Solutions to Dobson & Barnett’s An Introduction to Generalized Linear Models

Write the transformation in matrix form. With \(\mathbf{y} = [Y_1, Y_2]^T\) and

\[ \mathbf{w} = \begin{pmatrix} W_1 \\ W_2 \end{pmatrix} = \mathbf{A}\mathbf{y}, \qquad \mathbf{A} = \begin{pmatrix} 1 & 2 \\ 4 & -1 \end{pmatrix}, \]

the two components of \(\mathbf{w}\) are linear combinations of independent Normal variables, so by Section 1.4.1 result 4 each of \(W_1, W_2\) is Normal; more than that, every linear combination \(a_1 W_1 + a_2 W_2\) is again a linear combination of \(Y_1\) and \(Y_2\) and hence Normal, which is exactly the condition for \(\mathbf{w}\) to be bivariate Normal. So only the mean vector and the variance-covariance matrix remain to be found.

Since the \(Y_i\) are independent,

\[ \boldsymbol{\mu}_y = \begin{pmatrix} 1 \\ 2 \end{pmatrix}, \qquad \mathbf{V}_y = \begin{pmatrix} 3 & 0 \\ 0 & 5 \end{pmatrix}. \]

By Equations (1.2) and (1.3), \(E(\mathbf{w}) = \mathbf{A}\boldsymbol{\mu}_y\) and \(\mathrm{var}(\mathbf{w}) = \mathbf{A}\mathbf{V}_y\mathbf{A}^T\). Componentwise,

\[ E(W_1) = 1 + 2(2) = 5, \qquad E(W_2) = 4(1) - 2 = 2, \] \[ \mathrm{var}(W_1) = 1^2(3) + 2^2(5) = 23, \qquad \mathrm{var}(W_2) = 4^2(3) + (-1)^2(5) = 53, \] \[ \mathrm{cov}(W_1, W_2) = \mathrm{cov}(Y_1 + 2Y_2,\; 4Y_1 - Y_2) = 4\,\mathrm{var}(Y_1) - 2\,\mathrm{var}(Y_2) = 4(3) - 2(5) = 2, \]

where the cross terms vanish because \(\mathrm{cov}(Y_1, Y_2) = 0\).

Hence

\[ \begin{pmatrix} W_1 \\ W_2 \end{pmatrix} \sim \mathrm{MVN}\left( \begin{pmatrix} 5 \\ 2 \end{pmatrix}, \begin{pmatrix} 23 & 2 \\ 2 & 53 \end{pmatrix} \right). \]

The matrix arithmetic, checked in R:

A  <- matrix(c(1, 2, 4, -1), nrow = 2, byrow = TRUE)
mu <- c(1, 2)
V  <- diag(c(3, 5))
A %*% mu
A %*% V %*% t(A)

     [,1]
[1,]    5
[2,]    2
     [,1] [,2]
[1,]   23    2
[2,]    2   53

Note that \(W_1\) and \(W_2\) are not independent even though \(Y_1\) and \(Y_2\) are: the correlation is \(2/\sqrt{23 \times 53} = 0.057\). It is small only by coincidence – the positive contribution \(4\,\mathrm{var}(Y_1) = 12\) nearly cancels the negative contribution \(-2\,\mathrm{var}(Y_2) = -10\). The covariance \(4\sigma_1^2 - 2\sigma_2^2\) vanishes exactly when \(\sigma_1^2 : \sigma_2^2 = 1 : 2\); here the ratio is \(3 : 5\), close to but not equal to \(1 : 2\). Had it been exactly \(1 : 2\), the two linear combinations would have been uncorrelated and, being jointly Normal, independent.

Part (a): \(Y_1\) already has the standard Normal distribution, so by the definition of the central chi-squared distribution in Section 1.4.2 result 1 with \(n = 1\),

\[ Y_1^2 \sim \chi^2(1). \]

Part (b): The second element standardises \(Y_2\): since \(Y_2 \sim N(3, 4)\) has \(\mu_2 = 3\) and \(\sigma_2 = 2\), the variable \(Z_2 = (Y_2 - 3)/2 \sim N(0,1)\). Therefore

\[ \mathbf{y}^T\mathbf{y} = Y_1^2 + \frac{(Y_2 - 3)^2}{4}. \]

This is the sum of squares of two independent standard Normal variables, which is precisely the form in Equation (1.4) with \(n = 2\):

\[ \mathbf{y}^T\mathbf{y} \sim \chi^2(2). \]

Part (c): Here \(\mathbf{y}\) is not centred. Its mean vector and variance-covariance matrix are

\[ \boldsymbol{\mu} = \begin{pmatrix} 0 \\ 3 \end{pmatrix}, \qquad \mathbf{V} = \begin{pmatrix} 1 & 0 \\ 0 & 4 \end{pmatrix}, \qquad \mathbf{V}^{-1} = \begin{pmatrix} 1 & 0 \\ 0 & 1/4 \end{pmatrix}, \]

the off-diagonal entries being zero because \(Y_1\) and \(Y_2\) are independent. Hence

\[ \mathbf{y}^T\mathbf{V}^{-1}\mathbf{y} = Y_1^2 + \frac{Y_2^2}{4}. \]

By Section 1.4.2 result 6, \(\mathbf{y}^T\mathbf{V}^{-1}\mathbf{y}\) has the non-central chi-squared distribution \(\chi^2(n, \lambda)\) with \(n = 2\) degrees of freedom and non-centrality parameter

\[ \lambda = \boldsymbol{\mu}^T\mathbf{V}^{-1}\boldsymbol{\mu} = 0^2 + \frac{3^2}{4} = \frac{9}{4} = 2.25, \]

so \(\mathbf{y}^T\mathbf{V}^{-1}\mathbf{y} \sim \chi^2(2, 2.25)\).

The contrast between (b) and (c) is the whole point of the exercise. Subtracting the mean before standardising gives a central chi-squared; failing to subtract it leaves the squared standardised mean behind as the non-centrality parameter. Note \(\lambda = 2.25 = (\mu_2/\sigma_2)^2\), and indeed \(\mathbf{y}^T\mathbf{y}\) from (b) equals \(\mathbf{y}^T\mathbf{V}^{-1}\mathbf{y}\) from (c) with \(Y_2\) replaced by \(Y_2 - 3\).

A simulation confirms the moments of (c), using \(E = n + \lambda\) and \(\mathrm{var} = 2n + 4\lambda\) from Section 1.4.2 result 4:

set.seed(1)
n  <- 2e6
Y1 <- rnorm(n, 0, 1)
Y2 <- rnorm(n, 3, 2)
Q  <- Y1^2 + Y2^2 / 4
c(mean = mean(Q), var = var(Q))
c(theory.mean = 2 + 2.25, theory.var = 2 * 2 + 4 * 2.25)
ks.test(sample(Q, 5000), "pchisq", df = 2, ncp = 2.25)$p.value

     mean       var
 4.251925 13.007515
theory.mean  theory.var
       4.25       13.00
[1] 0.4517404

The simulated mean and variance match the theoretical \(4.25\) and \(13\), and a Kolmogorov-Smirnov test against \(\chi^2(2, 2.25)\) gives \(p = 0.45\), so there is no evidence against the claimed distribution.

Both parts need \(\mathbf{V}^{-1}\). For a \(2 \times 2\) matrix, \(\det \mathbf{V} = 4(9) - 1(1) = 35\) and

\[ \mathbf{V}^{-1} = \frac{1}{35}\begin{pmatrix} 9 & -1 \\ -1 & 4 \end{pmatrix}. \]

Note \(\det \mathbf{V} > 0\) and \(v_{11} = 4 > 0\), so \(\mathbf{V}\) is positive definite by the determinant criterion of Section 1.5 result 3, and the inverse exists as required.

Part (a): Writing \(u_1 = y_1 - 2\) and \(u_2 = y_2 - 3\), the quadratic form expands as

\[ (\mathbf{y} - \boldsymbol{\mu})^T\mathbf{V}^{-1}(\mathbf{y} - \boldsymbol{\mu}) = \frac{1}{35}\left[ 9u_1^2 - 2u_1u_2 + 4u_2^2 \right], \]

that is,

\[ \frac{1}{35}\left[ 9(y_1 - 2)^2 - 2(y_1 - 2)(y_2 - 3) + 4(y_2 - 3)^2 \right]. \]

The cross term picks up a factor of \(2\) because the two off-diagonal entries \(-1/35\) are equal. Since \(\mathbf{V}\) is non-singular, Equation (1.5) applies directly with \(n = 2\):

\[ (\mathbf{y} - \boldsymbol{\mu})^T\mathbf{V}^{-1}(\mathbf{y} - \boldsymbol{\mu}) \sim \chi^2(2). \]

Part (b): The same algebra without centring gives

\[ \mathbf{y}^T\mathbf{V}^{-1}\mathbf{y} = \frac{1}{35}\left[ 9y_1^2 - 2y_1y_2 + 4y_2^2 \right]. \]

By Section 1.4.2 result 6 this has the non-central chi-squared distribution with \(2\) degrees of freedom and non-centrality parameter

\[ \lambda = \boldsymbol{\mu}^T\mathbf{V}^{-1}\boldsymbol{\mu} = \frac{1}{35}\left[ 9(2)^2 - 2(2)(3) + 4(3)^2 \right] = \frac{36 - 12 + 36}{35} = \frac{60}{35} = \frac{12}{7} \approx 1.714, \]

so \(\mathbf{y}^T\mathbf{V}^{-1}\mathbf{y} \sim \chi^2(2, 12/7)\).

Confirming the inverse and \(\lambda\) in R:

mu <- c(2, 3)
V  <- matrix(c(4, 1, 1, 9), nrow = 2)
det(V)
Vinv <- solve(V)
35 * Vinv
lambda <- drop(t(mu) %*% Vinv %*% mu)
c(lambda = lambda, as.fraction = 60 / 35)

[1] 35
     [,1] [,2]
[1,]    9   -1
[2,]   -1    4
     lambda as.fraction
   1.714286    1.714286

Interpretation: the degrees of freedom are the same in both parts because \(\mathbf{V}^{-1}\) has full rank \(2\) (Section 1.5 result 4). What differs is the location. In (a) the quadratic form is a squared Mahalanobis distance from the mean, and it has expected value \(2\); in (b) it is measured from the origin, and the extra distance of \(\boldsymbol{\mu}\) from \(\mathbf{0}\) inflates the expectation to \(n + \lambda = 2 + 12/7 \approx 3.71\). Quantities of exactly form (a) reappear in later chapters as goodness-of-fit statistics; form (b) is what arises under an alternative hypothesis, and \(\lambda\) is what drives the power of the test.

Part (a): \(\overline{Y} = \sum a_i Y_i\) with every \(a_i = 1/n\), so by Section 1.4.1 result 4,

\[ \overline{Y} \sim N\!\left( \sum \frac{1}{n}\mu, \; \sum \frac{1}{n^2}\sigma^2 \right) = N\!\left( \mu, \frac{\sigma^2}{n} \right). \]

Part (b): Insert and remove \(\mu\) inside the square, writing \(Y_i - \overline{Y} = (Y_i - \mu) - (\overline{Y} - \mu)\). Then

\[ \sum_{i=1}^{n}(Y_i - \overline{Y})^2 = \sum (Y_i - \mu)^2 - 2(\overline{Y} - \mu)\sum (Y_i - \mu) + n(\overline{Y} - \mu)^2. \]

The middle sum is \(\sum (Y_i - \mu) = n\overline{Y} - n\mu = n(\overline{Y} - \mu)\), so the cross term equals \(-2n(\overline{Y} - \mu)^2\) and

\[ \sum_{i=1}^{n}(Y_i - \overline{Y})^2 = \sum (Y_i - \mu)^2 - n(\overline{Y} - \mu)^2. \]

Dividing by \(n-1\) gives the required identity. This is an algebraic identity: it holds for any numbers, no distributional assumption is used.

Part (c): Divide the identity by \(\sigma^2\) to get the decomposition stated in the question,

\[ \underbrace{\sum \frac{(Y_i - \mu)^2}{\sigma^2}}_{Q} = \underbrace{\frac{(n-1)S^2}{\sigma^2}}_{Q_1} + \underbrace{\frac{n(\overline{Y} - \mu)^2}{\sigma^2}}_{Q_2}. \]

Each \(Y_i - \mu \sim N(0, \sigma^2)\) independently, so \(Q \sim \chi^2(n)\) by Equation (1.4), and \(Q_1, Q_2\) are quadratic forms in the \(Y_i - \mu\). Their degrees of freedom (ranks) are \(m_1 = n-1\) and \(m_2 = 1\): \(Q_2\) is a single squared linear combination, hence rank \(1\), and \(Q_1\) is the sum of squared deviations about their own mean, whose matrix \(\mathbf{I} - \frac{1}{n}\mathbf{J}\) is idempotent with trace \(n-1\). Since \(m_1 + m_2 = n\), Cochran’s theorem (Section 1.5 result 5) applies and gives that \(Q_1\) and \(Q_2\) are independent, with \(Q_1 \sim \chi^2(n-1)\) and \(Q_2 \sim \chi^2(1)\).

Strictly, Cochran’s theorem delivers independence of \(S^2\) and \((\overline{Y} - \mu)^2\), not of \(S^2\) and \(\overline{Y}\) itself. The gap is easily closed and worth stating explicitly rather than glossing over: \(\mathrm{cov}(\overline{Y},\, Y_i - \overline{Y}) = \sigma^2/n - \sigma^2/n = 0\) for every \(i\), and the vector \((\overline{Y}, Y_1 - \overline{Y}, \ldots, Y_n - \overline{Y})\) is multivariate Normal because all its entries are linear in the \(Y_i\). Zero covariance plus joint Normality gives independence of \(\overline{Y}\) from the whole residual vector, hence from \(S^2\), which is a function of that vector alone. With that supplement the deduction is rigorous.

Part (d): From (c),

\[ \frac{(n-1)S^2}{\sigma^2} \sim \chi^2(n-1). \]

The single degree of freedom lost is exactly the one spent estimating \(\mu\) by \(\overline{Y}\), which is what the decomposition in (c) makes visible. It also gives \(E(S^2) = \sigma^2\) immediately, since \(E[\chi^2(n-1)] = n-1\).

Part (e): Standardise the numerator using (a): \(Z = (\overline{Y} - \mu)/(\sigma/\sqrt{n}) \sim N(0,1)\). Let \(X^2 = (n-1)S^2/\sigma^2 \sim \chi^2(n-1)\), independent of \(Z\) by (c). Then

\[ \frac{\overline{Y} - \mu}{S/\sqrt{n}} = \frac{(\overline{Y} - \mu)/(\sigma/\sqrt{n})}{S/\sigma} = \frac{Z}{\left[ X^2/(n-1) \right]^{1/2}}, \]

since \(S/\sigma = [X^2/(n-1)]^{1/2}\). That is exactly the definition of the t-distribution in Equation (1.6), so

\[ \frac{\overline{Y} - \mu}{S/\sqrt{n}} \sim t(n-1). \]

The unknown \(\sigma\) has cancelled, which is why this statistic is usable in practice: it is the one-sample t-statistic.

A simulation with \(n = 8\), \(\mu = 10\), \(\sigma = 3\) checks (c), (d) and (e) numerically:

set.seed(42)
n <- 8; mu <- 10; sigma <- 3; B <- 2e5
sim  <- matrix(rnorm(n * B, mu, sigma), nrow = B)
ybar <- rowMeans(sim)
s2   <- apply(sim, 1, var)
c(cor.ybar.s2 = cor(ybar, s2))
Q1 <- (n - 1) * s2 / sigma^2
c(mean.Q1 = mean(Q1), df = n - 1, var.Q1 = var(Q1), two.df = 2 * (n - 1))
Tstat <- (ybar - mu) / sqrt(s2 / n)
c(mean.T = mean(Tstat), var.T = var(Tstat), t.var = (n - 1) / (n - 3))
ks.test(sample(Tstat, 5000), "pt", df = n - 1)$p.value

  cor.ybar.s2
-0.0006370237
  mean.Q1        df    var.Q1    two.df
 7.007965  7.000000 14.053002 14.000000
      mean.T        var.T        t.var
0.0002216423 1.4003904492 1.4000000000
[1] 0.4874567

The sample correlation between \(\overline{Y}\) and \(S^2\) is essentially zero, consistent with (c). \(Q_1\) has simulated mean \(7.01\) and variance \(14.05\) against the \(\chi^2(7)\) values \(7\) and \(14\), confirming (d). The t-statistic has simulated variance \(1.400\) against the \(t(7)\) value \((n-1)/(n-3) = 1.4\), and a Kolmogorov-Smirnov test against \(t(7)\) gives \(p = 0.49\), confirming (e). Note the variance exceeds \(1\): the heavier tails of \(t(7)\) relative to \(N(0,1)\) are the price of estimating \(\sigma\).

Part (a): With \(f(y_i; \theta) = \theta^{y_i}e^{-\theta}/y_i!\) for \(y_i = 0, 1, 2, \ldots\),

\[ E(Y_i) = \sum_{y=0}^{\infty} y\,\frac{\theta^{y}e^{-\theta}}{y!} = \sum_{y=1}^{\infty} \frac{\theta^{y}e^{-\theta}}{(y-1)!}, \]

the \(y = 0\) term contributing nothing. Substituting \(z = y - 1\) and taking one factor of \(\theta\) outside,

\[ E(Y_i) = \theta e^{-\theta}\sum_{z=0}^{\infty}\frac{\theta^{z}}{z!} = \theta e^{-\theta}e^{\theta} = \theta, \]

using the exponential series \(\sum_{z \ge 0} \theta^{z}/z! = e^{\theta}\).

Part (b): Reparameterising the log-likelihood of Section 1.6.2 by \(\theta = e^{\beta}\),

\[ l(\beta; \mathbf{y}) = \left( \sum y_i \right)\beta - n e^{\beta} - \sum \log y_i!, \]

since \(\log \theta = \beta\) and \(n\theta = ne^{\beta}\). Differentiating,

\[ \frac{dl}{d\beta} = \sum y_i - n e^{\beta}. \]

Setting this to zero gives \(e^{\beta} = \overline{y}\), so

\[ \widehat{\beta} = \log \overline{y}, \]

provided \(\overline{y} > 0\). The second derivative is \(d^2l/d\beta^2 = -ne^{\beta} < 0\) everywhere, so \(l\) is strictly concave in \(\beta\) and \(\widehat{\beta}\) is the unique maximum, as required by the check in Section 1.6.1.

This is the invariance property of maximum likelihood estimators in action: \(\beta = g(\theta) = \log\theta\) is a one-to-one function of \(\theta\), so \(\widehat{\beta} = g(\widehat{\theta}\,) = \log \overline{y}\) without redoing the calculus. Note also that \(\beta\) is the natural parameter of the Poisson distribution, and \(\theta = e^{\beta}\) is the log link that will define the Poisson generalized linear model in Chapter 9.

Part (c): Put \(u = e^{\beta}\), so \(S = \sum (Y_i - u)^2\) and, by the chain rule,

