# Mixed model specification with nlme in R

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. :-)

• (+1) Just a tip: not everything needs to be done using mixed models - you might be fine using "simpler" (more robust??) methods like gls! But I just started with it and have no experience so it's just a tip for you where to investigate... – Curious Jan 22 '13 at 18:52
• I agree with Tomas. You can also use standard multivariate regression, see section 3 of John Fox's appendix socserv.mcmaster.ca/jfox/Books/Companion/appendix/… (but use Type II sum of squares and not Type III) – Stéphane Laurent Jan 22 '13 at 21:47
• And be careful with your design from the practictal point of view. Smoking 2 grams of marijuana has not the same effect according to whether you have already smoked 1 gram yesterday or not. – Stéphane Laurent Jan 22 '13 at 21:51
• Thanks for your comments. I did not really think about using different methods. @Stephane Laurent: There is an adequate period between the runs with different doses. But thanks for your tip. – Neuhier Jan 23 '13 at 8:23

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.

Looks like you want to fit a mixed effects model with Dose-Time interaction with random intercepts for each animale (ID) e.g.

lmer(Change ~ Dose*Time + ( 1 | ID ), data=data)


A model that allows for random slopes (over Time) in addition to random intercepts would look like this:

lmer(Change ~ Dose*Time + ( Time | ID ), data=data)


Also have a look here

• Thanks for your answer. I thought about this models, you are talking about, too. But I got the impression, that these models do not really respect the "nestedness" of my data. I mean in this data I got kind of two repeated measure variables nested in each other: Each dose is used to each rat on the first level. On the second level I got the three measurements for each dose. – Neuhier Jan 23 '13 at 8:33