0
$\begingroup$

I have a model of the following structure:

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

and find a significant main effect for factorA.

I then compute a post hoc test using lsmeans::lsm

summary(glht(finalModel, lsm(pairwise ~ factorA)))

and get the following result:

Note: df set to 132

     Simultaneous Tests for General Linear Hypotheses

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

Linear Hypotheses:
               Estimate Std. Error t value Pr(>|t|)   
f1 - f2 == 0  -1.4293     0.6743  -2.120  0.21647   
f1 - f3 == 0  -2.3160     0.6743  -3.434  0.00694 **
f1 - f4 == 0  -1.9429     0.6750  -2.878  0.03661 * 
f1 - f5 == 0  -1.0871     0.7796  -1.394  0.63049   
f2 - f3 == 0  -0.8867     0.6373  -1.391  0.63240   
f2 - f4 == 0  -0.5136     0.6668  -0.770  0.93822   
f2 - f5 == 0   0.3423     0.7742   0.442  0.99192   
f3 - f4 == 0   0.3731     0.6668   0.560  0.98039   
f3 - f5 == 0   1.2289     0.7742   1.587  0.50619   
f4 - f5 == 0   0.8558     0.7327   1.168  0.76805

I actually just wanted to check, whether the statistics were really computed with df=132, so that t(132)=-3.434, p=0.007.

So I computed: 2*pt(3.434, 132, lower=F), but I get a much smaller p-value instead.

Any suggestions why this goes wrong?

Is it, because the df in a mixed model are not so straightforward as in simple models or even not appropriate, that the pt()-formula fails?

$\endgroup$

1 Answer 1

1
$\begingroup$

There are at least 2 issues here.

First, the df of 132 presumably represents the total degrees of freedom left in the model after parameter estimation. That is not the df for each of the pairwise comparisons among levels of factorA. The df for each of the pairwise comparisons should be based on the number of observations for the 2 factor levels being compared, as in a standard t-test.

Second, the p-values reported by the summary function are adjusted for the multiple post-hoc comparisons. If you did not specify a particular type of adjustment then a default was chosen. See the documentation for the lsmeans package or the related glht methods in the multcomp package. Thus the reported p-values can't be compared directly against the values for a t distribution in the way that you tried.

There may be additional issues in determining df in a mixed model, but those are beyond my expertise. See how far the two issues just noted help you and ask a follow-up question focused on df in mixed models if this doesn't help adequately.

$\endgroup$
3
  • $\begingroup$ Good points, thanks. It seems to me, after some further reading, to be a convolution of both your observations and the difficulties of asserting df in a mixed model. $\endgroup$ Dec 5, 2015 at 14:45
  • $\begingroup$ @evoked_potential : Don't forget that you can always inspect the code for these R packages to get a better handle on how they work. $\endgroup$
    – EdM
    Dec 5, 2015 at 14:59
  • $\begingroup$ I don't think glht can output the unadjusted P values -- at least that isn't it's principle purpose. If you use the lsmeans package directly, e.g. pairs(lsmeans(finalModel, ~ factorA), adjust ="none") you can get the K-R d.f. and the unadjusted P values. $\endgroup$
    – Russ Lenth
    Dec 5, 2015 at 16:24

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.