[R] unstructured correlation matrix in lme

From: David Hajage <dhajage.r_at_gmail.com>
Date: Wed, 25 Jun 2008 18:01:56 +0200


Hello R users,

I'm student and I'm actually having a lecture introducing repeated mesures analysis. Unfortunately, all examples use SAS system...

I'm working with lme function (package "nlme"), and I'm using extract.lme.cov (package "mgcv") to extract covariance structure of models.

One example is about weight evolution of rats (3 treatment groups) :

> summary(rats)

       id      trt          y         s
 1      :  5   1:50   Min.   : 46.0   0:27
 10     :  5   2:50   1st Qu.: 71.0   1:27
 11     :  5   3:35   Median :100.0   2:27
 12     :  5          Mean   :100.8   3:27
 13     :  5          3rd Qu.:122.5   4:27
 14     :  5          Max.   :189.0

 (Other):105

id : subjects
trt : treatment groups (1, 2 or 3)
y : weights
s : weeks (Week 0 to Week 4)

There are 5 mesearements of weight (on week 0, 1, 2, 3 and 4) for each rat.

I first fit a model with a compound symmetry correlation structure :

> library(nlme)
> library(mgcv)
> fm1 <- lme(y ~ as.factor(s) * as.factor(trt), data = rats, random = ~
1|id, correlation = corCompSymm())
> summary(fm1)
...
Fixed effects: y ~ as.factor(s) * as.factor(trt)
                                  Value Std.Error DF  t-value p-value
(Intercept)                    54.00000  3.469019 96 15.56636  0.0000
as.factor(s)1                  24.50000  3.125974 96  7.83756  0.0000
as.factor(s)2                  52.00000  3.125974 96 16.63481  0.0000
as.factor(s)3                  76.10000  3.125974 96 24.34441  0.0000
as.factor(s)4                 106.60000  3.125974 96 34.10137  0.0000
as.factor(trt)2                 0.70000  4.905934 24  0.14268  0.8877
as.factor(trt)3                 1.57143  5.406076 24  0.29068  0.7738
as.factor(s)1:as.factor(trt)2  -2.90000  4.420795 96 -0.65599  0.5134
as.factor(s)2:as.factor(trt)2 -10.90000  4.420795 96 -2.46562  0.0155
as.factor(s)3:as.factor(trt)2 -22.60000  4.420795 96 -5.11220  0.0000
as.factor(s)4:as.factor(trt)2 -37.30000  4.420795 96 -8.43740  0.0000
as.factor(s)1:as.factor(trt)3  -4.21429  4.871479 96 -0.86509  0.3891
as.factor(s)2:as.factor(trt)3  -2.71429  4.871479 96 -0.55718  0.5787
as.factor(s)3:as.factor(trt)3   1.04286  4.871479 96  0.21407  0.8309
as.factor(s)4:as.factor(trt)3 0.68571 4.871479 96 0.14076 0.8884 ...

> extract.lme.cov(fm1,rats)[1:5, 1:5]

          1 2 3 4 5

1 120.34095  71.48238  71.48238  71.48238  71.48238
2  71.48238 120.34095  71.48238  71.48238  71.48238
3  71.48238  71.48238 120.34095  71.48238  71.48238
4  71.48238  71.48238  71.48238 120.34095  71.48238
5  71.48238  71.48238  71.48238  71.48238 120.34095

As you can see, I obtain exactly the same solutions for fixed parameters and estimated covariance matrix in SAS :
proc mixed data = rat ;
class id trt s ;
format trt ftrt. s fs. ;
model y = trt s trt * s / solution ddfm = satterth ; repeated s / type = cs subject = id r rcorr; run ;

                                                   Erreur
    Effet        trt       s       Estimation    standard      DF    Valeur
du test t Pr > |t|
    Intercept                         54.0000      3.4690
49.8               15.57      <.0001
    trt          A.trt3                1.5714      5.4061
49.8                0.29      0.7725
    trt          B.trt2                0.7000      4.9059
49.8                0.14      0.8871
    trt          C.trt1                     0           .

. . .
s A.S4 106.60 3.1260 96 34.10 <.0001 s B.S3 76.1000 3.1260 96 24.34 <.0001 s C.S2 52.0000 3.1260 96 16.63 <.0001 s D.S1 24.5000 3.1260 96 7.84 <.0001 s E.S0 0 .
. . .
trt*s A.trt3 A.S4 0.6857 4.8715 96 0.14 0.8884 trt*s A.trt3 B.S3 1.0429 4.8715 96 0.21 0.8309 trt*s A.trt3 C.S2 -2.7143 4.8715 96 -0.56 0.5787 trt*s A.trt3 D.S1 -4.2143 4.8715 96 -0.87 0.3891 trt*s A.trt3 E.S0 0 .
. . .
trt*s B.trt2 A.S4 -37.3000 4.4208 96 -8.44 <.0001 trt*s B.trt2 B.S3 -22.6000 4.4208 96 -5.11 <.0001 trt*s B.trt2 C.S2 -10.9000 4.4208 96 -2.47 0.0155 trt*s B.trt2 D.S1 -2.9000 4.4208 96 -0.66 0.5134 trt*s B.trt2 E.S0 0 .
. . .
trt*s C.trt1 A.S4 0 .
. . .
trt*s C.trt1 B.S3 0 .
. . .
trt*s C.trt1 C.S2 0 .
. . .
trt*s C.trt1 D.S1 0 .
. . .
trt*s C.trt1 E.S0 0 .
. . .

Estimated R Matrix for id 1

1      120.34     71.4824     71.4824     71.4824     71.4824
2     71.4824      120.34     71.4824     71.4824     71.4824
3     71.4824     71.4824      120.34     71.4824     71.4824
4     71.4824     71.4824     71.4824      120.34     71.4824
5     71.4824     71.4824     71.4824     71.4824      120.34

