SAS and R

SAS Macro Simplifies SAS and R integration (Updated)

2012-01-26T16:59:00.003-05:00

Many of us feel very enthusiastic about R. It's free, it features cutting edge applications, it has a large community of users contributing for mutual benefit, and on and on. There are also many things to like about SAS, including stability, backwards compatibility, and professional support among them. The way to be the best analyst you can be is to be flexible and have as many tools at your disposal as you can manage. That's the main motivating principle behind our book and what we do here in this blog.

Today we call attention to a SAS macro that greatly eases integrating R from SAS. Published last month in the Journal of Statistical Software, the macro (written by Xin Wei of Roche Pharmaceuticals) is called Proc_R, and we discuss its installation and use today. For a fuller write-up, see the paper, here. For SAS users, the macro is a huge productivity booster, allowing one to easily complete data management and/or partial data analysis in SAS, skip out quickly to R for analyses that are awkward or impossible in SAS, then return to SAS for completion. For people in industry, this may also ease integrating R into documentation systems built for SAS code. See this post on DecisionStats for a review of other integration attempts.

Getting ready

1. Download the "SAS source code" and the "Replication code and instructions".

2. Move the macro somewhere you have write access.

3. Open the macro in a text editor and change line 46 so that the rpath option points to the location of your R executable.

(4. If you're running Windows 7 or Vista, and you has SAS 9.1 or above, follow instructions in a PDF in the second supplemental file you downloaded. This makes a shortcut for a special version of SAS. I'm not at all sure why you have to do this, though. I had the same results running in my usual SAS set-up.)

That's it! The way the macro works is to read in your R code as a SAS data set, write it out to a file, and call R to run it, then do a bunch of post-processing. The basic macro call looks like this:


%include "C:\ken\sasmacros\Proc_R.sas";
%Proc_R (SAS2R =, R2SAS =);
Cards4;

******************************
***Please Enter R Code Here***
******************************

;;;;
%Quit;

You just replace the starred lines with R code, and run-- the R results, if any, appear in your SAS output and/or results windows. The SAS2R value is a list of the names of SAS data sets you'd like to send to R; they're added into the R environment before your code is executed. The R2SAS value is a list of the names of R objects (that can be coerced to data frames) that you'd like to become SAS data sets.

Use
Here's a trivial example-- generate two data sets in SAS, send them to R to run linear regressions, and send the resulting parameter estimates back to SAS.


data test;
do i = 1 to 1000;
  x = normal(0);
  y = x + normal(0);
  output;
  end;
run;

data t2;
do i = 1 to 100;
  x = normal(0);
  y = x + uniform(0);
  output;
  end;
run;

%include "C:\Proc_R.sas";
%Proc_R (SAS2R =test t2, R2SAS =mylm mylm2);
Cards4;
setwd("c:/temp")
an.lm = with(test,lm(y ~x))
mylm = t(coef(an.lm))

an.lm2 = with(t2,lm(y~x))
mylm2 = t(coef(an.lm2))
;;;;
%Quit;

proc print data = mylm; run;
proc print data = mylm2; run;

And here's what you get in the SAS log.


[First, proc_r result]

******************R OUTPUT***********************    

R_OUTPUT_LOG

> setwd("c:/temp")
> library(grDevices)
> png("c:/temp/....png")
> test<- read.csv('c:/temp/test.csv')
> t2<- read.csv('c/temp/t2.csv')
> an.lm = with(test,lm(y ~x))
> mylm = t(coef(an.lm))
> summary(an.lm)

Call:
lm(formula = y ~ x)

Residuals:
Min      1Q  Median      3Q     Max
-2.8571 -0.6430 -0.0051  0.6713  3.5903

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.008568   0.031686    0.27    0.787
x           1.020640   0.033315   30.64   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.002 on 998 degrees of freedom
Multiple R-squared: 0.4846, Adjusted R-squared: 0.4841
F-statistic: 938.5 on 1 and 998 DF,  p-value: < 2.2e-16

>
> an.lm2 = with(t2,lm(y~x))
> mylm2 = t(coef(an.lm2))
> write.csv(mylm,'mylm.csv',row.names=F)
> write.csv(mylm2,'mylm2.csv',row.names=F)
> dev.off()
null device
1
> q()
> proc.time()
user  system elapsed
0.28    0.10    0.37


[Here are the proc print results]
Obs     _Intercept_               x
 1     0.0085676126    1.0206400545

Obs     _Intercept_               x
 1      0.528410053    0.9851225238

(Page breaks and some extraneous stuff removed.)

It's pretty magical for a SAS user to see R living in the SAS output like this. But there are some caveats. First, this is a windows-only macro. If you run SAS on *nix, you may not be able to get it to work. Second, while the article has examples of graphics from R neatly appearing in SAS, this failed for me. This may be due to the fact that I run SAS 9.3, while the author of the macro is still in earlier versions of SAS. I may try to diagnose and fix this problem, and will update this entry if I find a fix. (Fixed! See update below.)

(UPDATE: Reader Abhijit suggested a setwd() in the R code as a fix for the graphics problem. This works, and I now get R grapics in my SAS results viewer. Even more magical. Code and output above updated to show this. Thanks, Abhijit!)

However, these seem like minor problems, compared with the overall simplification offered by the macro. It's been of great use to me in the past few months, and I expect it will help others as well. Many thanks and congratulations to Xin Wei!

Example 9.19: Demonstrating the central limit theorem

2012-01-11T15:30:00.001-05:00

A colleague recently asked "why should the average get closer to the mean when we increase the sample size?" We should interpret this question as asking why the standard error of the mean gets smaller as n increases. The central limit theorem shows that (under certain conditions, of course) the standard error must do this, and that the mean approaches a normal distribution. But the question was why does it? This seems so natural that it may have gone unquestioned in the past.

The best simple rationale may be that there are more ways to get middle values than extreme values--for example, the mean of a die roll (uniform discrete distribution on 1, 2, ..., 6) is 3.5. With one die, you're equally likely to get an "average" of 3 or of 1. But with two dice there are five ways to get an average of 3, and only one way to get an average of 1. You're 5 times more likely to get the value that's closer to the mean than the one that's further away.

Here's a simple graphic to show that the standard error decreases with increasing n.

SAS
We begin by simulating some data-- normal, here, but of course that doesn't matter (assuming that the standard deviation exists for whatever distribution we pick and the sample size is appropriately large). Rather than simulate separate samples with n = 1 ... k, it's easier to add a random variate to a series and keep a running tally of the mean, which is easy with a little algebra. This approach also allows tracking the progress of the mean of each series, which could also be useful.


%let nsamp = 100;
data normal;
do sample = 1 to &nsamp;
  meanx = 0;
  do obs = 1 to &nsamp;
    x = normal(0);
 meanx = ((meanx * (obs -1)) + x)/obs;
 output;
  end;
end;
run;

We can now plot the means vs. the number of observations.


symbol1 i = none v = dot h = .2;
proc gplot data = normal;
plot meanx * obs;
run;
quit;

symbol1 i=join v=none r=&nsamp;
proc gplot data=normal;
  plot meanx * obs = sample / nolegend;
run; quit;

The graphic resulting from the first proc gplot is shown above, and demonstrates both how quickly the variability of the estimate of the mean decreases when n is small, and how little it changes when n is larger. A plot showing the means for each sequence converging can be generated with the second block of code. Note the use of the global macro variable nsamp assigned using the %let statement (section A.8.2).

R
We'll also generate sequences of variates in R. We'll do this by putting the random variates in a matrix, and treating each row as a sequence. We'll use the apply() function (sections 1.10.6 and B.5.3) to treat each row of the matrix separately.


numsim = 100
matx = matrix(rnorm(numsim^2), nrow=numsim)

runavg = function(x) { cumsum(x)/(1:length(x)) }
ramatx = t(apply(matx, 1, runavg))

The simple function runavg() calculates the running average of a vector and returns the a vector of equal length. By using it as the function in apply() we can get the running average of each row. The result must be transposed (with the t() function, section 1.9.2) to keep the sequences in rows. To plot the values, we'll use the type="n" option to plot(), specifying the first column of the running total as the y variable. While it's possible that the running average will surpass the average when n=1, we ignore that case in this simple demonstration.


plot(x=1:numsim, y = ramatx[,1], type="n",
  xlab="number of observations", ylab="running mean")
rapoints = function(x) points(x~seq(1:length(x)), pch=20, cex=0.2)
apply(ramatx,1,rapoints)

plot(x=1:numsim, y = ramatx[,1], type="n",
  xlab="number of observations", ylab="running mean")
ralines = function(x) lines(x~seq(1:length(x)))
apply(ramatx, 1, ralines)

Here we define another simple function to plot the values in a vector against the place number, then again use the apply() function to plot each row as a vector. The first set of code generates a plot resembling the SAS graphic presented above. The second set of code will connect the values in each sequence, with results shown below.

Example 9.18: Constructing the fastest relay team via enumeration

2012-01-05T14:51:00.009-05:00

In competitive swimming, the medley relay is a team event in which four different swimmers each swim one of the four strokes: freestyle, breaststroke, backstroke, and butterfly. In general, every swimmer might be able swim any given stroke. How can we choose the fastest relay team? Here we solve this by enumerating all possible teams, though a more efficient routine likely exists.

Some example practice times can be seen on this Google Spreadsheet. Also, using the steps outlined here, the same spreadsheet is available as a CSV file here. (FTR, these are actual practice times for 100 yards for mostly 12-13 year-old swimmers; the names have been changed.)

SAS
We first read the data from the URL, using the technique outlined in section 1.1.6. Note that if you cut-and-paste this, you'll need to get the whole URL onto one line-- we break it up here for display only.


filename swimurl url 'https://docs.google.com/spreadsheet
/pub?key=0AvJKgZUzMYLYdE5xTHlEWkNUM3NoOHB1ZTJoTFMzUUE&
single=true&gid=0&output=csv';

proc import datafile=swimurl out=swim dbms=csv;
run;

Next, we use the point= option in nested set statements to generate a single data set with all the possible combinations of names and times. Meanwhile we change the names of the variables so they don't get overwritten in the next set statement. Note the use of the nobs option to find the number of rows in the data set.


data permute 
    (keep=free freetime fly flytime back backtime  breast breasttime);
set swim (rename = (swimmer=free freestyle=freetime)) nobs=nobs;
do i = 1 to nobs;
  set swim(rename = (swimmer=fly butterfly=flytime)) point=i;
  do j = 1 to nobs;
    set swim(rename = (swimmer=back backstroke=backtime)) point=j;
    do k = 1 to nobs;
      set swim(rename = (swimmer=breast breaststroke=breasttime)) 
        point = k;
      output;
 end;
  end;
end;
run;

The resulting data set has 12^4 rows, and includes rosters with the same swimmer swimming all four legs. In fact, a quick glance will show that Anna has the best time in each stroke, and thus the best "team" based on these practice times would use her for each stroke. This is against the rules, and also probably isn't reflective of performance in a race. We'll remove illegal line-ups using a where statement (section 1.5.1) and also calculate the total team time.


data prep;
set permute;
where free ne back and free ne breast and free ne fly and 
  back ne breast and back ne fly and breast ne fly;
time = sum(freetime, flytime, backtime, breasttime);
run;

The resulting data set has (12 permute 4) lines. To find the best team, we just sort by the total time and look at the first line. Here the first 10 lines (10 best teams) are shown.


