22  Parametric Proportional Hazards Regression

22.1 The model

The parametric PH regression model specifies both the baseline hazard \(h_0(t)\) and the effect of explanatory variables: \[ h(t \mid \mathbf{x}_i) = h_0(t)\exp(\bbeta^T \mathbf{x}_i), \] where \(\mathbf{x}_i = (x_{i1}, \ldots, x_{ip})^T\) is the vector of explanatory variables for individual \(i\).

We use a Weibull baseline hazard \(h_0(t) = \lambda\gamma t^{\gamma-1}\), giving \[ h(t \mid \mathbf{x}_i) = \lambda\gamma t^{\gamma-1} \exp(\bbeta^T \mathbf{x}_i). \]

The survival function is \[ S(t \mid \mathbf{x}_i) = \exp\!\left(-\lambda t^\gamma \exp(\bbeta^T \mathbf{x}_i)\right). \]

22.2 Log-likelihood

Using the censored likelihood from Chapter 19, with \(h(t_i \mid \mathbf{x}_i) = \lambda\gamma t_i^{\gamma-1}\exp(\bbeta^T\mathbf{x}_i)\) and \(S(t_i \mid \mathbf{x}_i)\) as above:

\[ \ell(\lambda,\gamma,\bbeta) = \sum_{i=1}^n \Bigl[\delta_i\bigl\{\bbeta^T\mathbf{x}_i + \log(\lambda\gamma) + (\gamma-1)\log t_i\bigr\} - \lambda t_i^\gamma \exp(\bbeta^T\mathbf{x}_i)\Bigr]. \tag{22.1}\]

The MLEs \(\hat\lambda, \hat\gamma, \hat{\bbeta}\) are obtained numerically by maximising Equation 22.1.

22.3 Interpretation of coefficients

The hazard ratio for a one-unit increase in explanatory variable \(x_k\), with all other explanatory variables fixed, is \[ \psi_k = e^{\beta_k}. \]

  • \(\beta_k > 0\): higher \(x_k\) increases the hazard (worse prognosis).
  • \(\beta_k < 0\): higher \(x_k\) decreases the hazard (better prognosis).
  • \(\beta_k = 0\): \(x_k\) has no effect on survival.

For a binary explanatory variable (e.g., group indicator \(x_k \in \{0, 1\}\)): \(\psi_k = e^{\beta_k}\) is the hazard ratio between groups.

22.4 Fitting Weibull PH models in R

The survreg() function uses the accelerated failure time (AFT) parametrisation. For a Weibull model, we must convert from the AFT output to the PH parametrisation.

Code 
library(survival)

# Fit Weibull PH regression
fit <- survreg(Surv(t, delta) ~ x1 + x2, data = mydata,
               dist = "weibull")
summary(fit)

# Convert AFT to PH parametrisation:
gamma_hat  <- 1 / fit$scale            # shape parameter
lambda_hat <- exp(-gamma_hat * coef(fit)[1])  # intercept -> baseline
beta_hat   <- -gamma_hat * coef(fit)[-1]      # explanatory variables (note sign)
psi_hat    <- exp(beta_hat)            # hazard ratios

# 95% CI for hazard ratios (via delta method):
se_beta <- gamma_hat * sqrt(diag(vcov(fit)))[-1]
ci_lo   <- exp(beta_hat - 1.96 * se_beta)
ci_hi   <- exp(beta_hat + 1.96 * se_beta)
Note: Sign convention in survreg()

In survreg(), the explanatory-variable effects \(\tilde\beta_k\) are in the AFT parametrisation: \(\log T = \mu + \tilde\bbeta^T\mathbf{x} + \sigma\varepsilon\). The relationship to the PH coefficient is \(\beta_k^{PH} = -\hat\gamma \cdot \tilde\beta_k\), so the signs are reversed when converting.

Example: Breast cancer data

For the breast cancer data from Chapter 20, fitting a Weibull model with a group indicator (\(x = 0\): negative staining, \(x = 1\): positive): \[ \hat\gamma \approx 1.0 \quad (\text{essentially exponential}), \qquad \hat\psi = e^{\hat\beta} \approx 2.59. \]

The estimated hazard ratio matches the two-sample exponential result, confirming consistency. The 95% CI \([0.98, 6.87]\) is just consistent with \(\psi = 1\).

22.5 Multiple explanatory variables: Weibull regression

With multiple continuous or categorical explanatory variables:

Code 
library(survival)

# Weibull model with two explanatory variables
fit <- survreg(Surv(t, delta) ~ age + factor(treatment),
               data = mydata, dist = "weibull")
summary(fit)

# Extract hazard ratios for each explanatory variable
gamma_hat <- 1 / fit$scale
beta_ph   <- -gamma_hat * coef(fit)[-1]
exp(beta_ph)  # hazard ratios

22.6 Model checking

After fitting, check:

  1. Log-cumulative hazard plot: should be linear in \(\log t\) for each group or combination of explanatory variables (Weibull check).
  2. Parallel log-cumulative hazard curves: should be parallel for different groups (PH check).
  3. Residual plots: deviance or Cox–Snell residuals should be approximately exponentially distributed if the model fits.
Code 
# Log-cumulative hazard by group
km_g <- survfit(Surv(t, delta) ~ group, data = mydata)
plot(km_g, fun = "cloglog", col = c("blue", "red"),
     xlab = "log(t)", ylab = "log H(t)")

Exercises

22.1. Consider the survival model with hazard \(h(t) = \lambda e^t\) (from Exercise 19.1) extended to a proportional hazards regression model with a single binary explanatory variable \(x_i \in \{0,1\}\): \[h_i(t) = e^{\beta x_i} \cdot \lambda e^t, \qquad i = 1,\ldots,n.\] There are \(r\) uncensored event times and \(n-r\) right-censored observations.

Write the full log-likelihood \(\ell(\lambda,\beta) = \sum_i \delta_i\log h_i(t_i) + \sum_i \log S_i(t_i)\). Note that \(S_i(t) = \exp(-e^{\beta x_i}\lambda(e^t-1))\). Differentiate with respect to \(\lambda\) and \(\beta\) separately and set both partial derivatives to zero.

  1. Write down the log-likelihood \(\ell(\lambda, \beta)\).
  2. Derive the two score equations \(\partial\ell/\partial\lambda = 0\) and \(\partial\ell/\partial\beta = 0\).
  3. Show that the two equations can be reduced to a single equation in \(\beta\) alone, plus a formula expressing \(\hat\lambda\) in terms of \(\hat\beta\).