[R-sig-ME] Trouble Replicating Unstructured Mixed Procedure in R

David Atkins datkins at u.washington.edu
Wed Jan 25 20:39:56 CET 2012


Hi Charles--

So, I don't know SAS at all, so just providing the SAS code and anova 
table don't tell me (at least) all that much about what SAS is doing. 
Also keep in mind that the anova table summary will depend on the "type" 
of sums of squares (and there have been some diatribes about this and 
type III sums of squares in the past on this list).

I find it a bit more informative to look at the summary output from lme(r):

 > y=lme(distance~Sex*age, random=(~1|Subject), data=dental)
 >
 > summary(y)
Linear mixed-effects model fit by REML
  Data: dental
        AIC      BIC    logLik
   445.7572 461.6236 -216.8786

Random effects:
  Formula: ~1 | Subject
         (Intercept) Residual
StdDev:    1.816214 1.386382

Fixed effects: distance ~ Sex * age
                   Value Std.Error DF   t-value p-value
(Intercept)   16.340625 0.9813122 79 16.651810  0.0000
SexFemale      1.032102 1.5374208 25  0.671321  0.5082
age            0.784375 0.0775011 79 10.120823  0.0000
SexFemale:age -0.304830 0.1214209 79 -2.510520  0.0141
  Correlation:
               (Intr) SexFml age
SexFemale     -0.638
age           -0.869  0.555
SexFemale:age  0.555 -0.869 -0.638

Standardized Within-Group Residuals:
         Min          Q1         Med          Q3         Max
-3.59804400 -0.45461690  0.01578365  0.50244658  3.68620792

Number of Observations: 108
Number of Groups: 27

DAVE: so, we fit a random-intercept, residual error, and 4 fixed-effects 
(N = 27 people and 108 observations)

Here is your second model:

 > y=lme(distance~Sex*age, random=(~1|Subject), corr=corSymm(,~1|Subject),
+ data=dental)
 > summary(y)
Linear mixed-effects model fit by REML
  Data: dental
        AIC      BIC    logLik
   450.1706 481.9033 -213.0853

Random effects:
  Formula: ~1 | Subject
         (Intercept) Residual
StdDev:    1.827112 1.375353

Correlation Structure: General
  Formula: ~1 | Subject
  Parameter estimate(s):
  Correlation:
   1      2      3
2 -0.174
3 -0.002 -0.177
4 -0.342  0.306  0.227
Fixed effects: distance ~ Sex * age
                   Value Std.Error DF   t-value p-value
(Intercept)   15.932927 0.9978753 79 15.966852  0.0000
SexFemale      1.473662 1.5633701 25  0.942619  0.3549
age            0.824339 0.0824260 79 10.000962  0.0000
SexFemale:age -0.348100 0.1291367 79 -2.695594  0.0086
  Correlation:
               (Intr) SexFml age
SexFemale     -0.638
age           -0.874  0.558
SexFemale:age  0.558 -0.874 -0.638

Standardized Within-Group Residuals:
          Min           Q1          Med           Q3          Max
-3.180968424 -0.544002662  0.001195789  0.495264398  3.730242930

Number of Observations: 108
Number of Groups: 27

DAVE: Now, in addition to earlier parameters, we have an unstructured 
covariance matrix for repeated measures.  By all indications this is a 
*worse* fitting model...

Finally, I wonder whether SAS is fitting random-effects at all (no idea, 
just a guess).  If so, you might check it relative to a gls() fit, a la:

 > y2=gls(distance~Sex*age, corr=corSymm(,~1|Subject),
+ data=dental)
 > summary(y2)
Generalized least squares fit by REML
   Model: distance ~ Sex * age
   Data: dental
        AIC      BIC    logLik
   448.1706 477.2589 -213.0853

Correlation Structure: General
  Formula: ~1 | Subject
  Parameter estimate(s):
  Correlation:
   1     2     3
