Gamma regression

Topic of Underdispersion, positive responses, and quasi-likelihood, part of Generalised linear models · Also in Introduction to generalised linear models

Gamma regression

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

A GLM for positive continuous responses with a gamma random component, so that the standard deviation is proportional to the mean. The log link is the usual choice; the canonical inverse link can fail when the linear predictor goes negative.

In Underdispersion, positive responses, and quasi-likelihood

From Underdispersion, positive responses, and quasi-likelihood

A clotting time, or an insurance claim amount, is continuous and cannot be negative. Three links are in ordinary use: \begin{align*} \log\mu_i&={\boldsymbol{X}}_i^\top\boldsymbol{\beta}, & \mu_i&=\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta}),\\ \mu_i^{-1}&={\boldsymbol{X}}_i^\top\boldsymbol{\beta}, & \mu_i&=1/({\boldsymbol{X}}_i^\top\boldsymbol{\beta}),\\ \mu_i^{-2}&={\boldsymbol{X}}_i^\top\boldsymbol{\beta}, & \mu_i&=1/\sqrt{{\boldsymbol{X}}_i^\top\boldsymbol{\beta}}. \end{align*} The log link is the one with a direct reading, and it keeps \mu positive for every finite predictor. It is not the canonical link for the gamma or the inverse Gaussian, so R will not choose it unless you name it. The canonical link for the gamma is the inverse link. The canonical link for the inverse Gaussian is the inverse quadratic. Both of those are only positive when {\boldsymbol{X}}^\top\boldsymbol{\beta}>0. If the linear predictor changes sign, the inverse link produces a negative mean, and the square root in the inverse quadratic leaves the reals. It is sometimes said that the inverse quadratic “ensures” non-negativity because a complex mean is not a mean. That is a restriction on the domain of \boldsymbol{\beta}, not a guarantee for every coefficient vector. You still have to keep the linear predictor positive.

The two families are written by their mean and variance, so the shape-rate and shape-scale arguments do not have to be reconstructed: \begin{align*} Y_i\mid {\boldsymbol{X}}_i&\stackrel{\mathrm{ind}}{\sim}\mathrm{Gamma}\!\left(\mu_i=\mu({\boldsymbol{X}}_i^\top\boldsymbol{\beta}),\ \mathrm{Var}=\phi\mu_i^2\right),\\ Y_i\mid {\boldsymbol{X}}_i&\stackrel{\mathrm{ind}}{\sim}\mathrm{InverseGaussian}\!\left(\mu_i=\mu({\boldsymbol{X}}_i^\top\boldsymbol{\beta}),\ \mathrm{Var}=\phi\mu_i^3\right). \end{align*} The mean function \mu(\cdot) is whichever inverse link you chose. It is not tied to the family. The score for a gamma regression has the same shape as the Bernoulli and Poisson scores, (y-\mu)\mu'/\mathrm{Var}, and showing that is an exercise rather than a new idea. Intervals are again z intervals from the inverse Fisher information.

Clotting times, three ways

Main article: Clotting-time data

Return to the clotting data. Two lots of thromboplastin, log dose, clotting time in seconds. If a dose of zero is imagined to mean no clotting agent at all, an infinite mean time is the story the inverse link tells, and it is the default for family = Gamma. Fit an interaction, so each lot has its own intercept and its own slope:

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(clot, package = "faraway")
clot$Y <- clot$time
clot$logdose <- log(clot$conc)
clot$lot.f <- factor(clot$lot, levels = c("one", "two"), labels = c("1", "2"))
summary(glm(Y ~ lot.f * logdose, data = clot, family = Gamma))

                Estimate Std. Error t value Pr(>|t|)
(Intercept)     -0.01655   0.000866   -19.13  2.0e-11
lot.f2          -0.00735   0.001678    -4.38  0.00063
logdose          0.01534   0.000387    39.63  8.9e-16
lot.f2:logdose   0.00826   0.000735    11.23  2.2e-08

Dispersion \hat\phi=0.00213, three iterations, AIC 63.2. The first lot is the baseline factor level. Its mean is \mu_1=(-0.01655+0.01534\,\mathrm{logdose})^{-1}. The second lot adds both coefficients. The intercept becomes -0.01655-0.00735=-0.0239, and the slope becomes 0.01534+0.00826=0.0236, so \mu_2=(-0.0239+0.0236\,\mathrm{logdose})^{-1}. On the inverse scale the second lot starts lower and climbs faster. Because the mean is the reciprocal, that is a higher curve which then falls faster as dose increases. The interaction line has t=11.2 and a p-value of 2\times 10^{-8}. The slopes are not parallel. The dispersion is tiny: these points sit close to the curve, as a controlled clotting experiment should.

The default inverse-Gaussian fit uses the inverse-quadratic link. The coefficients are on a different scale, of order 10^{-3} rather than 10^{-2}, because they sit inside a square root. The AIC is 110.5, much worse than 63.

Figure 1: Clotting time against concentration, for the two lots. Solid curves are the gamma model with the inverse link. Dotted curves are the inverse Gaussian with the inverse-quadratic link. The vertical scale stops at 140 seconds; the dotted curves run off the top at low concentration and then flatten below the points.

