Refugee from SPSS having issues with fit.contrast in R

Question

I'm trying to wean myself off of SPSS and switch to R for psych statistics and I'm mostly getting there, but I've run into a real roadblock with fit.contrast from the gmodels package. As far as I can tell, even the provided example doesn't work properly--either that, or I grossly misunderstand contrasts in R (always a possibility).

First let's look at the example code that comes with the documentation for gmodels::fit.contrast():

set.seed(03215)
Genotype <- sample(c("WT","KO"), 1000, replace=TRUE)
Time <- factor(sample(1:3, 1000, replace=TRUE))
y <- rnorm(1000)
data <- data.frame(y, Genotype, Time)

model <- aov( y ~ Genotype + Time + Genotype:Time, data=data )
model.tables(model, "means")

fit.contrast( model, "Genotype", rbind("KO vs WT"=c(-1,1) ), conf=0.95 , df=T)

For me, this produces the following output:

                   Estimate Std. Error   t value Pr(>|t|)   lower CI  upper CI
GenotypeKO vs WT 0.01683876  0.1095764 0.1536714   0.8779 -0.1981888 0.2318664

Now consider the later example code:

model <- aov( y ~ Genotype + Time + Genotype:Time, data=data,
              contrasts=list(Genotype=make.contrasts(cm.G),
                             Time=make.contrasts(cm.T) )
)

summary(model, split=list( Genotype=list( "KO vs WT"=1 ),
                           Time = list( "1 vs 2" = 1,
                                        "2 vs 3" = 2 ) ) )

This produces this output for me:

                                  Df Sum Sq Mean Sq F value Pr(>F)
Genotype                           1    1.2  1.1687   1.121  0.290
  Genotype: KO vs WT               1    1.2  1.1687   1.121  0.290
Time                               2    0.7  0.3677   0.353  0.703
  Time: 1 vs 2                     1    0.2  0.1784   0.171  0.679
  Time: 2 vs 3                     1    0.6  0.5571   0.534  0.465
Genotype:Time                      2    1.2  0.5760   0.552  0.576
  Genotype:Time: KO vs WT.1 vs 2   1    0.3  0.3265   0.313  0.576
  Genotype:Time: KO vs WT.2 vs 3   1    0.8  0.8256   0.792  0.374
Residuals                        994 1036.6  1.0429  

I'm completely flummoxed here. The p-values for the two tests on Genotype should match, but they don't! The p-values are not the same and the F-value isn't the square of the t-value. Also, the value computed for the contrast estimate (0.0168) is NOT the difference in the two means... why not?

Any help in understanding what fit.contrast is doing here would be greatly appreciated!

Answer 1

My initial remark is: why do you need these functions from the gmodels package to accomplish this? These are basic tasks that can be accomplished straightforwardly in base R, which is what I would recommend.

But to answer your question directly, the issue here is that fit.contrast is using "type 3" tests of effects while summary is using "type 2" tests. This can be verified most easily using Anova from the car package, which lets you select the type of sums of squares that you desire:

library(car)

> Anova(model, type=2)
Anova Table (Type II tests)

Response: y
Sum Sq  Df F value Pr(>F)
Genotype         1.09   1  1.0495 0.3059
Time             0.74   2  0.3526 0.7029
Genotype:Time    1.15   2  0.5524 0.5758
Residuals     1036.64 994       

> Anova(model, type=3)
Anova Table (Type III tests)

Response: y
Sum Sq  Df F value Pr(>F)
(Intercept)      0.10   1  0.0935 0.7599
Genotype         0.02   1  0.0236 0.8779
Time             1.87   2  0.8954 0.4088
Genotype:Time    1.15   2  0.5524 0.5758
Residuals     1036.64 994

A description of the difference between these two methods can be found here: https://stat.ethz.ch/pipermail/r-help/2006-August/111854.html

Answer 2

Strongly recommended to not use reserved words as R objects names. For genotype contrast it is clear that the p-values from ANOVA and contrast statement from gmodels are identical and F=t^2.

library(gmodels)
set.seed(03215)
Genotype <- sample(c("WT","KO"), 1000, replace=TRUE)
Time <- factor(sample(1:3, 1000, replace=TRUE))
y <- rnorm(1000)
dat <- data.frame(y, Genotype, Time)

fit1 <- aov( y ~ Genotype + Time + Genotype:Time, data=dat)
summary(fit1)


              Df Sum Sq Mean Sq F value Pr(>F)
Genotype        1    1.2  1.1687   1.121  0.290
Time            2    0.7  0.3677   0.353  0.703
Genotype:Time   2    1.2  0.5760   0.552  0.576
Residuals     994 1036.6  1.0429

model.tables(fit1, "means")
Tables of means
Grand mean

0.01447773 

 Genotype 
           KO        WT
     -0.02154   0.04693
rep 474.00000 526.00000

 Time 
            1         2         3
      0.03267   0.03313  -0.02539
rep 350.00000 334.00000 316.00000

 Genotype:Time 
        Time
Genotype 1      2      3     
     KO    0.02   0.02  -0.11
     rep 160.00 155.00 159.00
     WT    0.04   0.04   0.06
     rep 190.00 179.00 157.00

As (-1)(-0.02154)+(1)(0.04693) = 0.06847

fit.contrast(fit1, "Genotype", rbind("KO vs WT"=c(1, -1)), conf=0.95, df=TRUE)

                  Estimate Std. Error  t value  Pr(>|t|)  DF   lower CI upper CI
GenotypeKO vs WT 0.06869178 0.06477589 1.060453 0.2891962 994 -0.0584214 0.195805