Model checking · added in 2026

How far can the 2023 models be trusted?

The assignments fitted and interpreted models; they did little to check them. This page wraps the original fits in the checks I would run today: residuals against a simulated envelope, influence, over-dispersion, profile-likelihood intervals, a test of the proportional-odds assumption and the sensitivity of the GEE analysis. The 2023 estimates are not changed; every simulation shows its seed.

Assignment 2 · Q1 · binomial GLM

Beetle dose–response: fit, influence and intervals

The straight-line logistic model is the one the assignment asked for. Its residual deviance (13.63 on 6 df, p = 0.0340) already hinted at a problem; these checks say what kind.

Half-normal plot with a simulated envelope

Observed (sorted)
Outside the envelope
Simulated median
95% pointwise envelope

3 of 8 points outside the envelope · B = 999 simulations from the fitted straight-line model, seed 20230501, no failed refits.

Points outside the envelope across five fixed seeds (B = 999)
Seed	Straight line	Quadratic
20230501	3 of 8	0 of 8
1	3 of 8	0 of 8
42	3 of 8	0 of 8
90139	3 of 8	0 of 8
271828	3 of 8	0 of 8

Pearson dispersion, straight line

2.02

X² / df on 6 df

Pearson dispersion, quadratic

1.00

on 5 df

Quasi-binomial comparison (straight line)

(Intercept)

Estimate: −60.10
SE, binomial · quasi: 5.16 · 7.34
95% CI (quasi, t₆): (−78.1, −42.1)

Dosage

Estimate: 33.93
SE, binomial · quasi: 2.90 · 4.12
95% CI (quasi, t₆): (23.8, 44.0)

Term	Estimate	SE (binomial)	SE (quasi)	95% CI (quasi, t₆)
(Intercept)	−60.10	5.16	7.34	(−78.1, −42.1)
Dosage	33.93	2.90	4.12	(23.8, 44.0)

Under the quasi-likelihood the quadratic term is still needed, though less emphatically: F = 8.55 on (1, 5) df, p = 0.0329, against p = 0.0035 for the binomial likelihood-ratio test.

Influence: Cook's distance and leverage

With only 8 groups every point matters. Dose 1.76 has Cook's D = 1.03 and dose 1.72 has 0.87; both sit where the straight line misses the S-shape. No group has leverage above 2p/n = 0.50 (largest 0.35). Influence here is a symptom of the curvature rather than of a single bad batch.

LD50: four intervals, then two other mean models

The first four rows are intervals for the 2023 straight-line estimate. The last two refit the LD50 under the mean models the checks above point to.

Wald (delta method)1.7712 (1.7636 to 1.7788)
Profile likelihood1.7712 (1.7634 to 1.7787)
Parametric bootstrapB = 2,000 · seed 202305011.7712 (1.7638 to 1.7790)
Quasi-binomial (delta, t₆)1.7712 (1.7577 to 1.7847)
Quadratic logitother mean model · profile1.7798 (1.7703 to 1.7884)
Complementary log-logother mean model · profile1.7782 (1.7699 to 1.7859)

LD50 (log₁₀ dosage)LD50 (log₁₀ dosage); dashed line: 2023 estimate

2023 straight-line logit model, estimate 1.7712
Method	95% interval	Width
Wald (delta method)	(1.7636, 1.7788)	0.0152
Profile likelihood	(1.7634, 1.7787)	0.0153
Parametric bootstrap (seed 20230501)	(1.7638, 1.7790)	0.0151
Quasi-binomial delta	(1.7577, 1.7847)	0.0270

Other mean models

Straight line, logit (2023)

Deviance (df) · AIC: 13.63 (6) · 43.8
LD50, 95% profile CI: 1.7712 (1.7634, 1.7787)

Quadratic, logit

Deviance (df) · AIC: 5.11 (5) · 37.3
LD50, 95% profile CI: 1.7798 (1.7703, 1.7884)

Straight line, complementary log-log

Deviance (df) · AIC: 5.63 (6) · 35.8
LD50, 95% profile CI: 1.7782 (1.7699, 1.7859)

