[R] translating sas proc mixed to lme()

Fri Apr 6 07:49:27 CEST 2007

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