Re: [R] mixed model results from SAS and R

From: Aldi Kraja <aldi_at_wustl.edu>
Date: Thu, 22 May 2008 19:21:56 -0500

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_at_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_at_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]]

______________________________________________
R-help_at_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.
Received on Fri 23 May 2008 - 00:33:57 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 23 May 2008 - 01:30:40 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive