6

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?

Community
  • 1
Nate
  • 768
  • 5
  • 10
  • 3
    You've made everything reproducible, which is excellent. But you're expecting people to run this to understand what is going on. Some will, but you'll get more readers if you actually show some of the results too. – Nick Cox Jul 10 '13 at 16:07
  • 5
    So, what specific estimates don't match? Your factors all have more than 2 levels, so the ANOVA table is going to show effects that take up more than 1 degree of freedom (pooled effects, in a way), whereas all parameter estimates shown in a regression output are going to be 1-df effects. In other words, in order to get the same F test and p-value in regression as in ANOVA for the effect of a factor with more than 2 levels, you would need to do a nested-model comparison between a regression model including all terms, and one excluding all dummy codes affiliated with the factor in question... – Patrick Coulombe Jul 10 '13 at 16:15
  • (cont.) ...One way to check that you've estimated the same models would be to compare the R² obtained in regression to the sum of your SSeffects divided by SStotal in ANOVA – Patrick Coulombe Jul 10 '13 at 16:16
  • 1
    In the anova() nutrient was sig; in the summary() one level of it was sig. Similarly for herbivore, and more or less the opposite for the interaction. anova() and summary() ask different questions, so they give different answers (at least sometimes). – Peter Flom Jul 10 '13 at 17:05
  • Someone was nice enough to add the output for me, just as I was trying to figure out how to do it. I have no idea how they did that, though. I was trying to use Markdown, but apparently using it incorrectly. – Nate Jul 10 '13 at 17:32
  • I understand Patrick's comments, but what I don't understand is how a series of significant interactions becomes non-significant when pooled together (albeit with non-significant interactions). What I understand even less is how a series of non-significant coefficients becomes significant when pooled. – Nate Jul 10 '13 at 17:35
  • 1
    @Nate you can use the brackets that say "code sample" – Peter Flom Jul 10 '13 at 17:46
  • Nate, when you're editing, investigate the functions of the tools across the top left, but also when you edit there's a little gold circle at the top right with a "?" in it. It's worth clicking. But if you're familiar with Markdown, the way to get code blocks formatted (including output) is the same as in vanilla Markdown - indent 4 spaces. The "**{}**" tool does that for blocks of code. – Glen_b Jul 11 '13 at 05:21

2 Answers2

4

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.

January
  • 6,999
  • 1
  • 32
  • 55
2

What is going on here is a multiple comparisons issue. You have 10 df for interaction so you can look at 10 independent interaction effects (although I suspect the 10 you are actually looking at, i.e. the 10 regression effects, are not independent).

The 10 df interaction test will be significant at some level if and only if Scheffe's multiple comparison test can find an interaction effect that is significant at that level. So using Scheffe's method you would not be able to find an interaction regression coefficient that is significant. What is being report as P values for the regression coefficients are equivalent to looking at Fisher's LSD multiple comparison methods, which is notoriously easier at declaring significance. So basically you have one method that is declaring no effects and another that finds a few but since they are different methods that is not surprising. You need to decide what standards you want to use. (A more sophisticated use of LSD would not look at individual coefficients unless the overall test was significant.)

Another way of thinking about this is that the 10df interaction test is an average of ten 1df tests and if the interaction effects are not very striking they can get lost in the process of averaging them. However, if you look at them individually, you can see their effect.

I will not get into the main effects issues. But I think what R is telling you most strongly about the interactions (P=.00513) is that the differential effect of using Nutrients a and d changes depending on whether you use the unnamed Herbivore or the Paired Herbivore. If the differential effect of a and d can change, then there has to be some effect for the pair a and d, however the regression coefficient for Nut. d (which really looks at their difference) SEEMS to be saying there is none -- but it only seems to be saying that because main effects in the presence of interaction get so convoluted that they are not worth trying to figure out.

Ronald Christensen
  • 21
  • 1