Mean model	Deviance (df)	AIC	LD50	95% profile CI
Straight line, logit (2023)	13.63 (6)	43.8	1.7712	(1.7634, 1.7787)
Quadratic, logit	5.11 (5)	37.3	1.7798	(1.7703, 1.7884)
Straight line, complementary log-log	5.63 (6)	35.8	1.7782	(1.7699, 1.7859)

Under the 2023 model the three binomial intervals agree to within 0.0005: with a slope this well determined the log-likelihood is close to quadratic. Across five seeds the bootstrap endpoints range over (1.7633, 1.7638) (lower) and (1.7781, 1.7790) (upper), with 0 failed refits.

All four are conditional on a mean model the diagnostics reject. Adding the quadratic term moves the LD50 to 1.7798, a shift of 0.0086, more than the Wald half-width of 0.0076; it lies outside all three binomial intervals. A complementary log-log line, which fits about as well as the quadratic with one parameter fewer (AIC 35.8 against 37.3), puts it at 1.7782. The quasi-binomial interval (1.8 times as wide as Wald) covers both, but by inflating the variance rather than fixing the curve. Model choice moves the estimate by up to 0.0086, which is the uncertainty the straight-line intervals leave out.

Profile likelihood: interval method against model

Profile, 2023 model
Wald, 2023 model
Profile, quadratic logit
Profile, complementary log-log
±1.96: 95% interval endpoints

Each curve crosses ±1.96 at its 95% profile interval. For the 2023 model the profile curve and the Wald line almost coincide, so the choice of interval method hardly matters. The curves for the other two mean models sit to the right: changing the model moves the whole interval, by more than any change of method.

Parameter	Wald 95%	Profile 95%
$\beta_0$	(−70.22, −49.98)	(−70.82, −50.51)
$\beta_1$	(28.24, 39.62)	(28.54, 39.96)
LD50	(1.7636, 1.7788)	(1.7634, 1.7787)

