Solved – ANOVA and Regression give opposite results in R

anovainteractionrregression

I conducted an experiment in a factorial design: I measured light (PAR) in three herbivore treatments as well as six nutrient treatments. The experiment was blocked.

I've run the linear model as follows (you can download the data from my website to replicate)

dat <- read.csv('http://www.natelemoine.com/testDat.csv')
mod1 <- lm(light ~ Nutrient*Herbivore + BlockID, dat)

The residual plots look pretty good

par(mfrow=c(2,2))
plot(mod1)

When I look at the ANOVA table, I see main effects of Nutrient and Herbivore. 

anova(mod1)

Analysis of Variance Table 

Response: light 
                    Df  Sum Sq Mean Sq F value    Pr(>F)     
Nutrient             5  4.5603 0.91206  7.1198 5.152e-06 *** 
Herbivore            2  2.1358 1.06791  8.3364 0.0003661 *** 
BlockID              9  5.6186 0.62429  4.8734 9.663e-06 *** 
Nutrient:Herbivore  10  1.7372 0.17372  1.3561 0.2058882     
Residuals          153 19.5996 0.12810                       
--- 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

However, the regression table shows non-significant main effects and significant interactions.

summary(mod1)

Call: 
lm(formula = light ~ Nutrient * Herbivore + BlockID, data = dat) 

Residuals: 
     Min       1Q   Median       3Q      Max  
-0.96084 -0.19573  0.01328  0.24176  0.74200  

Coefficients: 
                           Estimate Std. Error t value Pr(>|t|)     
(Intercept)                1.351669   0.138619   9.751  < 2e-16 *** 
Nutrientb                  0.170548   0.160064   1.066  0.28833     
Nutrientc                 -0.002172   0.160064  -0.014  0.98919     
Nutrientd                 -0.163537   0.160064  -1.022  0.30854     
Nutriente                 -0.392894   0.160064  -2.455  0.01522 *   
Nutrientf                  0.137610   0.160064   0.860  0.39129     
HerbivorePaired           -0.074901   0.160064  -0.468  0.64049     
HerbivoreZebra            -0.036931   0.160064  -0.231  0.81784     
... 
Nutrientb:HerbivorePaired  0.040539   0.226364   0.179  0.85811     
Nutrientc:HerbivorePaired  0.323127   0.226364   1.427  0.15548     
Nutrientd:HerbivorePaired  0.642734   0.226364   2.839  0.00513 **  
Nutriente:HerbivorePaired  0.454013   0.226364   2.006  0.04665 *   
Nutrientf:HerbivorePaired  0.384195   0.226364   1.697  0.09168 .   
Nutrientb:HerbivoreZebra   0.064540   0.226364   0.285  0.77594     
Nutrientc:HerbivoreZebra   0.279311   0.226364   1.234  0.21913     
Nutrientd:HerbivoreZebra   0.536160   0.226364   2.369  0.01911 *   
Nutriente:HerbivoreZebra   0.394504   0.226364   1.743  0.08338 .   
Nutrientf:HerbivoreZebra   0.324598   0.226364   1.434  0.15362     
--- 
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.3579 on 153 degrees of freedom 
Multiple R-squared:  0.4176,    Adjusted R-squared:  0.3186  
F-statistic: 4.219 on 26 and 153 DF,  p-value: 8.643e-09

I know that this question has been previously asked and answered in multiple posts. In the earlier posts, the issue revolved around the different types of SS used in anova() and lm(). However, I don't think that is the issue here. First of all, the design is balanced:

with(dat, tapply(light, list(Nutrient, Herbivore), length))

Second, using the Anova() option doesn't change the anova table. This isn't a surprise because the design is balanced.

Anova(mod1, type=2)
Anova(mod1, type=3)

Changing the contrast doesn't change the results (qualitatively). I still get pretty much backwards intepretations from anova() vs. summary().

options(contrasts=c("contr.sum","contr.poly"))
mod2 <- lm(light ~ Nutrient*Herbivore + BlockID, dat)
anova(mod2)
summary(mod2)

I'm confused because everything I've read on regression not agreeing with ANOVA implicates differences in the way R uses SS for summary() and anova() functions. However, in the balanced design, the SS types are equivalent, and the results here don't change. How can I have completely opposite interpretations depending on which output I use?

Best Answer

Essentially, the question is, how come that one coefficient in the linear model is significantly different from 0, but ANOVA shows no significant effect and vice versa.

For this, let's consider a simpler example.

set.seed( 123 )
data <- data.frame( x= rnorm( 100 ), g= rep( letters[1:10], each= 10 ) )
data$x[ data$g == "d" ] <- data$x[ data$g == "d" ] + 0.5
boxplot( x ~ g, data )
l <- lm( x ~ 0 + g, data )
summary( l )
anova( l )

You can see that there is only one group (d) that stands out of the line (has a coefficient significantly different from zero). However, given that the nine other groups do not show an effect, the anova returns $p > 0.1$. However, let us remove some of the groups:

data2 <- data[ data$g %in% c( "a", "d" ), ]
anova( lm( x ~ 0 + g, data2 )

returns

          Df  Sum Sq Mean Sq F value  Pr(>F)  
g          2  6.8133  3.4066  5.7363 0.01182 *
Residuals 18 10.6898  0.5939

ANOVA considers the overall variance within and between the groups. In the first case (10 groups) the variance between the groups is smaller because of the many groups with no effect. In the second, there are only two groups, and all the between groups variance comes from the difference between these two groups.

How about the reverse? This is easier: imagine three groups with means equal to -1, 0, 1. Total average is 0. Each group separately does not necessarily has a significant difference from 0, but there is enough difference between group 1 and 3 to account for significant total between group variance.

Related Question