proc sort data=prep; by time; run;
proc print data=prep (obs=10); run;

                                                        b
                                                        r
           f                             b              e
           r              f              a              a
           e              l              c      b       s
           e              y              k      r       t
   f       t              t      b       t      e       t      t
   r       i      f       i      a       i      a       i      i
   e       m      l       m      c       m      s       m      m
   e       e      y       e      k       e      t       e      e

 Kara    109.3  Dora    126.8  Lara    117.7  Anna    126.9  480.7
 Anna    102.8  Dora    126.8  Lara    117.7  Beth    134.6  481.9
 Kara    109.3  Anna    120.5  Lara    117.7  Beth    134.6  482.1
 Anna    102.8  Dora    126.8  Lara    117.7  Honora  136.4  483.7
 Kara    109.3  Jane    129.8  Lara    117.7  Anna    126.9  483.7
 Kara    109.3  Anna    120.5  Lara    117.7  Dora    136.4  483.9
 Kara    109.3  Anna    120.5  Lara    117.7  Honora  136.4  483.9
 Kara    109.3  Lara    123.1  Jane    124.7  Anna    126.9  484.0
 Anna    102.8  Dora    126.8  Lara    117.7  Inez    136.8  484.1
 Carrie  112.7  Dora    126.8  Lara    117.7  Anna    126.9  484.1

The best time shaves a whole second off the predicted time using the second-best team.

