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
[EMAIL PROTECTED] 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]]
______________________________________________
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.
