Report p-values from model fitting or glht?

Question

I am running a mixed model to compare two groups in a repeated measures design. The model is very simple:

> linM11 <- lme(values ~ grf, random = ~1|id, data=dat_trf, na.action=na.omit, method = "ML", control=lCtr )

where grf is a factor with two levels. I am asking if the two levels are significantly different. When I call summary() I obtain (among the other things):

> summary(linM11)

Fixed effects: values ~ grf
                   Value   Std.Error   DF   t-value   p-value
(Intercept)    -8.513064   0.9908567   16   -8.59162   0.0000
grf1            3.027705   1.4158346   15    2.13846   0.0493

which makes me think that the two groups are barely significantly different. But If I run this post-hoc test I get this other (quite different) result:

> ph_conditional <- c("grf1 = 0");
> linM.ph <- glht(linM11, linfct = ph_conditional);
> summary(linM.ph)

Linear Hypotheses:
          Estimate   Std. Error   z value      Pr(>|z|) 
grf1 == 0    3.028        1.372     2.206      0.0274 *

Which one should I trust? I am always confused on whether I should use the p-values of the model fit or those of the post-hoc tests. If I had multiple tests, I would certainly run the post-hoc and adjust for multiple comparisons (with glht). But here there is only one test, and I am not sure why I get so different results, and which one I should trust.

Answer 1

There are two things going on here:

1. the difference between t-tests and Z-tests (as pointed out by @vkehayas); t-tests account for the uncertainty in the estimate of the standard error, so should be preferred to Z-tests where available.
2. the fact that summary.lme by default adjusts the residual standard error for ML estimates (while glht doesn't); ML estimation in general gives a slightly downward-biased estimate of the standard error (by a factor $\sqrt{(n-p)/n}$), so this adjustment should be preferred where available. This is the adjustSigma parameter of summary.lme:

adjustSigma: an optional logical value. If ‘TRUE’ and the estimation method used to obtain ‘object’ was maximum likelihood, the residual standard error is multiplied by sqrt(nobs/(nobs - npar)), converting it to a REML-like estimate. ... Default is ‘TRUE’.

Both of these adjustments should in general make little difference unless your sample is small, but both adjustSigma=TRUE and t-tests rather than Z-test are technically more correct, so in a pinch you should probably accept the results of summary(.) rather than those of glht().

If you have a factor with more than two levels (so that you need to summarize the joint significance of multiple parameters), you can use anova(), which uses F tests (the analog of t-tests) and includes an adjustSigma option: if you want to do more complicated post hoc testing (e.g. Tukey pairwise comparisons), you will probably need to use glht() and accept that your answers will be slightly anticonservative/optimistic.

Try to keep in mind that $p=0.0274$ and $p=0.0493$ (from your example) are not very different from each other; in practice people behave as if there's a magic line at $p=0.05$, but there isn't.

Here's an example:

library(multcomp)
library(nlme)
data("sleepstudy",package="lme4")

m2 <- lme(Reaction~Days, random = ~Days|Subject,
          data=sleepstudy, method="ML")

Results (fancy code with printCoefmat() etc. is just to isolate the information we want from summary(m2)):

printCoefmat(summary(m2)$tTab["Days",,drop=FALSE])
##         Value Std.Error       DF t-value   p-value    
## Days  10.4673    1.5106 161.0000   6.929 9.651e-11 ***

With adjustSigma=FALSE, the standard error changes from 1.5106 to 1.5022:

printCoefmat(summary(m2,adjustSigma=FALSE)$tTab["Days",,drop=FALSE])
##         Value Std.Error       DF t-value   p-value    
## Days  10.4673    1.5022 161.0000  6.9678 7.811e-11 ***

A direct calculation of the p-value using the unadjusted sigma:

2*pt(6.9678,161,lower.tail=FALSE)
## [1]  7.811903e-11

If we instead use a Z-test:

2*pnorm(6.9678,lower.tail=FALSE)
## [1] 3.219354e-12

This agrees with the answer we get from glht:

    summary(glht(m2, linfct=c("Days=0")))
##           Estimate Std. Error z value Pr(>|z|)    
## Days == 0   10.467      1.502   6.968 3.22e-12 ***

In your case most of the difference is from the t- vs Z-test distinction; 2*pt(2.206,df=15,lower.tail=FALSE) (i.e. using the unadjusted standard error with a t test) gives $p=0.043$, most of the way from $p=0.027$ (summary(.) result) to $p=0.049$ (glht(.) result).

Answer 2

These results do not look different to me, in the sense that I would not change my interpretation based on so minuscule discrepancies. Any difference you think you detect can be considered rounding error. Perhaps the most important factor that needs addressing in the issue you are facing is the psychological bias towards an arbitrary threshold. Having said that, I see that one function reports t-values and the other z-values. Given you small sample size, I'd venture to guess that your sample displays minor deviations from the normal distribution that prevent the values from being more close to one another.