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

Question

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 ...

Answer 1

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.

Added comment (revised)

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.

Answer 2

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.

The documentation speaks about 'approximation'

For all other models, results rely on the normal approximation

• For a Gaussian family the normal distribution is an approximation in the sense that the t-distribution will approach a normal distribution for more degrees of freedom. If the degrees of freedom are low, then the approximation is less correct.

  So if you care about the small difference you will have to supply the degrees of freedom in order to force the function to use a t-distribution.

• In the case of a binomial or Poisson family it is appropriate to use a z-test instead of a t-test. In the case of a binomial family or Poisson family, the normal approximation for the distribution of the estimate does not depend on the degrees of freedom but on the number of counts.