SAS and R

Example 10.8: The upper 95% CI is 3.69

2012-12-10T10:08:00.000-05:00

Apologies for the long and unannounced break-- the longest since we started blogging, three and a half years ago. I was writing a 2-day course for SAS users to learn R. Contact me if you're interested. And Nick and I are beginning work on the second edition of our book-- look for it in the fall. Please let us know if you have ideas about what we omitted last time or would otherwise like to see added. In the mean time, we'll keep blogging, though likely at a reduced rate.

Today: what can you say about the probability of an event if the observed number of events is 0? It turns out that the upper 95% CI for the probability is 3.69/N. There's a sweet little paper with some rationale for this, but it's in my other office. And I couldn't recall the precise value-- so I used SAS and R to demonstrate it to myself.

R

The R code is remarkably concise. After generating some Ns, we write a little function to perform the test and extract the (exact) upper 95% confidence limit. This is facilitated by the "..." notation, which passes along unused arguments to functions. Then we use apply() to call the new function for each N, passing the numerator 0 each time. Note that apply() needs a matrix argument, so the simple vector of Ns is converted to a matrix before use. [The sapply() function will accept a vector input, but took about 8 times as long to run.] Finally, we plot the upper limit * N against N. showing the asymptote. A log scaled x-axis is useful here, and is achieved with the log='x' option. (Section 5.3.12.) the result is shown above.

bin.m = seq(10, 10000, by=5)
mybt = function(...) { binom.test(...)$conf.int[2] }
uci = apply(as.matrix(bin.m), 1, mybt, x=0)
plot(y=bin.m * uci, x=bin.m, ylim=c(0,4), type="l", 
     lwd=5, col="red", cex=5, log='x',  
     ylab="Exact upper CI", xlab="Sample size", 
     main="Upper CI when there are 0 cases observed")
abline(h=3.69)

SAS

In SAS, the data, really just the N and a numerator of 0, are generated in a data step. The CI are found using the binomial option in the proc freq tables statement and saved using the output statement. Note that the weight statement is used here to avoid having a row for each Bernoulli trial.

data binm;
do n = 10 to 10000 by 5;
  x=0;
  output;
  end;
run;

ods select none;
proc freq data=binm;
by n;
weight n;
tables x / binomial;
output out=bp binomial;
run;
ods select all;

To calculate the upper limit*N, another data step is needed-- note that in this setting SAS will only produce the lower limit against the probability that all observations share the same value, thus the subtraction from 1 shown below. The log scale x-axis is obtained with the logbase option to the axis statement. (Section 5.3.12.) The result is shown below.

data uci;
set bp;
limit = (1-xl_bin) * n;
run;

axis1 order = (0 to 4 by 1);
axis2 logbase=10 logstyle=expand;
symbol1 i = j v = none c = red w=5 l=1;
proc gplot data=uci;
plot limit * n / vref=3.69 vaxis=axis1 haxis=axis2;
label n="Sample size" limit="Exact upper CI";
run;
quit;

It's clear that the upper 95% limit on the number of successes asymptotes to about 3.69. Thus the upper limit on the binomial probability p is 3.69/N.

An unrelated note about aggregators: We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 10.7: Fisher vs. Pearson

2012-10-29T09:33:00.000-04:00

In the early days of the discipline of statistics, R.A. Fisher argued with great vehemence against Egon Pearson (and Jerzy Neyman) over the foundational notions supporting statistical inference. The personal invective recorded is somewhat amusing and also reminds us how very puerile even very smart people can be.

Today, we'll compare Fisher's exact test for 2*2 tables with the Pearson chi-square, developed by Karl Pearson, Egon's father and another early pioneer of statistics. This blog entry was inspired by a questioner on LinkedIn who asked when should the one be preferred over the other. One commenter gave the typical rule of thumb-- "If the expected count in any cell is less than 5, use the exact test, otherwise use the chi-square." My copy of "Statistical rules of thumb" is AWOL at the moment, so I don't know if this one is covered there. A quick googling did not reveal an answer either.

The rule of thumb dates back to the days before the exact test became computationally feasible outside of small problems. In contrast, today it can be performed quickly for all tables, either through a complete enumeration of the possible tables or through Monte Carlo hypothesis testing, which is simple to apply in either SAS or R. My default in recent years has been to take advantage of this capability and use the exact test all the time, ignoring the traditional approximate chi-square test. My idea was that if there were any small cells, I'd be covered, while allowing a simpler methods section.

Let's develop some code to see what happens. Is the rule of thumb accurate? What's the power cost of using the exact test instead of the Chi-square?

Our approach will be to set cell probabilities and the sample size, and simulate data under this model, then perform each test and evaluate the proportion of rejections under each test. One complication here is that simulated data might result in a null margin, i.e., there might be no observed values in a row or in a column. We'll calculate rejections of the null only among the tables where this does not happen. This means that the average observed cell counts among included tables may be different from the expect cell counts. This makes sense from a practical perspective-- we probably would not do the test if we observed 0 subjects in one of our planned categories.

SAS
In SAS, we'll do the dumb straightforward thing and simulate 100 pairs of dichotomous variables. Here we just code the null case of no association, with margins of 70% in one column and 80% in one row. The smallest cell has an expected count of 6%, so that a total sample size of 83 will have an expected count of 5 in that cell.

data test;
pdot1 = .7;
p1dot = .8;
do tablen = 20, 50, 100, 200, 500, 1000;
  do ds = 1 to 10000;
    do i = 1 to tablen;
      xnull = uniform(0) gt pdot1;
      ynull = uniform(0) gt p1dot;
      output;
      end;
    end;
  end;
run;

Then proc freq can be used to generate the two p-values, using the by statement to do the calculations for all the tables at once. The output statement extracts the p-values into a data set.

ods select none;
options nonotes;
proc freq data = test;
by tablen ds;
tables ynull * xnull / chisq fisher;
output out = kk1 chisq fisher;
run;
options notes;
ods select all;

To get the proportion of rejections, we first use a data step to calculate whether each test was rejected, then go back to proc freq to find the proportion of rejections and the CI on the probability of rejections.

data summ;
set kk1 (keep = tablen p_pchi xp2_fish);
rej_pchi = (p_pchi lt 0.05);
rej_fish = (xp2_fish lt .05);
run;

ods output binomialprop = kk2;
proc freq data = summ;
by tablen;
tables rej_pchi rej_fish / binomial(level='1');
run;

You may have noticed that we didn't do anything to account for the tables with empty rows or columns. When the initial proc freq encounters such a table, it performs neither test. Thus the second proc freq is calculating the proportion and CI with a denominator that might be smaller than the number of tables we simulated. Happily, they'll still be correct, though the CI may be wider than we'd intended. Finally, we're ready to plot the results, using the hilob interpolation described in Example 10.4. Using hiloc instead shows the "close" as a tick mark between the high and low values.

data kk2a;
set kk2;
if table eq "Table rej_pchi" then tablen = tablen + 1;
run;

symbol1 i = hiloc;
symbol2 i = hiloc;
proc gplot data = kk2a;
where name1 in ("_BIN_","XL_BIN","XU_BIN");
plot nvalue1 * tablen = table / vref = 0.05 href=83;
/* ref lines where the expected count in the smallest cell is > 5, 
and the nominal alpha */
run; quit;

The results are shown above. The confidence limits should include 0.05 for all numbers of subjects in order to be generally recommended. Both tests reach this standard, with these margins, even for tables with only 20 subjects, i.e., with expected cell counts of 11, 5, 3, and 1. The exact test appears conservative (rejects less than 5% of the time), probably due to small cell counts and the resulting ties in the list of possible tables.

R
In R, we'll simulate observations from a multinomial distribution with the desired cell probabilities, and assemble the result into a table to calculate the p-values. This will make it easier to simulate tables under the alternative, as we need to do to assess power. If there are empty rows or columns, the chisq.test() function produces a p-value of "NaN", which will create problems later. To avoid this, we'll put the table generation inside a while() function. This operates like the do while construction in SAS (and other programming languages). The condition we check for is whether there is a row or column with 0 observations; if so, try generating the data again. We begin by initializing the table with 0's.

makeitm = function(n.in.table, probs)  {
   myt = matrix(rep(0,4), ncol=2)
  while( (min(colSums(myt)) == 0) | (min(rowSums(myt)) == 0) )  { 
    myt = matrix(rmultinom(n=1, size=n.in.table, probs), ncol=2,byrow=TRUE)
}
  chisqp = chisq.test(myt, correct=FALSE)$p.value
  fishp = fisher.test(myt)$p.value
  return(c(chisqp, fishp))
}

With this basic building block in place, we can write a function to call it repeatedly (using the replicate() function, then calculate the proportion of rejections and get the CI for the probability of rejections.

many.ns = function(tablen, nds,probs) {
  res1 = t(replicate(nds,makeitm(tablen,probs)))
  res2 = res1 < 0.05
  pear = sum(res2[,1])/nds
  fish = sum(res2[,2])/nds
  pearci = binom.test(sum(res2[,1]),nds)$conf.int[1:2]
  fishci = binom.test(sum(res2[,2]),nds)$conf.int[1:2]
  return(c("N.ds" = nds, "N.table" = tablen, probs, 
          "Pearson.p" = pear, "PCIl"=pearci[1], "PCIu"=pearci[2],
          "Fisher.p" = fish, "FCIl" = fishci[1], "FCIu" = fishci[2]))
}

Finally, we're ready to actually do the experiment, using sapply() to call the function that calls replicate() to call the function that makes a table. The result is converted to a data frame to make the plotting simpler. The first call below replicates the SAS result shown above and has very similar estimates and CI, but is not displayed here. The second shows an odds ratio of 3, the third (plotted below) has an OR of 1.75, and the last an OR of 1.5.

#res3 = data.frame(t(sapply(c(20, 50, 100, 200, 500, 1000),many.ns, nds=10000, 
  probs = c(.56,.24,.14,.06))))
#res3 = data.frame(t(sapply(c(20, 50, 100, 200, 500, 1000),many.ns, nds=1000, 
  probs = c(.6,.2,.1,.1))))
res3 = data.frame(t(sapply(c(20, 50, 100, 200, 500, 1000),many.ns, nds=1000, 
  probs = c(.58,.22,.12,.08))))
#res3 = data.frame(t(sapply(c(20, 50, 100, 200, 500, 1000),many.ns, nds=1000, 
  probs = c(.57,.23,.13,.07))))

with(res3,plot(x = 1, y =1, type="n", ylim = c(0, max(c(PCIu,FCIu))), xlim=c(0,1000),
               xlab = "N in table", ylab="Rejections", main="Fisher vs. Pearson"))
with(res3,points(y=Pearson.p, x=N.table,col="blue"))
with(res3,segments(x0 = N.table, x1=N.table, y0 = PCIl, y1= PCIu, col = "blue"))
with(res3,points(y=Fisher.p, x=N.table + 2,col="red"))
with(res3,segments(x0 = N.table+2, x1=N.table+2, y0 = FCIl, y1= FCIu, col = "red"))
abline(v=83)
abline(h=0.05)

The plotting commands used above have been demonstrated in Examples 10.4 and 8.42.

Overall, the results show (in the SAS plot, at the top) that the Pearson chi-square test does perfectly well at protecting the alpha level under the null with these margins, even when the expected number of cases in one cell is as small as 1. In contrast, compared to the exact test, the Chi-square test has a bit more power, for these cell probabilities. The power difference is greatest when the N is smaller. Given this example, I would say that the rule of thumb may be too conservative, pushing people away from a more powerful test unnecessarily. In the future, I plan to be more positive about using the Pearson chi-square.

An unrelated note about aggregators: We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 10.6: Should Poisson regression ever be used? Negative binomial vs. Poisson regression

2012-10-15T10:30:00.000-04:00

In practice, we often find that count data is not well modeled by Poisson regression, though Poisson models are often presented as the natural approach for such data. In contrast, the negative binomial regression model is much more flexible and is therefore likely to fit better, if the data are not Poisson. In example 8.30 we compared the probability mass functions of the two distributions, and found that for a given mean, the negative binomial closely approximates the Poisson, as the scale parameter increases. But how does this affect the choice of regression model? How might another alternative, the overdispersed, or quasi-Poisson model compete with these? Today we generate a rudimentary toolkit for assessing the effects of Poisson, negative binomial, and quasi-Poisson models, assuming data are truly generated by one or the other process.

SAS
We'll begin by simulating Poisson and negative binomial data. Note that we also rely on the poismean_nb function that we created in example 8.30-- this is needed because SAS only accepts the natural parameters of the distribution, while the mean is a (simple) function of the two parameters. As is typical in such settings, we'll begin by generating data under the null of no association between, in this case, the normal covariate and the count outcome. The proportion of rejections should be no greater than alpha (5%, here). However, we'll include code to easily simulate data under the alternative as well. This will facilitate assessing the relative power of the models, later.

data nbp;
do ds = 1 to 10000;
  do i = 0 to 250;
    x = normal(0);
    mean = exp(-.25 + (0 * x));
    pois = rand("POISSON",mean);
    nb1 = rand("NEGBINOMIAL", poismean_nb(mean, 1), 1);
    output;
    end;
  end;
run;

The models will be fit in proc genmod. (See sections 4.1.3, 4.1.5, table 4.1.) It would be good to write a little macro to change the distribution and the output names, but it's not necessary. To save space here, the repetitive lines are omitted. The naming convention is that the true distribution (p or nb) is listed first, followed by the fit model (p, nb, or pod, for overdispersed).

 
ods select none;
ods output parameterestimates = ppests;
proc genmod data = nbp;
by ds;
model pois = x /dist=poisson;
run;

ods output parameterestimates = ppodests;
model pois = x /dist=poisson scale = p;

ods output parameterestimates = pnbests;
model pois = x /dist=negbin;

ods output parameterestimates = nbnbests;
model nb1 = x /dist=negbin;

ods output parameterestimates = nbpests;
model nb1 = x /dist=poisson;

ods output parameterestimates = nbpodests;
model nb1 = x /dist=poisson scale=p;
ods select all;

For analysis, we'll bring all the results together using the merge statement (section 1.5.7). Note that the output data sets contain the Wald CI limits as well as the estimates themselves; all have to be renamed in the merge, or they will overwrite each other.

data results;
merge
ppests (rename = (estimate = pp_est lowerwaldcl = pp_l 
  upperwaldcl = pp_u))
ppodests (rename = (estimate = ppod_est lowerwaldcl = ppod_l 
  upperwaldcl = ppod_u))
pnbests (rename = (estimate = pnb_est lowerwaldcl = pnb_l 
  upperwaldcl = pnb_u))
nbnbests (rename = (estimate = nbnb_est lowerwaldcl = nbnb_l 
  upperwaldcl = nbnb_u))
nbpests (rename = (estimate = nbp_est lowerwaldcl = nbp_l 
  upperwaldcl = nbp_u))
