[R] Help with lme and correlated residuals
Spencer Graves
spencer.graves at pdf.com
Sun Mar 5 01:40:43 CET 2006
I was able to replicate your error; without the data you provided,
it would have been much more difficult to comment.
Regarding "correlation" and "weights", have you consulted Pinheiro
and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer)? For
me this is an essential reference for "lme" and is also an excellent
book for these kinds of models more generally, whether you use S-Plus or
R or some other software. The standard distribution of R includes in
"~library\nlme\scripts" text files of R commands used in that book.
This makes it easy to study that book, because you can open these script
files and walk through them line by line while you read the book, trying
modifications, etc.
I won't comment more about "correlation" and "weights", but permit me
to ask some other questions about things I saw in your data and in your
call the "lme": First, I noticed that your model treats "subject" as a
numeric variable. Do you really intend to tell the algorithm that
subject 1111 is half way between subjects 1110 and 1112, and that the
difference between either and subject 1003 is 50 times the difference
between 1110 and 1112? I assume you want to make "subject" a factor.
Your data shows that subject is nested within Trial.
Second, I have another question about your "fixed" model: Are you
really trying to say that you want to force to 0 the b0 and b2 terms in
the more common model:
outcome = b0 + b1*endpoint + b2*trt + b12*endpoint*trt + error
Unless you really do want to foce b0 and b2 to 0, then I suggest your
consider a fixed model like the following:
outcome ~ (endpoint+trt)^2
This model will ask mle to estimate b0, b1, b2, and b12. Your
"random" model says you have the same thing, but "Trial" is a random
effect and b1 and b12 (or b0, b1, b2, and b12 if you use the alternative
I just suggested) assume different values for each trial. With this
more common model, I got the following:
> (bmr0 <- lme(outcome~(endpoint+trt)^2,
+ random=~1|Trial/subject,data=datt))
Linear mixed-effects model fit by REML
Data: datt
Log-restricted-likelihood: -55.83479
Fixed: outcome ~ (endpoint + trt)^2
(Intercept) endpoint trt endpoint:trt
-2.675 0.625 -3.425 -0.125
Random effects:
Formula: ~1 | Trial
(Intercept)
StdDev: 0.0003620097
Formula: ~1 | subject %in% Trial
(Intercept) Residual
StdDev: 7.398842 6.158618
When I tried to add "endpoint" to the random model, I got the
following:
> (bmr1 <- lme(outcome~(endpoint+trt)^2,
+ random=~endpoint|Trial/subject,data=datt))
Error in lme.formula(outcome ~ (endpoint + trt)^2, random = ~endpoint | :
iteration limit reached without convergence (9)
I got the same result when I tied to add "trt" to the random model.
I then tried to set b0 and b2 to 0 and got the following:
> (bmr0. <- lme(outcome~-1+endpoint+endpoint:trt,
+ random=~1|Trial/subject,data=datt))
Linear mixed-effects model fit by REML
Data: datt
Log-restricted-likelihood: -61.0509
Fixed: outcome ~ -1 + endpoint + endpoint:trt
endpoint endpoint:trt
0.625 -0.125
Random effects:
Formula: ~1 | Trial
(Intercept)
StdDev: 0.0004386636
Formula: ~1 | subject %in% Trial
(Intercept) Residual
StdDev: 7.699722 6.158618
Number of Observations: 18
Number of Groups:
Trial subject %in% Trial
3 9
> (bmr1. <- lme(outcome~-1+endpoint+endpoint:trt,
+ random=
+ ~-1+endpoint+endpoint:trt|Trial/subject,
+ data=datt))
Linear mixed-effects model fit by REML
Data: datt
Log-restricted-likelihood: -62.79024
Fixed: outcome ~ -1 + endpoint + endpoint:trt
endpoint endpoint:trt
0.625 -0.125
Random effects:
Formula: ~-1 + endpoint + endpoint:trt | Trial
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
endpoint 0.0001850424 endpnt
endpoint:trt 0.0002087138 0
Formula: ~-1 + endpoint + endpoint:trt | subject %in% Trial
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
endpoint 1.858120e-04 endpnt
endpoint:trt 1.861459e-04 0
Residual 1.022864e+01
Number of Observations: 18
Number of Groups:
Trial subject %in% Trial
3 9
hope this helps.
spencer graves
Pryseley Assam wrote:
>
> Dear R - Users
> I have some problems fitting a linear mixed effects model using the lme function (nlme library). A sample data is as shown at the bottom of this mail. I fit my linear mixed model
> using the following R code:
>
> bmr <-lme (outcome~ -1 + as.factor(endpoint)+ as.factor(endpoint):trt, data=datt,
> random=~-1 + as.factor(endpoint) + as.factor(endpoint):trt|as.factor(Trial),
> correlation = corSymm(form=~subject|as.factor(endpoint)),
> weights=varIdent (form=~subject|endpoint))
>
> With this code, i want to obtain random effects for each Trial. ALso my residuals are correlated, and the correlated residuals are heteroscedastic over endpoint. That is, each endpoint has a different constant variance.
>
> However, I recieve the following error message using the code above:
>
> Error in lme.formula(outcome ~ -1 + as.factor(endpoint) + as.factor(endpoint):trt, :
> Incompatible formulas for groups in "random" and "correlation"
>
> May someone kindly inform me of how to correct my code.
>
> Does this imply that the grouping factor in the correlation formula must be the same grouping factor in the random formula ? If so, is there a way of getting pass this restriction ?
>
> In essence, I wish to know how to fit a linear mixed model with correlated residuals (based on a variable in my model) and to obtain the correlation matrix.
>
> Best regards
> Pryseley
>
> Sample data:
>
> RowNames Trial subject VISUAL0 TRT VISUAL24 VISUAL52 TREAT outcome endpoint trt
> 4 1 1003 65 4 65 55 2 0 1 1
> 8 1 1007 67 1 64 68 2 -3 1 -1
> 12 2 1110 59 4 53 42 2 -6 1 1
> 14 2 1111 64 1 72 65 2 8 1 -1
> 16 2 1112 39 1 37 37 2 -2 1 -1
> 18 2 1115 59 4 54 58 2 -5 1 1
> 24 3 1806 46 4 27 24 2 -19 1 1
> 26 3 1813 31 4 33 48 2 2 1 1
> 28 3 1815 64 1 67 64 2 3 1 -1
>
> 4 1 1003 65 4 65 55 2 -10 -1 1
> 8 1 1007 67 1 64 68 2 1 -1 -1
> 12 2 1110 59 4 53 42 2 -17 -1 1
> 14 2 1111 64 1 72 65 2 1 -1 -1
> 16 2 1112 39 1 37 37 2 -2 -1 -1
> 18 2 1115 59 4 54 58 2 -1 -1 1
> 24 3 1806 46 4 27 24 2 -22 -1 1
> 26 3 1813 31 4 33 48 2 17 -1 1
> 28 3 1815 64 1 67 64 2 0 -1 -1
>
>
> Thanks
>
>
>
>
>
>
>
>
> ---------------------------------
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
More information about the R-help
mailing list