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.
