testing logistic regression coefficients using $t$ and residual deviance degrees of freedom

Question

Summary: Is there any statistical theory to support the use of the $t$-distribution (with degrees of freedom based on the residual deviance) for tests of logistic regression coefficients, rather than the standard normal distribution?

---

Some time ago I discovered that when fitting a logistic regression model in SAS PROC GLIMMIX, under the default settings, the logistic regression coefficients are tested using a $t$ distribution rather than the standard normal distribution.$^1$ That is, GLIMMIX reports a column with the ratio $\hat{\beta}_1/\sqrt{\text{var}(\hat{\beta}_1)}$ (which I will call $z$ in the rest of this question), but also reports a "degrees of freedom" column, as well as a $p$-value based on assuming a $t$ distribution for $z$ with degrees of freedom based on the residual deviance -- that is, degrees of freedom = total number of observations minus number of parameters. At the bottom of this question I provide some code and output in R and SAS for demonstration and comparison.$^2$

This confused me, since I thought that for generalized linear models such as logistic regression, there was no statistical theory to support the use of the $t$-distribution in this case. Instead I thought what we knew about this case was that

• $z$ is "approximately" normally distributed;
• this approximation might be poor for small sample sizes;
• nevertheless it cannot be assumed that $z$ has a $t$ distribution like we can assume in the case of normal regression.

Now, on an intuitive level, it does seem reasonable to me that if $z$ is approximately normally distributed, it might in fact have some distribution that is basically "$t$-like", even if it is not exactly $t$. So the use of the $t$ distribution here does not seem crazy. But what I want to know is the following:

1. Is there in fact statistical theory showing that $z$ really does follow a $t$ distribution in the case of logistic regression and/or other generalized linear models?
2. If there is no such theory, are there at least papers out there showing that assuming a $t$ distribution in this way works as well as, or maybe even better than, assuming a normal distribution?

More generally, is there any actual support for what GLIMMIX is doing here other than the intuition that it is probably basically sensible?

R code:

summary(glm(y ~ x, data=dat, family=binomial))

R output:

Call:
glm(formula = y ~ x, family = binomial, data = dat)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.352  -1.243   1.025   1.068   1.156  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.22800    0.06725   3.390 0.000698 ***
x           -0.17966    0.10841  -1.657 0.097462 .  
---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1235.6  on 899  degrees of freedom
Residual deviance: 1232.9  on 898  degrees of freedom
AIC: 1236.9

Number of Fisher Scoring iterations: 4

SAS code:

proc glimmix data=logitDat;
    model y(event='1') = x / dist=binomial solution;
run;

SAS output (edited/abbreviated):

The GLIMMIX Procedure

               Fit Statistics

-2 Log Likelihood            1232.87
AIC  (smaller is better)     1236.87
AICC (smaller is better)     1236.88
BIC  (smaller is better)     1246.47
CAIC (smaller is better)     1248.47
HQIC (smaller is better)     1240.54
Pearson Chi-Square            900.08
Pearson Chi-Square / DF         1.00


                       Parameter Estimates

                         Standard
Effect       Estimate       Error       DF    t Value    Pr > |t|

Intercept      0.2280     0.06725      898       3.39      0.0007
x             -0.1797      0.1084      898      -1.66      0.0978

$^1$Actually I first noticed this about mixed-effects logistic regression models in PROC GLIMMIX, and later discovered that GLIMMIX also does this with "vanilla" logistic regression.

$^2$I do understand that in the example shown below, with 900 observations, the distinction here probably makes no practical difference. That is not really my point. This is just data that I quickly made up and chose 900 because it is a handsome number. However I do wonder a little about the practical differences with small sample sizes, e.g. $n$ < 30.

Answer 1

Is there in fact statistical theory showing that z really does follow a t distribution in the case of logistic regression and/or other generalized linear models?

As far as I am aware, no such theory exists. I do regularly see hand-wavy arguments, and occasionally simulation experiments to support such an approach for some particular GLM family or another. The simulations are more convincing than the handwavy arguments.

If there is no such theory, are there at least papers out there showing that assuming a t distribution in this way works as well as, or maybe even better than, assuming a normal distribution?

Not that I recall seeing, but that's not saying much.

My own (limited) small-sample simulations suggest assuming a t-distribution in the logistic case may be substantially worse than assuming a normal:

Here, for example, are the results (as Q-Q plots) of 10000 simulations of the Wald statistic for an ordinary logistic regression (i.e. fixed-effects, not mixed) on 15 equispaced x-observations where the population parameters were both zero. The red line is the y=x line. As you see, in each case the normal is quite quite a good approximation over a good range in the middle - out to about the 5th and 95th percentiles (1.6-1.7ish), and then outside that the actual distribution of the test statistic is substantially lighter tailed than the normal.

So for the logistic case, I'd say any argument to use the t- rather than the z- seems to be unlikely to succeed on this basis, since simulations like these tend to suggest the results may tend to lie on the lighter-tailed side of the normal, rather than the heavier tailed.

[However, I recommend you don't trust my simulations any further than as a warning to beware - try some of your own, perhaps for circumstances more representative of your own situations typical of your IVs and models (of course, you need to simulate the case where some null is true to see what distribution to use under the null). I'd be interested to hear how they come out for you.]

Answer 2

Here are a few additional simulations just to expand a bit on what Glen_b already presented.

In these simulations I looked at the slope of a logistic regression where the predictor had a uniform distribution in $[-1,1]$. The true regression slope was always 0. I varied the total sample size ($N=10,20,40,80$) and the base rate of the binary response ($p=0.5,0.731,0.881,0.952$).

Here are Q-Q plots comparing the observed $z$ values (Wald statistics) to theoretical quantiles of the corresponding $t$ distribution ($df=N-2$). These are based on 1000 runs for each parameter combination. Notice that with small sample sizes and extreme base rates (i.e., the upper-right region of the figure), there were many cases the response only took on a single value, in which case $z=0$ and $p$-value $=1$.

Here are histograms showing the distributions of $p$-values for the logistic regression slopes based on those same $t$ distributions. These are based on 10,000 runs for each parameter combination. The $p$-values are grouped into bins of width 0.05 (20 bins in total). The dashed horizontal line shows the 5% mark, that is, frequency=500. Of course, one wants the distribution of $p$-values under the null hypothesis to be uniform, that is, all the bars should be right around the dashed line. Notice again the many degenerate cases in the upper-right part of the figure.

The conclusion seems to be that the use of $t$ distributions in this case can lead to severely conservative results when the sample size is small and/or when the base rate is approaching 0 or 1.

Answer 3

Nice work both of you. Bill Gould studied this in http://www.citeulike.org/user/harrelfe/article/13264166 making the same conclusions, in a standard fixed-effects binary logistic model.

Briefly, since the logistic model does not have an error term there is no residual variance to estimate hence the $t$ distribution does not apply [at least outside the context of multiple imputation adjustments].