[R] proc mixed vs. lme
Grathwohl,Dominik,LAUSANNE,NRC/NT
dominik.grathwohl at rdls.nestle.com
Wed Oct 9 10:28:44 CEST 2002
Dear All,
Comparing linear mixed effect models in SAS and R, I found the following
discrepancy:
SAS R
random statement random subj(program); random = ~ 1 |
Subj
-2*loglik 1420.8 1439.363
random effects
variance(Intercept) 9.6033 9.604662
variance(residual) 1.1969 1.187553
the first 3 fixed effects
intercept 83.0952 81.10544
ProgramCont -3.4952 -1.11526
ProgramRI -1.9702 -1.04517
... ... ...
Can somebody explain me this different results?
thanks in advance,
Dominik
The two examples can be found in the internet:
http://ftp.sas.com/samples and http://cran.r-project.org
#My work around:
> R.version
_
platform i386-pc-mingw32
arch i386
os mingw32
system i386, mingw32
status
major 1
minor 6.0
year 2002
month 10
day 01
language R
#The SAS code from "http://ftp.sas.com/samples/A55235":
data weights;
input subj program$ s1 s2 s3 s4 s5 s6 s7;
datalines;
1 CONT 85 85 86 85 87 86 87
2 CONT 80 79 79 78 78 79 78
3 CONT 78 77 77 77 76 76 77
4 CONT 84 84 85 84 83 84 85
5 CONT 80 81 80 80 79 79 80
6 CONT 76 78 77 78 78 77 74
7 CONT 79 79 80 79 80 79 81
8 CONT 76 76 76 75 75 74 74
9 CONT 77 78 78 80 80 81 80
10 CONT 79 79 79 79 77 78 79
11 CONT 81 81 80 80 80 81 82
12 CONT 77 76 77 78 77 77 77
13 CONT 82 83 83 83 84 83 83
14 CONT 84 84 83 82 81 79 78
15 CONT 79 81 81 82 82 82 80
16 CONT 79 79 78 77 77 78 78
17 CONT 83 82 83 85 84 83 82
18 CONT 78 78 79 79 78 77 77
19 CONT 80 80 79 79 80 80 80
20 CONT 78 79 80 81 80 79 80
1 RI 79 79 79 80 80 78 80
2 RI 83 83 85 85 86 87 87
3 RI 81 83 82 82 83 83 82
4 RI 81 81 81 82 82 83 81
5 RI 80 81 82 82 82 84 86
6 RI 76 76 76 76 76 76 75
7 RI 81 84 83 83 85 85 85
8 RI 77 78 79 79 81 82 81
9 RI 84 85 87 89 88 85 86
10 RI 74 75 78 78 79 78 78
11 RI 76 77 77 77 77 76 76
12 RI 84 84 86 85 86 86 86
13 RI 79 80 79 80 80 82 82
14 RI 78 78 77 76 75 75 76
15 RI 78 80 77 77 75 75 75
16 RI 84 85 85 85 85 83 82
1 WI 84 85 84 83 83 83 84
2 WI 74 75 75 76 75 76 76
3 WI 83 84 82 81 83 83 82
4 WI 86 87 87 87 87 87 86
5 WI 82 83 84 85 84 85 86
6 WI 79 80 79 79 80 79 80
7 WI 79 79 79 81 81 83 83
8 WI 87 89 91 90 91 92 92
9 WI 81 81 81 82 82 83 83
10 WI 82 82 82 84 86 85 87
11 WI 79 79 80 81 81 81 81
12 WI 79 80 81 82 83 82 82
13 WI 83 84 84 84 84 83 83
14 WI 81 81 82 84 83 82 85
15 WI 78 78 79 79 78 79 79
16 WI 83 82 82 84 84 83 84
17 WI 80 79 79 81 80 80 80
18 WI 80 82 82 82 81 81 81
19 WI 85 86 87 86 86 86 86
20 WI 77 78 80 81 82 82 82
21 WI 80 81 80 81 81 82 83
;
/*---Data Set 3.2(b)---*/
data weight2;
set weights;
time=1; strength=s1; output;
time=2; strength=s2; output;
time=3; strength=s3; output;
time=4; strength=s4; output;
time=5; strength=s5; output;
time=6; strength=s6; output;
time=7; strength=s7; output;
keep subj program time strength;
run;
proc sort data=weight2;
by program time;
run;
/*---Data Set 3.2(c)---*/
proc means data=weight2 noprint;
by program time;
var strength;
output out=avg mean=strength;
run;
/*---produces Output 3.1 on pages 91-92---*/
proc mixed data=weight2;
class program subj time;
model strength = program time program*time /s;
random subj(program);
run;
#The SAS results:
Covariance Parameter
Estimates
Cov Parm Estimate
subj(program) 9.6033
Residual 1.1969
The Mixed Procedure
Fit Statistics
-2 Res Log Likelihood 1420.8
AIC (smaller is better) 1424.8
AICC (smaller is better) 1424.9
BIC (smaller is better) 1428.9
Solution for Fixed Effects
Standard
Effect program time Estimate Error DF t
Value Pr > |t|
Intercept 83.0952 0.7171 54
115.87 <.0001
program CONT -3.4952 1.0268 54
-3.40 0.0013
program RI -1.9702 1.0906 54
-1.81 0.0764
program WI 0 . .
. .
time 1 -2.0476 0.3376 324
-6.06 <.0001
time 2 -1.4286 0.3376 324
-4.23 <.0001
time 3 -1.1905 0.3376 324
-3.53 0.0005
time 4 -0.5714 0.3376 324
-1.69 0.0915
time 5 -0.4762 0.3376 324
-1.41 0.1594
time 6 -0.3810 0.3376 324
-1.13 0.2600
time 7 0 . .
. .
program*time CONT 1 2.1976 0.4834 324
4.55 <.0001
program*time CONT 2 1.7786 0.4834 324
3.68 0.0003
program*time CONT 3 1.5905 0.4834 324
3.29 0.0011
program*time CONT 4 1.0214 0.4834 324
2.11 0.0354
program*time CONT 5 0.6762 0.4834 324
1.40 0.1628
program*time CONT 6 0.3810 0.4834 324
0.79 0.4312
program*time CONT 7 0 . .
. .
program*time RI 1 0.6101 0.5134 324
1.19 0.2356
program*time RI 2 0.8661 0.5134 324
1.69 0.0926
program*time RI 3 0.8780 0.5134 324
1.71 0.0882
program*time RI 4 0.4464 0.5134 324
0.87 0.3852
program*time RI 5 0.6012 0.5134 324
1.17 0.2425
program*time RI 6 0.3810 0.5134 324
0.74 0.4586
program*time RI 7 0 . .
. .
program*time WI 1 0 . .
. .
program*time WI 2 0 . .
. .
program*time WI 3 0 . .
. .
program*time WI 4 0 . .
. .
program*time WI 5 0 . .
. .
program*time WI 6 0 . .
. .
program*time WI 7 0 . .
. .
#The R code:
library(nlme)
library(SASmixed)
options( contrasts = c(unordered = "contr.SAS", ordered = contr.poly"))
data(Weights)
fm1Weight <- lme( strength ~ Program * Time, data = Weights, random = ~ 1 |
Subj)
summary( fm1Weight )
VarCorr( fm1Weight )
#The R results:
> summary( fm1Weight )
Linear mixed-effects model fit by REML
Data: Weights
AIC BIC logLik
1455.363 1487.153 -719.6815
Random effects:
Formula: ~1 | Subj
(Intercept) Residual
StdDev: 3.099139 1.089749
Fixed effects: strength ~ Program * Time
Value Std.Error DF t-value p-value
(Intercept) 81.10544 0.7001315 339 115.84315 <.0001
ProgramCONT -1.11526 1.0024358 54 -1.11255 0.2708
ProgramRI -1.04517 1.0646834 54 -0.98168 0.3306
Time 0.15986 0.0224702 339 7.11447 <.0001
ProgramCONT:Time -0.18397 0.0321725 339 -5.71827 <.0001
ProgramRI:Time -0.05495 0.0341703 339 -1.60822 0.1087
Correlation:
(Intr) PrCONT PrgrRI Time PCONT:
ProgramCONT -0.698
ProgramRI -0.658 0.459
Time -0.225 0.157 0.148
ProgramCONT:Time 0.157 -0.225 -0.103 -0.698
ProgramRI:Time 0.148 -0.103 -0.225 -0.658 0.459
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.11669139 -0.60963281 0.03963523 0.63074983 3.33520406
Number of Observations: 399
Number of Groups: 57
> VarCorr( fm1Weight )
Subj = pdLogChol(1)
Variance StdDev
(Intercept) 9.604662 3.099139
Residual 1.187553 1.089749
Dominik Grathwohl
Biostatistician
Nestlé Research Center
PO Box 44, CH-1000 Lausanne 26
Phone: + 41 21 785 8034
Fax: + 41 21 785 8556
e-mail: dominik.grathwohl at rdls.nestle.com
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list