R
Since published Google Spreadsheets use https rather than http, we use the RCurl package and its getURL() function. (Note that if you cut-and-paste this, you'll need to get the whole URL onto one line-- we break it up here for display only.) Then we can read the data with read.csv() and textConnection().


library(RCurl)
swim = getURL("https://docs.google.com/spreadsheet
/pub?key=0AvJKgZUzMYLYdE5xTHlEWkNUM3NoOHB1ZTJoTFMzUUE&
single=true&gid=0&output=csv")

swim2=read.csv(textConnection(swim))

To make an object with the combinations of names, we use the expand.grid() function highlighted in Example 7.22, providing as arguments the swimmers names four times. As in the SAS example, the result has has 12^4 rows. The combn() function might be a better fit here, but was more difficult to use.


test2 = expand.grid(swim2$Swimmer,swim2$Swimmer, swim2$Swimmer, swim2$Swimmer)

It'll be useful to assign these copies of the names to each of the strokes. We'll do that with the rename() function available in the reshape package. (This approach is mentioned in section 1.3.4.). Then we can remove the rows where the same name appears twice using some logic. The logic is nested in the with() function to save some keystrokes and is generally preferable to attach()ing the test2 object.


library(reshape)
test2 = rename(test2, c("Var1" = "free", "Var2" = "fly", 
  "Var3" = "back", "Var4" = "breast"))
test3 = with(test2, test2[(free != breast) & (free != fly) 
  & (free != back) & (breast != fly) & (breast != back) 
  & (fly != back) ,])

Finally, we can use the which.min() function to pick the best team.


> bestteam =   
+ test3[which.min(swim2$Freestyle[test3$free]+swim2$Breaststroke[test3$breast] +
+ swim2$Butterfly[test3$fly] + swim2$Backstroke[test3$back]),]
>bestteam
     free  fly back breast
1631 Kara Dora Lara   Anna

For new users of R, this may look very peculiar-- it uses elegant but powerful features of the R language that may be challenging for new users to grasp. Essentially, in swim2$Freestyle[test3$free] we say: from the "freestyle" times in the swim2 object, take the time from the row that has the name in a row of "free" names in the test3 object. The which.min() function replicates this request for every row in the test3 object (which has all of the permutations) in it, returning the row number with that minimum sum. The outer test3[rows,columns] syntax grabs the values in this row. (The number 1631 is the row number, for some reason showing the row in the test2 object created by expand.grid().)

Now, we might also want the actual times associated with the best team. We can find them by calling the correct rows (names from the best team) and columns (stroke associated with that name) from the original data set.


> times = c(swim2[swim2$Swimmer == bestteam$free,2], 
+      swim2[swim2$Swimmer == bestteam$fly,3], 
+      swim2[swim2$Swimmer == bestteam$back,4], 
+      swim2[swim2$Swimmer == bestteam$breast,5])
> times
[1] 109.3 126.8 117.7 126.9

If instead, one wanted to list the times in order, one approach would be to add columns to the test3 object with the time for each stroke, calculate their sum, and sort on the sum.

Example 9.17: (much) better pairs plots

2011-12-06T09:24:00.004-05:00

Pairs plots (section 5.1.17) are a useful way of displaying the pairwise relations between variables in a dataset. But the default display is unsatisfactory when the variables aren't all continuous. In this entry, we discuss ways to improve these displays that have been proposed by John Emerson, Walton Green, Barret Schloerke, Dianne Cook, Heike Hofmann, and Hadley Wickham in a manuscript under review entitled The Generalized Pairs Plot. http://www.blogger.com/img/blank.gif

Implementations of the methods in the paper are available in the gpairs and GGally packages; here we use the latter, which is based on the grammar of graphics and the ggplot2 package. This is an R-only entry: we are unaware of efforts to replicate this approach in SAS.

New users may find it easier to break process down into steps, rather than to do everything at once, as the R language allows. One way to do that is to make a smaller version of a dataset, with just the analysis variables included. here we use the HELP data set and choose two categorical variables (gender and housing status) and two continuous ones (the number of drinks per day and a measure of depressive symptoms). Once this new subset is created, the call to ggpairs() is straightforward.

R


library(GGally)
ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
ds$sex = as.factor(ifelse(ds$female==1, "female", "male"))
ds$housing = as.factor(ifelse(ds$homeless==1, "homeless", "housed"))
smallds = subset(ds, select=c("housing", "sex", "i1", "cesd"))
ggpairs(smallds, diag=list(continuous="density", discrete="bar"), axisLabels="show")

For users more comfortable with R, the ggpairs function allows you to select variables to include, via its columns option. The following line produces a plot identical to the above, without the subset().


ggpairs(ds, columns=c("housing", "sex", "i1", "cesd"),
    diag=list(continuous="density",   discrete="bar"), axisLabels="show")

Various options are available for the diagonal elements of the plot matrix, and the off-diagonals can be controlled with upper and lower options. The examples(ggpairs) command is very helpful for visualizing some of the possibilities.

Example 9.16: Small multiples

2011-11-29T09:09:00.010-05:00

Small multiples are one of the great ideas of graphics visionary Edward Tufte (e.g., in Envisioning Information). Briefly, the idea is that if many variations on a theme are presented, differences quickly become apparent. Today we offer general guidance on creating figures with small multiples.

As an example, we'll show graphics for the popularity, salary, and unemployment rates for college majors. This data was discussed here where a scatterplot graphic was presented. We draw on data and code presented there as well. The scatterplot does not show the unemployment rate, and the width of the field names and arbitrarily sized text makes it difficult to determine the popularity ranking. In contrast, the small multiples plot, demonstrated above, makes each of these features easy to see. (Click on the picture for a clearer image.)

R

The graphics options in R, particularly par("mfrow") or par("mfcol"), are well-suited to small multiples. The main tip here is to minimize the space reserved for margins and titles. In the example below, we do this with the mar, mgp, and oma settings. We'll begin by setting up the data in a process that relies heavily on the code here. (Note that the zip file referred to includes data already in the R format-- since our point today is not data management, we don't replicate the process used to make this out of the raw data.)


clean = function(x) as.numeric(gsub("\\$|, ", "", x))
clean2 = function(x) as.numeric(gsub("%", "", x))
load("table.rda")
X[,2] = clean2(X[,2])
for (i in 3:5) X[,i] = clean(X[,i])


X$cols = NA
X$cols[grep("BUSI|ACC|FINAN",X[,1])] = 1
X$cols[grep("ENGINEERING",X[,1])] = 2
X$cols[grep("STAT|COMPU",X[,1])] = 3
X$cols[grep("BIOL|CHEMI|PHYSICS|MATHEM",X[,1])] =  4
X$cols[grep("ENGLISH|HISTORY|LANG|FINE|MUSIC|PHILOS|DRAMA|LIBERAL|ARCH|THEO",X[,1])] = 5
X$cols[grep("SOCIO|PSYCH|POLI|ECONO|JOURNAL",X[,1])] = 6
X$cols[grep("COMMUN|MARKET|COMMER|MULTI|MASS|ADVERT",X[,1])] = 7
X$cols[grep("NURS|CRIM|EDU|PHYSICAL|FAMI|SOCIAL|TREAT|HOSP|HUMAN",X[,1])] = 8

This removes some funny characters and groups the fields together in a coherent manner. Then we write a function to set the par() values we need to change, and plot a pie for each row of the data set. Here a for loop is used-- we're not aware of a vectorized version of the pie() function. Colors for each pie are assigned via the colors.


sortx = X[order(X$Popularity),]

smajors = function(mt) {
  old.par = par(no.readonly=TRUE) 
  on.exit(par(old.par))
  nrows = sqrt(nrow(mt)) + (ceiling(sqrt(nrow(mt))) !=  sqrt(nrow(mt)))
  par(mfrow=c(nrows,nrows), mar=c(1,0,0,0), oma=c(0,0,0,0), mgp=c(0,0,0))
  colors = c("blue", "purple", "purple", "gray", "orange", "gray", "red", "black")
  for (i in 1:nrow(mt)) {
    pie(c(mt[i,2], 100 - mt[i,2]), labels=NA, radius=mt[i,4]/max(mt[,4]), 
        col = c("white",colors[mt[i,7]]))
    mtext(substr(mt[i,1],1,19), side=1, cex=.8)
  }
}

smajors(sortx[1:49,])

The resulting plot is shown above. There may be too many colors, though statistics and computing were assigned the same color as engineering. We can easily read from the plot that computing, statistics, and engineering (purple) are the largest circles, and thus the best paying. Similarly, the humanities (orange) are the worst paying. They are also not terribly popular-- the first orange pie appears in the second row. The "professions" (nursing, teaching, policing, therapy) don't pay well but are fairly popular. Most pies have roughly equal unemployment, though nursing and the professions generally are notable for near full employment. All in all, the rank, salary, unemployment percentage, and field are all clearer than in the scatterplot.

SAS

In SAS, one can use the greplay procedure to reproduce images in miniature. One has to define the size and shape allotted for each replayed plot, in a stored "template." This allows enormous control, but at the cost of some complexity. For example, one can create a scatterplot matrix using proc gplot instead of proc sgplot, as in this implementation by Michael Friendly. If you can generate your multiple images with a by statement, the coding for this is not too painful. However, in this case, it was not obvious how to change the color and radius for each pie using a by statement in proc gchart which includes pie charts, and would thus have been the obvious choice. In such cases, it may be easier to plot the figures directly using an annotate data set.

However, having demonstrated the use of annotate previously (e.g., Example 8.13), we show here an application using greplay, though the coding is a little bit involved. In outline, we use a macro to call for a pie to be made from each observation of the original data set. Then we use a template-making macro found here to generate the 7X7 grid template. Finally, we replay the pies into the grid.


/* read the data-- note that text file edited to remove spaces in 
    variable names */
proc import datafile = 
   'c:\sas-r dictionary\p\book\web\blog\majors\majors\table.txt'
   out = majors;
   getnames = yes;
   run;

/* set up the field categories and colors */
data m2;
set majors;
cols="      ";     /* make a missing character variable to hold them */
mf = majorfield;   /* just a copy to save keystrokes */
   /* the find() function is discussed in section 1.4.6 */
if sum(find(mf,"BUSI"), find(mf,"ACC"), find(mf,"FINAN")) ne 0 then cols = "blue";
if find(mf,"ENGINEERING") ne 0 the cols = "purple";
if sum(find(mf,"STAT"), find(mf,"COMPU"), find(mf,"SYSTEMS")) ne 0 then cols = "purple";
if sum(find(mf,"BIOL"), find(mf,"CHEMI"), find(mf,"PHYSICS"), find(mf,"MATH")) ne 0 then cols = "gray";
if sum(find(mf,"ENGL"), find(mf,"HIST"), find(mf,"FRENCH"), find(mf,"FINE"),
     find(mf,"MUSIC"), find(mf,"PHIL"), find(mf,"DRAMA"), find(mf,"LIBERAL"),
     find(mf,"ARCH"), find(mf,"THEO")) ne 0 then cols = "orange";
if sum(find(mf,"SOCI"), find(mf,"PSYCH"), find(mf,"POLI"), find(mf,"ECON"),
     find(mf,"JOURN"), find(mf,"LIBERAL")) ne 0 then cols = "gray";
if sum(find(mf,"COMMU"), find(mf,"MARKET"), find(mf,"COMMER"), find(mf,"MULTI"),
     find(mf,"MASS"), find(mf,"ADVERT")) ne 0 then cols = "red";
if sum(find(mf,"NURS"), find(mf,"CRIM"), find(mf,"EDU"), find(mf,"PHYSICAL"),
     find(mf,"FAMI"), find(mf,"SOCIAL"), find(mf,"TREAT"),
     find(mf,"HOSP"), find(mf,"HUMAN")) ne 0 then cols = "black";
drop MF;   /* get rid of that keystroke-saver */
run; 

/* order by popularity */
proc sort data = m2; by popularity; run;

The next macro takes a line number as input. A data step then reads that line from the real data set and uses call symput (section A.8.2) to extract as macro variables the median earnings used for the radius, the color, and the major name. It then produces two rows of data-- one with the unemployed percent and the other with the employed percent. We need this for the proc gchart input, as shown in the second half of the macro.


%macro smpie(obs);
data t1;
set m2 (firstobs = &obs obs = &obs);
call symput('rm', medianearnings);
call symput('color', cols);
call symput('name', strip(substr(majorfield,1,19)));
employed = "No"; percent = unemploymentpercent; output;
employed = "Yes"; percent = 100 - unemploymentpercent; output;
run;

pattern1 color= white v=solid;
pattern2 color= &color v=solid;  /* only pattern2 should be needed, I think, but */
pattern3 color= &color v=solid;  /* sometimes SAS required pattern3, in my trials */
title h=7pct "&name";
proc gchart data = t1 gout = kkpies;
pie percent / group = majorfield sumvar = percent value = none 
    noheading nogroupheading radius = %sysevalf((&rm * 45)/ 105000)
    name = "PIE&obs" ; 
run;quit;
%mend smpie;

Of particular note in the forgoing are the gout and name options. The former specifies a location where output plots should be stored. The latter assigns a name to this particular plot. In addition, the %sysevalf macro function is needed to perform mathematical functions on the radius variable.

Next, we write and call another macro to call the first repeatedly. Making the image small to begin with (using the goptions statement [section 5.3]) reduces quality loss when replaying as small multiples.


%macro pies;
goptions hsize=1in vsize=1in;
%do i = 1 %to 49;
  %smpie(&i);
%end;
%mend;  

%pies;

Finally, we can make the template to replay the images, and replay them, both using proc greplay.


/* key elements of the %makegridtemplate macro: 
     how many rows and columns (down and across). 
     name for the template.
   Note that the macro is called *inside* the proc greplay statement, 
   which allows the user access to  pro greplay statment options */
proc greplay nofs  tc=work.templates;
  %makeGridTemplate (across=7, down=7, ordering=LRTB, 
      hgap=0, vgap=0, gapT=0, gapL=0, gapr=0, name=sevensq,bordercolor=white)
quit;

/* this macro just types out text for us the sequence 1:pie1 2:pie2 ... 49:pie49 */
/* we need that to replay the figures in proc greplay */
%macro pielist;
%do i = 1 %to 49; &i:pie&i %end;;
%mend;

filename pies "c:\temp\pies2.png";
goptions hsize=7in vsize=7in gsfname=pies device=png;

/* now ready to replay the plots
   The proc greplay options say what template to use and where to find it, 
     and where to find the input and place the output */
/* The treplay statement plays the old plots to the locations specified 
     in the template
proc greplay template=sevensq tc=work.templates nofs gout=kkpies igout=kkpies;
treplay %pielist;
run;
quit;

The net outcome of this is shown below. The image is pretty disappointing-- the circles are not round,and the text is pretty blurry. However, the message is as clear as in the prettier R version.

Example 9.15: Bar chart with error bars ("Dynamite plot")

2011-11-22T10:10:00.000-05:00

The "dynamite plot", a bar chart plotting the a mean with a error bar, is one of the most reviled types of image among statisticians. Reasons to dislike them are numerous, and are nicely summarized here. (Edward Tufte also suggests they be avoided.)

Nonetheless, as consulting statisticians, we're often required to meet the needs of our collaborators, or of the reviewers who will judge their submissions. If you need to do it, here's how. We demonstrate with the HELP data set.

SAS

The plot can be created easily with proc gchart (section 5.1.3).


proc gchart data = "c:\book\help.sas7bdat";
vbar substance /
   group=female
   type=mean 
   sumvar=sexrisk
   errorbar=top;
run;

The syntax requests the basic chart variable be substance, in groups defined by female, where the mean of sexrisk is plotted. The errorbar option can be used to make the full CI or just the top lines; the default is a 95% standard error of the mean, but that can also be adjusted. The result is shown below.

R

There are several easily googlable examples of how to do a basic dynamite plot in R, including one in example(barchart). However, a grouped example with a nice legend is a bit harder to find. The bargraph.CI() function in the sciplot package fits the bill.


library(sciplot)
bp = with(ds, bargraph.CI(x.factor=female, group=substance, response=sexrisk,
  lc=FALSE, xlab="Female",
  legend=TRUE, x.leg=3.3, cex.leg=1.3, cex.names=1.5, cex.lab = 1.5,
  ci.fun=function(x) {c(mean(x) - 1.96*se(x), mean(x) + 1.96*se(x))}))

The results are shown at the top. The first row of options describe the variables to be used; the default plot statistic is the mean. The second and third rows suppress the bottom of the confidence interval and customize the location, label, and font size in the legend and the x-axis label. The only tricky bit is the last row, which provides the 95% CI limits to be plotted, as opposed to the default 1 standard error.

Example 9.14: confidence intervals for logistic regression models

2011-11-15T09:53:00.001-05:00

Recently a student asked about the difference between confint() and confint.default() functions, both available in the MASS library to calculate confidence intervals from logistic regression models. The following example demonstrates that they yield different results.

R


ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
library(MASS)
glmmod = glm(homeless ~ age + female, binomial, data=ds)

> summary(glmmod)
Call:
glm(formula = homeless ~ age + female, family = binomial, data = ds)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3600  -1.1231  -0.9185   1.2020   1.5466  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.89262    0.45366  -1.968   0.0491 *
age          0.02386    0.01242   1.921   0.0548 .
female      -0.49198    0.22822  -2.156   0.0311 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 625.28  on 452  degrees of freedom
Residual deviance: 617.19  on 450  degrees of freedom
AIC: 623.19

Number of Fisher Scoring iterations: 4

> exp(confint(glmmod))
Waiting for profiling to be done...
                2.5 %    97.5 %
(Intercept) 0.1669932 0.9920023
age         0.9996431 1.0496390
female      0.3885283 0.9522567
> library(MASS)
> exp(confint.default(glmmod))
                2.5 %    97.5 %
(Intercept) 0.1683396 0.9965331
age         0.9995114 1.0493877
female      0.3909104 0.9563045

Why are they different? Which one is correct?

SAS

Fortunately the detailed documentation in SAS can help resolve this. The logistic procedure (section 4.1.1) offers the clodds option to the model statement. Setting this option to both produces two sets of CL, based on the Wald test and on the profile-likelihood approach. (Venzon, D. J. and Moolgavkar, S. H. (1988), “A Method for Computing Profile-Likelihood Based Confidence Intervals,” Applied Statistics, 37, 87–94.)


ods output cloddswald = waldcl cloddspl = plcl;
proc logistic data = "c:\book\help.sas7bdat"  plots=none;
class female (param=ref ref='0');
model homeless(event='1') = age female / clodds = both;
run;

 Odds Ratio Estimates and Profile-Likelihood Confidence Intervals

 Effect                Unit     Estimate     95% Confidence Limits

 AGE                 1.0000        1.024        1.000        1.050
 FEMALE 1 vs 0       1.0000        0.611        0.389        0.952


        Odds Ratio Estimates and Wald Confidence Intervals

 Effect                Unit     Estimate     95% Confidence Limits

 AGE                 1.0000        1.024        1.000        1.049
 FEMALE 1 vs 0       1.0000        0.611        0.391        0.956

Unfortunately, the default precision of the printout isn't quite sufficient to identify whether this distinction aligns with the differences seen in the two R methods. We get around this by using the ODS system to save the output as data sets (section A.7.1). Then we can print the data sets, removing the default rounding formats to find all of the available precision.


title "Wald CL";
proc print data=waldcl; format _all_; run;
title "PL CL";
proc print data=plcl; format _all_; run;

                              Wald CL  
                                     Odds
   Obs    Effect           Unit    RatioEst    LowerCL    UpperCL

    1     AGE                1      1.02415    0.99951    1.04939
    2     FEMALE 1 vs 0      1      0.61143    0.39092    0.95633
 
  
                               PL CL                            
                                     Odds
   Obs    Effect           Unit    RatioEst    LowerCL    UpperCL

    1     AGE                1      1.02415    0.99964    1.04964
    2     FEMALE 1 vs 0      1      0.61143    0.38853    0.95226

With this added precision, we can see that the confint.default() function in the MASS library generates the Wald confidence limits, while the confint() function produces the profile-likelihood limits. This also explains the confint() comment "Waiting for profiling to be done..." Thus neither CI from the MASS library is incorrect, though the profile-likelihood method is thought to be superior, especially for small sample sizes. Little practical difference is seen here.

Example 9.13: Negative binomial regression with proc mcmc

2011-11-08T09:10:00.010-05:00

In practice, data that derive from counts rarely seem to be fit well by a Poisson model; one more flexible alternative is a negative binomial model. In this SAS-only entry, we discuss how proc mcmc can be used for estimation. An overview of support for Bayesian methods in R can be found in the Bayesian Task View.

SAS

As noted in example 8.30, the SAS rand function lacks the option to input the mean directly, instead using the basic parameters of the probability of success and the number of successes k. (Though note the negative binomial has several formulations, which can cause problems when using multiple software systems.) As developed in that example, we use the the proc fcmp function to instead work with the mean.


proc fcmp outlib=sasuser.funcs.test;
function poismean_nb(mean, size);
  return(size/(mean+size));
  endsub;
run;

options cmplib=sasuser.funcs;
run;

With that preparation out of the way, we simulate some data--here an intercept of 0 and a slope of 1.


data test;
  do i = 1 to 10000;
    x = normal(0);
 mu = exp(0 + x);
 k = 2;
 y = rand("NEGBINOMIAL", poismean_nb(mu, k),k);
 output;
 end;
run;

The proc mcmc code presents a slight difficulty: the k successes before the random number of failures ought to be an integer, and proc mcmc appears to lack an integer-valued distribution. The model will run with continuous values of k, but its behavior is strange. Instead, we put a prior on a new parameter, kstar and take k as the rounded value (section 1.8.4) of kstar; since the values must be > 0, we also add 1 to the rounded value.


proc mcmc data=test nmc=1000 thin=1 seed=10061966;
parms beta0 1 beta1 1 kstar 10;

prior b: ~ normal(0, var = 10000);
prior kstar ~ igamma(.01, scale=0.01);

k=round(kstar+1, 1);
mu = exp(beta0 + beta1 * x);

model y ~ negbin(k, poismean_nb(mu, k));
run;

The way the kstar and k business works is that SAS actually processes the programming statements in each iteration of the chain. Posterior summaries just below, sample diagnostic plot above.


                       Posterior Summaries

Parameter        N  Mean  Standard              Percentiles
                          Deviation         25%      50%     75%
beta0        10000 0.00712 0.0131        -0.00171  0.00721 0.0156
beta1        10000 0.9818  0.0128         0.9732   0.9814  0.9905
kstar        10000 0.9648  0.2855         0.7112   0.9481  1.1974

                       Posterior Intervals
Parameter Alpha Equal-Tail Interval   HPD Interval
beta0    0.050 -0.0195 0.0321       -0.0182 0.0328
beta1    0.050  0.9569 1.0074        0.9562 1.0063
kstar    0.050  0.5208 1.4709        0.5001 1.4348

If a simple model like the one shown here is all you need, proc genmod's bayes statement can work for you. But the formulation demonstrated above would be useful for a generalized linear mixed model, for example.

Example 9.12: simpler ways to carry out permutation tests

2011-10-31T10:23:00.012-04:00

In a previous entry, as well as section 2.4.3 of the book, we describe how to carry out a 2 group permutation test in SAS as well as with the coin package in R. We demonstrate with comparing the ages of the female and male subjects in the HELP study.

In this entry, we revisit the permutation test using other functions.

R

We describe a simpler interface to carry out and visualize permutation tests using the functions from the mosaic package.

We begin by replicating our previous example (section 2.6.4, p. 87).


ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
library(coin)
numsim = 1000
oneway_test(age ~ as.factor(female), 
  distribution=approximate(B=numsim-1), data=ds)

which yields the following output:


 Approximative 2-Sample Permutation Test

data:  age by as.factor(female) (0, 1) 
Z = -0.9194, p-value = 0.3894
alternative hypothesis: true mu is not equal to 0

We conclude that there is minimal evidence to contradict the null hypothesis that the two groups have the same ages back in their respective populations.

Now we demonstrate another way to run this test in a more general form, using the mosaic package's do() function combined with the * operator to repeatedly carry out fitting a linear model with a parameter for female which will calculate our test statistic (difference in means between females and males) repeatedly after shuffling the group indicators. The shuffle() function permutes the group labels, and then the summary statistic is calculated.


> library(mosaic)
> obsdiff = with(ds, mean(age[female==1]) - mean(age[female==0]))
> obsdiff
     mean 
0.7841284 
> summary(age ~ female, data=ds, fun=mean)
age    N=453

+-------+---+---+--------+
|       |   |N  |mean    |
+-------+---+---+--------+
|female |No |346|35.46821|
|       |Yes|107|36.25234|
+-------+---+---+--------+
|Overall|   |453|35.65342|
+-------+---+---+--------+

Now we can run the permutation test, then display the results on a souped-up histogram with different shading for values larger in magnitude than the observed statistic (see above).


res = do(numsim) * lm(age ~ shuffle(female), data=ds)
pvalue = sum(abs(res$female) > abs(obsdiff)) / numsim
xhistogram(~ female, groups = abs(female) > abs(obsdiff), 
  n=50, density=TRUE, data=res, xlab="difference between groups",
  main=paste("Permutation test result: p=", round(pvalue, 3)))

The results are similar to those from the previous test: there is little evidence to contradict the null hypothesis.

SAS

In SAS, we'll take another approach, delving into the capabilities of proc iml to make a manual permutation test. We begin by reading the data and replicating the example in the book.


libname k 'c:\book';
proc npar1way data = k.help;
class female;
var age;
exact scores=data / mc n= 9999 alpha = .05;
run;

                   Data Scores One-Way Analysis

                  Chi-Square               0.8453
                  DF                            1
                  Pr > Chi-Square          0.3579

Permuting data is a very awkward thing to do in data steps. But it turns out to be easy in proc iml (the built-in SAS matrix language). Here we read in the key variables from the data set (use and read). Then we generate the permutations (ranperm). However, this generates row for each permuted data set, while we need a column for each, so we transpose the matrix (t) before saving it. Then we save the resulting data in a SAS data set with the female variable. Note that we permuted the ages only, as opposed to the R example-- it doesn't matter which is permuted, of course. Much of the proc iml code used here can be found in section 1.9 of the book-- however, note that curly braces are required in the read statement, as shown below.


proc iml;
use k.help;
read all var{female age} into x;
p = t(ranperm(x[,2],1000));
paf = x[,1]||p;
create newds from paf;
append from paf;
quit;

With the permuted data in hand, we use proc ttest (section 2.4.1) with the ODS system to generate and save the differences. Note that the default variable names from proc iml are fairly nondescript. With the 1000 permuted statistics in hand, we can generate a histogram of the statistics and a p-value with proc univariate.


ods output conflimits=diff;
proc ttest data=newds plots=none;
  class col1;
  var col2 - col1001;
run;

proc univariate data=diff;
  where method = "Pooled";
  var mean;
  histogram mean / normal;
run;

data diff2;
set diff;
absdiff = abs(mean);
run;

proc univariate data=diff2
  loccount mu0 = 0.7841284;
  where method = "Pooled";
  var absdiff;
run;

                     Location Counts: Mu0=0.78

                     Count                Value

                     Num Obs > Mu0          357
                     Num Obs ^= Mu0        1000
                     Num Obs < Mu0          643

Proc tabulate for simple statistics (corrected)

2011-10-30T08:43:00.007-04:00

Ken Beath, of Macquarie University, commented on an earlier entry that the best way to generate summary statistics is using proc tabulate. While the best tools might differ, depending on the purpose, we wanted to share Ken's code demonstrating how to replicate the R mosaic package tables using proc tabulate.

SAS

Ken's fully annotated code is appended below; we highlight the key syntax elements here. Reading in the data is shown in many examples.


proc tabulate data=help;
 class substance;
 var cesd;
 table (substance all),cesd*(n nmiss mean median);
run;

The class and var statements serve their usual purpose of identifying the categorical and analysis variables. The table statement does the work of the procedure, specifying a table with rows (requested before the comma) for each level of substance and overall, with columns (requested after the the *) to include the listed statistics for the analysis variable cesd. The resulting table is shown below.


    ------------------------------------------------------------------------
    |                  |                       cesd                        |
    |                  |---------------------------------------------------|
    |                  |     N      |   NMiss    |    Mean    |   Median   |
    |------------------+------------+------------+------------+------------|
    |substance         |            |            |            |            |
    |------------------|            |            |            |            |
    |alcohol           |      177.00|        0.00|       34.37|       36.00|
    |------------------+------------+------------+------------+------------|
    |cocaine           |      152.00|        0.00|       29.42|       30.00|
    |------------------+------------+------------+------------+------------|
    |heroin            |      124.00|        0.00|       34.87|       35.00|
    |------------------+------------+------------+------------+------------|
    |All               |      453.00|        0.00|       32.85|       34.00|
    ------------------------------------------------------------------------

Below we show Ken's code, complete with his helpful annotations. Note his use of formats. We're not fond of formats, but for presentation, as opposed to analysis, they can be very useful.

(Note: a prior version of this entry inexplicably referred to proc report, rather than proc tabulate.)


PROC IMPORT OUT= help
  DATAFILE= "C:\Users\kbeath\Documents\tabulate\help.csv" 
  DBMS=CSV REPLACE;
  GETNAMES=YES;
  DATAROW=2; 
RUN;

/* missing option create a category missing for each categorical 
    variable, always a good idea;
  the table statement specifies row then column;
  so for this example we have substance defining the rows, and 
    cesd statistics the columns */

proc tabulate data=help missing;
 class substance;
 var cesd;
 table substance,cesd*(mean n nmiss);
run;

/* formchar specifies the characters used to form the borders 
      - we set them all to blank to have no borders;
   mean*f=8.2 specifies that the mean is formatted using 
      an 8.2 format, etc*/

proc tabulate data=help missing formchar='           ';
 class substance;
 var cesd;
 table substance,cesd*(mean*f=8.2 n*f=8. nmiss*f=8.);
run;


/* substance="Substance" causes change in label to Substance */

proc format;
 value $subf "alcohol"="Alcohol" 
          "cocaine"="Cocaine" "heroin"="Heroin";

proc tabulate data=help missing formchar='           ';
 class substance;
 var cesd;
 format substance $subf.;
 table substance="Substance",
            cesd="CESD"*(mean*f=8.2 n*f=8. nmiss*f=8.);
run;

/* all the statistics. I've changed the format to 7.2 
     so they all fit on a line */

proc tabulate data=help missing formchar='           ';
 class substance;
 var cesd;
 format substance $subf.;
 table substance="Substance",
            cesd="CESD"*(n*f=8. nmiss*f=8. 
            (mean std min q1 median q3 max)*f=7.2);
run;

/* add a line for all */

proc tabulate data=help missing formchar='           ';
 class substance;
 var cesd;
 format substance $subf.;
 table (substance="Substance" all),
            cesd="CESD"*(n*f=8. nmiss*f=8. 
            (mean std min q1 median q3 max)*f=7.2);
run;

/* to show how easy it is, further subdivide 
   by racial group */

proc tabulate data=help missing formchar='           ';
 class substance racegrp;
 var cesd;
 format substance $subf.;
 table (racegrp all)*(substance="Substance" all),
            cesd="CESD"*(n*f=8. nmiss*f=8. 
            (mean std min q1 median q3 max)*f=7.2);
run;

Example 9.11: Employment plot

2011-10-25T14:30:00.027-04:00

A facebook friend posted the picture reproduced above-- it makes the case that President Obama has been a successful creator of jobs, and also paints GW Bush as a president who lost jobs. Another friend pointed out that to be fair, all of Bush's presidency ought to be included. Let's make a fair plot of job growth and loss. Data can be retrieved from the Bureau of Labor Statistics, where Nick will be spending his next sabbatical. The extract we use below is also available from the book website. This particular table reports the cumulative change over the past three months, adjusting for seasonal trends. This tends to smooth out the line.

SAS

The first job is to get the data into SAS. Here we demonstrate reading it directly from a URL, as outlined in section 1.1.6.


filename myurl
   url "http://www.math.smith.edu/sasr/datasets/bls.csv";
       
data q_change;
   infile myurl delimiter=',';
input  Year Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Annual;
run;

The raw data are in a pretty inconvenient format for plotting. To make a long, narrow data set with a row for each month, we'll use proc transpose (section 1.5.3) to flip each year on its side. Then, to attach a date to each measure, we'll use the compress function. First we add "01" (the first of the month) to the month name, which is in a variable created by proc transpose with the default name "_name_". Then we tack on the year variable, and input the string in the date format. The resulting variable is a SAS date (number of days since 12/31/1959, see section 1.6.1).


proc transpose data=q_change out=q2;
  by year;
run;

data q3;
set q2;
  date1 = input(compress("01"||_name_||year),date11.);
run;

Now the data are ready to plot. It would probably be possible to use proc sgplot but proc gplot is more flexible and allows better control for presentation graphics.


title "3-month change in private-sector jobs, seasonally adjusted";
axis1 minor = none label = (h=2 angle = 90 "Thousands of jobs")
  value = (h = 2);
axis2 minor = none value = (h=2)label = none 
  offset = (1cm, -5cm)
  reflabel = (h=1.5 "Truman" "Eisenhower" "Kennedy/Johnson" 
    "Nixon/Ford" "Carter" "Reagan" "GHW Bush" "Clinton" "GW Bush" "Obama" );
symbol1 i=j v=none w=3;

proc gplot data=q3;
plot col1 * date1 / vaxis=axis1 haxis=axis2 vref=0 
  href = '12apr1945'd '21jan1953'd '20jan1961'd '20jan1969'd 
    '20jan1977'd '21jan1981'd '20jan1989'd '20jan1993'd 
    '20jan2001'd '20jan2009'd;
format date1 monyy6.;
run;
quit;

Much of the syntax above has been demonstrated in our book examples and blog entries. What may be unfamiliar is the use of the href option in the plot statement and the reflabel option in the axis statement. The former draws reference lines at the listed values in the plot, while the latter adds titles to these lines. The resulting plot is shown here.

Looking fairly across postwar presidencies, only the Kennedy/Johnson and Clinton years were mostly unmarred by periods with large losses in jobs. The Carter years were also times jobs were consistently added. While the graphic shared on facebook overstates the case against GW Bush, it fairly shows Obama as a job creator thus far, to the extent a president can be credited with jobs created on his watch.

R

The main trick in R is loading the data and getting it into the correct format.
here we use cbind() to grab the appropriate columns, then transpose that matrix and turn it into a vector which serves as input for making a time series object with the ts() command (as in section 4.2.8). Once this is created, the default plot for a time series object is close to what we have in mind.


ds = read.csv("http://www.math.smith.edu/sasr/datasets/bls.csv", 
              header=FALSE)
jobs = with(ds, cbind(V2, V3, V4, V5, V6, V7, V8, V9, V10,
                      V11, V12, V13))
jobsts = ts(as.vector(t(jobs)), start=c(1945, 1), 
            frequency=12)
plot(jobsts, plot.type="single", col=4,
     ylab="number of jobs (in thousands)")

All that remains is to add the reference lines for 0 jobs and the presidencies. The lines are most easily added with the abline() function (section 5.2.1). Easier than adding labels for the lines within the plot function will be to use the mtext() function to place the labels in the margins. We'll write a little function to save a few keystrokes by plotting the line and adding the label together.


abline(h=0)
presline = function(date,line,name){
  mtext(at=date,text= name, line=line)
  abline(v = date)
}
presline(1946,1,"Truman")
presline(1953,2,"Eisenhower")
presline(1961,1,"Kennedy/Johnson")
presline(1969,2,"Nixon/Ford")
presline(1977,1,"Carter")
presline(1981,2,"Reagan")
presline(1989,1,"GHW Bush")
presline(1993,2,"Clinton")
presline(2001,1,"GW Bush")
presline(2009,2,"Obama")

It might be worthwhile to standardize the number of jobs to the population size, since the dramatic loss of jobs due to demobilization after the Second World War during a single month in 1945 (2.4 million) represented 1.7% of the population, while the recent loss of 2.3 million jobs in 2009 represented only 0.8% of the population.

Example 9.10: more regression trees and recursive partitioning with "partykit"

2011-10-17T11:08:00.000-04:00

We discuss recursive partitioning, a technique for classification and regression using a decision tree in section 6.7.3 of the book. Support for these methods is available within the rpart package. Torsten Hothorn and Achim Zeileis have extended the support for these methods with the partykit package, which provides a toolkit with infrastructure for representing, summarizing, and visualizing tree-structured regression and classification models.

In this entry, we revisit the example from the book, which worked to classify predictors of homelessness in the HELP study.

R


ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
library(rpart); library(partykit)
ds$sub = as.factor(ds$substance)
homeless.rpart = rpart(homeless ~ female + i1 + sub + sexrisk + mcs +
  pcs, method="class", data=ds)
plot(homeless.rpart)
text(homeless.rpart)

This reproduces Figure 6.2 (p. 236) from the book, while we can display the output from the classification tree using the printcp() command.


> printcp(homeless.rpart)
Classification tree:
rpart(formula = home ~ female + i1 + sub + sexrisk + mcs + pcs, 
    data = ds, method = "class")
Variables actually used in tree construction:
[1] female  i1      mcs     pcs     sexrisk

Root node error: 209/453 = 0.5
n= 453 
    CP nsplit rel error xerror xstd
1 0.10      0       1.0    1.0 0.05
2 0.05      1       0.9    1.1 0.05
3 0.03      4       0.8    1.1 0.05
4 0.02      5       0.7    1.0 0.05
5 0.01      7       0.7    0.9 0.05
6 0.01      9       0.7    0.9 0.05

Using the partykit package, we can make a nice graphic describing these results. We'll use the plot.party() function on a party object (coerced from the rpart object generated above using as.party()). This provides more information about the tree (as seen in the Figure above).


plot(as.party(homeless.rpart), type="simple")

More information as well as a lovely vignette can be found here.

SAS

Recursive partitioning is available through SAS Enterprise Miner, a module not always included in SAS installations.

Example 9.9: Simplifying R using the mosaic package (part 1)

2011-10-13T11:35:00.009-04:00

While both SAS and R are powerful systems for statistical analysis, they can be frustrating to new users or those learning statistics for the first time.

R
The mosaic package is designed to help simplify the interface for such new users, while allowing them to undertake sophisticated analyses.

As an example of how the package simplifies life for the novice user, consider calculating summary statistics and displaying a densityplot for the CESD (measure of depressive symptom) scores by substance abuse group in the HELP dataset. Doing this in R without the package would require mastering a package such as plyr to replicate results by substance or a typing-intensive use of syntax to select rows corresponding to each substance.


ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
library(mosaic)
options(digits=3)

After loading the data and the package, and setting the number of digits to a more reasonable default, we can call the mean() function to easily calculate this statistic (denoted by S in the result) for each of the three substance abuse groups alcohol, cocaine or heroin.


> mean(cesd ~ substance, data=ds)
  substance    S   N Missing
1   alcohol 34.4 177       0
2   cocaine 29.4 152       0
3    heroin 34.9 124       0

Similar results are seen when we calculate the standard deviations per group:


> sd(cesd ~ substance, data=ds)
  substance    S   N Missing
1   alcohol 12.1 177       0
2   cocaine 13.4 152       0
3    heroin 11.2 124       0

Another function can calculate a raft of summary statistics for each group that are nicely formatted.


> summary(cesd ~ substance, data=ds, fun=favstats)
cesd    N=453
+---------+-------+---+----+---+-------+---+----+-----+----+---+--------+
|         |       |N  |min |Q1 |median |Q3 |max |mean |sd  |n  |missing |
+---------+-------+---+----+---+-------+---+----+-----+----+---+--------+
|substance|alcohol|177|4   |26 |36     |42 |58  |34.4 |12.1|177|0       |
|         |cocaine|152|1   |19 |30     |39 |60  |29.4 |13.4|152|0       |
|         |heroin |124|4   |28 |35     |43 |56  |34.9 |11.2|124|0       |
+---------+-------+---+----+---+-------+---+----+-----+----+---+--------+
|Overall  |       |453|1   |25 |34     |41 |60  |32.8 |12.5|453|0       |
+---------+-------+---+----+---+-------+---+----+-----+----+---+--------+

These commands allow quick review of the data to ensure, for example, that assumptions of equal variance are justified, or that coding errors or missing values haven't crept in.

A graphical depiction using a set of densityplots (shown above) can be created using the command:


densityplot(~ cesd, group=substance, data=ds, auto.key=TRUE)

SAS
We're unaware of any similar program that attempts to simplify SAS syntax for educational use. To replicate the above results, we would use the means and sgpanel procedures.


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

options ls=80;
proc means data=ds fw=4
  min q1 median q3 max mean std nmiss n;
  class substance;
  var cesd;
run;
                     Analysis Variable : CESD

              N           Lower             Upper               Std
 SUBSTANCE  Obs   Min  Quartile  Median  Quartile   Max  Mean   Dev
 ------------------------------------------------------------------
 alcohol    177  4.00      26.0    36.0      42.0  58.0  34.4  12.1
 cocaine    152  1.00      19.0    30.0      39.0  60.0  29.4  13.4
 heroin     124  4.00      28.0    35.0      43.0  56.0  34.9  11.2
 ------------------------------------------------------------------

                                 N     N
                    SUBSTANCE  Obs  Miss      N
                    ---------------------------
                    alcohol    177     0    177
                    cocaine    152     0    152
                    heroin     124     0    124
                    ---------------------------

After reading the data in, the meansprocedure can produce any of the desired statistics (plus may others) directly. To replicate the mosaic package in printing a single statistic, list only that statistic in the proc means statement. Note that the overall statistic in the R table is not included. To replicate that row, you would re-run the above code, omitting the class statement.

To the best of our knowledge, there still does not exist an easy way to plot multiple densities in a single SAS plot. In example 2.6.4 we show how it can be done using proc kde, saving the density estimates, and plotting separately. (Code for this is included at the book web site.) But in the interest of simple code, we show a simpler method using proc sgpanel. The result, show below, is less useful than the R plot from the the mosaic package, but still gets the point across.


proc sgpanel data=ds;
  panelby substance / columns=1;
  density cesd / type=kernel;
run;

Example 9.8: New stuff in SAS 9.3-- Bayesian random effects models in Proc MCMC

2011-10-04T10:16:00.006-04:00

Rounding off our reports on major new developments in SAS 9.3, today we'll talk about proc mcmc and the random statement.

Stand-alone packages for fitting very general Bayesian models using Markov chain Monte Carlo (MCMC) methods have been available for quite some time now. The best known of these are BUGS and its derivatives WinBUGS (last updated in 2007) and OpenBUGS . There are also some packages available that call these tools from R.

Today we'll consider a relatively simple model: Clustered Poisson data where cluster means are a constant plus a cluster-specific exponentially-distributed random effect. To be clear:
y_ij ~ Poisson(mu_i)
log(mu_i) = B_0 + r_i
r_i ~ Exponential(lambda)
Of course in Bayesian thinking all effects are random-- here we use the term in the sense of cluster-specific effects.

SAS
Several SAS procedures have a bayes statement that allow some specific models to be fit. For example, in Section 6.6 and example 8.17, we show Bayesian Poisson and logistic regression, respectively, using proc genmod. But our example today is a little unusual, and we could not find a canned procedure for it. For these more general problems, SAS has proc mcmc, which in SAS 9.3 allows random effects to be easily modeled.

We begin by generating the data, and fitting the naive (unclustered) model. We set B_0 = 1 and lambda = 0.4. There are 200 clusters of 10 observations each, which we might imagine represent 10 students from each of 200 classrooms.


data test2;
truebeta0 = 1;
randscale = .4;
call streaminit(1944);
  do i = 1 to 200;
    randint = rand("EXPONENTIAL") * randscale;
    do ni = 1 to 10;
      mu = exp(truebeta0  + randint); 
      y = rand("POISSON", mu);
      output;
    end;
  end;
run;

proc genmod data = test2;
model y = / dist=poisson;
run;

                      Standard       Wald 95%     
Parameter  Estimate     Error   Confidence Limits

Intercept    1.4983    0.0106    1.4776    1.5190

Note the inelegant SAS syntax for fitting an intercept-only model. The result is pretty awful-- 50% bias with respect to the global mean. Perhaps we'll do better by acknowledging the clustering. We might try that with normally distributed random effects in proc glimmix.


proc glimmix data = test2 method=laplace;
class i;
model y = / dist = poisson solution;
random int / subject = i type = un;
run;

   Cov                               Standard
   Parm       Subject    Estimate       Error
   UN(1,1)    i            0.1682     0.01841

                       Standard
Effect     Estimate     Error  t Value  Pr > |t|
Intercept    1.3805   0.03124    44.20    <.0001

No joy-- still a 40% bias in the estimated mean. And the variance of the random effects is biased by more than 50%! Let's try fitting the model that generated the data.


proc mcmc data=test2 nmc=10000 thin=10 seed=2011;
parms fixedint 1 gscale 0.4;

prior fixedint ~ normal(0, var = 10000);
prior gscale ~ igamma(.01 , scale = .01 ) ;

random rint ~ gamma(shape=1, scale=gscale) subject = i initial=0.0001;
mu = exp(fixedint + rint);
model y ~ poisson(mu);
run;

The key points of the proc mcmc statement are nmc, the total number of Monte Carlo iterations to perform, and thin, which includes only every nth sample for inference. The prior and model statements are fairly obvious; we note that in more complex models, parameters that are listed within a single prior statement are sampled as a block. We're placing priors on the fixed (shared) intercept and the scale of the exponential. The mu line is actually just a programming statement-- it uses the same syntax as data step programming.
The newly available statement is random. The syntax here is similar to those for the other priors, with the addition of the subject option, which generates a unique parameter for each level of the subject variable. The random effects themselves can be used in later statements, as shown, to enter into data distributions. A final note here is that the exponential distribution isn't explicitly available, but since the gamma distribution with shape fixed at 1 defines the exponential, this is not a problem. Here are the key results.


           Posterior Summaries

                               Standard        
  Parameter        N     Mean Deviation  
  fixedint      1000   1.0346    0.0244  
  gscale        1000   0.3541    0.0314  

           Posterior Intervals

 Parameter    Alpha        HPD Interval
 fixedint     0.050      0.9834      1.0791
 gscale       0.050      0.2937      0.4163

The 95% HPD regions include the true values of the parameters and the posterior means are much less biased than in the model assuming normal random effects.

As usual, MCMC models should be evaluated carefully for convergence and coverage. In this example, I have some concerns (see default diagnostic figure above) and if it were real data I would want to do more.

R
The CRAN task view on Bayesian Inference includes a summary of tools for general and model-specific MCMC tools. However, there is nothing like proc mcmc in terms of being a general and easy to use tool that is native to R. The nearest options are to use R front ends to WinBUGS/OpenBUGS (R2WinBUGS) or JAGS (rjags). (A brief worked example of using rjags was posted last year by John Myles White.) Alternatively, with some math and a little sweat, the mcmc package would also work. We'll explore an approach through one or more of these packages in a later entry, and would welcome a collaboration from anyone who would like to take that on.

Example 9.7: New stuff in SAS 9.3-- Frailty models

2011-09-27T10:14:00.006-04:00

Shared frailty models are a way of allowing correlated observations into proportional hazards models. Briefly, instead of l_i(t) = l_0(t)e^(x_iB), we allow l_ij(t) = l_0(t)e^(x_ijB + g_i), where observations j are in clusters i, g_i is typically normal with mean 0, and g_i is uncorrelated with g_i'. The nomenclature frailty comes from examining the logs of the g_i and rewriting the model as l_ij(t) = l_0(t)u_i*e^(xB) where the u_i are now lognormal with median 0. Observations j within cluster i share the frailty u_i, and fail faster (are frailer) than average if u_i > 1.

In SAS 9.2, this model could not be fit, though it is included in the survival package in R. (Section 4.3.2) With SAS 9.3, it can now be fit. We explore here through simulation, extending the approach shown in example 7.30.

SAS
To include frailties in the model, we loop across the clusters to first generate the frailties, then insert the loop from example 7.30, which now represents the observations within cluster, adding the frailty to the survival time model. There's no need to adjust the censoring time.


data simfrail;
  beta1 = 2;
  beta2 = -1;
  lambdat = 0.002; *baseline hazard;
  lambdac = 0.004; *censoring hazard;
  do i = 1 to 250; *new frailty loop;
    frailty = normal(1999) * sqrt(.5);
    do j = 1 to 5; *original loop;
      x1 = normal(0);
      x2 = normal(0);
      * new model of event time, with frailty added;
      linpred = exp(-beta1*x1 - beta2*x2 + frailty);
      t = rand("WEIBULL", 1, lambdaT * linpred);
        * time of event;
      c = rand("WEIBULL", 1, lambdaC);
        * time of censoring;
      time = min(t, c);    * which came first?;
      censored = (c lt t);
    output;
 end;
  end;
run;

For comparison's sake, we replicate the naive model assuming independence:


proc phreg data=simfrail;
  model time*censored(1) = x1 x2;
run;

               Parameter   Standard                         Hazard
 Parameter DF   Estimate      Error Chi-Square Pr > ChiSq    Ratio

 x1         1    1.68211    0.05859   824.1463     <.0001    5.377
 x2         1   -0.88585    0.04388   407.4942     <.0001    0.412

The parameter estimates are rather biased. In contrast, here is the correct frailty model.


proc phreg data=simfrail;
  class i;
  model time*censored(1) = x1 x2;
  random i / noclprint;
run;
                    Cov         REML    Standard
                    Parm    Estimate       Error

                    i         0.5329     0.07995

               Parameter   Standard                         Hazard
 Parameter DF   Estimate      Error Chi-Square Pr > ChiSq    Ratio

 x1         1    2.03324    0.06965   852.2179     <.0001    7.639
 x2         1   -1.00966    0.05071   396.4935     <.0001    0.364

This returns estimates gratifyingly close to the truth. The syntax of the random statement is fairly straightforward-- the noclprint option prevents printing all the values of i. The clustering variable must be specified in the class statement. The output shows the estimated variance of the g_i.

R
In our book (section 4.16.14) we show an example of fitting the uncorrelated data model, but we don't display a frailty model. Here, we use the data generated in SAS, so we omit the data simulation in R. As in SAS, it would be a trivial extension of the code presented in example 7.30. For parallelism, we show the results of ignoring the correlation, first.


> library(survival)
> with(simfrail, coxph(formula = Surv(time, 1-censored) ~ x1 + x2))

     coef exp(coef) se(coef)     z p
x1  1.682     5.378   0.0586  28.7 0
x2 -0.886     0.412   0.0439 -20.2 0

with identical results to above. Note that the Surv function expects an indicator of the event, vs. SAS expecting a censoring indicator.

As with SAS, the syntax for incorporating the frailty is simple.


> with(simfrail, coxph(formula = Surv(time, 1-censored) ~ x1 + x2 
    + frailty(i)))

            coef se(coef) se2    Chisq DF  p
x1          2.02 0.0692   0.0662 850     1 0
x2         -1.00 0.0506   0.0484 393     1 0
frailty(i)                       332   141 0

Variance of random effect= 0.436

Here, the results differ slightly from the SAS model, but still return parameter estimates that are quite similar. We're not familiar enough with the computational methods to diagnose the differences.

Example 9.6: Model comparison plots (Completed)

2011-09-21T09:30:00.001-04:00

We often work in settings where the data set has a lot of missing data-- some missingness in the (many) covariates, some in the main exposure of interest, and still more in the outcome. (Nick describes this as "job security for statisticians").

Some analysts are leery of imputing anything at all, preferring to rely on the assumption that the data are missing completely at random. Others will use multiple imputation for covariates, but feel they should use "real" data for the exposure and outcome. Still others will impute the exposure but not the outcome. Theory and experiments suggest (Moons et al JCE 2006) that all missing data should be imputed. Depending on the imputation method, this may offer some protection against missing data that missing at random, more general than missing completely at random.

In one analysis, we decided to use each of these approaches and demonstrate the results that would be obtained. The data are shown below. The first column denotes the data used, the second has the effect on the mean and CI limits for the effect. How can we present these results clearly? We designed a graphic that requires some customization using either SAS or R but which makes the point elegantly.

SAS
The SAS version is shown below. (Click on it for a larger image.) To generate it, add a final column to the data, where the effect estimate is repeated but the other values are not. Then a basic plot can be created in proc gplot with the hiloc interpolation in the symbol statement and the overlay option in the plot statement. (See book section 5.3 and other blog entries for details.) Try the code without the axis statements to see what happens.


data ke1;
input datatype estimate meanval;
cards;
1 .11 .11 
1 -.05 .
1 .28 .
2 .07 .07
2 .21 .
2 -.07 .
3 .06 .06
3 -.08 .
3 .2 .
4 0 0
4 -.13 .
4 .12 .
;;
cards;
run;

symbol1 i=hiloc c=black v=none;
symbol2 i=none v=dot h=1  c=black;
axis1 minor=none order = (1 to 4 by 1)
  value = (tick = 1 "Complete" 
             justify=c  "Case" justify = c "(N = 2055)"
           tick=2 "MI" justify=c  "Covariates only" 
             justify=c "(N = 2961)"
           tick=3 "MI" justify=c  "Covariates and exposure" 
             justify=c "(N = 3994)"
           tick=4 "MI" justify=c  "All variables" 
             justify=c "(N = 6782)"
           )
  label = none
  offset = (2 cm, 2 cm)
;
axis2 minor=none order = (-.2 to .3 by .1) 
  label = (angle=90 "Effect of exposure on outcome");
title "Compare missingness approaches";
proc gplot data = ke1;
plot (estimate meanval) * datatype / 
  overlay haxis=axis1 vref=0 vaxis=axis2;
run;
quit;

The two axis statements make the plot work. The axis1 statement uses the value option to hand-write the labels describing the data sets. Note that the justify = c causes a new line to be started. The offset option adds a little space to the left and right of the data. The axis2 statement specifies the range and label for the vertical axis. The extra symbol statement and the overlay option just plot the dots that call attention to the effect estimates-- otherwise they would show just a small crossbar at the effect.

The plot suggests that as more observations are included and the multiple imputation gains accuracy the effect attenuates and the standard errors decrease.

R

In R we create the equivalent plot in multiple steps, first by creating an empty plot of the correct size then iterating through each of the lines. As with the SAS approach, a little manipulation of the raw data is required.


n = c(2055, 2961, 3994, 6782)
labels = c("Complete Case", "MI\ncovariates only",
           "MI\ncovariates and exposure",
           "MI\nall variables")
est =    c(0.11,  0.07,  0.06,  0)
lower = c(-0.05, -0.07, -0.08, -0.13)
upper =  c(0.28,  0.21,  0.20,  0.12)

plot(c(0.5, 4.5), c(min(lower)-.10, max(upper)), type="n", xlab="", 
  xaxt="n", ylab="Effect of exposure on outcome")
title("Compare missingness approaches")
for (i in 1:length(n)) {
  points(i, est[i])
  lines(c(i,i), c(lower[i], upper[i]))
  stringval = paste(labels[i],"\n(N=",n[i],")")
  text(i, min(lower) - .05, stringval, cex=.6)
}
abline(h=0, lty=2)

The resulting plot is shown at the top. As opposed to the SAS approach, more of the figure can be defined using the data. For example, the y-axis values are determined from the minimum and maximum values to plot.

Note: a draft of this entry was published accidentally. Many apologies. --Ken

Example 9.5: New stuff in SAS 9.3-- proc FMM

2011-09-13T14:13:00.005-04:00

Finite mixture models (FMMs) can be used in settings where some unmeasured classification separates the observed data into groups with different exposure/outcome relationships. One familiar example of this is a zero-inflated model, where some observations come from a degenerate distribution with all mass at 0. In that case the exposure/outcome relationship is less interesting in the degenerate distribution group, but there would be considerable interest in the estimated probability of group membership. Another possibly familiar setting is the estimation of a continuous density as a mixture of normal distributions.

More generally, there could be several groups, with "concomitant" covariates predicting group membership. Each group might have different sets of predictors and outcomes from different distribution families. On the other hand, in a "homogenous" mixture setting, all groups have the same distributional form, but with different parameter values. If the covariates in the model are the same, this setting is similar to an ordinary regression model where every observed covariate is interacted with the (unobserved) group membership variable.

SAS

SAS 9.3 includes the "experimental" FMM procedure to fit these models. We're unsure what criteria SAS uses to decide when a procedure is experimental and when it becomes "production", but experimental procedures in SAS/STAT usually do become production eventually.

As hinted at above, the generality implies by the models is fairly vast, and the FMM procedure includes a lot of generality. The most obvious limitation is that it requires independent observations.

We'll demonstrate with a simulated data set. We create a variable x that predicts both group membership and an outcome y with different linear regression parameters depending on group. The mixing probability follows a logistic regression with intercept=-2 and slope=1. The intercept and slope for the outcome are (0, 1) and (3, 1.2) for groups 0 and 1, respectively. The resulting data is plotted above using a default plot from proc glm using the actual group membership. A histogram and nearest normal density for the residuals from a simple OLS are shown, with R code, below-- they appear more or less normal.


data fmmtest;
do i = 1 to 5000;
  x = normal(0);
  group = (exp(-1 + 2*x)/(1 + exp(-1 + 2*x))) 
    gt uniform(0);
  y = (group * 3) + ((1 + group/5) * x) + 
    normal(0) * sqrt(1);
  output;
  end;
run;

title "Fixed k = 2";
proc fmm data = fmmtest;
  model y = x / k=2 equate=scale;
  probmodel x;
run;

In the model statement, the option k defines how many groups (or "components") are to be included. There are also options kmin and kmax which will fit models with various numbers of components and report results of one of them based on some model criterion. The equate option allows the user to force some elements of the component distributions to be equal. Here we force the residuals to be equal, since that's the model that generated our data. The probmodel statement is how the concomitant variables enter the model. The default settings model a logistic (or generalized logit) regression using the listed covariates.

We show only the parameter estimates.


                                Standard
Component  Effect     Estimate     Error  z Value  Pr > |z|
        1  Intercept    2.9856   0.04856    61.48    <.0001
        1  x            1.2060   0.03927    30.71    <.0001
        2  Intercept  -0.01978   0.02512    -0.79    0.4309
        2  x            0.9889   0.02453    40.32    <.0001
        1  Variance     1.0279   0.02531
        2  Variance     1.0279   0.02531


            Parameter Estimates for Mixing Probabilities
                         Standard
Effect       Estimate       Error    z Value    Pr > |z|
Intercept     -1.0018     0.05665     -17.68      <.0001
x              1.9403     0.07413      26.17      <.0001

The recovery of the true parameters is quite good. This may be unsurprising given the sample size and the distinctness of the components, but we think it's pretty impressive.

R
The CRAN Task View on Cluster Analysis & Finite Mixture Models provides an overview of packages available for R.

First, we'll make the data and take a look at the simple model diagnostics.


> x = rnorm(5000)
> probgroup1 = exp(-1 + 2*x)/(1 + exp(-1 + 2*x))
> group = ifelse(probgroup1 > runif(5000),1,0)
> y = (group * 3) + ((1 + group/5) * x) + rnorm(5000);

> resids = residuals(lm(y~x))
> hist(resids, freq = FALSE)
> dvals = seq(from=min(resids), to=max(resids),length=100)
> lines(dvals, dnorm(dvals, mean(resids), sd(resids)))

The histogram of the residuals from the simple OLS model is shown below. They appear reasonably normal.

Ron Pearson recently showed a simple example using the mixtools package. We tried to replicate the example above using mixtools but were unable to find functionality for concomitant variables. If you fit the closest model, assuming the mixing probabilities depend on no covariates, this is what happens:


> library(mixtools)
> mixout.mt = regmixEM(y,x,k=2,arbvar = FALSE)> summary(mixout.mt)
summary of regmixEM object:
  comp 1   comp 2
lambda  0.5275915 0.472408
sigma   1.0553877 1.055388
beta1  -0.0035015 2.237736
beta2   1.5068146 2.004008

which is pretty awful.

Fortunately, the flexmix package, while a bit more difficult to use, offers the needed flexibility. As a side note, the generality of the package is pretty awe-inspiring. Nice work!


library(flexmix)
mixout.fm=flexmix(y~x, k=2, model=FLXMRglmfix(y~x, varFix=TRUE),
   concomitant=FLXPmultinom(~ x))

The flexmix function uses a variety of special objects that are created by other functions the package provides. Here we use the FLXMRglmfix function to force equal variances across the components and the FLXPmultinom function to define the logistic regression on the covariate x. The results are:


> parameters(mixout5)
                   Comp.1       Comp.2
coef.(Intercept) 3.000412 -0.008360171
coef.x           1.190454  1.047660263
sigma            0.972306  0.972306015
> parameters(mixout5, which = "concomitant")
            1          2
(Intercept) 0  0.9605367
x           0 -1.9109429

which also reproduces reality with impressive accuracy. Note that the concomitant variable model apparently predicts membership in the second component, so the signs are reversed from the generating model.

Example 9.4: New stuff in SAS 9.3-- MI FCS

2011-09-06T14:12:00.007-04:00

We begin the new academic year with a series of entries exploring new capabilities of SAS 9.3, and some functionality we haven't previously written about.

We'll begin with multiple imputation. Here, SAS has previously been limited to multivariate normal data or to monotonic missing data patterns.

SAS

SAS 9.3 adds the FCS statement to proc mi. This implements a fully conditional specification imputation method (e.g., van Buuren, S. (2007), "Multiple Imputation of Discrete and Continuous Data by Fully Conditional Specification," Statistical Methods in Medical Research, 16, 219–242.) Briefly, we begin by imputing all the missing data with a simple method. Then missing values for each variable are imputed using a model created with the real and current imputed values for the other variables, iterating across the variables several times.

We replicate the multiple imputation example from the book, section 6.5. In that example, we used the mcmc statement for imputation: at the time, this was the only method available in SAS when a non-monotonic missingness pattern was present. We noted at the time that this was not "strictly appropriate" since mcmc method assumes multivariate normality, and two of our missing variables were dichotomous.

filename myhm url "http://www.math.smith.edu/sasr/datasets/helpmiss.csv" lrecl=704;

proc import replace datafile=myhm out=help dbms=dlm;
delimiter=',';
getnames=yes;
run;

proc mi data = help nimpute=20 out=helpmi20fcs;
class homeless female;
var i1 homeless female sexrisk indtot mcs pcs;
fcs
 logistic (female)
 logistic (homeless);
run;

In the fcs statement, you list the method (logistic, discrim, reg, regpmm) to be used, naming the variable for which the method is to be used in parentheses following the method. (You can also specify a subset of covariates to be used in the method, using the usual SAS model-building syntax.) Omitted covariates are imputed using the default reg method.

ods output parameterestimates=helpmipefcs
covb = helpmicovbfcs;
proc logistic data=helpmi20fcs descending;
by _imputation_;
model homeless=female i1 sexrisk indtot /covb;
run;

proc mianalyze parms=helpmipefcs covb=helpmicovbfcs;
  modeleffects intercept female i1 sexrisk indtot;
run;

with the following primary result:

Parameter    Estimate   Std Error  95% Conf. Limits

intercept   -2.492733    0.591241  -3.65157  -1.33390 
female      -0.245103    0.244029  -0.72339   0.23319
i1           0.023207    0.005610   0.01221   0.03420
sexrisk      0.058642    0.035803  -0.01153   0.12882
indtot       0.047971    0.015745   0.01711   0.07883

which is quite similar to our previous results. Given the small proportion of missing values, this isn't very surprising.

R
Several R packages allow imputation for a general pattern of missingness and missing outcome distribution. A brief summary of missing data tools in R can be found in the CRAN Task view on Multivariate Statistics. We'll return to this topic from the R perspective in a future entry.

Taking August off!

2011-07-31T17:32:00.004-04:00

We'll be back with recharged batteries and lots of new entries in September. Have a great summer*!

As usual, please send any questions you have about using SAS or R.

*Not valid in the southern hemisphere.

Really useful R package: sas7bdat

2011-07-25T15:24:00.005-04:00

For SAS users, one hassle in trying things in R, let alone migrating, is the difficulty of getting data out of SAS and into R. In our book (section 1.2.2) and in a blog entry we've covered getting data out of SAS native data sets. Unfortunately, for all of these methods, you need a working, licensed version of SAS.

However Matt Shotwell has reverse-engineered the sas7bdat file format. This means that you can now read a SAS data set without a working copy of SAS. This is a wonderful thing, and in fact SAS Institute ought to have provided this ability long ago. The package is experimental, but it worked fine for two small data sets. Matt tells me that as of 7/2011, the package only works for sas7bdat files generated on 32-bit Windows systems.

R
Install the package sas7bdat. The use the read.sas7bdat() function.


library(sas7bdat)
helpfromSAS = read.sas7bdat("http://www.math.smith.edu
/sasr/datasets/help.sas7bdat")

(Note that newlines are not allowed in the URL in practice, but formatting for the blog required it.)


> is.data.frame(helpfromSAS)
[1] TRUE
> summary(helpfromSAS$MCS)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  6.763  21.680  28.600  31.680  40.940  62.180 
> with(helpfromSAS, summary(SUBSTANCE))
alcohol cocaine  heroin 
    177     152     124

It's unclear why all the variable names are all capitalized. That didn't happen in another trial, so it must be something about the way the help.sas7bdat data set was constructed.

Example 9.3: augmented display of contingency table

2011-07-18T11:10:00.017-04:00

SAS and R often provide different levels of details from output. This is particularly true for the descriptive analysis of contingency tables, where SAS makes it easy to display tables with additional quantities (such as the observed cell count).

The mosaic package has added functionality to calculate these quantities in R. We demonstrate using an example from the HELP dataset.

R


ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
library(mosaic)
ds$gender = ifelse(ds$female==1, "female", "male")
ds$homeless = ifelse(ds$homeless==1, "homeless", "housed")
tab = xtabs(~ gender + homeless, data=ds)
> tab
        homeless
gender   homeless housed
  female       40     67
  male        169    177
> xchisq.test(tab)

 Pearson's Chi-squared test with Yates' continuity correction

data:  tab 
X-squared = 3.8708, df = 1, p-value = 0.04913

  40.00    67.00 
( 49.37) ( 57.63)
 [1.78]   [1.52] 
<-1.33>  < 1.23> 
   
 169.00   177.00 
(159.63) (186.37)
 [0.55]   [0.47] 
< 0.74>  <-0.69> 
   
key:
 observed
 (expected)
 [contribution to X-squared]

We see that there is a borderline statistically significant association between gender and homeless status in the HELP study. We interpret that we see fewer than expected females who are homeless, and more males who are homeless.

Another idea is to use graphical depictions of the association in this table. One approach is a mosaic plot (note: no relation to Project MOSAIC and the mosaic package). A mosaic plot starts as a square with area equal to one. It is divided into columns based on the prevalence in each of the values for the column variable (in this case, gender). Then each bar is divided vertically based on the conditional probability of the other variable within that category.

Another graphical display of a table is the association plot. In an association plot, there is also a box for each cell of the table. The area of the box is proportional to the difference between the observed and expected (assuming no association) frequencies. In a typical presentation, excess observed counts are black and above the line, while deficient counts are red and below the line.

Above, we show the mosaic plot (on the left) and association plot (on the right). Both of these displays demonstrate that there is an association. The mosaic plot indicates that only about a quarter of the sample is female (indicated by the width of the columns), and that homelessness is present in about half the subjects (area shaded in light grey). The slight association demonstrated is that there are fewer homeless women than expected (since the horizontal line moves down between the first and second column). Similarly, for the association plot we note that the expected cell count is less than the observed (indicated in red with values below the line) for the female homeless group.


par(mfrow=c(1,2))
mosaicplot(tab, color=TRUE, main="mosaic plot")
assocplot(tab)
title("association plot")

SAS
As in Example 8.32, we find SAS macros for mosaic plots among the contributions of Michael Friendly. In this complex case, they are somewhat more difficult to access than others. The code for the plots themselves can be downloaded here, while it's useful to also run a wrapper macro. After downloading the files, the following code can be used to make the figure below.


title 'Install mosaic modules';
* location of the zipped files;
filename mosaic  'c:\ken\sasmacros\mosaics';
* storage location of compiled macros;
libname  mosaic   'c:\ken\sasmacros\mosaics';

* Code to read in, compile and store the macros;
proc iml ;
   reset storage=mosaic.mosaic;
   %include mosaic(mosaics) ;
   store module=_all_;
   show storage;
quit;

* Prep: create the table, save the cell counts;
proc freq data = "c:\book\help.sas7bdat";
tables homeless * female / out=outhelp;
run;

* Read in the wrapper macro;
%include "c:\ken\sasmacros\mosaics\mosaic.sas";

* Make the plot;
%mosaic(data=outhelp,var = female homeless, 
        sort=homeless descending female, space = 1 1);

The sort and space options make the results more similar to those shown for mosaicplot(). In this version, the colors reflect the signs of the residuals.

Example 9.2: Transparency and bivariate KDE

2011-07-11T13:57:00.012-04:00

In Example 9.1, we showed a binning approach to plotting bivariate relationships in a large data set. Here we show more sophisticated approaches: transparent overplotting and formal two-dimensional kernel density estimation. We use the 10,000 simulated bivariate normals shown in Example 9.1.

SAS
In SAS, transparency can be found in proc sgplot, with results shown above. The options here are fairly self-explanatory.


proc sgplot data=mvnorms;
  scatter x=x1 y=x2 / markerattrs=(symbol=CircleFilled size = .05in) 
          transparency=0.85;
run;

The image gives a good sense of the overall density, with the darker (overplotted) areas reflecting more observations. Overplotting was the problem we sought to avoid with the binning, but here it becomes an advantage.

Another approach is to use bivariate kernel density estimation. This is perhaps more similar to the binning shown previously, but without the stricture of regular polygons. It also offers some default values for smoothing, though whether or not these are good default values could be debated.


proc kde data=mvnorms;
  bivar x1 x2 / plots=contour;
run;

R

In R, the basic plot() function appears to include transparency, though you must select a suitably pale color to see it. The pch, col, and cex parameters govern the shape, color, and size of the plotted symbols, respectively.


plot(xvals[,1], xvals[,2], pch=19, col="#00000022", cex=0.1)

Bivariate kernel density estimation is available in the smoothScatter() function, which is in included in the R distribution as part of the graphics package.


smoothScatter(xvals[,1], xvals[,2])

Example 9.1: Scatterplots with binning for large datasets

2011-07-05T14:02:00.006-04:00

Scatterplots can get very hard to interpret when displaying large datasets, as points inevitably overplot and can't be individually discerned. A number of approaches have been crafted to help with this problem. One approach uses binning. This approach is also sometimes called a heat map, and can be though of as a two-dimensional histogram, where shades of the bins take the place of the heights of the bars. Any regular tesselation of the plane can be used, but there is some attraction to using hexagons. Why? In the vignettes for the hexbin package author Nicholas Lewin-Koh notes:

There are many reasons for using hexagons, at least over squares. Hexagons have symmetry of nearest neighbors which is lacking in square bins. Hexagons are the maximum number of sides a polygon can have for a regular tesselation of the plane, so in terms of packing a hexagon is 13% more efficient for covering the plane than squares. This property translates into better sampling efficiency at least for elliptical shapes. Lastly hexagons are visually less biased for displaying densities than other regular tesselations.

On the other hand, it's unclear whether these advantages are relevant here or whether they outweigh the simplicity of the square and the constant x and y values accompanying it.

In this entry, we demonstrate the use of a binned scatterplot for data from a sample of 10,000 generated bivariate normal random variables (section 1.10.6).

R

In R, we use the hexbin package to generate our plot, after generating our bivariate normals with correlation approximately 0.52.


library(MASS)
library(hexbin)
mu = c(1, -1)
Sigma = matrix(c(3, 2,
                 2, 5), nrow=2)
xvals = mvrnorm(10000, mu, Sigma)
Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])    # correlation
plot(hexbin(xvals[,1], xvals[,2]), xlab="X1", ylab="X2")

