[R-sig-ME] lme() vs aov()
Andrew Robinson
A.Robinson at ms.unimelb.edu.au
Thu May 22 23:12:37 CEST 2008
Hi all,
well spotted, Peter!
Federico, another tool that will help diagnose this problem is
VarCorr, which provides an estimate of the variance components.
> VarCorr(mod.lme)
source = pdLogChol(1)
Variance StdDev
(Intercept) 6.675305e-10 2.583661e-05
Residual 1.166667e+00 1.080123e+00
Andrew
On Thu, May 22, 2008 at 01:54:42PM -0700, Kingsford Jones wrote:
> 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
> >
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
Andrew Robinson
Department of Mathematics and Statistics Tel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/
More information about the R-sig-mixed-models
mailing list