26

I am trying to move from using the ez package to lme for repeated measures ANOVA (as I hope I will be able to use custom contrasts on with lme).

Following the advice from this blog post I was able to set up the same model using both aov (as does ez, when requested) and lme. However, whereas in the example given in that post the F-values do perfectly agree between aov and lme (I checked it, and they do), this is not the case for my data. Although the F-values are similar, they are not the same.

aov returns a f-value of 1.3399, lme returns 1.36264. I am willing to accept the aov result as the "correct" one as this is also what SPSS returns (and this is what counts for my field/supervisor).

Questions:

  1. It would be great if someone could explain why this difference exists and how I can use lme to provide credible results. (I would also be willing to use lmer instead of lme for this type of stuff, if it gives the "correct" result. However, I haven't used it so far.)

  2. After solving this problem I would like to run a contrast analysis. Especially I would be interested in the contrast of pooling the first two levels of factor (i.e., c("MP", "MT")) and compare this with the third level of factor (i.e., "AC"). Furthermore, testing the third versus the fourth level of factor (i.e., "AC" versus "DA").

Data:

tau.base <- structure(list(id = structure(c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 
9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 
22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 
14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 20L, 21L, 22L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L), .Label = c("A18K", 
"D21C", "F25E", "G25D", "H05M", "H07A", "H08H", "H25C", "H28E", 
"H30D", "J10G", "J22J", "K20U", "M09M", "P20E", "P26G", "P28G", 
"R03C", "U21S", "W08A", "W15V", "W18R"), class = "factor"), factor = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("MP", "MT", "AC", "DA"
), class = "factor"), value = c(0.9648092876, 0.2128662077, 1, 
0.0607615485, 0.9912814024, 3.22e-08, 0.8073856412, 0.1465590332, 
0.9981672618, 1, 1, 1, 0.9794401938, 0.6102546108, 0.428651501, 
1, 0.1710644881, 1, 0.7639763913, 1, 0.5298989196, 1, 1, 0.7162733447, 
0.7871177434, 1, 1, 1, 0.8560509327, 0.3096989662, 1, 8.51e-08, 
0.3278862311, 0.0953598576, 1, 1.38e-08, 1.07e-08, 0.545290432, 
0.1305621416, 2.61e-08, 1, 0.9834051136, 0.8044114935, 0.7938839461, 
0.9910112678, 2.58e-08, 0.5762677121, 0.4750002288, 1e-08, 0.8584252623, 
1, 1, 0.6020385797, 8.51e-08, 0.7964935271, 0.2238374288, 0.263377904, 
1, 1.07e-08, 0.3160751898, 5.8e-08, 0.3460325565, 0.6842217296, 
1.01e-08, 0.9438301877, 0.5578367224, 2.18e-08, 1, 0.9161424562, 
0.2924856039, 1e-08, 0.8672987992, 0.9266688748, 0.8356425464, 
0.9988463913, 0.2960361777, 0.0285680426, 0.0969063841, 0.6947998266, 
0.0138254805, 1, 0.3494775301, 1, 2.61e-08, 1.52e-08, 0.5393467752, 
1, 0.9069223275)), .Names = c("id", "factor", "value"), class = "data.frame", row.names = c(1L, 
6L, 10L, 13L, 16L, 17L, 18L, 22L, 23L, 24L, 27L, 29L, 31L, 33L, 
42L, 43L, 44L, 45L, 54L, 56L, 58L, 61L, 64L, 69L, 73L, 76L, 79L, 
80L, 81L, 85L, 86L, 87L, 90L, 92L, 94L, 96L, 105L, 106L, 107L, 
108L, 117L, 119L, 121L, 124L, 127L, 132L, 136L, 139L, 142L, 143L, 
144L, 148L, 149L, 150L, 153L, 155L, 157L, 159L, 168L, 169L, 170L, 
171L, 180L, 182L, 184L, 187L, 190L, 195L, 199L, 202L, 205L, 206L, 
207L, 211L, 212L, 213L, 216L, 218L, 220L, 222L, 231L, 232L, 233L, 
234L, 243L, 245L, 247L, 250L))

And the code:

require(nlme)

summary(aov(value ~ factor+Error(id/factor), data = tau.base))

anova(lme(value ~ factor, data = tau.base, random = ~1|id))
amoeba
  • 104,745
