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 <[EMAIL PROTECTED]>:
> 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 <[EMAIL PROTECTED]>:
>
> 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@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.