The profile interval for the slope is shifted up by about 0.3 relative to Wald, the mild skew a binomial likelihood usually has. Profile endpoints are solved exactly by root-finding on constrained refits (MASS's confint() interpolates a grid and agrees to about 10⁻³).

Assignment 3 · Q1–Q2 · ordinal response

Pneumoconiosis: is proportional odds justified?

In 2023 the proportional-odds model was preferred on AIC without testing its key assumption: that years of exposure shift every cumulative split by the same amount. Two tests and a picture of the splits (see DR-002).

Cumulative logits by exposure group

Mild or severe (vs normal)
Severe (vs normal or mild)
Proportional odds (common slope)

Solid lines: proportional-odds fit. Dashed: separate logistic fit for each split. Whiskers: Wilson 95% intervals mapped to the logit scale (cut at the plot edge where a group has no cases above the split).

Brant test of parallel slopes

p = 0.82

X² = 0.051 on 1 df

LR test vs separate slopes

p = 0.99

G² = 0.0002 on 1 df

Slope per year, fitted separately for each split
Split	Slope	95% CI
Mild or severe (vs normal)	0.0963	(0.0720, 0.1205)
Severe (vs normal or mild)	0.0935	(0.0632, 0.1237)
Proportional odds (common)	0.0959	(0.0725, 0.1193)

Assignment 3 · Q3 · correlated binary data

Wheeze: does the working correlation matter?

GEE estimates stay consistent whatever working correlation is assumed, and the sandwich standard errors stay valid; the choice affects efficiency. Here are the three usual structures, a bootstrap over children, and the random-intercept GLMM discussed in DR-003.

Working-correlation sensitivity

Independence

α: —
Smoking β (robust SE): 0.272 (0.178)
Robust / model SE: 1.44
Smoking OR, 95% CI: 1.31 (0.93, 1.86)
Age β (robust SE): −0.113 (0.044)
QIC · CIC: 1829.49 · 4.80

Exchangeable(2023 choice)

α: 0.354
Smoking β (robust SE): 0.265 (0.178)
Robust / model SE: 1.00
Smoking OR, 95% CI: 1.30 (0.92, 1.85)
Age β (robust SE): −0.113 (0.044)
QIC · CIC: 1829.48 · 4.80

AR(1)

α: 0.491
Smoking β (robust SE): 0.234 (0.181)
Robust / model SE: 1.02
Smoking OR, 95% CI: 1.26 (0.89, 1.80)
Age β (robust SE): −0.115 (0.045)
QIC · CIC: 1830.26 · 4.99

Working correlation	α	Smoking β (robust SE)	Robust / model SE	Smoking OR, 95% CI	Age β (robust SE)	QIC	CIC
Independence	—	0.272 (0.178)	1.44	1.31 (0.93, 1.86)	−0.113 (0.044)	1829.49	4.80
Exchangeable(2023 choice)	0.354	0.265 (0.178)	1.00	1.30 (0.92, 1.85)	−0.113 (0.044)	1829.48	4.80
AR(1)	0.491	0.234 (0.181)	1.02	1.26 (0.89, 1.80)	−0.115 (0.045)	1830.26	4.99

QIC differs by less than one unit across structures, so the data do not prefer one. The smoking estimate moves from 0.272 to 0.234, well inside its standard error, and every interval for the odds ratio includes 1. QIC and CIC follow geepack's QIC(); smaller is better.

Maternal smoking: one effect, several analyses

Naive GLMignores clustering1.31 (1.03 to 1.67)
GEE, independencerobust SE1.31 (0.93 to 1.86)
GEE, exchangeablerobust SE1.30 (0.92 to 1.85)
GEE, AR(1)robust SE1.26 (0.89 to 1.80)
Cluster bootstrapB = 1,000 · seed 202306011.30 (0.91 to 1.83)
GLMM (child-specific)different scale1.49 (0.87 to 2.54)

Odds ratio (log scale)Odds ratio, smoking vs non-smoking mother (log scale)

Ignoring clustering gives a misleadingly narrow interval (p = 0.0275). Every analysis that respects it includes an odds ratio of 1.

Cluster bootstrap across seeds

Whole children are resampled with replacement (B = 1,000 per seed) and the exchangeable GEE is refitted. Sandwich interval for comparison: (0.92, 1.85).

Seed	Smoking OR, 95% percentile CI	Bootstrap SE (log OR)
20230601	(0.91, 1.83)	0.180
7	(0.92, 1.87)	0.183
2023	(0.90, 1.88)	0.186
90139	(0.90, 1.86)	0.180
314159	(0.91, 1.82)	0.181

Robust (sandwich) SE 0.178; 0 failed refits.

GEE or GLMM? Two scales for the same data

(Intercept)

GEE, population-averaged (robust SE): −1.880 (0.114)
GLMM, child-specific (SE): −3.101 (0.219)
GLMM attenuated: −1.916

age

GEE, population-averaged (robust SE): −0.113 (0.044)
GLMM, child-specific (SE): −0.176 (0.068)
GLMM attenuated: −0.108

smoke

GEE, population-averaged (robust SE): 0.265 (0.178)
GLMM, child-specific (SE): 0.399 (0.273)
GLMM attenuated: 0.246

Between-child spread

GEE working α: 0.354
GLMM σ (SE): 2.16 (0.18)
Latent ICC: 0.59

Coefficient	GEE, population-averaged (robust SE)	GLMM, child-specific (SE)	GLMM attenuated
(Intercept)	−1.880 (0.114)	−3.101 (0.219)	−1.916
age	−0.113 (0.044)	−0.176 (0.068)	−0.108
smoke	0.265 (0.178)	0.399 (0.273)	0.246
Between-child spread	working α = 0.354	σ = 2.16 (SE 0.18)	ICC 0.59

The GLMM gives each child a normal random intercept. Its coefficients compare two children with the same underlying propensity to wheeze, so they are larger than the population-averaged GEE coefficients whenever children differ a lot (here σ = 2.16, a latent intraclass correlation of 0.59).

Attenuating by $1/\sqrt{1 + 0.346\,\sigma^2}$ (Zeger, Liang and Albert 1988) brings the smoking effect from 0.399 to 0.246, close to the GEE value of 0.265: the two models tell the same story on different scales. The 2023 analysis wanted the population-level contrast, so GEE stays the main analysis (DR-003).

Maximum likelihood with 100-point Gauss–Hermite quadrature; agrees with lme4::glmer(nAGQ = 25) to about 10⁻³ (glmer's optimiser stops earlier).

Where these checks come from

Every check on this page is computed by TypeScript code tested against R (base stats, MASS, VGAM, brant, geepack and lme4) by scripts/diagnostics-reference.R. See Methods for the evaluation design and Verification for the parity table.