[R] mixed model resuts from SAS and R
aldi at dsgmail.wustl.edu
aldi at dsgmail.wustl.edu
Thu May 22 14:08:40 CEST 2008
Hi,
I was wondering if there is a way to figure out why in SAS random beta
coefficients are 0 vs. in R the beta-s are non zero.
The variables of the data are nidl, time, and sub (for subject). Time and
nidl are continuous variables. I am applying random coefficients model.
Any input is greatly appreciated, Thanks, Aldi
1. mixed model in SAS:
======================
ods output SolutionR = out1.randomnidltest2;
proc mixed data = a1 ;
class sub ;
model nidl = time / solution ;
random int time / sub = sub solution;
run;
ods output close;
2. mixed model in R:
====================
a1<-read.table(file="c:\\aldi\\a1.txt",sep=",",header=T)
library(nlme)
fm1nidl.lme<-lme(nidl~time,data=a1,random=~time | sub)
plot(coef(fm1nidl.lme))
3. SAS output:
==============
Plot of nidl*time. Symbol used is '*'.
40
*
*
* *
* *
* * *
20 * * *
* *****
* ** **
* *** **** *
* *** ** *
* * ****** *
* ****** * *
* * ******* * *
0 * * ******** *
------------------------------------
0 25 50 75
time
NOTE: 684 obs hidden.
Dimensions
Covariance Parameters 3
Columns in X 2
Columns in Z Per Subject 2
Subjects 425
Max Obs Per Subject 2
Number of Observations
Number of Observations Read 763
Number of Observations Used 763
Number of Observations Not Used 0
Covariance Parameter Estimates
Cov Parm Subject Estimate
Intercept sub 17.1773
time sub 0
Residual 27.0023
The Mixed Procedure
Fit Statistics
-2 Res Log Likelihood 5005.5
AIC (smaller is better) 5009.5
AICC (smaller is better) 5009.5
BIC (smaller is better) 5017.6
Solution for Fixed Effects
Standard
Effect Estimate Error DF t Value Pr > |t|
Intercept 7.7608 0.3214 424 24.15 <.0001
time -0.08148 0.01605 337 -5.08 <.0001
Solution for Random Effects
Std Err
Effect sub Estimate Pred DF t Value Pr > |t|
Intercept 1 5.6722 3.2426 0 1.75 .
time 1 0 . . . .
Intercept 2 3.6722 3.2426 0 1.13 .
time 2 0 . . . .
Intercept 3 -2.0774 2.7539 0 -0.75 .
time 3 0 . . . .
... ... .... ... .... ...
R output:
(Intercept) time
1 17.432680 -0.3155496745
2 14.024527 -0.2345787274
3 3.105323 0.0469759240
4 23.047062 -0.5565796200
5 10.708267 -0.1557909941
... ... ...
> summary(fm1nidl.lme)
Linear mixed-effects model fit by REML
Data: a1
AIC BIC logLik
4982.15 5009.958 -2485.075
Random effects:
Formula: ~time | sub
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 6.0090637 (Intr)
time 0.1717771 -0.831
Residual 4.2885993
Fixed effects: nidl ~ time
Value Std.Error DF t-value p-value
(Intercept) 7.779379 0.3582945 424 21.71225 0
time -0.086206 0.0158615 337 -5.43494 0
Correlation:
(Intr)
time -0.675
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.0234047 -0.5132952 -0.2246854 0.4249250 3.5611259
Number of Observations: 763
Number of Groups: 425
3. data:
========
see the text file (comma delimited) attached.
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: a1.txt
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20080522/db9c1399/attachment.txt>
More information about the R-help
mailing list