Example: hospital readmissions after colorectal cancer surgery
The readmission data (Gonzalez et al., 2005) contain repeated
hospital readmissions of 403 patients after surgery for colorectal
cancer, with death as competing terminal event.
readmission_prep() returns one row per readmission (status
1) and one final row per patient, either death (status 2) or censoring
(status 0). Time is measured in days since surgery.
rdata <- readmission_prep(tau = 1460)
head(rdata)
#> id time status chemo sex dukes charlson
#> 1 1 24 1 Treated Female D 3
#> 2 1 457 1 Treated Female D 0
#> 3 1 1037 0 Treated Female D 0
#> 4 2 489 1 NonTreated Male C 0
#> 5 2 1182 0 NonTreated Male C 0
#> 6 3 15 1 NonTreated Male C 3
table(rdata$status)
#>
#> 0 1 2
#> 296 447 107The covariates are chemotherapy, sex, Dukes’ tumor stage and the Charlson comorbidity index. For our analysis we set \(\tau = 4\) years (1460 days): readmissions after \(\tau\) are removed and patients still under observation are censored at \(\tau\). At \(\tau\), 25% of the patients are still followed and 98% of the readmissions have occurred, while beyond \(\tau\) the number of patients at risk drops rapidly (21 patients at 5 years).
Choosing the transformation
plot_loglik() plots the profile log-likelihood over a
grid of transformation parameters for both model classes, with \(r\) (logarithmic) on the left and \(\rho\) (Box-Cox) on the right. The open
circle marks the Ghosh–Lin model (\(\rho =
1\)), the filled circle the proportional odds model (\(r = 1\)).
plot_loglik(
Surv(time, status) ~ chemo + sex + dukes + charlson,
data = rdata,
id = "id",
tau = 1460,
rho_grid = seq(0.05, 1.2, by = 0.05),
r_grid = seq(0.05, 1.2, by = 0.05)
)
The maxima are at \(\hat\rho \approx
0.85\) and \(\hat r \approx
0.15\), both close to the Ghosh–Lin model. The exact optima can
be found with optimize():
Fitting the selected model
fit_rd <- wnpmle_fit(
Surv(time, status) ~ chemo + sex + dukes + charlson,
data = rdata,
id = "id",
model = "boxcox",
rho = 0.85,
tau = 1460,
se = "sandwich_adj"
)
summary(fit_rd)
#>
#> Weighted NPMLE - Recurrent Events with Competing Terminal Event
#> Type : recurrent
#> Model : BOXCOX transformation ( rho = 0.85 )
#> Subjects : 403
#> Events : recurrent = 447 terminal = 107 censored = 296
#> Log-lik : -3016.694
#> Convergence: relative convergence (4)
#>
#> Coefficients:
#> Estimate SE z value Pr(>|z|)
#> chemoTreated -0.6079 0.1861 -3.267 0.0010880
#> sexFemale -0.5159 0.1792 -2.879 0.0039920
#> dukesC 0.2823 0.1991 1.418 0.1562000
#> dukesD 0.9322 0.2765 3.371 0.0007483
#> charlson1-2 0.7336 0.4075 1.800 0.0718000
#> charlson3 -0.0805 0.1957 -0.411 0.6808000
#>
#> Cumulative baseline mean at time grid:
#> Lambda SE
#> A(tau/4) = 365 0.7016 0.1371
#> A(tau/2) = 730 1.0874 0.1994
#> A(tau) = 1460 1.5716 0.2711Comparison with the Ghosh–Lin and proportional odds models:
fit_gl <- wnpmle_fit(Surv(time, status) ~ chemo + sex + dukes + charlson,
data = rdata, id = "id", model = "boxcox", rho = 1,
tau = 1460, se = "none")
fit_po <- wnpmle_fit(Surv(time, status) ~ chemo + sex + dukes + charlson,
data = rdata, id = "id", model = "log", rho = 1,
tau = 1460, se = "none")
#> Note: Using Makevars in C:\Users\abell\AppData\Local\Temp\RtmpMZz5ds\file93cc2ded7de9
#> using C++ compiler: 'G__~1.EXE (GCC) 13.3.0'
sapply(list(boxcox_0.85 = fit_rd, ghosh_lin = fit_gl, prop_odds = fit_po), AIC)
#> boxcox_0.85 ghosh_lin prop_odds
#> 6045.388 6045.461 6054.987In the selected model, chemotherapy, sex and Dukes’ stage D have a significant effect on the expected number of readmissions.
A positive coefficient increases and a negative coefficient decreases
the expected number of recurrences at all times. In the Ghosh–Lin model
(\(\rho = 1\)), \(\exp(\beta)\) is the ratio of the marginal
means; in the proportional odds model (\(r =
1\)), \(\exp(\beta)\) is the
ratio of \(\exp(\mu(t)) - 1\), where
\(\mu(t)\) is the marginal mean. For
other transformation models, the size of the effect depends on time and
is best illustrated with predict(), as shown below.
Cumulative baseline mean
baseline() returns the cumulative baseline mean with
pointwise 95% confidence limits.
bl <- baseline(fit_rd)
plot(bl$time, bl$Lambda, type = "s", lwd = 2,
xlab = "Days since surgery", ylab = expression(hat(Lambda)(t)),
ylim = range(c(bl$lower, bl$upper), na.rm = TRUE))
lines(bl$time, bl$lower, type = "s", lty = 2, col = "grey50")
lines(bl$time, bl$upper, type = "s", lty = 2, col = "grey50")
Prediction
predict() gives the marginal mean with pointwise 95%
confidence limits at new covariate values, here for a man with Dukes’
stage C and Charlson index 0, with and without chemotherapy.
newdat <- data.frame(chemo = c("NonTreated", "Treated"),
sex = "Male",
dukes = "C",
charlson = "0")
pred <- predict(fit_rd, newdata = newdat, times = seq(0, 1460, by = 7))
head(pred)
#> time mu_1 lower_1 upper_1 mu_2 lower_2 upper_2
#> 1 0 0.00000000 0.000000000 0.00000000 0.00000000 0.000000000 0.00000000
#> 2 7 0.02214344 0.009756028 0.05025939 0.01206602 0.004954678 0.02938414
#> 3 14 0.05192327 0.028269828 0.09536760 0.02832098 0.014213503 0.05643071
#> 4 21 0.08541735 0.051230223 0.14241836 0.04664011 0.025176942 0.08640049
#> 5 28 0.10400615 0.064370436 0.16804732 0.05682323 0.031239694 0.10335824
#> 6 35 0.11891761 0.075519166 0.18725575 0.06500002 0.036408368 0.11604481
plot(pred$time, pred$mu_1, type = "n",
xlab = "Days since surgery",
ylab = "Marginal mean number of readmissions",
ylim = range(pred[, -1]))
polygon(c(pred$time, rev(pred$time)), c(pred$lower_1, rev(pred$upper_1)),
col = adjustcolor("black", 0.12), border = NA)
polygon(c(pred$time, rev(pred$time)), c(pred$lower_2, rev(pred$upper_2)),
col = adjustcolor("firebrick", 0.15), border = NA)
lines(pred$time, pred$mu_1, type = "s", lwd = 2)
lines(pred$time, pred$mu_2, type = "s", lwd = 2, lty = 2, col = "firebrick")
legend("topleft", legend = c("No chemotherapy", "Chemotherapy"),
lty = c(1, 2), col = c("black", "firebrick"), lwd = 2, bty = "n")