2 0.575
3 0.638 0.574
4 0.515 0.749 0.721

Coefficients:
                   Value Std.Error   t-value p-value
(Intercept)   15.932911 0.9978729 15.966874  0.0000
SexFemale      1.473679 1.5633664  0.942632  0.3481
age            0.824340 0.0824258 10.001000  0.0000
SexFemale:age -0.348102 0.1291364 -2.695611  0.0082

  Correlation:
               (Intr) SexFml age
SexFemale     -0.638
age           -0.874  0.558
SexFemale:age  0.558 -0.874 -0.638

Standardized residuals:
         Min          Q1         Med          Q3         Max
-2.41707752 -0.64439706 -0.07388955  0.58480683  2.26288228

Residual standard error: 2.286908
Degrees of freedom: 108 total; 104 residual

So, not sure if this is helpful, but perhaps you could provide 
additional output / commentary on what exactly SAS is doing (beyond 
anova table).

cheers, Dave


Greetings,

I have been working on R for some time now and I have begun the endeavor of
trying to replicate some SAS code in R.  I have scoured the forums but
haven't been able to find an answer.  I hope one of you could be so kind as
to enlighten me.

I am attempting to replicate a repeated measures experiment using some
standard data.  I have posted the SAS code and output directly from a
publication as well as my attempts in R to replicate it.  My main issue
comes with the 'unstructured' component.

The 'dental' dataset from 'mixedQF' package,
equivalent to formixed data in SAS

     distance age Subject    Sex
1       26.0   8     M01   Male
2       25.0  10     M01   Male
3       29.0  12     M01   Male
4       31.0  14     M01   Male
5       21.5   8     M02   Male
6       22.5  10     M02   Male
7       23.0  12     M02   Male
8       26.5  14     M02   Male
9       23.0   8     M03   Male
10      22.5  10     M03   Male
11      24.0  12     M03   Male
12      27.5  14     M03   Male
13      25.5   8     M04   Male
14      27.5  10     M04   Male
15      26.5  12     M04   Male
16      27.0  14     M04   Male
17      20.0   8     M05   Male
18      23.5  10     M05   Male
19      22.5  12     M05   Male
20      26.0  14     M05   Male
21      24.5   8     M06   Male
22      25.5  10     M06   Male
23      27.0  12     M06   Male
24      28.5  14     M06   Male
25      22.0   8     M07   Male
26      22.0  10     M07   Male
27      24.5  12     M07   Male
28      26.5  14     M07   Male
29      24.0   8     M08   Male
30      21.5  10     M08   Male
31      24.5  12     M08   Male
32      25.5  14     M08   Male
33      23.0   8     M09   Male
34      20.5  10     M09   Male
35      31.0  12     M09   Male
36      26.0  14     M09   Male
37      27.5   8     M10   Male
38      28.0  10     M10   Male
39      31.0  12     M10   Male
40      31.5  14     M10   Male
41      23.0   8     M11   Male
42      23.0  10     M11   Male
43      23.5  12     M11   Male
44      25.0  14     M11   Male
45      21.5   8     M12   Male
46      23.5  10     M12   Male
47      24.0  12     M12   Male
48      28.0  14     M12   Male
49      17.0   8     M13   Male
50      24.5  10     M13   Male
51      26.0  12     M13   Male
52      29.5  14     M13   Male
53      22.5   8     M14   Male
54      25.5  10     M14   Male
55      25.5  12     M14   Male
56      26.0  14     M14   Male
57      23.0   8     M15   Male
58      24.5  10     M15   Male
59      26.0  12     M15   Male
60      30.0  14     M15   Male
61      22.0   8     M16   Male
62      21.5  10     M16   Male
63      23.5  12     M16   Male
64      25.0  14     M16   Male
65      21.0   8     F01 Female
66      20.0  10     F01 Female
67      21.5  12     F01 Female
68      23.0  14     F01 Female
69      21.0   8     F02 Female
70      21.5  10     F02 Female
71      24.0  12     F02 Female
72      25.5  14     F02 Female
73      20.5   8     F03 Female
74      24.0  10     F03 Female
75      24.5  12     F03 Female
76      26.0  14     F03 Female
77      23.5   8     F04 Female
78      24.5  10     F04 Female
79      25.0  12     F04 Female
80      26.5  14     F04 Female
81      21.5   8     F05 Female
82      23.0  10     F05 Female
83      22.5  12     F05 Female
84      23.5  14     F05 Female
85      20.0   8     F06 Female
86      21.0  10     F06 Female
87      21.0  12     F06 Female
88      22.5  14     F06 Female
89      21.5   8     F07 Female
90      22.5  10     F07 Female
91      23.0  12     F07 Female
92      25.0  14     F07 Female
93      23.0   8     F08 Female
94      23.0  10     F08 Female
95      23.5  12     F08 Female
96      24.0  14     F08 Female
97      20.0   8     F09 Female
98      21.0  10     F09 Female
99      22.0  12     F09 Female
100     21.5  14     F09 Female
101     16.5   8     F10 Female
102     19.0  10     F10 Female
103     19.0  12     F10 Female
104     19.5  14     F10 Female
105     24.5   8     F11 Female
106     25.0  10     F11 Female
107     28.0  12     F11 Female
108     28.0  14     F11 Female

