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

\begin{equation*} \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). \end{equation*}

Indeed \(\mathbf{w} = \mathbf{A}\mathbf{y}\) with \(\mathbf{y} = [Y_1, Y_2]^T\) and \(\mathbf{A} = \begin{pmatrix} 1 & 2 \\ 4 & -1 \end{pmatrix}\), so every linear combination of \(W_1, W_2\) is a linear combination of \(Y_1, Y_2\) and hence Normal (Section 1.4.1 result 4): \(\mathbf{w}\) is bivariate Normal. Independence gives \(\boldsymbol{\mu}_y = (1, 2)^T\) and \(\mathbf{V}_y = \mathrm{diag}(3, 5)\), so by Equations (1.2) and (1.3),

\begin{align*} E(\mathbf{w}) &= \mathbf{A}\boldsymbol{\mu}_y = (1 + 4,\; 4 - 2)^T = (5, 2)^T, \\ \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) &= 4\,\mathrm{var}(Y_1) - 2\,\mathrm{var}(Y_2) = 4(3) - 2(5) = 2, \end{align*}

the cross terms vanishing because \(\mathrm{cov}(Y_1, Y_2) = 0\). In particular \(W_1\) and \(W_2\) are not independent, the correlation being \(2/\sqrt{23 \times 53} = 0.057\).

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

(a) \(Y_1^2 \sim \chi^2(1)\), since \(Y_1\) is already standard Normal (Section 1.4.2 result 1, \(n = 1\)).

(b) Since \(Y_2 \sim N(3,4)\) gives \((Y_2 - 3)/2 \sim N(0,1)\) independent of \(Y_1\),

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

a sum of squares of two independent standard Normals, as in Equation (1.4) with \(n = 2\).

(c) Independence makes \(\mathbf{V} = \mathrm{diag}(1, 4)\) and \(\mathbf{V}^{-1} = \mathrm{diag}(1, 1/4)\), with \(\boldsymbol{\mu} = (0, 3)^T\), so

\[ \mathbf{y}^T\mathbf{V}^{-1}\mathbf{y} = Y_1^2 + \frac{Y_2^2}{4} \sim \chi^2(2,\, 2.25) \]

by Section 1.4.2 result 6 (\(\mathbf{V}\) non-singular), the non-centrality parameter being

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

Not centring \(\mathbf{y}\) is exactly what leaves \(\lambda = (\mu_2/\sigma_2)^2\) behind.

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.

Since \(\det \mathbf{V} = 4(9) - 1 = 35 > 0\) with \(v_{11} = 4 > 0\), \(\mathbf{V}\) is positive definite (Section 1.5 result 3) and

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

(a) Expanding the quadratic form, the cross term doubling because the off-diagonal entries agree,

\begin{equation*} (\mathbf{y} - \boldsymbol{\mu})^T\mathbf{V}^{-1}(\mathbf{y} - \boldsymbol{\mu}) = \frac{9(y_1 - 2)^2 - 2(y_1 - 2)(y_2 - 3) + 4(y_2 - 3)^2}{35}, \end{equation*}

which is \(\chi^2(2)\) by Equation (1.5), \(\mathbf{V}\) being non-singular.

(b) The same algebra without centring gives

\begin{equation*} \mathbf{y}^T\mathbf{V}^{-1}\mathbf{y} = \frac{1}{35}\left[ 9y_1^2 - 2y_1y_2 + 4y_2^2 \right] \sim \chi^2(2,\, 12/7), \end{equation*}

by Section 1.4.2 result 6, with

\begin{equation*} \lambda = \boldsymbol{\mu}^T\mathbf{V}^{-1}\boldsymbol{\mu} = \frac{9(2)^2 - 2(2)(3) + 4(3)^2}{35} = \frac{60}{35} = \frac{12}{7} \approx 1.714. \end{equation*}

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

Both forms carry \(2\) degrees of freedom because \(\mathbf{V}^{-1}\) has full rank (Section 1.5 result 4); only the location differs, (a) having expectation \(2\) and (b) expectation \(n + \lambda = 2 + 12/7 \approx 3.71\).

(a) \(\overline{Y} = \sum a_i Y_i\) with every \(a_i = 1/n\), so \(\overline{Y} \sim N(\mu, \sigma^2/n)\) by Section 1.4.1 result 4.

(b) Write \(Y_i - \overline{Y} = (Y_i - \mu) - (\overline{Y} - \mu)\) and expand; since \(\sum (Y_i - \mu) = n(\overline{Y} - \mu)\) the cross term is \(-2n(\overline{Y} - \mu)^2\), so

\begin{align*} \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 \\ &= \sum (Y_i - \mu)^2 - n(\overline{Y} - \mu)^2, \end{align*}

and dividing by \(n - 1\) gives the identity (purely algebraic; no distributional assumption).

(c) Dividing by \(\sigma^2\),

\begin{equation*} \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}, \end{equation*}

where \(Q \sim \chi^2(n)\) by Equation (1.4) and \(Q_1, Q_2\) are quadratic forms in the independent \(Y_i - \mu \sim N(0, \sigma^2)\) of ranks \(m_1 = n - 1\) (matrix \(\mathbf{I} - \frac{1}{n}\mathbf{J}\), idempotent with trace \(n-1\)) and \(m_2 = 1\). As \(m_1 + m_2 = n\), Cochran’s theorem (Section 1.5 result 5) makes \(Q_1 \sim \chi^2(n-1)\) and \(Q_2 \sim \chi^2(1)\) independent. That delivers independence of \(S^2\) and \((\overline{Y} - \mu)^2\); for \(\overline{Y}\) itself, \(\mathrm{cov}(\overline{Y}, Y_i - \overline{Y}) = \sigma^2/n - \sigma^2/n = 0\) for every \(i\) and \((\overline{Y}, Y_1 - \overline{Y}, \ldots, Y_n - \overline{Y})\) is multivariate Normal (all entries linear in the \(Y_i\)), so \(\overline{Y}\) is independent of the residual vector and hence of \(S^2\).

(d) \((n-1)S^2/\sigma^2 \sim \chi^2(n-1)\), by (c); in particular \(E(S^2) = \sigma^2\).

(e) With \(Z = (\overline{Y} - \mu)/(\sigma/\sqrt{n}) \sim N(0,1)\) from (a) and \(X^2 = (n-1)S^2/\sigma^2 \sim \chi^2(n-1)\) independent of \(Z\) by (c), so that \(S/\sigma = [X^2/(n-1)]^{1/2}\),

\begin{equation*} \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}} \sim t(n-1) \end{equation*}

by the definition of the t-distribution in Equation (1.6).

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); and the t-statistic has simulated variance \(1.400\) against \((n-1)/(n-3) = 1.4\) with Kolmogorov-Smirnov \(p = 0.49\) against \(t(7)\), confirming (e).

(a) With \(f(y_i; \theta) = \theta^{y_i}e^{-\theta}/y_i!\), the \(y = 0\) term dropping and \(z = y - 1\),

\begin{equation*} E(Y_i) = \sum_{y=1}^{\infty} \frac{\theta^{y}e^{-\theta}}{(y-1)!} = \theta e^{-\theta}\sum_{z=0}^{\infty}\frac{\theta^{z}}{z!} = \theta e^{-\theta}e^{\theta} = \theta. \end{equation*}

(b) \(\widehat{\beta} = \log \overline{y}\) (for \(\overline{y} > 0\)). Reparameterising the log-likelihood of Section 1.6.2 by \(\theta = e^{\beta}\),

\begin{equation*} l(\beta; \mathbf{y}) = \left( \sum y_i \right)\beta - n e^{\beta} - \sum \log y_i!, \qquad \frac{dl}{d\beta} = \sum y_i - n e^{\beta}, \end{equation*}

which vanishes at \(e^{\beta} = \overline{y}\); and \(d^2l/d\beta^2 = -ne^{\beta} < 0\) everywhere makes \(l\) strictly concave, so the root is the unique maximum. (Equivalently, invariance of maximum likelihood under the one-to-one \(\beta = \log\theta\).)

(c) \(\widetilde{\beta} = \log \overline{y}\) also. Putting \(u = e^{\beta}\), the chain rule gives

\begin{equation*} \frac{dS}{d\beta} = \left[ -2\sum (Y_i - u) \right] e^{\beta} = -2e^{\beta}\left( \sum Y_i - n e^{\beta} \right), \end{equation*}

which, since \(e^{\beta} > 0\), vanishes only at \(\sum Y_i = ne^{\beta}\); there \(dS/du = 0\), so \(d^2S/d\beta^2 = e^{2\beta}(2n) > 0\) and the point is a minimum. The two estimators coincide because both estimating equations reduce to \(\sum(Y_i - e^{\beta}) = 0\) (comment 2 of Section 1.6.4) – an agreement that depends on the \(Y_i\) sharing one \(\theta\), so that unweighted least squares is already correctly weighted (Section 1.6.3).

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).

