Quasi-likelihood in GLMs

To fit a GLM you choose a distribution — Poisson, binomial, Gamma — and maximise the likelihood. But look closely at the equations the algorithm actually solves and the distribution almost disappears: only the mean \(\mu_i\) and the variance function \(V(\mu_i)\) are left. Robert Wedderburn noticed this in 1974 and drew the conclusion: if only the mean and the variance matter, you can specify only the mean and the variance and throw away the rest of the distribution. The result is quasi-likelihood, the basis of quasi-Poisson, quasi-binomial and of models for data that follow no known distribution at all. This article tells where it came from, proves why it works, shows the estimators and the algebraic properties that make it so convenient and — with three examples and four simulations — when it is the right tool and when it is not. Every number comes from a reproducible script.

§1What a GLM really uses from the distribution

In a GLM the response belongs to the exponential family:

\(\theta\) is the canonical parameter, \(\phi\) the dispersion parameter, and \(V(\mu) = b''(\theta)\), written as a function of \(\mu\), is the variance function: \(V(\mu) = \mu\) for the Poisson, \(\mu(1-\mu)\) for the binomial (for proportions), \(\mu^2\) for the Gamma, \(1\) for the normal. Now differentiate the log-likelihood of one observation with respect to a coefficient \(\beta_j\), using the chain rule \(\beta_j \to \eta \to \mu \to \theta\):

(The middle term comes from \(d\mu/d\theta = b''(\theta) = V(\mu)\).) Summing over observations, the score equations that define the maximum likelihood estimator are

Read (3) carefully. There is no \(b(\cdot)\), no \(c(y,\phi)\), no shape of the density. Only the mean \(\mu_i(\boldsymbol\beta)\) and the variance \(V(\mu_i)\) appear. Two distributions with the same mean and the same variance function would give exactly the same \(\hat{\boldsymbol\beta}\). This is the observation everything else follows from.

§2The origin: Wedderburn and leaf blotch

GLMs were proposed by John Nelder and Robert Wedderburn in 1972. Two years later Wedderburn — a young statistician at Rothamsted Experimental Station — published “Quasi-likelihood functions, generalized linear models, and the Gauss–Newton method” in Biometrika. The concrete problem was a dataset on leaf blotch (Rhynchosporium secalis) in barley: the percentage of leaf area affected, for 10 varieties at 9 sites.

These data are proportions, but they are not binomial: there are no “\(m\) trials” and “\(y\) successes”, just a continuous area. There is no natural distribution for them. Even so, a lot can be said without choosing one — first about the mean, then about the variance.

The mean. The experiment crosses two factors, variety \(i\) and site \(j\), with one observation per combination. The natural model is