*Mixed modeling and fixed effect test*
SAS
proc mixed data=formixed;
class gender age person;
model y = gender|age;
repeated / type=cs sub=person;
run;

output of interest to me
           Tests of Fixed Effects
Source             NDF   DDF    Type III F    Pr > F
GENDER           1        25        9.29        0.0054
AGE                  3        75       35.35       0.0001
GENDER*AGE   3        75        2.36        0.0781

R (nlme package)
y=lme(distance~Sex*age, random=(~1|Subject), data=dental)
anova(y)

             numDF denDF  F-value p-value
(Intercept)     1    75 4123.156  <.0001
Sex              1    25    9.292  0.0054
age               3    75   40.032  <.0001
Sex:age        3    75    2.362  0.0781

Now this isn't exact but it is extremely close, however when I try to
replicate the unstructured,

SAS
proc mixed data=formixed;
class gender age person;
model y = gender|age;
repeated / type=un sub=person;
run;

              Tests of Fixed Effects
Source          NDF DDF Type III F Pr > F
GENDER         1    25     9.29    0.0054
AGE                3    25    34.45   0.0001
GENDER*AGE 3    25     2.93    0.0532

R
either
y=lme(distance~Sex*age, random=(~1|Subject), corr=corSymm(,~1|Subject),
data=dental)
anova(y)
or
z=lme(distance~Sex*age, random=(~1|Subject), corr=corSymm(), data=dental)
anova(z)

gives the output

             numDF denDF  F-value    p-value
(Intercept)     1    75     4052.028  <.0001
Sex              1    25       8.462      0.0075
age               3    75      39.022    <.0001
Sex:age        3    75       2.868      0.0421

What am I doing wrong to replicate the unstructured linear mixed model from
SAS?

Regards,

Charles

	[[alternative HTML version deleted]]

-- 
Dave Atkins, PhD
Research Associate Professor
Department of Psychiatry and Behavioral Science
University of Washington
datkins at u.washington.edu

Center for the Study of Health and Risk Behaviors (CSHRB)		
1100 NE 45th Street, Suite 300 	
Seattle, WA  98105 	
206-616-3879 	
http://depts.washington.edu/cshrb/
(Mon-Wed)	

Center for Healthcare Improvement, for Addictions, Mental Illness,
   Medically Vulnerable Populations (CHAMMP)
325 9th Avenue, 2HH-15
Box 359911
Seattle, WA 98104
http://www.chammp.org
(Thurs)




More information about the R-sig-mixed-models mailing list