nbpodests (rename = (estimate = nbpod_est lowerwaldcl = nbpod_l 
  upperwaldcl = nbpod_u));
where parameter eq "x";
target = 0;
pprej =   ((pp_l gt target) or (pp_u lt target));
ppodrej =   ((ppod_l gt target) or (ppod_u lt target));
pnbrej =  ((pnb_l gt target) or (pnb_u lt target));
nbnbrej = ((nbnb_l gt target) or (nbnb_u lt target));
nbprej =  ((nbp_l gt target) or (nbp_u lt target));
nbpodrej =  ((nbpod_l gt target) or (nbpod_u lt target));
run;

The indicators of CI that exclude the null are calculated with appropriate names using logical tests that are 1 if true (rejections) and 0 if false. (See, e.g., section 1.4.9.) The final results can be obtained from proc means

proc means data = results; 
var pp_est ppod_est pnb_est nbnb_est nbp_est nbpod_est 
    pprej ppodrej pnbrej nbnbrej nbprej nbpodrej; 
run;

                     Variable             Mean
                     -------------------------
                     pp_est       -0.000349738
                     ppod_est     -0.000349738
                     pnb_est      -0.000344668
                     nbnb_est        0.0013738
                     nbp_est         0.0013588
                     nbpod_est       0.0013588
                     pprej           0.0505000
                     ppodrej         0.0501000
                     pnbrej          0.0468000
                     nbnbrej         0.0535000
                     nbprej          0.1427000
                     nbpodrej        0.0555000
                     -------------------------

All of the estimates appear to be unbiased. However, Poisson regression, when applied to the truly negative binomial data, appears to be dramatically anticonservative, rejecting the null (i.e., with CI excluding the null value) 14% of the time. The overdispersed model may be slightly biased as well. The estimated proportion of rejections is 5.55%, or 555 of 10,000 experiments. An exact CI for the proportion excludes 5%, here, although the anticonservative bias appears to be slight. To test other effect sizes, we'd change the mean, set in the first data step and the target in the results data. It would also be valuable to change the scale parameter for the negative binomial.

R
We begin by defining two simple functions: one to extract the standard errors from a model, and the second to assess whether Wald-type CI for parameter estimates exclude some value. It's a bit confusing that a standard error extracting function is not part of R. Or perhaps it is, and someone will point out the obvious function in the comments. It's useful to use the standard errors and construct the Wald CI in the current setting because the obvious alternative for constructing CI, the confint() function, uses profile likelihoods, which would be too time-consuming in a simulation setting. The second function accepts the parameter estimate, its standard error, and a fixed value which we want to know is in or out of the CI. Both functions are actually single expressions, but having them in hand will reduce the typing in the main function.

# this will work for any model object that works with vcov()
# the test for positive variance should be unnecessary but can't hurt
stderrs = function(model) {
  ifelse(min(diag(vcov(model))) > 0, sqrt(diag(vcov(model))), NA)  
}

# short and sweet: 1 if target is out of Wald CI, 0 if in
ciout = function(est, se, target){
  ifelse( (est - 1.96*se > target) | (est + 1.96*se < target), 1,0)
}

With these ingredients prepared, we're ready to write a function to fit the three models to the two sets of observed data. The function will accept a number of observations per data set and a true beta. The Poisson and overdispersed Poisson are fit with the glm() function (section 4.1.3, table 4.1) but the negative binomial uses the glm.nb() function found in the MASS package (section 4.1.5).

testpnb = function(n, beta) {
# make data
n = 250
x = rnorm(n)
mean = exp(-.25 + (beta * x))
pois = rpois(n,mean)
nb1 = rnbinom(n, size=1, mu=mean)

# fit models to Poisson data
pp = glm(pois ~x, family="poisson")
ppod = glm(pois ~x, family="quasipoisson")
pnb = glm.nb(pois~x)

# fit models to nb data
nbnb = glm.nb(nb1 ~x)
nbp = glm(nb1 ~x, family="poisson")
nbpod = glm(nb1 ~x, family="quasipoisson")

# extract parameter estimates using the coef() function
est = as.numeric(c(coef(pp)[2], coef(ppod)[2], coef(pnb)[2], coef(nbnb)[2], coef(nbp)[2], coef(nbpod)[2]))
# use our two new functions to get the SE and the CI indicator
se = c(stderrs(pp), stderrs(ppod), stderrs(pnb), stderrs(nbnb), stderrs(nbp), stderrs(nbpod))
ci = ciout(est, se, rep(beta,6))
return(matrix(c(est,se,ci),ncol=3))
}

Now we can use the convenient replicate() function to call the experiment many times. Since the output of testnb() is a matrix, the result of replicate() is a three-dimensional matrix, R * C * sheet, where sheet here corresponds to each experimental replicate. To summarize the results, we can use the rowMeans() function to get the proportion of rejections or the mean of the estimates.

mainout = replicate(10000,testpnb(250,0))

# the [,3,] below means "all rows, column 3, all sheets"
> rowMeans(mainout[,3,])
[1] 0.0490 0.0514 0.0463 0.0490 0.1403 0.0493

> rowMeans(mainout[,1,])
[1]  0.0003482834  0.0003482834  0.0003558526 -0.0004123949 -0.0003972441 -0.0003972441

The results agree completely with the SAS results discussed above. The naive Poisson regression would appear a bad idea--if the data are negative binomial, tests don't have the nominal size. It would be valuable to replicate the experiment with some other distribution for the real data as well. One approach to modeling count data would be to fit the Poisson and assess the quality of the fit, which can be done in several ways. However, this iterative fitting also jeopardizes the size of the test, in theory. Perhaps we'll explore the practical impact of this in a future entry. Fortunately, at least in this limited example, a nice alternative exists: We can just fit the negative binomial by default. The costs of this in terms of power could be assessed with a thorough simulation study, but are likely to be small, since only one additional parameter is estimated. And the size of the test is hardly affected at all. The quasi-Poisson model could also be recommended, but has the drawback of relying on what is actually not a viable distribution for the data. Some sources suggest that it may be even more flexible than the negative binomial, however.

An unrelated note about aggregators: We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 10.5: Convert a character-valued categorical variable to numeric

2012-10-08T10:58:00.000-04:00

In some settings it may be necessary to recode a categorical variable with character values into a variable with numeric values. For example, the matching macro we discussed in example 7.35 will only match on numeric variables. One way to convert character variables to numeric values is to determine which values exist, then write a possibly long series of conditional tests to assign numbers to the values. Surely there's a better way?

SAS
In SAS, Rick Wicklin offers an IML solution and links to a macro with the same function. But if you're not an IML coder, and you don't want to investigate a macro solution, it's simple enough to do with data steps. We'll begin by making some fake data.

data test;
  do i = 1 to 100;
  cat = "meow";
  if i gt 30 then cat = "Purr";
  if i gt 70 then cat = "Hiss";
  output;
  end;
run;

To make the new variable, we'll just sort (section 1.5.6) the data on the categorical variable we want to convert, then use the set ds; by x; syntax to keep track of when a new value is encountered in the data. It's hard to believe that we've never demonstrated this useful syntax before-- perhaps we just can't find it today. The set ds; by x; syntax makes new temporary variables first.x and last.x that are equal to 1 for the first and last observations of each new level of x, respectively, and 0 otherwise. When we find a new value, we'll increase a counter by 1; the counter is our new numeric-valued variable.

proc sort data = test; by cat; run;

data catize;
set test;
by cat;
retain catnum 0;
if first.cat then catnum = catnum + 1;
run;

/* check the result */
proc freq data = catize;
tables cat * catnum;
run;

The table also shows the recoding values.

                             Table of cat by catnum

                  cat       catnum

                  Frequency|
                  Percent  |
                  Row Pct  |
                  Col Pct  |       1|       2|       3|  Total
                  ---------+--------+--------+--------+
                  Hiss     |     30 |      0 |      0 |     30
                           |  30.00 |   0.00 |   0.00 |  30.00
                           | 100.00 |   0.00 |   0.00 |
                           | 100.00 |   0.00 |   0.00 |
                  ---------+--------+--------+--------+
                  Purr     |      0 |     40 |      0 |     40
                           |   0.00 |  40.00 |   0.00 |  40.00
                           |   0.00 | 100.00 |   0.00 |
                           |   0.00 | 100.00 |   0.00 |
                  ---------+--------+--------+--------+
                  meow     |      0 |      0 |     30 |     30
                           |   0.00 |   0.00 |  30.00 |  30.00
                           |   0.00 |   0.00 | 100.00 |
                           |   0.00 |   0.00 | 100.00 |
                  ---------+--------+--------+--------+
                  Total          30       40       30      100
                              30.00    40.00    30.00   100.00

R
We begin by making the data. To convert to numbers, we use the labels option to the factor() function, feeding it the sequences of numbers between 1 and however many different values there are. Note that we find this using the factor() function again. There's probably a better way of doing this, but it's a little bit amusing to code it this way. Then we have numbers, but they're store as a factor. We can get them out with a call to as.numeric().

cat = c(rep("meow",30),rep("Hiss",30), rep("Purr", 40))
catn1 = factor(cat, labels=(1:length(levels(factor(cat)))))
catn = as.numeric(catn1)
table(catn,cat)

    cat
catn Hiss meow Purr
   1   30    0    0
   2    0   30    0
   3    0    0   40

There's a warning in the documentation for factor() that the values are assigned in location-specific fashion, so the table should be used to establish how the codes were assigned. For the record, the use cases for this kind of recoding in R may be more strained than the SAS example given above.

An unrelated note about aggregators: We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 10.4: Multiple comparisons and confidence limits

2012-10-01T12:18:00.000-04:00

A colleague is a devotee of confidence intervals. To him, the CI have the magical property that they are immune to the multiple comparison problem-- in other words, he feels its OK to look at a bunch of 95% CI and focus on the ones that appear to exclude the null. This though he knows well the one-to-one relationship between 95% CIs that exclude the null and p-values below 0.05.

Today, we'll create a Monte Carlo experiment to demonstrate that fishing by CI is just as dangerous as fishing by p-value; generating the image above. We'll do this by replicating a bivariate experiment 100 times. Later, we'll examine the results of a single experiment with many predictors.

R
To begin with, we'll write a function to generate a single experiment, using a logistic regression. This is a simple modification of one of our first and most popular entries.

simci = function(){
  intercept = 0
  beta = 0
# beta = 0 because we're simulating under the null 
# make the variance of x in this experiment vary a bit
  xtest = rnorm(1000) * runif(1,.6,1.4)
  linpred = intercept + xtest*beta
  prob = exp(linpred)/(1 + exp(linpred))
  runis = runif(1000)
  ytest = ifelse(runis < prob,1,0)
# now, fit the model
  fit = glm(ytest~xtest,family=binomial(link="logit"))
# the standard error of the estimates is easiest to find in the
  pe = summary(fit)$coefficients
# calculate the Wald CI; an alternative would be confint(), but
# that calculated profile CI, which take longer to generate
  ci = exp(c(pe[2,1] - 1.96*pe[2,2], pe[2,1] + 1.96*pe[2,2] ))
  return(ci)
}

Then we can use the replicate() function to repeat the experiment 100 times. The t() function (section 1.9.2) transposes the resulting matrix to have one row per experiment.

sim100 = t(replicate(100,simci()))

plot(x = sim100[,1], y = 1:100, 
  xlim = c(min(sim100),max(sim100)), type="n")
segments(y0=1:100,x0=sim100[,1],y1 = 1:100,x1=sim100[,2], 
  col = ifelse(sim100[,1]>1 | sim100[,2]<1,"red","black"))
abline(v=1)

The result is shown at the top. In the code, we set the limits of the x-axis by finding the max and min across the whole matrix-- this is a little wasteful of CPU cycles, but saves some typing. The segments() function (see example 8.42) is a vector-enabled line-drawer. Here we draw a line from the lower CI limit to the upper, giving the experiment number as the x value for each. We assign a red plot line if the CI excludes 1, using the ifelse() function (section 1.11.2), a vectorwise logic test. Finally, a reference line helps the viewer see for far the end of the CI is from the null. We omit prettying up the axis labels.

SAS
In SAS, considerably more lines are required. We begin by simulating the data, as in example 7.2. The modifications are to generate 100 examples with an outside do loop (section 1.11.1) and the random element added to the variance.

data simci;
beta = 0;
intercept = 0;
do sim = 1 to 100;   /* outer loop */
  xvar = (uniform(0) *.8) + .6;  /* variance != 1 */
  do i = 1 to 1000;
    xtest = normal(0) * xvar;
    linpred = intercept + (xtest * beta);
    prob =  exp(linpred)/(1 + exp(linpred));
    ytest = (uniform(0) < prob);
    output;
  end;
end;
run;

Then we fit the logistic regression. We leave in the ods trace commands (section A.7.1) to remind you how to find the SAS names of the output elements, needed to save the results in the ods output statement. The CI for the odds ratios are requested in the clodds statement, which accepts a pl value for a profile likelihood based interval.

*ods trace on/listing; 
ods select none;
ods output cloddswald = lrci;
proc logistic data = simci;
by sim;
model ytest(event='1')=xtest / clodds=wald;
run;
*ods trace off;
ods select all;

Our plotting approach will require the "long" data set style, with two rows for each experiment. We'll generate that while checking whether the null is excluded from the CI.

data lrp2;
set lrci;
red = 0;
if lowercl > 1 or uppercl < 1 then red = 1;
point = lowercl; output;
point = uppercl; output;
run;

Finally, we're ready to make the graphic. We use the hilob interpolation to connect the upper and lower CI for each experiment; the b requests bars instead of a line, and the bwidth option specifies a very narrow bar. These options prevent the default plotting of the "mean" value with a tick. The a*b=c syntax (section 5.2.2) allows the different line colors.

symbol1 i=hilob bwidth=.05 c=black;
symbol2 i=hilob bwidth=.05 c=red;
proc gplot data = lrp2;
plot point * sim = red / vref = 1;
run;quit;

The result is just below. The vertical alignment seen in the R plot seems more natural, but this would not be possible with the hilo interpolation. As theory and logic would suggest, quite a few of the hundred simulated CI exclude the null, sometimes by a large proportion of the CI width.

An unrelated note about aggregators:We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 10.3: Enhanced scatterplot with marginal histograms

2012-09-24T12:32:00.001-04:00

Back in example 8.41 we showed how to make a graphic combining a scatterplot with histograms of each variable. A commenter suggested we change the R graphic to allow post-hoc plotting of, for example, lowess lines. In addition, there are further refinements to be made.

In this R-only entry, we'll make the figure more flexible and a bit more robust. See the example linked above for SAS code, or check out Rick Wicklin discussing the same subject-- Rick gives some additional resources.

