Dear Charles,

First, I hope you are not put off by the tone of some the responses in this
email thread.  Sometimes people have strong opinions and it might come
across aggressively.

Second, be sure to understand that reproducing a SAS analysis with lme in
no way violates any legal agreements that SAS may have, if for no other
reason than you never signed an agreement with SAS!  That bit in the EULA
about decompiling and reverse engineering means that people are prohibited
from creating a new version of PROC MIXED that does the same thing.  The
nlme package uses different methods than SAS. E.g. different optimizers,
even uses a log-parameterization deep in the code so that negative variance
components cannot happen.

Third, by now you've probably figured out that PROC MIXED and lme have very
different ideas about degrees of freedom.  Also, the loglikelihoods are on
different scales.  For that reason, when I try to reproduce an analysis, I
find the best way to compare is to look at the variance components.

Here is what PROC MIXED says:

    Cov Parm Estimate     Std Error       Z  Pr > |Z|
UN(1,1)     5.41545455    1.53172185    3.54    0.0004
UN(2,1)     2.71681818    1.09623989    2.48    0.0132
UN(2,2)     4.18477273    1.18363247    3.54    0.0004
UN(3,1)     3.91022727    1.41775367    2.76    0.0058
UN(3,2)     2.92715909    1.19304751    2.45    0.0141
UN(3,3)     6.45573864    1.82595863    3.54    0.0004
UN(4,1)     2.71022727    1.17209851    2.31    0.0208
UN(4,2)     3.31715909    1.12903016    2.94    0.0033
UN(4,3)     4.13073864    1.40356157    2.94    0.0033
UN(4,4)     4.98573864    1.41017984    3.54    0.0004
Residual         1.00000000 .       .         .

That's our target, and here is the lme code to get us there.

library(nlme)
orth <- as.data.frame(Orthodont)
m1 <- gls(distance ~ Sex*age,
          correlation=corSymm(form = ~ 1 | Subject),
          weights = varIdent(form = ~ 1 | age),
          data = orth)
cors <- corMatrix(m1$modelStruct$corStruct)[[1]]
vars <- c(1.0000000, 0.8788793, 1.0744596, 0.9586878)^2 * 2.329213^2
covs <- outer(vars,vars, function(x,y) sqrt(x)*sqrt(y))
cors*covs

Now, some explanations will surely help, so let's step through the code
with comments.

The Orthodont data is a groupedData object, which can be helpful, but
sometimes confusing, so coercing to a data.frame removes the formula--you
can see it with 'formula(Orthodont)'.

orth <- as.data.frame(Orthodont)

Since there is no "random" line in the PROC MIXED code, we don't want to
use 'lme' (also why I removed the formula from the data), but instead use
'gls' for generalized least squares.  The 'corSymm' part specifies a
symmetric correlation matrix with 1 on the diagonal.  Looking at the MIXED
parameters, the results are given a covariances, not correlations.  Note
how the variances UN(1,1), UN(2,2,), etc are different, not constant along
the diagonal.  We use the 'weights' statement to specify different stratum
variances.

m1 <- gls(distance ~ Sex*age,
          correlation=corSymm(form = ~ 1 | Subject),
          weights = varIdent(form = ~ 1 | age),
          data = orth)

Now we have to convert the correlations and standard deviations of lme into
covariance parameters of MIXED.

First, extract the correlation matrix from lme.

cors <- corMatrix(m1$modelStruct$corStruct)[[1]]

