# How does one compare two nested quasibinomial GLMs?

Lets say I have two models: Model 1 and Model 2, both of which are used to fit a quasibinomial GLM on some 0/1 response data (that I believe has overdispersion, hence quasibinomial GLM instead of logistic regression).

Model 1 is a subset of Model 2. So there are k explanatory variables in Model 2 that are not in Model 1, and I want to test whether the addition of these explanatory variables provides an improvement (in-sample).

In regular logistic regression, my understanding is that I could use deviance tests to compare nested models. However, from what I've read, this is connected to the likelihood of the model, which is not available in the same way with a quasibinomial regression, since we don't have a complete log likelihood. As such, I'm unsure of whether a deviance test has the same distributional properties that it would in the general GLM case.

So as a result, my questions are as follows:

1. Is a deviance test still valid when using quasibinomial regression? If so, how is it calculated, with the addition of estimated dispersion parameters?

2. Are there other ways to formally test whether the addition of the k variables in Model 2 improve performance. Is simply using a Wald test on each of the added variables (with some sort of correction to account for the multiple tests) be valid?

I've looked around here and on other places for an answer to this question, but I haven't seen anything conclusive.

Yes, a deviance test is still valid. Some more details: Since the general theory is not specific for binomial models, I will start out with some general theory, but use binomial examples (and R.)

GLM's is based on the exponential dispersion model $$f(y_i;\theta_i,\phi)= \exp\left\{ w_i [y_i \theta_i -\gamma(\theta_i)]/\phi +\tau(y_i,\phi/w_i)\right\}$$ where $$y_i$$ is the observation, $$\theta_i$$ parameter which depends on a linear predictor $$\eta_i=x^T\beta$$, $$\phi$$ a scale parameter and $$w_i$$ a prior weight. To understand the notation, think about a normal theory model, which glm's generalize. There $$\phi$$ is the variance $$\sigma^2$$, and if $$y_i$$ is the mean of a group of $$w_i$$ independent observations with the same covariables, then the variance is $$\phi/w_i$$. The last term $$\tau(y_i,\phi/w_i)$$ is often of little interest since it does not depend on the interest parameters $$\theta_i$$ (or $$\beta$$,) so we will treat it cavalierly.

So the binomial case. If we have an observation $$y_i^* \sim \mathcal{Binom}(w_i,p_i)$$ then we will treat $$y_i=y_i^*/w_i$$ as the observation, so that the expectation of $$y_i$$ is $$p_i$$ and its variance $$\frac{p_i(1-p_i)}{w_i}$$. The binomial pmf can then be written as $$f(y_i;\theta_i,\phi)=\exp\left\{ w_i[y_i\theta_i-\log(1+e^{\theta_i})]/\phi + \log\binom{w_i/\phi}{y_i w_i/\phi} \right\}$$ where $$\phi=1$$ and $$\theta_i=\log\frac{p_i}{1-p_i}$$. We can identify $$\gamma(\theta_i)=\log(1+e^{\theta_i})$$ and $$\tau(y_i,\phi/w_i) = \log\binom{w_i/\phi}{y_i w_i/\phi}$$. This form is chosen such that we can get the quasi-model just by allowing $$\phi>0$$ to vary freely.

The quasi-likelihood we then get from this model, is constructed *to function as a likelihood for the $$\theta_i$$ (or $$\beta$$) parameters, it will not work as a likelihood for $$\phi$$. This means that the quasi-likelihood function shares enough of the properties of a true likelihood function that the usual likelihood asymptotic theory still goes through, see also Idea and intuition behind quasi maximum likelihood estimation (QMLE). Since it does not have this properties as a function of $$\phi$$, inference about $$\phi$$ must be treated outside that framework. Specifically, there is no reason to hope that maximizing the qlikelihood in $$\phi$$ to give good results.

Now, the analysis of deviance. We define the saturated model S by giving each observation its own parameter, so setting $$\hat{\mu}_i=\gamma'(\hat{\theta}_i)=y_i$$. Then by assuming for the moment that $$\phi=1$$ we get $$D_M=2\sum_i \left\{ w_i[( y_i \theta(y_i)-\gamma(\theta(y_i)))-( y_i\hat{\theta}_i-\gamma(\hat{\theta_i }) ) ]\right\}$$ which is twice the loglikelihood-ratio for testing the reduced model M within the saturated model S. Note that this does not depend on the function $$\tau$$ at all. For the case of normal-theory models, this is the residual sum of squares (RSS), which is not a function of the scale parameter $$\phi=\sigma^2$$ either.