R
The R code relies heavily on the layout() function. We discussed it last time in a simpler setting with only one column of plots. The goal for the current plot is to enable a title for the whole figure-- this ought to be centered over the whole graphic-- and x- and y-axis labels. In the previous version, there was no title to the page at all and the axis titles would occasionally fail. To do this, we need a layout with a single cell at the top for the whole width of the graphic, a tall narrow cell at the left for the y-axis title, only in the bottom row, and a thin cell at the bottom, only on the left, for the x-axis title. This turns out to be fairly simple with layout() and the results can be checked with layout.show().

zones <- matrix(c(1,1,1, 
                  0,5,0, 
                  2,6,4, 
                  0,3,0), ncol = 3, byrow = TRUE)
layout(zones, widths=c(0.3,4,1), heights = c(1,3,10,.75))
layout.show(n=6)

The matrix input tells R to make the whole first row a single plot area, and that this will be the first thing plotted. The corners of the remaining 3*3 plot cells will be empty. The numbers in the matrix give the order in which the plot cells will be filled. This matrix is the key input to layout(), where we use the remaining options to give the relative widths and heights of the cells. It's possible to do this in the abstract, but is helpful to draw the intended layout first, then test whether the intended design was a achieved using the layout.show() function. The result is shown below. Putting the scatterplot in last will be useful for adding to it post hoc.

With that in hand, it's time to make a function. In generating last week's example, we noted that the layout persists-- that is, the graphics area retains the layout until you shut the graphics device or restore the old parameters. In the new plot, we'll add an option to revert to the old parameters (by default) or retain them. The latter option would be desirable, if, as suggested by a commenter, you wanted to add items to the scatterplot after generating the plot. We also add an option to allow different sized plot symbols.

scatterhist <- function(x, y, xlab = "", ylab = "", plottitle="", 
                        xsize=1, cleanup=TRUE,...){
  # save the old graphics settings-- they may be needed
  def.par <- par(no.readonly = TRUE)
  
  zones <- matrix(c(1,1,1, 0,5,0, 2,6,4, 0,3,0), ncol = 3, byrow = TRUE)
  layout(zones, widths=c(0.3,4,1), heights = c(1,3,10,.75))
  
  # tuning to plot histograms nicely
  xhist <- hist(x, plot = FALSE)
  yhist <- hist(y, plot = FALSE)
  top <- max(c(xhist$counts, yhist$counts))
  
  # for all three titles: 
  #   drop the axis titles and omit boxes, set up margins
  par(xaxt="n", yaxt="n",bty="n",  mar = c(.3,2,.3,0) +.05)
  # fig 1 from the layout
  plot(x=1,y=1,type="n",ylim=c(-1,1), xlim=c(-1,1))
  text(0,0,paste(plottitle), cex=2)
  # fig 2
  plot(x=1,y=1,type="n",ylim=c(-1,1), xlim=c(-1,1))
  text(0,0,paste(ylab), cex=1.5, srt=90)
  # fig 3
  plot(x=1,y=1,type="n",ylim=c(-1,1), xlim=c(-1,1))
  text(0,0,paste(xlab), cex=1.5)
  
  # fig 4, the first histogram, needs different margins
  # no margin on the left
  par(mar = c(2,0,1,1))
  barplot(yhist$counts, axes = FALSE, xlim = c(0, top),
          space = 0, horiz = TRUE)
  # fig 5, other histogram needs no margin on the bottom
  par(mar = c(0,2,1,1))
  barplot(xhist$counts, axes = FALSE, ylim = c(0, top), space = 0)
  # fig 6, finally, the scatterplot-- needs regular axes, different margins
  par(mar = c(2,2,.5,.5), xaxt="s", yaxt="s", bty="n")
  # this color allows traparency & overplotting-- useful if a lot of points
  plot(x, y , pch=19, col="#00000022", cex=xsize, ...)
  
  # reset the graphics, if desired 
  if(cleanup) {par(def.par)}
}

To test this, we'll generate some data and try it out. The results are immediately below; I like this example to help demonstrate that it's not the marginal normality of the data that matter.

x=rexp(1000)
y = x^2 + rnorm(1000)
scatterhist(x[x<4], y[x<4], ylab="This is x", xlab="This is y",   
  "Revised scatterhist", xsize =2)

But let's take advantage of the ability to add curves to the scatterplot.

x=rexp(1000)
y = x^2 + rnorm(1000)
scatterhist(x[x<4], y[x<4], ylab="This is x", xlab="This is y",   
  "Revised scatterhist", xsize =2, cleanup=FALSE)
abline(lm(y~x))
lines(lowess(x,y))

The results are shown at the top-- we can do anything with the scatterplot that we'd be able to do if there were no layout() in effect.

An unrelated note about aggregators:We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 10.2: Custom graphic layouts

2012-09-17T11:17:00.004-04:00

In example 10.1 we introduced data from a CPAP machine. In brief, it's hard to tell exactly what's being recorded in the data set, but it seems to be related to the pattern of breathing. Measurements are taken five times a second, leading to on the order of 100,000 data points in a typical night. To get a visual sense of what a night's breathing looks like is therefore non-trivial.

Today, we'll make the graphic shown above, which presents an hour of data.

SAS
In SAS, the sgpanel procedure (section 5.1.11) will produce a similar graphic pretty easily. But we need to make a data set with indicators of the hour, and of ten-minute blocks within the hour. This we'll do with the ceil function (section 1.8.4).


data cycles2;
set cycles;
hour = ceil(time_min/60);
tenmin = ceil(time_min/10);
time_in_ten = mod(time_min - 1/300,10);
/* 1/300 adjustment keeps last measure in the correct 
        10-min block */
run;

title "Hour 4 of pressure";
proc sgpanel data = cycles2;
where hour eq 4;
panelby tenmin / layout=rowlattice rows=6 spacing = 4;
colaxis display=none;
rowaxis display = (nolabel);
series x = time_in_ten y = byte;
run; quit;

The resulting plot is shown below. It would be nicer to omit the labels on the right of each plot, but this does not appear to be an option. It would likely only be possible with a fair amount of effort.

R
In R, we'll use the layout() function to make a 7-row layout-- one for the title and 6 for the 10-minute blocks of time. Before we get there, though, we'll construct a function to fill the time block plots with input data. The function accepts a data vector and plots only 3,000 values from it, choosing the values based on an input hour and 10-minute block within the hour. To ensure an equal y-axis range for each call, we'll also send minimum and maximum values as input to the function. All of this will be fed into plot() with the type="l" option to make a line plot.


plot10 = function(hour, tenmins, miny, maxy, data=cycles){
   start = hour*18000 + tenmins* 3000 +1 
   plot((1:3000)/300, cycles[(start + 1):(start +3000)], 
            ylim = c(miny,maxy),type="l", xaxs="i", yaxs="i")
}

The documentation for layout() is rather opaque, so we'll review it separately.


oldpar = par(no.readonly = TRUE)
# revert to this later 

layout(matrix(1:7), widths=1, heights=c(3,8,8,8,8,8,8), respect=FALSE)

The layout() function divides the plot area into a matrix of cells, some of which will be filled by the next output plots. The first argument says where in the matrix the next N objects will go. All the integers 1...N must appear in the matrix; cells that will be left empty have a 0 instead. Here, we have no empty cells, and only one column, so the "matrix" is really just a vector with 1...7 in order. The widths option specifies the relative widths of the columns-- here we have only one column so any constant will result in the use of the whole width of the output area. Similarly, the heightsoption gives the relative height of the cells. Here the title will get 3/51 of the height, while each 10-minute block will get 8/51. This unequal shape of the plot regions is one reason to prefer layout() to some other ways to plot multiple images on a page. The respect option, when "TRUE" makes the otherwise relative widths and heights conform, so that a unit of height is equal to a unit of width. We also use layout() in example 8.41.

With the layout in hand, we're ready to fill it.


par(xaxt="n",  mar = c(.3,2,.3,0) +.05)
# drop the x-axis, change the spacing around the plot
plot(x=1,y=1,type="n",ylim=c(-1,1), xlim=c(-1,1), yaxt="n",bty="n")
# the first (narrow) plot is just empty
hour=3
text(0,0,paste("Hour ", (hour + 1), " of pressure data"), cex=2)
# text to put in the first plot
miny = min(cycles[(hour * 18000 + 1):((hour + 1) * 18000)])
maxy = max(cycles[(hour * 18000 + 1):((hour + 1) * 18000)])
# find min and max across the whole hour, to keep range 
# of y-axis constant across the plots
for (x in 0:5) plot10(hour, x, miny, maxy)
# plot the 6 ten-minute blocks
par(oldpar)
# reset the graphics options

The resulting plot is shown at the top of the entry. There's clearly something odd going on around 11-15 minutes into the hour-- this could be a misadjusted mask, or a real problem with the breathing. There's also a period around 58 minutes when it looks like breathing stops. That's what the machine is meant to stop.

An unrelated note about aggregatorsWe love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 10.1: Read a file byte by byte

2012-09-10T11:16:00.003-04:00

More and more makers of electronic devices use standard storage media to record data. Sometimes this is central to the device's function, as in a camera, so that the data must be easy to recover. Other times, it's effectively incidental, and the device maker may not provide easy access to the stored data.

For example, I recently was prescribed a constant positive air pressure (CPAP) machine for sleep apnea. My machine, made by Philips, records data onto a SD card. The card's file system is readable, but Philips provides neither software to read the files, nor a data dictionary to explain what they contain. (I believe Philips sells software that does read it, for ludicrous prices, to physicians who prescribe their machines. Nice racket.)

If you open the files on the card as ASCII files, you get a bunch of gobbledygook. But the data they contain is mine, in the most fundamental sense! I want to be able to read it. Fortunately, some folks have done a fair amount of work to reverse engineer the files. Through them, I was able to find some guidance for data from a related machine. Now I know what's in the file, more or less: a header of 25 bytes, 1200 bytes of data-- representing 4 minutes of recording--, and a one-byte footer, repeated ad nauseum (ad somnum? somno contingit?).

Today we show how to read bytes stored in a file as signed integers. For this file, (download,) we also trim out the header and footer, and make a simple line plot of the recorded data, which appear to be some function of the variable pressure with which the CPAP machine outputs air. Next time we'll make a more useful plot.

SAS
In SAS, we can use the infile statement to read in the data (section 1.1.2).


data test;
infile "c:\ken\cpap\0000000007.005"  recfm=n;
input byte ib1. @@;
run;

The recfm=n option tells SAS (for Windows, may differ in other OS) to read the file in binary. The ib1. informat tells SAS to read the bytes in the native format. (We cover reading in various formats in section 1.1.3, A.6.4, and several examples.) The @@ tells SAS to hold this line of input, rather than skipping to a new line, when the data is read. (See Example 8.1) . SAS will read bytes until there are no more to read.

Now I have a file with 128,680 observations, each being a signed integer. Some of these are actually nonsense, since the header and footer contain data stored in a variety of formats. To get rid of the header, we'll use the _n_ implied variable (section 1.4.15) which is effectively the line number, in conjunction with the mod function. While we're processing the data set anyhow, we'll also figure out the total elapsed time, which will be useful for plotting.


data cycles;
set test;
if mod(_n_,1226) ge 25 and mod(_n_,1226) lt 1225;
   /* otherwise it's a header or the footer */
