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.
Putting my comment into an answer:
The estimates from glht are log odds ratios because they give you differences between logits. You can see this easily, if you write down the maths:
$\ln{\frac{p_{2014}}{1-p_{2014}}} - \ln{\frac{p_{2013}}{1-p_{2013}}}
= \ln{\frac{p_{2014}(1-p_{2013})}{p_{2013}(1-p_{2014})}}$
The estimate of 0.6131 means that the odds (i.e., $\frac{p}{1-p}$) of a chick surviving in 2014 were almost twice as high as in 2013: $\exp(0.6131) = 1.846146$
(Assuming I understand correctly, how you specified the dependent. Possible those are the odds that they died instead.)
Best Answer
By default,
lmertreats the reference level of a categorical predictor as the baseline and estimates parameters for the other levels. So you get some pairwise comparisons in the default output and you can get the others by usingrelevelto define a new reference level and re-fitting the model. This has the advantage of letting you use model comparisons or MCMC to get the p-values, but does not correct for multiple comparisons (though you could apply your own correction afterwards).To use
multcomp, you need to define the contrast matrix. Each row in the contrast matrix represents weights for the effects you get in the default model output, starting with the Intercept. So if you want an effect that's already included in the basic output, you just put a "1" in the position corresponding to that effect. Since the parameter estimates are relative to a common reference level, you can get comparisons between any two other levels by setting the weight of one to "-1" and of the other "1". Here's how that would work for theTime:Dietterms in theChickWeightexample:Caveat emptor: This approach seems to use the normal approximation to get p-values, which is somewhat anti-conservative, and then applies some correction for multiple comparisons. The upshot is that this method gives you easy access to as many pairwise parameter estimates and standard errors as you want, but the p-values may or may not be what you want.
(Thanks to Scott Jackson from r-ling-lang-L for help with this)