[R] Planned comparison ANOVA

Greg Snow Greg.Snow at imail.org
Mon Sep 19 18:37:46 CEST 2011 

Here is one example of doing it, first create a matrix with all your contrasts (I added some to make it full), then invert it to get the dummy variable encodings, use the 'contrast' function (or 'C') to set the contrasts, then use 'aov' or 'lm' and those contrasts will be used:

> my.contrasts <- rbind( 1/5, '1 v 5' = c(1,0,0,0,-1), 
+ '12 v 45' = c(1,1,0,-1,-1)/2, '1 v 2' = c(1,-1,0,0,0),
+ '3 v 1245' = c(-1,-1,4,-1,-1)/4 )
> my.contrasts
          [,1]  [,2] [,3]  [,4]  [,5]
          0.20  0.20  0.2  0.20  0.20
1 v 5     1.00  0.00  0.0  0.00 -1.00
12 v 45   0.50  0.50  0.0 -0.50 -0.50
1 v 2     1.00 -1.00  0.0  0.00  0.00
3 v 1245 -0.25 -0.25  1.0 -0.25 -0.25
> 
> my.dv <- solve(my.contrasts)
> my.dv
       1 v 5 12 v 45 1 v 2 3 v 1245
[1,] 1     0     0.5   0.5     -0.2
[2,] 1     0     0.5  -0.5     -0.2
[3,] 1     0     0.0   0.0      0.8
[4,] 1     1    -1.5  -0.5     -0.2
[5,] 1    -1     0.5   0.5     -0.2
> 
> g <- factor(rep(1:5, each=20))
> y <- rnorm(100, 10, 2)
> 
> contrasts(g) <- my.dv[,-1]
> 
> fit <- aov( y ~ g )
> summary(fit, split=list(g=1:4))
            Df Sum Sq Mean Sq F value  Pr(>F)  
g            4  19.20  4.8007  1.3110 0.27151  
  g: C1      1   2.01  2.0054  0.5476 0.46111  
  g: C2      1  15.95 15.9470  4.3548 0.03958 *
  g: C3      1   0.28  0.2814  0.0768 0.78223  
  g: C4      1   0.97  0.9689  0.2646 0.60819  
Residuals   95 347.88  3.6619                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> summary.lm(fit)

Call:
aov(formula = y ~ g)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.4519 -1.1546  0.0749  1.2288  5.2639 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.7437     0.1914  50.918   <2e-16 ***
g1 v 5       -1.0330     0.6051  -1.707   0.0911 .  
g12 v 45     -0.8929     0.4279  -2.087   0.0396 *  
g1 v 2        0.1677     0.6051   0.277   0.7822    
g3 v 1245     0.2461     0.4784   0.514   0.6082    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 1.914 on 95 degrees of freedom
Multiple R-squared: 0.05231,    Adjusted R-squared: 0.01241 
F-statistic: 1.311 on 4 and 95 DF,  p-value: 0.2715 

> 
> mean( y[g==1] ) - mean(y[g==5]) # compare to estimate for 1st contrast
[1] -1.032981
> (mean( y[g==1] )+mean(y[g==2]) - mean(y[g==4]) - mean(y[g==5]) )/2
[1] -0.8929435
> mean( y[g==1] ) - mean(y[g==2])
[1] 0.1677452
> mean( y[g==3] ) - mean(y[g!=3]) # change this if unbalanced
[1] 0.2460775
>

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at imail.org
801.408.8111


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-
> project.org] On Behalf Of robcinm
> Sent: Saturday, September 17, 2011 10:59 PM
> To: r-help at r-project.org
> Subject: [R] Planned comparison ANOVA
> 
> I am trying to do a priori ANOVA analysis for a class assignment. The
> professor uses SPSS and does not know R. I want to do a simple planned
> comparison but have been unable to find a function or specific help.
> There
> is a grouping variable with five levels and a subsequent response
> variable.
> I also created a few columns that contain my group contrasts to see if
> anything could come of that. My first hypothesis, for example, is that
> group
> 1 and group 5 differ. My second is that groups 1 & 2 differ from 4 & 5
> and
> so on...
> 
> I searched and found people wanting similar help, but their projects
> were a
> little beyond my R comprehension right now and the help given did not
> apply.
> 
> --
> View this message in context: http://r.789695.n4.nabble.com/Planned-
> comparison-ANOVA-tp3821357p3821357.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.

More information about the R-help mailing list