[R] Help with lme basics
Greenberg, Jeff (J.A.)
jgreenb2 at ford.com
Fri Apr 19 15:31:58 CEST 2002
Renaud,
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list