\[ \frac{dS}{d\beta} = \frac{dS}{du}\cdot\frac{du}{d\beta} = \left[ -2\sum (Y_i - u) \right] e^{\beta} = -2e^{\beta}\left( \sum Y_i - n e^{\beta} \right). \]

Because \(e^{\beta} > 0\) for every real \(\beta\), the derivative vanishes only where \(\sum Y_i = ne^{\beta}\), that is at

\[ \widetilde{\beta} = \log \overline{y}. \]

To confirm a minimum, note \(d^2S/d\beta^2 = e^{2\beta}\,d^2S/du^2 + (dS/du)e^{\beta}\), and at the stationary point \(dS/du = 0\), leaving \(e^{2\beta}(2n) > 0\).

So the least squares and maximum likelihood estimators coincide here, illustrating comment 2 of Section 1.6.4. The agreement is not automatic: least squares needs only \(E(Y_i) = e^{\beta}\) whereas maximum likelihood needs the full Poisson specification. They agree because both estimating equations reduce to the same moment condition \(\sum(Y_i - e^{\beta}) = 0\). The agreement also relies on the \(Y_i\) being identically distributed here: all \(n\) variables share the single parameter \(\theta\), so \(\mathrm{var}(Y_i) = \theta\) is the same for every \(i\) and unweighted least squares is already the correctly weighted one (Section 1.6.3, weights \(w_i = 1/\sigma_i^2\)). Once explanatory variables enter and the means differ across observations, the Poisson variances differ too, since variance equals mean, the two methods separate, and it is weighted least squares that tracks maximum likelihood.

Applying this to the tropical cyclone data of Table 1.2, where Section 1.6.5 reports \(\widehat{\theta} = 72/13 = 5.538\):

library(dobson)
data(cyclones)
y <- cyclones$number
c(n = length(y), total = sum(y), ybar = mean(y))
betahat <- log(mean(y))
c(betahat = betahat, exp.betahat = exp(betahat))

negll <- function(b) -sum(y * b - exp(b) - lfactorial(y))
opt <- optimize(negll, interval = c(-5, 5), tol = 1e-10)
c(numerical.betahat = opt$minimum, closed.form = betahat)

ss  <- function(b) sum((y - exp(b))^2)
opt2 <- optimize(ss, interval = c(-5, 5), tol = 1e-10)
c(least.squares.betahat = opt2$minimum)

coef(glm(y ~ 1, family = poisson(link = "log")))

        n     total      ybar
13.000000 72.000000  5.538462
    betahat exp.betahat
   1.711717    5.538462
numerical.betahat       closed.form
         1.711717          1.711717
least.squares.betahat
             1.711717
(Intercept)
   1.711717

All four routes agree on \(\widehat{\beta} = 1.7117\), and \(e^{\widehat{\beta}} = 5.538\) recovers \(\widehat{\theta}\), confirming (b) and (c). The intercept-only Poisson generalized linear model with log link, which is what the last line fits, returns the identical estimate – this exercise is the simplest possible instance of the machinery developed in the rest of the book.

Part (a): The group size is \(n_i = (\text{females})_i + (\text{males})_i\) and the proportion female is \(p_i = y_i/n_i\).

library(dobson)
data(moths)
moths$n <- moths$females + moths$males
moths$p <- moths$females / moths$n
print(round(moths[, c("group", "females", "n", "p")], 3), row.names = FALSE)
c(total.females = sum(moths$females), total.progeny = sum(moths$n),
  theta.hat = sum(moths$females) / sum(moths$n))

 group females  n     p
     1      18 29 0.621
     2      31 53 0.585
     3      34 61 0.557
     4      33 62 0.532
     5      27 51 0.529
     6      33 62 0.532
     7      28 53 0.528
     8      23 49 0.469
     9      33 71 0.465
    10      12 26 0.462
    11      19 42 0.452
    12      25 56 0.446
    13      14 34 0.412
    14       4 10 0.400
    15      22 56 0.393
    16       7 19 0.368
total.females total.progeny     theta.hat
  363.0000000   734.0000000     0.4945504

The proportions range from \(0.368\) to \(0.621\), scattered around one half. The smallest groups carry the least information: group \(14\) has only \(10\) progeny, so its proportion of \(0.400\) is one female away from \(0.500\).

Part (b): The likelihood function for the whole sample is the product over the \(16\) independent groups,

\[ L(\theta; \mathbf{y}) = \prod_{i=1}^{16}\binom{n_i}{y_i}\theta^{y_i}(1 - \theta)^{n_i - y_i}, \]

so the log-likelihood is

\[ l(\theta; \mathbf{y}) = \sum_{i=1}^{16}\log\binom{n_i}{y_i} + \left( \sum y_i \right)\log\theta + \left( \sum n_i - \sum y_i \right)\log(1 - \theta). \]

The binomial coefficients do not involve \(\theta\) and so play no part in the optimisation, exactly as the \(\sum \log y_i!\) term did not in Section 1.6.5. Write \(N = \sum n_i\) and \(T = \sum y_i\). Then, following Equation (1.9),

\[ \frac{dl}{d\theta} = \frac{T}{\theta} - \frac{N - T}{1 - \theta} = \frac{T(1 - \theta) - (N - T)\theta}{\theta(1 - \theta)} = \frac{T - N\theta}{\theta(1 - \theta)}. \]

Setting the numerator to zero gives

\[ \widehat{\theta} = \frac{T}{N} = \frac{\sum_{i=1}^{16} y_i}{\sum_{i=1}^{16} n_i}. \]

For the maximum check,

\[ \frac{d^2l}{d\theta^2} = -\frac{T}{\theta^2} - \frac{N - T}{(1 - \theta)^2} < 0 \]

for all \(\theta \in (0,1)\) whenever \(0 < T < N\), so \(l\) is strictly concave and \(\widehat{\theta}\) is the unique interior maximum. From the output above, \(T = 363\) and \(N = 734\), so

\[ \widehat{\theta} = \frac{363}{734} = 0.4946. \]

Note that this is the pooled proportion, not the average of the \(16\) group proportions. The estimator weights each group by its size \(n_i\), which is correct because larger groups carry more information about \(\theta\).

Part (c): Three numerical routes, none of which uses the closed form: direct maximisation of \(l^*(\theta) = \sum [ y_i \log\theta + (n_i - y_i)\log(1 - \theta) ]\) with optimize, a bisection search in the spirit of Table 1.3, and the intercept-only binomial generalized linear model.

y <- moths$females; n <- moths$n

lstar <- function(th) sum(y * log(th) + (n - y) * log(1 - th))
opt <- optimize(lstar, interval = c(0.001, 0.999), maximum = TRUE, tol = 1e-12)
c(optimize = opt$maximum, closed.form = sum(y) / sum(n))

lo <- 0.30; hi <- 0.70
for (k in 1:20) {
  mid <- (lo + hi) / 2
  if (lstar(mid + 1e-8) > lstar(mid - 1e-8)) lo <- mid else hi <- mid
}
c(bisection = (lo + hi) / 2)

fit <- glm(cbind(females, males) ~ 1, family = binomial, data = moths)
c(glm.theta = plogis(coef(fit)))

   optimize closed.form
  0.4945504   0.4945504
bisection
0.4945505
glm.theta.(Intercept)
            0.4945504

All three numerical answers agree with the calculus answer of (b) to at least six decimal places; the bisection differs only in the seventh, because \(20\) halvings of an interval of width \(0.4\) leave a resolution of about \(4 \times 10^{-7}\). This is the same bisection idea used for the cyclone data in Section 1.6.5, and it makes the same point: when the closed form is unavailable – which is the normal situation for the generalized linear models of later chapters – a numerical search recovers the estimate without any loss of accuracy that matters.

Finally, inference and adequacy. The estimate is only useful with a standard error and a check that the single-\(\theta\) model is tenable.

th <- sum(y) / sum(n)
se <- sqrt(th * (1 - th) / sum(n))
c(theta.hat = th, se = se, lower = th - 1.96 * se, upper = th + 1.96 * se)
c(z.vs.half = (th - 0.5) / se, p.value = 2 * pnorm(-abs((th - 0.5) / se)))
c(resid.dev = deviance(fit), df = df.residual(fit),
  p.het = pchisq(deviance(fit), df.residual(fit), lower.tail = FALSE))

 theta.hat         se      lower      upper
0.49455041 0.01845424 0.45838010 0.53072072
 z.vs.half    p.value
-0.2953029  0.7677625
 resid.dev         df      p.het
11.9288136 15.0000000  0.6844091

The estimated sex ratio is \(\widehat{\theta} = 0.495\) with standard error \(0.018\) and approximate \(95\%\) confidence interval \((0.458, 0.531)\). This comfortably contains \(0.5\): the test of \(H_0: \theta = 0.5\) gives \(z = -0.30\), \(p = 0.77\), so there is no evidence that light brown apple moth progeny depart from an even sex ratio, which is what Mendelian sex determination predicts.

The residual deviance of \(11.93\) on \(15\) degrees of freedom (the \(16\) groups less the one fitted parameter) tests whether a common \(\theta\) suffices, against the alternative that each mother has her own sex ratio. With \(p = 0.68\) there is no evidence of heterogeneity between mothers – the spread of the \(16\) proportions in (a) is no more than binomial sampling variation would produce. The figure below makes this visual: the observed proportions fan out for small \(n_i\) and tighten for large \(n_i\) in the manner the binomial variance \(\theta(1-\theta)/n_i\) requires, and all \(16\) sit inside the pointwise \(95\%\) band.

p  <- y / n
th <- sum(y) / sum(n)
par(mfrow = c(1, 2))
plot(n, p, pch = 19, ylim = c(0.3, 0.7), xlab = "progeny in group (n_i)",
     ylab = "proportion female", main = "Proportions vs group size")
abline(h = th, lty = 2); abline(h = 0.5, lty = 3, col = "grey40")
curve(th + 1.96 * sqrt(th * (1 - th) / x), add = TRUE, col = "grey50")
curve(th - 1.96 * sqrt(th * (1 - th) / x), add = TRUE, col = "grey50")
legend("topright", c("MLE 0.4946", "0.5", "95% binomial band"),
       lty = c(2, 3, 1), col = c("black", "grey40", "grey50"), bty = "n")
tt <- seq(0.30, 0.70, length.out = 400)
ls <- sapply(tt, function(t) sum(y * log(t) + (n - y) * log(1 - t)))
plot(tt, ls, type = "l", xlab = expression(theta), ylab = "log-likelihood (l*)",
     main = "Log-likelihood function")
abline(v = th, lty = 2)

The right-hand panel is the analogue of Figure 1.2 for these data: a single smooth peak at \(\widehat{\theta} = 0.4946\), sharply curved, which is why the standard error is small. The curvature at the maximum is \(-d^2l/d\theta^2 = N/[\theta(1-\theta)] = 2936.3\), and \(1/\sqrt{2936.3} = 0.01845\) reproduces the standard error above – the connection between likelihood curvature and precision that Chapter 5 develops formally.

The data ship in the dobson package as plants, with one column per group; I stack them once and reuse the stacked frame throughout.

library(dobson)
data(plants)
plant <- data.frame(
  weight = c(plants$treatment, plants$control),
  group  = factor(rep(c("treatment", "control"), each = 20),
                  levels = c("treatment", "control")))
stats <- function(v) c(n = length(v), mean = mean(v), median = median(v),
                       sd = sd(v), min = min(v), max = max(v))
round(t(sapply(split(plant$weight, plant$group), stats)), 3)

           n  mean median    sd  min  max
treatment 20 4.860  4.850 0.791 3.48 6.34
control   20 4.726  4.635 0.864 3.05 6.11

(a) Exploratory analysis. Following the checklist of Section 2.3.1: the response is continuous, the explanatory variable is a two-level nominal factor, and the question is whether the two group distributions differ in location.

stem(plants$treatment)
stem(plants$control)

The decimal point is at the |

3 | 568
4 | 234477899
5 | 024589
6 | 03


The decimal point is at the |

3 | 1679
4 | 0125567
5 | 1223666
6 | 11

par(mfrow = c(2, 2))
stripchart(weight ~ group, data = plant, method = "stack", pch = 16,
           xlab = "dried weight (g)", main = "Dot plot by group")
boxplot(weight ~ group, data = plant, horizontal = TRUE,
        xlab = "dried weight (g)", main = "Box plots")
for (g in levels(plant$group)) {
  qqnorm(plant$weight[plant$group == g], main = paste("Normal Q-Q:", g))
  qqline(plant$weight[plant$group == g])
}

var.test(weight ~ group, data = plant)

	F test to compare two variances

data:  weight by group
F = 0.83891, num df = 19, denom df = 19, p-value = 0.7057
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.332052 2.119473
sample estimates:
ratio of variances
         0.8389132

All values lie between \(3.05\) and \(6.34\) g, which is plausible for dried plant material, so there are no obvious data-entry errors. Both distributions are single-peaked and roughly symmetric; the Normal probability plots are close to straight, so the Normality assumption is credible. The two spreads are very similar (\(s_1 = 0.791\), \(s_2 = 0.864\); the variance-ratio test gives \(F_{19,19} = 0.839\), \(p = 0.71\)), which supports the common-\(\sigma^2\) assumption built into both models. The treatment mean exceeds the control mean by only \(4.860 - 4.726 = 0.134\) g, an eighth of a standard deviation, and the two dot plots overlap almost completely. The exploratory analysis therefore suggests any treatment effect is small relative to the plant-to-plant variation, and we should not expect a significant test.

(b) The unpaired t-test. Since the model assumes a common \(\sigma^2\), the pooled-variance (equal-variance) two-sample \(t\)-test is the one that matches it.

t.test(weight ~ group, data = plant, var.equal = TRUE)

	Two Sample t-test

data:  weight by group
t = 0.50985, df = 38, p-value = 0.6131
alternative hypothesis: true difference in means between group treatment and group control is not equal to 0
95 percent confidence interval:
 -0.3965733  0.6635733
sample estimates:
mean in group treatment   mean in group control
                 4.8600                  4.7265

The estimated difference is \(\overline{y}_1 - \overline{y}_2 = 0.1335\) g with \(t = 0.510\) on \(38\) degrees of freedom and \(p = 0.61\): no evidence against \(\mathrm{H}_0\). The 95% confidence interval \((-0.397, 0.664)\) contains zero. Its width matters more than its \(p\)-value: the data are consistent with the enriched environment reducing mean dry weight by up to about \(0.40\) g or raising it by up to about \(0.66\) g. With \(20\) plants per group the experiment simply cannot resolve effects of this size, so “no significant difference” here means “inconclusive”, not “no effect”.

(c) Estimation. Write \(N = 40\) for the total number of observations. Under \(\mathrm{H}_1\) the log-likelihood, with \(\sigma^2\) a known constant as in Example 2.2.2, is \[ \ell_1(\mu_1,\mu_2;\mathbf{y}) = -\tfrac{1}{2}N\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{j=1}^{2}\sum_{k=1}^{20}(y_{jk}-\mu_j)^2 . \] The least squares criterion is \[ S_1 = \sum_{j=1}^{2}\sum_{k=1}^{20}(y_{jk}-\mu_j)^2 , \] so \(\ell_1 = -\tfrac{1}{2}N\log(2\pi\sigma^2) - S_1/(2\sigma^2)\) is a strictly decreasing linear function of \(S_1\). Maximizing \(\ell_1\) and minimizing \(S_1\) therefore give the same equations, exactly as for Equations (2.8) and (2.10) in Example 2.2.2. Differentiating, \[ \frac{\partial \ell_1}{\partial \mu_j} = \frac{1}{\sigma^2}\sum_{k=1}^{20}(y_{jk}-\mu_j) = 0 \quad\Longrightarrow\quad \widehat{\mu}_j = \frac{1}{20}\sum_{k=1}^{20} Y_{jk} = \overline{Y}_j , \] for \(j = 1, 2\). The second derivative is \(-20/\sigma^2 < 0\), so this is a maximum of \(\ell_1\) (equivalently \(\partial^2 S_1/\partial\mu_j^2 = 40 > 0\), a minimum of \(S_1\)).

Under \(\mathrm{H}_0\) the single mean \(\mu\) appears in all \(40\) terms, so \[ \frac{\partial \ell_0}{\partial \mu} = \frac{1}{\sigma^2}\sum_{j=1}^{2}\sum_{k=1}^{20}(y_{jk}-\mu) = 0 \quad\Longrightarrow\quad \widehat{\mu} = \frac{1}{40}\sum_{j=1}^{2}\sum_{k=1}^{20} Y_{jk} = \overline{Y}. \] So the maximum likelihood and least squares estimates coincide and are the obvious sample means: numerically \(\widehat{\mu}_1 = 4.8600\), \(\widehat{\mu}_2 = 4.7265\), \(\widehat{\mu} = 4.79325\).

(d) The minimized criteria. Substituting the estimates from (c) back into the criteria gives the minimum values directly: \[ \widehat{S}_1 = \sum_{j=1}^{2}\sum_{k=1}^{20}(Y_{jk}-\widehat{\mu}_j)^2 = \sum\sum (Y_{jk}-\overline{Y}_j)^2, \qquad \widehat{S}_0 = \sum\sum (Y_{jk}-\overline{Y})^2 . \] That these are minima and not merely stationary values follows because \(S_1\) is a positive-definite quadratic in \((\mu_1,\mu_2)\). Explicitly, for any \(\mu_j\), \[ \sum_k (Y_{jk}-\mu_j)^2 = \sum_k (Y_{jk}-\overline{Y}_j)^2 + 20(\overline{Y}_j - \mu_j)^2 \geq \sum_k (Y_{jk}-\overline{Y}_j)^2 , \] the cross-product vanishing because \(\sum_k (Y_{jk}-\overline{Y}_j) = 0\); equality holds only at \(\mu_j = \overline{Y}_j\). The same identity with a single mean over all \(40\) observations gives \(\widehat{S}_0\).

(e) Chi-squared distributions. Two printing slips in the exercise as set should be noted: in the first display the subtracted term is summed over \(j = 1, 2\) (not \(k = 1,\ldots,20\)), and in the third display the term is \(\frac{40}{\sigma^2}(\overline{Y}-\mu)^2\), the sum over \(j\) being spurious. With those corrections the identities are exactly Exercise 1.4(b) applied group-by-group.

