[R] gls() error
toby909 at gmail.com
toby909 at gmail.com
Fri May 18 04:02:32 CEST 2007
Hi All
How can I fit a repeated measures analysis using gls? I want to start with a
unstructured correlation structure, as if the the measures at the occations are
not longitudinal (no AR) but plainly multivariate (corSymm).
My data (ignore the prox_pup and gender, occ means occasion):
> head(dta,12)
teacher occ prox_self prox_pup gender
1 1 0 0.76 0.41 1
2 1 1 1.00 1.05 1
3 1 3 0.81 0.91 1
4 2 0 0.96 0.64 0
5 3 0 1.12 1.13 1
6 3 2 1.34 1.35 1
7 4 1 0.35 -0.40 0
8 4 2 0.25 0.27 0
9 4 3 0.54 0.26 0
10 5 0 0.65 1.02 1
11 5 1 0.68 0.87 1
12 5 2 1.04 0.98 1
x=factor(dta$occ)
Following gives me an error message:
gls(prox_pup~x-1, dta, corSymm(, ~x-1|teacher))
> gls(prox_pup~x-1, dta, corSymm(, ~x-1|teacher))
Error in Initialize.corSymm(X[[1]], ...) :
Covariate must have unique values within groups for corSymm objects
In addition: There were 50 or more warnings (use warnings() to see the first 50)
I checked that the covariate, occ, has unique values within each of the teachers.
Aside, lme actually gives me what I want, except that the residual variance is
not where I want to have it. I want the residuals being part of the covariance
matrix to be estimated rather than in the 'level one' residual, ie the residuals
included on the diagonal of "G" displayed below.
> a4 = lme(prox_pup~x-1, dta, ~x-1|teacher)
Linear mixed-effects model fit by REML
Data: dta
Log-restricted-likelihood: -53.91059
Fixed: prox_pup ~ x - 1
x0 x1 x2 x3
0.5739072 0.7169963 0.6503699 0.6567064
Random effects:
Formula: ~x - 1 | teacher
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
x0 0.5424187 x0 x1 x2
x1 0.4326164 0.739
x2 0.3343281 0.611 0.681
x3 0.3638630 0.569 0.752 0.900
Residual 0.0929820
Number of Observations: 153
Number of Groups: 51
> G = lapply(pdMatrix(a4$modelStruct$reStruct), "*", a4$sigma^2)
$teacher
x0 x1 x2 x3
x0 0.2942180 0.17330375 0.11089028 0.1123597
x1 0.1733037 0.18715693 0.09847681 0.1183089
x2 0.1108903 0.09847681 0.11177526 0.1094374
x3 0.1123597 0.11830892 0.10943738 0.1323963
Thanks for your help on this.
Toby
More information about the R-help
mailing list