time_min = (4 * int(_n_/1226)) +  (mod(_n_, 1226) - 24 )/(300);
   /* 4 minutes for each header-data-footer block + 
      number of measurements in this data block / 300 
      (measurements per minute).
run;

Now it's easy to plot the data-- a simple connected line plot across time makes sense, and can be made using the symbol statement with the i=j (j for join) syntax. The result is shown above.


symbol i = j v = none c = black;
proc gplot data = cycles;
where time_min le 4;
plot byte * time_min;
run;
quit;

R
In R, we'll use the readBin() function to actually read the file, but we need to do a little prep, first. The readBin() function requires we input the number of data elements to read. An overestimate is OK, but we can easily find the exact length of the file using the file.info function; the resulting object has a size constituent with the number of bytes. We'll also need a "connection" to the file, which is established in a call to the file() function.


finfo = file.info("0000000007.005")
toread= file("0000000007.005", "rb")
alldata = readBin(toread, integer(), size=1, n = finfo$size, endian="little")

The size option is the length of the elements, in bytes, and the endian option helps describe how the bytes should be read.

Analogous to SAS, the alldata vector has 128,680 integers. All that remains is to remove the headers and footers. We'll do that by making a logical test with the %% operator, saving the result as a vector, and selecting out the data from among the headers and footers using this logical. All that then remains is to plot the data-- we replicate the SAS plot of the first 1200 observations (4 minutes).


keep = 1:finfo$size %% 1226 > 24 & 1:finfo$size %% 1226 < 1225
cycles = alldata[keep]
# cycles gets only the elements of alldata when the corresponding
#    element of keep is TRUE

plot(1:1200, cycles[1:1200], type = "l")

The result is shown below. The first four minutes of my night's sleep were apparently characterized by generally lengthening breaths that became increasingly shallow.

An unrelated note about aggregatorsWe love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

The Statistical Sleuth (second edition) in R

2012-08-14T20:33:00.002-04:00

For those of you who teach, or are interested in seeing an illustrated series of analyses, there is a new compendium of files to help describe how to fit models for the extended case studies in the Second Edition of the Statistical Sleuth: A Course in Methods of Data Analysis (2002), the excellent text by Fred Ramsey and Dan Schafer. If you are using this book, or would like to see straightforward ways to undertake analyses in R for intro and intermediate statistics courses, these may be of interest.

These files can be found here. The site includes both formatted pdf files as well as the original knitr files which were used to generate the output. Knitr is an elegant, flexible and fast means to undertake reproducible analysis and dynamic report generation within R and RStudio.

This work leverages efforts undertaken by Project MOSAIC, an NSF-funded initiative to improve the teaching of statistics, calculus, science and computing in the undergraduate curriculum. In particular, we utilize the mosaic package, which was written to simplify the use of R for introductory statistics courses.

An unrelated note about aggregators:We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Third year wrap-up

2012-07-23T09:10:00.002-04:00

July marks the end of three years of blogging for us. By our count, we've posted 121 examples across the first three years. We aim to be helpful and interesting.

As always, it's hard to get a sense of our readership. At the time we wrote this, Feedburner reports about 1050 regular readers (up from 650 last year), but this (still) omits people who see us on R-bloggers or SAS Community Planet or SAS-X or statsblogs. As consumers of those aggregators, we assume there are many others who see us without subscribing directly. Google Analytics reports over 200,000 total pageviews (up from 100,000), while Feedburner claims 525,000, (up from 250,000).

As in previous years (2010, 2011) we report here on our most popular entries:

Feedburner
RStudio in the cloud, for dummies
To attach() or not attach(): that is the question
Example 9.17: (much) better pairs plot
Really useful R package: sas7bdat
Example 8.37: Read sheets from Excel
Example 9.20: Visualizing Simpson's paradox

Blogger
Example 7.35: Propensity score matching
Example 8.7: Hosmer and Lemeshow goodness-of-fit
Example 7.30: Simulate censored survival data
Example 7.38: Kaplan-Meier survival estimates
Example 7.2: Simulate data from a logistic regression

No overlap at all! This points to the difficulty of knowing what kinds of things we do may be useful to you, our readers. So, as usual, any feedback or suggestions would be most welcome.

In previous years we've slavishly turned over into a new set of numbered entries on our anniversary (July 1) and then taken a hiatus in August. This year we're going to rationalize and both conclude the chapter and take our break with this entry. We'll be back with example 10.1 in September.

SAS

An unrelated note about aggregators:We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.38: dynamite plots, revisited

2012-07-16T13:17:00.014-04:00

Dynamite plots are a somewhat pejorative term for a graphical display where the height of a bar indicates the mean, and the vertical line on top of it represents the standard deviation (or standard error). These displays are commonly found in many scientific disciplines, as a way of communicating group differences in means.

Many find these displays troubling. One post entitled them unmitigated evil.
The Vanderbilt University Department of Biostatistics has a formal policy discouraing use of these plots, stating that:

Dynamite plots often hide important information. This is particularly true of small or skewed data sets. Researchers are highly discouraged from using them, and department members have the option to decline participation in papers in which the lead author requires the use of these plots.

Despite the limitations of the display, we believe that there may be times when the display is helpful as a way to compare groups, assuming distributions that are approximately normal. Samuel Brown also described creation of these figures, as a way to encourage computing in R. We previously demonstrated how to create them in SAS and R, and today discuss code created by Randall Pruim to demonstrate how such graphics can be created using lattice graphics within the mosaic package.

R


library(mosaic)
dynamitePlot <- function(height, error, names = as.character(1:length(height)), 
                         significance = NA, ylim = c(0, maxLim), ...) {
  if (missing(error)) { error = 0 }
  maxLim <- 1.2* max(mapply(sum, height, error))
  mError <- min(c(error, na.rm=TRUE))
  barchart(height ~ names, ylim=ylim, panel=function(x,y,...) {
    panel.barchart(x,y,...)
    grid.polyline(c(x,x), c(y, y+error), id=rep(x,2), default.units='native',
      arrow=arrow(angle=45, length=unit(mError, 'native'))) 
    grid.polyline(c(x,x), c(y, y-error), id=rep(x,2), default.units='native',
      arrow=arrow(angle=45, length=unit(mError, 'native'))) 
    grid.text(x=x, y=y + error + .05*maxLim, label=significance, 
      default.units='native')
  }, ...)
}

Much of the code involves setting up the appropriate axis limits, then drawing the lines and adding the text. We can call this new function with an artificial example with 4 groups:


Values <- c(1,2,5,4)
Errors <- c(0.25, 0.5, 0.33, 0.12)
Names <- paste("Trial", 1:4)
Sig <- c("a", "a", "b", "b")
dynamitePlot(Values, Errors, names=Names, significance=Sig)

We still don't recommend frequent use of these plots (as other displays may be better (e.g. dotplots for very small sample sizes or violin plots), but having the capability to generate dynamite plots within the lattice framework can be handy.

An unrelated note about aggregators:We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.37: (Mis)behavior of binomial confidence intervals

2012-07-09T13:04:00.028-04:00

While traditional statistics courses teach students to calculate intervals and test for binomial proportions using a normal or t approximation, this method does not always work well. Agresti and Coull ("Approximate is better than "exact' for interval estimation of binomial proportions". The American Statistician, 1998; 52:119-126) demonstrated this and reintroduced an improved (Plus4) estimator originally due to Wilson (1927).

In this entry, we demonstrate how the coverage varies as a function of the underlying probability and compare four intervals: (1) t interval, (2) Clopper-Pearson, (3) Plus4 (Wilson/Agresti/Coull) and (4) Score, using code contributed by Albyn Jones from Reed College. Here, coverage probability is defined as the expected value that a CI based on an observed value will cover the true binomial parameter that generated that value. The code calculates the coverage probability as a function of a given binomial probability p and a sample size n. Intervals are created for each of the possible outcomes from 0, ..., n, then checked to see if the intervals include the true value. Finally, the sum of the probabilities of observing outcomes in which the binomial parameter is included in the interval determines the exact coverage. Note several distinct probabilities: (1) binomial parameter "p", the probability of success on a trial; (2) probability of observing x successes in N trials, P(X=x); (3) coverage probability as defined above. For distribution quantiles and probabilities, see section 1.10 and table 1.1.

R

We begin by defining the support functions which will be used to calculate the coverage probabilities for a specific probability and sample size.


CICoverage = function(n, p, alpha=.05) {
  # set up a table to hold the results
  Cover = matrix(0, nrow=n+1, ncol=5)
  colnames(Cover)=c("X","ClopperPearson","Plus4","Score","T")
  Cover[,1]=0:n
  zq = qnorm(1-alpha/2)
  tq = qt(1-alpha/2,n-1)
  for (i in 0:n) {
    Phat = i/n
    P4 = (i+2)/(n+4)
    # Calculate T and plus4 intervals manually, 
    # use canned functions for the other 
    TInt = Phat + c(-1,1)*tq*sqrt(Phat*(1-Phat)/n)
    P4Int = P4 + c(-1,1)*zq*sqrt(P4*(1-P4)/(n+4))
    CPInt= binom.test(i,n)$conf.int
    SInt = prop.test(i,n)$conf.int
    # check to see if the binomial p is in each CI 
    Cover[i+1,2] = InInt(p, CPInt)
    Cover[i+1,3] = InInt(p, P4Int)
    Cover[i+1,4] = InInt(p, SInt)
    Cover[i+1,5] = InInt(p, TInt)
  }
  # probability that X=x 
  p = dbinom(0:n, n, p)
  ProbCover=rep(0, 4)
  names(ProbCover) = c("ClopperPearson", "Plus4", "Score", "T")
  # sum P(X=x) * I(p in CI from x)
  for (i in 1:4){
    ProbCover[i] = sum(p*Cover[,i+1])
  }
  list(n=n, p=p, Cover=Cover, PC=ProbCover)
}

In addition, we define a function to determine whether something is in the interval.


InInt = function(p,interval){
  interval[1] <= p && interval[2] >= p
}

Finally, there's a function which summarizes the results.


CISummary = function(n, p) {
  M = matrix(0,nrow=length(n)*length(p),ncol=6)
  colnames(M) = c("n","p","ClopperPearson","Plus4","Score","T")
  k=0
  for (N in n) {
    for (P in p) {
      k=k+1
      M[k,]=c(N, P, CICoverage(N, P)$PC)
    }
  }
  data.frame(M)
}

We then generate the CI coverage plot provided at the start of the entry, which uses sample size n=50 across a variety of probabilities.


lwdval = 2
nvals = 50
probvals = seq(.01, .30, by=.001)
results = CISummary(nvals, probvals)
plot(range(probvals), c(0.85, 1), type="n", xlab="true binomial p",
     ylab="coverage probability")
abline(h=0.95, lty=2)
lines(results$p, results$ClopperPearson, col=1, lwd=lwdval)
lines(results$p, results$Plus4, col=2, lwd=lwdval)
lines(results$p, results$Score, col=3, lwd=lwdval)
lines(results$p, results$T, col=4, lwd=lwdval)
tests = c("ClopperPearson", "Plus4", "Score", "T")
legend("bottomright", legend=tests,
       col=1:4, lwd=lwdval, cex=0.70)

The resulting plot is quite interesting, and demonstrates how non-linear the coverage is for these methods, and how the t (almost equivalent to the normal, in this case) is anti-conservative in many cases. It also confirms the results of Agresti and Coull, who concluded that for interval estimation of a proportion, coverage probabilities from inverting the standard binomial and too small when inverting the Wald large-sample normal test, with the Plus 4 yielding coverage probabilities close to the desired, even for very small sample sizes.

SAS
Calculating the coverage probability for a given N and binomial p can be done in a single data step, summing the probability-weighted coverage indicators over the realized values of the random variate. Once this machinery is developed, we can call it repeatedly, using a macro, to find the results for different binomial p. We comment on the code internally.


%macro onep(n=,p=,alpha=.05);
data onep;
n = &n;
p = &p;
alpha = α
/* set up collectors of the weighted coverage indicators */
expcontrib_t = 0;
expcontrib_p4 = 0;
expcontrib_s = 0;
expcontrib_cp = 0;
/* loop through the possible observed successes x*/
do x = 0 to n;
  probobs = pdf('BINOM',x,p,n);  /* probability X=x */
  phat = x/n;
  zquant = quantile('NORMAl', 1 - alpha/2, 0, 1);
  p4 = (x+2)/(n + 4);

  /* calculate the half-width of the t and plus4 intervals */
  thalf = quantile('T', 1 - alpha/2,(n-1)) * sqrt(phat*(1-phat)/n);
  p4half = zquant * sqrt(p4*(1-p4)/(n+4));

  /* the score CI in R uses a Yates correction by default, and is 
     reproduced here */
  yates = min(0.5, abs(x - (n*p)));
  z22n = (zquant**2)/(2*n);
  yl = phat-yates/n;
  yu = phat+yates/n;
  slower = (yl + z22n - zquant * sqrt( (yl*(1-yl)/n) + z22n / (2*n) )) /
    (1 + 2 * z22n); 
  supper = (yu + z22n + zquant * sqrt( (yu*(1-yu)/n) + z22n / (2*n) )) /
    (1 + 2 * z22n); 

  /* cover = 1 if in the CI, 0 else */
  cover_t = ((phat - thalf) < p) and ((phat + thalf) > p);
  cover_p4 = ((p4 - p4half) < p) and ((p4 + p4half) > p);
  cover_s = (slower < p) and (supper > p);
  /* the Clopper-Pearson interval can be easily calculated on the fly */ 
  cover_cp = (quantile('BETA', alpha/2 ,x,n-x+1) < p) and 
    (quantile('BETA', 1 - alpha/2 ,x+1,n-x) > p); 

  /* cumulate the weighted probabilities */
  expcontrib_t = expcontrib_t + probobs * cover_t;
  expcontrib_p4 = expcontrib_p4 + probobs * cover_p4;
  expcontrib_s = expcontrib_s + probobs * cover_s;
  expcontrib_cp = expcontrib_cp + probobs * cover_cp;
  /* only save the last interation */
  if x = N then output;
  end;
run;
%mend onep;

The following macro calls the first one for a series of binomial p for a fixed N. Since the macro %do loop can only iterate through integers, we have to do a little division; the %sysevelf function will do this within the macro.


%macro repcicov(n=, lop=, hip=, byp=, alpha= .05);
/* need an empty data set to store the results */
data summ; set _null_; run;
%do stepp = %sysevalf(&lop / &byp, floor) %to %sysevalf(&hip / &byp,floor);
  /* note that the p sent to the %onep macro is a 
                           text string like "49 * .001" */
  %onep(n = &n, p = &stepp * &byp, alpha = &alpha);
  /* tack on the current results to the ones finished so far */
  /* this is a simple but inefficient way to add each binomial p into 
     the output data set */
  data summ; set summ onep; run;
%end;
%mend repcicov;

/* same parameters as in R */
%repcicov(n=50, lop = .01, hip = .3, byp = .001);

Finally, we can plot the results. One option shown here and not mentioned in the book are the mode=include option to the symbol statement, which allows the two distinct pieces of the T coverage to display correctly.


goptions reset=all;
legend1 label=none position=(bottom right inside)
        mode=share across=1 frame value = (h=2);
axis1 order = (0.85 to 1 by 0.05) minor=none 
   label = (a=90 h=2 "Coverage probability") value=(h=2);
axis2 order = (0 to 0.3 by 0.05) minor=none 
   label = (h=2 "True binomial p") value=(h=2);
symbol1 i = j v = none l =1 w=3 c=blue mode=include;
symbol2 i = j v = none l =1 w=3 c=red;
symbol3 i = j v = none l =1 w=3 c=lightgreen;
symbol4 i = j v = none l =1 w=3 c=black;
proc gplot data = summ;
plot (expcontrib_t expcontrib_p4  expcontrib_s  expcontrib_cp) * p 
  / overlay legend vaxis = axis1 haxis = axis2 vref = 0.95 legend = legend1;
label expcontrib_t = "T approximation" expcontrib_p4 = "P4 method"
      expcontrib_s = "Score method" expcontrib_cp = "Exact (CP)";
run; quit;

Citing R or SAS

2012-07-02T12:19:00.005-04:00

One of us recently read a colleague's first draft of a paper, in which she had written: "All analyses were done in R 2.14.0." We assume we're preaching to the converted here, when we say that the enormous amount of work that goes into R needs to be recognized as often as possible, and that R's creators deserve to reap some credit for their labors. In contrast to SAS, after all, most work on R is not compensated with a paycheck. As a reminder, the citation() function produces the correct citation for R in general and is good to use when citing R.

The project in question had used a negative binomial regression function from the MASS package, but colleague had omitted any reference to it. In this case, a citation would provide both credit to the authors and a useful guide to anyone wanting to replicate our approach. It would also allow readers to consider whether changes in the package might affect the results observed. A call to citation(package="MASS") will provide the preferred citation here. (Any package name can be inserted, or course, though some authors may not have provided a full citation.)

Similarly, while SAS authors are rarely identified by name and presumably get a salary from SAS, it's preferable to identify the version of the software and where it can be obtained. In medical research this is usually done by an in-text reference. For example: "Analyses were performed in SAS 9.3 (SAS Institute, Cary NC)."

For complex analyses, it is also best to mention the SAS procedure used. As with the R package, this can help readers plan similar analyses, and may inform interpretation.

So a multi-software analysis section might end with the following statement: Analyses were performed in R 2.14.2 [1] using the MASS package [2] glm.nb() function for negative binomial regression and in SAS 9.3 (SAS Institute, Cary NC) using the MCMC procedure for negative binomial mixture models." The references to [1] and [2] would be found using the citation() function.

An unrelated note about aggregators:We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.36: Levene's test for equal variances

2012-06-25T09:02:00.004-04:00

The assumption of equal variances among the groups in analysis of variance is an expression of the assumption of homoscedasticity for linear models more generally. For ANOVA, this assumption can be tested via Levene's test. The test is a function of the residuals and means within each group, though various modifications are used, including the Brown-Forsythe test. This uses the medians within group, rather than the mean, and is recommended when normality may be suspect.

We illustrate using the HELP data set available here, modeling depressive symptoms (assessed via CESD) as a function of abused substance.

SAS
In SAS, the tests are available as an option to the means statement in proc glm


data help;
set "C:\book\help.sas7bdat";
run;

proc glm data = help;
class substance;
model cesd = substance;
means substance / hovtest=levene(type=abs) hovtest=bf;
run;

