3
$\begingroup$

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 ...
$\endgroup$
5
  • 11
    $\begingroup$ 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. $\endgroup$ Nov 23 at 16:14
  • 6
    $\begingroup$ 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. $\endgroup$
    – EdM
    Nov 23 at 18:02
  • 1
    $\begingroup$ Looks like you are using multcomp::glht, not emmeans. So maybe you want to clarify the question? $\endgroup$
    – Russ Lenth
    Nov 24 at 3:45
  • $\begingroup$ @RussLenth thanks! had a load of packages open, should have checked! $\endgroup$ Nov 24 at 15:49
  • $\begingroup$ 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(). $\endgroup$
    – Russ Lenth
    Nov 25 at 19:04
8
$\begingroup$

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.

$\endgroup$
9
  • $\begingroup$ 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? $\endgroup$ Nov 25 at 7:22
  • $\begingroup$ 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). $\endgroup$ Nov 25 at 7:35
  • 1
    $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Nov 25 at 13:49
  • $\begingroup$ 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. $\endgroup$ Nov 25 at 17:16
  • $\begingroup$ Not true. All other things being equal, the approximation is better with a larger sample size. The tests are derived from asymptotic theory. $\endgroup$
    – Russ Lenth
    Nov 25 at 17:18
1
$\begingroup$

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.

$\endgroup$

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.