Solved – Mixed model specification with nlme in R

lme4-nlmemixed model

i'm having a hard time finding the right model for my data. The study is about the effect of a drug (more precise: the effect of the dose of a drug) and the time on the subject (rat). So i got the following variables:

  • Dependent Variable: Relative Change
  • Independent Variables: Dose (Levels: 1, 2, 4, 10) and Time (1,2,3)
  • Further we got the subject ID (1 to 6)

resulting in data like this:

ID Dose Time Change
1  1    1    0.342
1  1    2    0.212
1  1    3    0.121
1  2    1    0.423
1  2    2    0.321
1  2    3    ...
1  4    1    ...
....
6  10   1    0.121
6  10   2    0.221
6  10   3    0.111

As you can see the different doses were applied to all subjects. Three measurements were made for each subject+dose after 1, 2 and 3 minutes.

My idea for a model:

$$ y_{ijt} = \beta_0 + \beta_1 dose_j + \beta_2 time_t + \beta_3 dose_j*time_t + u_{i} + \epsilon_{ijt}$$

with Dose and Time being fixed effects and $u_i$ being the subject-specific random intercept. Of course this model does not respect the repeated measurements. So I would expect the residuals to be correlated (over the three time-points), modeling this with an adequate residual variance structure like the following:

$$ \epsilon_{i} \sim N(0, R_i)$$

with $\epsilon_i$ being a vector of all residuals of subject i and $R_i$ being a block-diagonal matrix with diagonal elements

$$ R_{i,diag} = \begin{pmatrix} \sigma^2 & \sigma_{1,2} & \sigma_{1,3} \\
\sigma_{1,2} & \sigma^2 & \sigma_{2,3} \\
\sigma_{1,3} & \sigma_{2,3}& \sigma^2
\end{pmatrix}$$.

$$ R_i = \begin{pmatrix}
R_{i, diag} & 0 & 0 & 0 \\
0 & R_{i, diag} & 0 & 0 \\
0 & 0 & R_{i, diag} &0 \\
0 & 0 & 0 & R_{i, diag} \\
\end{pmatrix}$$

I would like to know:

  1. Is this an adequate way to model this situation or is there a more simple way to do that (would be great! :-)).
  2. If this is a right way: Is there a way to estimate this using the lme-function in R (nlme-package)? The problem I have: how might I specify this residual variance-covariance-structure? I got the impression that I need to implement my own varFunc, which doesn't seem to be fun.

Concluding remark: At the moment I have no access to "Mixed-Effects Models in S (Pinheiro and Bates)". It will arrive next week. 🙂

Best Answer

If I understand your question correctly, you can specify your model with nested random effects like this:

fit.1 <- lme(Change ~ Dose*Time, random=~1|ID/Dose, data=mydata)

To specify the covariance structure, e.g. a simple compound symmetry form, try this:

fit.2 <- lme(Change ~ Dose*Time, random=~1|ID/Dose, data=mydata, cor=corCompSymm())

To look at the estimated parameters try:

summary(fit.1)

To get all estimated coefficients try:

coef(fit.1)

To get the p-values then use:

anova(fit.1)

Notice that if you need to specify the covariance structure of the residuals, you'll have to use nlme as although lme4 (i.e. the lmer function) is a more advanced package, currently it does not support that feature.

Related Question