I answer my own question.

I was disturbed by the syntax of the SAS code for "proc mixed", which uses
the same procedure for mixed or non mixed effect model. Moreover, the type
of variance-covariance structure is indicate by one parameter.

R is more intuitive and supple : there is two different functions (gls and
lme), and two parameters (one for variance structure, one for correlation
structure). I needed too much time to understand this.

My mistake.

david

2008/6/25 David Hajage <dhajage.r@gmail.com>:

> I would like to make clear that the SAS "unstructured" correlation matrix
> in the second model was :
> Estimated R Correlation Matrix for id 1
>
> 1      1.0000      0.8575      0.6984      0.4657      0.3165
> 2      0.8575      1.0000      0.8557      0.5670      0.4039
> 3      0.6984      0.8557      1.0000      0.8717      0.7551
> 4      0.4657      0.5670      0.8717      1.0000      0.9333
> 5      0.3165      0.4039      0.7551      0.9333      1.0000
>
> which seems to be equal to :
>
> > cov2cor(cov(datares))
>           [,1]      [,2]      [,3]      [,4]      [,5]
> [1,] 1.0000000 0.8574841 0.6984365 0.4657120 0.3165061
> [2,] 0.8574841 1.0000000 0.8557224 0.5670489 0.4038805
> [3,] 0.6984365 0.8557224 1.0000000 0.8716642 0.7550925
> [4,] 0.4657120 0.5670489 0.8716642 1.0000000 0.9333166
> [5,] 0.3165061 0.4038805 0.7550925 0.9333166 1.0000000
>
> but not to :
>
> > cov2cor(extract.lme.cov(fm2,rats)[1:5, 1:5]) # corClasses in fm2 :
> corSymm
>           1         2         3         4         5
> 1 1.0000000 0.9153199 0.8424135 0.6222813 0.2125583
> 2 0.9153199 1.0000000 0.9344902 0.6960017 0.3115286
> 3 0.8424135 0.9344902 1.0000000 0.8709635 0.5573176
> 4 0.6222813 0.6960017 0.8709635 1.0000000 0.8218556
> 5 0.2125583 0.3115286 0.5573176 0.8218556 1.0000000
>
>
> 2008/6/25 David Hajage <dhajage.r@gmail.com>:
>
> 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]]

