# Post-hoc tests in mixed model using lsmeans and glht

I have a mixed model with two fixed factors (Insecte and temperature) with two random factors (bloc and date) nested in the temperature treatment. I also include interactions between the fixed factor Insecte and the random factors. The response variable is aphid density. This give the following model:

m2.nlmer = lmer(aphid_density~ Insecte*amp_freq
+(1|bloc:(date:amp_freq))+(1|date:amp_freq)+(1|Insecte:(bloc:(date:amp_freq)))+(1|Insecte:(date:amp_freq)),
na.action=na.exclude,
data = fr2)


where amp_freq are the temperature treatments. I detected significant interaction between Insecte and amp_freq and want to analyse Insecte effect for each temperature levels as follows:

> lsmeans(m2.nlmer, pairwise ~ Insecte|amp_freq, adjustSigma = TRUE, adjust = "tukey")

$lsmeans amp_freq = 23-const: Insecte lsmean SE df lower.CL upper.CL COCCINELLE 568.17947 160.6389 5.39 164.0381 972.3208 PUCERON 1783.53300 160.0224 5.32 1379.5075 2187.5585 amp_freq = 30-jour: Insecte lsmean SE df lower.CL upper.CL COCCINELLE 570.64910 161.6708 5.48 165.8695 975.4287 PUCERON 1314.53794 160.0629 5.32 910.4669 1718.6089 amp_freq = 30-semaine: Insecte lsmean SE df lower.CL upper.CL COCCINELLE 738.12107 163.2045 5.74 334.3272 1141.9149 PUCERON 1867.47572 160.0204 5.32 1463.4669 2271.4845 amp_freq = 40-jour: Insecte lsmean SE df lower.CL upper.CL COCCINELLE 84.87312 162.5806 5.64 -319.1406 488.8869 PUCERON 256.55122 159.6785 5.27 -147.5851 660.6875 amp_freq = 40-semaine: Insecte lsmean SE df lower.CL upper.CL COCCINELLE 296.63947 159.6828 5.27 -107.5317 700.8107 PUCERON 837.65225 160.0241 5.32 433.6129 1241.6916 Confidence level used: 0.95$contrasts
amp_freq = 23-const:
contrast               estimate       SE   df t.ratio p.value
COCCINELLE - PUCERON -1215.3535 61.59491 4.44 -19.731  <.0001

amp_freq = 30-jour:
contrast               estimate       SE   df t.ratio p.value
COCCINELLE - PUCERON  -743.8888 64.11036 4.99 -11.603  0.0001

amp_freq = 30-semaine:
contrast               estimate       SE   df t.ratio p.value
COCCINELLE - PUCERON -1129.3546 68.03274 6.60 -16.600  <.0001

amp_freq = 40-jour:
contrast               estimate       SE   df t.ratio p.value
COCCINELLE - PUCERON  -171.6781 65.61395 5.69  -2.616  0.0418

amp_freq = 40-semaine:
contrast               estimate       SE   df t.ratio p.value
COCCINELLE - PUCERON  -541.0128 58.89646 3.77  -9.186  0.0010


However, I get a different results if I split the data set by temperature level and run again the model. For instance:

m40_jour = lmer(aphid_density~ Insecte+(1|bloc:(date:amp_freq))+(1|date:amp_freq)+(1|Insecte:(bloc:(date:amp_freq)))+(1|Insecte:(date:amp_freq)),
na.action=na.exclude,
data = fr2[which(fr2$amp_freq == "40-jour"),]) lsmeans(m40_jour, pairwise ~ Insecte, adjustSigma = TRUE, adjust = "tukey")$lsmeans
Insecte       lsmean       SE   df  lower.CL upper.CL
COCCINELLE  81.90721 29.60546 2.37 -28.08425 191.8987
PUCERON    258.15338 25.84196 1.58 113.00870 403.2981

Confidence level used: 0.95

$contrasts contrast estimate SE df t.ratio p.value COCCINELLE - PUCERON -176.2462 34.73064 1.02 -5.075 0.1202  Finally, I get a different result if I use the glht function with mcp:  summary(glht(m40_jour, mcp(Insecte="Tukey"))) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: Tukey Contrasts Fit: lmer(formula = aphid_density ~ Insecte + (1 | bloc:(date:amp_freq)) + (1 | date:amp_freq) + (1 | Insecte:(bloc:(date:amp_freq))) + (1 | Insecte:(date:amp_freq)), data = fr2[which(fr2$amp_freq ==
"40-jour"), ], na.action = na.exclude)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
PUCERON - COCCINELLE == 0   176.25      33.43   5.273 1.34e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


So I do not understand why (1) I get different results between lsmeans with the interaction between fixed factor and lsmeans with the split dataset and model m40-day and (2) lsmeans and glht produce different results?

What is the correct/more robust approach?

• For the 2nd part of your question, the main reason is degrees of freedom. The test from lsmeans uses 1.02 df and the one from glht uses infinite df (asymptotic z test). That makes a huge difference in the P value. On a lesser note, The estimates match but the SEs differ slightly. That's because lsmeans uses an adjusted covariance matrix from the pbkrtest package. If you do lsm.options(disable.pbkrtest = TRUE), the results will match. – rvl Jan 16 '16 at 18:10
• Thanks for your reply. I tried with lsm.options(disable.pbkrtest=TRUE) and indeed get exactly the same estimate, SEs, z and p values than with glht. So my question is then: is it more accurate to use an asymptotic z test or a t-test with a given df? – Arno Jan 16 '16 at 18:24
• Assumping you have infinite information about a variance instead of only 1 df worth is probably not a good idea! I expanded on my comment, as well as the other part of the question, in the answer I just posted. – rvl Jan 16 '16 at 18:36

Which is better? Well, that's not cut-and-dried. How well do you believe the underlying assumptions for the fuill model? They include the idea that the variances of each random effect is the same for each level of amp_freq. If that's believable, then the full-data model is preferable because you can pool the information about those effects from the whole dataset, making the tests more powerful and the confidence intervals shorter. But if you doubt those assumptions, then separate models may be better, or a more sophisticated model that allows for different variances across amp_freq.
For the 2nd part of your question, the main reason is degrees of freedom. The test from lsmeans uses 1.02 df and the one from glht uses infinite df (asymptotic z test). That makes a huge difference in the P value. Also, while the estimates match, the SEs differ slightly. That's because lsmeans uses an adjusted covariance matrix from the pbkrtest package. That package is also used to get the degrees of freedom. So, if you do lsm.options(disable.pbkrtest = TRUE), the lsmeans results will match those from glht. But I suggest that isn't a good idea, because that 1.02 df is a warning that you really don't have a lot of information about your estimate of that difference.