Fitting, prediction & simulation

Fitting and inference

A model is fitted from data with fit(Model, @formula(Surv(time, status) ~ ...), df) (or the two-formula variant for GeneralHazard); see Illustrative example for complete calls. The optimizer is seeded automatically from the data.

Fit statistics

GeneralHazardModel <: StatsAPI.StatisticalModel, so a fitted model supports the standard statistical-model accessors. aic, aicc, and bic follow from loglikelihood, dof, and nobs; stderror follows from vcov. coef and vcov are reported on the inference scale [log.(baseline parameters); active regression coefficients], so MvNormal(coef(m), vcov(m)) is a coherent parameter-uncertainty distribution.

For a fitted-model report, coeftable(m) gives the Wald table of covariate effects (with exp(coef) and its confidence interval) and confint(m) the coefficient intervals as a DataFrame; both exclude the baseline distribution parameters, which show reports as the fitted baseline instead.

StatsAPI.loglikelihoodMethod
loglikelihood(m::GeneralHazardModel)

Maximized log-likelihood of the fitted model, captured from the optimizer at fit time (NaN for models built by the direct, non-optimizing constructor).

source
StatsAPI.nobsMethod
nobs(m::GeneralHazardModel)

Number of observations (event/censoring times) the model was fit to.

source
StatsAPI.dofMethod
dof(m::GeneralHazardModel)

Number of free parameters consumed by the fit: the baseline distribution's parameters plus the regression coefficients that actually enter the hazard for this model's structure (β for PH/AFT, α for AH, both for GH).

source
StatsAPI.coefMethod
coef(m::GeneralHazardModel)

Identified parameters on the scale used for inference: [log.(baseline parameters); active regression coefficients], where the active coefficients are β (PH/AFT), α (AH), or [α; β] (GH). This ordering matches vcov, so MvNormal(coef(m), vcov(m)) is a coherent sampling distribution for the parameter uncertainty (exponentiate the baseline entries to recover the natural-scale baseline parameters).

source
StatsAPI.vcovMethod
vcov(m::GeneralHazardModel)

Covariance of coef: the inverse observed information, i.e. the inverse of the Hessian of the fitted negative log-likelihood at the optimum, restricted to the identified parameters (the inactive coefficient block makes the full Hessian singular). Computed on demand — LBFGS never forms the Hessian during the fit, so there is nothing cached to reuse. stderror follows from this via the generic StatisticalModel method.

source
StatsAPI.coeftableMethod
coeftable(m::GeneralHazardModel; level::Real=0.95)

Wald coefficient table for the covariate effects: estimate, standard error, z, p-value, exp(coef), and the confidence interval for exp(coef) at level level. exp(coef) is a hazard ratio for hazard-level effects and a time/acceleration ratio for time-scale effects. For a General Hazard model the rows are tagged with their role (hazard-level vs time-scale), since a covariate may enter both. Baseline distribution parameters are reported by show, not here: a null of zero is not the relevant hypothesis for a baseline shape/scale.

source
StatsAPI.confintMethod
confint(m::GeneralHazardModel; level::Real=0.95)

