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>
Hi,
To continue the test of the differences between mixed model SAS (9.1.3)
and lme (R 7.0), I removed the intercept of the random effect in the R
case. In that case lme is showing that beta coefficients are almost all
of them identical (-0.08424262) and somehow all of them are spread along
the plot of the coefficients axes (!!!) which surprises me. Unless the
authors of lme use jittering of the points. Is lme having a bug in this
case? I know that the "a1" data are not very good. If one looks at
the
distribution of the residuals they are not beautifully iid. But lme had
to detect this problem of random beta-s being almost identical and close
to zero. In addition is there any new development on mixed models in R
for the quantitative traits that is not covered by lme package, but can
help to understand this case?
fm2nidl.lme<-lme(nidl~time,data=a1,random=~time-1 | sub)
plot(coef(fm2nidl.lme))
print(coef(fm2nidl.lme))
anova(fm1nidl.lme,fm2nidl.lme)
Model df AIC BIC logLik Test L.Ratio p-value
fm1nidl.lme 1 6 4982.150 5009.958 -2485.075
fm2nidl.lme 2 4 5069.081 5087.619 -2530.540 1 vs 2 90.9309 <.0001
R output (contd.)
=========== (Intercept) time
1 7.782863 -0.08424262
2 7.782863 -0.08424262
3 7.782863 -0.08424262
4 7.782863 -0.08424262
5 7.782863 -0.08424262
6 7.782863 -0.08424262
7 7.782863 -0.08424262
8 7.782863 -0.08424262
9 7.782863 -0.08424262
10 7.782863 -0.08424262
...
Thank you in advance for any insights,
Aldi
aldi@dsgmail.wustl.edu wrote:> 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.
>
>
>
> ------------------------------------------------------------------------
>
> ______________________________________________
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
[[alternative HTML version deleted]]