My problem is that I can't find how to obtain an "unstructered" correlation structure with lme function (SAS option : type = un), ie to use this covariance matrix :
> datares <- t(unstack(data.frame(fm1$residuals[,1], rats$id))) > cov(datares)

         [,1] [,2] [,3] [,4] [,5]

[1,] 19.91593 30.47967  29.15275  25.79780  21.44505
[2,] 30.47967 63.44066  63.74835  56.06209  48.84066
[3,] 29.15275 63.74835  87.47912 101.19670 107.22527
[4,] 25.79780 56.06209 101.19670 154.07418 175.88901
[5,] 21.44505 48.84066 107.22527 175.88901 230.50989

My first thought was that it was corClasses "corSymm", but I don't have the same results in R and SAS :

With R :
> fm2 <- lme(y ~ as.factor(s) * as.factor(trt), data = rats, random = ~ 1|id, correlation = corSymm())
 > summary(fm2)
...

                                  Value Std.Error DF   t-value p-value
(Intercept)                    54.00000  3.767777 96 14.332057  0.0000
as.factor(s)1                  24.50000  1.550568 96 15.800659  0.0000
 as.factor(s)2                  52.00000  2.115240 96 24.583496  0.0000
as.factor(s)3                  76.10000  3.274798 96 23.238076  0.0000
as.factor(s)4                 106.60000  4.728348 96 22.544871  0.0000
as.factor(trt)2                 0.70000  5.328442 24  0.131370  0.8966
as.factor(trt)3                 1.57143  5.871657 24  0.267629  0.7913
as.factor(s)1:as.factor(trt)2  -2.90000  2.192835 96 -1.322489  0.1891
as.factor(s)2:as.factor(trt)2 -10.90000  2.991401 96 -3.643777  0.0004
as.factor(s)3:as.factor(trt)2 -22.60000  4.631263 96 -4.879878  0.0000
as.factor(s)4:as.factor(trt)2 -37.30000  6.686894 96 -5.578075  0.0000
as.factor(s)1:as.factor(trt)3  -4.21429  2.416386 96 -1.744045  0.0844
as.factor(s)2:as.factor(trt)3  -2.71429  3.296364 96 -0.823418  0.4123
as.factor(s)3:as.factor(trt)3   1.04286  5.103404 96  0.204345  0.8385
as.factor(s)4:as.factor(trt)3 0.68571 7.368598 96 0.093059 0.9261 ...

> extract.lme.cov(fm2, rats)[1:5, 1:5]

          1 2 3 4 5

1 141.96147 129.94016 119.59027  88.33997  30.17509
2 129.94016 141.96147 132.66161  98.80543  44.22506
3 119.59027 132.66161 141.96147 123.64326  79.11763
4  88.33997  98.80543 123.64326 141.96147 116.67183
5  30.17509  44.22506  79.11763 116.67183 141.96147

With SAS :
proc mixed data = rat ;
class id trt s ;
format trt ftrt. s fs. ;
model y = trt s trt * s / solution ddfm = satterth ; repeated s / type = un subject = id r rcorr; run ;

                                                   Erreur
     Effet        trt       s       Estimation    standard      DF    Valeur
du test t Pr > |t|
    Intercept                         54.0000      1.4689
24               36.76      <.0001
    trt          A.trt3                1.5714      2.2891
24                0.69      0.4990
    trt          B.trt2                0.7000      2.0773
24                0.34      0.7391
    trt          C.trt1                     0           .

. . .
s A.S4 106.60 4.7416 24 22.48 <.0001 s B.S3 76.1000 3.6413 24 20.90 <.0001 s C.S2 52.0000 2.3061 24 22.55 <.0001 s D.S1 24.5000 1.5577 24 15.73 <.0001 s E.S0 0 .
. . .
trt*s A.trt3 A.S4 0.6857 7.3893 24 0.09 0.9268 trt*s A.trt3 B.S3 1.0429 5.6746 24 0.18 0.8557 trt*s A.trt3 C.S2 -2.7143 3.5938 24 -0.76 0.4574 trt*s A.trt3 D.S1 -4.2143 2.4275 24 -1.74 0.0954 trt*s A.trt3 E.S0 0 .
. . .
trt*s B.trt2 A.S4 -37.3000 6.7057 24 -5.56 <.0001 trt*s B.trt2 B.S3 -22.6000 5.1496 24 -4.39 0.0002 trt*s B.trt2 C.S2 -10.9000 3.2613
24 -3.34 0.0027 trt*s B.trt2 D.S1 -2.9000 2.2029 24 -1.32 0.2005 trt*s B.trt2 E.S0 0 .
. . .
trt*s C.trt1 A.S4 0 .
. . .
trt*s C.trt1 B.S3 0 .
. . .
trt*s C.trt1 C.S2 0 .
. . .
trt*s C.trt1 D.S1 0 .
. . .
trt*s C.trt1 E.S0 0 .
. . .

Estimated R Matrix for id 1

1     21.5756     33.0196     31.5821     27.9476     23.2321
2     33.0196     68.7274     69.0607     60.7339     52.9107
3     31.5821     69.0607     94.7690      109.63      116.16
4     27.9476     60.7339      109.63      166.91      190.55
 5     23.2321     52.9107      116.16      190.55      249.72

I know that R is not SAS, and that the word "unstructured" is perhaps not appropriate and a bit confusing. But I would like to know if there is an equivalent code for this in R ?

Thank you very much.

        [[alternative HTML version deleted]]



R-help_at_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. Received on Wed 25 Jun 2008 - 16:10:10 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 25 Jun 2008 - 17:31:13 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.

list of date sections of archive