R help - Oct 2003 - mixed effects with nlme

Dear R users:

 I have some difficulties analizing data with mixed effects NLME and the
last version of R. More concretely, I have a repeated measures design with
a single group and 2 experimental factors (say A and B) and my interest is
to compare additive and nonadditive models. 

   suj  rv    A       B
1   s1   4   a1      b1
2   s1   5   a1      b2
3   s1   7   a1      b3
4   s1   1   a2      b1
5   s1   4   a2      b2
6   s1   2   a2      b3
7   s2   6   a1      b1
8   s2   8   a1      b2
9   s2  10   a1      b3
10  s2   3   a2      b1
11  s2   6   a2      b2
12  s2   6   a2      b3
13  s3   1   a1      b1
14  s3   6   a1      b2
15  s3   5   a1      b3
16  s3   3   a2      b1
17  s3   5   a2      b2
18  s3   4   a2      b3
19  s4   2   a1      b1
20  s4  10   a1      b2
21  s4  12   a1      b3
22  s4   1   a2      b1
23  s4   4   a2      b2
24  s4   7   a2      b3
25  s5   5   a1      b1
26  s5  10   a1      b2
27  s5  10   a1      b3
28  s5   5   a2      b1
29  s5   6   a2      b2
30  s5   5   a2      b3
31  s6   1   a1      b1
32  s6   7   a1      b2
33  s6   8   a1      b3
34  s6   2   a2      b1
35  s6   8   a2      b2
36  s6   7   a2      b3

It is very easy to fit these data with base R function AOV:

