[R] lme and mixed effects

[Doug Davidson]

Peter Dalgaard BSA writes:
 > Douglas Bates <bates at stat.wisc.edu> writes:
 > 
 > > Doug Davidson <dvdsn at printhost.beckman.uiuc.edu> writes:
 > > 
 > > > With lme, is there a way to specify multiple fixed factors under one
 > > > level of grouping?
 > > 
 > > I'm not sure I understand the question.  Can you tell us more about
 > > the model that you want to fit?
 > > 
 > > If you mean that you want a random interaction term then you might
 > > want to check the description in section 1.3 of Pinheiro and Bates
 > > (Springer, 2000) on different ways that interactions can be expressed
 > > in lme.
 > 
 > I read it as two crossed within-subject terms. lme() is not quite
 > built for that with its emphasis on nested random effects, but this
 > kind of stuff has worked for me before:
 > 
 > > lme((y~fact1*fact2,
 > +       random=list(subj=pdIdent(form=~fact1-1),
 > +                   subj=~1,
 > +                   fact2=~1),
 > 
 > (That implies a subj within subj grouping, but the algorithm treats
 > that just like subj.)
 > 

Yes, I have two crossed within-subject terms.  The design is for an
experiment where each subject sees all levels of two crossed factors.

Here is an example, and the results I get using aov and lme (note the
variable names have changed slightly from the above).  The y is a
continous response, x1 and x2 are fixed factors (both
within-subjects), and sbj is the random factor.

> mydata.gd
   sbj x1 x2    y
1    1  1  1  5.0
2    2  1  1  1.5
3    3  1  1  1.6
4    4  1  1  2.8
5    1  1  2  5.6
6    2  1  2  4.3
7    3  1  2 16.3
8    4  1  2  6.6
9    1  2  1 31.7
10   2  2  1 22.8
11   3  2  1 28.7
12   4  2  1 65.5
13   1  2  2 20.8
14   2  2  2 23.3
15   3  2  2 21.2
16   4  2  2 47.3

> fm1.aov <- aov(y~x1*x2+Error(sbj/(x1+x2)), data=mydata.gd)

> summary(fm1.aov)

Error: sbj
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  3    736     245               

Error: sbj:x1
          Df Sum Sq Mean Sq F value Pr(>F)  
x1         1   2958    2958    10.9  0.046 *
Residuals  3    814     271                 
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Error: sbj:x2
          Df Sum Sq Mean Sq F value Pr(>F)
x2         1   12.5    12.5    0.46   0.55
Residuals  3   81.3    27.1               

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)  
x1:x2      1  208.8   208.8    9.25  0.056 .
Residuals  3   67.7    22.6                 
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 



> fm1.lme <- lme(y~x1*x2, random=list(sbj=pdIdent(form=~x1-1), sbj=~1,
x2=~1), data=mydata.gd)

> anova(fm1.lme)
            numDF denDF F-value p-value
(Intercept)     1     6   0.970  0.3628
x1              1     6  32.545  0.0013
x2              1     3   0.219  0.6715
x1:x2           1     6   3.660  0.1043


I think I need to investigate pdIdent(form=...) further, but I'm unclear
why the degrees of freedom are allocated as they are in the lme anova.

In any case, thank you for your help!!

-Doug

Doug Davidson/2143 Beckman Institute/dvdsn at casper.beckman.uiuc.edu
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._