with 1 + 9 + 8 = 18 parameters for 90 observations. The logit link does two things. It keeps \(\mu_{ij}\) between 0 and 1 whatever \(\alpha_i\) and \(\beta_j\) are. And it gives the effects a multiplicative reading: \(e^{\alpha_i - \alpha_{i'}}\) is the ratio of the odds of infection of two varieties, the same at every site. A resistant variety cuts the odds by the same factor at a dry site and at a wet one — the no-interaction assumption, which with one observation per cell could not even be tested.

The variance: what is compulsory. A proportion lives in \([0, 1]\), and that alone bounds the variance. Since \(Y^2 \le Y\) on that interval, \(\mathbb E[Y^2] \le \mathbb E[Y] = \mu\) and

Every variable in \([0,1]\) obeys this inequality, binomial or not. Consequence: the variance must go to zero as the mean proportion goes to 0 or 1. A site where the disease barely appears cannot have widely scattered proportions. But that still leaves many possible shapes, and two are natural candidates: \(\mu(1-\mu)\), the binomial shape, and \(\mu^2(1-\mu)^2\).

The variance: which shape. The difference between them shows up in the coefficient of variation, the standard deviation relative to the mean, near zero:

With the binomial shape, the coefficient of variation grows like \(1/\sqrt\mu\): a variety losing 1% of its leaf on average would have relative fluctuations more than three times those of one losing 10% (\(\sqrt{99}/\sqrt{9} \approx 3.3\)). With \(\mu^2(1-\mu)^2\), the coefficient of variation is roughly constant for small proportions: a disease on 1% of the leaf fluctuates, proportionally, as much as one on 10%. That is how an epidemic, which grows multiplicatively, is expected to behave. And by symmetry the same holds near 1, swapping the diseased area for the healthy area \(1-\mu\).

There is a second reading, which ties the variance to the link. By the delta method, \(\operatorname{Var}(\operatorname{logit} Y) \approx \operatorname{Var}(Y)\,\big/\,[\mu(1-\mu)]^2\), since \(d\operatorname{logit}(\mu)/d\mu = 1/[\mu(1-\mu)]\). So

Wedderburn's variance is exactly the one that makes \(\operatorname{logit} Y\) homoscedastic: it translates to the original scale the model an agronomist would fit by transforming the data, \(\operatorname{logit} y_{ij} = \alpha_i + \beta_j + \varepsilon_{ij}\) with constant-variance errors. (The binomial shape does the same for the arcsine square-root transform.) The advantage of not transforming is that the model is still for \(\mathbb E[Y]\), not \(\mathbb E[\operatorname{logit} Y]\), and that proportions equal to 0 — common in disease data — break nothing.

Wedderburn fitted the model with the binomial variance, saw that it did not describe the residuals well and proposed

No known exponential-family distribution has this variance function. Curiously, the part that depends on \(\mu\) does have the right form: integrating \(d\theta/d\mu = 1/V(\mu)\) gives a “canonical parameter” \(\theta = 2\operatorname{logit}\mu - 1/\mu + 1/(1-\mu)\), and the function \(Q\) of §3.2 can be written as \(y\theta - b(\theta)\), with \(b = \operatorname{logit}\mu + 1/(1-\mu)\). What is missing is a term \(c(y,\phi)\) that makes \(\exp\{\cdot\}\) integrate to 1 — that is, an actual density. Wedderburn did not need it. The paper's question was: can \(\boldsymbol\beta\) be estimated by solving (3) with an arbitrary \(V\), even with no likelihood behind it? And does the resulting estimator have good properties? The answer, yes and yes, is the content of §3–§5. §8 redoes the example with simulated data in the same format.

Others completed the idea:

• Peter McCullagh (1983), in the Annals of Statistics, gave the asymptotic theory: consistency, normality and the covariance matrix of the estimator, including for dependent observations.
• Kung-Yee Liang and Scott Zeger (1986) extended the idea to longitudinal and clustered data: GEE (generalized estimating equations) are multivariate quasi-likelihood equations with a “working” correlation matrix.
• John Nelder and Daryl Pregibon (1987) proposed the extended quasi-likelihood, which makes it possible to compare different variance functions (§4.6).
• The book by McCullagh and Nelder (1989), ch. 9, consolidated the method — it is where most statisticians learned it, and where the original leaf-blotch data are printed.

Wedderburn died in 1975, a year after the paper, before turning 30. The quasi-Poisson and quasi-binomial that are now in every package — family = quasipoisson in R — are a direct consequence of it.

§3The mathematical foundation

3.1 The quasi-score

Assume only two things about an observation \(Y\):

with \(V\) a known function and \(\phi \gt 0\) possibly unknown. Nothing about the shape of the distribution. Define the quasi-score by imitating the term in (2):

A true likelihood score, \(\partial\ell/\partial\mu\), has three properties that support the whole asymptotic theory of maximum likelihood. The quasi-score has the same three, using only (5).

Property 1 — zero mean.

Property 2 — the variance. Since \(\mu\) and \(V(\mu)\) are constants,

Property 3 — the information identity. Differentiate \(U\) with respect to \(\mu\) (quotient rule):

The term with \(V'(\mu)\) vanishes because it multiplies \(Y - \mu\), which has mean zero. The equality “\(-\mathbb E[\text{derivative of the score}] = \operatorname{Var}(\text{score})\)” is the Bartlett identity that, for a likelihood, makes the Fisher information both the expected curvature and the variance of the score. It is what gives \(\operatorname{Cov}(\hat{\boldsymbol\beta}) \approx \mathcal I^{-1}\). Since (7)–(9) hold, every asymptotic argument of maximum likelihood that uses only these three properties also holds for the quasi-score.

3.2 The quasi-likelihood

If \(U\) behaves like the derivative of a log-likelihood, integrate it:

This is the quasi-(log-)likelihood. By construction \(\partial Q/\partial\mu = U\) and, by the chain rule, \(\partial Q/\partial\beta_j = \sum_i U_i\,\partial\mu_i/\partial\beta_j\): maximising \(Q\) is solving (3). The lower limit \(y\) only fixes the constant — it makes \(Q(y;y) = 0\) — and affects nothing that depends on \(\boldsymbol\beta\). For each variance function the integral gives (with \(\phi = 1\) and up to terms that depend only on \(y\)):

When an exponential-family distribution with that \(V\) exists, \(Q\) equals its log-likelihood (up to terms in \(y\)). The last row is Wedderburn's: a “likelihood” with no distribution.

Example: checking Wedderburn's row. Differentiate each term of \(Q = (2y-1)\log\frac{\mu}{1-\mu} - \frac{y}{\mu} - \frac{1-y}{1-\mu}\):

Expanding the numerator: \((2y\mu - 2y\mu^2 - \mu + \mu^2) + (y - 2y\mu + y\mu^2) - (\mu^2 - y\mu^2) = y - \mu\). So \(\partial Q/\partial\mu = (y-\mu)/[\mu^2(1-\mu)^2]\), as it should be. Any row of the table can be checked the same way.

3.3 The quasi-deviance

As with a likelihood, the measure of fit is the distance to the saturated model (\(\mu_i = y_i\)):

Multiplying by \(\phi\) removes the dispersion parameter: \(D\) depends only on \(V\). For \(V = \mu\) it is the Poisson deviance, \(2\sum[y\log(y/\mu) - (y - \mu)]\), and so on. It is used to compare nested models with the same \(V\) (§4.5) — and not to compare different variance functions, because each \(V\) defines a different scale.

§4The estimators

4.1 The coefficients: the quasi-score equations

In matrix notation, with \(\mathbf D = \partial\boldsymbol\mu/\partial\boldsymbol\beta^\top\) (an \(n\times p\) matrix, \(D_{ij} = \partial\mu_i/\partial\beta_j\)) and \(\mathbf V = \operatorname{diag}\{V(\mu_i)\}\), equations (3) become

\(\hat{\boldsymbol\beta}\) is the root of (13). Note that \(\phi\) multiplies the whole equation and can be cancelled: \(\hat{\boldsymbol\beta}\) does not depend on \(\phi\).

4.2 The algorithm: Fisher scoring = IRLS = Gauss–Newton

Fisher scoring replaces the derivative of \(\mathbf U\) by its expectation, \(-\mathbb E[\partial\mathbf U/\partial\boldsymbol\beta^\top] = \mathbf D^\top\mathbf V^{-1}\mathbf D/\phi\) (the matrix version of (9)), and iterates

with everything evaluated at \(\boldsymbol\beta^{(r)}\). Once again \(\phi\) cancels. To see that this is the GLM IRLS, write \(\mathbf D = \boldsymbol\Delta\mathbf X\), with \(\boldsymbol\Delta = \operatorname{diag}(d\mu_i/d\eta_i)\), and define the weights and the working response

Then \(\mathbf D^\top\mathbf V^{-1}\mathbf D = \mathbf X^\top\mathbf W\mathbf X\) and \(\mathbf D^\top\mathbf V^{-1}(\mathbf y-\boldsymbol\mu) = \mathbf X^\top\mathbf W\boldsymbol\Delta^{-1}(\mathbf y-\boldsymbol\mu)\), and (14) becomes

a weighted least-squares regression of \(\mathbf z\) on \(\mathbf X\), repeated until convergence. It is the same algorithm that fits any GLM; the only thing that changes from one quasi model to another is the function \(V\) inside the weights. The “Gauss–Newton” in Wedderburn's title comes from here: (14) is exactly the Gauss–Newton step for minimising \(\sum_i (y_i - \mu_i)^2/V_i\) with the \(V_i\) frozen at the current iteration — weighted nonlinear least squares.

Example. On the insect dataset of §6 (four sprays, log link, \(V(\mu) = \mu\)), the iterations of (16) from a naive start are:

Starting point: \(\mu_i^{(0)} = (y_i + \bar y)/2 + 0.1\). By the third iteration the coefficients are right to the 5th decimal; the last two only confirm convergence. No row depends on \(\phi\).

4.3 The dispersion \(\phi\): the Pearson estimator

By (5), \(\mathbb E[(Y_i - \mu_i)^2/V(\mu_i)] = \phi\) for every observation. This suggests the moment estimator

Pearson's \(X^2\) statistic divided by the degrees of freedom. The \(n - p\) instead of \(n\) compensates for the \(p\) estimated parameters, exactly as in \(s^2 = \text{RSS}/(n-p)\) in linear regression (which is the case \(V = 1\)). An alternative is \(D/(n-p)\), the deviance over the degrees of freedom, but it is biased when counts are small or proportions are near 0 and 1; Pearson's \(\hat\phi\) is consistent whenever (5) holds, and it is the default in R and in statsmodels.

4.4 The covariance of \(\hat{\boldsymbol\beta}\)

Expand (13) around the true value: \(\mathbf 0 = \mathbf U(\hat{\boldsymbol\beta}) \approx \mathbf U(\boldsymbol\beta) - \mathcal I(\hat{\boldsymbol\beta} - \boldsymbol\beta)\), with \(\mathcal I = \mathbf D^\top\mathbf V^{-1}\mathbf D/\phi\). So \(\hat{\boldsymbol\beta} - \boldsymbol\beta \approx \mathcal I^{-1}\mathbf U(\boldsymbol\beta)\) and, since the observations are independent,

The middle step uses \(\operatorname{Cov}(\mathbf U) = \mathcal I\), the matrix version of (8)–(9). In practice \(\phi\) is replaced by \(\hat\phi\) from (17). Result: the standard errors of the quasi model are those of the “nominal” GLM multiplied by \(\sqrt{\hat\phi}\). For quasi-Poisson, the Poisson ones times \(\sqrt{\hat\phi}\); for quasi-binomial, the binomial ones times \(\sqrt{\hat\phi}\).

And if the variance function is wrong — if the true variance is \(\boldsymbol\Sigma \ne \phi\mathbf V\)? Then \(\operatorname{Cov}(\mathbf U) \ne \mathcal I\) and the correct covariance is the sandwich

estimated by replacing \(\boldsymbol\Sigma\) with \(\operatorname{diag}\{(y_i - \hat\mu_i)^2\}\). It assumes no form for the variance, but it is noisy when \(n\) is small (§9 measures this).

4.5 Tests: \(t\) and \(F\) instead of \(z\) and \(\chi^2\)

Since \(\phi\) is estimated, inference follows the linear-regression template, not the likelihood one:

• One coefficient: \(t = \hat\beta_j/\text{SE}(\hat\beta_j)\), compared with \(t_{n-p}\).
• Nested models (the smaller with \(p_0\) parameters, the larger with \(p_1\)):

with \(\hat\phi_1\) estimated in the larger model. If \(\phi\) were known and equal to 1, \((D_0 - D_1)\) would be the likelihood-ratio test with \(\chi^2_{p_1-p_0}\). Dividing by \(\hat\phi\) “discounts” the overdispersion, and using \(F\) instead of \(\chi^2\) discounts the uncertainty of having estimated \(\hat\phi\).

4.6 When variance functions must be compared

\(Q\) cannot choose between \(V(\mu) = \mu(1-\mu)\) and \(V(\mu) = \mu^2(1-\mu)^2\): the two live on different scales. Nelder and Pregibon (1987) solved this with the extended quasi-likelihood

where \(d_i\) is observation \(i\)'s contribution to the deviance (12). The \(\log V(y_i)\) term is what allows different \(V\)s to be compared: \(Q^+\) is the saddlepoint approximation to the log-likelihood, and for the normal and inverse Gaussian it is exact. Maximising over \(\phi\) gives \(\hat\phi = D/n\). For model selection with overdispersion, ecologists also use the QAIC of Lebreton et al. (1992), \(\text{QAIC} = -2\ell/\hat c + 2k\), with \(\hat c = \hat\phi\) and \(k\) counting \(\phi\).

§5The algebraic advantages

The properties below are why quasi-likelihood is so widely used. Each one is a short calculation.

5.1 Consistency needs only the right mean

Take the expectation of (13) at the true value of \(\boldsymbol\beta\):

The calculation never used the variance. An estimating equation with zero expectation at the true value (an unbiased estimating equation) yields, under regularity conditions, a consistent estimator. So \(\hat{\boldsymbol\beta}\) is consistent even if \(V\) is wrong. A wrong variance spoils the standard errors (which (19) fixes) and the efficiency (§5.3), but it does not point the estimator at the wrong place.

5.2 \(\phi\) separates from \(\boldsymbol\beta\)

We have seen that \(\phi\) cancels in (13) and (14). Algebraically, this means the problem splits in two: first estimate \(\boldsymbol\beta\) as if \(\phi = 1\); then \(\phi\) by (17); finally, the covariance (18) is the nominal one multiplied by \(\hat\phi\). Example: in §6, quasi-Poisson returns exactly the same \(\hat\beta\) as Poisson (in the table, a single column of estimates serves both). That is why switching from Poisson to quasi-Poisson never changes a prediction; it only changes the uncertainty.

5.3 Optimality: the “Gauss–Markov” of estimating equations

Why \(\mathbf D^\top\mathbf V^{-1}(\mathbf y-\boldsymbol\mu)\) in particular? Consider any estimating equation linear in \(\mathbf y\),

with \(\mathbf A\) any \(n\times p\) matrix (which may depend on \(\boldsymbol\beta\)). By the same argument as in (18), the root has \(\hat{\boldsymbol\beta} - \boldsymbol\beta \approx (\mathbf A^\top\mathbf D)^{-1}\mathbf A^\top(\mathbf y-\boldsymbol\mu)\) and asymptotic covariance

The quasi-likelihood choice is \(\mathbf A^\star = \boldsymbol\Sigma^{-1}\mathbf D\), which gives \(\mathbf C(\mathbf A^\star) = (\mathbf D^\top\boldsymbol\Sigma^{-1}\mathbf D)^{-1}\). Claim: \(\mathbf C(\mathbf A) - \mathbf C(\mathbf A^\star)\) is positive semidefinite for every \(\mathbf A\). Proof. Define

\(\mathbf M\boldsymbol\Sigma\mathbf M^\top\) is positive semidefinite (it is a covariance). Expanding it, both cross terms equal

and the last term is also \(\mathbf C(\mathbf A^\star)\). Therefore

Among all estimators defined by equations linear in \(\mathbf y\), the quasi-likelihood one has the smallest asymptotic variance (Godambe, 1960; McCullagh, 1983). It is the generalisation of the Gauss–Markov theorem: with \(V = 1\) and the identity link, \(\mathbf A^\star \propto \mathbf X\) and (23) are the least-squares normal equations. Two practical consequences: (i) getting the shape of \(V\) right buys efficiency — §9 measures how much; (ii) the guarantee is only within the linear class: if the true distribution is known and is not in the exponential family, full maximum likelihood can be more efficient.

5.4 The exponential family is a special case

If there is an exponential-family distribution with variance function \(V\), then \(Q\) is its log-likelihood (up to terms in \(y\)) — compare (2) with (6) and integrate. In that case the quasi-likelihood estimator is the maximum likelihood estimator, with all its efficiency. Quasi-likelihood is not an approximate stand-in for the GLM: it is the general theory, and the classical GLM is the case in which there is a distribution behind it.

5.5 More freedom with the same algorithm

• Any positive \(V\), even with no distribution (like (4)), fits in the same IRLS (16).
• \(\phi\) free in both directions: \(\hat\phi \lt 1\) (underdispersion) is as admissible as \(\hat\phi \gt 1\). In the Poisson and binomial, \(\phi = 1\) is imposed.
• Responses outside the “official” support: quasi-Poisson accepts continuous non-negative \(y\) (spending, concentrations), and quasi-binomial accepts continuous proportions — \(Q\) only requires the integral (10) to exist.

§6Example 1 — quasi-Poisson: insects per plot

Four sprays (A, B, C, D) are applied to 12 plots each, and the number of live insects per plot is counted. The true means are 15, 12, 4 and 6, but plots differ from one another (microclimate, edges, vegetation), and the data were generated with \(\operatorname{Var}(Y) = 3\mu\) — exactly the case quasi-Poisson describes, with \(\phi = 3\). The first diagnostic needs no model: with groups, compare mean and variance.

12 plots per spray. Data generated with \(\operatorname{Var}(Y) = 3\mu\). Under Poisson the last column would be close to 1; here it ranges from 2 to 6 with no clear trend in the mean — the pattern quasi-Poisson assumes.

The model is \(\log\mu_{i} = \beta_0 + \beta_1[\text{B}] + \beta_2[\text{C}] + \beta_3[\text{D}]\). Fitting Poisson and quasi-Poisson:

n = 48, \(\hat\phi = X^2/(n-p) = 147.1/44 = 3.34\) (true: 3). The estimate is the same in all three columns; only the standard error changes. Quasi-Poisson = Poisson × \(\sqrt{\hat\phi} = 1.83\). The quasi interval uses the \(t_{44}\) quantile = 2.015.

Three readings:

• \(\hat\phi = 3.34\), close to the true 3. The Poisson standard errors are \(\sqrt{3.34} = 1.83\) times too small.
• The interval for B − A goes from \([-0.272,\ 0.170]\) to \([-0.467,\ 0.365]\). As a rate ratio, \(e^{\hat\beta_1} = 0.95\): the two sprays are indistinguishable in both models, but Poisson promises almost twice the precision the data have.
• The global test of the spray effect is still significant — \(F_{3,44} = 7.18\), p = 0.0005 — but the evidence is more than ten orders of magnitude less extreme than Poisson's (\(\chi^2_3 = 72.0\), p \(\approx 10^{-15}\)). With smaller effects, this gap is where false positives live.

The sandwich column lies between the other two for some coefficients and above quasi-Poisson for C − A. With 12 plots per group it estimates a separate variance per group and is much noisier than quasi-Poisson, which pools all 48 observations into a single \(\hat\phi\).

§7Example 2 — quasi-binomial: seed germination

Two seed types are germinated on plates of 30 seeds, with water or with a plant extract; eight plates per combination. The question is whether the extract acts differently on the two types — the interaction. Seeds on the same plate share moisture and temperature, so they germinate in a correlated way: the number germinated per plate is not binomial. If each plate's probability varies around the mean with intra-plate correlation \(\rho\), the variance is

that is, the binomial variance times \(\phi = 1 + (m-1)\rho\). With the same \(m\) on every plate this is exactly quasi-binomial. Here \(\rho = 0.06\) gives \(\phi = 2.74\), and the true interaction is zero.

8 plates of 30 seeds per combination; plate probabilities drawn from a Beta with intra-plate correlation \(\rho = 0.06\), so \(\phi = 1 + (m-1)\rho = 2.74\). The true interaction is zero. Standard errors in brackets.

This dataset is instructive precisely because it is “unlucky”. The type 1 · extract combination germinated 15% against 31% expected, which produced an estimated interaction of 1.2 where the true one is zero, and \(\hat\phi = 1.54\), well below the true 2.74. The binomial rejects with p \(\approx 3\times 10^{-5}\); quasi-binomial also rejects, with p = 0.002. A single sample cannot tell which method is calibrated. For that, we repeat the same design 2,000 times, always with no interaction, and count rejections at 5%:

• binomial (\(\chi^2_1\)): rejects 22.6% of the time — more than four times the nominal level;
• quasi-binomial (\(F_{1,28}\)): rejects 4.4% — within Monte Carlo error of 5%.

Across those 2,000 repetitions \(\hat\phi\) averaged 2.77 (true 2.74), with 80% of values between 1.92 and 3.68: our sample was in the lowest 10%. Quasi-binomial is calibrated in aggregate, and the variability of \(\hat\phi\) with 28 degrees of freedom is the price of not knowing \(\phi\) — the \(F\) test already accounts for it.

§8Example 3 — Wedderburn's variance: continuous proportions

Back to 1974. We simulate data in the leaf-blotch format: 10 varieties × 9 sites, the proportion of leaf area affected with \(\operatorname{logit}\mu_{ij} = \alpha_i + \beta_j\) — from almost zero to almost 1 — and a Beta distribution whose variance is exactly \(\mu^2(1-\mu)^2\). We fit the same mean model with two variance functions: the binomial \(\mu(1-\mu)\) (quasi-binomial, what many people would do) and Wedderburn's (4).

First five rows: mean squared standardised residual in each fifth of the observations, sorted by fitted value. If the variance function is right, all are close to 1. Larger \(Q^+\) = better; unlike the quasi-deviance it can compare different variance functions.

• The residual pattern. With the binomial variance, the mean squared residual rises from 0.68 in the lowest fifth to 1.55 in the fourth fifth: the assumed variance is too large for small proportions relative to those in the middle. With Wedderburn's, the five fifths stay between 0.84 and 1.14. It is the same diagnosis Wedderburn made in 1974. (The pattern is less dramatic than the theoretical ratio \(\mu(1-\mu)\) would suggest, because the wrong \(V\) itself distorts the fit: with weights \(\propto\mu(1-\mu)\), the extreme observations barely count and end up poorly fitted.)
• \(\hat\phi\) does not compare models. 0.158 versus 1.113 does not mean the first fits “better”; the two \(V\)s are on different scales. \(Q^+\) is what compares them: 121.8 against 97.1, in favour of \(\mu^2(1-\mu)^2\).
• Efficiency. The contrast variety 10 − variety 1 (true 3.5) comes out at 4.00 with SE 0.50 under the right variance, and 4.17 with SE 0.68 under the binomial one. The sandwich tells the same story (0.35 against 0.46): even with robust standard errors, the wrong variance wastes information — the theorem of §5.3 at work.

§9When quasi models are good

To measure rather than opine, four simulation studies with 2,000 datasets each. S1 and S3 generate counts with \(\operatorname{Var} = 3\mu\) (the quasi-Poisson assumption is true); S2 generates continuous positive data with \(\operatorname{Var} = 0.25\,\mu^2\) and fits them with the wrong \(V\) (\(\mu\)) and the right one (\(\mu^2\)).

2000 datasets per scenario; \(\log\mu = 1 + 0.5x\) (S1, S3) or \(1 + 1.5x\) (S2), \(x \sim U(0, 2)\). Each cell of the first table: share of 95% intervals for \(\beta_1\) that contain the true value. Monte Carlo error: ± 0.5 pp.

Quasi-likelihood is a good choice when:

1. The overdispersion is proportional to the nominal variance. In S1, Poisson covers 73.8% and quasi-Poisson 94.5%. This is the case for aggregated counts of homogeneous units, for proportions with a constant group size (28), and for any situation in which a plot of variance against mean (like the table in §6) shows a roughly constant ratio.
2. The interest is in \(\boldsymbol\beta\), not in the distribution. Effects, rate ratios, odds ratios, hypothesis tests: all of these need only the mean and honest standard errors.
3. There is no reasonable distribution. Continuous proportions, indices, counts with weak and unknown dependence — leaf blotch is the textbook case.
4. The sample is small and \(V\) is plausible. In S3, with \(n = 12\), quasi-Poisson with the \(t\) quantile covers 95.2%, while the sandwich drops to 85.8%: a single parameter \(\phi\) is far more stable than one variance per observation. Note the difference between the normal quantile (92.0%) and \(t\) (95.2%): with small \(n\), using \(t\) and \(F\) matters.
5. You want robustness without complication. The fit is that of the nominal GLM, and \(\hat{\boldsymbol\beta}\) stays consistent even if \(V\) is wrong (5.1).

§10When quasi models are bad

1. When the shape of the variance is wrong. In S2, using \(V = \mu\) for data with variance \(\propto\mu^2\) keeps \(\hat\beta_1\) unbiased, but the actual standard deviation of the estimator rises from 0.124 to 0.162 — a variance 1.7 times larger, which is like throwing away about 40% of the data. And the reported standard error (0.148 on average) understates the real one (0.162): coverage drops to 93.3%. Multiplying by \(\sqrt{\hat\phi}\) fixes the average size of the variance, not its shape. The same holds for counts with \(\operatorname{Var} = \mu + \mu^2/\theta\): quasi-Poisson gives weight \(\propto\mu\) precisely to the noisiest observations — the article From Poisson to the negative binomial shows its coverage falling to 88.6% in that case.
2. When you need probabilities. There is no distribution, so there is no \(\Pr(Y = 0)\), no predictive quantiles, no prediction intervals and no simulation of new data. To predict “how many policies will have zero claims”, you need a full model (negative binomial, beta-binomial).
3. When you need a likelihood. Without \(\ell\) there is no AIC, no BIC, no likelihood ratio between different families and no Bayesian inference. \(Q^+\) and QAIC are useful patches, not full substitutes.
4. With binary (0/1) data. For \(Y \in \{0, 1\}\), \(\operatorname{Var}(Y) = \mu(1-\mu)\) necessarily: there is no overdispersion in a Bernoulli, and \(\phi\) is not identifiable. We simulated 2,000 datasets of \(n = 200\) binary responses with strong unobserved heterogeneity (a \(\mathcal N(0,\ 1.5^2)\) random effect on the logit scale) and fitted quasi-binomial: \(\hat\phi\) averaged 1.009, with 95% of values between 0.979 and 1.046. The heterogeneity exists, but \(\hat\phi\) cannot see it. It only shows up with grouping (\(m_i \gt 1\)) or with a mixed model.
5. When the problem is the mean, not the variance. A forgotten covariate, a missing interaction or a nonlinear relationship inflate \(\hat\phi\). The quasi model “absorbs” the symptom with larger standard errors, but the coefficients stay biased. Get \(\mathbb E[Y]\) right first — it is the only assumption consistency depends on (5.1).
6. With dependence. Repeated measures on the same individual and clustered units violate the independence that (18) uses. The natural extension is GEE, which is quasi-likelihood with a working correlation matrix — or a mixed model.
7. With structure the mean and variance do not capture. Structural excess zeros, truncation, mixtures of populations: here the shape of the distribution is the phenomenon, and a model for it (zero-inflated, hurdle) is needed.

§11In practice: R and Python

# R — quasi-Poisson and quasi-binomial
qp <- glm(insects ~ spray, family = quasipoisson, data = d)
summary(qp)$dispersion                      # Pearson phi-hat (17)
anova(update(qp, . ~ 1), qp, test = "F")    # F test (20)

qb <- glm(cbind(germ, 30 - germ) ~ seed * extract, family = quasibinomial, data = g)

# Wedderburn's variance: a quasi family with its own V and deviance
wedd <- quasi(link = "logit", variance = "mu(1-mu)")
wedd$variance   <- function(mu) (mu * (1 - mu))^2
wedd$dev.resids <- function(y, mu, wt)
  wt * 2 * ((2*y - 1) * (qlogis(y) - qlogis(mu)) - 2 + y/mu + (1 - y)/(1 - mu))
fw <- glm(prop ~ site + variety, family = wedd, data = leaf)

import statsmodels.api as sm
import statsmodels.formula.api as smf

# quasi-Poisson: the Poisson fit with phi estimated by Pearson
qp = smf.glm("insects ~ spray", data=d, family=sm.families.Poisson()).fit(scale="X2")
print(qp.scale)                       # phi-hat (17)
sand = smf.glm("insects ~ spray", data=d, family=sm.families.Poisson()).fit(cov_type="HC0")  # (19)

# quasi-binomial
qb = smf.glm("germ + fail ~ seed * extract", data=g,          # two columns: successes and failures
             family=sm.families.Binomial()).fit(scale="X2")

For a \(V\) with no ready-made family, this article's script implements the IRLS (16), \(\hat\phi\) (17), the deviance (12), the sandwich (19) and \(Q^+\) (21) in about 60 lines of numpy.

§12Summary

1. The score equations of a GLM use only the mean \(\mu_i(\boldsymbol\beta)\) and the variance function \(V(\mu_i)\). Wedderburn (1974) drew the conclusion: specify only \(\mathbb E[Y] = \mu\) and \(\operatorname{Var}(Y) = \phi V(\mu)\).
2. The quasi-score \(U = (Y-\mu)/(\phi V(\mu))\) has mean zero, variance \(1/(\phi V)\) and satisfies the information identity — the three properties behind the asymptotic theory of likelihood. Its integral \(Q\) is the quasi-likelihood.
3. \(\hat{\boldsymbol\beta}\) solves \(\mathbf D^\top\mathbf V^{-1}(\mathbf y-\boldsymbol\mu) = \mathbf 0\) with the same IRLS as GLMs. \(\phi\) cancels; it is estimated afterwards by \(X^2/(n-p)\); the covariance is \(\hat\phi(\mathbf X^\top\mathbf W\mathbf X)^{-1}\). Tests use \(t\) and \(F\).
4. Advantages: consistency requiring only the right mean; separation of \(\boldsymbol\beta\) and \(\phi\); optimality among linear estimating equations (a generalised Gauss–Markov); the exponential family as a special case; any \(V\), even with no distribution, and \(\phi \lt 1\).
5. Examples: quasi-Poisson (\(\hat\phi = 3.34\), standard errors 1.8 times larger than Poisson's), quasi-binomial (type I error 22.6% → 4.4%) and Wedderburn's variance (patternless residuals, larger \(Q^+\), a more precise contrast).
6. They are good when the variance is proportional to the nominal one, the interest is in the coefficients, there is no natural distribution or the sample is small. They are bad when the shape of \(V\) is wrong (loss of efficiency and coverage), when probabilities or a likelihood are needed, with binary data, with dependence, or when the problem is in the mean.

§13Further reading

• R. W. M. Wedderburn (1974). Quasi-likelihood functions, generalized linear models, and the Gauss–Newton method. Biometrika, 61(3), 439–447. The original paper, short and readable.
• J. A. Nelder and R. W. M. Wedderburn (1972). Generalized linear models. Journal of the Royal Statistical Society A, 135(3), 370–384.
• P. McCullagh (1983). Quasi-likelihood functions. Annals of Statistics, 11(1), 59–67. The asymptotic theory.
• P. McCullagh and J. A. Nelder (1989). Generalized Linear Models, 2nd ed., ch. 9 and 10. Chapman & Hall. Includes the leaf-blotch data.
• J. A. Nelder and D. Pregibon (1987). An extended quasi-likelihood function. Biometrika, 74(2), 221–232.
• K.-Y. Liang and S. L. Zeger (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13–22. GEE.
• V. P. Godambe (1960). An optimum property of regular maximum likelihood estimation. Annals of Mathematical Statistics, 31(4), 1208–1211. The origin of optimal estimating-equation theory.
• D. A. Williams (1982). Extra-binomial variation in logistic linear models. Applied Statistics, 31(2), 144–148.
• J.-D. Lebreton, K. P. Burnham, J. Clobert and D. R. Anderson (1992). Modeling survival and testing biological hypotheses using marked animals. Ecological Monographs, 62(1), 67–118. Origin of QAIC.
• A. Agresti (2015). Foundations of Linear and Generalized Linear Models, ch. 8. Wiley.