[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