Simultaneous Inferences in R

Using package:multcomp

mice example

The multcomp does not come with the standard R distribution. You need to download it from CRAN. See Install R Packages for example on the installation of R add-on packages.

> require(multcomp) # load package:multcomp

Loading required package: multcomp
Loading required package: mvtnorm
[1] TRUE

> mice <- data.frame(group=as.factor(rep(c("C","1","2"),rep(6,3))),

+ score=scan())

1: 58 32 59 64 55 49 73 70 68 71 60 62 53 74 72 62 58 61

19:

Read 18 items

> levels(mice$group) # show factor levels

[1] "1" "2" "C"

> # We wish to use "C" as the reference level, so do:

> mice$group <- relevel(mice$group, ref="C")

> levels(mice$group) # now "C" becomes the reference level

[1] "C" "1" "2"

> mice.aov <- aov(score ~ group, data=mice) # fit aov model

> summary(mice.aov) # produce the summary of the fit

            Df  Sum Sq Mean Sq F value  Pr(>F)  
group        2     673   336.5   4.549 0.0286 *
Residuals   15    1110    74.0
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

>

>

> # The following shows Tukey's all pairwise comparisons

> # using glht (general linear hypothesis test)

> mice.glht <- glht(mice.aov, linfct = mcp(group="Tukey"))

> summary(mice.glht) # the summary of the tests

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = score ~ group, data = mice)

Linear Hypotheses:
           Estimate Std. Error t value p value  
1 - C == 0   14.500      4.965   2.920  0.0268 *
2 - C == 0   10.500      4.965   2.115  0.1204  
2 - 1 == 0   -4.000      4.965  -0.806  0.7054  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported -- single-step method)

>

> # The results can be displayed graphically (function confint

> # was silently called to produce confidence intervals)

> windows(width=5,height=3,pointsize=10)

> plot(mice.glht)

> title(sub="Mice Data",adj=0)

> mtext("Tukey Honest Significant Differences",side=3,line=0.5)

(mice05.png here)

>

>

> # The following demonstrates Dunnett's method

> mice.Dunnett <- glht(mice.aov, linfct=mcp(group="Dunnett"))

> summary(mice.Dunnett)

       Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = score ~ group, data = mice)

Linear Hypotheses:
           Estimate Std. Error t value p value  
1 - C == 0   14.500      4.965   2.920  0.0195 *
2 - C == 0   10.500      4.965   2.115  0.0916 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported -- single-step method)

> windows(width=5,height=3,pointsize=10)

> plot(mice.Dunnett,sub="Mice Data")

> mtext("Dunnet's Method",side=3,line=0.5)

(mice06.png here)

> confint(mice.Dunnett) # show the confidence intervals

         Simultaneous Confidence Intervals

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = score ~ group, data = mice)

Quantile = 2.4392
95% family-wise confidence level


Linear Hypotheses:
           Estimate lwr     upr    
1 - C == 0 14.5000   2.3884 26.6116
2 - C == 0 10.5000  -1.6116 22.6116

> confint(mice.Dunnett,level=0.90) # can change confidence level

         Simultaneous Confidence Intervals

Multiple Comparisons of Means: Dunnett Contrasts


Fit: aov(formula = score ~ group, data = mice)

Quantile = 2.0663
90% family-wise confidence level


Linear Hypotheses:
           Estimate lwr     upr    
1 - C == 0 14.5000   4.2399 24.7601
2 - C == 0 10.5000   0.2399 20.7601

>

> # The plot below displays 90\% Dunnett's intervals

> plot(confint(mice.Dunnett,level=0.90),sub="Mice Data",

+ ylim=c(0.5,2.5))

> mtext("Dunnett's Method",side=3,line=0.5)

(mice07.png here)

>

>

> # Can also do general contrasts

> contr <- rbind("avg1&2 - C"=c(-1,.5,.5),"2 - 1"=c(0,-1,1))

> confint(glht(mice.aov,linfct=mcp(group=contr)))

       Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = score ~ group, data = mice)

Estimated Quantile = 2.4742

Linear Hypotheses:
                Estimate lwr      upr     
avg1&2 - C == 0  12.5000   1.8605  23.1395
2 - 1 == 0       -4.0000 -16.2854   8.2854

95% family-wise confidence level