(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 being the least informative (group \(14\) has \(10\) progeny, so its \(0.400\) is one female from \(0.500\)).

(b) \(\widehat{\theta} = T/N = \sum y_i / \sum n_i\), the pooled proportion rather than the average of the \(16\) group proportions. Writing \(N = \sum n_i\), \(T = \sum y_i\), and dropping the \(\theta\)-free binomial coefficients from the log-likelihood of the product over the \(16\) independent groups,

\begin{align*} l(\theta; \mathbf{y}) &= \mathrm{const} + T\log\theta + (N - T)\log(1 - \theta), \\ \frac{dl}{d\theta} &= \frac{T}{\theta} - \frac{N - T}{1 - \theta} = \frac{T - N\theta}{\theta(1 - \theta)}, \end{align*}

following Equation (1.9), whose numerator vanishes at \(T/N\); and \(d^2l/d\theta^2 = -T/\theta^2 - (N-T)/(1-\theta)^2 < 0\) on \((0,1)\) whenever \(0 < T < N\), so this is the unique interior maximum. Here \(T = 363\), \(N = 734\) and \(\widehat{\theta} = 363/734 = 0.4946\).

(c) Three numerical routes, none using 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 agree with (b) to six decimal places; the bisection differs in the seventh, \(20\) halvings of an interval of width \(0.4\) leaving a resolution of about \(4 \times 10^{-7}\).

A standard error and a check on the single-\(\theta\) model:

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 sex ratio is \(\widehat{\theta} = 0.495\) with standard error \(0.018\) and approximate \(95\%\) interval \((0.458, 0.531)\), which contains \(0.5\); the test of \(H_0: \theta = 0.5\) gives \(z = -0.30\), \(p = 0.77\), so there is no evidence of departure from an even sex ratio. The residual deviance \(11.93\) on \(15\) degrees of freedom (\(16\) groups less one fitted parameter) gives \(p = 0.68\), so no evidence of heterogeneity between mothers either: the spread of the \(16\) proportions is no more than binomial sampling variation. The figure shows this, the proportions fanning out for small \(n_i\) as the variance \(\theta(1-\theta)/n_i\) requires, all \(16\) 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 shows a single sharply curved peak at \(\widehat{\theta} = 0.4946\); the curvature there is \(-d^2l/d\theta^2 = N/[\theta(1-\theta)] = 2936.3\), and \(1/\sqrt{2936.3} = 0.01845\) reproduces the standard error above.

There is no detectable treatment effect: \(F = 0.260\) on \((1,38)\) degrees of freedom, \(p = 0.61\). The data ship in the dobson package as plants, stacked once and reused 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.

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 in \(3.05\)–\(6.34\) g, so no obvious data-entry errors; both distributions are single-peaked, roughly symmetric, with near-straight Normal probability plots, so Normality is credible. The spreads match (\(s_1 = 0.791\), \(s_2 = 0.864\); variance-ratio \(F_{19,19} = 0.839\), \(p = 0.71\)), supporting the common-\(\sigma^2\) assumption. The treatment mean exceeds the control mean by \(4.860 - 4.726 = 0.134\) g, an eighth of a standard deviation, with near-complete overlap of the dot plots: any treatment effect is small relative to plant-to-plant variation.

(b) The unpaired t-test, pooled-variance to match the common-\(\sigma^2\) model.

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 difference is \(\overline{y}_1 - \overline{y}_2 = 0.1335\) g with \(t = 0.510\) on \(38\) degrees of freedom, \(p = 0.61\): no evidence against \(\mathrm{H}_0\). The \(95\%\) interval \((-0.397, 0.664)\) contains zero and is wide – the data are consistent with the enriched environment lowering mean dry weight by \(0.40\) g or raising it by \(0.66\) g – so the result is inconclusive rather than a demonstration of no effect.

(c) \(\widehat{\mu}_j = \overline{Y}_j\) and \(\widehat{\mu} = \overline{Y}\), by maximum likelihood and least squares alike. With \(\sigma^2\) known, \(\ell_1 = -\tfrac{1}{2}N\log(2\pi\sigma^2) - S_1/(2\sigma^2)\) where \(S_1 = \sum_j\sum_k (y_{jk}-\mu_j)^2\) and \(N = 40\), so the two criteria give the same equations (Equations (2.8) and (2.10)); then

\begin{align*} \frac{\partial \ell_1}{\partial \mu_j} &= \frac{1}{\sigma^2}\sum_{k=1}^{20}(y_{jk}-\mu_j) = 0 \;\Longrightarrow\; \widehat{\mu}_j = \overline{Y}_j, \\ \frac{\partial \ell_0}{\partial \mu} &= \frac{1}{\sigma^2}\sum_{j}\sum_{k}(y_{jk}-\mu) = 0 \;\Longrightarrow\; \widehat{\mu} = \overline{Y}, \end{align*}

with \(\partial^2\ell_1/\partial\mu_j^2 = -20/\sigma^2 < 0\) and \(\partial^2\ell_0/\partial\mu^2 = -40/\sigma^2 < 0\), so both are maxima. Numerically \(\widehat{\mu}_1 = 4.8600\), \(\widehat{\mu}_2 = 4.7265\), \(\widehat{\mu} = 4.79325\).

(d) Substituting (c) into the criteria gives \(\widehat{S}_1 = \sum\sum (Y_{jk}-\overline{Y}_j)^2\) and \(\widehat{S}_0 = \sum\sum (Y_{jk}-\overline{Y})^2\); these are global minima because, the cross-product vanishing by \(\sum_k (Y_{jk}-\overline{Y}_j) = 0\),

\begin{equation*} \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 \end{equation*}

with equality only at \(\mu_j = \overline{Y}_j\), and the same identity over all \(40\) observations gives \(\widehat{S}_0\).

(e) Two printing slips: in the first display the subtracted term is summed over \(j = 1,2\) (not \(k\)), and in the third it is \(\frac{40}{\sigma^2}(\overline{Y}-\mu)^2\) with no sum over \(j\). Applying Exercise 1.4(b) within group \(j\) (\(n = 20\), mean \(\mu_j\)), summing over \(j\) and dividing by \(\sigma^2\),

\begin{equation*} \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 . \end{equation*}

Under \(\mathrm{H}_1\) the first term is \(\chi^2(40)\) and, as \(\overline{Y}_j \sim \mathrm{N}(\mu_j, \sigma^2/20)\) independently across \(j\), the second is \(\chi^2(2)\), independent of \(\widehat{S}_1/\sigma^2\) by Exercise 1.4(c); comparing moment generating functions in \(\chi^2(40) = \widehat{S}_1/\sigma^2 + \chi^2(2)\) gives \(\widehat{S}_1/\sigma^2 \sim \chi^2(38) = \chi^2(40-2)\). Likewise with \(n = 40\) and the single mean \(\mu\),

\begin{equation*} \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 , \end{equation*}

where the second term is \(\chi^2(1)\) (since \(\overline{Y}\sim \mathrm{N}(\mu,\sigma^2/40)\)) and independent of \(\widehat{S}_0\), so \(\widehat{S}_0/\sigma^2 \sim \chi^2(39)\) under \(\mathrm{H}_0\).

(f) Expanding \(Y_{jk}-\overline{Y} = (Y_{jk}-\overline{Y}_j) + (\overline{Y}_j-\overline{Y})\) with \(\sum_k (Y_{jk}-\overline{Y}_j) = 0\), and using \(\overline{Y} = (\overline{Y}_1+\overline{Y}_2)/2\) for equal group sizes,

\begin{equation*} \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 . \end{equation*}

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)\) (Section 2.2.2 with \(J = 2\)), independent of \(\widehat{S}_1\) because the former depends on the group means alone and the latter on the within-group deviations alone. Hence by the definition of the \(F\)-distribution (Section 1.4.4),

\begin{equation*} 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). \end{equation*}

If \(\mathrm{H}_0\) is false then \(\mathrm{E}(\overline{Y}_1-\overline{Y}_2) = \mu_1-\mu_2 \neq 0\), the numerator is non-central \(\chi^2\) with \(\lambda = 10(\mu_1-\mu_2)^2/\sigma^2\), and \(F\) is non-central (Figure 2.5).

(g) Computing the two minima and the statistic 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 from the nested model comparison:

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

Separate group means reduce the sum of squares from \(26.232\) to \(26.053\), an improvement of \(0.178\) or \(0.7\)% of the total; \(F = 0.260\) on \((1,38)\), \(p = 0.61\), far below \(F_{0.95}(1,38) = 4.10\). There is no evidence against \(\mathrm{H}_0\), so the single-mean model is preferable and the nutritional enrichment has no demonstrable effect on mean dry weight.

(h) By Section 1.4.4, \(T \sim t(n)\) implies \(T^2 \sim F(1,n)\); here \(t = 0.50985\) gives \(t^2 = 0.25994 = F\), and the \(p\)-values agree exactly at \(0.6131\) because the two-sided \(t\)-test and the upper-tail \(F\)-test reject on the same event \(\{|T| > c\}\). The conclusions of (b) and (g) are identical, the \(t\) formulation additionally supplying the signed estimate and interval that \(F\) discards.

(i) Under \(\mathrm{H}_0\) every fitted value is \(\overline{y}\), so the residuals are \(y_{jk}-\overline{y}\), standardized by \(s = \sqrt{\widehat{S}_0/39}\).

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 Normal probability plot is near-straight with no residual beyond \(\pm 2.13\) and Shapiro-Wilk \(p = 0.85\), and the spread is the same in the two groups, supporting Normality and constant variance. The residual group means are \(+0.081\) and \(-0.081\) against a residual standard deviation of \(1\): no group structure remains, so the \(\mathrm{H}_0\) model is adequate.

The paired analysis rejects \(\mathrm{H}_0\) and the unpaired one does not; the paired analysis is the correct one, and the men retain a mean loss of \(2.645\) kg. The data ship in the dobson package as waist, one row per man.

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 points lie almost on \(y = x\) and nearly all just below it: the men differ enormously from each other (\(80\) to \(135\) kg) but each changes only a little.

