Three examples the linear model does not fit
Section of Introduction to generalised linear models, part of Generalised linear models
Three examples the linear model does not fit
| Type | Section |
| Taught in | Introduction to generalised linear models |
| Part | Generalised linear models |
Runs in the past week
Main article: Runs in the past week
A hundred people at the gym were asked how many times they had run in the past week, and that count was plotted against age. The responses are integers, 0,1,2,\ldots, with plenty of ties. Jitter in a scatterplot is only there so that those ties do not sit on top of each other. A person cannot have run 1.37 times in a given week.
The mean is the first problem. A straight line through this cloud may not look absurd between the ages we actually saw. Extend it. By age 50 or 60 the line has crossed below zero, and a mean number of runs cannot be negative. Any function that is allowed to go to \pm\infty will do that for some covariate, whether or not that covariate appeared in the sample. A model has to stay valid for data we could have seen.
The repair is to put the straight line on the log of the mean, not on the mean: \mu=\exp({\boldsymbol{X}}^\top\boldsymbol{\beta}) \qquad\text{equivalently}\qquad {\boldsymbol{X}}^\top\boldsymbol{\beta}=\log\mu, where \mu=\mathrm{E}(Y\mid {\boldsymbol{X}}). The exponential of anything real is positive, so the mean cannot cross zero. The transformation that takes the mean across to the linear predictor is the link. The map the other way, from {\boldsymbol{X}}^\top\boldsymbol{\beta} back to \mu, is the inverse link. We will often interpret on the mean scale. The literature names the link in the other direction, and we will follow the literature so that we can talk to it.
Do not take logarithms of the counts. Zero runs is a common and perfectly sensible observation, and \log 0 is not a number. Logging the observations is not what a generalised linear model does. The data stay on the original scale. Only the mean curve moves.
The distribution is the next problem. Integer counts with ties are not normal draws. A reasonable starting model is Y_i\mid {\boldsymbol{X}}_i\stackrel{\mathrm{ind}}{\sim} \mathrm{Poisson}\!\left(\mu_i=\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta})\right). There are other distributions for counts. The Poisson is the default, not the only choice.
The variance comes with that choice, and it is not constant. For a Poisson random variable the variance equals the mean, \mathrm{Var}(Y_i\mid {\boldsymbol{X}}_i)=\mu_i=\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta}). That is what the scatterplot is already saying. Near a mean of zero the counts pile up on 0 and 1, because they cannot go negative, so the spread is small. Where the mean is larger, the spread is larger. A Poisson model does not have to be told this. It is the variance formula.
Clotting times
Main article: Clotting-time data
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.
Disease and cholesterol
Main article: Cholesterol and heart disease
A hundred patients have a cholesterol measurement and a binary disease indicator, with sex marked on the plot. Vertical jitter stops a pile of zeroes, or a pile of ones, looking like a fractional response. The response really is 0 or 1. As cholesterol increases, particularly once it is above about 6.5, disease becomes the typical outcome. Down near 4.5, absence is more common.
The mean of a binary response is the success probability, p=\mathrm{E}(Y\mid {\boldsymbol{X}})=\mathrm{P}(Y=1\mid {\boldsymbol{X}}). A straight line cannot represent it. The line will eventually leave (0,1), and a probability cannot. The exponential \exp({\boldsymbol{X}}^\top\boldsymbol{\beta}) stays positive, but it is not bounded by 1. Dividing by one plus itself fixes that: p=\frac{\exp({\boldsymbol{X}}^\top\boldsymbol{\beta})}{1+\exp({\boldsymbol{X}}^\top\boldsymbol{\beta})} \qquad\text{equivalently}\qquad {\boldsymbol{X}}^\top\boldsymbol{\beta}=\log\frac{p}{1-p}. The right-hand equation is the logistic link, also called the logit: the log of the odds p/(1-p). The left-hand equation is the inverse link, sometimes called the expit. Log-odds run over the whole real line, which is what a linear predictor needs. The probability stays in (0,1).
Again, the link is applied to the mean. The logistic function of a zero or a one is not defined in any useful way, and any function of a binary observation still takes at most two values. You cannot manufacture a curve by transforming zeroes and ones.
The distribution is Bernoulli, a single coin toss, not a normal: Y_i\mid {\boldsymbol{X}}_i\stackrel{\mathrm{ind}}{\sim} \mathrm{Bernoulli}\!\left(p_i=\frac{\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta})}{1+\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta})}\right). The variance of a coin toss is p(1-p), so \mathrm{Var}(Y_i\mid {\boldsymbol{X}}_i)=p_i(1-p_i)=\frac{\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta})}{\bigl(1+\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta})\bigr)^2}. This variance is not monotone in p. It is largest at p=1/2 and it shrinks towards 0 as p approaches 0 or 1. If almost everyone is diseased, there is little disagreement left in the sample. The same is true if almost nobody is.
In all three examples the mean is a function of a linear predictor, and the observations are not transformed. The link is chosen so that the mean stays inside the set of values the response can actually have. The distribution is chosen to match the kind of observation: counts, positive times, or zeroes and ones. The variance formula is then a consequence, not a separate constant \sigma^2.
Put the four models next to each other: \begin{align*} Y_i\mid {\boldsymbol{X}}_i &\stackrel{\mathrm{ind}}{\sim} \mathrm{N}({\boldsymbol{X}}_i^\top\boldsymbol{\beta},\sigma^2),\\ Y_i\mid {\boldsymbol{X}}_i &\stackrel{\mathrm{ind}}{\sim} \mathrm{Poisson}\!\left(\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta})\right),\\ Y_i\mid {\boldsymbol{X}}_i &\stackrel{\mathrm{ind}}{\sim} \mathrm{Gamma}\!\left(\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta}),\ \mathrm{Var}=\phi\mu_i^2\right),\\ Y_i\mid {\boldsymbol{X}}_i &\stackrel{\mathrm{ind}}{\sim} \mathrm{Bernoulli}\!\left(\frac{\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta})}{1+\exp({\boldsymbol{X}}_i^\top\boldsymbol{\beta})}\right). \end{align*} Before Nelder’s work in the early 1970s, each new distribution arrived with its own theory and its own fitting algorithm. These four, and others like them, are exponential families. One piece of theory covers them, and one numerical algorithm fits them.
Used in
See also
Exponential family · Mean and variance of an exponential family · Generalised linear model · Canonical link · What a regression statement is hiding · Runs in the past week · Clotting-time data · Cholesterol and heart disease