Wald confidence intervals for the covariate coefficients at confidence level level. Returns a DataFrame with columns component (the coefficient's role, "hazard-level" or "time-scale"), term, lower, and upper. Baseline distribution parameters are not included (see coeftable).

source

Brier score

Inverse-probability-of-censoring-weighted Brier score (Graf et al. 1999) and its integrated form work for GeneralHazardModel through the same brier_score(model, ...) / integrated_brier_score(model, ...) API used for Cox. See the Model Evaluation: Brier Score section of the Cox documentation for the mathematical definition and signature list.

Prediction

Once a GeneralHazardModel is fitted (or directly constructed), you can evaluate per-subject cumulative hazards and survival probabilities at user-supplied times. The four hazard structures share the same closed-form expression via the unified representation

\[H(t \,|\, x) = H_0\!\left(t\, c_1(x)\right) c_2(x)\]

where $H_0$ is the cumulative hazard of the baseline distribution and $(c_1, c_2)$ are the method-specific time- and hazard-scale multipliers (c1/c2 in the code). The survival is $S(t \,|\, x) = \exp(-H(t \,|\, x))$.

predict(model, :survival)              # length-n vector, each subject at own Tᵢ
predict(model, :expected)              # length-n vector of Λᵢ(Tᵢ)
predict(model, :survival, t)           # length-n vector at scalar t
predict(model, :expected, t)
predict(model, :survival, ts)          # n × length(ts) matrix
predict(model, :expected, ts)

The default no-arg form (predict(model) or predict(model, :survival)) evaluates each subject at their own observed time $T_i$, matching the convention used by the Cox interface.

Predict on new data

Each prediction also accepts a newdata::DataFrame argument. The fit's stored formula(s) are re-applied to newdata to rebuild the design matrices $X_1$, $X_2$, so newdata must contain every predictor column referenced in the original @formula(...). For models fit with fit(GHM, formula, df) (one formula), the same formula is stored twice and used for both $X_1$ and $X_2$; for fit(GeneralHazard, formula1, formula2, df) the two are stored separately.

predict(model, :survival, newdata, t)        # length-n_new at scalar t
predict(model, :expected, newdata, t)
predict(model, :survival, newdata, ts)       # n_new × length(ts) matrix
predict(model, :expected, newdata, ts)

Newdata predict requires an explicit time argument — there is no "own time" default for arbitrary new subjects. Models built directly via the positional constructor (without the formula1 / formula2 keyword arguments) do not have stored formulas and will error on newdata predict.

SurvivalModels.predict_expectedMethod
predict_expected(m::GeneralHazardModel)
predict_expected(m::GeneralHazardModel, t::Real)
predict_expected(m::GeneralHazardModel, ts::AbstractVector)
predict_expected(m::GeneralHazardModel, newdata::DataFrame, t::Real)
predict_expected(m::GeneralHazardModel, newdata::DataFrame, ts::AbstractVector)

Per-subject cumulative hazard $\Lambda_i(t) = H_0(t\, c_{1i})\, c_{2i}$, where $H_0$ is the cumulative hazard of the baseline distribution and $(c_{1i}, c_{2i})$ are the method-specific time- and hazard-scale multipliers (PH, AFT, AH, GH share the same closed form via the unified $H(t \mid x) = H_0(t\, c_1)\, c_2$ representation).

Output shape:

  • no time argument → length-n vector with each subject evaluated at their own observed time $T_i$;
  • t::Real → length-n vector at the scalar time;
  • ts::AbstractVectorn × length(ts) matrix.

With newdata::DataFrame the design matrices are rebuilt by applying the fit's stored formula(s) — newdata must contain every predictor column referenced in the original @formula(...). Newdata predict requires an explicit time argument (no "own time" default).

source
SurvivalModels.predict_survivalMethod
predict_survival(m::GeneralHazardModel)
predict_survival(m::GeneralHazardModel, t::Real)
predict_survival(m::GeneralHazardModel, ts::AbstractVector)
predict_survival(m::GeneralHazardModel, newdata::DataFrame, t::Real)
predict_survival(m::GeneralHazardModel, newdata::DataFrame, ts::AbstractVector)

Per-subject survival probability $S_i(t) = \exp(-\Lambda_i(t))$ derived from predict_expected. Shapes match predict_expected; newdata variants require an explicit time argument.

source

Simulation

The simulate command (ported from HazReg.jl; formerly simGH, now a deprecated alias) allows one to simulate times to event from the following models:

  • General Hazard (GH) model [5] [6].
  • Accelerated Failure Time (AFT) model [7].
  • Proportional Hazards (PH) model [8].
  • Accelerated Hazards (AH) model [9].
SurvivalModels.simulateFunction
simulate(n, method, baseline, X1, X2, α, β; rng=Random.default_rng())
simulate(n, method, baseline, formula[, formula2], df, α, β; rng=Random.default_rng())
simulate(n, model::GeneralHazardModel; rng=Random.default_rng())

Simulate n times to event from a general-hazard model with hazard structure method (PHMethod()/AFTMethod()/AHMethod()/GHMethod()) and baseline distribution. Each of the n rows of the design (X1/X2, or built from formula/df via modelcols exactly as in fit) yields one event time by inverting that subject's survival Sᵢ(t)=U, i.e. Tᵢ = quantile(baseline, 1 − (1−Uᵢ)^(1/c2ᵢ)) / c1ᵢ. α are the X2 coefficients and β the X1 coefficients (only the ones the structure uses).

The fitted-model form delegates here using the model's own design/coefficients.

References:

source