Variance function

Topic of Underdispersion, positive responses, and quasi-likelihood, part of Generalised linear models

Variance function

Type Topic
Taught in Underdispersion, positive responses, and quasi-likelihood
Part Generalised linear models
Questions using it 1

The function V(\mu) in \operatorname{Var}(Y)=\phi V(\mu): 1 for the normal, \mu for the Poisson, \mu(1-\mu) for the Bernoulli, \mu^2 for the gamma, \mu^3 for the inverse Gaussian. Quasi-likelihood needs only a mean and a variance function, not a full distribution.

In Underdispersion, positive responses, and quasi-likelihood

From Underdispersion, positive responses, and quasi-likelihood

The normal, the gamma, and the inverse Gaussian already have a dispersion \phi. The Poisson and the binomial do not: the variance is a function of the mean alone. Putting a free \phi in front, \mathrm{Var}(Y)=\phi\mu \qquad\text{or}\qquad \mathrm{Var}(Y)=\phi\, np(1-p), is the quasi-Poisson and quasi-binomial we have already used for dead time, spurts, and pseudo-replication. In general, any non-negative function V gives a variance \phi V(\mu). Most of those choices are not the variance of an exponential family. They describe the data. They do not, by themselves, prescribe a distribution you can simulate from. That is the whole of the name quasi.

Definition 1 If \mathrm{Var}(Y)=\phi V(\mu), the quasi-score in the mean is S_\mu=\frac{Y-\mu}{\phi V(\mu)}. If the mean depends on covariates through \mu({\boldsymbol{X}}^\top\boldsymbol{\beta}), the quasi-score in \boldsymbol{\beta} is {\boldsymbol{S}}_{\boldsymbol{\beta}} =\frac{Y-\mu}{\phi V(\mu)}\,\mu'({\boldsymbol{X}}^\top\boldsymbol{\beta})\,{\boldsymbol{X}}.

Article: Quasi-score and quasi-likelihood →

Two properties make S_\mu behave as a score. The expectation is zero, because \mathrm{E}(Y-\mu)=0 and the rest does not depend on Y. And the variance equals the negative expected derivative: \mathrm{Var}(S_\mu)=\frac{1}{\phi^2 V(\mu)^2}\mathrm{Var}(Y)=\frac{1}{\phi V(\mu)} =-\mathrm{E}\!\left(\frac{\partial S_\mu}{\partial\mu}\right). Every GLM score has the shape of {\boldsymbol{S}}_{\boldsymbol{\beta}}. Not every choice of V comes from a GLM. Inside an exponential family the variance function determines the family, so V(\mu)=\mu is the Poisson and V(\mu)=np(1-p) is the binomial. A free \phi takes you outside that class. You can still invent a probability model with the same moments, for example a normal with variance \phi\mu, but that normal is not a distribution on the counts, and it is not an exponential family in \mu because \mu appears in both the mean and the variance. Fitting it needs its own algorithm. The quasi-score does not: glm only solves {\boldsymbol{S}}_{\boldsymbol{\beta}}={\boldsymbol{0}}.

The quasi-log-likelihood is the integral of that score, Q(\mu;y)=\int_y^{\mu}\frac{y-t}{\phi V(t)}\,\mathrm{d}t, and the sample version is the sum of Q(\mu_i;y_i). The integral is rarely needed. Estimation and standard errors use the score and its variance. The integral is what you need if you want something to play the role of a likelihood in an AIC comparison. Because \phi factors out of {\boldsymbol{S}}_{\boldsymbol{\beta}}, the estimate of \boldsymbol{\beta} does not depend on \phi. Estimate the mean first. Then estimate the dispersion by moments, exactly as \sigma^2 is estimated after a normal line: \hat\phi=\frac{1}{n-p}\sum_{i=1}^n\frac{(y_i-\hat\mu_i)^2}{V(\hat\mu_i)}. The covariance of \hat{\boldsymbol{\beta}} is \hat\phi times the inverse of the weight matrix built from (\mu')^2/V(\mu).

Leaf blotch

Main article: Leaf-blotch data

Wedderburn (1974) is the paper that introduced this. Ten barley varieties were grown at nine sites, one proportion per cell: the percentage of leaf area covered by blotch, Rhynchosporium secalis. The variety and site numbers in the table are labels, ordered roughly from less blotch to more. They are not numeric covariates. Treating “variety 2” as twice “variety 1” is the same mistake as averaging student ID numbers.

There is one observation per variety–site pair, so a full interaction cannot be fitted. The mean model is main effects only, on the logit of the proportion p_{ij}: \log\frac{p_{ij}}{1-p_{ij}}=\alpha_i+\beta_j, \qquad p_{ij}=\frac{\exp(\alpha_i+\beta_j)}{1+\exp(\alpha_i+\beta_j)}. A structured interaction would be allowed if a reason for that particular combination existed before looking at the table. It does not. The variance starts as quasi-binomial, \phi p(1-p). The proportions are averages of many leaves from the same plot, so the trials are not independent Bernoulli draws. Pseudo-replication makes \phi much smaller than 1.

need <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg, repos = "https://cloud.r-project.org")
  }
  suppressPackageStartupMessages(library(pkg, character.only = TRUE))
}
need("faraway")
data(leafblotch, package = "faraway")
leafblotch$proportions <- leafblotch$blotch
summary(glm(proportions ~ factor(site) + factor(variety),
            data = leafblotch, family = quasibinomial))