The two requested tests are a version of Levene's test that is produced in R, below, and the aforementioned Brown-Forsythe test. The relevant results are shown below.


           Levene's Test for Homogeneity of CESD Variance
           ANOVA of Absolute Deviations from Group Means

                           Sum of        Mean
  Source           DF     Squares      Square    F Value    Pr > F

  SUBSTANCE         2       272.4       136.2       2.61    0.0747
  Error           450     23480.7     52.1793


     Brown and Forsythe's Test for Homogeneity of CESD Variance
          ANOVA of Absolute Deviations from Group Medians

                           Sum of        Mean
  Source           DF     Squares      Square    F Value    Pr > F

  SUBSTANCE         2       266.0       133.0       2.46    0.0864
  Error           450     24310.9     54.0243

There's some suggestion of a lack of homoscedasticity; it might be wise to consider methods robust to violations of this assumption.

R
In R, the test can be found in the levene.test() function in the lawstat package.


help = read.csv("http://www.math.smith.edu/r/data/help.csv")
library(lawstat)
with(help, levene.test(cesd, as.factor(substance), location="mean"))

 classical Levene's test based on the absolute deviations from the mean 
          ( none not applied because the location is not set to median )

data:  cesd 
Test Statistic = 2.6099, p-value = 0.07465

with(help, levene.test(cesd, as.factor(substance),location="median"))

 modified robust Brown-Forsythe Levene-type test based on the absolute 
           deviations from the median

data:  cesd 
Test Statistic = 2.462, p-value = 0.08641

Example 9.35: Discrete randomization and formatted output

2012-06-18T10:11:00.005-04:00

A colleague asked for help with randomly choosing a kid within a family. This is for a trial in which families are recruited at well-child visits, but in each family only one of the children having a well-child visit that day can be in the study. The idea is that after recruiting the family, the research assistant needs to choose one child, but if they make that choice themselves, the children are unlikely to be representative. Instead, we'll allow them to make a random decision through an easily used slip that can be put into sealed envelopes. The envisioned process is that the RA will recruit the family, determine the number of eligible children, then open the envelope to find out which child was randomly selected.

One thought here would be to generate separate stacks of envelopes for each given family size, and have the research assistant open an envelope from the appropriate stack. However, this could be logistically challenging, especially since the RAs will spend weeks away from the home office. Instead, we'll include all plausible family sizes on each slip of paper. It seems unlikely that more than 5 children in a family will have well-child visits on the same day.

SAS
We'll use the SAS example to demonstrate using SAS macros to write SAS code, as well as showing a plausible use for SAS formats (section 1.4.12) and making use of proc print.


/* the following macro will write out equal probabilities for selecting 
each integer between 1 and the argument, in the format needed for the 
rand function.  E.g., if the argument is 3, 
it will write out
1/3,1/3,1/3
*/

%macro tbls(n);
%do i = 1 %to &n;
1/&n %if &i < &n %then ,
%end;
%mend tbls;

/* then we can use the %tbls macro to create the randomization
via rand("TABLE") (section 1.10.4). */ 
data kids;
do family = 1 to 10000;
  nkids = 2; chosen = rand("TABLE",%tbls(2)); output;
  nkids = 3; chosen = rand("TABLE",%tbls(3)); output;
  nkids = 4; chosen = rand("TABLE",%tbls(4)); output;
  nkids = 5; chosen = rand("TABLE",%tbls(5)); output;
end;
run;

/* check randomization */
proc freq data = kids;
table nkids * chosen / nocol nopercent;
run; 

   nkids     chosen

   Frequency|
   Row Pct  |       1|       2|       3|       4|       5|  Total
   ---------+--------+--------+--------+--------+--------+
          2 |  50256 |  49744 |      0 |      0 |      0 | 100000
            |  50.26 |  49.74 |   0.00 |   0.00 |   0.00 |
   ---------+--------+--------+--------+--------+--------+
          3 |  33429 |  33292 |  33279 |      0 |      0 | 100000
            |  33.43 |  33.29 |  33.28 |   0.00 |   0.00 |
   ---------+--------+--------+--------+--------+--------+
          4 |  25039 |  24839 |  25245 |  24877 |      0 | 100000
            |  25.04 |  24.84 |  25.25 |  24.88 |   0.00 |
   ---------+--------+--------+--------+--------+--------+
          5 |  19930 |  20074 |  20188 |  20036 |  19772 | 100000
            |  19.93 |  20.07 |  20.19 |  20.04 |  19.77 |
   ---------+--------+--------+--------+--------+--------+
   Total      128654   127949    78712    44913    19772   400000

Looks pretty good. Now we need to make the output usable to the research assistants, by formatting the results into English. We'll use the same format for each number of kids. This saves some keystrokes now, but may possibly cause the RAs some confusion-- it means that we might refer to the "4th oldest" of 4 children, rather than the "youngest". We could fix this using a different format for each number of children, analogous to the R version below.


proc format;
value chosen
1 = "oldest"
2 = '2nd oldest'
3 = '3rd oldest'
4 = '4th oldest'
5 = '5th oldest';
run;
 
/* now, make a text variable the concatenates (section 1.4.5) the variables 
and some explanatory text */
data k2;
set kids;
if nkids eq 2 then
  t1 = "If there are " || strip(nkids) ||" children then choose the " ||
       strip(put(chosen,chosen.)) || " child.";
else
  t1 = "             " || strip(nkids) ||" ________________________ " ||
       strip(put(chosen,chosen.));
run;

/* then we print.  Notice the options to print in plain text, shorten the 
page length and width, and remove the date and page number from the SAS output, as
well as in the proc print statement to remove the observation number and
show the line number, with a few other tricks */
options nonumber nodate ps = 60 ls = 68;
OPTIONS FORMCHAR="|----|+|---+=|-/\<>*";
proc print data = k2 (obs = 3) noobs label sumlabel;
by family;
var t1;
label t1 = '00'x family = "Envelope";
run;

---------------------------- Envelope=1 ----------------------------



     If there are 2 children then choose the 2nd oldest child.
                  3 ________________________ 3rd oldest
                  4 ________________________ 4th oldest
                  5 ________________________ 5th oldest


---------------------------- Envelope=2 ----------------------------



     If there are 2 children then choose the 2nd oldest child.
                  3 ________________________ oldest
                  4 ________________________ oldest
                  5 ________________________ 3rd oldest


---------------------------- Envelope=3 ----------------------------



     If there are 2 children then choose the 2nd oldest child.
                  3 ________________________ 2nd oldest
                  4 ________________________ 3rd oldest
                  5 ________________________ 2nd oldest

R
For R, we leave some trial code in place, to demonstrate how one might discover, test, and build R code in this setting. Most results have been omitted.


sample(5, size = 1)   
# choose a (discrete uniform) random integer between 1 and 5

apply(matrix(2:5),1,sample,size=1)   
# choose a random integer between 1 and 2, then between 1 and 3, etc., 
# using apply() to repeat the call to sample() with different maximum number
# apply() needs a matrix or array input
# result of this is the raw data needed for one family

replicate(3,apply(matrix(2:5),1,sample,size=1))
# replicate() is in the apply() family and just repeats the 
# function n times

     [,1] [,2] [,3]
[1,]    2    1    2
[2,]    2    1    2
[3,]    2    2    2
[4,]    3    5    4

Now we have the raw data for the envelopes. Before formatting it for printing, let's check it to make sure it works correctly.


test=replicate(100000, apply(matrix(2:5), 1, sample, size=1))
apply(test, 1, summary)
        [,1] [,2]  [,3]  [,4]
Min.     1.0    1 1.000 1.000
1st Qu.  1.0    1 1.000 2.000
Median   1.0    2 2.000 3.000
Mean     1.5    2 2.492 3.003
3rd Qu.  2.0    3 3.000 4.000
Max.     2.0    3 4.000 5.000
# this is not so helpful-- need the count or percent for each number
# this would be the default if the data were factors, but they aren't
# check to see if we can trick summary() into treating these integers
# as if they were factors
methods(summary)
# yes, there's a summary() method for factors-- let's apply it
# there's also apply(test,1,table) which might be better, if you remember it
apply(test, 1, summary.factor)
[[1]]
    1     2 
50025 49975 

[[2]]
    1     2     3 
33329 33366 33305 

[[3]]
    1     2     3     4 
25231 25134 24849 24786 

[[4]]
    1     2     3     4     5 
19836 20068 20065 20022 20009 
# apply(test,1,table) will give similar results, if you remember it

Well, that's not too pretty, but it's clear that the randomization is working. Now it's time to work on formatting the output.


mylist=replicate(5, apply(matrix(2:5), 1, sample, size=1))
# brief example data set

# We'll need to use some formatted values (section 1.14.12), as in SAS. 
# Here, we'll make new value labels for each number of children,
# which will make the output easier to read.  We add in an envelope 
# number and wrap it all into a data frame.
df = data.frame(envelope = 1:5,
   twokids=factor(mylist[1,],1:2,labels=c("youngest","oldest")),
  threekids=factor(mylist[2,],1:3,labels=c("youngest", "middle", "oldest")),
  fourkids=factor(mylist[3,],1:4,labels=c("youngest", "second youngest", 
      "second oldest", "oldest")),
  fivekids=factor(mylist[4,],1:5,labels=c("youngest", "second youngest", 
      "middle", "second oldest", "oldest"))
)

# now we need a function to take a row of the data frame and make a single slip
# the paste() function (section 1.4.5) puts together the fixed and variable 
# content of each row, while the cat() function will print it without quotes
slip = function(kidvec) {
  cat(paste("------------- Envelope", kidvec[1], "------------------"))
  cat(paste("\nIf there are", 2:5, " children, select the", kidvec[2:5],"child"))
  cat("\n \n \n")
}

# test it on one row
slip(df[1,])

# looks good-- now we can apply() it to each row of the data frame
apply(df, 1, slip)

------------- Envelope 1 ------------------
If there are 2  children, select the youngest child 
If there are 3  children, select the youngest child 
If there are 4  children, select the second youngest child 
If there are 5  children, select the youngest child
 
 
------------- Envelope 2 ------------------
If there are 2  children, select the youngest child 
If there are 3  children, select the youngest child 
If there are 4  children, select the second oldest child 
If there are 5  children, select the middle child


------------- Envelope 3 ------------------
If there are 2  children, select the youngest child 
If there are 3  children, select the youngest child 
If there are 4  children, select the youngest child 
If there are 5  children, select the second youngest child

# and so forth

# finally, we can save the result in a file with
# capture.output()
capture.output(apply(df,1,slip), file="testslip.txt")

An unrelated note about aggregators:
We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.34: Bland-Altman type plot

2012-06-05T15:59:00.003-04:00

The Bland-Altman plot is a visual aid for assessing differences between two ways of measuring something. For example, one might compare two scales this way, or two devices for measuring particulate matter.

The plot simply displays the difference between the measures against their average. Rather than a statistical test, it is intended to demonstrate both typical differences between the measures and any patterns such differences may take. The utility of the plot, as compared with linear regression or sample correlation is that the plot is not affected by the range, while the sample correlation will typically increase with the range. In contrast, linear regression shows the strength of the linear association but not how closely the two measures agree. The Bland-Altman plot allows the user to focus on differences between the measures, perhaps focusing on the clinical relevance of these differences.

A peer reviewer recently asked a colleague to consider a Bland-Altman plot for two methods of assessing fatness: the familiar BMI (kg/m^2) and the actual fat mass measured by a sophisticated DXA machine. These are obviously not measures of the same thing, so a Bland-Altman plot is not exactly appropriate. But since the BMI is so simple to calculate and the DXA machine is so expensive, it would be nice if the BMI could be substituted for DXA fat mass.

For this purpose, we'll generate a modified Bland-Altman plot in which each measure is first standardized to have mean 0 and standard deviation 1. The resulting plot should be assessed for pattern as usual, but typical differences must be considered on the standardized scale-- that is, differences of a unit should be considered large, and good agreement might require typical differences of 0.2 or less.

SAS
Since this is a job we might want to repeat, we'll build a SAS macro to do it. This will also demonstrate some useful features. The macro accepts a data set name and the names of two variables as input. We'll comment on interesting features in code comments. If you're an R coder, note that SAS macro variables are merely text, not objects. We have to manually assign "values" (i.e., numbers represented as text strings) to newly created macro variables.


%macro baplot(datain=,x=x,y=y);

/* proc standard  standardizes the variables and saves the results in the
   same variable names in the output data set.  This means we can continue 
   using the input variable names throughout. */
proc standard data = &datain out=ba_std mean=0 std=1;
var &x &y;
run;

/* calculate differences and averages */
data ba;
set ba_std;
bamean = (&x + &y)/2;;
badiff = &y-&x;
run;

ods output summary=basumm;
ods select none;
proc means data = ba mean std;
var badiff;
run;
ods select all;

/* In the following, we take values calculated from a data set for the 
   confidence limits and store them in macro variables.  That's the 
   only way to use them later in code.
   The syntax is: call symput('varname', value);
   Note that 'bias' is purely nominal, as the standardization means that 
   the mean difference is 0. */
data lines;
set basumm;
call symput('bias',badiff_mean);
call symput('hici',badiff_mean+(1.96 * badiff_stddev));
call symput('loci',badiff_mean-(1.96 * badiff_stddev));
run;

/* We use the macro variables just created in the vref= option below;
   vref draws reference line(s) on the vertical axis.  lvref specifies
   a line type. */
symbol1 i = none v = dot h = .5;
title "Bland-Altman type plot of &x and &y";
title2 "&x and &y standardized";
proc gplot data=ba;
plot badiff * bamean / vref = &bias &hici &loci lvref=3;
label badiff = "difference" bamean="mean";
run;
%mend baplot;

Here is a fake sample data set, with the plot resulting from the macro shown above. An analysis would suggest that despite the correlation of 0.59 and p-value for the linear association < .0001, that these two measures don't agree too well.


data fake;
do i = 1 to 50;
/* the "42" in the code below sets the seed for the pseudo-RNG
   for this and later calls.  See section 1.10.9. */
  x = normal(42);
  y = x + normal(0);
  output;
  end;
run;

%baplot(datain=fake, x=x, y=y);

R
Paralleling SAS, we'll write a small function to draw the plot, annotating within to highlight some details. If you're primarily a SAS coder, note the syntax needed to find the name of an object submitted to a function. In contrast, assigning values to new objects created with the function is entirely natural. The resulting plot is shown below.


# set seed, for replicability
set.seed(42)
x = rnorm(50)
y = x + rnorm(50) 

baplot = function(x,y){
  xstd = (x - mean(x))/sd(x)
  ystd = (y - mean(y))/sd(y)
  
  bamean = (xstd+ystd)/2
  badiff = (ystd-xstd)
  
  plot(badiff~bamean, pch=20, xlab="mean", ylab="difference")
# in the following, the deparse(substitute(varname)) is what retrieves the
# name of the argument as data
  title(main=paste("Bland-Altman plot of x and y\n",
    deparse(substitute(x)), "and", deparse(substitute(y)),
    "standardized"), adj=".5")
#construct the reference lines on the fly: no need to save the values in new 
# variable names
  abline(h = c(mean(badiff), mean(badiff)+1.96 * sd(badiff),
    mean(badiff)-1.96 * sd(badiff)), lty=2)
} 