Exercise 1.4(b) states that for \(n\) independent \(\mathrm{N}(\mu,\sigma^2)\) variables, \(\sum_i (Y_i-\overline{Y})^2 = \sum_i (Y_i-\mu)^2 - n(\overline{Y}-\mu)^2\). Applying it within group \(j\) with \(n = 20\) and mean \(\mu_j\), then summing over \(j\) and dividing by \(\sigma^2\), \[ \frac{1}{\sigma^2}\widehat{S}_1 = \frac{1}{\sigma^2}\sum_{j=1}^{2}\sum_{k=1}^{20}(Y_{jk}-\mu_j)^2 - \frac{20}{\sigma^2}\sum_{j=1}^{2}(\overline{Y}_j-\mu_j)^2 . \] If \(\mathrm{H}_1\) is true then \((Y_{jk}-\mu_j)/\sigma \sim \mathrm{N}(0,1)\) independently, so the first term is \(\chi^2(40)\); and \(\overline{Y}_j \sim \mathrm{N}(\mu_j, \sigma^2/20)\) independently across \(j\), so \(\sqrt{20}(\overline{Y}_j-\mu_j)/\sigma \sim \mathrm{N}(0,1)\) and the second term is \(\chi^2(2)\). By Exercise 1.4(c) the sample mean and the within-group sum of squares are independent, so \(\widehat{S}_1/\sigma^2\) is independent of the second term. Writing \(\chi^2(40) = \widehat{S}_1/\sigma^2 + \chi^2(2)\) as a sum of two independent terms and comparing moment generating functions gives \[ \frac{1}{\sigma^2}\widehat{S}_1 \sim \chi^2(38), \] which is \(N\) minus the number of parameters estimated, \(40 - 2 = 38\).

Under \(\mathrm{H}_0\) all \(40\) observations have the same mean \(\mu\), so applying Exercise 1.4(b) once with \(n = 40\), \[ \frac{1}{\sigma^2}\widehat{S}_0 = \frac{1}{\sigma^2}\sum_{j=1}^{2}\sum_{k=1}^{20}(Y_{jk}-\mu)^2 - \frac{40}{\sigma^2}(\overline{Y}-\mu)^2 , \] where now the first term is \(\chi^2(40)\) and the second, since \(\overline{Y}\sim \mathrm{N}(\mu,\sigma^2/40)\), is \(\chi^2(1)\) and independent of \(\widehat{S}_0\). Hence \[ \frac{1}{\sigma^2}\widehat{S}_0 \sim \chi^2(39) = \chi^2(40-1). \]

(f) The F-statistic. Repeat the argument of Example 2.2.2. Expanding \(Y_{jk}-\overline{Y} = (Y_{jk}-\overline{Y}_j) + (\overline{Y}_j-\overline{Y})\) and using \(\sum_k (Y_{jk}-\overline{Y}_j) = 0\), \[ \widehat{S}_0 - \widehat{S}_1 = 20\sum_{j=1}^{2}(\overline{Y}_j - \overline{Y})^2 = 10(\overline{Y}_1-\overline{Y}_2)^2 , \] the last step because with equal group sizes \(\overline{Y} = (\overline{Y}_1+\overline{Y}_2)/2\). Under \(\mathrm{H}_0\), \(\overline{Y}_1-\overline{Y}_2 \sim \mathrm{N}(0, \sigma^2/10)\), so \((\widehat{S}_0-\widehat{S}_1)/\sigma^2 \sim \chi^2(1)\), consistent with the general statement \((\widehat{S}_0-\widehat{S}_1)/\sigma^2 \sim \chi^2(J-1)\) in Section 2.2.2 with \(J = 2\). Moreover \(\widehat{S}_0 - \widehat{S}_1\) is a function of the group means alone while \(\widehat{S}_1\) is a function of the within-group deviations alone, and these are independent for Normal data; so the two chi-squared variables are independent. Since \(\sigma^2\) is unknown in practice we eliminate it by taking the ratio of the two independent chi-squared variables, each divided by its degrees of freedom: \[ F = \frac{(\widehat{S}_0-\widehat{S}_1)/\sigma^2}{1} \bigg/ \frac{\widehat{S}_1/\sigma^2}{38} = \frac{\widehat{S}_0-\widehat{S}_1}{\widehat{S}_1/38} \sim F(1,38) \] if \(\mathrm{H}_0\) is true, by the definition of the \(F\)-distribution in Section 1.4.4. If \(\mathrm{H}_0\) is false then \(\mathrm{E}(\overline{Y}_1-\overline{Y}_2) = \mu_1-\mu_2 \neq 0\), so \((\widehat{S}_0-\widehat{S}_1)/\sigma^2\) has a non-central chi-squared distribution with non-centrality parameter \(10(\mu_1-\mu_2)^2/\sigma^2\) and \(F\) has a non-central \(F\)-distribution, shifted to the right of the central one (Figure 2.5).

(g) Testing. Computing the two minima and the statistic directly from their definitions:

y <- plant$weight
S0 <- sum((y - mean(y))^2)
S1 <- sum((y - ave(y, plant$group))^2)
Fstat <- (S0 - S1) / (S1 / 38)
c(S0 = S0, S1 = S1, F = Fstat, p = pf(Fstat, 1, 38, lower.tail = FALSE))

        S0         S1          F          p
26.2316775 26.0534550  0.2599446  0.6131068

The same numbers come out of the model-comparison machinery, which is a useful check that \(\widehat{S}_0\) and \(\widehat{S}_1\) are the residual sums of squares of the two nested linear models:

anova(lm(weight ~ 1, data = plant), lm(weight ~ group, data = plant))

Analysis of Variance Table

Model 1: weight ~ 1
Model 2: weight ~ group
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     39 26.232
2     38 26.053  1   0.17822 0.2599 0.6131

Fitting separate group means reduces the sum of squares from \(26.232\) to \(26.053\), an improvement of only \(0.178\), or \(0.7\)% of the total. With \(F = 0.260\) on \((1, 38)\) degrees of freedom, \(p = 0.61\), this is far below the \(95\)th percentile \(F_{0.95}(1,38) = 4.10\). There is no evidence against \(\mathrm{H}_0\); on grounds of simplicity the single-mean model is preferable, and the nutritional enrichment has no demonstrable effect on mean dry weight.

(h) Relation to the t-test. Section 1.4.4 records that if \(T \sim t(n)\) then \(T^2 \sim F(1,n)\). Here \(t = 0.50985\) and \(t^2 = 0.25994 = F\) to five decimal places, and the two \(p\)-values agree exactly at \(0.6131\) because the two-sided \(t\) test and the one-sided upper-tail \(F\) test reject on the same event \(\{|T| > c\}\). The two procedures are algebraically the same test: the least squares comparison of nested models with a one-degree-of-freedom difference is the unpaired \(t\)-test rewritten. The conclusions in (b) and (g) are therefore identical, but the \(t\)-test formulation is more informative in practice because it also delivers the signed estimate and the confidence interval, whereas \(F\) discards the sign of the effect.

(i) Residuals. Under \(\mathrm{H}_0\) the fitted value for every observation is \(\widehat{\mu} = \overline{y}\), so the raw residuals are \(y_{jk}-\overline{y}\); standardizing by \(s = \sqrt{\widehat{S}_0/39}\) gives residuals comparable with Figures 2.3 and 2.4 of Example 2.2.2.

r <- (y - mean(y)) / sqrt(S0 / 39)
par(mfrow = c(1, 3))
plot(as.integer(plant$group), r, xaxt = "n", xlab = "group",
     ylab = "standardised residual", main = "Residuals by group", pch = 16)
axis(1, at = 1:2, labels = levels(plant$group)); abline(h = 0, lty = 2)
plot(seq_along(r), r, xlab = "observation", ylab = "standardised residual",
     main = "Residuals in order", pch = 16); abline(h = 0, lty = 2)
qqnorm(r, main = "Normal Q-Q of residuals"); qqline(r)

round(c(mean = mean(r), sd = sd(r), min = min(r), max = max(r)), 3)
tapply(r, plant$group, mean)
shapiro.test(r)

  mean     sd    min    max
 0.000  1.000 -2.126  1.886

 treatment    control
 0.0813899 -0.0813899

	Shapiro-Wilk normality test

data:  r
W = 0.9844, p-value = 0.8457

The residuals have mean zero and unit standard deviation by construction. The Normal probability plot is close to a straight line and no residual lies beyond \(\pm 2.13\), well inside what \(40\) standard Normal draws would produce; the Shapiro-Wilk test gives \(p = 0.85\). The spread is visually the same in the two groups, supporting the constant-variance assumption. Crucially the residuals show no group structure: the two group means of the residuals are \(+0.081\) and \(-0.081\), tiny compared with a residual standard deviation of \(1\). That is the diagnostic restatement of (g) — the omitted group term would explain almost none of the residual variation, so the \(\mathrm{H}_0\) model is adequate for these data.

The data ship in the dobson package as waist, one row per man with columns before and after.

library(dobson)
data(waist)
w <- data.frame(weight = c(waist$before, waist$after),
                time = factor(rep(c("before", "after"), each = 20),
                              levels = c("before", "after")))
round(t(sapply(split(w$weight, w$time),
       function(v) c(n = length(v), mean = mean(v), sd = sd(v),
                     min = min(v), max = max(v)))), 3)

        n    mean     sd min   max
before 20 103.245 13.517  80 135.0
after  20 100.600 12.475  80 127.5

d <- waist$before - waist$after
par(mfrow = c(1, 3))
plot(waist$before, waist$after, pch = 16, xlab = "weight before (kg)",
     ylab = "weight after (kg)", main = "Paired weights")
abline(0, 1, lty = 2)
stripchart(d, method = "stack", pch = 16,
           xlab = "difference before - after (kg)",
           main = "Individual weight losses")
abline(v = 0, lty = 2)
qqnorm(d, main = "Normal Q-Q of differences"); qqline(d)

The scatter plot is the whole story: the points lie almost on the line \(y = x\), and nearly all of them sit just below it. The men differ enormously from each other (weights run from \(80\) to \(135\) kg), but each man changes only a little.

(a) The unpaired t-test. Treating the two sets of \(20\) weights as independent samples, as the exercise’s stated model does:

t.test(weight ~ time, data = w, var.equal = TRUE)

	Two Sample t-test

data:  weight by time
t = 0.64309, df = 38, p-value = 0.524
alternative hypothesis: true difference in means between group before and group after is not equal to 0
95 percent confidence interval:
 -5.68128 10.97128
sample estimates:
mean in group before  mean in group after
             103.245              100.600

The estimated mean loss is \(103.245 - 100.600 = 2.645\) kg, but \(t = 0.643\) on \(38\) degrees of freedom gives \(p = 0.52\) and the 95% confidence interval \((-5.68, 10.97)\) is very wide. On this analysis there is no evidence of a retained weight loss.

(b) The paired analysis. Work with the \(20\) differences \(D_k = Y_{1k} - Y_{2k}\). Under the stated assumptions \(D_k \sim \mathrm{N}(\delta, \sigma_D^2)\) with \(\delta = \mu_1 - \mu_2\), and the \(D_k\) are independent across men. The two competing models are \[ \mathrm{H}_0 : \mathrm{E}(D_k) = 0; \quad D_k \sim \mathrm{N}(0,\sigma_D^2), \qquad \mathrm{H}_1 : \mathrm{E}(D_k) = \delta; \quad D_k \sim \mathrm{N}(\delta,\sigma_D^2), \] for \(k = 1,\ldots,20\), with \(\sigma_D^2\) treated as a known constant. Exactly as in Exercise 2.1(c), maximizing the log-likelihood \[ \ell_1(\delta;\mathbf{d}) = -\tfrac{1}{2}\cdot 20\log(2\pi\sigma_D^2) - \frac{1}{2\sigma_D^2}\sum_{k=1}^{20}(d_k-\delta)^2 \] is equivalent to minimizing \(S_1 = \sum_k (d_k-\delta)^2\), and gives \(\widehat{\delta} = \overline{D}\). The minimized criteria are \[ \widehat{S}_0 = \sum_{k=1}^{20} D_k^2 \quad (\text{no parameters estimated}), \qquad \widehat{S}_1 = \sum_{k=1}^{20}(D_k-\overline{D})^2 \quad (\text{one parameter}). \] By Exercise 1.4(b) with \(n = 20\), \(\widehat{S}_1/\sigma_D^2 = \sum_k (D_k-\delta)^2/\sigma_D^2 - 20(\overline{D}-\delta)^2/\sigma_D^2\), a \(\chi^2(20)\) minus an independent \(\chi^2(1)\), so \(\widehat{S}_1/\sigma_D^2 \sim \chi^2(19)\) whether or not \(\mathrm{H}_0\) holds. If \(\mathrm{H}_0\) is true then \(\widehat{S}_0/\sigma_D^2 \sim \chi^2(20)\) and \(\widehat{S}_0 - \widehat{S}_1 = 20\overline{D}^2\) with \(\overline{D} \sim \mathrm{N}(0,\sigma_D^2/20)\), so \((\widehat{S}_0-\widehat{S}_1)/\sigma_D^2 \sim \chi^2(1)\) independently of \(\widehat{S}_1\). Hence \[ F = \frac{\widehat{S}_0-\widehat{S}_1}{\widehat{S}_1/19} \sim F(1,19) \quad \text{if } \mathrm{H}_0 \text{ is true.} \]

S0 <- sum(d^2); S1 <- sum((d - mean(d))^2)
Fstat <- (S0 - S1) / (S1 / 19)
c(S0 = S0, S1 = S1, diff = S0 - S1, F = Fstat,
  p = pf(Fstat, 1, 19, lower.tail = FALSE))

          S0           S1         diff            F            p
4.619100e+02 3.219895e+02 1.399205e+02 8.256448e+00 9.730463e-03

\(F = 8.256\) on \((1,19)\) degrees of freedom gives \(p = 0.0097\), well beyond the \(95\)th percentile \(F_{0.95}(1,19) = 4.38\), so \(\mathrm{H}_0\) is rejected. Equivalently, as a \(t\)-statistic:

t.test(waist$before, waist$after, paired = TRUE)

	Paired t-test

data:  waist$before and waist$after
t = 2.8734, df = 19, p-value = 0.00973
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.718348 4.571652
sample estimates:
mean difference
          2.645

Note \(t^2 = 2.8734^2 = 8.256 = F\), the same \(t\)-versus-\(F\) identity used in Exercise 2.1(h). On average the men retain a loss of \(2.645\) kg twelve months on, with 95% confidence interval \((0.72, 4.57)\) kg.

(c) Comparison. No — the conclusions are opposite. Both analyses estimate the same quantity, \(2.645\) kg, but the unpaired analysis attaches a standard error of \(4.11\) kg to it while the paired analysis attaches \(0.92\) kg, a four-and-a-half-fold reduction. The reason is the correlation between a man’s two weights:

c(cor = cor(waist$before, waist$after),
  se_unpaired = sqrt(var(waist$before)/20 + var(waist$after)/20),
  se_paired = sd(d)/sqrt(20))

      cor se_unpaired   se_paired
0.9529734   4.1129735   0.9205112

With \(\mathrm{corr}(Y_{1k},Y_{2k}) = \rho = 0.953\), \(\mathrm{var}(D_k) = 2\sigma^2(1-\rho)\) is only about \(5\)% of the \(2\sigma^2\) the unpaired analysis assumes. Pairing removes the between-man variation — which is the dominant source of scatter and is irrelevant to the question — and leaves only the within-man change. The unpaired test is not wrong arithmetically; it is answering the question with the wrong yardstick, and its failure to reject is a Type II error.

(d) Assumptions. Assumptions for (a): the \(40\) observations are mutually independent; both sets are Normally distributed; the two variances are equal; the sample is representative of the population of interest. The independence assumption is the fatal one, since the same \(20\) men are measured twice.

Assumptions for (b): the \(20\) differences \(D_k\) are mutually independent (reasonable — different men); each \(D_k\) is Normal with common variance \(\sigma_D^2\) (the Q-Q plot above is acceptably straight, Shapiro-Wilk \(p = 0.16\), though there are several exact zeros where a man’s weight was recorded unchanged); and the \(D_k\) share a common mean \(\delta\), i.e. the retained loss does not depend on initial weight. Nothing is assumed about the between-man distribution of weights, and no equal-variance assumption across times is needed.

shapiro.test(d)
summary(lm(d ~ waist$before))$coefficients

	Shapiro-Wilk normality test

data:  d
W = 0.93148, p-value = 0.1649

               Estimate Std. Error   t value   Pr(>|t|)
(Intercept)  -9.8004283  6.8611113 -1.428402 0.17029887
waist$before  0.1205427  0.0659201  1.828618 0.08407696

Regressing the loss on initial weight gives a slope of \(0.121\) kg per kg of initial weight, \(p = 0.084\): a heavier man tends to retain a slightly larger loss, but the evidence is weak and the common-\(\delta\) assumption is not clearly contradicted. This is worth flagging rather than acting on with \(n = 20\).

The paired analysis in (b) is the appropriate one. The design is a before-and-after study on the same subjects, so the data are paired by construction and the independence assumed in (a) does not hold. Conclusion: twelve months after the program the men retain a mean weight loss of about \(2.6\) kg (95% CI \(0.7\) to \(4.6\) kg). Two caveats on interpretation remain, neither statistical: there is no control group, so regression to the mean and secular trends cannot be separated from the program effect; and the individual results are heterogeneous — \(14\) men lost weight, \(4\) were unchanged to the recorded precision, and \(2\) gained.

Throughout, fix a group \(j\) and write \[ S_{xx} = \sum_{k=1}^{K}(x_{jk}-\overline{x}_j)^2 = \sum_{k=1}^{K}x_{jk}^2 - K\overline{x}_j^2, \qquad S_{xY} = \sum_{k=1}^{K}(x_{jk}-\overline{x}_j)(Y_{jk}-\overline{Y}_j), \] \[ S_{YY} = \sum_{k=1}^{K}(Y_{jk}-\overline{Y}_j)^2 , \] so that the least squares estimates of Section 2.2.2 are \(b_j = S_{xY}/S_{xx}\) and \(a_j = \overline{Y}_j - b_j\overline{x}_j\). Assume \(S_{xx} > 0\), i.e. the gestational ages within a group are not all equal.

(1) The algebraic identity. Substituting \(a_j = \overline{Y}_j - b_j\overline{x}_j\) removes the intercept: \[ Y_{jk}-a_j-b_j x_{jk} = (Y_{jk}-\overline{Y}_j) - b_j(x_{jk}-\overline{x}_j), \] whence the group-\(j\) contribution to \(\widehat{S}_1\) is \[ \widehat{S}_{1j} = S_{YY} - 2b_j S_{xY} + b_j^2 S_{xx} = S_{YY} - b_j^2 S_{xx}, \tag{i} \] using \(S_{xY} = b_j S_{xx}\) twice.

Now work on the right-hand side. Put \(e_{jk} = Y_{jk} - (\alpha_j + \beta_j x_{jk})\), so that \(\overline{e}_j = \overline{Y}_j - \alpha_j - \beta_j\overline{x}_j\) is precisely the quantity squared in the second term. Exercise 1.4(b), applied to \(e_{j1},\ldots,e_{jK}\) with \(n = K\), gives \[ \sum_{k}e_{jk}^2 - K\overline{e}_j^2 = \sum_k (e_{jk}-\overline{e}_j)^2 . \] Since \(e_{jk}-\overline{e}_j = (Y_{jk}-\overline{Y}_j) - \beta_j(x_{jk}-\overline{x}_j)\), \[ \sum_k (e_{jk}-\overline{e}_j)^2 = S_{YY} - 2\beta_j S_{xY} + \beta_j^2 S_{xx} = S_{YY} - 2\beta_j b_j S_{xx} + \beta_j^2 S_{xx}. \] Subtracting the third term, and expanding \((b_j-\beta_j)^2 S_{xx} = b_j^2 S_{xx} - 2b_j\beta_j S_{xx} + \beta_j^2 S_{xx}\), everything involving \(\beta_j\) cancels: \[ \sum_k e_{jk}^2 - K\overline{e}_j^2 - (b_j-\beta_j)^2 S_{xx} = S_{YY} - b_j^2 S_{xx} = \widehat{S}_{1j} \] by (i). Summing over \(j = 1,\ldots,J\) gives the stated identity. Note it is an algebraic identity: it holds for any numbers \(\alpha_j,\beta_j\), not only the true parameter values. A numerical check on the birthweight data, using both the true-ish fitted values from Table 2.5 and a pair of arbitrary values, confirms this:

