8
$\begingroup$

After doing a model comparison with my mixed lmer model, I have a model with three main effects, no interaction, say signal ~ factor A + factor B + factor C + (1|subj).

Factor C has three levels, so I want to do a post-hoc test to see how the levels differ from each other. I tried two methods:

Method 1: mcp with Tukey (from multcomp package)

summary(glht(myModel, mcp(factorC="Tukey"))

where I get the following result:

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lme4::lmer(formula = signal ~ factorA + factorB + factor C + (1 | 
    subj), data = s)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)  
e1 - e2 == 0   0.8071     0.4681   1.724   0.1984    
e1 - e3 == 0   1.9926     0.4681   4.257   <1e-04 ***
e2 - e3 == 0   1.1855     0.4681   2.533   0.0321 *    ---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Method 2: lsm (from lsmeans package)

summary(glht(myModel, lsm(pairwise ~ factorC)))

giving me the following result:

     Simultaneous Tests for General Linear Hypotheses

Fit: lme4::lmer(formula = signal ~ factorA + factorB + factorC + 
    (1 | ID), data = s)

Linear Hypotheses:
             Estimate Std. Error t value Pr(>|t|)    
e1 - e2 == 0   0.8071     0.4681   1.724    0.198    
e1 - e3 == 0   1.9926     0.4681   4.257   <1e-04 ***
e2 - e3 == 0   1.1855     0.4681   2.533    0.032 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

The results are pretty similar, and I would guess that the lsm-results are more reliable, since lsmeans is explicitly suited for models with unequal observations. I still wonder, though, whether it is acceptable to do so and would appreciate any comment!

$\endgroup$

1 Answer 1

6
$\begingroup$

Are you sure that the results really differ at all? I see only one change in the 4th decimal place (of the p-value of the difference between groups 2 and 3), a relative difference of 0.3%, which could be a numerical difference due to doing equivalent computations in a different sequence.

?lsm says:

It works similarly to ‘mcp’ except with ‘specs’ (and optionally ‘by’ and ‘contr’ arguments) provided as in a call to ‘lsmeans’.

which suggests strongly to me (since lsmeans is generally well-documented) that this is only a different interface to the same functionality: if there were important statistical differences I think they would be mentioned ...

It would be helpful to tell us that lsm comes from the lsmeans package (and glht is from multcomp).

$\endgroup$
4
  • 1
    $\begingroup$ What Ben says is exactly correct. lsm is just another way of generating the same set of linear contrast coefficients, which are then passed to glht. There are other sets of contrasts that may not be available in both. $\endgroup$
    – Russ Lenth
    Commented Oct 24, 2014 at 2:26
  • $\begingroup$ Thanks you two, I conclude from your answers, that there is no difference in calculation between lsm and mcp (as the results suggest, but it is good to be sure!). $\endgroup$ Commented Oct 29, 2014 at 14:21
  • 1
    $\begingroup$ I don't think the methods are the same. The first method uses a "z" statistic, whereas the second one uses a "t" statistic (look also the outputs here). I don't know the details, but I guess the first method is an asymptotic one whereas the second one is more appropriate for small sample sizes. In this example the absence of difference may be due to a large degrees of freedom (a large sample size). cc @evoked_potential $\endgroup$ Commented Mar 31, 2016 at 14:36
  • $\begingroup$ Thanks for pointing out that the underlying distributions differ. In my example, the degrees of freedom were large indeed, so you may be right that this is the cause for the numbers being practically identical. $\endgroup$ Commented Apr 4, 2016 at 13:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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