baplot(x,y)

Example 9.33: Multiple imputation, rounding, and bias

2012-05-29T09:18:00.012-04:00

Nick has a paper in the American Statistician warning about bias in multiple imputation arising from rounding data imputed under a normal assumption. One example where you might run afoul of this is if the data are truly dichotomous or count variables, but you model it as normal (either because your software is unable to model dichotomous values directly or because you prefer the theoretical soundness of multivariate normal imputation to, e.g., chained equations). In such cases, one might impute assuming normality, then round the imputed values to plausible integers. The paper shows theoretically the bias that can result if this process is pursued, and also that allowing the "implausible values" will eliminate the bias. (Of course, modeling the missing variable using a logistic regression model will be most appropriate here).

In another paper, Nick and Stuart Lipsitz (TAS 2001) comment that the method of predictive mean matching (PMM) "ensures that imputed values are plausible, and may be more appropriate if the normality assumption is violated." Briefly, the PMM method predicts a value from a model for both missing and observed values. The imputation for a subject with a missing value is the observed value of the subject with the nearest predicted value (or random draw of observed values from among the subjects with the nearest predicted values).

How does this play out in practice? Can the PMM method overcome the theoretical rounding bias while still generating only plausible imputed values?

SAS
We begin by simulating dichotomous data, choosing the value of p (probability of 1) = .25, a value with a large absolute bias, according to the paper. We set values to missing with probability 0.5, using a MCAR mechanism. Then we use proc mi (section 6.5, example 9.4) to impute the missing values, assuming normality. The mean and standard error of the mean of y are calculated in proc summary (section 2.1.1) and combined in proc mianalyze. Then the values are rounded manually and the analysis repeated. Next, we impute separately with PMM. Finally, we impute again with a logistic imputation. We use 5 imputations throughout, though 50 would likely be preferable.

Note that a Poisson regression imputation is not yet available for proc mi, so that the exercise is not wholly academic--if you needed to impute count values, you'd have to choose among implausible values, rounding, and PMM. Also note our use of the fcs imputation method, though it is not needed here with an obviously monotone missingness pattern. Finally, note that proc mi here requires at least two variables, for no reason we know of. We generate a normally-distributed and uncorrelated covariate.


data testpmm;
do i = 1 to 5000;
  x = normal(0);
  y = rand('BINOMIAL', .25, 1);
  missprob = ranuni(0);
  if missprob le .5 then y = .;
  output;
  end;
run;

title "Normal imputation";
proc mi data=testpmm out=normal nimpute=5;
var x y;
fcs reg;
run;

title2 "Implausible values";
proc summary data = normal mean stderr;
by _imputation_;
var y;
output out=outnormal mean=meany stderr=stderry;
run;

proc mianalyze data = outnormal;
modeleffects meany;
stderr stderry;
run;


title2 "Rounded";
/* make the rounded data */
data normalrnd;
set normal;
if y lt .5 then y=0;
else y=1;
run;

proc summary data = normalrnd mean stderr;
by _imputation_;
var y;
output out=outnormalrnd mean=meany stderr=stderry;
run;

proc mianalyze data = outnormalrnd;
modeleffects meany;
stderr stderry;
run;

title "regpmm imputation";
proc mi data=testpmm out=pmm nimpute=5;
var x y;
fcs regpmm;
run;
...

title "logistic imputation";
proc mi data=testpmm out=logistic nimpute=5;
class y;
var x y;
fcs logistic;
run;
...

We omit the summary and mianalyze procedures for the latter imputations. Ordinarily, it would be easiest to do this kind of repetitive task with a macro, but we leave it in open code here for legibility.
The results are shown below


                          Normal imputation-- Implausible values

          Parameter        Estimate      Std Error    95% Confidence Limits  
          meany            0.249105       0.008634     0.230849     0.267362
    
                                Normal imputation-- Rounded 
          meany            0.265280       0.006408     0.252710     0.277850
    
                                       regpmm imputation      
          meany            0.246320       0.006642     0.233204     0.259436  

                                      logistic imputation     
          meany            0.255120       0.008428     0.237449     0.272791

As theory suggests, rounding the normally imputed values leads to bias, while using the normal imputations does not (though it results in implausible values). Nether PMM imputation nor direct logistic imputation appear to be biased.

R
We will use the mice package written by Stef van Buuren, one of the key developers of chained imputation. Stef also has a new book describing the package and demonstrating its use in many applied examples. We use 5 imputations throughout, though 50 would likely be preferable.

We begin by creating the data. Note that mice(), like proc mi, requires at least two columns of data. To do the logistic regression imputation, mice() wants the missing data to be a factor, so we make a copy of the data as a data frame object as well.


library(mice)
n = 5000  # number of observations
m = 5   # number of imputations (should be 25-50 in practice)
x = rnorm(n)
y = rbinom(n, 1, .25)   # interesting point according to Horton and Lipsitz (TAS 2004)
unif = runif(n)
y[unif < .5] = NA    # make half of the Y's be missing
ds = cbind(y, x)
ds2 = data.frame(factor(y), x)

The mice package works analogously to proc mi/proc mianalyze. The mice() function performs the imputation, while the pool() function summarizes the results across the completed data sets. The method option to mice() specifies an imputation method for each column in the input object. Here we fit the simplest linear regression model (intercept only).


# normal model with implausible values
impnorm = mice(ds, method="norm", m=m)
summary(pool(with(impnorm, lm(y ~ 1))))

Rounding could be done by tampering with the mids-type object that mice() produces, but there is a more direct way to do this through the post= option. It accepts text strings with R commands that will be applied to the imputed values. Here we use the ifelse() function to make the normal values equal to 0 or 1. The code for the predictive mean matching and logistic regression follow.


impnormround = mice(ds, method="norm", m=m, 
   post= c("imp[[j]][,i] = ifelse(imp[[j]][,i] < .5, 0, 1)",""))

imppmm = mice(ds, method="pmm", m=m)

implog = mice(ds2, method="logreg", m=m)

The results of summary(pool()) calls are shown below..


> summary(pool(with(impnorm, lm(y ~ 1))))
                 est          se     lo 95     hi 95 
(Intercept) 0.272912 0.007008458 0.2589915 0.2868325
> summary(pool(with(impnormround, lm(y ~ 1))))
                est         se     lo 95     hi 95 
(Intercept) 0.28544 0.00854905 0.2676263 0.3032537
> summary(pool(with(imppmm, lm(y ~ 1))))
                 est         se     lo 95     hi 95
(Intercept) 0.277636 0.03180604 0.2145564 0.3407156
> summary(pool(with(implog, lm(y ~ 1))))
                  est         se     lo 95     hi 95 
(Intercept) 0.2652899 0.00879988 0.2480342 0.2825457

The message on bias is similar, though there is some hint of trouble in the CI for the PMM method (it seems to have a bias towards 0.5). The default option of 3 donors may be too few (this can be tweaked by use of the donors = NUMBER option).

Example 9.32: Multiple testing simulation

2012-05-21T10:36:00.024-04:00

In examples 9.30 and 9.31 we explored corrections for multiple testing and then extracting p-values adjusted by the Benjamini and Hochberg (or FDR) procedure. In this post we'll develop a simulation to explore the impact of "strong" and "weak" control of the family-wise error rate offered in multiple comparison corrections. Loosely put, weak control procedures may fail when some of the null hypotheses are actually false, in that the remaining (true) nulls may be rejected more than the nominal proportion of times.

For our simulation, we'll develop flexible code to generate some p-values from false nulls and others from true nulls. We'll assume that the true nulls have p-values distributed uniform (0,1); the false nulls will have p-values distributed uniform with a user-determined maximum. We'll also allow the number of tests overall and the number of false nulls to be set.

SAS
In SAS, a macro does the job. It accepts the user parameters described above, then generates false and true nulls for each desired simulation. With the data created, we can use proc multtest to apply the FDR procedure, with the ODS system saving the results. Note how the by statement allows us to replicate the analysis for each simulated set of p-values without creating a separate data set for each one. (Also note that we do not use proc sort before that by statement-- this can be risky, but works fine here.)


%macro fdr(nsims=1, ntests = 20, nfalse=10, howfalse=.01);
ods select none;
data test;
do sim = 1 to &nsims;
  do i = 1 to &ntests;
    raw_p = uniform(0) * 
      ( ((i le &nfalse) * &howfalse ) + ((i gt &nfalse) * 1 ) );
    output;
  end;
end;
run;

ods output pvalues = __pv;
proc multtest inpvalues=test fdr;
by sim;
run;

With the results in hand, (still within the macro) we need to do some massaging to make the results usable. First we'll recode the rejections (assuming a 0.05 alpha level) so that non-rejections are 0 and rejections are 1/number of tests. That way we can just sum across the results to get the proportion of rejections. Next, we transform the data to get each simulation in a row (section 1.5.4). (The data output from proc multtest has nsims*ntests rows. After transposing, there are nsims rows.) Finally, we can sum across the rows to get the proportion of tests rejected in each simulated family of tests. The results are shown in a table made with proc freq.


data __pv1;
set __pv;
if falsediscoveryrate lt 0.05 then fdrprop = 1/&ntests;
else fdrprop =0;
run;

proc transpose data = __pv1 (keep =sim fdrprop) out = pvals_a;
by sim; run;

data pvals;
set pvals_a;
prop = sum(of col1 - col&ntests);
run;
ods select all;

proc freq data = pvals; tables prop; run;
%mend fdr;

%fdr(nsims = 1000, ntests = 20, nfalse = 10, howfalse=.001);

                                      Cumulative    Cumulative
     prop    Frequency     Percent     Frequency      Percent
     ---------------------------------------------------------
      0.5         758       75.80           758        75.80
     0.55         210       21.00           968        96.80
      0.6          27        2.70           995        99.50
     0.65           5        0.50          1000       100.00

So true nulls were rejected 24% of the time, which seems like a lot. Multiple comparison procedures with "strong" control of the familywise error rate will reject them only 5% of the time. Building this simulation as a macro facilitates exploring the effects of the multiple comparison procedures in a variety of settings.

R
As in example 9.31, the R code is rather simpler, though perhaps a bit opaque. To make the p-values, we make them first for all of tests with the false, then for all of the tests with the true nulls. The matrix function reads these in by column, by default, meaning that the first nfalse columns get the nsims*nfalse observations. The apply function generates the FDR p-values for each row of the data set. The t() function just transposes the resulting matrix so that we get back a row for each simulation. As in the SAS version, we'll count each rejection as 1/ntests, and non-rejections as 0; we do this with the ifelse() statement. Then we sum across the simulations with another call to apply() and show the results with a simple table.


checkfdr = function(nsims=1, ntests=100, nfalse=0, howfalse=0.001) {
  raw_p = matrix(c(runif(nfalse * nsims) * howfalse, 
                   runif((ntests-nfalse) * nsims)), nrow=nsims)
  fdr = t(apply(raw_p, 1, p.adjust, "fdr"))
  reject = ifelse(fdr<.05, 1/ntests,0)
  prop = apply(reject, 1, sum)
  prop.table(table(prop)) 
}

> checkfdr(nsims=1000, ntests=20, nfalse=10, howfalse=.001)
prop
  0.5  0.55   0.6  0.65 
0.755 0.210 0.032 0.003

The results are reassuringly similar to those from SAS. In this R code, it's particularly simple to try a different test-- just replace "fdr" in the p.adjust() call. Here's the result with the Hochberg test, which has strong control.


checkhoch = function(nsims=1, ntests=100, nfalse=0, howfalse=0.001) {
   pvals = matrix(c(runif(nfalse * nsims) * howfalse, 
                    runif((ntests-nfalse) * nsims)), nrow=nsims)
   hochberg = t(apply(pvals, 1, p.adjust,"hochberg"))
   reject = ifelse(hochberg<.05,1/ntests,0)
   prop = apply(reject, 1, sum)
   prop.table(table(prop)) 
}
 
> checkhoch(nsims=1000, ntests=20, nfalse=10, howfalse=.001)
prop
  0.5  0.55   0.6 
0.951 0.046 0.003

With this procedure one or more of the true nulls is rejected an appropriate 4.9% of the time. For the most part, we feel more comfortable using multiple testing procedures with "strong control".

An unrelated note about aggregators
We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.31: Exploring multiple testing procedures

2012-05-14T09:51:00.011-04:00

In example 9.30 we explored the effects of adjusting for multiple testing using the Bonferroni and Benjamini-Hochberg (or false discovery rate, FDR) procedures. At the time we claimed that it would probably be inappropriate to extract the adjusted p-values from the FDR method from their context. In this entry we attempt to explain our misgivings about this practice.

The FDR procedure is described in Benjamini and Hochberg (JRSSB, 1995) as a "step-down" procedure. Put simply, the procedure has the following steps:


0. Choose the familywise alpha
1. Rank order the unadjusted p-values
2. Beginning with the Mth of the ordered p-values p(m), 
2a.    if p(m) < alpha*(m/M), then reject all tests 1 ... m, 
2b.    if not, m = m-1
3. Repeat steps 2a and 2b until the condition is met 
               or p(1) > alpha/M

where M is the number of tests. The "adjusted p-value" based on this procedure is the smallest familywise alpha under which the current test would have been rejected. To calculate this, we can modify the routine above:


1. Rank order the unadjusted p-values
2. For ordered p-values p(m) M to 1, 
2a.    candidate ap(m) = p(m) *(M/m) 
2b.    if candidate ap(m) > ap(m+1) then ap(m) = ap(m+1)
2c.    else ap(m) = candidate ap(m)

where ap(m) refers to the adjusted p-value corresponding to the mth ordered unadjusted p-value. It's interesting to note that the adjusted p-value for the Mth ordered test is the same as the unadjusted p-value, while the candidate adjusted p-value for the smallest test is the Bonferroni adjusted p-value. The primary difficulty with taking these p-values (as opposed to the test results) out of context is captured in steps 2b and 2c. They imply that the p-value for a given test may be lowered by other observed p-values in the family of tests. It's also true that the adjusted p-value depends on the number of tests included in the family, but this seems somewhat less troubling.

To examine the impact of the procedure on the adjusted p-values for the individual tests, we'll compare the candidate ap(m) from step 2a against the actual ap(m). Our sense is that to the degree these are different, the adjusted p-value should not be extracted from the context of the observed family of tests.

