32
$\begingroup$

I am trying to analyse some data using a mixed effect model. The data I collected represent the weight of some young animals of different genotype over time.

I am using the approach proposed here: https://gribblelab.wordpress.com/2009/03/09/repeated-measures-anova-using-r/

In particular I'm using solution #2

So I have something like

require(nlme)
model <- lme(weight ~ time * Genotype, random = ~1|Animal/time, 
         data=weights)    
av <- anova(model)

Now, I would like to have some multiple comparisons. Using multcomp I can do:

require(multcomp)
comp.geno <- glht(model, linfct=mcp(Genotype="Tukey"))
print(summary(comp.geno))

And, of course, I could do the same with time.

I have two questions:

  1. How do I use mcp to see the interaction between Time and Genotype?
  2. When I run glht I get this warning:

    covariate interactions found -- default contrast might be inappropriate

    What does it mean? Can I safely ignore it? Or what should I do to avoid it?

EDIT: I found this PDF that says:

Because it is impossible to determine the parameters of interest automatically in this case, mcp() in multcomp will by default generate comparisons for the main effects only, ignoring covariates and interactions. Since version 1.1-2, one can specify to average over interaction terms and covariates using arguments interaction_average = TRUE and covariate_average = TRUE respectively, whereas versions older than 1.0-0 automatically averaged over interaction terms. We suggest to the users, however, that they write out, manually, the set of contrasts they want. One should do this whenever there is doubt about what the default contrasts measure, which typically happens in models with higher order interaction terms. We refer to Hsu (1996), Chapter~7, and Searle (1971), Chapter~7.3, for further discussions and examples on this issue.

I do not have access to those books, but maybe someone here has?

$\endgroup$
  • $\begingroup$ Have a look at InvivoStat (invivostat.co.uk), it should do the type of analysis that you are looking for $\endgroup$ – user2988 Jan 31 '11 at 12:19
29
$\begingroup$

If time and Genotype are both categorical predictors as they appear to be, and you are interested in comparing all time/Genotype pairs to each other, then you can just create one interaction variable, and use Tukey contrasts on it:

weights$TimeGeno <- interaction(weigths$Time, weights$Geno)
model <- lme(weight ~ TimeGeno, random = ~1|Animal/time, data=weights) 
comp.timegeno <- glht(model, linfct=mcp(TimeGeno="Tukey")) 

If you are interested in other contrasts, then you can use the fact that the linfct argument can take a matrix of coefficients for the contrasts - this way you can set up exactly the comparisons you want.

EDIT

There appears some concern in the comments that the model fitted with the TimeGeno predictor is different from the original model fitted with the Time * Genotype predictor. This is not the case, the models are equivalent. The only difference is in the parametrization of the fixed effects, which is set up to make it easier to use the glht function.

I have used one of the built-in datasets (it has Diet instead of Genotype) to demonstrate that the two approaches have the same likelihood, predicted values, etc:

> # extract a subset of a built-in dataset for the example
> data(BodyWeight)
> ex <- as.data.frame(subset(BodyWeight, Time %in% c(1, 22, 44)))
> ex$Time <- factor(ex$Time)
> 
> #create interaction variable
> ex$TimeDiet <- interaction(ex$Time, ex$Diet)
    > 
    > model1 <- lme(weight ~ Time * Diet, random = ~1|Rat/Time,  data=ex)    
    > model2 <- lme(weight ~ TimeDiet, random = ~1|Rat/Time, data=ex)    
    > 
    > # the degrees of freedom, AIC, BIC, log-likelihood are all the same 
    > anova(model1, model2)
           Model df      AIC      BIC    logLik
    model1     1 12 367.4266 387.3893 -171.7133
    model2     2 12 367.4266 387.3893 -171.7133
    Warning message:
    In anova.lme(model1, model2) :
      fitted objects with different fixed effects. REML comparisons are not meaningful.
    > 
    > # the second model collapses the main and interaction effects of the first model
    > anova(model1)
                numDF denDF   F-value p-value
    (Intercept)     1    26 1719.5059  <.0001
    Time            2    26   28.9986  <.0001
    Diet            2    13   85.3659  <.0001
    Time:Diet       4    26    1.7610  0.1671
    > anova(model2)
                numDF denDF   F-value p-value
    (Intercept)     1    24 1719.5059  <.0001
    TimeDiet        8    24   29.4716  <.0001
    > 
    > # they give the same predicted values
    > newdata <- expand.grid(Time=levels(ex$Time), Diet=levels(ex$Diet))
    > newdata$TimeDiet <- interaction(newdata$Time, newdata$Diet)
> newdata$pred1 <- predict(model1, newdata=newdata, level=0)
    > newdata$pred2 <- predict(model2, newdata=newdata, level=0)
> newdata
  Time Diet TimeDiet   pred1   pred2
1    1    1      1.1 250.625 250.625
2   22    1     22.1 261.875 261.875
3   44    1     44.1 267.250 267.250
4    1    2      1.2 453.750 453.750
5   22    2     22.2 475.000 475.000
6   44    2     44.2 488.750 488.750
7    1    3      1.3 508.750 508.750
8   22    3     22.3 518.250 518.250
9   44    3     44.3 530.000 530.000

The only difference is that what hypotheses are easy to test. For example, in the first model it is easy to test whether the two predictors interact, in the second model there is no explicit test for this. On the other hand, the joint effect of the two predictors is easy to test in the second model, but not the first one. The other hypotheses are testable, it is just more work to set those up.

| cite | improve this answer | |
$\endgroup$
  • 3
    $\begingroup$ glht uses the degrees of freedom given in the lme model. I'm not sure these degrees of freedom are appropriate... ? $\endgroup$ – Stéphane Laurent Feb 4 '13 at 13:12
  • 2
    $\begingroup$ I am also curious how this is best done. This approach, however, is giving effects from a different model - one that essentially only tests for an interaction. The second model does not include the two potential main effects at all. This does not seem to be an appropriate method for checking the effects in the first model. $\endgroup$ – Marcus Morrisey Oct 21 '14 at 4:27
  • $\begingroup$ @Aniko, I was thinking of combining 2 categorical variable into one as you just did, but I hesited because just one of the variables is within subject, the other is between. Can you confirm that this doesn't mather? I noticed that in the example you keep Animal/time which now is not one of the factors. Does lme really understand this? $\endgroup$ – toto_tico Oct 27 '15 at 21:54
  • $\begingroup$ @toto_tico, I have edited the response to show that the second model is equivalent to the first one. $\endgroup$ – Aniko Oct 28 '15 at 21:21
  • 3
    $\begingroup$ @toto_tico, I gave you a reproducible example. Why don't you try all.equal(resid(model1), resid(model2)) and see that they are the same before guessing otherwise? Only the fixed effects parametrization is different. TimeDiet is not a pure interaction term, and it is not equivalent to Time:Diet, but rather to Time + Diet + Time:Diet. $\endgroup$ – Aniko Oct 29 '15 at 21:52

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.