Should degrees of freedom corrections be used for inference on GLM parameters?

Question

This question is inspired by Martijn's answer here.

Suppose we fit a GLM for a one parameter family like a binomial or Poisson model and that it is a full likelihood procedure (as opposed to say, quasipoisson). Then, the variance is a function of the mean. With binomial: $\text{var}[X] = E[X]E[1-X]$ and with Poisson $\text{var}[X] = E[X]$.

Unlike linear regression when the residuals are normally distributed, the finite, exact sampling distribution of these coefficients is not known, it is a possibly complicated combination of the outcomes and covariates. Also, using the GLM's estimate of the mean, that be used as a plugin estimate for the variance of the outcome.

Like linear regression, however, the coefficients have an asymptotic normal distribution, and so in finite sample inference we can approximate their sampling distribution with the normal curve.

My question is: do we gain anything by using the T-distribution approximation to the sampling distribution of the coefficients in finite samples? On one hand, we know the variance yet we don't know the exact distribution, so a T approximation seems like the wrong choice when a bootstrap or jackknife estimator could properly account for these discrepancies. On the other hand, perhaps the slight conservatism of the T-distribution is simply preferred in practice.

Answer 1

Short answer: Not a full answer yet, but you might be interested in the following distributions related to the linked question: It compares z-test (as also used by glm) and t-test

    layout(matrix(1:2,1,byrow=TRUE))

    # trying all 100 possible outcomes if the true value is p=0.7
    px <- dbinom(0:100,100,0.7)
    p_model = rep(0,101)
    p_model2 = rep(0,101)
    for (i in 0:100) {
      xi = c(rep(1,i),rep(0,100-i))
      model = glm(xi ~ 1, offset=rep(qlogis(0.7),100), family="binomial")
      p_model[i+1] = 1-summary(model)$coefficients[4]
      model2 <- glm(xi ~ 1, family = "binomial")
      coef <- summary(model2)$coefficients
      p_model2[i+1] = 1-2*pt(-abs((qlogis(0.7)-coef[1])/coef[2]),99,ncp=0)
    }


    # plotting cumulative distribution of outcomes z-test
    outcomes <- p_model[order(p_model)]
    cdf <- cumsum(px[order(p_model)])
    plot(1-outcomes,1-cdf, 
         ylab="cumulative probability", 
         xlab= "calculated glm p-value",
         xlim=c(10^-4,1),ylim=c(10^-4,1),col=2,cex=0.5,log="xy")
    lines(c(0.00001,1),c(0.00001,1))
    for (i in 1:100) {
      lines(1-c(outcomes[i],outcomes[i+1]),1-c(cdf[i+1],cdf[i+1]),col=2)
    #  lines(1-c(outcomes[i],outcomes[i]),1-c(cdf[i],cdf[i+1]),col=2)
    }

    title("probability for rejection with z-test \n as function of set alpha level")


    # plotting cumulative distribution of outcomes t-test
    outcomes <- p_model2[order(p_model2)]
    cdf <- cumsum(px[order(p_model2)])
    plot(1-outcomes,1-cdf, 
         ylab="cumulative probability", 
         xlab= "calculated glm p-value",
         xlim=c(10^-4,1),ylim=c(10^-4,1),col=2,cex=0.5,log="xy")
    lines(c(0.00001,1),c(0.00001,1))
    for (i in 1:100) {
      lines(1-c(outcomes[i],outcomes[i+1]),1-c(cdf[i+1],cdf[i+1]),col=2)
      #  lines(1-c(outcomes[i],outcomes[i]),1-c(cdf[i],cdf[i+1]),col=2)
    }

    title("probability for rejection with t-test \n as function of set alpha level")
    [![p-test vs t-test][1]][1]

And there is only a small difference. And also the z-test is actually better (but this might be because both the t-test and z-test are "wrong" and possibly the error of the z-test compensates this error).

Long Answer: ...