NonAdditive model:
 aov(rv ~ A*B + Error(suj+suj/A+suj/B)

Additive model:
 aov(rv ~ A*B + Error(suj)

and also easy with SAS MIXED (I missed some obvious lines):

NonAdditive model
 model vr = A B A*B;
 random suj A*suj B*suj;
 repeated / type=cs subj=suj;

Additive model;
 model vr = A B A*B /ddfm=satterth;
 repeated / type=cs subj=suj;
 
Using LME I do not find any problems to fit the additive model with

 lme(vr~A*B, random=~1|suj, cor=corCompSymm())

but I have found some difficulties fitting the nonadditive model.

 Can anyone help me?
 
 Thanks in advance.

	Manuel Ato
	Dpto. Psic.B?sica y Metodolog?a
	Apartado 4021
	30080 MURCIA (Spain)
	e-mail: matogar at um.es

1.  Have you reviewed Pinhiero and Bates (2000) Mixed Effects Models in 
S and S-Plus (Springer)?  I think the answer is there. 

2.  Observing the philosophy that a quick response that may contain an 
error is often better than a perfect but late response, I will offer a 
suggestion I have not checked:  Have you considered "random = 
~1+A+B|subj"? 

hope this helps.  spencer graves

Manuel Ato Garcia wrote:
>Dear R users:
>
> I have some difficulties analizing data with mixed effects NLME and the
>last version of R. More concretely, I have a repeated measures design with
>a single group and 2 experimental factors (say A and B) and my interest is
>to compare additive and nonadditive models. 
>
>   suj  rv    A       B
>1   s1   4   a1      b1
>2   s1   5   a1      b2
>3   s1   7   a1      b3
>4   s1   1   a2      b1
>5   s1   4   a2      b2
>6   s1   2   a2      b3
>7   s2   6   a1      b1
>8   s2   8   a1      b2
>9   s2  10   a1      b3
>10  s2   3   a2      b1
>11  s2   6   a2      b2
>12  s2   6   a2      b3
>13  s3   1   a1      b1
>14  s3   6   a1      b2
>15  s3   5   a1      b3
>16  s3   3   a2      b1
>17  s3   5   a2      b2
>18  s3   4   a2      b3
>19  s4   2   a1      b1
>20  s4  10   a1      b2
>21  s4  12   a1      b3
>22  s4   1   a2      b1
>23  s4   4   a2      b2
>24  s4   7   a2      b3
>25  s5   5   a1      b1
>26  s5  10   a1      b2
>27  s5  10   a1      b3
>28  s5   5   a2      b1
>29  s5   6   a2      b2
>30  s5   5   a2      b3
>31  s6   1   a1      b1
>32  s6   7   a1      b2
>33  s6   8   a1      b3
>34  s6   2   a2      b1
>35  s6   8   a2      b2
>36  s6   7   a2      b3
>
>It is very easy to fit these data with base R function AOV:
>
>NonAdditive model:
> aov(rv ~ A*B + Error(suj+suj/A+suj/B)
>
>Additive model:
> aov(rv ~ A*B + Error(suj)
>
>and also easy with SAS MIXED (I missed some obvious lines):
>
>NonAdditive model
> model vr = A B A*B;
> random suj A*suj B*suj;
> repeated / type=cs subj=suj;
>
>Additive model;
> model vr = A B A*B /ddfm=satterth;
> repeated / type=cs subj=suj;
> 
>Using LME I do not find any problems to fit the additive model with
>
> lme(vr~A*B, random=~1|suj, cor=corCompSymm())
>
>but I have found some difficulties fitting the nonadditive model.
>
> Can anyone help me?
> 
> Thanks in advance.
>
>	Manuel Ato
>	Dpto. Psic.B?sica y Metodolog?a
>	Apartado 4021
>	30080 MURCIA (Spain)
>	e-mail: matogar at um.es
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>  
>

Hi, R-users:

 Last week I send a request for help to this list. I have receive until now
two kindly responses from Spencer Graves and Renauld Lancelot. They both 
point interesting things to fit an adequate model to my data but
unfortunately 
it persists without a satisfactory solution. 

 I propose again the same problem but using with a different dataset
(Assay), taken from Pinheiro and Bates' book on page 163, that is relevant
with crossed 
random effects. I have fitted the same model (p. 165)
>fmAssay <- lme(logDens ~ sample * dilut, Assay, random=pdBlocked(list(,
         pdIdent(~1), pdIdent(~sample-1),pdIdent(~dilut-1))) )

and the results with "anova" function (p. 166) are
 
             numDF denDF  F-value p-value
(Intercept)      1    29 537.6294  <.0001
sample           5    29  11.2103  <.0001
dilut            4    29 420.5458  <.0001
sample:dilut    20    29   1.6072  0.1192

 The problem is that with this approach one obtains correct F-values, but 
using a common residual term for DenDF that is a combination of (Block +
Block:sample + Block:dilut). Then probability values are different to that
obtained when we used the classical AOV funtion to fit the same model,
because in this case each term is mapped with a error term (so
"sample"
uses "Block:sample", "dilut" uses "Block:dilut",
and "sample:dilut" uses
"Block:sample:dilut"):
>mod<-aov(logDens ~ sample*dilut + Error(Block+Block/sample+Block/dilut),
data=Assay)
>summary(mod)
Error: Block
          Df    Sum Sq   Mean Sq F value Pr(>F)
Residuals  1 0.0083115 0.0083115               

Error: Block:sample
          Df   Sum Sq  Mean Sq F value   Pr(>F)
sample     5 0.276153 0.055231  11.213 0.009522
Residuals  5 0.024627 0.004925                 

Error: Block:dilut
          Df Sum Sq Mean Sq F value    Pr(>F)
dilut      4 3.7491  0.9373  420.79 1.684e-05
Residuals  4 0.0089  0.0022                  

Error: Within
             Df   Sum Sq  Mean Sq F value Pr(>F)
sample:dilut 20 0.055525 0.002776  1.6069 0.1486
Residuals    20 0.034555 0.001728  


 Obviously, the interest on linear mixed effects is open with the
possibility of modeling with correlated and/or heterocedastic errors, and
this end cannot
be pursued with AOV function.

 Summarizing, the problem is that I have not found a way to obtain with
NLME the same results (DF, F-ratios and probabilities) that I get with AOV and
multistratum errors. Is this an inconvenience of program?, probably due
to the impossibility to use multiple nested arguments as 

 random(~1|Block/sample+dilut) or  random(~1|Block/sample*dilut)
 
SAS MIXED can also fit these data and obtain correct results by means of a
combination of "random" and "repeated" arguments:

 model = sample dilut sample*dilut;
 random = Block sample*Block dilut*Block;
 repeated /type=cs Sub=Block;


              Type 3 Tests of Fixed Effects

                      Num     Den
            Effect     DF      DF    F Value    Pr > F
            sample      5       5      11.21    0.0095
             dilut      4       4     420.79    <.0001
      sample*dilut     20      20       1.61    0.1486


May be possible obtain the same results with NLME function?

 I would appreciate any kind of help.

 Best regards,


			Manuel Ato
			University of Murcia
			Spain
			e-mail: matogar at um.es