library(dobson)
data(birthweight)
bw <- data.frame(y   = c(birthweight[[2]], birthweight[[4]]),
                 x   = c(birthweight[[1]], birthweight[[3]]),
                 sex = rep(c("boys", "girls"), each = 12))
K <- 12
check <- function(alpha, beta) {
  S1 <- 0; A <- 0; B <- 0; C <- 0
  for (j in unique(bw$sex)) {
    yj <- bw$y[bw$sex == j]; xj <- bw$x[bw$sex == j]
    i  <- match(j, unique(bw$sex))
    bj <- cov(xj, yj) / var(xj); aj <- mean(yj) - bj * mean(xj)
    Sxx <- sum(xj^2) - K * mean(xj)^2
    S1 <- S1 + sum((yj - aj - bj * xj)^2)
    A  <- A + sum((yj - (alpha[i] + beta[i] * xj))^2)
    B  <- B + K * (mean(yj) - alpha[i] - beta[i] * mean(xj))^2
    C  <- C + (bj - beta[i])^2 * Sxx
  }
  c(S1 = S1, rhs = A - B - C, A = A, B = B, C = C)
}
round(check(c(-1268.672, -2141.667), c(111.983, 130.400)), 4)
round(check(c(0, 500), c(100, 90)), 4)

         S1         rhs           A           B           C
652424.5218 652424.5218 652424.5230      0.0011      0.0000

         S1         rhs           A           B           C
  652424.52   652424.52 22475004.00 21757861.67    64717.81

The value \(\widehat{S}_1 = 652424.5\) agrees with Table 2.5, and rhs reproduces it in both cases even though the three components differ by orders of magnitude.

(2) The distributions. \(Y_{jk} \sim \mathrm{N}(\alpha_j+\beta_j x_{jk},\sigma^2)\) is the model assumption in (2.7) itself. Both \(\overline{Y}_j\) and \(b_j\) are linear combinations of independent Normal variables, hence Normal, and it remains to compute the two moments in each case. For the mean, \[ \mathrm{E}(\overline{Y}_j) = \frac{1}{K}\sum_k (\alpha_j+\beta_j x_{jk}) = \alpha_j + \beta_j \overline{x}_j, \qquad \mathrm{var}(\overline{Y}_j) = \frac{1}{K^2}\sum_k \sigma^2 = \frac{\sigma^2}{K}. \] For the slope, write it in the equivalent centred form \[ b_j = \frac{S_{xY}}{S_{xx}} = \frac{1}{S_{xx}}\sum_k (x_{jk}-\overline{x}_j)Y_{jk}, \] the replacement of \(Y_{jk}-\overline{Y}_j\) by \(Y_{jk}\) being legitimate because \(\sum_k (x_{jk}-\overline{x}_j) = 0\). Then, using \(\sum_k (x_{jk}-\overline{x}_j)x_{jk} = S_{xx}\), \[ \mathrm{E}(b_j) = \frac{1}{S_{xx}}\sum_k (x_{jk}-\overline{x}_j)(\alpha_j+\beta_j x_{jk}) = \frac{\beta_j S_{xx}}{S_{xx}} = \beta_j , \] \[ \mathrm{var}(b_j) = \frac{\sigma^2}{S_{xx}^2}\sum_k (x_{jk}-\overline{x}_j)^2 = \frac{\sigma^2}{S_{xx}} = \frac{\sigma^2}{\sum_k x_{jk}^2 - K\overline{x}_j^2}, \] as required.

(3) The independence. The exercise’s phrase “the random variables \(Y_{jk}\), \(\overline{Y}_j\) and \(b_j\) are all independent” cannot be taken literally, since \(\overline{Y}_j\) is a function of the \(Y_{jk}\) in its own group. What the degrees-of-freedom argument of Section 2.2.2 actually needs, and what is true, is that within each group the three terms of the decomposition are mutually independent; across groups everything is independent because the groups involve disjoint sets of \(Y_{jk}\).

Rearranged, the identity reads \[ \frac{1}{\sigma^2}\sum_j\sum_k e_{jk}^2 = \frac{\widehat{S}_1}{\sigma^2} + \frac{K}{\sigma^2}\sum_j \overline{e}_j^2 + \frac{1}{\sigma^2}\sum_j (b_j-\beta_j)^2 S_{xx,j}. \] Fix \(j\) and regard \(\mathbf{e}_j = (e_{j1},\ldots,e_{jK})^{\mathrm{T}} \sim \mathrm{N}_K(\mathbf{0},\sigma^2 \mathbf{I})\). Let \(V_j\) be the two-dimensional subspace of \(\mathbb{R}^K\) spanned by the orthogonal vectors \(\mathbf{u}_1 = \mathbf{1}/\sqrt{K}\) and \(\mathbf{u}_2 = (\mathbf{x}_j - \overline{x}_j\mathbf{1})/\sqrt{S_{xx}}\). Then \[ \mathbf{u}_1^{\mathrm{T}}\mathbf{e}_j = \sqrt{K}\,\overline{e}_j, \qquad \mathbf{u}_2^{\mathrm{T}}\mathbf{e}_j = \frac{S_{xY}-\beta_j S_{xx}}{\sqrt{S_{xx}}} = (b_j-\beta_j)\sqrt{S_{xx}}, \] so the three terms of the decomposition are exactly \(\|(\mathbf{I}-\mathbf{P}_j)\mathbf{e}_j\|^2\), \((\mathbf{u}_1^{\mathrm{T}}\mathbf{e}_j)^2\) and \((\mathbf{u}_2^{\mathrm{T}}\mathbf{e}_j)^2\), where \(\mathbf{P}_j\) projects onto \(V_j\). Because \(\mathbf{e}_j\) is spherically Normal, its components along mutually orthogonal directions — and its component in the orthogonal complement — are independent. Dividing by \(\sigma^2\), \[ \frac{\widehat{S}_{1j}}{\sigma^2}\sim\chi^2(K-2), \qquad \frac{K\overline{e}_j^2}{\sigma^2}\sim\chi^2(1), \qquad \frac{(b_j-\beta_j)^2 S_{xx}}{\sigma^2}\sim\chi^2(1), \] independently, and these add to \(\chi^2(K)\) as they must. In particular \(\overline{Y}_j\) and \(b_j\) are uncorrelated, \[ \mathrm{cov}(\overline{Y}_j, b_j) = \frac{\sigma^2}{K S_{xx}}\sum_k (x_{jk}-\overline{x}_j) = 0, \] and jointly Normal, hence independent — which is the substantive content of the exercise’s claim. Summing the independent contributions over \(j = 1,\ldots,J\) gives \(\widehat{S}_1/\sigma^2 \sim \chi^2(JK-2J)\), the result quoted in Section 2.2.2, the degrees of freedom being the number of observations \(JK\) minus the number of parameters estimated, \(2J\).

The task is to identify the link. The mean is given as a function of the linear component rather than the other way round, so invert it: applying \(\exp\) to both sides of \(\mathrm{E}(Y_i) = \ln(\beta_0+\beta_1 x_i+\beta_2 x_i^2)\) gives \[ \exp[\mathrm{E}(Y_i)] = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 . \] This is of the form \(g[\mathrm{E}(Y_i)] = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}\) required by Section 2.4, with link function \[ g(\mu) = e^{\mu} \] and with two explanatory terms \(x_{i1} = x_i\) and \(x_{i2} = x_i^2\). Note that although the model is quadratic in \(x\), it is still linear in the parameters, which is all that (2.13) requires: \(x\) and \(x^2\) are simply two columns of measured constants.

There are \(N = 6\) observations and \(p = 3\) parameters, so \[ \mathbf{y} = \begin{bmatrix} 3.15 \\ 4.85 \\ 6.50 \\ 7.20 \\ 8.25 \\ 16.50 \end{bmatrix}, \qquad \boldsymbol{\beta} = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix}, \qquad \mathbf{X} = \begin{bmatrix} 1 & 1.0 & 1.00 \\ 1 & 1.2 & 1.44 \\ 1 & 1.4 & 1.96 \\ 1 & 1.6 & 2.56 \\ 1 & 1.8 & 3.24 \\ 1 & 2.0 & 4.00 \end{bmatrix}, \] and the model in the form of (2.13) is \[ g[\mathrm{E}(\mathbf{y})] = \begin{bmatrix} \exp[\mathrm{E}(Y_1)] \\ \vdots \\ \exp[\mathrm{E}(Y_6)] \end{bmatrix} = \mathbf{X}\boldsymbol{\beta}. \]

The design matrix is just a Vandermonde-style construction:

x <- c(1.0, 1.2, 1.4, 1.6, 1.8, 2.0)
y <- c(3.15, 4.85, 6.50, 7.20, 8.25, 16.50)
X <- cbind(1, x, x^2)
colnames(X) <- c("beta0", "beta1", "beta2")
X

     beta0 beta1 beta2
[1,]     1   1.0  1.00
[2,]     1   1.2  1.44
[3,]     1   1.4  1.96
[4,]     1   1.6  2.56
[5,]     1   1.8  3.24
[6,]     1   2.0  4.00

Two remarks that the exercise does not ask for but that matter if the model is ever fitted. First, the parameter space is restricted: the logarithm requires \(\beta_0+\beta_1 x+\beta_2 x^2 > 0\) over the range of \(x\), so unlike an ordinary linear model \(\boldsymbol{\beta}\) is not free to roam \(\mathbb{R}^3\). Second, this link makes the fit numerically awkward here, because \(\exp(y)\) ranges from \(\exp(3.15) \approx 23\) to \(\exp(16.5) \approx 1.5\times 10^7\) and a quadratic cannot follow that. Fitting by least squares with a custom link confirms it:

explink <- structure(list(
  linkfun  = function(mu)  exp(mu),
  linkinv  = function(eta) log(pmax(eta, .Machine$double.eps)),
  mu.eta   = function(eta) 1 / pmax(eta, .Machine$double.eps),
  valideta = function(eta) all(eta > 0),
  name     = "exp"), class = "link-glm")
fit <- glm(y ~ x + I(x^2), family = gaussian(link = explink),
           start = c(exp(3), 1, 1))
round(coef(fit), 2)
round(cbind(y = y, fitted = fitted(fit)), 3)

(Intercept)           x      I(x^2)
   60843.03  -111879.21    51059.77

      y fitted
1  3.15  3.161
2  4.85  4.737
3  6.50  8.364
4  7.20  9.437
5  8.25 10.122
6 16.50 10.629

The residual sum of squares is \(46.47\), and a multi-start Nelder-Mead search over the constrained parameter space finds no better optimum, so this is the least squares solution rather than a convergence failure. It is simply a bad fit: the last observation \(y = 16.50\) is far above what the model can reach. The exercise is about notation, but it also illustrates that writing a model in the form of (2.13) says nothing about whether that model suits the data.

There are \(N = 6\) observations but the unconstrained expression \(\mu+\alpha_j+\beta_k\) carries \(1+2+3 = 6\) parameters, and they are not all identifiable: adding a constant to \(\mu\) and subtracting it from every \(\alpha_j\) leaves all six means unchanged. This is the situation flagged in Example 2.4.3(b) — “too many parameters” — and the sum-to-zero constraints of Example 2.4.3(d) are the fix. They remove one parameter from each factor, leaving \(p = 4\) free parameters \[ \boldsymbol{\beta} = \begin{bmatrix} \mu \\ \alpha_1 \\ \beta_1 \\ \beta_2 \end{bmatrix}. \]

Following the hint, substitute \(\alpha_2 = -\alpha_1\) and \(\beta_3 = -\beta_1-\beta_2\) into the six means: \[ \mathrm{E}(Y_{11}) = \mu+\alpha_1+\beta_1, \qquad \mathrm{E}(Y_{12}) = \mu+\alpha_1+\beta_2, \qquad \mathrm{E}(Y_{13}) = \mu+\alpha_1-\beta_1-\beta_2, \] \[ \mathrm{E}(Y_{21}) = \mu-\alpha_1+\beta_1, \qquad \mathrm{E}(Y_{22}) = \mu-\alpha_1+\beta_2, \qquad \mathrm{E}(Y_{23}) = \mu-\alpha_1-\beta_1-\beta_2 . \] Reading off the coefficients of \((\mu,\alpha_1,\beta_1,\beta_2)\) row by row gives the model in the form of (2.13), with \(g\) the identity function: \[ \mathrm{E}(\mathbf{y}) = \mathbf{X}\boldsymbol{\beta}, \qquad \mathbf{y} = \begin{bmatrix} Y_{11} \\ Y_{12} \\ Y_{13} \\ Y_{21} \\ Y_{22} \\ Y_{23} \end{bmatrix}, \qquad \mathbf{X} = \begin{bmatrix} 1 & 1 & 1 & 0 \\ 1 & 1 & 0 & 1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & 1 & 0 \\ 1 & -1 & 0 & 1 \\ 1 & -1 & -1 & -1 \end{bmatrix}. \]

The second column is the sum-to-zero (effect) coding for the two-level factor, exactly the rows \([\,1\ \ 1\,]\) and \([\,1\ \ -1\,]\) of Example 2.4.3(d) repeated three times each; the last two columns are the corresponding coding for the three-level factor, with the third level coded as minus the sum of the other two so that the level effects sum to zero.

This is what R produces with contr.sum applied to both factors:

d <- expand.grid(k = factor(1:3), j = factor(1:2))
d <- d[order(d$j, d$k), ]
X <- model.matrix(~ j + k, data = d,
                  contrasts.arg = list(j = "contr.sum", k = "contr.sum"))
colnames(X) <- c("mu", "alpha1", "beta1", "beta2")
rownames(X) <- paste0("Y", d$j, d$k)
X[, ]
c(rank = qr(X)$rank, ncol = ncol(X))

    mu alpha1 beta1 beta2
Y11  1      1     1     0
Y12  1      1     0     1
Y13  1      1    -1    -1
Y21  1     -1     1     0
Y22  1     -1     0     1
Y23  1     -1    -1    -1

rank ncol
   4    4

\(\mathbf{X}\) has full column rank \(4\), so the four parameters are estimable — which the unconstrained six-parameter version was not. With \(N = 6\) observations and \(p = 4\) parameters only \(2\) degrees of freedom remain for estimating \(\sigma^2\), the price of having no replication; and because there is no replication no interaction term can be included, since a \(j \times k\) interaction would use up exactly those remaining \(2\) degrees of freedom and leave nothing to estimate the error variance with.

Finally, the design is balanced and the two factors are orthogonal to each other and to the intercept:

t(X) %*% X

       mu alpha1 beta1 beta2
mu      6      0     0     0
alpha1  0      6     0     0
beta1   0      0     4     2
beta2   0      0     2     4

The block-diagonal structure of \(\mathbf{X}^{\mathrm{T}}\mathbf{X}\) means \(\widehat{\mu}\), \(\widehat{\alpha}_1\) and the pair \((\widehat{\beta}_1,\widehat{\beta}_2)\) are mutually uncorrelated, so the estimated effect of one factor does not change when the other factor is added to or dropped from the model. The off-diagonal \(2\) in the lower block is not a defect of the design: it reflects only that the two columns coding a three-level factor cannot both be orthogonal to each other under the sum-to-zero parameterization.

Each part is answered with the three components of a generalized linear model set out in Section 3.4: the response distribution from the exponential family, the linear component \(\mathbf{x}_i^T\boldsymbol\beta\), and the link \(g\) with \(g(\mu_i)=\mathbf{x}_i^T\boldsymbol\beta\).

Part (a). The response is \(Y_i\), the weight (kg) of person \(i\). Weight is a continuous measurement, roughly symmetric within a homogeneous group and with no hard upper bound, so \(Y_i \sim N(\mu_i,\sigma^2)\) is the natural first choice; the Normal is in the exponential family (Section 3.2.2) and gives the normal linear model of Section 3.5.1 with the identity link.

The explanatory variables are \(x_{i1}\) age (years), \(x_{i2}\) sex (a dummy, \(1\) for male and \(0\) for female), \(x_{i3}\) height (cm), \(x_{i4}\) mean daily food intake (kJ) and \(x_{i5}\) mean daily energy expenditure (kJ). The linear component is \[ \mathbf{x}_i^T\boldsymbol\beta = \beta_1 + \beta_2 x_{i1} + \beta_3 x_{i2} + \beta_4 x_{i3} + \beta_5 x_{i4} + \beta_6 x_{i5}, \qquad g(\mu_i)=\mu_i . \] Two remarks. Weight is strictly positive and right-skewed in the general population; if the fitted model produced negative predictions or funnel-shaped residuals one would move to \(\log(\text{weight})\) as the response, or to a Gamma distribution with a log link, which keeps \(\mu_i>0\). Also, intake and expenditure enter as a balance, so in practice one would consider the single covariate \(x_{i4}-x_{i5}\).

Part (b). Here the mice at a given exposure level form a set of \(n=20\) independent binary trials with a common infection probability, so the response is \(Y_i\), the number infected out of \(n_i = 20\) at exposure level \(i\) (\(i=1,\dots,5\)), and \(Y_i \sim \mathrm{Bin}(20,\pi_i)\). The Binomial is in the exponential family with natural parameter \(\log[\pi/(1-\pi)]\) (Section 3.2.3), and it is the model of first choice for a count of successes out of a known number of trials. The proportion \(Y_i/n_i\) is just a rescaling of \(Y_i\); modelling the count keeps the correct binomial variance \(n_i\pi_i(1-\pi_i)\), which the proportion alone would not convey.

The explanatory variable is \(x_i\), the exposure level (or \(\log\) dose, which is usual for dose-response work). With the natural link, the logit, \[ g(\pi_i)=\log\!\left(\frac{\pi_i}{1-\pi_i}\right) = \beta_1 + \beta_2 x_i . \] This is a logistic dose-response model; it keeps \(\pi_i \in (0,1)\) for every value of \(\beta_1,\beta_2\) and \(x_i\), unlike a linear model for \(\pi_i\) itself.

Part (c). The response is \(Y_i\), the number of shopping trips made by household \(i\) in a week. This is a count of events in a fixed time period with no natural upper bound, so \(Y_i \sim \mathrm{Po}(\mu_i)\) (Section 3.2.1).