SAS
We're not aware of a SAS procedure to generate a binned scatterplot or of previously existing macros to do it. Ken wrote a relatively simple macro to do it, which can be found here. The macro uses proc gmap, and we hope that someone will develop an approach using proc template and proc sgrender, as demonstrated in an example from SAS Institute.

After running the macro, the following code generates the image shown below.


data Sigma (type=cov);
infile cards;
input _type_ $ _Name_ $ x1 x2;
cards;
cov x1 3 2
cov x2 2 5
;
run;

proc simnormal data=Sigma out=mvnorms numreal = 10000;
  var x1 x2;
run;

%twodhist(data=mvnorms,x=x1,y=x2,nbinsx=30,nbinsy=30,nshades=9);

We note that the default number of shades shown in R, and the number chosen here for SAS, seem to exceed the eye's ability to differentiate, especially for the darker shades.

Update

An anonymous commenter reported that the SAS code bombed when run. I (Ken) added a new version of the code at the link listed above. I note it here only to emphasize that in either SAS or R, settings or objects in the environment can affect the performance of code. If your plan to share code, an item to add to your checklist is to run the code in a fresh session.

A third year of entries!

2011-07-01T09:50:00.006-04:00

Contrary to previous reports, we started blogging after our book was published, with the conceit that we were adding examples to the book. Today marks the second anniversary of the book's appearance and of the blog. To celebrate, we're turning over our calendar, and starting a new volume of entries-- Example 9.1 will appear on July 5th.

