# Report p-values from model fitting or glht?

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.

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)
##   7.811903e-11


If we instead use a Z-test:

2*pnorm(6.9678,lower.tail=FALSE)
##  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).

• This is incredibly useful, thanks! I have two issues though: (1) why is the unadjusted-SE/t-tests more correct than the z-test (glht)? In a more general case with more factors in the model, I would use glht because of the adjustment for multiple comparisons, right? (2) In my case, given the numbers I get (p=0.0274 and p=0.0493), what would you conclude? – Cristiano Jun 7 '18 at 21:12

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.

• Thanks! I understand the issue of rounding due to small sample size, and also that of the arbitrary threshold. In practical terms though, what would you conclude from this analysis, that grf1 is significant or that is not? – Cristiano Jun 7 '18 at 20:56
• Practically, you should ignore both the p-value and the idea of statistical significance. Both estimate 3.0 +/- 1.4, for a (rounded) 95% confidence interval of 0.2 to 5.8. Is an average difference somewhere between 0.2 and 5.8 practically meaningful? – Aaron - Reinstate Monica Jun 7 '18 at 21:05
• This is one side of the coin. The flip side would be that since this is a repeated measure design, I want to know whether grf1 affects the outcome variable coherently across subjects. My understanding, from a non statistician, is that for this I should look at the significance of the test. Correct? If so, what would you conclude given my numbers? – Cristiano Jun 7 '18 at 21:21
• @Cristiano No, the fixed-effects coefficients are the estimates of the marginal (average) effect. If you want to investigate the per subject influence, you need to look at the distribution of the residuals per subject and, depending on exactly what you want to answer, the BLUPs. – vkehayas Jun 8 '18 at 9:20