$$D_M/\phi$$ is the scaled deviance while $$D_M$$ often is called the residual deviance, since in normal models it corresponds to the RSS. In normal models we have $$D_M/\phi \sim \chi^2_{n-p}$$ so an unbiased estimator of the variance parameter $$\phi$$ in this case is $$\hat{\phi}=D_M/(n-p)$$ and this might hold as an approximation also in other cases, but often better is $$\tilde{\phi}=\frac1{n-p}\sum_i \frac{(y_i-\hat{\mu}_i)^2}{V(\hat{\mu_i})/w_i}$$ where $$V$$ is the variance function, in the binomial case $$V(\mu)=\mu(1-\mu)$$. In the binomial case, this is considered to be better, and is the scale estimate used by R.

If we are interested in a submodel $$M_0 \subset M$$, with $$q < p$$ regression parameters, then the likelihood ratio test is $$\frac{D_{M_0}-D_M}{\phi} \stackrel{\text{approx}}{\sim} \chi^2_{p-q}$$ and with estimated scale we might use $$\frac{D_{M_0}-D_M}{\hat{\phi}(p-q)} \stackrel{\text{approx}}{\sim} \mathcal{F}_{p-q,n-p}$$ in analogy with the normal theory.

So, let us look at a simulated example.

    set.seed(7*11*13) # My public seed
n <- 200
k <- 5
N <- n*k
intercept <- rnorm(n, 0, 1)
x <- rnorm(n, 1, 1.5)
beta <- 0.1
expit <- function(x) 1/(1+exp(-x))
eta <- intercept + beta*x
p <- expit(eta)
Y <- rbinom(n, k, p)


This creates overdispersion by simulating a random intercept for each of the $$n=200$$ groups of size $$k=5$$. Then we will estimate a simple model two ways, by using a binomial likelihood, and then a quasibinomial likelihood:

    mod0 <- glm( cbind(Y, k-Y) ~ x, family=binomial)
modq <- glm( cbind(Y, k-Y) ~ x, family=quasibinomial)


Then the model summaries:

     summary(mod0)

Call:
glm(formula = cbind(Y, k - Y) ~ x, family = binomial)

Deviance Residuals:
Min      1Q  Median      3Q     Max
-3.053  -1.180  -0.103   1.180   2.836

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.02787    0.07632  -0.365  0.71496
x            0.12941    0.04170   3.103  0.00192 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 441.41  on 199  degrees of freedom
Residual deviance: 431.62  on 198  degrees of freedom
AIC: 749.1

Number of Fisher Scoring iterations: 3

> summary(modq)

Call:
glm(formula = cbind(Y, k - Y) ~ x, family = quasibinomial)

Deviance Residuals:
Min      1Q  Median      3Q     Max
-3.053  -1.180  -0.103   1.180   2.836

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.02787    0.10117  -0.275   0.7832
x            0.12941    0.05529   2.341   0.0202 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 1.757479)

Null deviance: 441.41  on 199  degrees of freedom
Residual deviance: 431.62  on 198  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 3


Compare the two summaries. They are very similar, the differences are in coefficient standard errors, and the printed scale parameter estimate, and lacking AIC of the modq summary. Check that you can calculate, "by hand", the standard errors for the quasimodel modq from the standard errors for mod0 and the estimated scale.

The printed deviances, and deviance residuals, are identical. This is because the residual deviance is defined by taking $$\phi=1$$ in both cases. The null deviance is the residual deviance for the null model, the model with only an intercept. The scaled deviance is not printed, but can be calculated from the output.

The analysis of deviance is calculated by the anova() function. Here we will see differences. First the model based on a binomial likelihood:

     anova(mod0, test="Chisq")
Analysis of Deviance Table

Response: cbind(Y, k - Y)

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL                   199     441.41
x     1   9.7883       198     431.62 0.001756 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

so here the regression seems significant. Then for the quasi-model:

anova(modq, test="F")
Analysis of Deviance Table

Response: cbind(Y, k - Y)

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev      F  Pr(>F)
NULL                   199     441.41
x     1   9.7883       198     431.62 5.5695 0.01925 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


What is printed as F here is (in this case) the scaled deviance (since $$p-q=1$$.)

(I will come back to your second question)

• I appreciate the very detailed response. Thank you. Jul 16 '19 at 20:50