SAS
Our SAS code relies heavily on the array statement (section 1.11.5). We loop through the p-values from largest to smallest, calculating the candidate fdr p-value as above, before arriving at the final adjusted p-value. To compare the values conveniently, we make a new data set with two copies of the original data set, renaming first the candidate and then the adjusted p-values to have the same names. The in = data set option creates a temporary variable which identifies which data set an observation was read from; here it denotes which version of the same data set (and which set of p-values) was used.


data fdr;
array pvals [10] pval1 - pval10 
     (.001 .001 .001 .001 .001 .03 .035 .04 .05 .05);
array cfdrpvals [10] cfdr1 - cfdr10;
array fdrpvals [10] fdr1 - fdr10;
fdrpvals[10] = pvals[10];
do i = 9 to 1 by -1;
  cfdrpvals[i] = pvals[i] * 10/i;
  if cfdrpvals[i] > fdrpvals[i+1] then fdrpvals[i] = fdrpvals[i+1];
  else fdrpvals[i] = cfdrpvals[i];
  end;
run;

data compare;
set fdr (in = cfdr rename=(cfdr1=c1 cfdr2=c2 cfdr3=c3 cfdr4=c4 
           cfdr5=c5 cfdr6=c6 cfdr7=c7 cfdr8=c8 cfdr9=c9)) 
    fdr (in = fdr rename=(fdr1=c1 fdr2=c2 fdr3=c3 fdr4=c4 fdr5=c5 
           fdr6=c6 fdr7=c7 fdr8=c8 fdr9=c9));
if cfdr then adjustment = "Candidate fdr";
if fdr then adjustment = "Final fdr";
run;

proc print data = compare; var adjustment c1-c9; run;

adjustment       c1    c2    c3     c4    c5     c6    c7    c8    c9

Candidate fdr   0.010  .005  .0033  .0025  .002  .05   .05   .05   .055
Final fdr       0.002  .002  .0020  .0020  .002  .05   .05   .05   .050

(We omit the last p-value because the adjustment does not affect it.) The result shows that for many of the tests in this family, a substantially smaller p-value is obtained with the final FDR p-value than the candidate. To this degree, the FDR p-value is dependent on the observed values of the p-values in the tests in the family, and ought not to be removed from the context of these other tests. We would recommend caution in displaying the FDR p-values in such settings, given readers' propensity to use them as if they were ordinary p-values, safely adjusted for multiple testing.

R
Comparison of the R and SAS code may make SAS programmers weep. The candidate values are easily calculated, and can be presented with the final p-values in one step using the p.adjust() function. Three lines of code, albeit incorporating multiple functions in each line. (And it could sensibly be done in two, calculating the candidate p-values within the rbind() function call.) Note especially the line calculating the candidate p-values, in which vectorization allows a for loop to be avoided in a very natural fashion.


fakeps = c(rep(.2, 5), 6, 7, 8, 10, 10)/200
cfdr = fakeps * 10/(1:10)
rbind(cfdr, fdr=p.adjust(fakeps, "fdr"))[,1:9]

      [,1]  [,2]   [,3]   [,4]  [,5] [,6] [,7] [,8]   [,9] [,10]
cfdr 0.010 0.005 0.0033 0.0025 0.002 0.05 0.05 0.05 0.0556  0.05
fdr  0.002 0.002 0.0020 0.0020 0.002 0.05 0.05 0.05 0.0500  0.05

An unrelated note about aggregatorsWe love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers and PROC-X with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.30: addressing multiple comparisons

2012-05-07T09:55:00.015-04:00

We've been more sensitive to accounting for multiple comparisons recently, in part due to work that Nick and colleagues published on the topic.

In this entry, we consider results from a randomized trial (Kypri et al., 2009) to reduce problem drinking in Australian university students.
Seven outcomes were pre-specified: three designated as primary and four as secondary. No adjustment for multiple comparisons was undertaken. The p-values were given as 0.001, 0.001 for the primary outcomes and 0.02 and .001, .22, .59 and .87 for the secondary outcomes.
In this entry, we detail how to adjust for multiplicity using R and SAS.

R

The p.adjust() function in R calculates a variety of different approaches for multiplicity adjustments given a vector of p-values. These include the Bonferroni procedure (where the alpha is divided by the number of tests or equivalently the p-value is multiplied by that number, and truncated back to 1 if the result is not a probability). Other, less conservative corrections are also included (these are Holm (1979), Hochberg (1988), Hommel (1988), Benjamini and Hochberg (1995) and Benjamini and Yekutieli (2001)). The first four methods provide strong control for the family-wise error rate and all dominate the Bonferroni procedure. Here we compare the results from the unadjusted, Benjamini and Hochberg method="BH" and Bonferroni procedure for the Kypri et al. study.


pvals = c(.001, .001, .001, .02, .22, .59, .87)
BONF = p.adjust(pvals, "bonferroni")
BH = p.adjust(pvals, "BH")
res = cbind(pvals, BH=round(BH, 3), BONF=round(BONF, 3))

This yields the following results:


    pvals    BH  BONF
[1,] 0.001 0.002 0.007
[2,] 0.001 0.002 0.007
[3,] 0.001 0.002 0.007
[4,] 0.020 0.035 0.140
[5,] 0.220 0.308 1.000
[6,] 0.590 0.688 1.000
[7,] 0.870 0.870 1.000

The only substantive difference between the three sets of unadjusted and adjusted p-values is seen for the 4th most significant outcome, which remains statistically significant at the alpha=0.05 level for all but the Bonferroni procedure.

It is straightforward to graphically display these results (as seen above):


matplot(res, ylab="p-values", xlab="sorted outcomes")
abline(h=0.05, lty=2)
matlines(res)
legend(1, .9, legend=c("Bonferroni", "Benjamini-Hochberg", "Unadjusted"), 
  col=c(3, 2, 1), lty=c(3, 2, 1), cex=0.7)

It bears mentioning here that the Benjamini-Hochberg procedure really only make sense in the gestalt. That is, it would probably be incorrect to take the adjusted p-values from above and remove them from the context of the 7 tests performed here. The correct use (as with all tests) is to pre-specify the alpha level, and reject tests with p-values that are smaller. What p.adjust() reports is the smallest family-wise alpha error under which each of the tests would result in a rejection of the null hypothesis. But the nature of the Benjamini-Hochberg procedure is that this value may well depend on the other observed p-values. We will explore this further in a later entry.

SAS
The multtest procedure will perform a number of multiple testing procedures. It works with raw data for ANOVA models, and can also accept a list of p-values as shown here. (Note that "FDR" (false discovery rate) is the name that Benjamini and Hochberg give to their procedure and that this nomenclature is used by SAS.) Various other procedures can do some adjustment through, e.g., the estimate statement, but multtest is the most flexible. A plot similar to that created in R is shown below.


data a;
   input Test$ Raw_P @@;
   datalines;
test01  0.001    test02  0.001    test03  0.001
test04  0.02    test05  0.22    test06  0.59
test07  0.87
;

proc multtest inpvalues=a bon fdr plots=adjusted(unpack);
run;
                                                     False
                                                 Discovery
            Test           Raw    Bonferroni          Rate

               1        0.0010        0.0070        0.0023
               2        0.0010        0.0070        0.0023
               3        0.0010        0.0070        0.0023
               4        0.0200        0.1400        0.0350
               5        0.2200        1.0000        0.3080
               6        0.5900        1.0000        0.6883
               7        0.8700        1.0000        0.8700

An unrelated note about aggregators:We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers and PROC-X with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.29: the perils of for loops

2012-04-30T10:34:00.000-04:00

A recent exchange on the R-sig-teaching list featured a discussion of how best to teach new students R. The initial post included an exercise to write a function, that given a n, will draw n rows of a triangle made up of "*", noting that for a beginner, this may require two for loops. For example, in pseudo-code:


for i = 1 to n
  for j = 1 to i
    print "*"

Unfortunately, as several folks (including Richard M. Heiberger and R. Michael Weylandt) noted, for loops in general are not the best way to take full advantage of R. In this entry, we review two solutions they proposed which fit within the R philosophy.

Richard's solution uses the outer() function to generate a 5x5 matrix of logical values indicating whether the column number is bigger than the row number. Next the ifelse() function is used to replace TRUE with *.


> ifelse(outer(1:5, 1:5, `>=`), "*", " ")  
     [,1] [,2] [,3] [,4] [,5]
[1,] "*"  " "  " "  " "  " " 
[2,] "*"  "*"  " "  " "  " " 
[3,] "*"  "*"  "*"  " "  " " 
[4,] "*"  "*"  "*"  "*"  " " 
[5,] "*"  "*"  "*"  "*"  "*"

Michael's solution uses the lapply() function to call a function repeatedly for different values of n. This returns a list rather than a matrix, but accomplishes the same task.


> lapply(1:5, function(x) cat(rep("*", x), "\n"))
* 
* * 
* * * 
* * * * 
* * * * *

While this exercise is of little practical value, it does illustrate some important points, and provides a far more efficient as well as elegant way of accomplishing the tasks. For those interested in more, another resource is the R Inferno project of Patric Burns.

SAS
We demonstrate a SAS data step solution mainly to call out some useful features and cautions. In all likelihood a proc iml matrix-based solution would be more elegant;


data test;
array star [5] $ star1 - star5;
do i = 1 to 5;
  star[i] = "*";
  output;
  end;
run;

proc print noobs; var star1 - star5; run;

             star1    star2    star3    star4    star5

               *
               *        *
               *        *        *
               *        *        *        *
               *        *        *        *        *

In particular, note the $ in the array statement, which allows the variables to contain characters; by default variables created by an array statement are numeric. In addition, note the reference to a sequentially suffixed list of variables using the single hyphen shortcut; this would help in generalizing to n rows. Finally, note that we were able to avoid a second do loop (SAS' primary iterative looping syntax) mainly by luck-- the most recently generated value of a variable is saved by default. This can cause trouble, in general, but here it keeps all the previous "*"s when moving on to the next row.

An unrelated note about aggregatorsWe love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers and PROC-X with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.28: creating datasets from tables

2012-04-23T10:14:00.000-04:00

R
There are often times when it is useful to create an individual level dataset from aggregated data (such as a table). While this can be done using the expand.table() function within the epitools package, it is also straightforward to do directly within R.

Imagine that instead of the individual level data, we had only the 2x2 table for the association between homeless status and gender within the HELP RCT:


> HELPrct = read.csv("http://www.math.smith.edu/r/data/help.csv")
> xtabs(~ homeless + female, data=HELPrct)
        female
homeless   0   1
       0 177  67
       1 169  40

We can use this information to create an analytic dataset using just the four rows of a new dataset:


> female = c(0, 1, 0, 1)
> homeless = c(1, 1, 0, 0)
> count = c(169, 40, 177, 67)
> ds=data.frame(cbind(female, homeless, count))
> ds
  female homeless count
1      0        1   169
2      1        1    40
3      0        0   177
4      1        0    67

Next we use the rep() function to generate a vector of indices to repeat. The index object repeats each row number count times.


> index = rep(seq_len(nrow(ds)), times=ds$count)
> newds = ds[index,]
> newds$count = NULL
> xtabs(~ homeless + female, data=newds)
        female
homeless   0   1
       0 177  67
       1 169  40

The resulting data set is identical to the summarized input data set.

SAS
Many SAS procedures offer a weight varname option (as a statement within the proc) which will duplicate each observation varname times. So, for example, we can make a data set such as that shown above, then use, e.g., proc freq to produce a table.


data ds;
female = 0; homeless = 1; count = 169; output;
female = 1; homeless = 1; count = 40; output;
female = 0; homeless = 0; count = 177; output;
female = 1; homeless = 0; count = 67; output;
run;

proc freq data = ds;
table homeless * female;
weight count;
run;
                homeless     female

                Frequency|
                Percent  |
                Row Pct  |
                Col Pct  |       0|       1|  Total
                ---------+--------+--------+
                       0 |    177 |     67 |    244
                         |  39.07 |  14.79 |  53.86
                         |  72.54 |  27.46 |
                         |  51.16 |  62.62 |
                ---------+--------+--------+
                       1 |    169 |     40 |    209
                         |  37.31 |   8.83 |  46.14
                         |  80.86 |  19.14 |
                         |  48.84 |  37.38 |
                ---------+--------+--------+
                Total         346      107      453
                            76.38    23.62   100.00

However, some procedures lack this option, and/or it may be difficult to arrange your data appropriately to take advantage of it. In such cases, it's useful to be able to expand the data manually, as we show for R above. We demonstrate this below, assuming the count variable can be constructed. The explicit output statement puts a line into the newds data set count times.


data newds;
set ds;
do i = 1 to count;
  output;
  end;
run;

proc freq data = newds;
table homeless * female;
run;
                homeless     female

                Frequency|
                Percent  |
                Row Pct  |
                Col Pct  |       0|       1|  Total
                ---------+--------+--------+
                       0 |    177 |     67 |    244
                         |  39.07 |  14.79 |  53.86
                         |  72.54 |  27.46 |
                         |  51.16 |  62.62 |
                ---------+--------+--------+
                       1 |    169 |     40 |    209
                         |  37.31 |   8.83 |  46.14
                         |  80.86 |  19.14 |
                         |  48.84 |  37.38 |
                ---------+--------+--------+
                Total         346      107      453
                            76.38    23.62   100.00

Example 9.27: Baseball and shrinkage

2012-04-16T10:15:00.001-04:00

To celebrate the beginning of the professional baseball season here in the US and Canada, we revisit a famous example of using baseball data to demonstrate statistical properties.

In 1977, Bradley Efron and Carl Morris published a paper about the James-Stein estimator-- the shrinkage estimator that has better mean squared error than the simple average. Their prime example was the batting averages of 18 player in the 1970 season: they considered trying to estimate the players' average over the remainder of the season, based on their first 45 at-bats. The paper is a pleasure to read, and can be downloaded here. The data are available here, on the pages of Statistician Phil Everson, of Swarthmore College.

Today we'll review plotting the data, and intend to look at some other shrinkage estimators in a later entry.

SAS
We begin by reading in the data for Everson's page. (Note the long address would need to be on one line, or you could could use a URL shortener like TinyURL.com. To read the data, we use the infile statement to indicate a tab-delimited file and to say that the data begin in row 2. The informat statement helps read in the variable-length last names.


filename bb url "http://www.swarthmore.edu/NatSci/peverso1/Sports%20Data/
    JamesSteinData/Efron-Morris%20Baseball/EfronMorrisBB.txt";

data bball;
infile bb delimiter='09'x MISSOVER DSD lrecl=32767 firstobs=2 ;
informat firstname $7. lastname $10.;
input FirstName $ LastName $ AtBats Hits BattingAverage RemainingAtBats
   RemainingAverage SeasonAtBats SeasonHits SeasonAverage;
run;

data bballjs;
set bball;
js = .212 * battingaverage + .788 * .265;

avg = battingaverage; time = 1; 
  if lastname not in("Scott","Williams", "Rodriguez", "Unser","Swaboda","Spencer") 
    then name = lastname; else name = ''; 
output;
avg = seasonaverage; name = ''; time = 2; output;
avg = js; time = 3; name = ''; output;
run;

In the second data step, we calculate the James-Stein estimator according to the values reported in the paper. Then, to facilitate plotting, we convert the data to the "long" format, with three rows for each player, using the explicit output statement. The average in the first 45 at-bats, the average in the remainder of the season, and the James-Stein estimator are recorded in the same variable in each of the three rows, respectively. To distinguish between the rows, we assign a different value of time: this will be used to order the values on the graphic. We also record the last name of (most of) the players in a new variable, but only in one of the rows. This will be plotted in the graphic-- some players' names can't be shown without plotting over the data or other players' names.

Now we can generate the plot. Many features shown here have been demonstrated in several entries. We call out 1) the h option, which increases the text size in the titles and labels, 2) the offset option, which moves the data away from the edge of the plot frame, 3) the value option in the axis statement, which replaces the values of "time" with descriptive labels, and 4) the handy a*b=c syntax which replicates the plot for each player.