(a) Treating the two sets of \(20\) weights as independent samples, as the 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\) with wide interval \((-5.68, 10.97)\): no evidence of a retained loss.

(b) With \(D_k = Y_{1k} - Y_{2k} \sim \mathrm{N}(\delta, \sigma_D^2)\) independent across men, \(\delta = \mu_1 - \mu_2\), and \(\sigma_D^2\) known, the models are

\begin{equation*} \mathrm{H}_0 : \mathrm{E}(D_k) = 0, \qquad \mathrm{H}_1 : \mathrm{E}(D_k) = \delta, \qquad k = 1,\ldots,20. \end{equation*}

As in Exercise 2.1(c), maximizing \(\ell_1\) is minimizing \(S_1 = \sum_k (d_k-\delta)^2\) and gives \(\widehat{\delta} = \overline{D}\), so the minimized criteria are \(\widehat{S}_0 = \sum_k D_k^2\) (no parameters) and \(\widehat{S}_1 = \sum_k (D_k-\overline{D})^2\) (one parameter). Exercise 1.4(b) with \(n = 20\) writes \(\widehat{S}_1/\sigma_D^2\) as a \(\chi^2(20)\) minus an independent \(\chi^2(1)\), so \(\widehat{S}_1/\sigma_D^2 \sim \chi^2(19)\) regardless of \(\mathrm{H}_0\); and under \(\mathrm{H}_0\), \(\widehat{S}_0 - \widehat{S}_1 = 20\overline{D}^2\) with \(\overline{D} \sim \mathrm{N}(0,\sigma_D^2/20)\) is \(\sigma_D^2\chi^2(1)\) independent of \(\widehat{S}_1\). Hence

\begin{equation*} 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.} \end{equation*}

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)\) gives \(p = 0.0097\), beyond \(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

Here \(t^2 = 2.8734^2 = 8.256 = F\) by the identity of Exercise 2.1(h), and the men retain a mean loss of \(2.645\) kg with \(95\%\) interval \((0.72, 4.57)\) kg.

(c) No — the conclusions are opposite. Both estimate the same \(2.645\) kg, but the unpaired standard error is \(4.11\) kg against the paired \(0.92\) kg, because of 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 \(\rho = 0.953\), \(\mathrm{var}(D_k) = 2\sigma^2(1-\rho)\) is about \(5\)% of the \(2\sigma^2\) the unpaired analysis assumes: pairing removes the dominant between-man variation, which is irrelevant to the question, so the unpaired test’s failure to reject is a Type II error.

(d) For (a): the \(40\) observations mutually independent, both sets Normal, equal variances, representative sample. Independence is the fatal one, the same \(20\) men being measured twice. For (b): the \(20\) differences mutually independent (different men); each \(D_k\) Normal with common \(\sigma_D^2\) (Q-Q acceptably straight, Shapiro-Wilk \(p = 0.16\), though with several exact zeros); and a common mean \(\delta\), i.e. retained loss independent of initial weight. Nothing is assumed about the between-man weight distribution or about equal variances across times.

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 loss on initial weight gives slope \(0.121\) kg per kg, \(p = 0.084\): weak evidence that heavier men retain a larger loss, so the common-\(\delta\) assumption is not clearly contradicted. The paired analysis of (b) is the appropriate one, the design being before-and-after on the same subjects, and the men retain a mean loss of about \(2.6\) kg (\(95\%\) CI \(0.7\) to \(4.6\) kg) – though with no control group the program effect cannot be separated from regression to the mean.

Fix a group \(j\), assume \(S_{xx} > 0\), and write \(S_{xx} = \sum_k (x_{jk}-\overline{x}_j)^2 = \sum_k x_{jk}^2 - K\overline{x}_j^2\), \(S_{xY} = \sum_k (x_{jk}-\overline{x}_j)(Y_{jk}-\overline{Y}_j)\), \(S_{YY} = \sum_k (Y_{jk}-\overline{Y}_j)^2\), so that \(b_j = S_{xY}/S_{xx}\) and \(a_j = \overline{Y}_j - b_j\overline{x}_j\) (Section 2.2.2).

(1) Substituting \(a_j = \overline{Y}_j - b_j\overline{x}_j\) gives \(Y_{jk}-a_j-b_j x_{jk} = (Y_{jk}-\overline{Y}_j) - b_j(x_{jk}-\overline{x}_j)\), so the group-\(j\) contribution to \(\widehat{S}_1\) is, using \(S_{xY} = b_j S_{xx}\),

\[ \widehat{S}_{1j} = S_{YY} - 2b_j S_{xY} + b_j^2 S_{xx} = S_{YY} - b_j^2 S_{xx}. \tag{i} \]

For the right-hand side put \(e_{jk} = Y_{jk} - (\alpha_j + \beta_j x_{jk})\), so \(\overline{e}_j = \overline{Y}_j - \alpha_j - \beta_j\overline{x}_j\) is the quantity squared in the second term. Exercise 1.4(b) with \(n = K\) gives \(\sum_k e_{jk}^2 - K\overline{e}_j^2 = \sum_k (e_{jk}-\overline{e}_j)^2\), and since \(e_{jk}-\overline{e}_j = (Y_{jk}-\overline{Y}_j) - \beta_j(x_{jk}-\overline{x}_j)\),

\begin{align*} \sum_k (e_{jk}-\overline{e}_j)^2 &= S_{YY} - 2\beta_j b_j S_{xx} + \beta_j^2 S_{xx}, \\ (b_j-\beta_j)^2 S_{xx} &= b_j^2 S_{xx} - 2b_j\beta_j S_{xx} + \beta_j^2 S_{xx}, \end{align*}

so subtracting cancels every \(\beta_j\) term and leaves \(S_{YY} - b_j^2 S_{xx} = \widehat{S}_{1j}\) by (i). Summing over \(j\) gives the stated identity, which is purely algebraic: it holds for any \(\alpha_j, \beta_j\), as a numerical check on the birthweight data with both the Table 2.5 values and arbitrary ones confirms:

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 though the three components differ by orders of magnitude.

(2) \(Y_{jk} \sim \mathrm{N}(\alpha_j+\beta_j x_{jk},\sigma^2)\) is the model assumption (2.7) itself; \(\overline{Y}_j\) and \(b_j\) are linear in independent Normals, hence Normal, so only the moments are wanted. Immediately

\begin{equation*} \mathrm{E}(\overline{Y}_j) = \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}, \end{equation*}

and for the slope, \(\sum_k (x_{jk}-\overline{x}_j) = 0\) lets \(Y_{jk}-\overline{Y}_j\) be replaced by \(Y_{jk}\) in \(b_j = S_{xY}/S_{xx} = S_{xx}^{-1}\sum_k (x_{jk}-\overline{x}_j)Y_{jk}\), whence, with \(\sum_k (x_{jk}-\overline{x}_j)x_{jk} = S_{xx}\),

\begin{align*} \mathrm{E}(b_j) &= \frac{1}{S_{xx}}\sum_k (x_{jk}-\overline{x}_j)(\alpha_j+\beta_j x_{jk}) = \beta_j, \\ \mathrm{var}(b_j) &= \frac{\sigma^2}{S_{xx}^2}\sum_k (x_{jk}-\overline{x}_j)^2 = \frac{\sigma^2}{\sum_k x_{jk}^2 - K\overline{x}_j^2}. \end{align*}

(3) The stated independence cannot be read literally, \(\overline{Y}_j\) being a function of the \(Y_{jk}\) of its own group; what holds, and what Section 2.2.2 needs, is that within each group the three terms of the decomposition are mutually independent, with independence across groups from disjointness. Rearranged, the identity reads

\begin{equation*} \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}. \end{equation*}

Fix \(j\), take \(\mathbf{e}_j \sim \mathrm{N}_K(\mathbf{0},\sigma^2 \mathbf{I})\), and let \(V_j \subset \mathbb{R}^K\) be spanned by the orthonormal \(\mathbf{u}_1 = \mathbf{1}/\sqrt{K}\) and \(\mathbf{u}_2 = (\mathbf{x}_j - \overline{x}_j\mathbf{1})/\sqrt{S_{xx}}\). Then

\begin{equation*} \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}}, \end{equation*}

so the three terms are \(\|(\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\) with \(\mathbf{P}_j\) the projection onto \(V_j\); \(\mathbf{e}_j\) being spherically Normal, its components along orthogonal directions and in the orthogonal complement are independent, so dividing by \(\sigma^2\),

\begin{equation*} \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), \end{equation*}

independently, adding to \(\chi^2(K)\) as they must. In particular \(\overline{Y}_j\) and \(b_j\) are jointly Normal with

\begin{equation*} \mathrm{cov}(\overline{Y}_j, b_j) = \frac{\sigma^2}{K S_{xx}}\sum_k (x_{jk}-\overline{x}_j) = 0, \end{equation*}

hence independent, which is the substantive content of the claim. Summing over \(j\) gives \(\widehat{S}_1/\sigma^2 \sim \chi^2(JK-2J)\) as quoted in Section 2.2.2.

The link is \(g(\mu) = e^{\mu}\): exponentiating \(\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 , \]