Now, hard-code the stratum std deviations and square them to get
variances.  Multiply by the square of the residual std dev.  (There's a way
to extract all these values from the fitted model instead of hand-typing
them, but I can't find them at the moment.)

vars <- c(1.0000000, 0.8788793, 1.0744596, 0.9586878)^2 * 2.329213^2

Now create a matrix of variances and covariances.

covs <- outer(vars,vars, function(x,y) sqrt(x)*sqrt(y))

Finally, multiply the correlations by the covariances.  Let's compare
results:

PROC MIXED:

Row Col1 Col2 Col3 Col4
1 5.415 2.716 3.910 2.710
2 2.716 4.184 2.927 3.317
3 3.910 2.927 6.455 4.130
4 2.710 3.317 4.130 4.985

R> round(cors*covs,3)
      [,1]  [,2]  [,3]  [,4]
[1,] 5.425 2.709 3.841 2.715
[2,] 2.709 4.191 2.975 3.314
[3,] 3.841 2.975 6.263 4.133
[4,] 2.715 3.314 4.133 4.986

In my experience, the small differences in the results are typical for
unstructured models.

Hope this helps.

Kevin Wright


On Tue, Jan 24, 2012 at 8:32 PM, Charles Determan Jr <deter088@umn.edu>wrote:

> Greetings,
>
> I have been working on R for some time now and I have begun the endeavor of
> trying to replicate some SAS code in R.  I have scoured the forums but
> haven't been able to find an answer.  I hope one of you could be so kind as
> to enlighten me.
>
> I am attempting to replicate a repeated measures experiment using some
> standard data.  I have posted the SAS code and output directly from a
> publication as well as my attempts in R to replicate it.  My main issue
> comes with the 'unstructured' component.
>
> The 'dental' dataset from 'mixedQF' package,
> equivalent to formixed data in SAS
>
>    distance age Subject    Sex
> 1       26.0   8     M01   Male
> 2       25.0  10     M01   Male
> 3       29.0  12     M01   Male
> 4       31.0  14     M01   Male
> 5       21.5   8     M02   Male
> 6       22.5  10     M02   Male
> 7       23.0  12     M02   Male
> 8       26.5  14     M02   Male
> 9       23.0   8     M03   Male
> 10      22.5  10     M03   Male
> 11      24.0  12     M03   Male
> 12      27.5  14     M03   Male
> 13      25.5   8     M04   Male
> 14      27.5  10     M04   Male
> 15      26.5  12     M04   Male
> 16      27.0  14     M04   Male
> 17      20.0   8     M05   Male
> 18      23.5  10     M05   Male
> 19      22.5  12     M05   Male
> 20      26.0  14     M05   Male
> 21      24.5   8     M06   Male
> 22      25.5  10     M06   Male
> 23      27.0  12     M06   Male
> 24      28.5  14     M06   Male
> 25      22.0   8     M07   Male
> 26      22.0  10     M07   Male
> 27      24.5  12     M07   Male
> 28      26.5  14     M07   Male
> 29      24.0   8     M08   Male
> 30      21.5  10     M08   Male
> 31      24.5  12     M08   Male
> 32      25.5  14     M08   Male
> 33      23.0   8     M09   Male
> 34      20.5  10     M09   Male
> 35      31.0  12     M09   Male
> 36      26.0  14     M09   Male
> 37      27.5   8     M10   Male
> 38      28.0  10     M10   Male
> 39      31.0  12     M10   Male
> 40      31.5  14     M10   Male
> 41      23.0   8     M11   Male
> 42      23.0  10     M11   Male
> 43      23.5  12     M11   Male
> 44      25.0  14     M11   Male
> 45      21.5   8     M12   Male
> 46      23.5  10     M12   Male
> 47      24.0  12     M12   Male
> 48      28.0  14     M12   Male
> 49      17.0   8     M13   Male
> 50      24.5  10     M13   Male
> 51      26.0  12     M13   Male
> 52      29.5  14     M13   Male
> 53      22.5   8     M14   Male
> 54      25.5  10     M14   Male
> 55      25.5  12     M14   Male
> 56      26.0  14     M14   Male
> 57      23.0   8     M15   Male
> 58      24.5  10     M15   Male
> 59      26.0  12     M15   Male
> 60      30.0  14     M15   Male
> 61      22.0   8     M16   Male
> 62      21.5  10     M16   Male
> 63      23.5  12     M16   Male
> 64      25.0  14     M16   Male
> 65      21.0   8     F01 Female
> 66      20.0  10     F01 Female
> 67      21.5  12     F01 Female
> 68      23.0  14     F01 Female
> 69      21.0   8     F02 Female
> 70      21.5  10     F02 Female
> 71      24.0  12     F02 Female
> 72      25.5  14     F02 Female
> 73      20.5   8     F03 Female
> 74      24.0  10     F03 Female
> 75      24.5  12     F03 Female
> 76      26.0  14     F03 Female
> 77      23.5   8     F04 Female
> 78      24.5  10     F04 Female
> 79      25.0  12     F04 Female
> 80      26.5  14     F04 Female
> 81      21.5   8     F05 Female
> 82      23.0  10     F05 Female
> 83      22.5  12     F05 Female
> 84      23.5  14     F05 Female
> 85      20.0   8     F06 Female
> 86      21.0  10     F06 Female
> 87      21.0  12     F06 Female
> 88      22.5  14     F06 Female
> 89      21.5   8     F07 Female
> 90      22.5  10     F07 Female
> 91      23.0  12     F07 Female
> 92      25.0  14     F07 Female
> 93      23.0   8     F08 Female
> 94      23.0  10     F08 Female
> 95      23.5  12     F08 Female
> 96      24.0  14     F08 Female
> 97      20.0   8     F09 Female
> 98      21.0  10     F09 Female
> 99      22.0  12     F09 Female
> 100     21.5  14     F09 Female
> 101     16.5   8     F10 Female
> 102     19.0  10     F10 Female
> 103     19.0  12     F10 Female
> 104     19.5  14     F10 Female
> 105     24.5   8     F11 Female
> 106     25.0  10     F11 Female
> 107     28.0  12     F11 Female
> 108     28.0  14     F11 Female
>
> *Mixed modeling and fixed effect test*
> SAS
> proc mixed data=formixed;
> class gender age person;
> model y = gender|age;
> repeated / type=cs sub=person;
> run;
>
> output of interest to me
>          Tests of Fixed Effects
> Source             NDF   DDF    Type III F    Pr > F
> GENDER           1        25        9.29        0.0054
> AGE                  3        75       35.35       0.0001
> GENDER*AGE   3        75        2.36        0.0781
>
> R (nlme package)
> y=lme(distance~Sex*age, random=(~1|Subject), data=dental)
> anova(y)
>
>            numDF denDF  F-value p-value
> (Intercept)     1    75 4123.156  <.0001
> Sex              1    25    9.292  0.0054
> age               3    75   40.032  <.0001
> Sex:age        3    75    2.362  0.0781
>
> Now this isn't exact but it is extremely close, however when I try to
> replicate the unstructured,
>
> SAS
> proc mixed data=formixed;
> class gender age person;
> model y = gender|age;
> repeated / type=un sub=person;
> run;
>
>             Tests of Fixed Effects
> Source          NDF DDF Type III F Pr > F
> GENDER         1    25     9.29    0.0054
> AGE                3    25    34.45   0.0001
> GENDER*AGE 3    25     2.93    0.0532
>
> R
> either
> y=lme(distance~Sex*age, random=(~1|Subject), corr=corSymm(,~1|Subject),
> data=dental)
> anova(y)
> or
> z=lme(distance~Sex*age, random=(~1|Subject), corr=corSymm(), data=dental)
> anova(z)
>
> gives the output
>
>            numDF denDF  F-value    p-value
> (Intercept)     1    75     4052.028  <.0001
> Sex              1    25       8.462      0.0075
> age               3    75      39.022    <.0001
> Sex:age        3    75       2.868      0.0421
>
> What am I doing wrong to replicate the unstructured linear mixed model from
> SAS?
>
> Regards,
>
> Charles
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>



-- 
Kevin Wright

	[[alternative HTML version deleted]]


