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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._