which is the form \(g[\mathrm{E}(Y_i)] = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}\) of Section 2.4 with \(x_{i1} = x_i\), \(x_{i2} = x_i^2\) – quadratic in \(x\) but linear in the parameters, which is all (2.13) requires. With \(N = 6\) and \(p = 3\),

\begin{equation*} \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}, \end{equation*}

and the model in the form of (2.13) is

\begin{equation*} g[\mathrm{E}(\mathbf{y})] = \begin{bmatrix} \exp[\mathrm{E}(Y_1)] \\ \vdots \\ \exp[\mathrm{E}(Y_6)] \end{bmatrix} = \mathbf{X}\boldsymbol{\beta}. \end{equation*}

The design matrix:

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

The logarithm restricts \(\boldsymbol{\beta}\) to the region where \(\beta_0+\beta_1 x+\beta_2 x^2 > 0\) over the range of \(x\). Fitting by least squares with a custom link:

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 and not a convergence failure: the fit is simply poor, the last observation \(y = 16.50\) lying far above what the model can reach.

With \(g\) the identity, the model in the form of (2.13) is \(\mathrm{E}(\mathbf{y}) = \mathbf{X}\boldsymbol{\beta}\) with

\begin{equation*} \mathbf{y} = \begin{bmatrix} Y_{11} \\ Y_{12} \\ Y_{13} \\ Y_{21} \\ Y_{22} \\ Y_{23} \end{bmatrix}, \qquad \boldsymbol{\beta} = \begin{bmatrix} \mu \\ \alpha_1 \\ \beta_1 \\ \beta_2 \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}. \end{equation*}

The unconstrained \(\mu+\alpha_j+\beta_k\) carries \(1+2+3 = 6\) parameters for \(N = 6\) observations and is not identifiable (Example 2.4.3(b)); the sum-to-zero constraints of Example 2.4.3(d) remove one from each factor, leaving \(p = 4\). Substituting \(\alpha_2 = -\alpha_1\) and \(\beta_3 = -\beta_1-\beta_2\),

\begin{align*} \mathrm{E}(Y_{11}) &= \mu+\alpha_1+\beta_1, & \mathrm{E}(Y_{12}) &= \mu+\alpha_1+\beta_2, & \mathrm{E}(Y_{13}) &= \mu+\alpha_1-\beta_1-\beta_2, \\ \mathrm{E}(Y_{21}) &= \mu-\alpha_1+\beta_1, & \mathrm{E}(Y_{22}) &= \mu-\alpha_1+\beta_2, & \mathrm{E}(Y_{23}) &= \mu-\alpha_1-\beta_1-\beta_2, \end{align*}

whose coefficients of \((\mu,\alpha_1,\beta_1,\beta_2)\) are the rows of \(\mathbf{X}\) above. This is contr.sum coding on 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, unlike the unconstrained six; \(N - p = 2\) degrees of freedom remain for \(\sigma^2\), which a \(j \times k\) interaction would consume entirely.

The design is balanced, the two factors 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 \(\mathbf{X}^{\mathrm{T}}\mathbf{X}\) makes \(\widehat{\mu}\), \(\widehat{\alpha}_1\) and the pair \((\widehat{\beta}_1,\widehat{\beta}_2)\) mutually uncorrelated, so one factor’s estimated effect is unchanged by adding or dropping the other; the off-diagonal \(2\) reflects only that the two columns coding a three-level factor are not orthogonal to each other under sum-to-zero coding.

Each part gives the three components of Section 3.4: response distribution, linear component, and link.

(a) \(Y_i\) is the weight (kg) of person \(i\): a continuous, roughly symmetric measurement with no hard upper bound, so \(Y_i \sim N(\mu_i,\sigma^2)\) (exponential family, Section 3.2.2) with the identity link, giving the normal linear model of Section 3.5.1. With \(x_{i1}\) age, \(x_{i2}\) sex (dummy), \(x_{i3}\) height, \(x_{i4}\) mean daily intake and \(x_{i5}\) mean daily expenditure,

\begin{equation*} g(\mu_i) = \mu_i = \beta_1 + \beta_2 x_{i1} + \beta_3 x_{i2} + \beta_4 x_{i3}

• \beta_5 x_{i4} + \beta_6 x_{i5}. \end{equation*}

(b) \(Y_i\) is the number infected out of \(n_i = 20\) at exposure level \(i\) (\(i=1,\dots,5\)), the mice at a level being independent binary trials with common probability, so \(Y_i \sim \mathrm{Bin}(20,\pi_i)\) (exponential family with natural parameter \(\log[\pi/(1-\pi)]\), Section 3.2.3). Modelling the count rather than the proportion keeps the variance \(n_i\pi_i(1-\pi_i)\). With \(x_i\) the exposure level (or log dose) and the natural logit link,

\[ g(\pi_i)=\log\!\left(\frac{\pi_i}{1-\pi_i}\right) = \beta_1 + \beta_2 x_i , \]

a logistic dose-response model, which keeps \(\pi_i \in (0,1)\) for all \(\beta_1, \beta_2, x_i\).

(c) \(Y_i\) is the number of trips made by household \(i\) in a week: an unbounded count of events in a fixed period, so \(Y_i \sim \mathrm{Po}(\mu_i)\) (Section 3.2.1). With \(x_{i1}\) household size, \(x_{i2}\) income and \(x_{i3}\) distance, the log link keeps \(\mu_i>0\) and makes effects multiplicative on the rate:

\[ g(\mu_i)=\log \mu_i = \beta_1 + \beta_2 x_{i1} + \beta_3 x_{i2} + \beta_4 x_{i3}, \]

with \(\beta_2 > 0\) and \(\beta_4 < 0\) expected. The Poisson assumption \(\mathrm{E}(Y)=\mathrm{var}(Y)\) should be checked, a negative binomial (Chapter 9) being preferred if the counts are overdispersed.

The natural parameter is \(b(\beta) = -\beta\), with \(\mathrm{E}(Y)=\alpha/\beta\) and \(\mathrm{var}(Y)=\alpha/\beta^2\). Taking logarithms, for \(y>0\),

\begin{equation*} f(y;\beta) = \exp\left[ \alpha\log\beta - \log\Gamma(\alpha)

• (\alpha-1)\log y - y\beta \right], \end{equation*}

which is the form (3.3), \(f=\exp[a(y)b(\theta)+c(\theta)+d(y)]\), with

\begin{equation*} a(y)=y, \quad b(\beta) = -\beta, \quad c(\beta) = \alpha\log\beta - \log\Gamma(\alpha), \quad d(y) = (\alpha-1)\log y , \end{equation*}

so the family is exponential and, as \(a(y)=y\), in canonical form; the known shape \(\alpha\) enters only through \(c\) and \(d\), as Section 3.2 permits. Then \(b’ = -1\), \(b’’ = 0\), \(c’ = \alpha/\beta\), \(c’’ = -\alpha/\beta^2\), so by (3.9) and (3.12),

\begin{align*} \mathrm{E}(Y) &= -\frac{c’(\beta)}{b’(\beta)} = -\frac{\alpha/\beta}{-1} = \frac{\alpha}{\beta}, \\ \mathrm{var}(Y) &= \frac{b’‘c’ - c’‘b’}{(b’)^3} = \frac{0 - (-\alpha/\beta^2)(-1)}{(-1)^3} = \frac{\alpha}{\beta^2}. \end{align*}

A numerical check against R’s 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 each case take logarithms and read off (3.3), \(f(y;\theta) = \exp[a(y)b(\theta) + c(\theta) + d(y)]\); the form is canonical only if \(a(y)=y\) (Section 3.2).

(a) Pareto, \(y>1\), \(\theta>0\):

\begin{equation*} f(y;\theta) = \theta y^{-\theta-1} = \exp\left[ (\log y)(-\theta) + \log\theta - \log y \right], \end{equation*}

so (3.3) holds with \(a(y) = \log y\), \(b(\theta) = -\theta\), \(c(\theta) = \log\theta\), \(d(y) = -\log y\). Exponential family, but not canonical, since \(a(y) = \log y \neq y\); equivalently \(\log Y \sim \mathrm{Exp}(\theta)\), which is used in Exercises 3.9 and 3.11.

(b) Exponential, \(y>0\):

\[ f(y;\theta) = \theta e^{-y\theta} = \exp\left[ -y\theta + \log\theta \right], \]

so \(a(y)=y\), \(b(\theta)=-\theta\), \(c(\theta)=\log\theta\), \(d(y)=0\): canonical, with natural parameter \(-\theta\). It is the \(\alpha=1\) case of Exercise 3.2’s Gamma.

(c) Negative binomial, \(y=0,1,2,\dots\) with \(r\) known, so the binomial coefficient depends on \(y\) alone:

\begin{equation*} f(y;\theta) = \exp\left[ y\log(1-\theta) + r\log\theta

• \log\binom{y+r-1}{r-1} \right], \end{equation*}

giving \(a(y)=y\), \(b(\theta)=\log(1-\theta)\), \(c(\theta)=r\log\theta\), \(d(y)=\log\binom{y+r-1}{r-1}\): canonical, with natural parameter \(\log(1-\theta)\). Knowing \(r\) is essential, since otherwise \(d(y)\) would depend on a second parameter and (3.3) would fail.

All three are in canonical form, so \(a(Y)=Y\) and (3.9), \(\mathrm{E}[a(Y)] = -c’/b’\), and (3.12), \(\mathrm{var}[a(Y)] = (b’‘c’-c’‘b’)/(b’)^3\), give the moments of \(Y\) directly from the \(b, c\) of Table 3.1.