The explanatory variables are \(x_{i1}\) the number of people in the household, \(x_{i2}\) the household income and \(x_{i3}\) the distance to the supermarket. The log link is used because it keeps \(\mu_i>0\) and makes the covariate effects multiplicative on the trip rate, \[ g(\mu_i)=\log \mu_i = \beta_1 + \beta_2 x_{i1} + \beta_3 x_{i2} + \beta_4 x_{i3}. \] One would expect \(\beta_2>0\) (more people, more trips) and \(\beta_4<0\) (further away, fewer trips). If the observed variance of the counts exceeds their mean the data are overdispersed and a negative binomial model (Chapter 9) would be preferred; the Poisson assumption \(\mathrm{E}(Y)=\mathrm{var}(Y)\) should be checked rather than assumed.

Take logarithms of the density and collect terms in \(y\) and in \(\beta\) separately, for \(y>0\): \[ f(y;\beta) = \exp\left[ \alpha\log\beta - \log\Gamma(\alpha) + (\alpha-1)\log y - y\beta \right]. \] This is exactly the canonical form (3.3), \(f(y;\theta)=\exp[a(y)b(\theta)+c(\theta)+d(y)]\), with \[ a(y)=y, \qquad b(\beta) = -\beta, \qquad c(\beta) = \alpha\log\beta - \log\Gamma(\alpha), \qquad d(y) = (\alpha-1)\log y . \] Since \(a(y)=y\) the distribution is in canonical form, and by the definition in Section 3.2 the natural parameter is \[ b(\beta) = -\beta . \] The known shape \(\alpha\) is a nuisance parameter and appears only inside the known functions \(c\) and \(d\), as Section 3.2 permits.

For the moments use (3.9) and (3.12). The required derivatives are \[ b^{\prime}(\beta) = -1, \quad b^{\prime\prime}(\beta) = 0, \quad c^{\prime}(\beta) = \frac{\alpha}{\beta}, \quad c^{\prime\prime}(\beta) = -\frac{\alpha}{\beta^2}. \] From (3.9), \(\mathrm{E}[a(Y)]=-c^{\prime}(\beta)/b^{\prime}(\beta)\), and since \(a(Y)=Y\), \[ \mathrm{E}(Y) = -\frac{\alpha/\beta}{-1} = \frac{\alpha}{\beta}. \] From (3.12), \(\mathrm{var}[a(Y)] = \left[ b^{\prime\prime}(\beta)c^{\prime}(\beta) - c^{\prime\prime}(\beta)b^{\prime}(\beta) \right] / [b^{\prime}(\beta)]^3\), so \[ \mathrm{var}(Y) = \frac{0\cdot(\alpha/\beta) - (-\alpha/\beta^2)(-1)}{(-1)^3} = \frac{-\alpha/\beta^2}{-1} = \frac{\alpha}{\beta^2}. \] These are the familiar Gamma moments. Note the mean-variance relationship \(\mathrm{var}(Y) = [\mathrm{E}(Y)]^2/\alpha\): the standard deviation is proportional to the mean, which is why Gamma models are used for positive, right-skewed responses whose scatter grows with their level. With \(\alpha=1\) this reduces to the exponential distribution of Exercise 3.3(b), where \(\mathrm{var}(Y)=[\mathrm{E}(Y)]^2\).

A numerical check that the algebra above matches R’s own parametrisation (R’s rate is \(\beta\)):

set.seed(1); a <- 3.5; b <- 2
y <- rgamma(2e6, shape = a, rate = b)
c(mean = mean(y), theory_mean = a/b, var = var(y), theory_var = a/b^2)

     mean theory_mean         var  theory_var
1.7500115   1.7500000   0.8750181   0.8750000

In every case the method is the same: take logarithms and try to write the result in the form (3.3), \[ f(y;\theta) = \exp\left[ a(y)b(\theta) + c(\theta) + d(y) \right], \] which requires that \(y\) and \(\theta\) appear together in exactly one product term. Recall from Section 3.2 that the distribution is in canonical form only if \(a(y)=y\).

Part (a): Pareto. Here \(y>1\) and \(\theta>0\). Then \[ f(y;\theta) = \theta y^{-\theta-1} = \exp\left[ \log\theta - (\theta+1)\log y \right] = \exp\left[ (\log y)(-\theta) + \log\theta - \log y \right], \] which is of the form (3.3) with \[ a(y) = \log y, \qquad b(\theta) = -\theta, \qquad c(\theta) = \log\theta, \qquad d(y) = -\log y . \] So the Pareto belongs to the exponential family, but it is not in canonical form, because \(a(y)=\log y \neq y\). Equivalently, \(\log Y\) rather than \(Y\) is the quantity with a canonical exponential-family distribution: from the display above, \(\log Y \sim \mathrm{Exp}(\theta)\). This observation is used again in Exercises 3.9 and 3.11.

Part (b): Exponential. For \(y>0\), \[ f(y;\theta) = \theta e^{-y\theta} = \exp\left[ -y\theta + \log\theta \right], \] so (3.3) holds with \[ a(y)=y, \qquad b(\theta)=-\theta, \qquad c(\theta)=\log\theta, \qquad d(y)=0 . \] Since \(a(y)=y\) this is in canonical form with natural parameter \(-\theta\). It is the \(\alpha=1\) special case of the Gamma density of Exercise 3.2, as confirmed in Exercise 3.12(a).

Part (c): Negative Binomial. Here \(y=0,1,2,\dots\) counts the failures before the \(r\)th success, and \(r\) is known, so the binomial coefficient is a function of \(y\) alone. Taking logarithms, \[ f(y;\theta) = \exp\left[ y\log(1-\theta) + r\log\theta + \log\binom{y+r-1}{r-1} \right], \] so (3.3) holds with \[ a(y)=y, \qquad b(\theta)=\log(1-\theta), \qquad c(\theta)=r\log\theta, \qquad d(y)=\log\binom{y+r-1}{r-1}. \] Again \(a(y)=y\), so the negative binomial is in canonical form with natural parameter \(\log(1-\theta)\). Note that \(r\) being known is essential: if \(r\) were an unknown parameter, then \(d(y)\) would depend on it and the density would no longer be a one-parameter exponential family member of the form (3.3).

The two results to be used are, from Section 3.3, \[ \mathrm{E}[a(Y)] = -\frac{c^{\prime}(\theta)}{b^{\prime}(\theta)} \tag{3.9} \] and \[ \mathrm{var}[a(Y)] = \frac{b^{\prime\prime}(\theta)c^{\prime}(\theta)-c^{\prime\prime}(\theta)b^{\prime}(\theta)}{[b^{\prime}(\theta)]^3}. \tag{3.12} \] All three distributions are in canonical form, so \(a(Y)=Y\) throughout and the left-hand sides are the mean and variance of \(Y\) itself. The functions \(b\), \(c\), \(d\) are read off Table 3.1.

Part (a): Poisson. From Section 3.2.1, \(b(\theta)=\log\theta\), \(c(\theta)=-\theta\), \(d(y)=-\log y!\), so \[ b^{\prime}(\theta)=\frac{1}{\theta}, \quad b^{\prime\prime}(\theta)=-\frac{1}{\theta^2}, \quad c^{\prime}(\theta)=-1, \quad c^{\prime\prime}(\theta)=0. \] By (3.9), \[ \mathrm{E}(Y) = -\frac{-1}{1/\theta} = \theta, \] and by (3.12), \[ \mathrm{var}(Y) = \frac{(-1/\theta^2)(-1) - 0\cdot(1/\theta)}{(1/\theta)^3} = \frac{1/\theta^2}{1/\theta^3} = \theta. \]

Part (b): Normal. Here \(\mu\) is the parameter of interest and \(\sigma^2\) is a known nuisance parameter. From Section 3.2.2, \(b(\mu)=\mu/\sigma^2\) and \(c(\mu) = -\mu^2/(2\sigma^2) - \tfrac12\log(2\pi\sigma^2)\), so \[ b^{\prime}(\mu)=\frac{1}{\sigma^2}, \quad b^{\prime\prime}(\mu)=0, \quad c^{\prime}(\mu)=-\frac{\mu}{\sigma^2}, \quad c^{\prime\prime}(\mu)=-\frac{1}{\sigma^2}. \] By (3.9), \[ \mathrm{E}(Y) = -\frac{-\mu/\sigma^2}{1/\sigma^2} = \mu, \] and by (3.12), \[ \mathrm{var}(Y) = \frac{0\cdot(-\mu/\sigma^2) - (-1/\sigma^2)(1/\sigma^2)}{(1/\sigma^2)^3} = \frac{1/\sigma^4}{1/\sigma^6} = \sigma^2. \]

Part (c): Binomial. From Section 3.2.3, \(b(\pi)=\log[\pi/(1-\pi)] = \log\pi - \log(1-\pi)\) and \(c(\pi)=n\log(1-\pi)\). Differentiating, \[ b^{\prime}(\pi) = \frac{1}{\pi} + \frac{1}{1-\pi} = \frac{1}{\pi(1-\pi)}, \qquad c^{\prime}(\pi) = -\frac{n}{1-\pi}. \] By (3.9), \[ \mathrm{E}(Y) = -\frac{-n/(1-\pi)}{1/[\pi(1-\pi)]} = \frac{n}{1-\pi}\,\pi(1-\pi) = n\pi. \] For the variance, differentiate again. Writing \(b^{\prime}(\pi)=[\pi(1-\pi)]^{-1}\), \[ b^{\prime\prime}(\pi) = -\frac{1-2\pi}{\pi^2(1-\pi)^2}, \qquad c^{\prime\prime}(\pi) = -\frac{n}{(1-\pi)^2}. \] Substituting into (3.12), the numerator is \[ \left(-\frac{1-2\pi}{\pi^2(1-\pi)^2}\right)\!\left(-\frac{n}{1-\pi}\right) - \left(-\frac{n}{(1-\pi)^2}\right)\!\left(\frac{1}{\pi(1-\pi)}\right) = \frac{n(1-2\pi)}{\pi^2(1-\pi)^3} + \frac{n}{\pi(1-\pi)^3}, \] and the denominator is \([b^{\prime}(\pi)]^3 = [\pi(1-\pi)]^{-3}\). Hence \[ \mathrm{var}(Y) = \pi^3(1-\pi)^3\left[\frac{n(1-2\pi)}{\pi^2(1-\pi)^3} + \frac{n}{\pi(1-\pi)^3}\right] = n\pi(1-2\pi) + n\pi^2 = n\pi(1-\pi). \]

As an arithmetic check, (3.9) and (3.12) are applied numerically to the same \(b\) and \(c\) functions, using central differences for the derivatives, and compared with the known moments.

ef <- function(b, cc, theta, h = 1e-5) {
  d1 <- function(f, t) (f(t + h) - f(t - h)) / (2 * h)
  d2 <- function(f, t) (f(t + h) - 2 * f(t) + f(t - h)) / h^2
  bp <- d1(b, theta); bpp <- d2(b, theta)
  cp <- d1(cc, theta); cpp <- d2(cc, theta)
  c(mean = -cp / bp, var = (bpp * cp - cpp * bp) / bp^3)
}
th <- 4.3
poi <- ef(function(t) log(t), function(t) -t, th)
sig <- 1.7; mu <- 2.4
nor <- ef(function(t) t / sig^2,
          function(t) -t^2 / (2 * sig^2) - 0.5 * log(2 * pi * sig^2), mu)
n <- 12; pi0 <- 0.3
bin <- ef(function(t) log(t / (1 - t)), function(t) n * log(1 - t), pi0)
round(rbind(Poisson  = c(poi, truth.mean = th,      truth.var = th),
            Normal   = c(nor, truth.mean = mu,      truth.var = sig^2),
            Binomial = c(bin, truth.mean = n * pi0, truth.var = n * pi0 * (1 - pi0))), 6)

         mean      var truth.mean truth.var
Poisson   4.3 4.300009        4.3      4.30
Normal    2.4 2.889965        2.4      2.89
Binomial  3.6 2.520000        3.6      2.52

The tiny discrepancy in the Poisson and Normal variances is the truncation error of the numerical second derivative, not a failure of the identities.

Part (a). From Exercise 3.3(c) the negative binomial is in canonical form with \[ a(y)=y, \qquad b(\theta)=\log(1-\theta), \qquad c(\theta)=r\log\theta, \] so \[ b^{\prime}(\theta) = -\frac{1}{1-\theta}, \quad b^{\prime\prime}(\theta) = -\frac{1}{(1-\theta)^2}, \quad c^{\prime}(\theta) = \frac{r}{\theta}, \quad c^{\prime\prime}(\theta) = -\frac{r}{\theta^2}. \] By (3.9), \[ \mathrm{E}(Y) = -\frac{r/\theta}{-1/(1-\theta)} = \frac{r(1-\theta)}{\theta}. \] By (3.12) the numerator is \[ b^{\prime\prime}c^{\prime} - c^{\prime\prime}b^{\prime} = \left(-\frac{1}{(1-\theta)^2}\right)\frac{r}{\theta} - \left(-\frac{r}{\theta^2}\right)\left(-\frac{1}{1-\theta}\right) = -\frac{r}{\theta(1-\theta)^2} - \frac{r}{\theta^2(1-\theta)}, \] and the denominator is \([b^{\prime}(\theta)]^3 = -1/(1-\theta)^3\). Dividing, \[ \mathrm{var}(Y) = (1-\theta)^3\left[\frac{r}{\theta(1-\theta)^2} + \frac{r}{\theta^2(1-\theta)}\right] = \frac{r(1-\theta)}{\theta} + \frac{r(1-\theta)^2}{\theta^2} = \frac{r(1-\theta)}{\theta^2}. \] So \[ \mathrm{E}(Y) = \frac{r(1-\theta)}{\theta}, \qquad \mathrm{var}(Y) = \frac{r(1-\theta)}{\theta^2} = \frac{\mathrm{E}(Y)}{\theta}. \] Since \(0<\theta<1\) we have \(\mathrm{var}(Y) > \mathrm{E}(Y)\), with the ratio \(1/\theta\) measuring how far the distribution is from Poisson.

r <- 5; th <- 0.35
h <- 1e-5
b  <- function(t) log(1 - t); cc <- function(t) r * log(t)
d1 <- function(f, t) (f(t + h) - f(t - h)) / (2 * h)
d2 <- function(f, t) (f(t + h) - 2 * f(t) + f(t - h)) / h^2
bp <- d1(b, th); bpp <- d2(b, th); cp <- d1(cc, th); cpp <- d2(cc, th)
set.seed(2); y <- rnbinom(2e6, size = r, prob = th)
round(c(mean.3.9 = -cp / bp, formula.mean = r * (1 - th) / th, sim.mean = mean(y),
        var.3.12 = (bpp * cp - cpp * bp) / bp^3, formula.var = r * (1 - th) / th^2,
        sim.var = var(y)), 4)

mean.3.9 formula.mean     sim.mean     var.3.12  formula.var      sim.var
  9.2857       9.2857       9.2908      26.5306      26.5306      26.5283

Part (b). The three distributions are all models for counts, and they differ in exactly one respect that matters for practice: the mean-variance relationship they impose. Writing \(\mu = \mathrm{E}(Y)\), \[ \text{Binomial: } \mathrm{var}(Y) = \mu(1-\pi) < \mu, \qquad \text{Poisson: } \mathrm{var}(Y) = \mu, \qquad \text{Negative Binomial: } \mathrm{var}(Y) = \mu/\theta > \mu. \] Because a generalized linear model has no free variance parameter once the distribution is chosen (there is no \(\sigma^2\) to absorb misspecification), the choice of distribution is a choice of variance function. This has several practical consequences.

First, the choice should follow the sampling mechanism where one exists. If the count has a known upper limit \(n\) arising from \(n\) independent binary trials, use the Binomial; underdispersion relative to Poisson is then a feature of the design, not something to correct.

Second, for unbounded counts the Poisson is the default, but its equality \(\mathrm{E}(Y)=\mathrm{var}(Y)\) is a strong and testable assumption. It should be checked, for example by comparing the residual deviance (or Pearson \(X^2\)) with its residual degrees of freedom, since both have approximate expectation equal to the degrees of freedom when the model is right.

Third, real count data are frequently overdispersed, as noted in Section 3.2.1: unmeasured heterogeneity between units, clustering, or an excess of zeros all inflate the variance above the mean. The negative binomial is the natural remedy, since it has the same log link and the same interpretation of \(\boldsymbol\beta\) as multiplicative rate effects, but a variance \(\mu/\theta\) that exceeds the mean; the extra parameter is estimated from the data. Chapter 9 develops these models.

Finally, the practical cost of ignoring overdispersion is not in the point estimates, which usually remain reasonable, but in the standard errors: fitting a Poisson model to overdispersed data understates \(\mathrm{var}(\hat{\boldsymbol\beta})\) and so produces confidence intervals that are too narrow and \(p\)-values that are too small. Underdispersion has the reverse effect.

The key feature of the model of Example 3.5.3 is that \(\log n_i\) enters the linear predictor with its coefficient fixed at one (that is, as an offset) and that there is no intercept: the only free parameter is \(\theta\). Forcing the fitted line through the origin at \(i=0\) is a strong constraint, and it is the reason the model fails here.

library(dobson)
data(mortality)
d <- as.data.frame(mortality)
names(d) <- c("agegroup", "deaths", "population")
d$i <- 1:8
d$rate <- d$deaths / d$population * 1e5
d

  agegroup deaths population i       rate
1    30-34      1      17742 1   5.636343
2    35-39      5      16554 2  30.204180
3    40-44      5      16059 3  31.135189
4    45-49     12      13083 4  91.722082
5    50-54     25      10784 5 231.824926
6    55-59     38       9645 6 393.986522
7    60-64     54      10706 7 504.390062
8    65-69     65       9933 8 654.384375

Note that the poisson family object is masked by a data set of the same name in the dobson package, so it is referred to explicitly as stats::poisson below.

m1 <- glm(deaths ~ i - 1 + offset(log(population)), family = stats::poisson, data = d)
summary(m1)

Call:
glm(formula = deaths ~ i - 1 + offset(log(population)), family = stats::poisson,
    data = d)

Coefficients:
  Estimate Std. Error z value Pr(>|z|)
i -2.72150    0.02583  -105.4   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 206303.6  on 8  degrees of freedom
Residual deviance:   6990.4  on 7  degrees of freedom
AIC: 7028.1

Number of Fisher Scoring iterations: 8

The model is not adequate. The residual deviance is \(6990.4\) on \(7\) degrees of freedom; if the model were correct this would be an observation from approximately \(\chi^2(7)\), whose mean is \(7\). The estimate \(\hat\theta = -2.72\) is negative, so the fitted death rate \(e^{\hat\theta i}\) decreases with age, contradicting Figure 3.2 entirely. The reason is visible in the table of expected counts below: with no intercept the fitted rate at \(i=1\) is forced to be \(e^{\theta}\), and no single \(\theta\) can make \(e^{\theta i}\) pass through rates that run from \(5.6\) to \(654\) per \(100{,}000\) (that is, from \(5.6\times10^{-5}\) to \(6.5\times10^{-3}\)) while also starting near \(e^{\theta}\). The model is badly misspecified because it lacks an intercept; it is not that the exponential shape is wrong.

