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:


*

*Is this an adequate way to model this situation or is there a more simple way to do that (would be great! :-)).

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