(a) Poisson (Section 3.2.1): \(b(\theta)=\log\theta\), \(c(\theta)=-\theta\), so \(b’=1/\theta\), \(b’’=-1/\theta^2\), \(c’=-1\), \(c’’=0\) and

\begin{equation*} \mathrm{E}(Y) = -\frac{-1}{1/\theta} = \theta, \qquad \mathrm{var}(Y) = \frac{(-1/\theta^2)(-1)}{(1/\theta)^3} = \theta. \end{equation*}

(b) Normal (Section 3.2.2), \(\sigma^2\) a known nuisance parameter: \(b(\mu)=\mu/\sigma^2\), \(c(\mu) = -\mu^2/(2\sigma^2) - \tfrac12\log(2\pi\sigma^2)\), so \(b’=1/\sigma^2\), \(b’’=0\), \(c’=-\mu/\sigma^2\), \(c’’=-1/\sigma^2\) and

\begin{equation*} \mathrm{E}(Y) = -\frac{-\mu/\sigma^2}{1/\sigma^2} = \mu, \qquad \mathrm{var}(Y) = \frac{-(-1/\sigma^2)(1/\sigma^2)}{(1/\sigma^2)^3} = \sigma^2. \end{equation*}

(c) Binomial (Section 3.2.3): \(b(\pi)=\log\pi - \log(1-\pi)\), \(c(\pi)=n\log(1-\pi)\), so

\begin{equation*} b’(\pi) = \frac{1}{\pi(1-\pi)}, \quad b’’(\pi) = -\frac{1-2\pi}{\pi^2(1-\pi)^2}, \quad c’(\pi) = -\frac{n}{1-\pi}, \quad c’’(\pi) = -\frac{n}{(1-\pi)^2}, \end{equation*}

whence \(\mathrm{E}(Y) = [n/(1-\pi)]\pi(1-\pi) = n\pi\) and, with denominator \((b’)^3 = [\pi(1-\pi)]^{-3}\),

\begin{align*} \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). \end{align*}

Applying (3.9) and (3.12) numerically to the same \(b, c\) by central differences:

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 discrepancy in the Poisson and Normal variances is truncation error in the numerical second derivative.

(a) \(\mathrm{E}(Y) = r(1-\theta)/\theta\) and \(\mathrm{var}(Y) = r(1-\theta)/\theta^2 = \mathrm{E}(Y)/\theta\), so \(\mathrm{var}(Y) > \mathrm{E}(Y)\) for \(0<\theta<1\). From Exercise 3.3(c) the distribution is canonical with \(b(\theta)=\log(1-\theta)\), \(c(\theta)=r\log\theta\), giving

\begin{equation*} b’= -\frac{1}{1-\theta}, \quad b’’= -\frac{1}{(1-\theta)^2}, \quad c’= \frac{r}{\theta}, \quad c’’= -\frac{r}{\theta^2}, \end{equation*}

so by (3.9) \(\mathrm{E}(Y) = -(r/\theta)\big/(-1/(1-\theta)) = r(1-\theta)/\theta\), while (3.12) has numerator \(-r/[\theta(1-\theta)^2] - r/[\theta^2(1-\theta)]\) over \((b’)^3 = -1/(1-\theta)^3\), whence

\begin{align*} \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}. \end{align*}

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

(b) Choosing the distribution is choosing the variance function, since a generalized linear model has no free \(\sigma^2\) to absorb misspecification. Writing \(\mu = \mathrm{E}(Y)\),

\begin{equation*} \text{Binomial: } \mathrm{var}(Y) = \mu(1-\pi) < \mu, \quad \text{Poisson: } \mathrm{var}(Y) = \mu, \quad \text{NBin: } \mathrm{var}(Y) = \mu/\theta > \mu. \end{equation*}

So: follow the sampling mechanism where there is one – a count out of \(n\) independent binary trials is Binomial, its underdispersion being a feature of the design. For unbounded counts the Poisson is the default, but \(\mathrm{E}(Y)=\mathrm{var}(Y)\) is testable, e.g. by comparing residual deviance or Pearson \(X^2\) with the residual degrees of freedom. Counts are frequently overdispersed (Section 3.2.1) through unmeasured heterogeneity, clustering or excess zeros, and the negative binomial is the remedy, keeping the log link and the multiplicative interpretation of \(\boldsymbol\beta\) while estimating the extra parameter (Chapter 9). Ignoring overdispersion damages not the point estimates but the standard errors, understating \(\mathrm{var}(\hat{\boldsymbol\beta})\) and so giving intervals that are too narrow.

No: the model of Example 3.5.3 is inadequate, because \(\log n_i\) enters as an offset with no intercept, forcing the fitted log rate through the origin at \(i=0\).

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

The poisson family object is masked by a dobson data set of the same name, hence 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 residual deviance is \(6990.4\) on \(7\) degrees of freedom, against a \(\chi^2(7)\) mean of \(7\), and \(\hat\theta = -2.72 < 0\) makes the fitted rate \(e^{\hat\theta i}\) decrease with age, contradicting Figure 3.2. Without an intercept no single \(\theta\) can make \(e^{\theta i}\) both start near \(e^{\theta}\) and pass through rates running from \(5.6\times10^{-5}\) to \(6.5\times10^{-3}\); the exponential shape is not at fault.