It's very difficult to get an accurate measure of our viewership. As I write this, Feedburner reports about 650 regular readers, but this omits people who see us on R-bloggers and SAS Community Planet or SAS-X. As consumers of those aggregators, we assume there are many others who see us without subscribing directly.

We also get a fair number of views that derive directly from Google searches, which means we must be doing something right. Google Analytics reports over 100,000 total pageviews, while Feedburner claims 250,000.

So far, our most popular entries, according to feedburner are:
Example 7.38: Kaplan-Meier survival estimates
Example 7.39: Nelson-Aalen estimate of survivorship
Example 8.3: Pyramid plots
Example 7.30: Simulate censored survival data
Example 8.5: Bubble plots part 3

Blogger's internal metrics vary somewhat:
Example 7.35: Propensity score matching
Example 8.7: Hosmer and Lemeshow goodness-of-fit
Example 7.38: Kaplan-Meier survival estimates
Example 7.30: Simulate censored survival data
Example 7.25: Compare draws with distribution

All of your attention is really gratifying, and we hope we're useful to people-- it's hard work finding new material and collaborating on a new post every week!

Example 8.42: skewness and kurtosis and more moments (oh my!)

2011-06-27T10:45:00.000-04:00

While skewness and kurtosis are not as often calculated and reported as mean and standard deviation, they can be useful at times. Skewness is the 3rd moment around the mean, and characterizes whether the distribution is symmetric (skewness=0). Kurtosis is a function of the 4th central moment, and characterizes peakedness, where the normal distribution has a value of 3 and smaller values correspond to thinner tails (less peakedness).

