[R-sig-ME] lme() vs aov()

[Kingsford Jones]

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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>