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