Henrik
  • 14,198
  • 11
  • 69
  • 130
  • It looks like you just answered the part about the contrasts yourself in your answer here; if not, please edit this question so we know what difficulty remains. – Aaron left Stack Overflow Aug 11 '11 at 12:47
  • 2
    @Aaron, as long as there are differneces in the lme results from standard textbook ANOVA (given by aov, and which is what I need), this is not an option for me. In my paper I want to report an ANOVA, not something like an ANOVA. Interestingly Venables & Ripley (2002, p. 285) show that both approaches lead to identical estimates. But the differences in F values leave me with a bad feeling. Furthermore, Anova() (from car) returns only Chi²-values for lme objects. Therefore for me, my first question is not answered yet. – Henrik Aug 11 '11 at 13:00
  • I understand (but don't share) your wariness of lme; but for contrasts, glht works on lm fits too, not just lme fits. (Also, the lme results are standard textbook results too.) – Aaron left Stack Overflow Aug 11 '11 at 22:04
  • Unfortunately you cannot specify lm for a repeated measure analysis. Only aov can deal with repeated measures but will return an object of class aovlist which is unfortunately not handled by glht. – Henrik Aug 12 '11 at 08:56
  • Is this just a smaller example then? Because you don't need aov here as there are only within subject effects. – Aaron left Stack Overflow Aug 12 '11 at 14:09
  • @Aaron can you elaborate more on this. Why no aov? In general I am interested in a general solution with both between and within factors. But if there is another neat way for parts of this problem, I would be interested, too. – Henrik Aug 12 '11 at 14:29
  • 3
    lm uses the residual error as the error term for all effects; when there are effects that should use a different error term, aov is necessary (or instead, using the results from lm to compute the F-stats manually). In your example, the error term for factor is the id:factor interaction, which is the residual error term in an additive model. Compare your results to anova(lm(value~factor+id)). – Aaron left Stack Overflow Aug 12 '11 at 14:33
  • I've answered your first question below. As the answer is quite long, I would suggest starting a new question for your second question. – Aaron left Stack Overflow Aug 12 '11 at 16:00
  • Why do you write " random = ~1|id" and not " random = ~1|(id/factor)"? – skan Feb 22 '16 at 20:59

2 Answers2

32

They are different because the lme model is forcing the variance component of id to be greater than zero. Looking at the raw anova table for all terms, we see that the mean squared error for id is less than that for the residuals. 

> anova(lm1 <- lm(value~ factor+id, data=tau.base))

          Df  Sum Sq Mean Sq F value Pr(>F)
factor     3  0.6484 0.21614  1.3399 0.2694
id        21  3.1609 0.15052  0.9331 0.5526
Residuals 63 10.1628 0.16131

When we compute the variance components, this means that the variance due to id will be negative. My memory of expected mean squares memory is shaky, but the calculation is something like

(0.15052-0.16131)/3 = -0.003597.

This sounds odd but can happen. What it means is that the averages for each id are closer to each other than you would expect to each other given the amount of residual variation in the model.

In contrast, using lme forces this variance to be greater than zero.

> summary(lme1 <- lme(value ~ factor, data = tau.base, random = ~1|id))
...
Random effects:
 Formula: ~1 | id
        (Intercept)  Residual
StdDev: 3.09076e-05 0.3982667

This reports standard deviations, squaring to get the variance yields 9.553e-10 for the id variance and 0.1586164 for the residual variance.

Now, you should know that using aov for repeated measures is only appropriate if you believe that the correlation between all pairs of repeated measures is identical; this is called compound symmetry. (Technically, sphericity is required but this is sufficient for now.) One reason to use lme over aov is that it can handle different kinds of correlation structures.

In this particular data set, the estimate for this correlation is negative; this helps explain how the mean squared error for id was less than the residual squared error. A negative correlation means that if an individual's first measurement was below average, on average, their second would be above average, making the total averages for the individuals less variable than we would expect if there was a zero correlation or a positive correlation.

Using lme with a random effect is equivalent to fitting a compound symmetry model where that correlation is forced to be non-negative; we can fit a model where the correlation is allowed to be negative using gls:

> anova(gls1 <- gls(value ~ factor, correlation=corCompSymm(form=~1|id),
                    data=tau.base))
Denom. DF: 84 
            numDF   F-value p-value
(Intercept)     1 199.55223  <.0001
factor          3   1.33985   0.267

This ANOVA table agrees with the table from the aov fit and from the lm fit.

OK, so what? Well, if you believe that the variance from id and the correlation between observations should be non-negative, the lme fit is actually more appropriate than the fit using aov or lm as its estimate of the residual variance is slightly better. However, if you believe the correlation between observations could be negative, aov or lm or gls is better.

You may also be interested in exploring the correlation structure further; to look at a general correlation structure, you'd do something like

gls2 <- gls(value ~ factor, correlation=corSymm(form=~unclass(factor)|id),
data=tau.base)

Here I only limit the output to the correlation structure. The values 1 to 4 represent the four levels of factor; we see that factor 1 and factor 4 have a fairly strong negative correlation:

> summary(gls2)
...
Correlation Structure: General
 Formula: ~unclass(factor) | id 
 Parameter estimate(s):
 Correlation: 
  1      2      3     
2  0.049              
3 -0.127  0.208       
4 -0.400  0.146 -0.024

One way to choose between these models is with a likelihood ratio test; this shows that the random effects model and the general correlation structure model aren't statistically significantly different; when that happens the simpler model is usually preferred.

> anova(lme1, gls2)
     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
lme1     1  6 108.0794 122.6643 -48.03972                        
gls2     2 11 111.9787 138.7177 -44.98936 1 vs 2 6.100725  0.2965
amoeba
  • 104,745
Aaron left Stack Overflow
  • 7,198
  • Great answer. Thanks. One quick follow up. Is there a way to allow negative variances in lme? – Henrik Aug 14 '11 at 14:37
  • Not using lme, no. But gls with compound symmetry is equivalent, as I showed above. – Aaron left Stack Overflow Aug 14 '11 at 17:32
  • 2
    It is actually possible to use compound symmetry with lme to obtain the same results as with aov (and thereby enabling lme for all ANOVAs), namely using the correlation argument in the call to lme: anova(lme(value ~ factor, data = tau.base, random = ~1|id, correlation = corCompSymm(form = ~1|id))) – Henrik Aug 31 '11 at 15:57
  • 2
    Nice find. But isn't there an extra parameter in that fit? It has three variance parameters; variance for id, residual variance, and the correlation, while the gls has only a residual variance and a correlation. – Aaron left Stack Overflow Aug 31 '11 at 18:51
  • 1
    Your argument sounds plausible, however, the results do not agree. All the anova tables (aov, lme without compound symmetry, and lme with compound symmetry) have exactly the same number of dfs. – Henrik Aug 31 '11 at 19:49
  • 1
    You'll have to convince me that those three parameters are really an overparametrization of the first two. Have you worked out how they're related? – Aaron left Stack Overflow Sep 01 '11 at 03:55
  • 1
    No. I am trusting the output of anova.lme(). From your answer I got that the relation between the ANOVA and the mixed models lies in their correlation structure. I then read that imposing a compund symmetric correlation structure leads to equality between the two approaches. Therefore I imposed it. I have no idea if this eats up another df. The output however disagrees with this interpretation. – Henrik Sep 01 '11 at 07:15
  • 1
    Well, you're more trusting than me then. I know Doug Bates wasn't happy with the df calculation (hence the changes in lme4) so I'd be wary unless I knew exactly what was going on. – Aaron left Stack Overflow Sep 01 '11 at 14:20
  • Why "random = ~1|id" and not "random = ~ factor|id"? – skan Feb 22 '16 at 21:01
  • 1
    From memory, that would make different variances for each level of the factor, and I think different covariances for each pair as well. That could be worth exploring too, depending on the particulars of your situation. One potential downside is that it could be a lot of parameters to fit, especially if there are lots of levels of the factor. – Aaron left Stack Overflow Feb 23 '16 at 02:30
  • What's the meaning of a negative variance? As far as I know the variance is the sum of squares. – skan Jun 12 '17 at 12:09
  • Good question! It's more than can be answered in a comment, though. This answer https://stats.stackexchange.com/a/253323/3601 suggests that there are lots of questions about it, but in my searching I didn't find one that was either posed well or with a good clear canonical answer. If you're interested in benefitting the site, it would be great to do a search to confirm that and then write a new clear question about it, referencing the others, and hope that yours becomes the canonical question about it. – Aaron left Stack Overflow Jun 15 '17 at 16:12
5

aov() fits the model via lm() using least squares, lme fits via maximum likelihood. That difference in how the parameters of the linear model are estimated likely accounts for the (very small) difference in your f-values.

In practice, (e.g. for hypothesis testing) these estimates are the same, so I don't see how one could be considered 'more credible' than the other. They come from different model fitting paradigms.

For contrasts, you need to set up a contrast matrix for your factors. Venebles and Ripley show how to do this on p 143, p.146 and p.293-294 of the 4th edition.

Chris
  • 480
  • Hmm, but why then are there sometimes differences and sometimes the results are exactly equal? Furthemrore, then it seems to be impossible to use lme or lmerfor calculating an ANOVA (strictly speaking) as it uses a method that is similar but not identical. So is there no way of calculating contrasts for repeated measure ANOVAs in R? – Henrik Aug 10 '11 at 16:58
  • If the system your modelling is truly linear than least squares and ML should give the same f statistic. Its only when there is other structure in the data that the two methods will give different results. Pinheiro and Bates cover this in their mixed-effects models book. Also they're likely not 'exactly' equal, if you were to go far enough in sig digits I'm sure you'd find some difference. But for all practical purposes they are the same. – Chris Aug 10 '11 at 17:22