In Figure the inverse-link gamma curves pass through the points, and the inverse-quadratic curves compress too quickly and then flatten too soon.

Now force the inverse Gaussian to use the inverse link, the same mean function as the gamma:

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(clot, package = "faraway")
clot$Y <- clot$time
clot$logdose <- log(clot$conc)
clot$lot.f <- factor(clot$lot, levels = c("one", "two"), labels = c("1", "2"))
summary(glm(Y ~ lot.f * logdose, data = clot,
            family = inverse.gaussian(link = "inverse")))

The coefficients come back to -0.0178, -0.0074, 0.0158, 0.0083, close to the gamma fit, with AIC 60.7. The two mean curves lie on top of each other. The standard errors move, because the variance function changed from \phi\mu^2 to \phi\mu^3.

Key idea

The link chooses the mean curve. The response distribution mainly chooses the standard errors. Those two pieces are essentially orthogonal: you may pair a scientifically sensible link with a distribution that matches the scatter, and neither choice is doing the other’s job.

That separation is what Huang and Rathouz (2017) give in general. It is why families and links can be mixed, and why a worse AIC for the inverse quadratic is a statement about the link, not a mystery of the inverse Gaussian. Other positive families exist and are not in the default glm list: inverse gamma, lognormal, Weibull with the second parameter fixed. A custom family needs a log-likelihood, or, as the next section shows, sometimes only a mean and a variance.

In Introduction to generalised linear models

From Introduction to generalised linear models

McCullagh and Nelder discuss an experiment on the clotting time of blood. Normal plasma was diluted to nine percentage concentrations and clotting was induced by two lots of thromboplastin. The response is a time in seconds, somewhere from about twenty seconds up to two minutes. The doses on the horizontal axis already sit on a doubling scale: 5, 10, 20, 50, 100.

Log dose, rather than dose, is the covariate to use. A doctor who raises a dose by 5\,\mathrm{mg} is doubling a 5\,\mathrm{mg} prescription and barely touching a 100\,\mathrm{mg} one. What matters is the multiplicative change. Whether the logarithm is base 2 or base e only rescales the slope: a unit of \log_2 is a doubling, and a unit of \log_e is a factor of about 2.7. The scientific point is the same.

A straight line for the mean clotting time will eventually become negative, and a waiting time cannot. Two links keep the mean positive. The log link is Equation (3) in Introduction to generalised linear models again. The inverse, or hyperbolic, link is \mu=\frac{1}{{\boldsymbol{X}}^\top\boldsymbol{\beta}} \qquad\text{equivalently}\qquad {\boldsymbol{X}}^\top\boldsymbol{\beta}=\frac{1}{\mu}, provided {\boldsymbol{X}}^\top\boldsymbol{\beta}>0. If the linear predictor is allowed to change sign, this expression goes negative and the model has broken.

The two links agree that a large dose of clotting agent can drive the mean time down towards zero. They disagree at the other end. As the linear predictor in Equation (4) in Introduction to generalised linear models approaches zero, the mean time goes to infinity. The log link, for any finite predictor, stays finite. So the modelling question is biological, not a default in the software. If there were truly no plasma, would the blood never clot, or would it clot after a long but finite time? Nobody in the experiment has zero plasma. A person with very little plasma does eventually clot, though a haemophiliac may bleed out first. That leans towards a finite mean, and so towards the log link. The inverse link is the honest model of the hypothetical with no plasma at all. Either is allowed. They can be paired later with whatever response distribution we choose. The link is not dictated by the distribution.

The times themselves are continuous and strictly positive, so they are not normal. A gamma distribution is a reasonable start: it is supported on (0,\infty) and it extends the exponential. A gamma law has a second parameter, a dispersion, which plays the role that \sigma^2 plays for the normal and can be chosen separately from the mean. In the shape-and-rate parameterisation, a \mathrm{Gamma}(\alpha,\alpha/\mu) random variable has mean \mu and variance \mu^2/\alpha. The model is the mean and the variance: Y_i\mid {\boldsymbol{X}}_i\stackrel{\mathrm{ind}}{\sim} \mathrm{Gamma}\!\left(\mu_i=\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta}),\ \mathrm{Var}=\phi\mu_i^2\right). Here \phi is the dispersion. The variance grows with the square of the mean, which is a stronger mean–variance relationship than the Poisson’s \mathrm{Var}=\mu. It is the same shape you see in the data: times near zero have little room to vary, and large mean times can vary a lot. Written out in full, this is shape 1/\phi and rate 1/(\phi\mu), because then the mean is \mu and the variance is \phi\mu^2.

Used in

Taught in Underdispersion, positive responses, and quasi-likelihood · Introduction to generalised linear models
Questions that use it (1) — Q55 (answer)
All 6 questions on Underdispersion, positive responses, and quasi-likelihood

In Underdispersion, positive responses, and quasi-likelihood

Clotting times, three ways

In Introduction to generalised linear models

Used in

See also