The dispersion printed for this fit is \hat\phi=0.089. The intercept, -8.05, is the log-odds for variety 1 at site 1. Every other coefficient is a contrast with that baseline. After adjusting for site, varieties 2, 3, and 4 are not distinguishable from variety 1 at the 5% level (variety 2 even has t=0.21). Varieties 5 through 10 are worse. Variety 1 is the most resilient in this fit, with 2, 3, and 4 in the same neighbourhood.

The Pearson residual is (y_i-\hat\mu_i)/\sqrt{\hat\phi V(\hat\mu_i)}. Plotted against the linear predictor, these residuals are not a flat band. The quasi-binomial variance compresses the small fitted values too hard and does not compress the middle enough. Wedderburn’s replacement is V(\mu)=\mu^2(1-\mu)^2, still multiplied by a free \phi. Squaring p(1-p) shrinks the variance even more near 0 and 1, and it leaves relatively more variance in the middle, which is what the residual plot was asking for.

Figure 1: Variance as a function of the mean. On the left, mu, mu squared and mu cubed are the Poisson, gamma and inverse-Gaussian choices. On the right, the binomial variance and Wedderburn’s replacement: smaller in the tails, and relatively larger through the middle.

No exponential family has this variance function. It is a description. It is not a quasi-binomial. There is nothing to call it except a quasi-likelihood with this V.

It is not a built-in family. You supply V and its derivatives to glm. A proportion of exactly zero makes the quasi-deviance integral diverge, so the block below prints the variety contrasts and their standard errors, which do not need that integral.

need <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg, repos = "https://cloud.r-project.org")
  }
  suppressPackageStartupMessages(library(pkg, character.only = TRUE))
}
need("faraway")
data(leafblotch, package = "faraway")
leafblotch$proportions <- leafblotch$blotch
wedderburn <- quasi(
  link = "logit",
  variance = list(
    varfun = function(mu) mu^2 * (1 - mu)^2,
    validmu = function(mu) all(mu > 0 & mu < 1),
    dev.resids = function(y, mu, wt) wt * (y - mu)^2 / (mu^2 * (1 - mu)^2),
    initialize = expression({
      n <- rep.int(1, nobs)
      mustart <- pmax(0.001, pmin(0.999, y))
    }),
    name = "mu^2(1-mu)^2"
  )
)
fit <- glm(proportions ~ factor(site) + factor(variety),
           data = leafblotch, family = wedderburn)
co <- coef(summary(fit))
round(co[grep("variety", rownames(co)), ], 3)

After that fit, the contrasts with variety 1, and a common standard error of 0.47 because the design is balanced, are

variety 1 2 3 4 5 6 7 8 9 10
contrast 0 -0.47 0.08 0.95 1.35 1.33 2.34 3.26 3.14 3.89

Variety 2 is now the most resilient, slightly ahead of variety 1. Varieties 1, 2, and 3 remain the resistant group: variety 4 is about two standard errors up, and from variety 5 the gap is clear. The ranking moved a little. The standard errors moved more, from about 0.6 down to 0.47, and the residual plot is closer to a constant-width scatter.

Figure 2: Pearson residuals against the linear predictor. On the left the quasi-binomial variance squeezes the small fitted values, at the left of the panel, onto the horizontal axis. On the right, Wedderburn’s variance gives those same points more room. Dashed lines are at plus or minus 2.

For a proportion that is actually continuous, a Pearson plot is a fair diagnostic. For a binary or count response the same plot bands, because the observation can only sit on a lattice, and that banding is not evidence against the mean. We will handle discrete residuals later.

Key idea

Once the mean is specified, any non-negative variance function V(\mu) gives a quasi-score that glm can solve. Whether that V is the variance of a named distribution is a separate question. A descriptive fit is allowed. It is not a generator of new data.

  1. From the gamma interaction fit, recompute the second lot’s intercept and slope by adding the baseline coefficients to the lot.f2 coefficients, and write \mu_2 as a reciprocal.

  2. Show \mathrm{E}(S_\mu)=0 and \mathrm{Var}(S_\mu)=1/(\phi V(\mu)) from the definition of S_\mu, using only \mathrm{E}(Y)=\mu and \mathrm{Var}(Y)=\phi V(\mu).

  3. Explain why \hat{\boldsymbol{\beta}} in a quasi-model does not depend on \phi, and why \hat\phi in Equation (1) in Underdispersion, positive responses, and quasi-likelihood is computed only after \hat{\boldsymbol{\beta}}.

  4. The leaf-blotch design has one observation per variety–site cell. Say why a product interaction site * variety is not available, and why that is different from refusing the interaction on scientific grounds.

Used in

See also

Quasi-score and quasi-likelihood · Quasi-likelihood · Inverse Gaussian distribution · Variance · Bids for a bankrupt firm · Positive continuous responses · Leaf-blotch data · Log link

Glossary and index