title h=3 "Efron and Morris example of James-Stein estimation";
title2 h=2 "Baseball players' 1970 performance estimated from first 45 at-bats";
axis1 offset = (4cm,1cm) minor=none label=none
  value = (h = 2 "Avg. of first 45" "Avg. of remainder" "J-S Estimator");
axis2 order = (.150 to .400 by .050) minor=none offset=(0.5cm,1.5cm) 
  label = (h =2 r=90 a = 270 "Batting Average");
symbol1 i = j v = none l = 1 c = black r = 20 w=3 
  pointlabel = (h=2 j=l position = middle "#name");

proc gplot data = bballjs;
  plot avg * time = lastname / haxis = axis1 vaxis = axis2 nolegend;
run; quit;

To read the plot (shown at the top) consider approaching the nominal true probability of a hit, as represented by the average over the remainder of the season, in the center. If you begin on the left, you see the difference associated with using the simple average of the first 45 at-bats as the estimator. Coming from the right, you see the difference associated withe using the James-Stein shrinkage estimator. The improvement associated with the James-Stein estimator is reflected in the generally shallower slopes when coming from the left. With the exception of Pirates great Roberto Clemente and declining third-baseman Max Alvis, most every line has a shallower slope from the left; James' and Stein's theoretical work proved that overall the lines must be shallower from the right.

R
A similar process is undertaken within R. Once the data are loaded, and a subset of the names are blanked out (to improve the readability of the figure), the matplot() and matlines() functions are used to create the lines.


bball = read.table("http://www.swarthmore.edu/NatSci/peverso1/Sports%20Data/JamesSteinData/Efron-Morris%20Baseball/EfronMorrisBB.txt",
                   header=TRUE, stringsAsFactors=FALSE)
bball$js = bball$BattingAverage * .212 + .788 * (0.265)
bball$LastName[!is.na(match(bball$LastName, 
  c("Scott","Williams", "Rodriguez", "Unser","Swaboda","Spencer")))] = ""

a = matrix(rep(1:3, nrow(bball)), 3, nrow(bball))
b = matrix(c(bball$BattingAverage, bball$SeasonAverage, bball$js), 
   3, nrow(bball), byrow=TRUE)
matplot(a, b, pch=" ", ylab="predicted average", xaxt="n", xlim=c(0.5, 3.1), ylim=c(0.13, 0.42))
matlines(a, b)
text(rep(0.7, nrow(bball)), bball$BattingAverage, bball$LastName, cex=0.6)
text(1, 0.14, "First 45\nat bats", cex=0.5)
text(2, 0.14, "Average\nof remainder", cex=0.5)
text(3, 0.14, "J-S\nestimator", cex=0.5)

Example 9.26: More circular plotting

2012-04-09T10:00:00.011-04:00

SAS's Rick Wicklin showed a simple loess smoother for the temperature data we showed here. Then he came back with a better approach that does away with edge effects. Rick's smoothing was calculated and plotted on a cartesian plane. In this entry we'll explore another option or two for smoothing, and plot the results on the same circular plot.

Since Rick is showing SAS code, and Robert Allison has done the circular plot (plot) (code), we'll stick to the R again for this one.

R
We'll start out by getting the data and setting it up as we did earlier. We add the year back into the matrix t3old because it'll be needed later.


temp1 = read.table("http://academic.udayton.edu/kissock/http/Weather/
           gsod95-current/NYALBANY.txt")
leap = c(0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1)
days = rep(365,18) + leap
monthdays = c(31,28,31,30,31,30,31,31,30,31,30,31)
temp1$V3 = temp1$V3 - 1994

yearpart = function(daytvec,yeardays,mdays=monthdays){
  part = (sum(mdays[1:(daytvec[1]-1)],(daytvec[1] > 2) * (yeardays[daytvec[3]]==366)) 
          + daytvec[2] - ((daytvec[1] == 1)*31)) / yeardays[daytvec[3]]
  return(part)
}

temp2 = as.matrix(temp1)

radians = 2* pi * apply(temp2, 1, yearpart, days, monthdays)

t3old = matrix(c(temp1$V4[temp1$V4 != -99 & ((temp1$V3 < 18) )],
                 radians[temp1$V4 != -99  &  ((temp1$V3 < 18) )],
                 temp1$V3[temp1$V4 != -99 & ((temp1$V3 < 18) )]), ncol=3)


t3now= matrix(c(temp1$V4[temp1$V4 != -99 & ((temp1$V3 == 18) | 
                           (temp1$V3 == 17 & temp1$V1 == 12))],
                radians[temp1$V4 != -99 & ((temp1$V3 == 18) | 
                           (temp1$V3 == 17 & temp1$V1 == 12))]), ncol=2)

library(plotrix)
radial.plot(t3old[,1],t3old[,2],rp.type="s", point.col = 2, point.symbols=46, 
    clockwise=TRUE, start = pi/2, label.pos = (1:12)/6 * (pi),
    radial.lim=c(-20,10,40,70,100),  labels=c("February 1","March 1",
      "April 1","May 1","June 1","July 1","August 1","September 1",
      "October 1","November 1","December 1","January 1"))

radial.plot(t3now[,1],t3now[,2],rp.type="s", point.col = 1, point.symbols='*', 
    clockwise=TRUE, start = pi/2, add=TRUE, radial.lim=c(-20,10,40,70,100))

If you didn't happen to see the update on the previous entry, note that the radial.lim option makes the axes for the added points match those for the initial plot. Otherwise, the added points plotted lower than they appeared, making the recent winter look cooler.

Rick started with a smoother, but often cyclic data can be fit well parametrically, using the sine and cosine of the cycle length as the covariates. With the data set up in radians already, this is trivially simple. The predicted values for the data can be retrieved with the fitted() function (e.g., section 3.7.3), which works with many model objects. These can then be fed into the radial.plot() function with rp.type="p" to make a line plot. The result is shown at the top-- the parametric fit appears to do a good job. Of course, you can fit on a square plot very easily with the plot() function, with result shown below.


simple = lm(t3old[,1] ~ sin(t3old[,2]) + cos(t3old[,2]))

radial.plot(fitted(simple),t3old[,2],rp.type="p", clockwise=TRUE,
            start = pi/2, add=TRUE, radial.lim=c(-20,10,40,70,100))

plot(t3old[,1] ~ t3old[,2], pch='.')
lines(t3old[,2],fitted(simple))

I didn't change the order of the data, so the line comes back to the beginning of the plot at the end of the year.

Adding a smoothed fit is nearly as easy. Just replace the lm() call with a loess() (section 5.2.6) call. The new line is added on top of the old one, to see just how they differ. The result is show below.


simploess = loess(t3old[,1] ~ t3old[,2])
radial.plot(fitted(simploess),t3old[,2],rp.type="p", line.col="blue", 
          clockwise=TRUE, start = pi/2, add=TRUE, radial.lim=c(-20,10,40,70,100))

The parametric fit is pretty good, but misses the sharp dip seen in January, and the fit in the late fall and early spring appear to be slightly affected.

But this approach stacks up all the data from 18 years. It might be more appropriate to stretch the data across the calendar time, fit the smoothed line to that, and then wrap it around the circular plot. To do this, we'll need to add the year back into the radians. Finding an acceptable smoother was a challenge-- the smooth.spline() function used here was adequate, but as the second plot below shows, it misses some highs and lows. Adding the smoothed curve to the plot is as easy as before, however. The plot with smoothing by year is immediately below.


radyear = t3old[,2] + (2 * pi * t3old[,3])
better = smooth.spline(y=t3old[,1],x= radyear, all.knots=TRUE,spar=1.1)

radial.plot(fitted(better),t3old[,2],rp.type="p", line.col="green", 
        clockwise=TRUE, start = pi/2, add=TRUE, radial.lim=c(-20,10,40,70,100))

plot(t3old[,1] ~ radyear, pch = '.')
lines(better)

The relatively poor fit seen below makes the new (green) line at least as poor as the parametric fit. The extra variability in the winter is reflected in distinct lines in the winter. Rick's approach, to fit the data lumping across years, seems to be the best for fitting, though it's easier to see the heteroscedaticity in the ciruclar plot. But however you slice it, this winter has had an unusual number of very warm days.

An unrelated note about aggregators
We love aggregators! Aggregators are meta-blogs that collect content from blogs that have similar coverage, for the convenience of readers. For blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers and PROC-X with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.25: It's been a mighty warm winter? (Plot on a circular axis)

2012-04-02T09:33:00.009-04:00

Updated (see below)

People here in the northeast US consider this to have been an unusually warm winter. Was it?

The University of Dayton and the US Environmental Protection Agency maintain an archive of daily average temperatures that's reasonably current. In the case of Albany, NY (the most similar of their records to our homes in the Massachusetts' Pioneer Valley), the data set as of this writing includes daily records from 1995 through March 12, 2012.

In this entry, we show how to use R to plot these temperatures on a circular axis, that is, where January first follows December 31st. We'll color the current winter differently to see how it compares. We're not aware of a tool to enable this in SAS. It would most likely require a bit of algebra and manual plotting to make it work.

R
The work of plotting is done by the radian.plot() function in the plotrix package. But there are a number of data management tasks to be employed first. Most notably, we need to calculate the relative portion of the year that's elapsed through each day. This is trickier than it might be, because of leap years. We'll read the data directly via URL, which we demonstrate in Example 8.31. That way, when the unseasonably warm weather of last week is posted, we can update the plot with trivial ease.


library(plotrix)
temp1 = read.table("http://academic.udayton.edu/kissock/http/
            Weather/gsod95-current/NYALBANY.txt")
leap = c(0,1,0,0,0,1,0,0,0,1,0,0,0,1,0,0,0,1)
days = rep(365, 18) + leap
monthdays = c(31,28,31,30,31,30,31,31,30,31,30,31)
temp1$V3 = temp1$V3 - 1994

The leap, days, and monthdays vectors identify leap years, count the corrrect number of days in each year, and have the number of days in the month in non-leap years, respectively. We need each of these to get the elapsed time in the year for each day. The columns in the data set are the month, day, year, and average temperature (in Fahrenheit). The years are renumbered, since we'll use them as indexes later.

The yearpart() function, below, counts the proportion of days elapsed.


yearpart = function(daytvec,yeardays,mdays=monthdays){
  part = (sum(mdays[1:(daytvec[1]-1)],
           (daytvec[1] > 2) * (yeardays[daytvec[3]]==366)) 
          + daytvec[2] - ((daytvec[1] == 1)*31)) / yeardays[daytvec[3]]
  return(part)
}

The daytvec argument to the function will be a row from the data set. The function works by first summing the days in the months that have passed (,sum(mdays[1:(daytvec[1]-1)]) adding one if it's February and a leap year ((daytvec[1] > 2) * (yeardays[daytvec[3]]==366))). Then the days passed so far in the current month are added. Finally, we subtract the length of January, if it's January. This is needed, because sum(1:0) = 1, the result of which is that that January is counted as a month that has "passed" when the sum() function quoted above is calculated for January days. Finally, we just divide by the number of days in the current year.

The rest is fairy simple. We calculate the radians as the portion of the year passed * 2 * pi, using the apply() function to repeat across the rows of the data set. Then we make matrices with time before and time since this winter started, admittedly with some ugly logical expressions (section 1.14.11), and use the radian.plot() function to make the plots. The options to the function are fairly self-explanatory.


temp2 = as.matrix(temp1)
radians = 2* pi * apply(temp2,1,yearpart,days,monthdays)

t3old = matrix(c(temp1$V4[temp1$V4 != -99 & ((temp1$V3 < 18) | (temp1$V2 < 12))],
          radians[temp1$V4 != -99  &  ((temp1$V3 < 18) | (temp1$V2 < 2))]),ncol=2)

t3now= matrix(c(temp1$V4[temp1$V4 != -99 & 
          ((temp1$V3 == 18) | (temp1$V3 == 17 & temp1$V1 == 12))],
          radians[temp1$V4 != -99 & ((temp1$V3 == 18) | 
          (temp1$V3 == 17 & temp1$V1 == 12))]),ncol=2)
# from plottrix library
radial.plot(t3old[,1],t3old[,2],rp.type="s", point.col = 2, point.symbols=46,
            clockwise=TRUE, start = pi/2, label.pos = (1:12)/6 * (pi), 
            labels=c("February 1","March 1","April 1","May 1","June 1",
            "July 1","August 1","September 1","October 1","November 1",
            "December 1","January 1"), radial.lim=c(-20,10,40,70,100))

radial.plot(t3now[,1],t3now[,2],rp.type="s", point.col = 1, point.symbols='*', 
            clockwise=TRUE, start = pi/2, add=TRUE, radial.lim=c(-20,10,40,70,100))

The result is shown at the top. The dots (point.symbol is like pch so 20 is a point (section 5.2.2) show the older data, while the asterisks are the current winter. An alternate plot can be created with the rp.type="p" option, which makes a line plot. The result is shown below, but the lines connecting the dots get most of the ink and are not what we care about today.

Either plot demonstrates clearly that a typical average temperature in Albany is about 60 to 80 in August and about 10 to 35 in January, the coldest monthttp://www.blogger.com/img/blank.gifh.

Update
The top figure shows that it has in fact been quite a warm winter-- most of the black asterisks are near the outside of the range of red dots. Updating with more recent weeks will likely increase this impression. In the first edition of this post, the radial.lim option was omitted, which resulted in different axes in the original and "add" calls to radial.plot. This made the winter look much cooler. Many thanks to Robert Allison for noticing the problem in the main plot. Robert has made many hundreds of beautiful graphics in SAS, which can be found here. He also has a book. Robert also created a version of the plot above in SAS, which you can find here, with code here. Both SAS and R (not to mention a host of other environments) are sufficiently general and flexible that you can do whatever you want to do-- but varying amounts of expertise might be required.

An unrelated note about aggregators
We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers and PROC-X with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.