Why does GLM with Gaussian family give different results to LM in R?

From what I understand GLM with a gaussian family should give the same results as LM in R, because they're essentially the same thing (from reading other posts).

When I run both on my data I get differences in significance levels; could someone explain why? The regression output appears the same, but when I do pairwise comparisons the results are different.

I am using the multcomp package for significance levels (I know there is a strong argument that p-values should not be used for regression outputs.) Is this discrepancy to do with multcomp or am I missing something?

lm <- lm(gene ~ condition, data = data)
summary(lm)
Call:
lm(formula = gene ~ condition, data = data)

Residuals:
Min       1Q   Median       3Q      Max
-1.49604 -0.54119  0.02927  0.61199  1.26961

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  -5.0178     0.3395 -14.781 1.84e-14 ***
groupCS       0.4658     0.4383   1.063   0.2973
groupSC      -0.1508     0.4383  -0.344   0.7335
groupSS       1.2337     0.4626   2.667   0.0128 *
---
Residual standard error: 0.8316 on 27 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared:  0.3145,    Adjusted R-squared:  0.2384
F-statistic:  4.13 on 3 and 27 DF,  p-value: 0.01562

> ph <- glht(lm, mcp(group="Tukey"))
> summary(ph)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lm(formula = gene ~ condition, data = data)

Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
CS - CC == 0   0.4658     0.4383   1.063   0.7139
SC - CC == 0  -0.1508     0.4383  -0.344   0.9856
SS - CC == 0   1.2337     0.4626   2.667   0.0578 .
SC - CS == 0  -0.6165     0.3920  -1.573   0.4096
SS - CS == 0   0.7680     0.4191   1.833   0.2798
SS - SC == 0   1.3845     0.4191   3.304   0.0137 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


and running the automatic glm (with Gaussian distribution);

Call:
glm(formula = gene ~ condition, data = data)

Deviance Residuals:
Min        1Q    Median        3Q       Max
-1.49604  -0.54119   0.02927   0.61199   1.26961

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  -5.0178     0.3395 -14.781 1.84e-14 ***
groupCS       0.4658     0.4383   1.063   0.2973
groupSC      -0.1508     0.4383  -0.344   0.7335
groupSS       1.2337     0.4626   2.667   0.0128 *
---
(Dispersion parameter for gaussian family taken to be 0.6914902)

Null deviance: 27.237  on 30  degrees of freedom
Residual deviance: 18.670  on 27  degrees of freedom
(1 observation deleted due to missingness)
AIC: 82.255

Number of Fisher Scoring iterations: 2

> ph <- glht(lm, mcp(group="Tukey"))
> summary(ph)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glm(formula = gene ~ condition, data = data)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
CS - CC == 0   0.4658     0.4383   1.063  0.71155
SC - CC == 0  -0.1508     0.4383  -0.344  0.98596
SS - CC == 0   1.2337     0.4626   2.667  0.03803 *
SC - CS == 0  -0.6165     0.3920  -1.573  0.39337
SS - CS == 0   0.7680     0.4191   1.833  0.25721
SS - SC == 0   1.3845     0.4191   3.304  0.00517 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

> str(data)
'data.frame':   127 obs. of  23 variables:
$$sample.id : chr "1" "3" "4" "5" ...$$ group          : Factor w/ 4 levels "CC","CS","SC",..: 1 1 1 1 1 1 2 2 2 3 ...
\$ gene            : num  -4.57 -4.91 -4.48 -5.32 -5.24 ...

• Look at the test statistics in the outputs: For lm, it uses a t statistic whereas with glm, it uses a z statistic. I suspect the slight differences are due to that. Nov 23 '21 at 16:14
• You will note that the models themselves, in terms of the regression coefficients, are identical, as are all of the contrast estimates and their standard errors.
– EdM
Nov 23 '21 at 18:02
• Looks like you are using multcomp::glht, not emmeans. So maybe you want to clarify the question? Nov 24 '21 at 3:45
• @RussLenth thanks! had a load of packages open, should have checked! Nov 24 '21 at 15:49
• FWIW, starting with the next update of emmeans, there will be no difference between these models when you use emmeans(). I assume this difference will remain for some time when using glht(). Nov 25 '21 at 19:04

Note that in the model summaries, the regression coefficients and standard errors are the same. In the comparisons outputs, the estimates and SEs are the identical. Also the t ratios in the first are identical to the z ratios in the second. The difference is in the degrees of freedom, which are taken to be 27 in the first model and infinite in the second. Infinite d.f.in a t distribution is the same as the standard normal. That affects the P values.

So the is a whole lot that is exactly the same. But the difference between these models is that the glm model does not include an estimate of degrees of freedom, and instead, asymptotic tests and CIs result.

The asymptotic results are somewhat too optimistic because they don't take into account the uncertainty in the variance estimates. Some GLMs do not involve a variance estimate at all, but the usual inferential framework is asymptotic, and both glht() and emmeans() just rely on asymptotic results.

That said, I find there are df.residual() and sigma() methods for all glm objects, and the question is when should we use that information, and when should we ignore it? As emmeans developer, I had made the choice to ignore it for glm models but not for lm models. I am considering modifying the code. One option would be to just use the residual d.f. in tests for all these models, even the ones where there is a fixed dispersion parameter. There is an argument for doing so, but I don't think this would go over well with users. I suppose what users would like is for me to use it for Gaussian and Gamma families, and ignore it otherwise.

• Isn't it confusing and complicating to say that the degrees of freedom are taken to be infinite for the second model, in order to say that it just applies a normal distribution instead of a t-distribution? Nov 25 '21 at 7:22
• In the end you write that the theory for the uncertainty is not straightforward for GLM's but for many it is, depending on whether the dispersion parameter is estimated based on the residuals (as with the Gaussian distribution, Gamma distribution or overdispersed models) or whether it is a fixed constant (as with Poisson, Binomial anddistribution). Nov 25 '21 at 7:35
• I don't think it's at all confusing to talk about infinite df; it fits right in with "asymptotic." The last paragraph isn't well stated, but the fact is that for those models with fixed scale, the tests we get are asymptotic, not exact. Nov 25 '21 at 13:49
• For Poisson and Binomial distributions the use of a z-test is an approximation but the sample size or degrees of freedom has nothing to do with that approximation. Nov 25 '21 at 17:16
• Not true. All other things being equal, the approximation is better with a larger sample size. The tests are derived from asymptotic theory. Nov 25 '21 at 17:18

From the the package documentation

glht extracts the number of degrees of freedom for models of class lm (via modelparm) and the exact multivariate t distribution is evaluated. For all other models, results rely on the normal approximation. Alternatively, the degrees of freedom to be used for the evaluation of multivariate t distributions can be given by the additional df argument to modelparm specified via ….

So, the difference is in computing the p-values with a Gaussian distribution instead of a t-distribution (z-test instead of t-test) as the sample distribution of the parameters when the null hypothesis is true.

The rest is the same. All the other values, the estimates, the standard errors, the z-statistic/t-statistic are the same.

For a generalized linear model with Gaussian family you should actually use a t-distribution/t-test for a hypothesis test instead of a z-test. Using a generalized linear model instead of ordinary least squares (linear model) does not matter in this.