Following the hint, regress the log rate on the age-group index by least squares, Figure 3.2 being approximately linear in \(i\) (equally spaced indices being 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 line explains \(94\%\) of the variation, with \(\hat\beta_2 = 0.665\) multiplying the death rate by \(e^{0.665}=1.95\) per five-year band. But least squares assumes constant variance on the log-rate scale, whereas \(\mathrm{var}[\log(y_i/n_i)] \approx 1/y_i\) ranges from \(1\) to \(0.015\) across the groups. The correct version keeps the offset and adds an intercept, \[ \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 smaller than the OLS \(0.665\) because the poorly determined early groups no longer dominate the weighting; the rate ratio per five-year band is \(e^{0.522}=1.69\), \(95\%\) interval \(\exp(0.522 \pm 1.96\times 0.039) = (1.56, 1.82)\).

Comparing observed and expected deaths group by group, as the hint asks:

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 Example 3.5.3 model predicts \(1167\) deaths in the youngest group, where one occurred, and essentially zero in the oldest, where \(65\) occurred. The OLS expectations are better but poor at the extremes (\(107\) against \(65\) observed in the oldest group) and sum to \(230.7\) against \(205\) observed, whereas the Poisson fit with an intercept reproduces the total exactly, its intercept score equation forcing \(\sum \hat\mu_i = \sum y_i\). That fit has deviance \(14.69\) on \(6\) degrees of freedom, \(p = 0.023\), and Pearson residuals negative at both ends and positive in the middle – the signature of curvature in \(\log(\text{rate})\) against \(i\), which a quadratic term removes.

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)

So the Example 3.5.3 model as written is inadequate – off by three orders of magnitude, with a slope of the wrong sign – though its exponential shape is right: with an intercept, \(\log\mu_i = \log n_i + \beta_1 + \beta_2 i\) describes the data reasonably, the coronary death rate rising by about \(1.69\) per five-year band. The quadratic model, deviance \(3.09\) on \(5\) df (\(p=0.69\)) and AIC lower by nearly \(10\), is preferred.

(a) Taking logarithms,

\begin{equation*} f(y_i;\pi_i) = \pi_i^{y_i}(1-\pi_i)^{1-y_i} = \exp\left[ y_i\log\left(\frac{\pi_i}{1-\pi_i}\right) + \log(1-\pi_i) \right], \end{equation*}

which is (3.3) with \(a(y_i)=y_i\), \(b(\pi_i)=\log[\pi_i/(1-\pi_i)]\), \(c(\pi_i)=\log(1-\pi_i)\), \(d(y_i)=0\): the Bernoulli is in the exponential family, being the \(n=1\) case of Section 3.2.3.

(b) \(a(y_i)=y_i\), so the form is canonical and the natural parameter is \(b(\pi_i) = \log[\pi_i/(1-\pi_i)]\), the logit – whence the logit is the natural link for binary data, the linear predictor then being the natural parameter itself.

(c) \(\mathrm{E}(Y_i)=1\cdot\pi_i + 0\cdot(1-\pi_i)=\pi_i\); by (3.9), with \(b’(\pi_i) = [\pi_i(1-\pi_i)]^{-1}\) and \(c’(\pi_i) = -(1-\pi_i)^{-1}\),

\begin{equation*} \mathrm{E}(Y_i) = -\frac{c’(\pi_i)}{b’(\pi_i)} = \frac{1/(1-\pi_i)}{1/[\pi_i(1-\pi_i)]} = \pi_i , \end{equation*}

which is Exercise 3.4(c) with \(n=1\).

(d) Exponentiating \(\log[\pi/(1-\pi)] = \mathbf{x}^T\boldsymbol\beta\) gives \(\pi = (1-\pi)e^{\mathbf{x}^T\boldsymbol\beta}\), so

\begin{equation*} \pi = \frac{e^{\mathbf{x}^T\boldsymbol\beta}}{1+e^{\mathbf{x}^T\boldsymbol\beta}} = \frac{1}{1+e^{-\mathbf{x}^T\boldsymbol\beta}} , \end{equation*}

the logit being a strictly increasing bijection \((0,1)\to\mathbb{R}\) with this inverse; whatever real value \(\mathbf{x}^T\boldsymbol\beta\) takes, \(\pi\) lies in \((0,1)\).

(e) The graph of \(\pi(x) = e^{\beta_1+\beta_2 x}/(1+e^{\beta_1+\beta_2 x})\) is the sigmoid logistic curve. For \(\beta_2>0\) it is strictly increasing, \(d\pi/dx = \beta_2\pi(1-\pi) > 0\), with asymptotes \(\pi\to 0\) and \(\pi\to 1\), passing through \(\pi = 1/2\) at \(x = -\beta_1/\beta_2\), symmetric about that point, where the slope is steepest at \(\beta_2/4\). For \(\beta_2<0\) it decreases; \(\beta_2=0\) gives the horizontal line \(\pi = e^{\beta_1}/(1+e^{\beta_1})\); larger \(|\beta_2|\) sharpens the transition.

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")

With \(x\) an insecticide dose and \(\pi\) the probability of death, this is the classical dose-response curve: almost no insects die at low doses and almost all at high ones, neither extreme reached exactly, with a steep rise over a narrow middle band beyond which extra dose buys little. The dose \(x = -\beta_1/\beta_2\) at which \(\pi = 1/2\) is the median lethal dose LD50, the standard measure of potency (\(5\) in the left panel, with maximum slope \(\beta_2/4 = 0.25\) per unit dose). The slope \(\beta_2\) multiplies the odds of death by \(e^{\beta_2}\) per unit dose, constantly across the range, whereas the effect on the probability is largest near LD50.

Yes, though not in canonical form. Taking logarithms with \(\phi\) known,

\begin{equation*} \log f(y;\theta) = -\log\phi + \frac{y}{\phi} - \frac{\theta}{\phi}

• e^{y/\phi}e^{-\theta/\phi}, \end{equation*}

whose last term is the only one mixing \(y\) and \(\theta\) and already factorises, so (3.3) holds with

\begin{equation*} a(y) = e^{y/\phi}, \quad b(\theta) = -e^{-\theta/\phi}, \quad c(\theta) = -\frac{\theta}{\phi} - \log\phi, \quad d(y) = \frac{y}{\phi}. \end{equation*}

The support is all of \(\mathbb{R}\) and free of \(\theta\), as Section 3.3 requires; and \(\phi\) must genuinely be known, since otherwise \(a\) and \(d\) depend on it and the family is two-parameter. Since \(a(y) = e^{y/\phi} \neq y\) the form is not canonical – it is \(e^{Y/\phi}\) that is canonical, the display showing \(e^{(Y-\theta)/\phi} \sim \mathrm{Exp}(1)\) – so (3.9) and (3.12) give the moments of \(a(Y)\), not of \(Y\). With \(b’ = \phi^{-1}e^{-\theta/\phi}\), \(b’’ = -\phi^{-2}e^{-\theta/\phi}\), \(c’ = -1/\phi\), \(c’’=0\),

\begin{align*} \mathrm{E}\left(e^{Y/\phi}\right) &= -\frac{c’}{b’} = \frac{1/\phi}{e^{-\theta/\phi}/\phi} = e^{\theta/\phi}, \\ \mathrm{var}\left(e^{Y/\phi}\right) &= \frac{b’‘c’ - c’‘b’}{(b’)^3} = \frac{\left(-\phi^{-2}e^{-\theta/\phi}\right)\left(-\phi^{-1}\right)}{\phi^{-3}e^{-3\theta/\phi}} = e^{2\theta/\phi}. \end{align*}

Confirmed by simulation, generating via \(e^{(Y-\theta)/\phi}\sim\mathrm{Exp}(1)\):

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

Yes, with link \(g(\mu)=\mu^{1/2}\) and linear predictor \(\beta_0+\beta_1x_i\), checking the three components of Section 3.4.

(i) Distribution. The \(Y_i\) are independent Pareto, which is in the exponential family with \(a(y)=\log y\), \(b(\theta_i)=-\theta_i\), \(c(\theta_i)=\log\theta_i\), \(d(y)=-\log y\) (Exercise 3.3(a)), each depending on one parameter. Strictly, Section 3.4 asks for the canonical form \(\exp[y_ib(\theta_i)+c(\theta_i)+d(y_i)]\), which the Pareto misses since \(a(y)=\log y\); but \(a\) is a fixed known transformation and \(\log Y_i \sim \mathrm{Exp}(\theta_i)\) is canonical, so the model is a generalized linear model for the transformed response.

(ii) Linear component. \(\beta_0+\beta_1x_i = \mathbf{x}_i^T\boldsymbol\beta\) with \(\mathbf{x}_i^T = [1\;\;x_i]\); linearity in \(\boldsymbol\beta\), not in \(x_i\), is what is required.

(iii) Link. \(\mu_i = (\beta_0+\beta_1x_i)^2\) inverts to \(g(\mu) = \mu^{1/2}\), strictly increasing and differentiable on \(\mu>0\).

Two restrictions on the parameter space follow. The square root inverts \(\mu \mapsto \mu^2\) on one branch only, so \((\beta_0,\beta_1)\) and \((-\beta_0,-\beta_1)\) give identical means and identifiability needs \(\beta_0+\beta_1x_i>0\) throughout. And \(\mathrm{E}(Y)=\theta/(\theta-1) > 1\) for every admissible \(\theta>1\) (the mean not existing for \(\theta\le 1\)), so coherence needs \(\beta_0+\beta_1x_i>1\) at every observed \(x_i\).

No: it is a nonlinear regression model. The distribution is fine – independent Normals, canonical exponential family (Section 3.2.2) – but the systematic part is not \(g(\mu_i)=\mathbf{x}_i^T\boldsymbol\beta\) for any link. Rearranging,

\[ e^{\mu_i - \beta_0} = \beta_1 + \beta_2 x_i , \]

whose right side is linear in \((\beta_1,\beta_2)\); so if \(\beta_0\) were known this would be a generalized linear model with link \(g(\mu)=e^{\mu-\beta_0}\). But \(\beta_0\) must be estimated, and a link is a fixed known function of \(\mu\) alone. Equivalently, \(\beta_0\) enters outside the logarithm and \(\beta_1, \beta_2\) inside it, so no \(g\) puts all three into one linear form. (The three are in any case unidentified without a constraint, since adding \(\log k\) to \(\beta_0\) and dividing \(\beta_1, \beta_2\) by \(k\) leaves \(\mu_i\) unchanged; imposing \(\beta_0=0\) does make it a generalized linear model with \(g(\mu)=e^{\mu}\).) As written the model needs nonlinear least squares, and \(\beta_1+\beta_2x_i>0\) at every observation.

\(U = 1/\theta - \log Y\) and \(\mathfrak{I} = 1/\theta^2\). From Exercise 3.3(a) the Pareto has \(a(y)=\log y\), \(b(\theta)=-\theta\), \(c(\theta)=\log\theta\), \(d(y)=-\log y\), so \(l(\theta;y) = -\theta\log y + \log\theta - \log y\) and, by (3.13) with \(b’=-1\), \(c’=1/\theta\),

\[ U = \frac{dl}{d\theta} = \frac{1}{\theta} - \log Y . \]

By (3.9), \(\mathrm{E}(\log Y) = -c’/b’ = 1/\theta\) (directly, \(\log Y\sim\mathrm{Exp}(\theta)\)), so \(\mathrm{E}(U) = 1/\theta - 1/\theta = 0\), the general result (3.14). For the information, Section 3.3 gives \(\mathfrak{I} = (b’)^2\mathrm{var}[a(Y)] = \mathrm{var}(\log Y) = 1/\theta^2\), the same as (3.15) with \(b’’=0\) and \(c’’=-1/\theta^2\),

\begin{equation*} \mathrm{var}(U) = \frac{b’’(\theta)c’(\theta)}{b’(\theta)} - c’’(\theta) = 0 + \frac{1}{\theta^2}, \end{equation*}

and as (3.16), \(U’ = -1/\theta^2\) being non-random. For \(N\) independent observations the log-likelihoods add, so \(U = N/\theta - \sum_i \log Y_i\) and \(\mathfrak{I} = N/\theta^2\), giving \(\hat\theta = N/\sum\log Y_i\) with asymptotic standard error \(\theta/\sqrt{N}\).

Simulation 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

(a) \(\mathrm{Exp}(\theta) = G(1,\theta)\): putting \(\alpha=1\), \(\beta=\theta\) in the Gamma density of Exercise 3.2, with \(\Gamma(1)=1\),

\[ f(y;1,\theta) = \frac{\theta^{1}}{\Gamma(1)}y^{0}e^{-y\theta} = \theta e^{-y\theta}, \]

the exponential density of Exercise 3.3(b); the moments agree, \(\alpha/\beta = 1/\theta\) and \(\alpha/\beta^2 = 1/\theta^2\). (Figure 3.3’s label \(\alpha=\theta\), \(\beta=1\) interchanges the substitutions: \(G(\theta,1)\) is exponential only at \(\theta=1\).)

(b) For \(X\sim U[0,1]\) and \(Y=-\theta\log X\), we have \(Y>0\) and, for \(y>0\),

\[ P(Y>y) = P\!\left(X < e^{-y/\theta}\right) = e^{-y/\theta}, \]

since \(P(X<x)=x\) on \((0,1)\); differentiating, \(f_Y(y) = \theta^{-1}e^{-y/\theta}\), exponential with mean \(\theta\). Under the book’s rate parametrisation \(f(y;\theta)=\theta e^{-y\theta}\) the printed statement reads \(-\theta\log X \sim \mathrm{Exp}(1/\theta)\); the transformation giving \(\mathrm{Exp}(\theta)\) exactly is \(Y = -\theta^{-1}\log X\).

(c)(i) For \(Y\sim\mathrm{Bin}(n,\pi)\), \(M_Y(t) = (1-\pi+\pi e^t)^n\), so with \(\lambda = n\pi\) fixed and \(n\to\infty\),

\begin{equation*} M_Y(t) = \left[1 + \frac{\lambda(e^t-1)}{n}\right]^n \longrightarrow \exp\left[\lambda\left(e^t-1\right)\right], \end{equation*}

which is the moment generating function of \(\mathrm{Po}(\lambda)\); the continuity theorem gives \(\mathrm{Bin}(n,\lambda/n)\to\mathrm{Po}(\lambda)\).

(ii) For \(Y\sim\mathrm{NBin}(r,\theta)\), the negative binomial series \(\sum_y \binom{y+r-1}{r-1}u^y = (1-u)^{-r}\) gives \(M_Y(t) = [\theta/(1-(1-\theta)e^t)]^r\) for \((1-\theta)e^t<1\), so with \(\lambda = r(1-\theta)\) fixed and \(r\to\infty\),

\begin{align*} \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[-\frac{\lambda}{r} + \frac{\lambda e^t}{r} + O(r^{-2})\right] \longrightarrow \lambda(e^t-1), \end{align*}

so \(\mathrm{NBin}(r,\theta)\to\mathrm{Po}(r(1-\theta))\) – consistent with Exercise 3.5, where \(\mathrm{var}(Y)/\mathrm{E}(Y) = 1/\theta \to 1\) as \(\theta\to1\).

Both limits checked by the largest absolute discrepancy between probability functions, \(\lambda=3\):

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

(d) Each family is closed under convolution, so each distribution is a sum of independent identically distributed pieces and the Central Limit Theorem applies.

(i) Poisson. For integer \(\mu\), \(Y\sim\mathrm{Po}(\mu)\) is the sum of \(\mu\) independent \(\mathrm{Po}(1)\) variables, each of mean and variance \(1\) (Exercise 3.4(a)), so \(Y\) is asymptotically \(N(\mu,\mu)\); for non-integer \(\mu\), write \(Y\) as a sum of \(n\) independent \(\mathrm{Po}(\mu/n)\) variables. Figure 3.3’s rule of thumb is \(\mu>15\).

(ii) Binomial. \(Y\sim\mathrm{Bin}(n,\pi)\) is a sum of \(n\) independent Bernoulli variables with mean \(\pi\) and variance \(\pi(1-\pi)\) (Exercise 3.7), so \(Y\) is asymptotically \(N(n\pi, n\pi(1-\pi))\). The summands are skew unless \(\pi=1/2\), the standardised skewness of \(Y\) being \((1-2\pi)/\sqrt{n\pi(1-\pi)}\), small only when \(n\pi(1-\pi)\) is large; for very small \(\pi\) the boundary at \(0\) truncates the Normal shape and the Poisson limit of (c)(i) applies instead. Hence \(n\pi>5\) and \(n\pi(1-\pi)>5\).

(iii) Gamma. For integer \(\alpha\), \(Y\sim G(\alpha,\beta)\) is a sum of \(\alpha\) independent \(\mathrm{Exp}(\beta)\) variables by (a), each of mean \(1/\beta\) and variance \(1/\beta^2\), so \(Y\) is asymptotically \(N(\alpha/\beta, \alpha/\beta^2)\), matching Exercise 3.2; for real \(\alpha>0\) use \(G(\alpha,\beta) = \sum_{j=1}^n G(\alpha/n,\beta)\) independently.

Checked by the largest absolute difference between exact and Normal distribution functions, with a continuity correction in the 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 parameters are \(\mu = 5, 20, 100\); \(n = 10, 50, 500\) with \(\pi=0.3\); and \(\alpha = 2, 10, 100\) with \(\beta=1\). Every row falls roughly as the square root of the parameter, the Berry-Esseen rate, with maximum error at or below about \(1.5\%\) by \(\mu=20\), \(n=50\), \(\alpha=100\) – which is where Figure 3.3’s rules of thumb sit.

The fitted model is \(\widehat\lambda_i = e^{0.996}i^{1.327}\), obtained identically by hand-coded iterative weighted least squares and by glm.

(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 wobbles at \(i = 7,8\) and \(i = 15\), and not linearly: the increment per quarter is about \(5\) cases early and \(15\) at the end. The scatter widens with the level, the mean-variance behaviour \(\operatorname{var}(Y_i) = \operatorname{E}(Y_i)\) that motivates the Poisson model.

(b) The log-log plot. Under \(\lambda_i = i^{\theta}\), \(\log \lambda_i = \theta \log i\), so \(\log y_i\) against \(\log i\) should be 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 line, so the power law is a reasonable description; least squares through the log-log points gives slope \(1.454\) and intercept \(0.686\). The intercept is not zero, which is why (c) generalises to \(\log \lambda_i = \beta_1 + \beta_2 x_i\), \(x_i = \log i\), of which \(\lambda_i = i^{\theta}\) is the case \(\beta_1 = 0\).

(c) Iterative weighted least squares from first principles. With \(Y_i \sim \operatorname{Po}(\mu_i)\) independent and \(\eta_i = \log \mu_i = \mathbf{x}_i^T \boldsymbol\beta\),

\begin{equation*} \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 . \end{equation*}

From equation (4.23) the diagonal weights are

\begin{equation*} 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 , \end{equation*}

so \(\mathbf{W} = \operatorname{diag}(\mu_1, \ldots, \mu_N)\), the weights reducing to the variance function because the log link is canonical for the Poisson. From equation (4.24) the working response is

\begin{equation*} 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} , \end{equation*}

with \(\eta_i, \mu_i\) at \(\mathbf{b}^{(m-1)}\). With \(\mathbf{X}\) having rows \((1, \log i)\), the information matrix \(\mathfrak{I} = \mathbf{X}^T \mathbf{W} \mathbf{X}\) has entries

\begin{equation*} \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 \left[ \mu_i \eta_i + (y_i - \mu_i) \right] \\ \sum_i x_i \left[ \mu_i \eta_i + (y_i - \mu_i) \right] \end{bmatrix} . \end{equation*}

Solving \(\mathbf{X}^T \mathbf{W} \mathbf{X} \mathbf{b}^{(m)} = \mathbf{X}^T \mathbf{W} \mathbf{z}\) repeatedly, from the log-log fit of (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 complete after three iterations, as in Table 4.4, giving \(\widehat\beta_1 = 0.99600\) and \(\widehat\beta_2 = 1.32661\). The inverse information matrix at the estimates gives the standard errors, by (4.12) in vector form:

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

(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

(The explicit poisson(link = "log") is needed because after library(dobson) the bare name poisson is the package’s Table 4.3 data set.)

Estimates and standard errors agree with (c) to every printed digit, \(b_1 = 0.99600\,(0.16971)\) and \(b_2 = 1.32661\,(0.06463)\), since glm runs the algorithm of (4.25) and Fisher scoring coincides with Newton-Raphson for the canonical link.

The fitted mean \(\widehat\lambda_i = 2.71\,i^{1.327}\) grows as a power of time with exponent about \(1.33\), not exponentially in \(i\); the approximate \(95\%\) interval \(1.3266 \pm 1.96 \times 0.0646 = (1.200, 1.453)\) excludes \(1\), so growth is faster than linear but far short of doubling behaviour, and the one-parameter model \(\beta_1 = 0\) of (b) is rejected (\(z = 5.87\)). The residual deviance \(21.755\) on \(18\) degrees of freedom gives \(p = 0.243\): no evidence of lack of fit or overdispersion. Against the literal exponential-growth model \(\log \lambda_i = \beta_1 + \beta_2 i\):

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\) against \(53.0\) on the same \(18\) degrees of freedom, the exponential model having \(p = 2.6 \times 10^{-5}\) for lack of fit and AIC \(31\) units worse. 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 or funnelling, running from \(-2.06\) (at \(i = 8\), where the series briefly flattens) to \(1.73\). The model is adequate.

The fitted exponential model is \(\widehat y_i = \exp(8.478 - 1.109x_i)\), which is adequate. Each \(y_i\) is a leukemia patient’s time to death in weeks; the \(x_i\) are \(\log_{10}\) initial white blood cell counts.

(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 trend is clearly negative: survival falls as the initial white cell count rises. The spread of \(y\) also shrinks with \(x\) – at \(x \approx 3\) survival ranges from \(16\) to \(156\) weeks, while at \(x = 5\) the three patients survived \(1\), \(1\) and \(65\) weeks – so variability scales with level, ruling out constant-variance least squares; and the decline is curved, consistent with a multiplicative effect of \(x\).

(b) The log link, \(g(\mu) = \log \mu\), since

\begin{equation*} \operatorname{E}(Y_i) = \exp(\beta_1 + \beta_2 x_i) \iff \log \operatorname{E}(Y_i) = \beta_1 + \beta_2 x_i , \end{equation*}

which is the form \(g(\mu_i) = \mathbf{x}_i^T \boldsymbol\beta\) of equation (4.16). It keeps \(\operatorname{E}(Y) > 0\) for every \(\beta\) and \(x\), accommodates the mean-dependent variance, and makes a unit increase in \(x_i\) multiply \(\operatorname{E}(Y)\) by \(e^{\beta_2}\); the trend in (a) implies \(\beta_2 < 0\).

(c) For \(f(y; \theta) = \theta e^{-y\theta}\), \(y>0\), integrating by parts with \(u=y\),

\begin{equation*} \operatorname{E}(Y) = \left[ -y e^{-y\theta} \right]_0^{\infty}

• \int_0^{\infty} e^{-y\theta}\,\mathrm{d}y = 0 + \left[ -\theta^{-1} e^{-y\theta} \right]_0^{\infty} = \frac{1}{\theta}, \end{equation*}

while the Gamma integral \(\int_0^{\infty} y^n e^{-y\theta}\,\mathrm{d}y = n!/\theta^{n+1}\) gives \(\operatorname{E}(Y^{2}) = \theta\cdot 2!/\theta^3 = 2/\theta^2\), so

\begin{equation*} \operatorname{var}(Y) = \frac{2}{\theta^{2}} - \frac{1}{\theta^{2}} = \frac{1}{\theta^{2}} = \left[ \operatorname{E}(Y) \right]^{2} , \end{equation*}

the Gamma variance \(\phi[\operatorname{E}(Y)]^2\) with \(\phi=1\): the standard deviation equals the mean, which (e) exploits.

(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

Gamma(link = "log") is used because the exponential is the Gamma with \(\phi = 1\). The estimates of \(\boldsymbol\beta\) are the exponential-model ones, \(\phi\) cancelling from the score equations (4.18) as a constant multiplier of \(\operatorname{var}(Y_i)\); only the standard errors are affected, and glm estimates \(\phi\) by default. Imposing \(\phi = 1\):

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 little (\(1.603\) to \(1.655\), \(0.387\) to \(0.400\)) because \(\widehat\phi = 0.939\) is already near \(1\). Here \(\widehat\beta_2 = -1.109\), so a tenfold increase in white cell count multiplies expected survival by \(e^{-1.109} = 0.33\); the \(95\%\) interval \(-1.109 \pm 1.96 \times 0.400 = (-1.893, -0.326)\) gives a factor between \(0.15\) and \(0.72\), the direction clear but the size imprecise with \(17\) patients. Fitted mean survival runs from \(197\) weeks at \(x = 2.88\) to \(19\) weeks at \(x = 5.00\).

(e) By (c), \(\operatorname{sd}(Y_i) = \operatorname{E}(Y_i)\) under the exponential model, so dividing by \(\widehat y_i\) divides by the estimated standard deviation and

\begin{equation*} r_i = \frac{y_i - \widehat y_i}{\widehat y_i}, \qquad \widehat y_i = \exp(\widehat\beta_1 + \widehat\beta_2 x_i) \end{equation*}

is genuinely standardized – these are 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

The model is adequate. The residuals are centred near zero with \(\operatorname{sd}(r_i) = 0.938\), and since the model implies \(\operatorname{sd}(Y_i) = \operatorname{E}(Y_i)\) a typical magnitude near \(1\) is expected; all but one fall in \((-0.95, 1.13)\), and the right skew is a property of \(r_i\) being bounded below by \(-1\) and unbounded above. Being Pearson residuals, \(\widehat\phi = \frac{1}{15}\sum r_i^{2} = 0.939\) matches part (d) and supports \(\phi = 1\). The residual deviance \(D = 19.457\) on \(15\) degrees of freedom is already the scaled deviance when \(\phi=1\), giving \(p = 0.194\) against \(\chi^2(15)\): no lack of fit. Neither plot shows funnelling or curvature, so the log link and mean function are appropriate. The one outlier, \(r_{17} = 2.47\) at \(x = 5.00\), \(y = 65\), reflects three patients sharing the highest count with survivals \(1\), \(1\) and \(65\) against a common \(\widehat y = 18.75\); since \(\Pr(Y > 3.47\mu) = e^{-3.47} = 0.031\), a residual this large among \(17\) is unusual but not extraordinary.

\(\widehat\beta = \exp(\overline{Y})\), with \(\overline{Y} = N^{-1}\sum_i Y_i\). From first principles:

\begin{align*} l(\beta; \mathbf{y}) &= -\frac{1}{2\sigma^{2}} \sum_{i=1}^{N} (y_i - \log \beta)^{2}

• N \log\left(\sigma\sqrt{2\pi}\right), \\ U = \frac{\mathrm{d}l}{\mathrm{d}\beta} &= \frac{1}{\beta\sigma^{2}} \sum_{i=1}^{N} (Y_i - \log\beta) = \frac{N}{\beta\sigma^{2}}\left(\overline{Y} - \log\beta\right), \\ U^{\prime} &= -\frac{N}{\beta^{2}\sigma^{2}}\left[1 + \overline{Y} - \log\beta\right]. \end{align*}

Since \(\beta > 0\) the factor \(N/(\beta\sigma^{2})\) never vanishes, so \(U = 0\) forces \(\log\beta = \overline{y}\); and \(U^{\prime}(\widehat\beta) = -N/(\widehat\beta^{2}\sigma^{2}) < 0\), a maximum. (Equivalently, by invariance under reparametrisation: \(\mu = \log\beta\) has \(\widehat\mu = \overline{Y}\).) Since \(\operatorname{E}(\overline{Y}) = \log\beta\),

\begin{equation*} \mathfrak{I} = \operatorname{E}(-U^{\prime}) = \frac{N}{\beta^{2}\sigma^{2}}, \qquad s.e.(\widehat\beta) = \sqrt{1/\mathfrak{I}} = \widehat\beta\sigma/\sqrt{N} \end{equation*}

by equation (4.12).

Equation (4.18). As a generalized linear model: \(p = 1\), \(\mathbf{X} = \mathbf{1}_N\), \(\eta_i = \beta\), \(\mu_i = \log\beta\), so the link required by (4.16) is \(g(\mu) = e^{\mu}\), the inverse of the mean function \(\mu = \log\eta\). With \(\operatorname{var}(Y_i) = \sigma^{2}\) and \(\partial\mu_i/\partial\eta_i = 1/\beta\),

\begin{align*} U_1 &= \sum_{i=1}^{N} \frac{(Y_i - \mu_i)}{\operatorname{var}(Y_i)} x_{i1} \frac{\partial\mu_i}{\partial\eta_i} = \frac{N}{\beta\sigma^{2}}\left(\overline{Y} - \log\beta\right), \\ \mathfrak{I}_{11} &= \sum_{i=1}^{N} \frac{x_{i1}^{2}}{\operatorname{var}(Y_i)} \left(\frac{\partial\mu_i}{\partial\eta_i}\right)^{2} = \frac{N}{\beta^{2}\sigma^{2}}, \end{align*}

matching the score and information above, which verifies (4.18) and (4.20).

Equation (4.25). The weights of (4.23) are \(w_{ii} = 1/(\sigma^{2}\beta^{2})\), so \(\mathbf{W} = (\sigma^{2}\beta^{2})^{-1}\mathbf{I}_N\) and \(\mathbf{X}^T\mathbf{W}\mathbf{X} = N/(\sigma^{2}\beta^{2}) = \mathfrak{I}\). Writing \(b = b^{(m-1)}\), the working response (4.24) is \(z_i = b + (y_i - \log b)b\), since \(\partial\eta_i/\partial\mu_i = \beta\) evaluated at \(b\). Hence

\begin{align*} \mathbf{X}^T\mathbf{W}\mathbf{z} &= \frac{1}{\sigma^{2}b^{2}} \sum_{i=1}^{N} \left[b + b(y_i - \log b)\right] = \frac{N}{\sigma^{2}b}\left[1 + \overline{y} - \log b\right], \\ b^{(m)} &= b^{(m-1)}\left[1 + \overline{y} - \log b^{(m-1)}\right], \end{align*}

the second line solving \(\mathbf{X}^T\mathbf{W}\mathbf{X}b^{(m)} = \mathbf{X}^T\mathbf{W}\mathbf{z}\). This is exactly the method of scoring (4.21), \(b + \mathfrak{I}^{-1}U = b[1 + \overline{y} - \log b]\), and its fixed point is \(\log b = \overline{y}\), i.e. \(b = \widehat\beta\). Note the step is genuinely iterative despite the Normal response, because \(g\) is not the identity, so \(\mathbf{W}\) and \(\mathbf{z}\) depend on \(b\).

Numerically, with \(N = 40\), \(\beta = 3\), \(\sigma = 0.5\):

set.seed(4)
N <- 40; beta.true <- 3; sigma <- 0.5
y <- rnorm(N, mean = log(beta.true), sd = sigma)

cat("closed form  exp(ybar) =", exp(mean(y)), "\n")
loglik <- function(b) sum(dnorm(y, mean = log(b), sd = sigma, log = TRUE))
cat("optimise     =", optimise(loglik, c(0.1, 20), maximum = TRUE)$maximum, "\n")

b <- 1                                     # IRLS recursion from (4.25)
for (m in 1:6) {
  b <- b * (1 + mean(y) - log(b))
  cat(sprintf("m = %d  b = %.8f\n", m, b))
}
cat("score (4.18) at b:", sum((y - log(b))/sigma^2) * (1/b), "\n")
cat("s.e. = sqrt(1/I):", sqrt(b^2 * sigma^2 / N), "\n")