3
$\begingroup$

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

$\endgroup$
  • 2
    $\begingroup$ (+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... $\endgroup$ – Curious Jan 22 '13 at 18:52
  • $\begingroup$ 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) $\endgroup$ – Stéphane Laurent Jan 22 '13 at 21:47
  • $\begingroup$ 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. $\endgroup$ – Stéphane Laurent Jan 22 '13 at 21:51
  • $\begingroup$ 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. $\endgroup$ – Neuhier Jan 23 '13 at 8:23
2
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Neuhier Jan 23 '13 at 8:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.