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. 🙂
Best Answer
If I understand your question correctly, you can specify your model with nested random effects like this:
To specify the covariance structure, e.g. a simple compound symmetry form, try this:
To look at the estimated parameters try:
To get all estimated coefficients try:
To get the p-values then use:
Notice that if you need to specify the covariance structure of the residuals, you'll have to use
nlme
as althoughlme4
(i.e. thelmer
function) is a more advanced package, currently it does not support that feature.