25  Model Selection and Comparison

25.1 Nested models and the change in deviance

Recall from Part I (Section 4.3) that nested models can be compared via the change in deviance. The same idea applies directly to survival models.

Two models \(M_0 \subset M_1\) are nested if \(M_0\) is obtained from \(M_1\) by fixing one or more parameters at specified values, not necessarily zero. (For example, the exponential model is nested within the Weibull model by fixing the shape parameter \(\gamma = 1\).)

For a fitted model \(M\) with maximised log-likelihood \(\hat\ell\), its deviance is \(D = -2\hat\ell\) (up to an additive constant that cancels whenever we compare two nested models). Then, for nested models \(M_0\) (with \(p\) parameters) and \(M_1\) (with \(q > p\) parameters),

\[ \Delta D = D_0 - D_1 \;\dot\sim\; \chi^2_{q - p} \quad \text{under } H_0, \]

exactly as in Section 4.3. We reject \(H_0\) (prefer \(M_1\)) if \(\Delta D > \chi^2_{q-p,\,1-\alpha}\).

For Cox models, the same test applies using the deviance based on the maximised partial log-likelihood. The asymptotic \(\chi^2\) distribution holds under the same conditions.

No separate dispersion step needed. Unlike the Normal/Gamma case in Section 4.3, where an unknown dispersion \(\phi\) required an \(F\)-test with \(\hat\phi\) estimated from a larger model, no such step is needed here. A parametric survival model’s shape parameter (e.g. \(\gamma\) for the Weibull) is a genuine parameter of the likelihood, not a separate dispersion term scaling an otherwise \(\phi\)-free deviance; it is already fully accounted for in \(\hat\ell_0\) and \(\hat\ell_1\) once each model is fitted. The Cox partial likelihood has no dispersion parameter at all. In both cases the plain asymptotic \(\chi^2\) result applies directly, exactly as it does for the known-\(\phi\) GLM families (Poisson, Binomial).

25.2 Step-by-step model selection strategy

A practical approach for building a survival regression model:

1. Univariate screening: fit each explanatory variable separately; identify those that are significant at a liberal threshold (e.g., \(p < 0.10\)).

2. Multivariable model: include all significant explanatory variables from step 1. Drop non-significant ones using backward elimination.

3. Add borderline variables: consider adding explanatory variables not significant alone but that may be confounders; use LRT to decide.

4. Final check: ensure no term can be significantly added or removed.

5. Check interactions: for key explanatory variables, test interaction terms.

25.3 Case study: Nephrectomy data

The hypernephroma.dat dataset contains survival data for 36 patients with kidney cancer. Explanatory variables are age group (three levels: \(<60\), 60–70, \(>70\)) and whether nephrectomy (kidney removal) was performed.

library(survival)

h <- read.table("https://www.richardpmann.com/MATH3701/Datasets/hypernephroma.dat", header = TRUE)
h$age         <- as.factor(h$age)
levels(h$age)  <- c("<60", "60-70", ">70")
h$nephrectomy <- as.factor(h$nephrectomy)
levels(h$nephrectomy) <- c("no", "yes")

# Fit sequence of nested Cox models
M0 <- coxph(Surv(time, status) ~ 1,                    data = h, ties = "breslow")
M1 <- coxph(Surv(time, status) ~ age,                  data = h, ties = "breslow")
M2 <- coxph(Surv(time, status) ~ nephrectomy,          data = h, ties = "breslow")
M3 <- coxph(Surv(time, status) ~ age + nephrectomy,    data = h, ties = "breslow")
M4 <- coxph(Surv(time, status) ~ age * nephrectomy,    data = h, ties = "breslow")

# Extract -2 log partial likelihoods
minus2logL <- sapply(list(M0, M1, M2, M3, M4),
                     function(m) -2 * m$loglik[2])
round(minus2logL, 3)

The key comparisons:

Comparison	\(\Delta\)	df	\(p\)-value
\(M_1\) vs \(M_0\) (age)	\(177.7 - 172.2 = 5.5\)	2	0.064
\(M_2\) vs \(M_0\) (nephrectomy)	\(177.7 - 170.2 = 7.4\)	1	0.006
\(M_3\) vs \(M_1\) (nephrectomy | age)	\(172.2 - 165.5 = 6.7\)	1	0.010
\(M_4\) vs \(M_3\) (interaction)	\(165.5 - 162.5 = 3.0\)	2	0.220

Nephrectomy is significant, both alone and after adjusting for age (nephrectomy \(\mid\) age, \(p=0.010\)). Age on its own is only marginal at the 5% level (\(p=0.064\) for the joint, 2 df test across all three groups) – this joint test is diluted by the middle (60–70) group, which is not clearly different from the youngest group, even though the oldest group (\(>70\)) does show a clearly elevated hazard (see the confidence interval in the next section). The interaction is clearly non-significant (\(p=0.220\)); we retain the main-effects model \(M_3\) as the final model. This illustrates that even when an explanatory variable’s overall significance is marginal, individual category effects can still be clinically important, and that we should always verify the added value of an interaction term.

# Final model
summary(M3)
# Hazard ratios and 95% CIs
exp(cbind(coef(M3), confint(M3)))

25.4 Parametric vs Cox: which to use?

Criterion	Parametric (Weibull)	Cox
Baseline hazard	Specified	Unspecified
Efficiency	Higher if model correct	Lower
Robustness	Biased if \(h_0\) wrong	Robust
Prediction	Full \(\hat{S}(t)\) easy	\(\hat{S}(t)\) requires estimate of \(H_0(t)\)
Interpretability	Hazard ratios + shape	Hazard ratios only

Recommendation: start with a graphical check (log-cumulative hazard plot). If the parametric form looks appropriate, use the parametric model for efficiency. When in doubt, use the Cox model.

Exercises

25.1. Using the Drug A/B Cox partial likelihood from Exercise 23.1, perform a likelihood ratio test of \(H_0: \beta = 0\) against \(H_1: \beta \neq 0\).

Simplify \(\ell(\hat\beta) - \ell(0)\) by cancelling the \(\log(1/24)\) terms. Use \(\log(4/3) = \log 4 - \log 3\) and collect terms.

The partial log-likelihoods are: \[\ell(0) = \log\!\left(\tfrac{1}{24}\right) - 4\log 2, \qquad \ell(\hat\beta) = \log\!\left(\tfrac{1}{24}\right) - \log 3 - 4\log\!\left(\tfrac{4}{3}\right).\]

1. Compute the likelihood ratio test statistic \(T = -2(\ell(0) - \ell(\hat\beta))\).
2. State the null distribution of \(T\) and the critical value at the 5% level.
3. State your conclusion. Is it consistent with the Wald test in Exercise 23.2?

25.2. When analysing survival data, one must often choose between a fully parametric model (e.g., Weibull PH) and the semi-parametric Cox model.

1. Which of the following is an advantage of the Cox model over a Weibull PH model?

2. Advantage of the Cox model: No assumption about the shape of the baseline hazard is neededThe Cox model always gives more precise parameter estimatesThe Cox model can be fitted without any censored observationsThe Cox model gives explicit estimates of the survival function at arbitrary future times without additional estimation

3. Which of the following is an advantage of the Weibull PH model over the Cox model?

4. Advantage of the Weibull model: The Weibull model makes no distributional assumptionsWhen the model is correctly specified, parametric estimates are more efficient and allow direct survival prediction at any future timeThe Weibull model is always preferred when the sample size is largeThe Weibull model cannot be used when there are tied event times

5. Suppose you fit both a Weibull PH model and an exponential PH model to the same data. Since the exponential model is the Weibull with shape \(\gamma = 1\), the two models are nested. The deviances (\(-2 \times\) log-likelihood) are: Weibull \(= 83.4\), Exponential \(= 89.2\). Use a likelihood ratio test (i.e. the change in deviance) to decide whether the extra shape parameter is justified.