Showing posts with label sequences. Show all posts
Showing posts with label sequences. Show all posts

Monday, April 5, 2010

Example 7.31: Contour plot of BMI by weight and height

A contour plot is a simple way to plot a surface in two dimensions. Lines with a constant Z value are plotted on the X-Y plane.

Typical uses include weather maps displaying "isobars" (lines of constant pressure), and maps displaying lines of constant elevation useful in, e.g., hiking. Unusual examples include maps of constant travel time, for example this map, which transforms space so that constant travel time can be presented as circles, or Francis Galton's travel map, below, in which the lines of constant travel time are reflected in the borders between colors.




As an example, here, we plot iso-BMI lines showing regions of underweight, normal weight, overweight, and obesity. As presented by National Heart, Lung, and Blood Institute, the borders between these categories are at BMI of 18.5, 25, and 30. Leaving aside questions of whether risk changes precipitously at these conveniently round numbers, it would be helpful to see what the categories imply across a range of heights and weights.




SAS

In SAS, we loop (section 1.11.1) through a range of heights (English) and weights (avoirdupois), calculating the BMI for each.


data bmi;
do ht = 48 to 84 by ((84-48)/1000);
do wt = 90 to 300 by ((300 -90)/1000);
bmi = (wt*703)/(ht**2);
output;
end;
end;
run;


We then prepare some labels for the plot, using the annotate macros (section 5.2, 5.2.11).


%annomac;
data annobmi;
%system(2,2,2);
%label(75,120,"Underweight",black,18,0,17);
%label(70,150,"Normal",black,33,0,17);
%label(65,170,"Overweight",black,37,0,17);
%label(55,200,"Obese",black,40,0,17);
run;


Finally, we use proc contour (section 5.1.12) to make the contour plot, using the anno option to plot the labels.


title h=1.25 "BMI categories by height and weight";
proc gcontour data = bmi anno=annobmi;
plot wt * ht = bmi / levels = 18.5,25,30 nolegend;
label ht="Height (inches)" wt="Weight (lbs)";
run;


R

In R, we build the data using the seq() function (section 1.11.3) for height and weight separately, and generate a data frame with all combinations of height and weight with the expand.grid() function (see example 7.22). After defining a function to calculate the BMI once, we apply it to the whole data set, storing the results in the matrix form we'll need using the matrix() function (section 1.9.1).


ht = seq(48,84, length.out=1000)
wt = seq(90,300, length.out=1000)
wtht = expand.grid(x=ht, y=wt)
bmi = function(h,w) {(w * 703)/(h*h)}
bmiwtht = matrix(bmi(wtht$x,wtht$y),length(ht),length(wt))


Now we're ready to use the contour() function (section 5.1.12) to make the plot and to add the labels with the text() function (section 5.2.11).


contour(ht,wt,bmiwtht,levels = c(18.5,25,30), drawlabels=FALSE,
xlab="Height (inches)",ylab="Weight (lbs)",
main="BMI categories by height and weight")

text(55,200,"Obese",cex=2,srt=45)
text(65,165,"Overweight",cex=2,srt=40)
text(70,150,"Normal",cex=2,srt=35)
text(75,120,"Underweight",cex=2,srt=18)



Saturday, July 25, 2009

Example 7.7: Tabulate binomial probabilities

Suppose we wanted to assess the probability P(X=x) for a binomial random variate with n = 10 and with p = .81, .84, ..., .99. This could be helpful, for example, in various game settings.

In SAS, we find the probability that X=x using differences in the CDF calculated via the cdf function (1.10.1). We loop through the various binomial probabilities and values of x using the do ... end structure (section 1.11.1). Finally, we store the probabilities in implicitly named variables via an array statement (section 1.11.5).

SAS

data table (drop=j);
array probs [11] prob0 prob1 - prob10;
do p = .81 to .99 by .03;
do j = 0 to 10;
if j eq 1 then probs[j+1] = cdf("BINOMIAL", 0, p, 10);
else probs[j+1] = cdf("BINOMIAL", j, p, 10)
- cdf("BINOMIAL", j-1, p, 10);
end;
output;
end;
run;


The results are printed with two decimal places using the format statement (section 1.2.4). The options statement (section A.4) uses the ls option to narrow the output width to be compatible with Blogger.


options ls=64;
proc print data=table noobs;
var p prob0 prob1 - prob10;
format _numeric_ 3.2;
run;


And the results are:

p
p p p p p p p p p p r
r r r r r r r r r r o
o o o o o o o o o o b
b b b b b b b b b b 1
p 0 1 2 3 4 5 6 7 8 9 0
.81 .00 .00 .00 .00 .00 .02 .08 .19 .30 .29 .12
.84 .00 .00 .00 .00 .00 .01 .05 .15 .29 .33 .17
.87 .00 .00 .00 .00 .00 .00 .03 .10 .25 .37 .25
.90 .00 .00 .00 .00 .00 .00 .01 .06 .19 .39 .35
.93 .00 .00 .00 .00 .00 .00 .00 .02 .12 .36 .48
.96 .00 .00 .00 .00 .00 .00 .00 .01 .05 .28 .66
.99 .00 .00 .00 .00 .00 .00 .00 .00 .00 .09 .90


R
In R we start by making a vector of the binomial probabilities, using the : operator (section 1.11.3) to generate a sequence of integers. After creating a matrix (section B.4.4) to hold the table results, we loop (section 1.11.1) through the binomial probabilities, calling the dbinom() function (section 1.1) to find the probability that X=x. This calculation is nested within the round() function (section 1.2.5) to reduce the digits displayed. Finally, we glue the vector of binomial probabilities to the results.


p <- .78 + (3 * 1:7)/100
allprobs <- matrix(nrow=length(p), ncol=11)
for (i in 1:length(p)) {
allprobs[i,] <- round(dbinom(0:10, 10, p[i]),2)
}
table <- cbind(p, allprobs)
table


With the result:

p
[1,] 0.81 0 0 0 0 0 0.02 0.08 0.19 0.30 0.29 0.12
[2,] 0.84 0 0 0 0 0 0.01 0.05 0.15 0.29 0.33 0.17
[3,] 0.87 0 0 0 0 0 0.00 0.03 0.10 0.25 0.37 0.25
[4,] 0.90 0 0 0 0 0 0.00 0.01 0.06 0.19 0.39 0.35
[5,] 0.93 0 0 0 0 0 0.00 0.00 0.02 0.12 0.36 0.48
[6,] 0.96 0 0 0 0 0 0.00 0.00 0.01 0.05 0.28 0.66
[7,] 0.99 0 0 0 0 0 0.00 0.00 0.00 0.00 0.09 0.90


As with the example of jittered scatterplots for dichotomous y by continuous x, (Example 7.3, Example 7.4, and Example 7.5) we will revisit this example for improvement in later entries.