Following the hint, the log rate is regressed on the age-group index by ordinary least squares. Figure 3.2 plots \(\log(y_i/n_i)\) against \(i\) and is approximately linear, so the transformation is a log of the response rate with the explanatory variable left untransformed (an equally spaced age-group index is equivalent to using the group midpoints up to a linear rescaling).

d$lograte <- log(d$deaths / d$population)
ols <- lm(lograte ~ i, data = d)
summary(ols)

Call:
lm(formula = lograte ~ i, data = d)

Residuals:
     Min       1Q   Median       3Q      Max
-0.59459 -0.28691  0.05263  0.34854  0.46029

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.85457    0.34698 -28.401 1.26e-07 ***
i            0.66547    0.06871   9.685 6.95e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4453 on 6 degrees of freedom
Multiple R-squared:  0.9399,	Adjusted R-squared:  0.9299
F-statistic:  93.8 on 1 and 6 DF,  p-value: 6.95e-05

The straight line explains \(94\%\) of the variation in the log rate, and \(\hat\beta_2 = 0.665\) means the death rate multiplies by \(e^{0.665}=1.95\) for each five-year increase in age, that is, it roughly doubles every five years. This is the familiar Gompertz-type law of mortality.

Simple linear regression is nevertheless the wrong tool for these data: it assumes constant variance on the log-rate scale, whereas the sampling variability of \(\log(y_i/n_i)\) is approximately \(1/y_i\), which is \(1\) for the first group and \(0.015\) for the last, a 65-fold range. It also gives the group with one death the same weight as the group with 65. The correct version of the same model is the Poisson generalized linear model with the offset retained and an intercept added, \[ \log \mu_i = \log n_i + \beta_1 + \beta_2 i . \]

m2 <- glm(deaths ~ i + offset(log(population)), family = stats::poisson, data = d)
summary(m2)

Call:
glm(formula = deaths ~ i + offset(log(population)), family = stats::poisson,
    data = d)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -9.01688    0.26188  -34.43   <2e-16 ***
i            0.52219    0.03905   13.37   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 257.51  on 7  degrees of freedom
Residual deviance:  14.69  on 6  degrees of freedom
AIC: 54.376

Number of Fisher Scoring iterations: 4

The Poisson slope \(\hat\beta_2 = 0.522\) (standard error \(0.039\)) is appreciably smaller than the OLS slope \(0.665\), precisely because the weighting has changed: the poorly determined early groups no longer dominate. The rate ratio per five-year age band is \(e^{0.522}=1.69\) with 95\% confidence interval \(\exp(0.522 \pm 1.96\times 0.039) = (1.56, 1.82)\).

Following the hint, observed and expected deaths are compared group by group.

m3 <- glm(deaths ~ i + I(i^2) + offset(log(population)), family = stats::poisson, data = d)
data.frame(agegroup = d$agegroup, observed = d$deaths,
           exp.Ex353 = round(fitted(m1), 2),
           exp.OLS   = round(d$population * exp(predict(ols)), 2),
           exp.m2    = round(fitted(m2), 2),
           exp.m3    = round(fitted(m3), 2),
           pearson.m2 = round(residuals(m2, type = "pearson"), 2))

  agegroup observed exp.Ex353 exp.OLS exp.m2 exp.m3 pearson.m2
1    30-34        1   1167.00    1.81   3.63   1.01      -1.38
2    35-39        5     71.62    3.29   5.71   2.94      -0.30
3    40-44        5      4.57    6.21   9.33   7.61      -1.42
4    45-49       12      0.24    9.84  12.82  14.23      -0.23
5    50-54       25      0.01   15.78  17.81  23.09       1.70
6    55-59       38      0.00   27.45  26.86  34.91       2.15
7    60-64       54      0.00   59.28  50.25  56.23       0.53
8    65-69       65      0.00  107.00  78.59  64.99      -1.53

The first column of expectations shows how absurd the Example 3.5.3 model is: it predicts 1167 deaths in the youngest group, where one occurred, and essentially zero in the oldest, where 65 occurred. The OLS-based expectations are much better but still poor at the extremes (107 expected against 65 observed in the oldest group), and they do not even reproduce the total: they sum to \(230.7\) against the \(205\) deaths actually observed, whereas the Poisson fit with an intercept reproduces the total exactly, because the score equation for the intercept forces \(\sum \hat\mu_i = \sum y_i\).

The Poisson model with an intercept fits acceptably but not perfectly: deviance \(14.69\) on \(6\) degrees of freedom, \(p = 0.023\). The Pearson residuals show the pattern responsible, negative at both ends and positive in the middle, which is exactly the signature of curvature in \(\log(\text{rate})\) against \(i\). Adding a quadratic term removes it.

anova(m2, m3, test = "Chisq")

Analysis of Deviance Table

Model 1: deaths ~ i + offset(log(population))
Model 2: deaths ~ i + I(i^2) + offset(log(population))
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)
1         6    14.6899
2         5     3.0938  1   11.596 0.0006609 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

gof <- function(m) c(deviance = deviance(m), df = df.residual(m),
                     p = pchisq(deviance(m), df.residual(m), lower.tail = FALSE),
                     AIC = AIC(m))
round(rbind(Ex353 = gof(m1), linear = gof(m2), quadratic = gof(m3)), 4)

           deviance df      p       AIC
Ex353     6990.4444  7 0.0000 7028.1303
linear      14.6899  6 0.0228   54.3758
quadratic    3.0938  5 0.6855   44.7797

par(mfrow = c(1, 2))
plot(d$i, log(d$rate), pch = 16, xaxt = "n",
     xlab = "Age group", ylab = "log(death rate per 100,000)",
     main = "Log rate vs age group")
axis(1, at = d$i, labels = d$agegroup, cex.axis = 0.7)
g <- seq(1, 8, length = 100)
lines(g, coef(m1)[1] * g + log(1e5), lty = 3, lwd = 2)
lines(g, coef(m2)[1] + coef(m2)[2] * g + log(1e5), lty = 1, lwd = 2)
lines(g, coef(m3)[1] + coef(m3)[2] * g + coef(m3)[3] * g^2 + log(1e5), lty = 2, lwd = 2)
legend("topleft", c("Example 3.5.3 (no intercept)", "linear in i", "quadratic in i"),
       lty = c(3, 1, 2), lwd = 2, bty = "n", cex = 0.8)
plot(d$i, residuals(m2, type = "pearson"), pch = 16, xaxt = "n", ylim = c(-3, 3),
     xlab = "Age group", ylab = "Pearson residual", main = "Residuals, linear model")
axis(1, at = d$i, labels = d$agegroup, cex.axis = 0.7)
abline(h = 0, lty = 2); points(d$i, residuals(m3, type = "pearson"), pch = 1)
legend("topleft", c("linear in i", "quadratic in i"), pch = c(16, 1), bty = "n", cex = 0.8)

Conclusions. The model of Example 3.5.3 as literally written, with no intercept, is not adequate: it is off by three orders of magnitude and gives a slope of the wrong sign. Its exponential shape is right, and once an intercept is admitted the model \(\log\mu_i = \log n_i + \beta_1 + \beta_2 i\) describes the data reasonably, with the coronary death rate rising by a factor of about \(1.69\) per five-year age band. Even so, the deviance of \(14.69\) on \(6\) df indicates mild lack of fit, and the residual pattern shows the log rate is slightly concave in age; the quadratic model has deviance \(3.09\) on \(5\) df (\(p=0.69\)) and an AIC lower by nearly \(10\), so it is the preferred description. Simple least squares on \(\log(y_i/n_i)\) recovers the same qualitative story but exaggerates the slope and mismatches the observed totals, because it ignores the Poisson mean-variance relationship and the very different amounts of information in the eight groups.

Part (a). Take logarithms of the probability function and collect terms: \[ f(y_i;\pi_i) = \pi_i^{y_i}(1-\pi_i)^{1-y_i} = \exp\left[ y_i\log\pi_i + (1-y_i)\log(1-\pi_i) \right] = \exp\left[ y_i\log\left(\frac{\pi_i}{1-\pi_i}\right) + \log(1-\pi_i) \right]. \] This is of the form (3.3), \(f(y;\theta)=\exp[a(y)b(\theta)+c(\theta)+d(y)]\), with \[ a(y_i)=y_i, \qquad b(\pi_i)=\log\left(\frac{\pi_i}{1-\pi_i}\right), \qquad c(\pi_i)=\log(1-\pi_i), \qquad d(y_i)=0, \] so the Bernoulli distribution belongs to the exponential family. This is the \(n=1\) case of the Binomial calculation in Section 3.2.3, since \(\binom{1}{y}=1\) and \(n\log(1-\pi)=\log(1-\pi)\).

Part (b). Because \(a(y_i)=y_i\), the distribution is in canonical form, and by the definition in Section 3.2 the natural parameter is \(b(\pi_i)\), that is, \[ b(\pi_i) = \log\left(\frac{\pi_i}{1-\pi_i}\right), \] the logit. This is why the logit is the natural (canonical) link for binary data: with \(g\) equal to the logit, the linear predictor \(\mathbf{x}_i^T\boldsymbol\beta\) is the natural parameter itself.

Part (c). Directly, \(\mathrm{E}(Y_i)=1\cdot\pi_i + 0\cdot(1-\pi_i)=\pi_i\). To use the results of the chapter instead, apply (3.9) with \(a(Y)=Y\): \[ b^{\prime}(\pi_i) = \frac{d}{d\pi_i}\left[\log\pi_i-\log(1-\pi_i)\right] = \frac{1}{\pi_i(1-\pi_i)}, \qquad c^{\prime}(\pi_i) = -\frac{1}{1-\pi_i}, \] so \[ \mathrm{E}(Y_i) = -\frac{c^{\prime}(\pi_i)}{b^{\prime}(\pi_i)} = \frac{1/(1-\pi_i)}{1/[\pi_i(1-\pi_i)]} = \pi_i . \] This is Exercise 3.4(c) with \(n=1\).

Part (d). Exponentiate both sides of \(\log[\pi/(1-\pi)] = \mathbf{x}^T\boldsymbol\beta\): \[ \frac{\pi}{1-\pi} = e^{\mathbf{x}^T\boldsymbol\beta} \implies \pi = (1-\pi)e^{\mathbf{x}^T\boldsymbol\beta} \implies \pi\left(1+e^{\mathbf{x}^T\boldsymbol\beta}\right) = e^{\mathbf{x}^T\boldsymbol\beta}, \] hence \[ \pi = \frac{e^{\mathbf{x}^T\boldsymbol\beta}}{1+e^{\mathbf{x}^T\boldsymbol\beta}} = \frac{1}{1+e^{-\mathbf{x}^T\boldsymbol\beta}} . \] The two statements are equivalent because the logit is a strictly increasing bijection from \((0,1)\) onto \(\mathbb{R}\), with the expression above as its inverse. The practical point is that whatever real value \(\mathbf{x}^T\boldsymbol\beta\) takes, the fitted \(\pi\) automatically lies strictly between \(0\) and \(1\).

Part (e). The graph of \[ \pi(x) = \frac{e^{\beta_1+\beta_2 x}}{1+e^{\beta_1+\beta_2 x}} \] is the sigmoid (“S-shaped”) logistic curve. Its properties, taking \(\beta_2>0\), are as follows. It is strictly increasing, since \[ \frac{d\pi}{dx} = \beta_2\,\pi(1-\pi) > 0 ; \] it has horizontal asymptotes \(\pi \to 0\) as \(x\to-\infty\) and \(\pi \to 1\) as \(x\to+\infty\); it passes through \(\pi = 1/2\) at \(x = -\beta_1/\beta_2\); and it is symmetric about that point, which is also the point of steepest ascent, where the slope equals \(\beta_2/4\). If \(\beta_2<0\) the curve is decreasing, and \(\beta_2=0\) gives a horizontal line at \(\pi = e^{\beta_1}/(1+e^{\beta_1})\). Larger \(|\beta_2|\) makes the transition from \(0\) to \(1\) sharper.

x <- seq(-2, 12, length = 400)
logistic <- function(x, b1, b2) exp(b1 + b2 * x) / (1 + exp(b1 + b2 * x))
par(mfrow = c(1, 2))
plot(x, logistic(x, -5, 1), type = "l", lwd = 2, ylim = c(0, 1),
     xlab = "dose x", ylab = expression(pi),
     main = expression(paste("Logistic curve, ", beta[1] == -5, ", ", beta[2] == 1)))
abline(h = c(0, 0.5, 1), lty = 3); abline(v = 5, lty = 2)
text(5, 0.05, "LD50 = -b1/b2", pos = 4, cex = 0.9)
plot(x, logistic(x, -5, 1), type = "l", lwd = 2, ylim = c(0, 1),
     xlab = "dose x", ylab = expression(pi), main = "Effect of the slope")
lines(x, logistic(x, -2.5, 0.5), lwd = 2, lty = 2)
lines(x, logistic(x, -10, 2), lwd = 2, lty = 3)
legend("topleft", c("b2 = 0.5", "b2 = 1", "b2 = 2"), lty = c(2, 1, 3), lwd = 2, bty = "n")

Interpretation as a dose-response curve. If \(x\) is the dose of an insecticide and \(\pi\) the probability that an insect dies, the curve is the classical dose-response relationship.

At very low doses almost no insects die and \(\pi \approx 0\); at very high doses almost all die and \(\pi\approx 1\). Neither extreme is reached exactly, which is realistic: no dose is perfectly safe and no dose is guaranteed lethal.

Over the middle range the response rises steeply, so within a fairly narrow band of doses the kill rate goes from low to high. Outside that band, raising the dose further buys very little extra mortality, which is why the curve flattens.

The dose \(x = -\beta_1/\beta_2\) at which \(\pi = 1/2\) is the median lethal dose, LD50, the standard summary of potency; a more potent insecticide has a smaller LD50. In the left panel, with \(\beta_1=-5\) and \(\beta_2=1\), the LD50 is \(5\) and the maximum slope is \(\beta_2/4 = 0.25\) per unit dose.

The slope \(\beta_2\) has an exact odds interpretation: increasing the dose by one unit adds \(\beta_2\) to the log odds of dying, so it multiplies the odds of death by \(e^{\beta_2}\). This interpretation is constant across the dose range, whereas the effect on the probability is largest near the LD50 and negligible in the tails.

A natural underlying mechanism is a tolerance distribution: suppose insect \(j\) dies if the dose exceeds its individual tolerance \(T_j\), and that \(T_j\) has a logistic distribution across insects. Then \(\pi(x)=P(T_j \le x)\) is exactly the curve above, and the sigmoid shape simply reflects the accumulation of a unimodal spread of tolerances in the population.

Yes. Take logarithms and expand the two occurrences of \((y-\theta)/\phi\), remembering that \(\phi\) is treated as known: \[ \log f(y;\theta) = -\log\phi + \frac{y-\theta}{\phi} - \exp\left(\frac{y-\theta}{\phi}\right) = -\log\phi + \frac{y}{\phi} - \frac{\theta}{\phi} - e^{y/\phi}e^{-\theta/\phi}. \] The last term is the only one in which \(y\) and \(\theta\) appear together, and it already factorises as a function of \(y\) times a function of \(\theta\). So \[ f(y;\theta) = \exp\left[ a(y)b(\theta) + c(\theta) + d(y) \right] \] with \[ a(y) = e^{y/\phi}, \qquad b(\theta) = -e^{-\theta/\phi}, \qquad c(\theta) = -\frac{\theta}{\phi} - \log\phi, \qquad d(y) = \frac{y}{\phi}. \] This is form (3.3), so the Gumbel distribution is a member of the exponential family. It is not in canonical form, because \(a(y) = e^{y/\phi} \neq y\). Equivalently, it is \(e^{Y/\phi}\), not \(Y\) itself, that has the canonical exponential-family structure; indeed the calculation above shows \(e^{(Y-\theta)/\phi} \sim \mathrm{Exp}(1)\).

Two conditions were used and are worth stating explicitly. The support of \(f\) is the whole real line and does not depend on \(\theta\), as required for the differentiation-under-the-integral arguments of Section 3.3. And \(\phi\) must genuinely be known: if \(\phi\) were also unknown, \(d(y)=y/\phi\) and \(a(y)=e^{y/\phi}\) would depend on it and the family would be a two-parameter one, outside the one-parameter form (3.3).

Since the distribution is a member of the family, (3.9) and (3.12) apply to \(a(Y)=e^{Y/\phi}\). The derivatives are \[ b^{\prime}(\theta) = \frac{1}{\phi}e^{-\theta/\phi}, \quad b^{\prime\prime}(\theta) = -\frac{1}{\phi^2}e^{-\theta/\phi}, \quad c^{\prime}(\theta) = -\frac{1}{\phi}, \quad c^{\prime\prime}(\theta)=0, \] so from (3.9) \[ \mathrm{E}\left(e^{Y/\phi}\right) = -\frac{c^{\prime}(\theta)}{b^{\prime}(\theta)} = \frac{1/\phi}{e^{-\theta/\phi}/\phi} = e^{\theta/\phi}, \] and from (3.12) \[ \mathrm{var}\left(e^{Y/\phi}\right) = \frac{b^{\prime\prime}c^{\prime} - c^{\prime\prime}b^{\prime}}{[b^{\prime}]^3} = \frac{\left(-\phi^{-2}e^{-\theta/\phi}\right)\left(-\phi^{-1}\right)}{\phi^{-3}e^{-3\theta/\phi}} = e^{2\theta/\phi}. \] Note that these are moments of \(e^{Y/\phi}\), not of \(Y\); the exponential-family machinery gives the mean and variance of \(a(Y)\), and here \(a\) is not the identity. (For the record, \(\mathrm{E}(Y) = \theta - \gamma\phi\) with \(\gamma\) Euler’s constant, which the chapter’s results do not deliver directly.)

Both results are confirmed by simulation, using \(e^{(Y-\theta)/\phi}\sim\mathrm{Exp}(1)\) to generate the variates.

set.seed(7); theta <- 1.3; phi <- 0.7
u <- rexp(2e6)
y <- theta + phi * log(u)
a <- exp(y / phi)
round(c(mean.a = mean(a), theory.mean = exp(theta / phi),
        var.a = var(a),  theory.var = exp(2 * theta / phi)), 4)

mean.a theory.mean       var.a  theory.var
6.4100      6.4054     41.0751     41.0293

Check the three components of a generalized linear model listed at the end of Section 3.4.

Component 1: the response distribution. The \(Y_i\) are independent and, by Exercise 3.3(a), each Pareto density \[ f(y;\theta_i) = \theta_i y^{-\theta_i-1}, \qquad y>1, \] belongs to the exponential family, with \(a(y)=\log y\), \(b(\theta_i)=-\theta_i\), \(c(\theta_i)=\log\theta_i\) and \(d(y)=-\log y\). All the \(Y_i\) have the same form of distribution and each depends on a single parameter \(\theta_i\), as required.

