5

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

Neuhier
  • 133
  • 1
  • 1
  • 6
  • 2
    (+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... – Tomas 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 http://socserv.mcmaster.ca/jfox/Books/Companion/appendix/Appendix-Multivariate-Linear-Models.pdf (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 08:23

2 Answers2

4

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

Community
  • 1
user12719
  • 1,009
  • 1
  • 8
  • 10
  • 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 08:33
4

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.

AlefSin
  • 945
  • 7
  • 13