[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

```