R help - Mar 2011 - Pairwise T-Tests and Dunnett's Test (possibly using multcomp)

Hello Everyone,
 
I've been learning to use R in my spare time over the past several months.
I've read about 7-8 books on the subject. Lately I've been testing what
I've learned by trying to replicate the analyses from some of my SAS books.
This helps me make sure I know how to use R properly and also helps me to
understand how the two programs are similar and different.
 
Below is an attempt to replicate some SAS analyses. This was a success in that I
was able to reproduce the results using R. I have a feeling though that it may
not represent an ideal approach. So I was hoping a few people might be willing
to look at what I've done and to suggest improvements or alternatives.
 
One thing I'm struggling with is setting the reference level. I inserted a
comment about this in the R code.
 
I've also pasted in the original SAS code in case some people are also SAS
users and might find this helpful or interesting.
 
Thanks,
 
Paul
 
  
R Code:
 
##################################
#### Chapter 6: One-Way ANOVA ####
##################################
 
#### Example 6.1: One-Way ANOVA using HAM-A Scores in GAD #### 
 
connection <- textConnection("
101 LO 21  104 LO 18
106 LO 19  110 LO  .
112 LO 28  116 LO 22
120 LO 30  121 LO 27
124 LO 28  125 LO 19
130 LO 23  136 LO 22
137 LO 20  141 LO 19
143 LO 26  148 LO 35
152 LO  .  103 HI 16
105 HI 21  109 HI 31
111 HI 25  113 HI 23
119 HI 25  123 HI 18
127 HI 20  128 HI 18
131 HI 16  135 HI 24
138 HI 22  140 HI 21
142 HI 16  146 HI 33
150 HI 21  151 HI 17
102 PB 22  107 PB 26
108 PB 29  114 PB 19
115 PB  .  117 PB 33
118 PB 37  122 PB 25
126 PB 28  129 PB 26
132 PB  .  133 PB 31
134 PB 27  139 PB 30
144 PB 25  145 PB 22
147 PB 36  149 PB 32
")
 
gad <- data.frame(scan(connection,
   list(patno=0,dosegrp="",hama=0),na.strings="."))
 
gad
 
#### Summary Statistics ####
 
summaryfn <- function(x) data.frame(
   mean=mean(x,na.rm=T),std.dev=sd(x,na.rm=T),n=sum(!is.na(x)))
 
by(gad$hama,gad$dosegrp,summaryfn)
 
#### ANOVA ####
 
model1 <- aov(hama~dosegrp,data=gad)
summary.aov(model1)
 
#### Multiple Comparisons ####
 
library(multcomp)
 
#### Pairwise T-Tests ####
 
cht <- glht(model1, linfct = mcp(dosegrp = "Tukey"))
summary(cht, test = univariate())
 
#### Dunnett's Test ####
 
#Not sure how to set the reference group to placebo
#cht <- glht(model1, linfct = mcp(dosegrp = "Dunnett"))
#summary(cht, test = univariate())
 
model2 <- aov(hama ~ dosegrp - 1, data = gad)
K <- rbind(c(1, 0, -1), c(0, 1, -1))
rownames(K) <- c("HI vs PL", "LO vs PL")
colnames(K) <- names(coef(model2))
summary(glht(model2, linfct = K))
 
#### Treatment vs. Placebo ####
 
#SAS produces an F-test but this is equivalent
 
K <- rbind(c(0.5, 0.5, -1))
rownames(K) <- c("TR vs PL")
colnames(K) <- names(coef(model2))
summary(glht(model2, linfct = K))

SAS code:
 
*-----------------------------------------------------------------------------------*;
** Chapter 6: One-Way ANOVA;
 
/* Example 6.1: HAM-A Scores in GAD */
 
DATA GAD;
    INPUT PATNO DOSEGRP $ HAMA @@;
    DATALINES;
101 LO 21  104 LO 18
106 LO 19  110 LO  .
112 LO 28  116 LO 22
120 LO 30  121 LO 27
124 LO 28  125 LO 19
130 LO 23  136 LO 22
137 LO 20  141 LO 19
143 LO 26  148 LO 35
152 LO  .  103 HI 16
105 HI 21  109 HI 31
111 HI 25  113 HI 23
119 HI 25  123 HI 18
127 HI 20  128 HI 18
131 HI 16  135 HI 24
138 HI 22  140 HI 21
142 HI 16  146 HI 33
150 HI 21  151 HI 17
102 PB 22  107 PB 26
108 PB 29  114 PB 19
115 PB  .  117 PB 33
118 PB 37  122 PB 25
126 PB 28  129 PB 26
132 PB  .  133 PB 31
134 PB 27  139 PB 30
144 PB 25  145 PB 22
147 PB 36  149 PB 32
;
 
PROC SORT DATA = GAD; BY DOSEGRP;
PROC MEANS MEAN STD N DATA = GAD;
    BY DOSEGRP;
    VAR HAMA;
    TITLE1 'One-Way ANOVA';
    TITLE2 'EXAMPLE 6.1: HAM-A Scores in GAD';
RUN;
 
PROC GLM DATA = GAD;
    CLASS DOSEGRP;
    MODEL HAMA = DOSEGRP;
    MEANS DOSEGRP/T DUNNETT('PB');
    CONTRAST 'ACTIVE vs. PLACEBO' DOSEGRP 0.5 0.5 -1;
RUN;


	[[alternative HTML version deleted]]

Hello Everyone,
 
Figured out one part of the code. Setting the reference level for a factor is
accomplished using the relevel funtion (pg. 383 of MASS; pg. 70 of Data
Manipulation with R):
 
gad$dosegrp <- relevel(gad$dosegrp,3)
 
This works very well. Much better than using a format in SAS procedures that
don't allow the "ref=" option for instance.
 
Does anyone have suggestions about how to improve other aspects of my code?
 
Thanks,
 
Paul
 


	[[alternative HTML version deleted]]