Some packages (including SAS) subtract three from the kurtosis, so that the normal distribution has a kurtosis of 0 (this is sometimes called excess kurtosis.

R


library(moments)
library(lattice)
ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
ds$gender = ifelse(ds$female==1, "female", "male")
densityplot(~ cesd, data=ds, groups=gender, auto.key=TRUE)

We see that the distribution of CESD scores is skewed with a long left tail, and appears somewhat less peaked than a normal distribution. This is confirmed by the actual statistics:


> with(ds, tapply(cesd, gender, skewness))
    female       male 
-0.4906171 -0.2464390 
> with(ds, tapply(cesd, gender, kurtosis))    # kurtosis
  female     male 
2.748968 2.547061 
> with(ds, tapply(cesd, gender, kurtosis))-3  # excess kurtosis
    female       male 
-0.2510318 -0.4529394

SAS
SAS includes much detail on the moments and other statistics in the output from proc univariate. As usual, the quantity of output can be off-putting for new users and students. Here we extract the moments we need with the ODS system. We also generate kernel density estimates roughly analogous to the densityplot() results shown above.


ods output moments = cesdmoments;
proc univariate data="c:\book\help.sas7bdat";
  class female;
  var cesd;
  histogram cesd / kernel;
run;

proc print data=cesdmoments; 
  where label1 = "Skewness";
  var female label1 nvalue1 label2 nvalue2;
run;

With the result:


Obs   FEMALE    Label1         nValue1    Label2         nValue2

  4     0      Skewness      -0.247513   Kurtosis      -0.442010
 10     1      Skewness      -0.497620   Kurtosis      -0.204928

We note that the default is to produce unbiased (REML) estimates, rather than the biased method of moments estimator produced by the kurtosis() function (and that SAS presents the excess kurtosis).