[R] translating sas proc mixed to lme()

From: <toby909_at_gmail.com>
Date: Fri 06 Apr 2007 - 05:49:27 GMT


Hi All

I am trying to translate a proc mixed into a lme() syntax. It seems that I was able to do it for part of the model, but a few things are still different.

It is a 2-level bivariate model (some call it a pseudo-3-level model).

PROC MIXED DATA=psdata.bivar COVTEST METHOD = ml; CLASS cluster_ID individual_id variable_id ; MODEL y = Dp Dq / SOLUTION NOINT;
RANDOM Dp Dq / SUBJECT = cluster_ID TYPE=UN G GCORR; REPEATED variable_id / SUBJECT = individual_ID(cluster_ID) TYPE=UN R RCORR; RUN; Here is my try:

dta = sqlQuery(odbcConnect("sasodbc", believeNRows=FALSE), "select * from psdata.bivar")
v = lme(Y~Dp+Dq-1, data=dta, random=~Dp+Dq-1|cluster_ID, method="ML", weights=varIdent(form=~1|Dp), corr=corCompSymm(, ~1|cluster_ID/individual_id))

Below is the data if you want to try, and also the outputs that I got from sas and R. What is different is the denominator DF for the F test for significance of the fixed effects.

Thanks for your hints.

Toby

DATA bivar;
INPUT cluster_ID individual_id variable_id $ Y; CARDS;

1 1 P 11.1
1 1 Q 20.6
1 2 P 13.3
1 2 Q 16.0
1 3 P 14.3
1 3 Q 15.1
2 4 P 18.6

2 4 Q 8.9
2 5 P 19.2
2 5 Q 7.3
2 6 P 21.5
2 6 Q 3.6
3 7 P 10.4
3 7 Q 2.7
3 8 P 11.1
3 8 Q 3.2
3 9 P 12.8
3 9 Q 5.0
;
run;
data psdata.bivar;
set bivar;
IF variable_ID eq "P" then Dp = 1; else Dp = 0; IF variable_ID eq "Q" then Dq = 1; else Dq = 0; run;

The SAS System
16:26 Friday, March 30, 2007 1

The Mixed Procedure

                  Model Information

Data Set                     WORK.BIVAR
Dependent Variable           Y
Covariance Structure         Unstructured
Subject Effects              cluster_ID,
                             individua(cluster_I)
Estimation Method            ML
Residual Variance Method     None
Fixed Effects SE Method      Model-Based
Degrees of Freedom Method Containment
                 Class Level Information

Class            Levels    Values

cluster_ID            3    1 2 3
individual_id         9    1 2 3 4 5 6 7 8 9
variable_id           2    P Q


            Dimensions

Covariance Parameters             6
Columns in X                      2
Columns in Z Per Subject          2
Subjects                          3
Max Obs Per Subject               6


          Number of Observations

Number of Observations Read              18
Number of Observations Used              18
Number of Observations Not Used           0


                     Iteration History

Iteration    Evaluations        -2 Log Like       Criterion

        0              1       109.94855540
        1              1        87.30141251      0.00000000


                   Convergence criteria met.


   Estimated R Matrix for
  individua(cluster_I) 1 1

Row Col1 Col2

   1      2.1822     -2.4739
   2     -2.4739      5.8522


   Estimated R Correlation
         Matrix for

  individua(cluster_I) 1 1

Row Col1 Col2

   1      1.0000     -0.6923
   2     -0.6923      1.0000


                Estimated G Matrix

                  cluster_
Row    Effect    ID              Col1        Col2

   1    Dp        1            12.4667     -2.3250
   2    Dq        1            -2.3250     32.1414


          Estimated G Correlation Matrix

                  cluster_
Row    Effect    ID              Col1        Col2

   1    Dp        1             1.0000     -0.1161
   2    Dq        1            -0.1161      1.0000


                        Covariance Parameter Estimates

                                                Standard         Z
Cov Parm    Subject                 Estimate       Error
Value        Pr Z

UN(1,1)     cluster_ID               12.4667     10.7811      1.16
0.1238
UN(2,1)     cluster_ID               -2.3250     12.3933     -0.19
0.8512
UN(2,2)     cluster_ID               32.1414     27.8589      1.15
0.1243
UN(1,1)     individua(cluster_I)      2.1822      1.2599      1.73
0.0416
UN(2,1)     individua(cluster_I)     -2.4739      1.7744     -1.39
0.1633
UN(2,2)     individua(cluster_I)      5.8522      3.3788      1.73
0.0416

           Fit Statistics

-2 Log Likelihood                87.3
AIC (smaller is better)         103.3
AICC (smaller is better)        119.3
BIC (smaller is better)          96.1


  Null Model Likelihood Ratio Test

    DF Chi-Square Pr > ChiSq

     5 22.65 0.0004

                 Solution for Fixed Effects

                      Standard
Effect    Estimate       Error      DF    t Value    Pr > |t|

Dp         14.7000      2.0971       2       7.01      0.0198
Dq          9.1556      3.3711       2       2.72      0.1130


       Type 3 Tests of Fixed Effects

           Num     Den
Effect      DF      DF    F Value    Pr > F

Dp           1       2      49.13    0.0198
Dq           1       2       7.38    0.1130







> summary(v)

Linear mixed-effects model fit by maximum likelihood   Data: dta

        AIC BIC logLik
   103.3014 110.4244 -43.65071

Random effects:
  Formula: ~Dp + Dq - 1 | cluster_ID
  Structure: General positive-definite, Log-Cholesky parametrization

          StdDev   Corr
Dp       3.530817 Dp
Dq       5.669334 -0.116

Residual 1.477236

Correlation Structure: Compound symmetry   Formula: ~1 | cluster_ID/individual_id   Parameter estimate(s):

       Rho
-0.692262
Variance function:
  Structure: Different standard deviations per stratum   Formula: ~1 | Dp
  Parameter estimates:

        1 0
1.000000 1.637610
Fixed effects: Y ~ Dp + Dq - 1

        Value Std.Error DF t-value p-value Dp 14.700000 2.224360 14 6.608641 0.0000 Dq 9.155556 3.575547 14 2.560603 0.0226   Correlation:
    Dp
Dq -0.149

Standardized Within-Group Residuals:

        Min Q1 Med Q3 Max -1.4002840 -0.5467747 -0.2042444 0.7126606 1.6044986

Number of Observations: 18
Number of Groups: 3



R-help@stat.math.ethz.ch 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 Apr 06 15:58:22 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Fri 06 Apr 2007 - 07:30:55 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.