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]]
______________________________________________
[email protected] 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.
