Here are the books that I have taken the time to create metadata and/or notes for.
An Introduction to Generalized Linear Models
install.packages("glmnet")
install.packages("dobson")
library(glmnet)
library(dobson)
data(carbohydrate)
y = carbohydrate$carbohydrate
x = as.matrix(carbohydrate[,c('age', 'weight', 'protein')])
fit = glmnet(x, y)
plot(fit, xvar='lambda')
cvfit = cv.glmnet(x,y)
coef(cvfit, s = "lambda.1se")
Exercise 4.2
Each \(y_i\) is ’time to death’ for a particular leukemia patient (in weeks from diagnosis).
The \(x_i\) are \(\log_{10}\) initial white blood cell counts.
a.
x = c(3.36,2.88,3.63,3.41,3.78,4.02,4.00,4.23,3.73,3.85,3.97,4.51,4.54,5.00,5.00,4.72,5.00)
y = c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65)
plot(x,y)
The data shows a negative trend – as WBC increase, TTL decreases.
length(x) # 17 patients
[1] 17
b. \[ E(Y_i) = \exp(\beta_1 + \beta_2 x_i)\]
which ensures that \(E(Y)\) is non-negative for all values of the parameters, and all values of \(x\).
Which uses the log link. This is because survival time is strictly positive, and the exponential mean function guarantees \(E(Y) > 0 \quad\forall \beta, x\) (unlike OLS which can predict negative survival times).
Furthermore, the variance likely scales with the mean (since survival times are right-skewed) – a constant variance assumption (OLS) is inappropriate.
Finally, the log link linearises a multiplicative relationship – unit increase in \(x_i\) multiplies \(E(Y)\) by \(e^{\beta^2}\) which is more natural for this biological relationship.
Note, the negative trend implies \(\beta_2 <0\), along with the presence of non-constant variance.
c. \[f(y;\theta) =\theta e^{-y\theta}\]
This is the PDF of the Exponential Distribution, which is a special case of the Gamma distribution with shape parameter \(\phi = 1\)
Show that \(E(Y) = \cfrac1 \theta\) and \(\text{var}(Y)=\cfrac 1 {\theta^2}\)
Proof:
\(E(Y) = \int^\infty_0 y \cdot \theta e^{-y\theta} \mathrm{d}y\), then with IBP, we have \(u=y\) and \(dv=\theta e^{-y\theta}dy\), which yields: \[\left[-ye^{-y\theta} \right]^\infty_0 + \int^\infty_0 e^{-y\theta} dy = 0 + \left[-\cfrac1\theta e^{-y\theta}\right]^\infty_0=\cfrac1\theta\] Similarly, for the variance, we have:
\(E(Y^2) = \int_0^\infty y^2 \cdot \theta e^{-y\theta} \, dy\), but require 2 applications of IBP (or cheating with Gamma integral: \(\int_0^\infty y^n e^{-y\theta}dy = \frac{n!}{\theta^{n+1}}\)):
\[E(Y^2) = \theta \cdot \frac{2!}{\theta^3} = \frac{2}{\theta^2}\]
Whence, \[\text{var}(Y) = E(Y^2) - [E(Y)]^2 = \frac{2}{\theta^2} - \frac{1}{\theta^2} = \boxed{\frac{1}{\theta^2}}\] Realise that \(\text{var}(Y) = [E(Y)]^2\), confirming the Gamma variance structure \(\text{var}(Y) = \phi \cdot [E(Y)]^2\) with \(\phi=1\).
d. Fit a model with the equation for \(E(Y_i)\):
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
Here we use Gamma(link = "log") since the exponential distribution is a special case of the Gamma family with \(\phi = 1\), and the log link corresponds to the mean function \(E(Y_i) = \exp(\beta_1 + \beta_2 x_i)\).
We have \(\hat\beta_2 \approx -1.11 < 0\), confirming the negative trend observed in part (a) – increasing WBC reduces expected survival time multiplicatively.
e. For the exponential / Gamma model with \(\phi =1\): \(\text{var}(Y_i) = [E(Y_i)]^2\), so the standard deviation equals the mean: \(\text{sd}(Y_i) = E(Y_i)\).
Dividing the raw residual by \(\hat y_i\) therefore divides by the estimated standard deviation, which is what makes
\[r_i = \cfrac{y_i - \hat y_i} {\hat y_i}, \qquad \hat y_i = \exp(\hat\beta_1 + \hat\beta_2 x_i)\]
a standardised residual here (rather than an ad hoc scaling). Indeed, these are the Pearson residuals for the Gamma family.
b1 <- coef(fit)[1]
b2 <- coef(fit)[2]
yhat <- exp(b1 + b2*x) # we can also do this with fitted(fit)
r <- (y - yhat) / yhat # again; 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"); abline(h = 0, lty =2)
plot(x, r, xlab = "x (log10 WBC)", ylab=expression(r[i]),
main = "Residuals vs x"); 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
Adequacy of the model.
Centre and spread. The residuals are centred near zero (by construction) with \(\text{sd}(r_i) \approx 0.94\). Because the model implies \(\text{sd}(Y_i) = E(Y_i)\), a “typical” residual has magnitude \(\approx 1\), and indeed all but one observation fall in roughly \((-1, 1.13)\). This is consistent with the exponential assumption – observations routinely sit a full standard deviation from their fitted mean, which for this distribution is expected rather than alarming.
Link to the dispersion parameter. Since these are Pearson residuals, \(\hat\phi = \frac{1}{n-p}\sum r_i^2 = \frac{1}{15}\sum r_i^2 \approx 0.939\), matching the dispersion parameter reported in part (d). The value being close to \(1\) supports the choice \(\phi = 1\) (exponential), i.e. the Gamma variance structure with a fixed shape is reasonable.
Structure. Plotting \(r_i\) against \(\hat y_i\) and against \(x_i\) shows no obvious funnelling (non-constant spread) and no systematic curvature, so the log link and the mean function \(E(Y_i) = \exp(\beta_1 + \beta_2 x_i)\) appear appropriate. Modelling on the log scale has successfully removed the mean-variance dependence visible in the raw data of part (a).
One outlier. The last observation (\(x = 5.00\), \(y = 65\)) gives \(r_{17} = 2.47\), well outside the others. Three patients share the highest WBC \(x = 5.00\); two survived \(1\) week while the third survived \(65\), so the model (fitting \(\hat y = 18.75\) for all three) cannot accommodate that spread. This is a genuinely surprising survival, not a model failure, but it inflates the residual sum of squares and is worth flagging.
Overall the exponential model is adequate: the standardised residuals are unstructured, their spread is close to \(1\), and the estimated dispersion is close to the assumed \(\phi = 1\). The single high-leukocyte survivor is the only notable departure.