A smooth curve is not evidence of a faithful model. Once the optimizer has converged, we still have five decisions to make. Does the fitted shape follow the empirical survival? Is a continuous covariate sensible on the scale we gave it? Does predicted risk agree with observed experience across the cohort? Would variable selection or resampling tell a different story? And have we defined competing outcomes with the right risk set?
These checks form a sequence, not a collection of stamps to put on a model. Global misfit sends us back to the hazard shape. Curvature in a grouped functional-form check sends us back to the covariate transformation. Poor decile agreement sends us back to the linear predictor or baseline hazard. Unstable selection or bootstrap estimates tell us to simplify the model or collect more information. None of them turns an observational association into a causal effect.
Everything here uses public data bundled with TemporalHazard(Ehrlinger 2026). The chapter is independently runnable and requires no patient extract.
24.2 Refit the reviewed model
The AVC example follows 305 complete records after atrioventricular canal repair. Follow-up is in months. We refit the reviewed Weibull model so every diagnostic below refers to one known object and one analysis cohort.
data(avc, package ="TemporalHazard")avc <-na.omit(avc)fit <-hazard( survival::Surv(int_dead, dead) ~ age + status + mal + com_iv,data = avc,dist ="weibull",theta =c(mu =0.20, nu =1.0, rep(0, 4)),fit =TRUE,control =list(maxit =500))data.frame(records =nrow(avc),deaths =sum(avc$dead),maximum_follow_up_months =max(avc$int_dead),converged = fit$fit$converged,log_likelihood = fit$fit$objective)
Stop here if converged is not TRUE. A diagnostic cannot repair a fit that never reached its optimum.
24.3 Does the overall shape fit?
hzr_gof() places the fitted Weibull survival beside the Kaplan-Meier estimate and accumulates observed and expected events over follow-up. With covariates in the model, the parametric curve is evaluated at the mean design values. The Kaplan-Meier curve is marginal over the cohort. Thus this is a first-pass shape and conservation check, not patient-specific calibration.
records observed_events expected_events_at_mean_profile
1 305 68 68.00001
final_residual_expected_minus_observed
1 1.116152e-05
ggplot() +geom_step(data = gof,aes(time, km_surv *100, colour ="Kaplan-Meier, cohort",linetype ="Kaplan-Meier, cohort"),linewidth =0.75 ) +geom_line(data = gof,aes(time, par_surv *100, colour ="Weibull, mean profile",linetype ="Weibull, mean profile"),linewidth =0.75 ) +scale_colour_manual(values =c("Kaplan-Meier, cohort"="grey35","Weibull, mean profile"="steelblue"),name =NULL ) +scale_linetype_manual(values =c("Kaplan-Meier, cohort"="solid","Weibull, mean profile"="dashed"),name =NULL ) +scale_y_continuous(limits =c(0, 100), breaks =seq(0, 100, 20),labels =function(x) paste0(x, "%")) +labs(x ="Months after AVC repair", y ="Freedom from death") +theme_hv_manuscript(base_size =12) +theme(legend.position ="top")
Figure 24.1: Kaplan-Meier survival for the AVC cohort compared with Weibull survival evaluated at the mean covariate profile. The solid grey step line is Kaplan-Meier and the dashed steel-blue line the Weibull curve.
The curves differ, and the summary reports fewer expected events at the mean covariate profile than were observed in the cohort. Those quantities have different estimands. One is conditional on a single reference profile, while the other is marginal over the cohort. The discrepancy therefore reflects an estimand mismatch, not underprediction or a misspecified hazard shape. Before concluding that the model is miscalibrated, standardize predictions over the cohort or use a subject-specific check such as the decile observed-to-expected comparison below.
24.4 Is age sensible on the logit scale?
The name hzr_calibrate() comes from the historical grouped-screening workflow, but its current contract is narrower than a fitted survival-model calibration. It bins a continuous x, counts a binary event within each group, and transforms that observed proportion. It does not use event time or censoring. Here it asks whether eventual death is roughly linear in age on the logit scale, which is a functional-form screen for the next model fit.
cal <-hzr_calibrate(x = avc$age,event = avc$dead,groups =10,link ="logit")cal
ggplot(cal, aes(mean, link_value)) +geom_point(size =2.4, colour ="steelblue") +geom_line(colour ="steelblue", linewidth =0.6) +geom_smooth(method ="lm", formula = y ~ x, se =FALSE,colour ="grey35", linetype ="dashed", linewidth =0.65) +labs(x ="Mean age at repair within group (months)",y ="Logit of observed death proportion") +theme_hv_manuscript(base_size =12)
Figure 24.2: Observed binary death proportion by age group, shown on the logit scale as a functional-form check. Steel-blue filled circles mark the age groups, a steel-blue solid line joins them, and the dashed grey line is the linear fit.
Strong curvature would argue for a transformation, spline, or clinically defined age categories before refitting. A straight pattern would support a linear age term, but it would not prove that the full survival model is calibrated.
24.5 Does risk-group agreement hold?
hzr_deciles() asks a different question. It ranks all 305 patients by predicted survival at 120 months and divides them into ten nearly equal groups. The horizon defines the ranking only. Within each group, observed events use each patient’s full follow-up and expected events sum predicted cumulative hazard at each patient’s own follow-up time. No patient is excluded for having less than 120 months of observation.
dec <-hzr_deciles(fit, time =120, groups =10)dec_overall <-attr(dec, "overall")data.frame(risk_horizon_months = dec_overall$time,included = dec_overall$n_included,excluded = dec_overall$n_excluded,observed_events = dec_overall$total_events,expected_events = dec_overall$total_expected,chi_square = dec_overall$chi_sq,degrees_freedom = dec_overall$df,p_value = dec_overall$p_value)
Figure 24.3: Full-follow-up observed and expected event rates within groups ranked by predicted 120-month risk. Grey bars are the observed rate and steel-blue bars the expected rate, side by side in each group.
The omnibus result is compatible with agreement across these groups, but a non-significant test is not proof of calibration. Read the group sizes and the observed-versus-expected pattern with the statistic. Sparse events and the choice of groups both limit what this check can see.
24.6 Does automated selection tell a stable story?
For the selection demonstration, we start from the reviewed early-plus- constant model. The scope is phase-specific, so age can enter the early phase without being forced into the constant phase. We request two-way selection with the refit-based Wald criterion and keep the optimizer bounded.
base_mp <-hazard( survival::Surv(int_dead, dead) ~1,data = avc,dist ="multiphase",phases =list(early =hzr_phase("cdf", t_half =0.5, nu =1, m =1,fixed ="shapes" ),constant =hzr_phase("constant") ),fit =TRUE,control =list(n_starts =3, maxit =500))selection_warnings <-character()selected <-withCallingHandlers(hzr_stepwise( base_mp,scope =list(early =~ age + status + mal + com_iv,constant =~ age + status + mal + com_iv ),data = avc,direction ="both",criterion ="wald",trace =FALSE,control =list(n_starts =2, maxit =500) ),warning =function(w) { selection_warnings <<-c(selection_warnings, conditionMessage(w))invokeRestart("muffleWarning") })trace <-stepwise_trace(selected)cat(paste(trace, collapse ="\n"))
Stepwise selection (direction = both, criterion = wald, slentry = 0.30, slstay = 0.20)
Step 1: ENTER status into early (p = 0.000)
Step 2: ENTER com_iv into early (p = 0.000)
Step 3: ENTER status into constant (p = 0.065)
Step 4: ENTER mal into early (p = 0.109)
Step 5: ENTER age into early (p = 0.117)
(no further action after 5 steps)
Final model: 5 covariates, logLik = -192.10, AIC = 398.21
The trace names the criterion and records each accepted addition or removal. This run adds five phase-specific terms. The final model converges with a positive-definite covariance matrix, but intermediate candidate fits produce ill-conditioned Hessians and non-positive variances. That instability is part of the selection result. The entry and stay thresholds are a search rule, not a correction for repeated testing, and the selected path is conditional on this sample. Preserve clinically required covariates with force_in, examine the fit after every accepted step, and do not describe a selected coefficient as a causal effect.
24.7 How much do the estimates move under resampling?
hzr_bootstrap() resamples the analysis rows, refits the same Weibull formula, and summarizes the fitted parameters among successful replicates. The seed makes this teaching run reproducible. Thirty replicates are enough to exercise the real refit path in a book render; they are not enough for study inference.
Figure 24.4: Teaching-run percentile intervals for the four Weibull covariate coefficients across 30 bootstrap refits. Steel-blue filled circles mark the bootstrap mean, horizontal steel-blue bars the percentile interval, and the dashed grey vertical line zero.
The mal interval crosses zero in this small run, while the other displayed intervals do not. That is a stability signal, not a final confidence interval. For a study analysis, prespecify a substantially larger resample, review the failure count, and report the resampling method with the model specification.
24.8 Are competing events using the right denominator?
Competing-risk incidence is built from a shared risk set. At each event time, n_risk is the number still free of every modeled event just before the events at that time. A cause-specific cumulative-incidence curve adds the current event-free survival multiplied by that cause’s event proportion. The curves therefore share one probability budget: event-free survival plus every cause’s cumulative incidence equals 1.
The public omc extract contains death/last-contact follow-up and indicators that one or more thromboembolic events were recorded, but it does not contain the thromboembolism date. We use the required three-state construction to show the estimator mechanics. Do not interpret the resulting type 1 curve as a clinical time-to-thromboembolism estimate. For that analysis, supply the time to first thromboembolism or death and code the first observed event.
cif <-rbind(data.frame(time = cr$time, incidence = cr$incid_1 *100,cause ="Any thromboembolism recorded"),data.frame(time = cr$time, incidence = cr$incid_2 *100,cause ="Death without recorded thromboembolism"))ggplot(cif, aes(time, incidence, colour = cause, linetype = cause)) +geom_step(linewidth =0.75) +scale_colour_manual(values =c("Any thromboembolism recorded"="steelblue","Death without recorded thromboembolism"="firebrick"),name =NULL ) +scale_linetype_manual(values =c("Any thromboembolism recorded"="solid","Death without recorded thromboembolism"="dashed"),name =NULL ) +labs(x ="Months to death or last contact",y ="Cumulative incidence in constructed example (%)") +theme_hv_manuscript(base_size =12) +theme(legend.position ="top")
Figure 24.5: Estimator demonstration using OMC follow-up classification; the public extract does not contain the thromboembolism time required for clinical incidence. The solid steel-blue step line is any thromboembolism recorded and the dashed firebrick step line death without recorded thromboembolism.
If the probability sum is not 1 to numerical precision, the event coding or risk-set construction needs attention. In a real competing-risk analysis, also verify that every row carries the time to the first event of any modeled cause, not a later administrative time.
