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.