5 Generalized linear models
This chapter treats a generalization of the linear model that preserves much of its tractability, but can be adapted to cases where the linear model is inappropriate. This generalization is tightly linked to the exponential dispersion distributions treated in Chapter 4.
Section 5.1 introduces the fundamental models assumptions as generalizations of the model assumptions for the linear model. Section 5.2 covers how generalized linear models can be fitted to data via maximum-likelihood estimation. Section 5.3 covers the standard constructions of confidence intervals and test statistics and their distributional approximations that are used for statistical inference. The content of these sections is an extension of similar methods and results for the linear model as treated in Chapter 2.
5.1 Model assumptions
To motivate the general model assumptions we will consider three examples. In all three examples \(X\) denotes a vector of predictor variables, and we introduce a model of the outcome given the predictors via a linear combination \(X^T \beta\) of the predictors for some parameter vector \(\beta.\)
Example 5.1 For a binary outcome \(Y \in \{0, 1\}\) we know that \[ \mathbf{E}(Y \mid X) = \mathbf{P}(Y = 1 \mid X) \in [0, 1] \] and \[ \mathbf{V}(Y \mid X) = \mathbf{P}(Y = 1 \mid X)(1 - \mathbf{P}(Y = 1 \mid X)) = \mathcal{V}(\mathbf{E}(Y \mid X)) \] for the variance function \(\mathcal{V}(\mu) = \mu(1 - \mu).\)
A linear model, \(\mathbf{P}(Y = 1 \mid X) = X^T \beta,\) of the conditional expectation is problematic. The constraint that a probability must be in \([0,1]\) can only be enforced by restricting the parameter space of \(\beta.\) The restriction will depend on the values of the observed predictors, and it is generally impossible to ensure that all future predictions will be in \([0,1].\) Moreover, the variance is not constant but depends upon the expectation via the variance function.
One solution is to consider a model of the form \[ \mathbf{P}(Y = 1 \mid X) = \mu(X^T \beta) \] for \(\mu : \mathbb{R} \to [0, 1]\) a given mean value function1. One possible choice is the logistic function, which corresponds to the canonical link for the Bernoulli distribution, and which results in the logistic regression model. Thus the model of the expectation is given in terms of the linear combination \(X^T \beta,\) but it is (nonlinearly) transformed using \(\mu.\) By assumption, \(\mu\) takes values in \([0, 1].\) The variance is then \(\mathcal{V}(\mu(X^T \beta)).\)
1 Typically a continuous and monotone function—a distribution function for instance.
Example 5.2 If \(Y\) is Poisson distributed, and thus takes values in \(\mathbb{N}_0,\) we know that \[ \mathbf{E}(Y \mid X) \geq 0, \] and \[ \mathbf{V}(Y \mid X) = \mathbf{E}(Y \mid X) = \mathcal{V}(\mathbf{E}(Y \mid X)) \] for the variance function \(\mathcal{V}(\mu) = \mu.\) As for the Bernoulli model, a linear model of the mean can produce values outside the distribution’s range of possible mean values. One solution is the log-linear model \[\mathbf{E}(Y \mid X) = e^{X^T \beta}.\] That is, the mean value function is the exponential function, and the logarithm of the mean is linear in the parameters.
To accommodate models such as the log-linear Poisson model and the logistic Bernoulli model, we generalize the linear model assumptions A1–A3 by allowing for nonlinear relations between a linear combination of the predictors and the mean and variance. This also allows for a non-constant variance as a function of the mean. For exponential dispersion distributions, such a relation is a mathematical necessity as the two examples above show.
The assumptions GA1–GA3 below, that replace A1–A3 for generalized linear models, follow as abstractions of the structure found for exponential dispersion models. To ease notation we introduce \[ \eta = X^T \beta \] to denote the linear predictor for a vector of predictors \(X\) and a vector of parameters \(\beta.\)
GA1. The conditional expectation of \(Y\) given \(X\) is \[ \mathbf{E}(Y \mid X) = \mu(\eta), \] with2 \(\mu : \mathbb{R} \to J\) the mean value function. The set \(J \subseteq \mathbb{R}\) denotes the range of the mean value function.
2 Occasionally, models have a mean value function defined on a domain \(H \subseteq \mathbb{R}\) that is not the entire real line. This implicitly restricts the \(\beta\) parameter space to ensure that the linear predictor ends up in \(H\), which is a practical and theoretical nuisance. We will develop the theory assuming that \(H = \mathbb{R}\).
GA2. The conditional variance of \(Y\) given \(X\) is \[ \mathbf{V}(Y \mid X) = \varphi \mathcal{V}(\mu(\eta)), \] with \(\mathcal{V} : J \to (0, \infty)\) the variance function and \(\varphi > 0\) the dispersion parameter.
GA3. The conditional distribution of \(Y\) given \(X\) is the \((\theta(\eta), \nu_{\varphi})\)-exponential dispersion distribution, \[ Y \mid X \sim \mathcal{E}(\theta(\eta), \nu_{\varphi}). \]
The exponential dispersion distributions referred to in GA3 are introduced in Chapter 4. The normal distribution, as in Assumption A3, is a special case. As for the linear model, the assumption GA3 is a strong distributional assumption, which implies GA1 and GA2 for specific choices of \(\mu\) and \(\mathcal{V}\) that depend on the map \(\eta \mapsto \theta(\eta)\) and the chosen exponential dispersion model.
Whenever the mean value function \(\mu\) is one-to-one, its inverse \[ g \coloneqq \mu^{-1} : J \to \mathbb{R} \] is called the link function. Besides the identity link, common link functions are the log-link and logit-link functions corresponding to the exponential and logistic mean value functions, respectively. These link functions are also the canonical link functions for the Poisson model and the Bernoulli model, respectively.
Example 5.3 If \(Y\) is a positive random variable given by \[ \log(Y) = X^T\beta + \varepsilon \] for \(\varepsilon \sim \mathcal{N}(0,\sigma^2)\) then \[ \mathbf{E}(Y \mid X) = e^{X^T \beta} \mathbf{E}(e^{\varepsilon}). \] The distribution of \(e^{\varepsilon}\) is a log-normal distribution with mean \(e^{\sigma^2/2},\) hence \[ \mathbf{E}(Y \mid X) = e^{X^T \beta + \sigma^2/2}, \] which shows that the model is a log-linear model. The variance is \[ \begin{align*} \mathbf{V}(Y \mid X) & = e^{2 X^T \beta} \mathbf{V}(e^{\varepsilon}) = e^{2 X^T \beta} (e^{\sigma^2} - 1) e^{\sigma^2} \\ & = (e^{\sigma^2} - 1)(\mathbf{E}(Y \mid X))^2. \end{align*} \] Thus by introducing the dispersion parameter \(\varphi = e^{\sigma^2} - 1\) and the variance function \(\mathcal{V}(\mu) = \mu^2\) we see that \[ \mathbf{V}(Y \mid X) = \varphi \mathcal{V}(\mathbf{E}(Y \mid X)). \] In conclusion, the log-normal model is a log-linear model with a quadratic variance function. It follows from Example 4.4 that the exponential dispersion model with quadratic variance function is the Gamma model. Thus an alternative to the log-normal model with the same mean and variance structure is a Gamma model with the log-link.
To specify a generalized linear model, we can specify the exponential dispersion distribution, and through the strong distributional assumption GA3 we implicitly specify the mean value function via the parametrization \(\eta \mapsto \theta(\eta)\) and the variance function. We can also just through GA1 and GA2 specify the mean value function and the variance function directly, but the situation is a little more complicated than for the linear model. With arbitrary choices of \(\mu,\) \(\mathcal{V}\) and \(\varphi\) we cannot always find a corresponding exponential dispersion model. Moreover, Example 5.3 shows that there exist models fulfilling GA1 and GA2, which are not exponential dispersion models. For the log-normal model there is an exponential dispersion model with an equivalent mean and variance structure, but the outcome distribution is different. In this light it is noteworthy that the methodology developed does not hinge on GA3 to any great extent. Generalized linear models is largely a semiparametric modeling framework of mean and variance.
5.2 Estimation theory
In this section we cover the theory behind estimation in generalized linear models. We introduce maximum likelihood estimation for generalized linear models based on the exponential dispersion distributions. This includes the derivation of the nonlinear estimating equation, known as the score equation, and the iterative weighted least squares (IWLS) algorithm that is used in practice to fit the models to data.
5.2.1 Maximum likelihood estimation
We first consider the simple case where \(Y \sim \mathcal{E}(\theta(\eta), \nu_{\varphi}),\) that is, the distribution of \(Y\) is given by the exponential dispersion model with canonical parameter \(\theta(\eta)\) and structure measure \(\nu_{\varphi}.\) Derivations of the score equation and Fisher information in this case can then be used to derive the score equation and Fisher information in the general case when we have observations \(Y_1, \ldots, Y_n\) that are conditionally independent given the predictors with \(Y_i \mid \mathbf{X} \sim \mathcal{E}(\theta(\eta_i), \nu_{\varphi})\) and \(\eta_i = X_i^T \beta.\)
Definition 5.1 With \(\ell\) denoting a generic \(C^2\) log-likelihood function the score statistic is \[ U(\eta) \coloneqq \nabla_{\eta} \ell(\eta). \] The Fisher information, \[ \mathcal{J}(\eta) \coloneqq - \mathbf{E}(D_{\eta} U(\eta)), \] is minus the expectation of the Jacobian of the score statistic, or, equivalently, the expectation of the Hessian of the negative log-likelihood.
The score equation is obtained by equating the score statistic to \(0.\) In all that follows, the dispersion parameter \(\varphi\) is regarded as fixed.
Lemma 5.1 If \(Y \sim \mathcal{E}(\theta(\eta), \nu_{\varphi})\) then the log-likelihood function is \[ \ell_Y(\eta) = \frac{\theta(\eta) Y - \zeta(\eta)}{\varphi}, \] the score statistic is \(U(\eta) = \theta'(\eta)( Y - \mu(\eta)) / \varphi,\) and the Fisher information is \[ \mathcal{J}(\eta) = \frac{\theta'(\eta) \mu'(\eta)}{\varphi}. \]
Proof. The density for the distribution of \(Y\) w.r.t. \(\nu_{\varphi}\) is by definition \[ e^{\frac{\theta(\eta)y - \zeta(\eta)}{\varphi}}, \] and it follows that \(\ell_Y(\eta)\) has the stated form. Differentiation of \(\varphi \ell_Y(\eta)\) yields \[ \varphi U(\eta) = \varphi \ell'(\eta) = \theta'(\eta) Y - \zeta'(\eta) = \theta'(\eta)\Big( Y - \underbrace{\frac{\zeta'(\eta)}{\theta'(\eta)}}_{\mu(\eta)}\Big), \] where we have used Corollary 4.1. Furthermore, we find that \[ \varphi U'(\eta) = \theta''(\eta) (Y - \mu(\eta)) - \theta'(\eta) \mu'(\eta), \] and since \(\mathbf{E}_{\eta} (Y) = \mu(\eta)\) it follows that \[ \mathcal{J}(\eta) = - \mathbf{E}_{\eta}(U'(\eta)) = \frac{\theta'(\eta) \mu'(\eta)}{\varphi}. \]
With only a single observation, the score equation is equivalent to \(\mu(\eta) = Y,\) and it follows that there is a solution to the score equation if \(Y \in J = \mu(I).\) However, the situation with a single observation is not relevant for practical purposes. The result is only given as an intermediate step towards the next result.
The score statistic \(U\) above is a function of the univariate parameter \(\eta.\) We adapt in the following the convention that for a vector \(\boldsymbol{\eta} = (\eta_1, \ldots, \eta_n)^T\) \[ U(\boldsymbol{\eta}) = (U(\eta_1), \ldots, U(\eta_n))^T. \] That is, as a map the score statistic is applied coordinatewise to the vector \(\boldsymbol{\eta}.\) Note that the derivative (the Jacobian) of \(\boldsymbol{\eta} \mapsto U(\boldsymbol{\eta})\) is an \(n \times n\) diagonal matrix.
Theorem 5.1 Assume that \(Y_1, \ldots, Y_n\) are conditionally independent given \(\mathbf{X}\) and that \(Y_i \mid \mathbf{X} \sim \mathcal{E}(\theta(\eta_i), \nu_{\varphi})\) where \(\eta_i = X_i^T \beta.\) Then with \(\boldsymbol{\eta} = \mathbf{X}\beta\) the score statistic expressed in the \(\beta\)-parameter is \[ \mathcal{U}(\beta) = \mathbf{X}^T U(\boldsymbol{\eta}). \] The Fisher information is
\[ \mathcal{J}(\beta) = \mathbf{X}^T \mathbf{W} \mathbf{X}, \] with the entries in the diagonal weight matrix \(\mathbf{W}\) being \[ w_{ii} = \frac{(\mu_i')^2}{\varphi \mathcal{V}(\mu_i)} = \frac{\theta'_i\mu'_i}{\varphi} \] where \(\mu_i = \mu(\eta_i),\) \(\mu'_i = \mu'(\eta_i)\) and \(\theta'_i = \theta'(\eta_i).\)
Proof. By the independence assumption the log-likelihood is \[ \ell_{\mathbf{Y}} (\beta) = \sum_{i=1}^n \ell_{Y_i}(\eta_i) \] where \(\eta_i = X_i^T \beta.\) By the chain rule, \[ \mathcal{U}(\beta) = \nabla_{\beta} \ell_{\mathbf{Y}} (\beta) = \sum_{i=1}^n X_i U(\eta_i) = \mathbf{X}^T U(\boldsymbol{\eta}). \]
As argued above, \(- D_{\boldsymbol{\eta}} U(\boldsymbol{\eta})\) is diagonal, and the expectation of the diagonal entries are according to Lemma 5.1 \[w_{ii} = \frac{\theta'(\eta_i) \mu'(\eta_i)}{\varphi}.\] The alternative formula for the weights follows from Corollary 4.1. Thus \[ \mathcal{J}(\beta) = - \mathbf{E} (D_{\beta} \mathcal{U}(\beta) \mid \mathbf{X}) = - \mathbf{X}^T \mathbf{E} ( D_{\boldsymbol{\eta}} U(\boldsymbol{\eta}) \mid \mathbf{X}) \mathbf{X} = \mathbf{X}^T \mathbf{W} \mathbf{X}. \]
The score equation is the equation \[ \mathcal{U}(\beta) = 0, \] and if the maximum-likelihood estimator exists in \(\mathbb{R}^p\) it solves this equation. Theorem 5.1 shows that the score equation can be written as \[ \mathbf{X}^T U(\boldsymbol{\eta}) = \mathbf{X}^T U(\mathbf{X} \beta) = 0. \] By observing that \[ U(\boldsymbol{\eta})_i = \frac{\theta'(\eta_i) ( Y_i - \mu_i)}{\varphi}, \] it follows that the score equation is equivalent to the system of equations \[ \sum_{i = 1}^n \theta'_i ( Y_i - \mu_i) X_{ij} = 0, \qquad j = 1, \ldots, p. \] Note that the score equation, and thus its solution, do not depend upon the dispersion parameter. For the canonical link, the equations simplify further as \(\theta'_i = 1.\)
Chapter 6 explores existence and uniqueness of the solution to the score equation. For an arbitrary link function it can be difficult to determine if there is a solution and whether it is unique, but for certain choices the negative log-likelihood becomes convex and some guarantees can be given. The canonical link function is one such choice, and Chapter 6 completely characterizes existence and uniqueness for the canonical link function.
Example 5.4 For the normal distribution \(\mathcal{N}(\mu, \sigma^2)\) and with the canonical link function the log-likelihood function becomes \[ \begin{align*} \ell(\beta) & = \frac{1}{\varphi} \sum_{i=1}^n Y_i X_i^T \beta - \frac{(X_i^T \beta)^2}{2} \\ & = \frac{1}{2\varphi} \left( 2 \mathbf{Y}^T \mathbf{X} \beta - \beta^T \mathbf{X}^T \mathbf{X} \beta \right) \\ & = \frac{1}{2\varphi} \left( \|\mathbf{Y}\|^2 - \|\mathbf{Y} - \mathbf{X}\beta\|^2 \right). \end{align*} \] Up to the term \(\|\mathbf{Y}\|^2\)—that doesn’t depend upon the unknown \(\beta\)-vector—the log-likelihood function is proportional to the squared error loss with proportionality constant \(- 1/(2 \varphi).\) The maximum likelihood estimator is thus equal to the least squares estimator.
5.2.2 Algorithms
The score equation is nonlinear, and it does not in general have a closed form solution. Solutions must thus be found by iterative methods. Newton’s algorithm is based on a first order Taylor approximation of the score statistic. The resulting approximation of the score equation is a linear equation. Newton’s algorithm consists of iteratively computing the first order Taylor approximation and solving the resulting linear equation. A slight modification of Newton’s algorithm is often used, where the derivative of the score is replaced by its expectation, that is, by the Fisher information. To present the idea we consider a simple example of estimation in the exponential distribution with i.i.d. observations but with a nonlinear and slightly non-standard parametrization of the mean.
Example 5.5 Consider the parametrization \(\theta(\eta) = - \eta^{-k}\) for \(\eta > 0\) (and a fixed \(k > 0\)) of the canonical parameter in the exponential distribution. That is, the density is \[ e^{\theta(\eta) y - k \log \eta} \] w.r.t. Lebesgue measure on \((0, \infty).\) The mean value function is \[ \mu(\eta) = -\frac{1}{\theta(\eta)} = \eta^k. \] With \(Y_1, \ldots, Y_n\) i.i.d. observations from this distribution and with \[ S = \sum_{i=1}^n Y_i \] the score statistic amounts to \[ U(\eta) = \sum_{i=1}^n \theta'(\eta) (Y_i - \mu(\eta)) = k\eta^{- k-1} S - nk\eta^{-1}, \] where we have used that \(\theta'(\eta) = k \eta^{- k - 1}.\)
then \(Y = Z^k\) is exponentially distributed with scale parameter \(\eta^k.\) This explains the interest in this particular parametrization as it allows us to fit models with Weibull distributed outcomes.
We also find that \[ U'(\eta) = nk\eta^{-2} - k(k+1)\eta^{- k-2} S. \] We get this expression directly from the explicit formula for \(U\) above, but it is instructive to see how it also follows from the general formula, first derived in the proof of Lemma 5.1, \[ \begin{align*} U'(\eta) &= \theta''(\eta)(S - n\mu(\eta)) - n\theta'(\eta) \mu'(\eta) \\ & = - k (k + 1) \eta^{- k - 2} S + n k (k + 1) \eta^{-2} - n k^2 \eta^{- k - 1} \eta^{k -1} \\ & = n k \eta^{-2} - k (k + 1) \eta^{- k - 2} S \end{align*} \]
The score equation is the nonlinear equation \(U(\eta) = 0.\) To illustrate the general techniques for solving a nonlinear equation, we Taylor expand \(U\) around \(\eta_m\) to first order to get \[ U(\eta) \approx U(\eta_m) + U'(\eta_m)(\eta - \eta_m). \] We then swap this approximation into the score equation and solve the resulting linear equation \[ U(\eta_m) + U'(\eta_m)(\eta - \eta_m) = 0. \] Provided that \(U'(\eta_m) \neq 0,\) the solution is \[ \begin{align*} \eta_{m+1} &= \eta_m - \frac{U(\eta_m)}{U'(\eta_m)} \\ & = \eta_m - \frac{k\eta_m^{- k-1} S - nk\eta_m^{-1}}{n k \eta_m^{-2} - k (k + 1) \eta_m^{- k - 2} S } \\ & = \eta_m + \frac{\eta_m(\eta_m^{-k}S - n)}{(k+1)\eta_{m}^{-k}S - n}. \end{align*} \] This is Newton’s algorithm. With a suitable choice of starting value \(\eta_1\) we iteratively compute \(\eta_{m}\) until convergence.
Now by Lemma 5.1, the Fisher information for \(n\) independent observations equals \[ \mathcal{J}(\eta) = n \theta'(\eta) \mu'(\eta) = n k^2 \eta^{- 2}, \] which is somewhat simpler than \(U'(\eta).\) Swapping \(\mathcal{J}(\eta_m)\) for \(U'(\eta_m)\) in Newton’s algorithm gives the Fisher scoring algorithm \[ \begin{align*} \eta_{m+1} &= \eta_m + \frac{U(\eta_m)}{\mathcal{J}(\eta_m)} \\ & = \eta_m + \frac{k\eta_m^{- k-1} S - nk\eta_m^{-1}}{n k^2 \eta_m^{- 2}} \\ & = \eta_m + \frac{\eta_m(\eta_m^{-k}S - n)}{nk}. \end{align*} \] For the Fisher scoring algorithm we simply replace the denominator \((k+1)\eta^{-k}S - n\) in the Newton algorithm by its expectation \[ \mathbf{E}((k+1)\eta^{-k}S - n) = (k + 1)\eta^{-k} n \eta^{k} - n = kn. \] In this simple example, the score equation is equivalent to \[ \eta^{- k} S - n = 0, \] which is straightforward to solve analytically3. The solution is \[ \hat{\eta} = \left(\frac{1}{n} S \right)^{1/k}, \] but the example illustrates the general iterative algorithms for solving the nonlinear score equation.
3 Which shows that if \(Z_1, \ldots, Z_n\) are i.i.d. Weibull distributed with known shape parameter \(k\) and scale parameter \(\eta\) the MLE of \(\eta\) is \[\hat{\eta} = \left(\frac{1}{n} \sum_{i=1}^n Z_i^k\right)^{1/k}.\]
The general Fisher scoring algorithm for solving the score equation follows the one-dimensional example above closely. The Newton algorithm is based on the first order Taylor approximation of the score function, \[ \mathcal{U}(\beta) \approx U(\beta_m) + D_\beta \mathcal{U}(\beta_m)(\beta - \beta_m), \] and equating this approximation equal to zero gives the update \[ \beta_{m + 1} = \beta_m - D_\beta \mathcal{U}(\beta_m)^{-1}\mathcal{U}(\beta_m), \] provided that the derivative is invertible. The Fisher scoring algorithm is obtained by replacing the negative derivative by the Fisher information matrix. That is, \[ \beta_{m + 1} = \beta_m + \mathcal{J}(\beta_m)^{-1}\mathcal{U}(\beta_m). \]
Since the dispersion parameter enters as a multiplicative constant in the log-likelihood, its value does not affect the maximum likelihood estimate of \(\beta.\) We take it to be equal to \(1\) for the subsequent computations.
The derivative of the negative score statistic4 is found as in the proof of Theorem 5.1 to be \[ - D_{\beta} \mathcal{U}(\beta) = \mathbf{X}^T \mathbf{W}^{\text{obs}} \mathbf{X} \] where \[ \mathbf{W}^{\text{obs}} = - \left(\begin{array}{ccc} U_1'(\eta_{1}) & \ldots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \ldots & U_n'(\eta_{n}) \end{array}\right). \tag{5.1}\]
4 Often called the observed Fisher information.
The linearized score equation around \(\beta_m\) can be written as \[ \mathbf{X}^T U(\boldsymbol{\eta}_m) - \mathbf{X}^T \mathbf{W}^{\text{obs}}_m \mathbf{X} (\beta - \beta_m) = 0, \] and its solution is \[ \beta_{m+1} = \beta_m + (\mathbf{X}^T \mathbf{W}^{\text{obs}}_m \mathbf{X})^{-1}\mathbf{X}^T U(\boldsymbol{\eta}_m), \] provided that \(\mathbf{X}^T \mathbf{W}^{\text{obs}}_m \mathbf{X}\) has full rank \(p.\) Theorem 5.1 shows that \[ \mathcal{J}(\beta_m) = \mathbf{X}^T \mathbf{W}_m \mathbf{X}, \] where \(\mathbf{W}_m\) is the diagonal weight matrix with diagonal entries \[ w_{m,ii} = \frac{(\mu'_{m,i})^2}{\mathcal{V}(\mu_{m,i})} = \frac{(\mu'(X_i^T \beta_m))^2}{\mathcal{V}(\mu(X_i^T \beta_m))}. \]
The Fisher scoring algorithm simply amounts to replacing \(\mathbf{W}^{\text{obs}}_m\) from the Newton algorithm by \(\mathbf{W}_m.\)
We may note that the diagonal entries in \(\mathbf{W}_m\) are always strictly positive if the mean value function is strictly monotone, which implies that \(\mathbf{X}^T \mathbf{W}_m \mathbf{X}\) is positive definite and has rank \(p\) if and only if \(\mathbf{X}\) has rank \(p.\) By contrast, the diagonal weights in \(\mathbf{W}^{\text{obs}}_m\) may be negative.
We can rewrite the update formula for the Fisher scoring algorithm as follows \[ \begin{align*} \beta_{m + 1} & =\beta_{m} + (\mathbf{X}^T \mathbf{W}_m \mathbf{X})^{-1} \mathbf{X}^T U(\boldsymbol{\eta}_{m}) \\ & = (\mathbf{X}^T \mathbf{W}_m \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W}_m \Big(\underbrace{ \mathbf{X} \beta_m + \mathbf{W}^{-1}_m U(\boldsymbol{\eta}_{m})}_{\mathbf{Z}_m}\Big). \end{align*} \] The vector \(\mathbf{Z}_m\) is known as the working outcome, and its coordinates can be written out as \[ Z_{m,i} = X_i^T \beta_{m} + \frac{Y_i - \mu_{m, i}}{\mu'_{m,i}}. \tag{5.2}\] In terms of the working outcome, the vector \(\beta_{m + 1}\) is the minimizer of the weighted squared error loss \[ (\mathbf{Z}_m- \mathbf{X}\beta)^T \mathbf{W}_m (\mathbf{Z}_m - \mathbf{X} \beta) = \|\mathbf{Z}_m- \mathbf{X}\beta \|_{\mathbf{W}_m}^2 , \tag{5.3}\] see Theorem 2.1. The Fisher scoring algorithm for generalized linear models is known as iterative weighted least squares (IWLS), since it can be understood as iteratively solving a weighted least squares problem. To implement the algorithm we can also rely on general solvers of weighted least squares problems. This results in the following version of IWLS. Given \(\beta_1\) we iterate over the steps 1–3 until convergence:
- Compute the working outcome vector \(\mathbf{Z}_m\) based on \(\beta_m\) using (5.2).
- Compute the weights \[ w_{m,ii} = \frac{(\mu'_{m,i})^2}{\mathcal{V}(\mu_{m,i})}. \]
- Compute \(\beta_{m+1}\) by minimizing the weighted sum of squares (5.3).
It is noteworthy that the computations only rely on the mean value function \(\mu,\) its derivative \(\mu'\) and the variance function \(\mathcal{V}.\) Thus the IWLS algorithm depends on the mean and variance structure, as specified in the assumptions GA1 and GA2, and not on any other aspects of the outcome distribution.
Exercise 5.5 explores how the IWLS algorithm can be derived for general outcome distributions. From a general viewpoint of maximum-likelihood estimation in a regression model there is nothing special about the exponential dispersion distributions. However, for the exponential dispersion distributions the algorithm becomes—as demonstrated above—an algorithm that only depends on the mean and variance structure of the model.
5.3 Sampling distributions
It is not possible to derive exact distributional results about the estimators or the various test statistics that are used for generalized linear models—except for the special case of the linear model with normal errors. We will need to rely on approximations. By introducing an oracle least squares estimator and reinterpreting the IWLS algorithm in this context, it is possible to introduce sensible approximate moment results under the weak GA1 and GA2 assumptions. Approximate confidence intervals are subsequently introduced based on \(Z\)-scores, and tests of linear hypotheses on \(\beta\) under the GA3 assumption are discussed using likelihood ratio (deviance) tests. The operational procedures for using the tests are the same as for the linear model—the only difference being that the distributions used to compute \(p\)-values are approximations. The formal asymptotic justifications will only be treated briefly.
5.3.1 An oracle weighted least squares estimator
We introduce in this section a so-called oracle weighted least squares estimator with weights and outcomes depending on the unknown parameters. Even though we cannot compute this estimator in practice, we can derive exact moment results for it, and we can interpret the IWLS algorithm as computing an approximation. The oracle estimator thus gives us an idea about the sampling properties of the estimator computed by IWLS. Moreover, the oracle estimator is a solution to an oracle weighted least squares problem, and its derivation underscores the fact that the MLE for a generalized linear model only draws on the model assumptions GA1 and GA2.
Pretending that we know the mean values \(\mu_i = \mu(X_i^T \beta^0)\) of \(Y_i\) for \(i = 1, \ldots n,\) we can form the weighted quadratic loss \[ \mathrm{Loss}_0(\beta) \coloneqq \sum_{i=1}^n \mathcal{V}(\mu_i)^{-1} (Y_i - \mu(X_i^T \beta))^2. \] We refer to \(\mathrm{Loss}_0\) as an oracle loss function since it contains components that only an oracle would know—here the values \(\mu_i^0.\) Thus we cannot in practice compute or minimize \(\mathrm{Loss}_0\) to estimate \(\beta^0.\)
Since the parameter \(\beta\) enters nonlinearly in the mean value above through \(\mu(X_i^T \beta),\) we cannot compute the minimizer of \(\mathrm{Loss}_0\) analytically either. By a first order Taylor expansion5 around \(\beta^0\) we get \[ \mu(X_i^T \beta) \approx \mu_i + \mu_i'X_i^T (\beta - \beta^0). \]
5 This is the key step in the Gauss-Newton algorithm for nonlinear least squares estimation.
Here \(\mu_i' = \mu'(X_i^T \beta^0).\) Plugging this approximation into the oracle loss yields the quadratic oracle loss \[ \begin{align*} \mathrm{Loss}_1(\beta) & \coloneqq \sum_{i=1}^n \mathcal{V}(\mu_i)^{-1} \Big(Y_i - \mu_i - \mu_i'X_i^T (\beta - \beta^0)\Big)^2 \\ & =\sum_{i=1}^n \frac{(\mu_i')^2}{\mathcal{V}(\mu_i)} \left(X_i^T \beta^0 + \frac{Y_i - \mu_i}{\mu_i'} - X_i^T \beta\right)^2. \end{align*} \]
By introducing the oracle outcomes \[ Z_i \coloneqq X_i^T \beta^0 + \frac{Y_i - \mu_i}{\mu_i'} \tag{5.4}\] and the oracle weights \[ w_{ii} \coloneqq \frac{(\mu_i')^2}{\mathcal{V}(\mu_i)} \geq 0, \] the quadratic oracle loss can be written as \[ \mathrm{Loss}_1(\beta) = (\mathbf{Z} - \mathbf{X} \beta)^T \mathbf{W} (\mathbf{Z} - \mathbf{X} \beta) = \|\mathbf{Z} - \mathbf{X} \beta\|^2_{\mathbf{W}}, \] with \(\mathbf{W}\) the diagonal weight matrix with the \(w_{ii}\)-s in the diagonal. By Theorem 2.1, the oracle weighted least squares estimator—the minimizer of \(\mathrm{Loss}_1(\beta)\)—equals \[ \hat{\beta}^{\mathrm{oracle}} = (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^T \mathbf{W} \mathbf{Z} \] whenever \(\mathbf{X}^T \mathbf{W} \mathbf{X}\) has full rank.
Under the assumptions GA1 and GA2 we find that \[ \mathbf{E} (Z_i \mid X_i) = X_i^T \beta^0 \] and \[ \mathbf{V} (Z_i \mid X_i) = \varphi w_{ii}^{-1}. \]
If A4 additionally holds, we have the following exact distributional results about \(\hat{\beta}^{\mathrm{oracle}}\): \[ \begin{align*} \mathbf{E}(\hat{\beta}^{\mathrm{oracle}} \mid \mathbf{X}) & = \beta^0, \\ \mathbf{V}(\hat{\beta}^{\mathrm{oracle}} \mid \mathbf{X}) & = \varphi (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1}, \\ \mathbf{E}(\|\mathbf{Z} - \mathbf{X} \hat{\beta}^{\mathrm{oracle}}\|^2_{\mathbf{W}} \mid \mathbf{X}) & = (n-p) \varphi. \end{align*} \]
Since the weights as well as \(\mathbf{Z}\) depend upon the unknown \(\beta^0,\) we cannot compute \(\mathrm{Loss}(\beta),\) but the IWLS algorithm can be seen as iteratively approximating the oracle estimator—by plugging in the current estimate of \(\beta^0\) in the computation of weights and the working outcomes. If we ultimately plug in the maximum likelihood estimator, \(\hat{\beta},\) in for \(\beta^0\) in the definition of the oracle quadratic loss, and thus replacing the oracle outcomes by \[ \hat{Z}_i = X_i^T \hat{\beta} + \frac{Y_i - \hat{\mu}_i}{\hat{\mu}_i'} \] we obtain a computable loss function \[ \beta \mapsto \|\hat{\mathbf{Z}} - \mathbf{X} \beta\|^2_{\hat{\mathbf{W}}}. \tag{5.5}\]
The fact that \(\hat{\beta}\) is a fixed point for the IWLS algorithm implies that \(\hat{\beta}\) is also the minimizer of (5.5). Provided that (5.5) is a good approximation of the oracle quadratic loss, we can expect that the distributional results on the oracle estimator will be good approximations for \(\hat{\beta}\) as well.
Finally, we observe that with \[ \mathcal{X}^2 \coloneqq \|\hat{\mathbf{Z}} - \mathbf{X} \hat{\beta}\|^2_{\hat{\mathbf{W}}} = \sum_{i=1}^n \frac{(Y_i - \hat{\mu}_i)^2}{\mathcal{V}(\hat{\mu}_i)}, \] where \(\mathcal{X}^2\) is known as the Pearson \(\chi^2\)-statistic, the results above suggest the estimator \[ \hat{\varphi} = \frac{1}{n-p} \mathcal{X}^2 \tag{5.6}\] of the dispersion parameter. This is the standard estimator of the dispersion parameter that we will use throughout.
5.3.2 Tests and confidence intervals
The distributional results for the oracle weighted least squares estimator suggest the approximation \[\sqrt{\varphi (\mathbf{X}^T \mathbf{W} \mathbf{X})^{-1}_{jj}}\] of the standard error of \(\hat{\beta}_j.\) By plugging in the estimate of the dispersion parameter and the estimate of the weights, we define the \(j\)-th \(Z\)-score as \[Z_j \coloneqq \frac{\hat{\beta}_j - \beta_j}{\sqrt{\hat{\varphi} ((\mathbf{X}^T \hat{\mathbf{W}} \mathbf{X})^{-1})_{jj}}}.\] More generally, we can introduce \[Z_a \coloneqq \frac{a^T \hat{\beta} - a^T \beta}{\sqrt{\hat{\varphi} a^T (\mathbf{X}^T \hat{\mathbf{W}} \mathbf{X})^{-1}a}}\] for \(a \in \mathbb{R}^p.\)
The \(Z\)-score is used exactly as it is used for linear models. First, it is used to test a one-dimensional restriction on the parameter vector—typically a hypothesis of the form \(H_0 : \beta_j = 0\) for some index \(j.\) The square of a \(Z\)-statistic, \(Z_a^2,\) is known as a univariate Wald statistic. It is possible to construct multivariate Wald statistics to test hypotheses about multivariate parameter constraints, but we will not pursue the construction here. Second, the \(Z\)-score is used to compute confidence intervals for \(a^T \beta\) of the form \[ a^T\hat{\beta} \pm z \cdot \sqrt{\hat{\varphi} a^T(\mathbf{X}^T\hat{\mathbf{W}}\mathbf{X})^{-1}a} \tag{5.7}\] where \(z\) is chosen suitably.
A distributional approximation of \(Z_a\) is needed for computing either a \(p\)-value associated with a hypothesis test or for choosing the quantile \(z\) in the construction of the confidence interval. Asymptotic theory supports using the \(\mathcal{N}(0, 1)\)-approximation. If we aim for a confidence interval with 95% nominal coverage, we choose \(z = 1.96\) as the 97.5% quantile for the standard normal distribution. The actual coverage depends upon how well the distribution of \(Z_a\) is approximated by \(\mathcal{N}(0, 1).\)
An alternative to using theoretical approximations is to use bootstrapping to compute approximating quantiles, see Section 13.3. Occasionally, \(z\) is chosen as the 97.5% quantile in the \(t_{n-p}\)-distribution instead—still aiming for a 95% coverage. The theoretical support for this practice is weak, but the practical consequence is clear. Since the tail quantiles for the \(t\)-distribution are larger than the tail quantiles for the normal distribution, the use of the \(t\)-distribution results in wider confidence intervals and more conservative conclusions. The difference is, however, negligible unless \(n-p\) is small—less than 20, say.
The likelihood ratio tests are alternatives to \(z\)-tests and Wald tests. Their derivation is based on the stronger distributional assumption GA3. That is, the outcomes follow an exponential dispersion distribution conditionally on the predictors. In the framework of generalized linear models the likelihood ratio tests are usually formulated in terms of deviances as introduced in Section 4.3.
For a generalized linear model we have outcome observations \(Y_1, \ldots, Y_n,\) and we let \(\hat{\mu}_i\) denote the maximum likelihood estimate of the mean of \(Y_i.\) That is, \[ \hat{\mu}_i = \mu(X_i^T \hat{\beta}) \] with \(\hat{\beta}\) the maximum likelihood estimate of \(\beta \in\mathbb{R}^p.\) Likewise, if \(\hat{\beta}^0\) denotes the maximum likelihood estimate of \(\beta\) under the null hypothesis \[ H_0: \beta \in L, \] where \(L \subseteq \mathbb{R}^p\) is a \(p_0\) dimensional subspace of \(\mathbb{R}^p,\) we let \(\hat{\mu}_i^0\) denote the corresponding maximum likelihood estimate of the mean for the \(i\)-th observation under the hypothesis. Recall that according to Equation 2.12, all such hypotheses can be rephrased in terms of a \(p \times p_0\) matrix \(\mathbf{C}\) of rank \(p_0\) such that \[ H_0 : \mathbf{E}(Y_i \mid X_i) = \mu(X_i^T \mathbf{C} \beta^0). \]
Definition 5.2 The total deviances for the model and the null hypothesis are \[ D = \sum_{i=1}^n d(Y_i, \hat{\mu}_i) \quad \textrm{ and} \quad D_0 = \sum_{i=1}^n d(Y_i, \hat{\mu}_i^0), \] respectively. The deviance test statistic is \[ D_0 - D \] and the \(F\)-test statistic is \[ F = \frac{(D_0 - D)/(p - p_0)}{D/(n-p)}. \]
The deviance test statistic is simply the log-likelihood ratio test statistic for the null hypothesis. The \(F\)-test statistic is inspired by the \(F\)-test for linear models. In both cases, large values of the test statistic are critical, and thus evidence against the null hypothesis. In contrast to the linear model, it is not possible to derive the exact distributions of the deviance test statistic or the \(F\)-test statistic under the null hypothesis. To compute \(p\)-values we need to rely on approximations or simulations. The general theoretical approximation of the deviance test statistic is
\[
D_0 - D \overset{\mathrm{approx}}{\sim} \varphi \chi^2_{p - p_0},
\tag{5.8}\] which is exact for the linear model with normal errors (under the assumptions A3 + A5). A purpose of the \(F\)-test is, as for linear models, to (approximately) remove the dependence upon the unknown dispersion parameter. The approximation of the \(F\)-test statistic is \[
F \overset{\mathrm{approx}}{\sim} F_{p - p_0, n - p},
\] which is exact for the linear model with normal errors. The approximation is generally good when \(\varphi\) is small.
A formal justification of the approximations is quite involved, and the typical strategy is to consider asymptotic scenarios with \(n \to \infty\) combined with suitable conditions on the predictors \(X_i.\) We will explain the basis for these approximations later, but we will not consider formal asymptotic arguments.
5.4 Model assessment and diagnostic
The definition of the adjusted \(R^2\) can easily be extended to generalized linear models by replacing \(\mathrm{RSS}/(n-p)\) with \(\mathcal{X}^2/(n-p),\) where \[ \mathcal{X}^2 = \sum_{i=1}^n \frac{(Y_i - \hat{\mu}_i)^2}{\mathcal{V}(\hat{\mu}_i)} \] is the Pearson \(\chi^2\)-statistic, or by \(D/(n-p),\) where \[D = \sum_{i=1}^n d(Y_i, \hat{\mu}_i)\] is the total deviance.
As for the linear model, model diagnostics can be based on residuals, but for generalized linear models there are several possible choices. With observations \(Y_1, \ldots, Y_n\) and fitted mean values \(\hat{\mu}_i, \ldots, \hat{\mu}_n\) the \(i\)-th raw residual is \[ Y_i - \hat{\mu}_i. \] In the R terminology the raw residual is called the response residual.
Whenever the variance function is not constant, the raw residual is not particularly useful. To take a non-constant variance function into account, a natural choice of residual is the Pearson residual defined as \[ \frac{Y_i - \hat{\mu}_i}{\sqrt{\mathcal{V}(\hat{\mu}_i)}}. \] Pearson residuals are used in the same way as the raw or standardized residuals are used for the linear model. We plot the Pearson residuals against the fitted values or a predictor variable to visually check the model assumptions—the residuals should show no distributional dependence upon what we plot it against. Systematic trends show that GA1 is not fulfilled, and variance inhomogeneity shows that the variance function in GA2 is not appropriate.
The Pearson residual is based on the mean and variance assumptions GA1 and GA2 only. Based on the distributional assumption GA3 we can also introduce the deviance residual for the \(i\)-th observation as \[ \mathrm{sign}(Y_i - \hat{\mu}_i) \sqrt{d(Y_i, \hat{\mu}_i)}. \] The deviance residuals can be used like the Pearson residuals, and Theorem 4.2 gives an approximate relation between them.
Finally, the working residual should be mentioned. It is defined for the \(i\)-th observation as \[ \frac{Y_i - \hat{\mu}_i}{\mu'(\hat{\eta}_i)}, \] and is the residual from the weighted least squares problem (5.5).
Example 5.6 For the binomial case we find that the raw residual is \[ Y_i - m_i \hat{p}_i \] where \(\hat{p}_i\) is the estimate of the success probability for the \(i\)-th observation.
The Pearson residual is \[ \frac{Y_i - m_i \hat{p}_i}{\sqrt{m_i \hat{p}_i(1-\hat{p}_i)}}, \] and the deviance residual is \[ \mathrm{sign}(Y_i - m_i \hat{p}_i) \sqrt{2 Y_i \log\frac{Y_i}{m_i\hat{p}_i} + 2 (m_i - Y_i) \log \frac{m_i - Y_i}{m_i(1-\hat{p}_i)}}. \]
The distribution of the residual is difficult to characterize in general. We cannot expect that the Pearson residuals follow a normal distribution, say. For generalized linear models there is thus no simple way to check if the strong distributional assumption— that the outcomes follow an exponential dispersion distribution given the predictors—are fulfilled, and the residuals are primarily used to investigate the mean-variance structure.
If the outcomes have a continuous distribution with distribution function \(F_{\mu_i, \varphi}\) we can, however, construct a pp-plot by comparing \[ F_{\hat{\mu}_i, \hat{\varphi}}(Y_i) \] to the quantiles of the uniform distribution. This can be useful for gamma distributed outcomes but not for Poisson distributed reponses, say, unless the counts are quite large.
Exercises
Exercise 5.1 Define the probability distribution on \(\{0,1, \ldots, m\}\) by the probability mass function \[ p_{\eta}(k) = \frac{1}{c(\eta)}\eta^k \] with \(c(\eta) = \sum_{k=0}^m \eta^k\) for \(\eta > 0.\)
Show that the family of probability distributions given by \(p_{\eta}\) for \(\eta > 0\) is an exponential family. Find the structure measure and the unit cumulant function.
Compute the mean value function \(\mu(\eta)\) and show that \[ \mathcal{V}(\mu(\eta)) = \frac{\eta}{(1-\eta)^2} - \frac{(m+1)^2 \eta^{m+1}}{(1-\eta^{m+1})^2} \] for \(\eta \neq 1.\) Plot the link function for \(m = 5,\) \(10,\) \(20.\)
Hint: You can plot the link function without a closed form expression.
Consider the following data set with one explanatory variable and \(n=10\) observations.
| \(i\) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| \(X_i\) | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 | 1.1 | 1.2 |
| \(Y_i\) | 1 | 5 | 2 | 8 | 7 | 5 | 9 | 7 | 10 | 10 |
Consider the generalized linear model with outcome distribution as given above, with the linear predictor \[ \eta_i = \beta_0 + \beta_1 x_i. \]
Implement the Fisher scoring algorithm to estimate \(\beta_0\) and \(\beta_1.\) Use as initial values \(\beta_0 = 0\) and \(\beta_1 = 1\) and report the estimated parameters for each iteration of the algorithm with a total of 10 iterations. Compute also an estimate of the standard error of the MLE \(\hat{\beta}_1\) (assuming that the algorithm has converged to the MLE after 10 iterations).
Compute and plot the Pearson residuals against the fitted values.
Introduce a dispersion parameter \(\varphi > 0\) and estimate it using the Pearson residuals. Recompute an estimate of the standard error of \(\hat{\beta}_1\) taking the estimated dispersion parameter into account. Comment on the result.
Exercise 5.2 Define a probability distribution on the non-negative integers, \(0, 1, 2, \ldots,\) as having probability mass function \[ p_{\eta}(k) = G(k, \rho) \left(\frac{1}{\rho+\eta}\right)^{\rho} \left(\frac{\eta}{\rho + \eta}\right)^k \tag{5.9}\] for \(\eta, \rho > 0.\) The explicit expression for \(G(k, \rho)\) is not needed—only that \(\sum_{k=0}^{\infty} p_{\eta}(k) = 1.\) In the following, \(\rho\) is considered a fixed nuisance parameter.
Show that the family of probability distributions given by \(p_{\eta}\) for \(\eta > 0\) is an exponential family for any fixed \(\rho > 0.\) Find the structure measure and the unit cumulant function.
Show that the mean value function is \(\mu(\eta) = \eta\) and that \[ \mathcal{V}(\mu(\eta)) = \eta + \eta^2/\rho. \]
Consider the following data set with one explanatory variable and \(n=10\) observations.
| \(i\) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| \(X_i\) | 5 | 10 | 15 | 20 | 25 | 30 | 35 | 40 | 45 | 50 |
| \(Y_i\) | 6 | 11 | 16 | 38 | 82 | 22 | 108 | 61 | 55 | 64 |
Consider the generalized linear model with outcome distribution as given above and with the linear predictor \[ \eta_i = \beta_0 + \beta_1 X_i. \]
Take \(\rho = 1\) and compute the MLE of \(\beta_0\) and \(\beta_1.\) You can use the Fisher-scoring algorithm. How will it affect the result if \(\rho\) is changed?
Still taking \(\rho = 1\) compute and plot the Pearson residuals against the fitted values and compute the \(\chi^2\)-statistic.
Write an R function that computes the \(\chi^2\)-statistic as a function of \(\rho\) using the general variance formula. Plot this \(\chi^2\)-statistic as a function of \(\rho\) and use it to argue that the data suggests a value of \(\rho = 4.\)
Compute a nominal 95% confidence interval of \(\beta_1\) based on the MLE when assuming \(\rho = 1.\) Then compute a nominal 95% confidence interval assuming that \(\rho = 4\) is the true value of \(\rho\) while the MLE is still based on assuming \(\rho = 1.\) Comment on the result.
Exercise 5.3 Recall from Example 5.3 that for the log-normal random variable, \(Y,\) whose distribution is given as \[ \log(Y) = \eta + \varepsilon \tag{5.10}\] for \(\eta \in \mathbb{R}\) the linear predictor and \(\varepsilon \sim \mathcal{N}(0,\sigma^2),\) GA1 holds with mean value function \[ \mu(\eta) =e^{\eta + \frac{\sigma^2}{2}} \] and GA2 holds with variance function \(\mathcal{V}(\mu) = \mu^2.\) The dispersion parameter is given as \[ \varphi = e^{\sigma^2} - 1 \] in terms of \(\sigma^2.\)
Introduce the parameter \(\tilde{\eta} = \eta + \frac{\sigma^2}{2}.\) Whenever the linear predictor includes an intercept, the term \(\frac{\sigma^2}{2}\) is simply absorbed into the intercept parameter, and we can regard \(\tilde{\eta}\) as the linear predictor with this reparametrization of the intercept.
- Identify the exponential dispersion model with mean value function \(\mu(\tilde{\eta}) = e^{\tilde{\eta}}\) and variance function \(\mathcal{V}(\mu) = \mu^2.\) Show that the unit deviance for this model is \[ d(y, \mu) = 2 \left(\log\left(\frac{\mu}{y} \right) + \frac{y}{\mu} - 1\right).\]
In general, the coefficient of variation of \(Y\) is defined as \(\sqrt{\mathbf{V}(Y)}/\mu.\) For the log-normal model as well as the exponential dispersion model above, this coefficient is constantly equal to \(\sqrt{\varphi}\) and independent of \(\mu.\) Both models are thus models using a log-link function and with a constant coefficient of variation.
Deviations in the model given by (5.10) are naturally measured in terms of the squared error \((\log(y) - \eta)^2,\) which is the unit deviance for the normal distribution, and the model is typically fitted on this log-scale. The exponential dispersion model is typically fitted using the IWLS algorithm.
- First fix \(y = 1\) and plot \(\tilde{\eta} \mapsto d(1, e^{\tilde{\eta}})\) for \(\tilde{\eta} \in (-2, 2).\) Compare with \(\tilde{\eta} \mapsto \tilde{\eta}^2.\) Then show that in general \[ d(y, e^{\tilde{\eta}}) = (\log(Y) - \tilde{\eta})^2 + o((\log(y) -\tilde{\eta})^2). \] What does this result imply about the relation between the least squares fit based on \(\log(y)\) and the IWLS fit based on \(y\)?
Exercise 5.4 This exercise is on the Danish fire insurance data introduced in Chapter 2. You should remove claims larger that 10 millions to avoid numerical problems. Recall that the outcome variable \(Y\) is the claim size and there are two potential predictors in the data set.
- Refit the linear, additive model \[ \mathbf{E} (\log(Y_i)) = \beta_0 + \beta_{X_{i,\texttt{grp}}} + \beta_{\texttt{sum}} \log X_{i, \texttt{sum}} \tag{5.11}\] as considered in Chapter 2, but using the reduced dataset. Construct model diagnostic plots for this model.
As shown in Example 5.3, and further elaborated on in Exercise 5.3, the log-normal regression model is also a generalized linear model with a log-link and a quadratic variance function. That is, \[ \log \mathbf{E} (Y_i) = \beta_0 + \beta_{X_{i,\texttt{grp}}} + \beta_{\texttt{sum}} \log X_{i, \texttt{sum}}. \tag{5.12}\]
Show that with a log-link and a quadratic variance function, \(\mathcal{V}(\mu) = \mu^2,\) the weights in the IWLS algorithm are always 1.
Fit the model given by (5.12). Investigate if the model fits the data. Compare the fitted values with those from the log-normal model.
An alternative to the log-normal model is the log-gamma model, where (5.11) still specifies the mean value, but where the log-outcome is assumed to have a gamma distribution. Such a model can be fitted as a generalized linear model with a gamma outcome, but with the identity link.
Fit the log-gamma model and discuss the model fit.
Construct pp-plots for both the log-normal and log-gamma model and discuss which outcome distribution appears most appropriate.
Exercise 5.5 Assume in this exercise that \(Y_i \mid X_i\) has distribution with density \(f_i(y_i, \eta_i)\) for \(i = 1, \ldots, n,\) where \(\eta_i = X_i^T \beta\) is the linear predictor. The generalized linear model with an exponential dispersion distribution is a special case with \[ \log f_i(y_i, \eta_i) = \frac{\theta(\eta_i) y_i - \zeta(\eta_i)}{\varphi}. \] The purpose of the exercise is to investigate the computations of the maximum-likelihood estimator in this more general setup, and to understand which results are specific to exponential dispersion models. It is assumed that \(f_i(y_i, \eta_i) > 0,\) and we define \[ U_i(\eta_i) = (\log f_i(Y_i, \eta_i))' \] and \[ w_i = -\mathbf{E} \left( (\log f_i(Y_i, \eta_i))'' \right) \] where the differentiation is w.r.t. \(\eta_i.\) Introduce also the score vector \(U(\boldsymbol{\eta}) = (U_1(\eta_1), \ldots, U_n(\eta_n))^T\) and the diagonal matrix \(\mathbf{W}\) with \(\mathbf{W}_{ii} = w_i.\)
Prove that the score statistic equals \(\mathbf{X}^TU(\boldsymbol{\eta})\) and that the Fisher information equals \(\mathbf{X}^T \mathbf{W} \mathbf{X}.\)
Prove that the \(m\)-th Fisher scoring step for solving the score equation amounts to solving the linear equation \[ \mathbf{X}^TU(\boldsymbol{\eta}_m) - (\mathbf{X}^T \mathbf{W}_m \mathbf{X}) (\beta - \beta_m) = 0. \] Here \(\mathbf{W}_m\) is defined in terms of \(\boldsymbol{\eta}_m = \mathbf{X} \beta_m.\)
Show that if \(w_{m, i} \geq 0\) then \(\beta_{m + 1}\) is the solution of a weighted least squares problem with diagonal weight matrix \(\mathbf{W}_m\) and working outcomes \(X_i^T \beta_m + U_i(\eta_m) / w_{m,i}.\)
Since \(w_i\) is the variance of the score statistic (and thus nonnegative) under some regularity conditions, estimation can be carried out by iteratively weighted least squares also in the more general case considered in this exercise. However, this requires formulas for the score statistic and the weights, which cannot generally be obtained from the first two moments.