lme() vs aov()

VarCorr(mod.lme)

Hi Federico,

Note also that the intervals function is a good tool for checking
reliability of lme objects:

mod.lme = lme(response ~ drug, random = ~1|source, dat)
intervals(mod.lme)

Error in intervals.lme(mod.lme) :
  Cannot get confidence intervals on var-cov components: Non-positive
definite approximate variance-covariance


Kingsford





On Thu, May 22, 2008 at 1:19 PM, Peter Dalgaard
<p.dalgaard at biostat.ku.dk> wrote:

Federico Calboli wrote:

Hi All,

I was playing with a small dataset of 12 observations, a very basic nested
model, with 3 drugs, 2 sources for each drug and two response counts for
each source (the response is some medical parameter, of no real interest
here). The data is:

drug source response
d1 a 102
d1 a 104
d1 q 103
d1 q 104
d2 d 108
d2 d 110
d2 b 109
d2 b 108
d3 l 104
d3 l 106
d3 s 105
d3 s 107

For kicks, and because the data is balanced I thought that I could use it
to compare the results of aov() with those of lme() -- I know the library
lme4  and lmer() should be preferred, but the stuff I am ultimately testig
was done with lme.

In any case I fit 2 models and got 2 different answers:

mod.lme = lme(response ~ drug, random = ~1|source, dat)
mod.aov = aov(response ~ drug + Error(source), dat)

summary(mod.aov)

Error: source
         Df Sum Sq Mean Sq F value   Pr(>F)
drug       2 61.167  30.583  61.167 0.003703 **
Residuals  3  1.500   0.500
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
         Df Sum Sq Mean Sq F value Pr(>F)
Residuals  6    9.0     1.5

anova(mod.lme) # I use anova here to directly compare the F-test

           numDF denDF   F-value p-value
(Intercept)     1     6 115207.14  <.0001
drug            2     3     26.21  0.0126

(incidentally the 3 denDF here make me think the F-test is exactly what
I'd expect)

Because the results look different, I thought the possibilities are:

1) I fit 2 different models without realising it
2) one model is more conservative than the other
3) I'm completely missing some point (despite searching the archives of
R-help and R-ME)

Just to be pesky, if I check the calculations against the book I got the
data from (Zar 4th ed, pgg 304-305) they agree with the aov() results.

Any illumination is gratefully asked for. I apologise in advance for any
annoyance past/present/future my question will cause.

The gut reaction is that you shouldn't trust lme() in low-df cases, but in
this particular case the issue is different:

summary(mod.aov)

Error: source
        Df Sum Sq Mean Sq F value   Pr(>F)  drug       2 61.167  30.583
 61.167 0.003703 **
Residuals  3  1.500   0.500                   ---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
        Df Sum Sq Mean Sq F value Pr(>F)
Residuals  6    9.0     1.5
Notice that the Residuals Mean Sq is larger in the Within stratum than in
the source stratum. In terms of a mixed-effects model, this implies a
negative estimate for the variance of the source effect. lme() will have
nothing of that and sets it to zero instead. If you drop the Error(source)
you get the same F as in lme() although the df differ.

(The  "negative variance" can be interpreted as negative within-source
correlation, but that only works properly for balanced designs. Long
story...)

--
 O__  ---- Peter Dalgaard             ?ster Farimagsgade 5, Entr.B
 c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen   Denmark      Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)              FAX: (+45) 35327907