# [R] unstructured correlation matrix in lme

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

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.