There is one caveat, and it should be stated plainly. Dobson’s definition in Section 3.4 asks that each \(Y_i\) have the canonical form \(f(y_i;\theta_i)=\exp[y_ib(\theta_i)+c(\theta_i)+d(y_i)]\), and the Pareto does not: its \(a(y)=\log y\). Taken literally, the model therefore falls outside the definition as written. In practice this is a mild objection, because \(a\) is a fixed known transformation and \(\log Y_i \sim \mathrm{Exp}(\theta_i)\) is canonical, so the same model can be written as a generalized linear model for the transformed response. Most treatments, and the intent of the exercise, accept the Pareto here.

Component 2: the linear component. The systematic part involves \(\beta_0 + \beta_1 x_i\), which is linear in the parameters, with \(\mathbf{x}_i^T = [1\;\;x_i]\) and \(\boldsymbol\beta^T = [\beta_0\;\;\beta_1]\). Note that it is linearity in \(\boldsymbol\beta\) that matters, not linearity in \(x_i\).

Component 3: the link function. We need a monotone differentiable \(g\) with \(g(\mu_i)=\mathbf{x}_i^T\boldsymbol\beta\). Since \(\mu_i = (\beta_0+\beta_1x_i)^2\), take \[ g(\mu) = \sqrt{\mu} = \mu^{1/2}. \] On \(\mu>0\) this is strictly increasing and differentiable, so it is a valid link. (In modern terminology it is the power link with exponent \(1/2\).)

Conclusion. Yes, this is a generalized linear model, with response from the exponential family, linear predictor \(\beta_0+\beta_1x_i\), and link \(g(\mu)=\mu^{1/2}\) – subject to the canonical-form caveat above and to two restrictions on the parameter space that are worth making explicit.

First, \(\sqrt{\cdot}\) inverts \(\mu \mapsto \mu^2\) only on one branch, so \((\beta_0,\beta_1)\) and \((-\beta_0,-\beta_1)\) give identical means. Identifiability requires restricting to \(\beta_0+\beta_1x_i>0\) for all \(i\) in the data.

Second, integrating the Pareto density gives \(\mathrm{E}(Y)=\theta/(\theta-1)\) for \(\theta>1\), which is greater than \(1\) for every admissible \(\theta\) (and the mean does not exist for \(\theta \le 1\)). So the model is coherent only if \((\beta_0+\beta_1x_i)^2>1\), that is \(\beta_0+\beta_1x_i>1\), at every observed \(x_i\). This is a genuine restriction on the parameter space, not a defect of the GLM structure; it is the same kind of restriction that a Binomial model with an identity link would impose.

No, this is not a generalized linear model. It is a nonlinear regression model.

The response distribution is not the problem: the \(Y_i\) are independent and Normal, which is in the exponential family in canonical form (Section 3.2.2). The failure is in the systematic part, which cannot be written as \(g(\mu_i)=\mathbf{x}_i^T\boldsymbol\beta\) for any link \(g\).

To see this, note what such a representation would require. A link function is a fixed known function of \(\mu_i\) alone; it may not involve unknown parameters. Rearranging the model gives \[ \mu_i - \beta_0 = \log(\beta_1+\beta_2x_i), \qquad\text{that is}\qquad e^{\mu_i - \beta_0} = \beta_1 + \beta_2 x_i . \] The right-hand side is linear in \((\beta_1,\beta_2)\), so if \(\beta_0\) were known the model would be a generalized linear model, with link \(g(\mu)=e^{\mu-\beta_0}\) and linear predictor \(\beta_1+\beta_2x_i\). But \(\beta_0\) is unknown and must be estimated, so the candidate link \(g(\mu)=e^{\mu-\beta_0}\) depends on an unknown parameter and is therefore not a link function in the sense of Section 3.4.

Equivalently: in a generalized linear model, all the parameters enter through the single linear combination \(\mathbf{x}_i^T\boldsymbol\beta\). Here \(\beta_0\) enters outside the logarithm while \(\beta_1\) and \(\beta_2\) enter inside it, and there is no transformation \(g\) that puts all three into one linear form. The parameters are not separable in the way the definition requires.

Note also that the three parameters are not all identified without a constraint even in the nonlinear model, because \(\beta_0 + \log(\beta_1+\beta_2x)\) is unchanged if one adds \(\log k\) to \(\beta_0\) and divides both \(\beta_1\) and \(\beta_2\) by \(k\) for any \(k>0\). One would fix this by setting \(\beta_0=0\) or \(\beta_1=1\); with \(\beta_0=0\) the model does become a generalized linear model with link \(g(\mu)=e^{\mu}\).

As written, therefore, the model must be fitted by general nonlinear least squares (for instance R’s nls), not by the iterative weighted least squares of Chapter 4. The model also requires \(\beta_1+\beta_2x_i>0\) at every observation for the logarithm to be defined.

From Exercise 3.3(a) the Pareto density \(f(y;\theta)=\theta y^{-\theta-1}\), \(y>1\), is in the exponential family with \[ a(y)=\log y, \qquad b(\theta)=-\theta, \qquad c(\theta)=\log\theta, \qquad d(y)=-\log y, \] so the log-likelihood from a single observation is \[ l(\theta;y) = a(y)b(\theta) + c(\theta) + d(y) = -\theta\log y + \log\theta - \log y . \]

The score statistic. Differentiating with respect to \(\theta\), or equivalently using (3.13), \(U = a(Y)b^{\prime}(\theta)+c^{\prime}(\theta)\) with \(b^{\prime}(\theta)=-1\) and \(c^{\prime}(\theta)=1/\theta\), \[ U = \frac{dl}{d\theta} = \frac{1}{\theta} - \log Y . \]

Verification that \(\mathrm{E}(U)=0\). By (3.9), \(\mathrm{E}[a(Y)]=-c^{\prime}(\theta)/b^{\prime}(\theta) = -(1/\theta)/(-1) = 1/\theta\), that is \(\mathrm{E}(\log Y)=1/\theta\). Hence \[ \mathrm{E}(U) = \frac{1}{\theta} - \mathrm{E}(\log Y) = \frac{1}{\theta} - \frac{1}{\theta} = 0, \] which is the general result (3.14). Directly: substituting \(t=\log y\) gives \(\int_1^\infty (\log y)\theta y^{-\theta-1}dy = \int_0^\infty t\,\theta e^{-\theta t}dt = 1/\theta\), since \(\log Y\sim\mathrm{Exp}(\theta)\).

The information. Using \(\mathfrak{I}=\mathrm{var}(U) = [b^{\prime}(\theta)]^2\,\mathrm{var}[a(Y)]\) from Section 3.3, and \(b^{\prime}(\theta)=-1\), \[ \mathfrak{I} = \mathrm{var}(\log Y) = \frac{1}{\theta^2}, \] the variance of an \(\mathrm{Exp}(\theta)\) variate. The same value follows from (3.15) with \(b^{\prime\prime}(\theta)=0\) and \(c^{\prime\prime}(\theta)=-1/\theta^2\): \[ \mathrm{var}(U) = \frac{b^{\prime\prime}(\theta)c^{\prime}(\theta)}{b^{\prime}(\theta)} - c^{\prime\prime}(\theta) = 0 + \frac{1}{\theta^2} = \frac{1}{\theta^2}, \] and from the third form (3.16), \(\mathrm{var}(U)=-\mathrm{E}(U^{\prime})\), since \(U^{\prime} = dU/d\theta = -1/\theta^2\) is non-random here so \(-\mathrm{E}(U^{\prime})=1/\theta^2\).

For a sample. If \(Y_1,\dots,Y_N\) are independent Pareto variables with the same \(\theta\), the log-likelihoods add, so \[ U = \sum_{i=1}^{N}\left(\frac{1}{\theta}-\log Y_i\right) = \frac{N}{\theta} - \sum_{i=1}^{N}\log Y_i, \qquad \mathfrak{I} = \frac{N}{\theta^2}. \] Setting \(U=0\) gives the maximum likelihood estimator \(\hat\theta = N/\sum\log Y_i\), and the asymptotic standard error is \(\mathfrak{I}^{-1/2}=\theta/\sqrt{N}\).

Simulation confirms the two results for a single observation.

set.seed(11); theta <- 2.5
y <- exp(rexp(2e6, rate = theta))
U <- 1 / theta - log(y)
round(c(mean.U = mean(U), theory = 0,
        var.U = var(U), information = 1 / theta^2,
        mean.logY = mean(log(y)), theory.logY = 1 / theta), 5)

 mean.U      theory       var.U information   mean.logY theory.logY
-0.0002      0.0000      0.1603      0.1600      0.4002      0.4000

Part (a). From Exercise 3.2 the Gamma density with shape \(\alpha\) and rate \(\beta\) is \[ f(y;\alpha,\beta) = \frac{\beta^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}e^{-y\beta}, \qquad y>0. \] Put \(\alpha=1\) and \(\beta=\theta\). Since \(\Gamma(1)=\int_0^\infty e^{-t}dt = 1\) and \(y^{0}=1\), \[ f(y;1,\theta) = \frac{\theta^{1}}{\Gamma(1)}y^{0}e^{-y\theta} = \theta e^{-y\theta}, \] which is exactly the exponential density of Exercise 3.3(b). So \(\mathrm{Exp}(\theta) = G(1,\theta)\). The moments agree too: Exercise 3.2 gives \(\mathrm{E}(Y)=\alpha/\beta = 1/\theta\) and \(\mathrm{var}(Y)=\alpha/\beta^2 = 1/\theta^2\), the exponential moments. Note that the arrow from Gamma to Exponential in Figure 3.3 is labelled \(\alpha=\theta\), \(\beta=1\), which is the two substitutions interchanged: \(G(\theta,1)\) has density \(y^{\theta-1}e^{-y}/\Gamma(\theta)\), which is exponential only when \(\theta=1\). The correct specialisation is \(\alpha=1\), \(\beta=\theta\), as derived above.

Part (b). Work with the survivor function. Let \(X\sim U[0,1]\) and put \(Y=-\theta\log X\). Since \(0<X<1\) we have \(\log X<0\), so \(Y>0\) as required, and for \(y>0\) \[ P(Y>y) = P(-\theta\log X > y) = P\!\left(\log X < -\frac{y}{\theta}\right) = P\!\left(X < e^{-y/\theta}\right) = e^{-y/\theta}, \] the last step because \(P(X<x)=x\) on \((0,1)\) and \(e^{-y/\theta}\in(0,1)\). Differentiating, the density of \(Y\) is \[ f_Y(y) = -\frac{d}{dy}P(Y>y) = \frac{1}{\theta}e^{-y/\theta}, \qquad y>0, \] so \(Y\) is exponential with mean \(\theta\).

A caution about parametrisation is needed here, and it is the reason to state the result this way. The book’s own exponential density, in Exercise 3.3(b) and Figure 3.3, is \(f(y;\theta)=\theta e^{-y\theta}\), which has mean \(1/\theta\), whereas the density just derived has mean \(\theta\). The two agree only if \(\mathrm{Exp}(\theta)\) is read with \(\theta\) as the scale (the mean). Under the book’s rate parametrisation the correct transformation is \(Y = -\theta^{-1}\log X\), for which the same argument gives \(P(Y>y) = P(X<e^{-\theta y}) = e^{-\theta y}\) and hence density \(\theta e^{-y\theta}\). Equivalently, the statement as printed shows \(-\theta\log X \sim \mathrm{Exp}(1/\theta)\) in the book’s notation. I flag this rather than paper over it; the substance of the result – that \(-\log X\) is a standard exponential variate, which is the basis of inversion sampling – is unaffected.

Part (c)(i): Binomial to Poisson. The moment generating function of \(Y\sim\mathrm{Bin}(n,\pi)\) is \[ M_Y(t) = \mathrm{E}\left(e^{tY}\right) = \sum_{y=0}^n \binom{n}{y}(\pi e^t)^y(1-\pi)^{n-y} = \left(1-\pi+\pi e^t\right)^n . \] Let \(n\to\infty\) with the mean held fixed at \(\lambda = n\pi\), so \(\pi=\lambda/n \to 0\). Then \[ M_Y(t) = \left[1 + \frac{\lambda(e^t-1)}{n}\right]^n \longrightarrow \exp\left[\lambda\left(e^t-1\right)\right], \] using \((1+a/n)^n \to e^a\). The limit is the moment generating function of \(\mathrm{Po}(\lambda)\), since for \(Z\sim\mathrm{Po}(\lambda)\) \[ M_Z(t) = \sum_{z=0}^\infty e^{tz}\frac{\lambda^z e^{-\lambda}}{z!} = e^{-\lambda}\sum_{z=0}^\infty \frac{(\lambda e^t)^z}{z!} = \exp\left[\lambda(e^t-1)\right]. \] By the continuity theorem for moment generating functions, \(\mathrm{Bin}(n,\lambda/n)\to\mathrm{Po}(\lambda)\). This is the “law of small numbers”: many trials, each very unlikely to succeed, with a stable expected number of successes.

Part (c)(ii): Negative Binomial to Poisson. With \(Y\sim\mathrm{NBin}(r,\theta)\) counting failures before the \(r\)th success, as in Exercise 3.3(c), \[ M_Y(t) = \sum_{y=0}^\infty e^{ty}\binom{y+r-1}{r-1}\theta^r(1-\theta)^y = \theta^r\sum_{y=0}^\infty\binom{y+r-1}{r-1}\left[(1-\theta)e^t\right]^y = \left[\frac{\theta}{1-(1-\theta)e^t}\right]^r, \] valid for \((1-\theta)e^t<1\), using the negative binomial series \(\sum_y \binom{y+r-1}{r-1}u^y = (1-u)^{-r}\). Now let \(r\to\infty\) with \(\lambda = r(1-\theta)\) fixed, so \(1-\theta = \lambda/r \to 0\). Taking logarithms, \[ \log M_Y(t) = r\left[\log\left(1-\frac{\lambda}{r}\right) - \log\left(1-\frac{\lambda e^t}{r}\right)\right] = r\left[\left(-\frac{\lambda}{r}+O(r^{-2})\right) - \left(-\frac{\lambda e^t}{r}+O(r^{-2})\right)\right], \] which tends to \(\lambda(e^t-1)\). Hence \(M_Y(t)\to\exp[\lambda(e^t-1)]\) and \(\mathrm{NBin}(r,\theta)\to\mathrm{Po}(r(1-\theta))\). This is consistent with the moments found in Exercise 3.5: as \(\theta\to1\) the ratio \(\mathrm{var}(Y)/\mathrm{E}(Y) = 1/\theta \to 1\), the Poisson value, so the overdispersion vanishes in the limit.

Both limits are checked numerically by the largest absolute discrepancy between the probability functions, with \(\lambda=3\) throughout.

lam <- 3
k <- 0:30
bin.err <- sapply(c(10, 50, 200, 5000),
                  function(n) max(abs(dbinom(k, n, lam / n) - dpois(k, lam))))
nb.err  <- sapply(c(5, 25, 100, 2000),
                  function(r) max(abs(dnbinom(k, size = r, prob = 1 - lam / r) - dpois(k, lam))))
out <- rbind(Bin.to.Po = bin.err, NBin.to.Po = nb.err)
colnames(out) <- c("n or r = 10/5", "50/25", "200/100", "5000/2000")
signif(out, 3)

           n or r = 10/5   50/25 200/100 5000/2000
Bin.to.Po         0.0428 0.00702 0.00170  6.72e-05
NBin.to.Po        0.1690 0.03250 0.00792  3.92e-04

Part (d): the Normal limits. The Central Limit Theorem states that if \(X_1,\dots,X_n\) are independent and identically distributed with mean \(m\) and finite variance \(v\), then \(\sum X_j\) is asymptotically \(N(nm, nv)\). In each case below the trick is to write the distribution as a sum of independent identically distributed pieces, which is possible because all three families are closed under convolution (reproductive).

(i) Poisson. If \(\mu\) is a positive integer, \(Y\sim\mathrm{Po}(\mu)\) has the same distribution as \(X_1+\dots+X_\mu\) where the \(X_j\) are independent \(\mathrm{Po}(1)\) variables, each with mean \(1\) and variance \(1\) (Exercise 3.4(a)). By the Central Limit Theorem \(Y\) is asymptotically \(N(\mu,\mu)\) as \(\mu\to\infty\). For non-integer \(\mu\) the same conclusion follows by writing \(Y\) as a sum of \(n\) independent \(\mathrm{Po}(\mu/n)\) variables and letting \(n\to\infty\) suitably, or by a direct characteristic-function argument; the essential point is that \(\mu\) large means many independent contributions. The rule of thumb printed in Figure 3.3 is \(\mu>15\).

(ii) Binomial. Here the decomposition is exact and immediate: \(Y\sim\mathrm{Bin}(n,\pi)\) is the sum \(X_1+\dots+X_n\) of independent Bernoulli variables with \(\mathrm{E}(X_j)=\pi\) and \(\mathrm{var}(X_j)=\pi(1-\pi)\) (Exercise 3.7). The Central Limit Theorem gives \(Y\) asymptotically \(N(n\pi, n\pi(1-\pi))\). The proviso matters: the Bernoulli summands are skew unless \(\pi=1/2\), and the standardised skewness of \(Y\) is \[ \frac{1-2\pi}{\sqrt{n\pi(1-\pi)}} , \] which is only small when \(n\pi(1-\pi)\) is large. If \(\pi\) is very small the boundary at \(0\) truncates the would-be Normal shape and the Poisson limit of part (c)(i) is the appropriate approximation instead. Hence the conditions \(n\pi>5\) and \(n\pi(1-\pi)>5\) in Figure 3.3.

(iii) Gamma. If \(\alpha\) is a positive integer, \(Y\sim G(\alpha,\beta)\) is the sum of \(\alpha\) independent \(\mathrm{Exp}(\beta) = G(1,\beta)\) variables by part (a), each with mean \(1/\beta\) and variance \(1/\beta^2\). The Central Limit Theorem then gives \(Y\) asymptotically \(N(\alpha/\beta, \alpha/\beta^2)\) as \(\alpha\to\infty\), matching the moments found in Exercise 3.2. For general real \(\alpha>0\), write \(\alpha = n(\alpha/n)\) and use \(G(\alpha,\beta) = \sum_{j=1}^n G(\alpha/n,\beta)\) independently, which is legitimate because the Gamma family is closed under convolution in the shape parameter.

In all three cases the underlying reason is the same: the coefficient of variation is \(O(\alpha^{-1/2})\), so as the shape or count parameter grows, the distribution concentrates and its standardised form flattens out to the Normal. The approximations are checked by the largest absolute difference between the exact and Normal cumulative distribution functions, with a continuity correction in the two discrete cases.

cc <- function(q, m, v) pnorm(q + 0.5, m, sqrt(v))
po <- sapply(c(5, 20, 100), function(mu) { k <- 0:ceiling(mu + 10 * sqrt(mu))
                 max(abs(ppois(k, mu) - cc(k, mu, mu))) })
bi <- sapply(c(10, 50, 500), function(n) { k <- 0:n
                 max(abs(pbinom(k, n, 0.3) - cc(k, n * 0.3, n * 0.3 * 0.7))) })
ga <- sapply(c(2, 10, 100), function(a) { q <- seq(0, a + 10 * sqrt(a), length = 2000)
                 max(abs(pgamma(q, shape = a, rate = 1) - pnorm(q, a, sqrt(a)))) })
