[R] Help with lme basics
Douglas Bates
bates at stat.wisc.edu
Fri Apr 19 17:34:49 CEST 2002
The lme model divides the terms in the fixed effects into those that
change within subjects and those that do not. Apparently both dose
and drug change within subject and only gp does not change within
subject. Thus the only term whose denominator degrees of freedom is
determined by the number of subjects is the gp term. All other terms
and interactions are estimated by all the data on all the subjects
hence they have 70 denominator degrees of freedom.
We explain the reasoning behind the lme model in chapter 3 of our book
@Book{pinh:bate:2000,
author = {Jos\'{e} C. Pinheiro and Douglas M. Bates},
title = {Mixed-Effects Models in \textsf{S} and \textsf{S-PLUS}},
publisher = {Springer},
year = 2000,
series = {Statistics and Computing}
}
It is derived from a hierarchical linear model with random effects at
the subject level.
I think this comparison shows that lme evaluates the significance of
terms differently from models that employ error strata. I don't think
it shows that the lme model is "wrong".
I will leave it to someone who understands error strata better than I
do to explain the rationale of the results based on error strata.
"Greenberg, Jeff (J.A.)" <jgreenb2 at ford.com> writes:
> Thanks. Im not a statistician so maybe I have this all wrong but here's what I was expecting. When I model this data set with:
>
> summary(aov(effect ~ gp * drug * dose + Error(subj/(dose+drug)), data=Ela.uni))
>
> what I mean is that subj is a random factor nested within gp and that dose and drug are fixed effects. So I want to construct tests that have the fixed effects crossed with random effects in the denominator. The output of aov is exactly what I would get from fitting this model in SYSTAT and the error strata and error df's are exactly what I would compute by hand:
>
> Error: subj
> Df Sum Sq Mean Sq F value Pr(>F)
> gp 1 270.01 270.01 7.0925 0.01855 *
> Residuals 14 532.98 38.07
> ---
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
>
> Error: subj:dose
> Df Sum Sq Mean Sq F value Pr(>F)
> dose 2 758.77 379.39 36.5097 1.580e-08 ***
> gp:dose 2 42.27 21.14 2.0339 0.1497
> Residuals 28 290.96 10.39
> ---
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
>
> Error: subj:drug
> Df Sum Sq Mean Sq F value Pr(>F)
> drug 1 348.84 348.84 13.001 0.002866 **
> gp:drug 1 326.34 326.34 12.163 0.003624 **
> Residuals 14 375.65 26.83
> ---
> Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
>
> Error: Within
> Df Sum Sq Mean Sq F value Pr(>F)
> drug:dose 2 12.062 6.031 0.6815 0.5140
> gp:drug:dose 2 14.812 7.406 0.8369 0.4436
> Residuals 28 247.792 8.850
>
> If I then follow your suggestion for fitting with lme the results I get are confusing:
>
> > fm <- lme(fixed = effect ~ gp * drug * dose,
> + random = ~ 1 | subj,
> + data = Ela.uni)
> > anova(fm)
> numDF denDF F-value p-value
> (Intercept) 1 70 1461.2979 <.0001
> gp 1 14 7.0924 0.0185
> drug 1 70 26.7052 <.0001
> dose 2 70 29.0433 <.0001
> gp:drug 1 70 24.9828 <.0001
> gp:dose 2 70 1.6180 0.2056
> drug:dose 2 70 0.4617 0.6321
> gp:drug:dose 2 70 0.5670 0.5698
>
> Here the gp term is similar to what I expect (at least the number of degrees of freedom in the denominator is right) but the rest of the table makes no sense to me (70 df's in the denominator???).
>
> I have been assuming that for a balanced case like this I should be able to get lme to report results that are very similar to the results of a 'traditional' analysis like aov. Until I can figure this out, I have no confidence that I can use lme for the unbalanced data I'm really interested in.
>
>
> -----Original Message-----
> From: Renaud Lancelot [mailto:lancelot at sentoo.sn]
> Sent: Thursday, April 18, 2002 2:52 PM
> To: Greenberg, Jeff (J.A.)
> Subject: Re: [R] Help with lme basics
>
>
> Hi Jeff,
>
> "Greenberg, Jeff (J.A.)" wrote:
> >
> > In Baron and Li's "Notes on the use of R for psychology experiments and questionnaires" http://cran.r-project.org/doc/contrib/rpsych.htm they describe a balanced data set for a drug experiment:
> >
> > "... a test of drug treatment effect by one between-subject factor: group (two groups of 8 subjects each) and two within-subject factors: drug (2 levels) and dose (3 levels). "
> >
> > This design is identical to an unbalanced one that I am interested in analyzing in R. I have both missing cells and unequal number of observations per cell. This is a repeated-measures design with test subject as a random factor.
> >
> > Baron and Li show how to analyze this design using aov:
> >
> > summary(aov(effect ~ gp * drug * dose + Error(subj/(dose+drug)), data=Ela.uni))
> >
> > When I analyze my design I get additional terms in each error stratum because of the missing cells. The results appear sensible, but I've been told that lme is a better way to model this data.
> >
> > My question is: how do I set up an equivalent model using lme? Using the balanced data set in Baron and Li I've tried this model:
> >
> > ela.lme<-lme(effect~gp*dose*drug,random=~1|subj/drug/dose)
>
> I don't thing dose and drugs can be considered as random effects like
> you specified it in the above model. They are both within-subject fixed
> effects. So:
>
> fm <- lme(fixed = effect ~ gp * drug * dose,
> random = ~ 1 | subj,
> data = Ela.uni)
>
> will give you population means for gp, drug and dose and interaction
> terms, and a subject-specific baseline (subject random effect).
>
> If you have strong evidence (e.g. graphical exploration of residuals) of
> an heterogeneous residuals variance according to strata defined by drug
> and dose, this can be specified with the weights argument:
>
> fm <- lme(fixed = effect ~ gp * drug * dose,
> random = ~ 1 | subj,
> weights = varIdent(form = ~ 1 | drug * dose),
> data = Ela.uni)
>
> This will give you one residuals variance for each combination of drug
> by dose.
>
> To examine the numeric results, you can use the summary, anova and
> intervals methods for lme objects. For example:
> anova(fm, type = "marginal")
> will produce an anova table for the fixed effects;
> intervals(fm)
> will produce confidence intervals for fixed and random parameters
> (random effects and variance parameters).
>
> A very useful and nice reference for mixed-effects models with lme is
> Pinheiro, J.C., Bates, D.M., 2000. Mixed-effect models in S and S-Plus.
> New York (USA), Springer, 598 p.
> which was written by the authors of nlme.
>
> Hope this helps,
>
> Renaud
>
> --
> Dr Renaud Lancelot, vétérinaire
> CIRAD, Département Elevage et Médecine Vétérinaire (CIRAD-Emvt)
> Programme Productions Animales
> http://www.cirad.fr/presentation/programmes/prod-ani.shtml (Français)
> http://www.cirad.fr/presentation/en/program-eng/prod-ani.shtml (English)
>
> ISRA-LNERV tel (221) 832 49 02
> BP 2057 Dakar-Hann fax (221) 821 18 79 (CIRAD)
> Senegal e-mail renaud.lancelot at cirad.fr
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
--
Douglas Bates bates at stat.wisc.edu
Statistics Department 608/262-2598
University of Wisconsin - Madison http://www.stat.wisc.edu/~bates/
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list