out <- rbind(Po.mu = po, Bin.n = bi, Gamma.alpha = ga)
colnames(out) <- c("small", "medium", "large")
round(out, 4)

             small medium  large
Po.mu       0.0290 0.0148 0.0066
Bin.n       0.0177 0.0081 0.0026
Gamma.alpha 0.0945 0.0421 0.0133

The parameter values are \(\mu = 5, 20, 100\) for the Poisson, \(n = 10, 50, 500\) with \(\pi=0.3\) for the Binomial, and \(\alpha = 2, 10, 100\) with \(\beta=1\) for the Gamma. In every row the error falls roughly as the square root of the parameter, the rate the Berry-Esseen bound predicts, and at \(\mu=20\), \(n=50\) and \(\alpha=100\) the maximum error is at or below about \(1.5\%\), which is why the rules of thumb quoted in Figure 3.3 are set where they are.

The practical relevance for this chapter is that these limits are what license the large-sample approximations used throughout the book: the Normal approximations to the sampling distributions of estimators in Chapter 4, the \(\chi^2\) approximations to deviances in Chapter 5, and the interchangeable use of Poisson and Binomial models for rare events.

Part (a) — the raw series.

library(dobson)
data(aids)
y <- aids$cases
i <- seq_along(y)
plot(i, y, pch = 19, xlab = "time period i (quarters from 1984)",
     ylab = expression(y[i]~"(cases of AIDS)"),
     main = "AIDS cases in Australia by quarter, 1984-1988")
lines(i, y, lty = 3)

The counts rise monotonically apart from small wobbles at \(i = 7,8\) (1985 quarters 3–4) and \(i = 15\) (1987 quarter 3). The growth is clearly not linear: the increment per quarter is about 5 cases early on and about 15 cases at the end. The scatter about any smooth trend also widens with the level of the series, which is exactly the mean–variance behaviour \(\operatorname{var}(Y_i) = \operatorname{E}(Y_i)\) that motivates a Poisson model here.

Part (b) — the log–log plot.

Under \(\lambda_i = i^{\theta}\) we have \(\log \lambda_i = \theta \log i\), so a plot of \(\log y_i\) against \(\log i\) should be roughly a straight line through the origin with slope \(\theta\).

plot(log(i), log(y), pch = 19, xlab = expression(log~i),
     ylab = expression(log~y[i]), main = "log-log plot of AIDS counts")
abline(lm(log(y) ~ log(i)), lty = 2)

The points lie close to a straight line, so the power-law mean \(\lambda_i = i^{\theta}\) is a reasonable description; an ordinary least-squares line through the log–log points has slope \(1.454\) and intercept \(0.686\). The intercept is not zero, which is why part (c) generalises the model to \(\log \lambda_i = \beta_1 + \beta_2 x_i\) with \(x_i = \log i\): the one-parameter model \(\lambda_i = i^{\theta}\) is the special case \(\beta_1 = 0\). Note also that the first observation (\(i = 1\), \(y_1 = 1\)) sits at the origin of this plot and is the point furthest below the line; the log transformation exaggerates its influence, which is another reason to fit the Poisson GLM properly rather than to regress \(\log y_i\) on \(\log i\).

Part (c) — iterative weighted least squares from first principles.

The model is \(Y_i \sim \operatorname{Po}(\lambda_i)\) independently, with \(\eta_i = \log \lambda_i = \beta_1 + \beta_2 x_i = \mathbf{x}_i^T \boldsymbol\beta\) and \(x_i = \log i\). Write \(\mu_i = \lambda_i\) for the mean. Then

\[ \mu_i = e^{\eta_i}, \qquad \frac{\partial \mu_i}{\partial \eta_i} = e^{\eta_i} = \mu_i, \qquad \operatorname{var}(Y_i) = \mu_i . \]

From equation (4.23) the diagonal weights are

\[ w_{ii} = \frac{1}{\operatorname{var}(Y_i)} \left( \frac{\partial \mu_i}{\partial \eta_i} \right)^{2} = \frac{\mu_i^{2}}{\mu_i} = \mu_i , \]

so \(\mathbf{W} = \operatorname{diag}(\mu_1, \ldots, \mu_N)\) — a consequence of the log link being the canonical link for the Poisson distribution, for which the weights reduce to the variance function itself. From equation (4.24) the working response is

\[ z_i = \sum_{k=1}^{p} x_{ik} b_k^{(m-1)} + (y_i - \mu_i) \frac{\partial \eta_i}{\partial \mu_i} = \eta_i + \frac{y_i - \mu_i}{\mu_i} , \]

with \(\eta_i\) and \(\mu_i\) evaluated at \(\mathbf{b}^{(m-1)}\). The design matrix is \(\mathbf{X}\) with rows \((1, \log i)\), and the information matrix is \(\mathfrak{I} = \mathbf{X}^T \mathbf{W} \mathbf{X}\) with entries

\[ \mathfrak{I} = \begin{bmatrix} \sum_i \mu_i & \sum_i \mu_i x_i \\ \sum_i \mu_i x_i & \sum_i \mu_i x_i^{2} \end{bmatrix}, \qquad \mathbf{X}^T \mathbf{W} \mathbf{z} = \begin{bmatrix} \sum_i \mu_i \eta_i + (y_i - \mu_i) \\ \sum_i x_i \left[ \mu_i \eta_i + (y_i - \mu_i) \right] \end{bmatrix} . \]

Solving \(\mathbf{X}^T \mathbf{W} \mathbf{X} \mathbf{b}^{(m)} = \mathbf{X}^T \mathbf{W} \mathbf{z}\) repeatedly, starting from the log–log fit of part (b):

x <- log(i)
X <- cbind(1, x)
b <- c(0.5, 1.5)                          # rough start suggested by the log-log plot
for (m in 1:6) {
  eta <- X %*% b
  mu  <- exp(eta)                         # inverse log link
  W   <- diag(as.vector(mu))              # w_ii = (dmu/deta)^2 / var(Y) = mu_i
  z   <- eta + (y - mu)/mu                # z_i = eta_i + (y_i - mu_i) deta/dmu
  b   <- solve(t(X) %*% W %*% X, t(X) %*% W %*% z)
  cat(sprintf("m = %d   b1 = %.6f   b2 = %.6f\n", m, b[1], b[2]))
}

m = 1   b1 = 1.054285   b2 = 1.305697
m = 2   b1 = 0.996735   b2 = 1.326344
m = 3   b1 = 0.995998   b2 = 1.326610
m = 4   b1 = 0.995998   b2 = 1.326610
m = 5   b1 = 0.995998   b2 = 1.326610
m = 6   b1 = 0.995998   b2 = 1.326610

Convergence is essentially complete after three iterations, as in Table 4.4 for the worked example of Section 4.4. The maximum likelihood estimates are \(\widehat\beta_1 = 0.99600\) and \(\widehat\beta_2 = 1.32661\). The inverse of the information matrix at the estimates gives the standard errors, by equation (4.12) generalised to the vector case:

I <- t(X) %*% diag(as.vector(exp(X %*% b))) %*% X
I
solve(I)
sqrt(diag(solve(I)))

                  x
  1311.000 3396.379
x 3396.379 9038.301
                        x
   0.02880067 -0.01082261
x -0.01082261  0.00417752
                    x
0.16970761 0.06463374

Part (d) — the same model via Poisson regression software.

aids$i <- seq_len(nrow(aids))
res.p <- glm(cases ~ log(i), family = poisson(link = "log"), data = aids)
summary(res.p)

Call:
glm(formula = cases ~ log(i), family = poisson(link = "log"),
    data = aids)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.99600    0.16971   5.869 4.39e-09 ***
log(i)       1.32661    0.06463  20.525  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 677.264  on 19  degrees of freedom
Residual deviance:  21.755  on 18  degrees of freedom
AIC: 138.05

Number of Fisher Scoring iterations: 4

(Note the explicit poisson(link = "log"): after library(dobson) the bare name poisson refers to the package’s data set of Table 4.3, not to the family function.)

The estimates and standard errors agree with part (c) to every printed digit — \(b_1 = 0.99600\,(0.16971)\), \(b_2 = 1.32661\,(0.06463)\) — which is what we expect, since glm uses precisely the algorithm of equation (4.25) and, for the canonical link, Fisher scoring and Newton–Raphson coincide.

On interpretation and adequacy: the fitted mean is \(\widehat\lambda_i = e^{0.996} i^{1.327} = 2.71\, i^{1.327}\): the quarterly incidence grows as a power of time with exponent about \(1.33\), not exponentially in \(i\). An approximate 95% confidence interval for \(\beta_2\) is \(1.3266 \pm 1.96 \times 0.0646 = (1.200, 1.453)\), which excludes \(1\), so the growth is faster than linear in \(i\) but far short of the doubling behaviour that “increasing exponentially” would suggest. The one-parameter model of part (b), \(\beta_1 = 0\), is decisively rejected (\(z = 5.87\)).

The residual deviance is \(21.755\) on \(18\) degrees of freedom, \(p = 0.243\) against \(\chi^2(18)\), so there is no evidence of lack of fit or of overdispersion. Compare the alternative log-linear-in-time model \(\log \lambda_i = \beta_1 + \beta_2 i\), which is the literal “exponential growth” hypothesis:

m.pow <- glm(cases ~ log(i), family = poisson(link = "log"), data = aids)
m.exp <- glm(cases ~ i,      family = poisson(link = "log"), data = aids)
c(power.dev = deviance(m.pow), power.AIC = AIC(m.pow),
  exp.dev = deviance(m.exp), exp.AIC = AIC(m.exp))

power.dev power.AIC   exp.dev   exp.AIC
 21.75511 138.05303  53.02000 169.31792

The power-law model is far better (deviance \(21.8\) versus \(53.0\) on the same \(18\) degrees of freedom; the exponential-growth model has \(p = 2.6 \times 10^{-5}\) for lack of fit, and its AIC is \(31\) units worse). Substantively, the epidemic was already decelerating on the log scale by 1988.

par(mfrow = c(1, 2))
plot(aids$i, aids$cases, pch = 19, xlab = "quarter i", ylab = "cases",
     main = "Fitted Poisson power-law mean")
lines(aids$i, fitted(m.pow), lwd = 2)
plot(aids$i, residuals(m.pow, type = "deviance"), pch = 19,
     xlab = "quarter i", ylab = "deviance residual", main = "Deviance residuals")
abline(h = 0, lty = 2)

The fitted curve tracks the counts closely and the deviance residuals show no trend and no funnelling; they run from \(-2.06\) (at \(i = 8\), the 1985 fourth quarter, where the series briefly flattens) to \(1.73\), so a single residual just crosses \(-2\) and nothing else is close to the boundary. The model is adequate.

Each \(y_i\) is time to death for a particular leukemia patient, in weeks from diagnosis; the \(x_i\) are \(\log_{10}\) initial white blood cell counts.

Part (a) — scatter plot.

library(dobson)
data(leukemia)
x <- leukemia$wbc
y <- leukemia$time
plot(x, y, pch = 19, xlab = "x (log10 initial white blood cell count)",
     ylab = "y (survival time, weeks)",
     main = "Survival time against log10 WBC, 17 leukemia patients")
lines(lowess(x, y), lty = 2)

length(x)  # 17 patients

[1] 17

The data show a clear negative trend: as the initial white cell count rises, survival time falls. Two further features matter for the modelling. First, the spread of \(y\) shrinks with \(x\) — at \(x \approx 3\) survival ranges over roughly \(16\) to \(156\) weeks, while at \(x = 5\) the three patients survived \(1\), \(1\) and \(65\) weeks. Variability therefore scales with the level, ruling out constant-variance least squares. Second, the decline looks curved rather than straight, consistent with a multiplicative (exponential) rather than additive effect of \(x\).

Part (b) — link function.

\[ \operatorname{E}(Y_i) = \exp(\beta_1 + \beta_2 x_i) \iff \log \operatorname{E}(Y_i) = \beta_1 + \beta_2 x_i , \]

so the appropriate link is the log link, \(g(\mu) = \log \mu\), in the notation \(g(\mu_i) = \mathbf{x}_i^T \boldsymbol\beta\) of equation (4.16).

This is the natural choice here for three reasons. Survival time is strictly positive, and the exponential mean function guarantees \(\operatorname{E}(Y) > 0\) for all \(\beta\) and all \(x\), unlike an identity link which can predict negative survival times. The variance scales with the mean (survival times are right-skewed), so a constant-variance assumption is inappropriate; the Gamma/exponential family with a log link handles this directly. Finally the log link linearises a multiplicative relationship — a unit increase in \(x_i\) multiplies \(\operatorname{E}(Y)\) by \(e^{\beta_2}\) — which is the more natural form for this biological relationship. The negative trend in (a) implies \(\beta_2 < 0\).

Part (c) — moments of the exponential distribution.

\[ f(y; \theta) = \theta e^{-y\theta}, \qquad y > 0, \]

which is the Gamma density with shape parameter \(\phi = 1\). For the mean, integrate by parts with \(u = y\) and \(\mathrm{d}v = \theta e^{-y\theta}\,\mathrm{d}y\):

\[ \operatorname{E}(Y) = \int_0^{\infty} y \cdot \theta e^{-y\theta}\,\mathrm{d}y = \left[ -y e^{-y\theta} \right]_0^{\infty} + \int_0^{\infty} e^{-y\theta}\,\mathrm{d}y = 0 + \left[ -\frac{1}{\theta} e^{-y\theta} \right]_0^{\infty} = \frac{1}{\theta}. \]

For the second moment, two applications of integration by parts (or the Gamma integral \(\int_0^{\infty} y^n e^{-y\theta}\,\mathrm{d}y = n! / \theta^{n+1}\)) give

\[ \operatorname{E}(Y^{2}) = \theta \cdot \frac{2!}{\theta^{3}} = \frac{2}{\theta^{2}}, \]

whence

\[ \operatorname{var}(Y) = \operatorname{E}(Y^{2}) - \left[ \operatorname{E}(Y) \right]^{2} = \frac{2}{\theta^{2}} - \frac{1}{\theta^{2}} = \frac{1}{\theta^{2}} . \]

Note that \(\operatorname{var}(Y) = \left[ \operatorname{E}(Y) \right]^{2}\), confirming the Gamma variance structure \(\operatorname{var}(Y) = \phi \left[ \operatorname{E}(Y) \right]^{2}\) with \(\phi = 1\). This is the fact that part (e) exploits: the standard deviation of \(Y_i\) equals its mean.

Part (d) — fitting the model.

fit <- glm(y ~ x, family = Gamma(link = "log"), data = data.frame(x, y))
summary(fit)

Call:
glm(formula = y ~ x, family = Gamma(link = "log"), data = data.frame(x,
    y))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   8.4775     1.6034   5.287 9.13e-05 ***
x            -1.1093     0.3872  -2.865   0.0118 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.9388638)

    Null deviance: 26.282  on 16  degrees of freedom
Residual deviance: 19.457  on 15  degrees of freedom
AIC: 173.97

Number of Fisher Scoring iterations: 8

We use Gamma(link = "log") because the exponential distribution is the special case of the Gamma family with \(\phi = 1\), and the log link is the mean function \(\operatorname{E}(Y_i) = \exp(\beta_1 + \beta_2 x_i)\) of part (b). The maximum likelihood estimates of \(\boldsymbol\beta\) do not depend on \(\phi\) at all — \(\phi\) cancels out of the score equations (4.18), since it enters \(\operatorname{var}(Y_i)\) as a constant multiplier — so these are exactly the exponential-model estimates. Only the standard errors are affected, and by default glm estimates \(\phi\) from the data rather than fixing it at \(1\). Imposing \(\phi = 1\), as the exponential assumption requires:

summary(fit, dispersion = 1)$coefficients

             Estimate Std. Error   z value     Pr(>|z|)
(Intercept)  8.477494  1.6548077  5.122948 3.007955e-07
x           -1.109297  0.3996545 -2.775640 5.509317e-03

The standard errors change only slightly (from \(1.603\) to \(1.655\) and from \(0.387\) to \(0.400\)), because \(\widehat\phi = 0.939\) is already close to \(1\).

On interpretation: \(\widehat\beta_2 = -1.109\), confirming the negative trend of part (a). On the original scale, an increase of one unit in \(\log_{10}\) (white cell count) — that is, a tenfold increase in the count — multiplies expected survival by \(e^{-1.109} = 0.33\), so it cuts expected survival to about a third. A 95% confidence interval for \(\beta_2\) using the exponential (\(\phi = 1\)) standard error is \(-1.109 \pm 1.96 \times 0.400 = (-1.893, -0.326)\), giving a multiplicative factor between \(0.15\) and \(0.72\); the effect is clearly established in direction, though its size is imprecisely determined with only 17 patients. The fitted mean survival runs from \(\widehat y = 197\) weeks at \(x = 2.88\) down to \(\widehat y = 19\) weeks at \(x = 5.00\).

Part (e) — standardized residuals.

For the exponential (Gamma with \(\phi = 1\)) model, part (c) gives \(\operatorname{var}(Y_i) = \left[ \operatorname{E}(Y_i) \right]^{2}\), so the standard deviation equals the mean: \(\operatorname{sd}(Y_i) = \operatorname{E}(Y_i)\). Dividing the raw residual by \(\widehat y_i\) therefore divides by the estimated standard deviation, which is what makes

\[ r_i = \frac{y_i - \widehat y_i}{\widehat y_i}, \qquad \widehat y_i = \exp(\widehat\beta_1 + \widehat\beta_2 x_i) \]

a genuinely standardized residual rather than an ad hoc scaling. These are exactly the Pearson residuals for the Gamma family.

b1 <- coef(fit)[1]
b2 <- coef(fit)[2]

yhat <- exp(b1 + b2*x)     # equivalently fitted(fit)
r    <- (y - yhat) / yhat  # equivalently residuals(fit, type = "pearson")

par(mfrow = c(1, 2))
plot(yhat, r, xlab = expression(hat(y)[i]), ylab = expression(r[i]),
     main = "Residuals vs fitted", pch = 19); abline(h = 0, lty = 2)
plot(x, r, xlab = "x (log10 WBC)", ylab = expression(r[i]),
     main = "Residuals vs x", pch = 19); abline(h = 0, lty = 2)

data.frame(x, y, yhat = round(yhat, 2), r = round(r, 3))

      x   y   yhat      r
1  3.36  65 115.61 -0.438
2  2.88 156 196.90 -0.208
3  3.63 100  85.69  0.167
4  3.41 134 109.38  0.225
5  3.78  16  72.56 -0.779
6  4.02 108  55.60  0.943
7  4.00 121  56.84  1.129
8  4.23   4  44.04 -0.909
9  3.73  39  76.69 -0.491
10 3.85 143  67.13  1.130
11 3.97  56  58.77 -0.047
12 4.51  26  32.28 -0.195
13 4.54  22  31.23 -0.295
14 5.00   1  18.75 -0.947
15 5.00   1  18.75 -0.947
16 4.72   5  25.57 -0.804
17 5.00  65  18.75  2.467