Illustrative example

In this example, we simulate $n=100$ times to event from the GH, PH, AFT, and AH models with PGW baseline hazards, using the simulate() function. This functionality was ported from HazReg.jl.

PGW-GH model

The fitted model displays the baseline distribution and the two coefficient sets grouped by role (time-scale for X_2, hazard-level for X_1):

using SurvivalModels, Distributions, DataFrames, Random, SurvivalDistributions
using SurvivalModels: simulate

# Simulate design matrices
n = 100
Random.seed!(123)
des = randn(n, 2)
des_t = randn(n, 2)

# True parameters
theta0 = [0.1, 2.0, 5.0]
alpha0 = [0.5, 0.8]
beta0 = [-0.5, 0.75]

# Construct the model directly (no optimization)
model = GeneralHazard(zeros(n), trues(n),
    PowerGeneralizedWeibull(theta0...),
    des, des_t, alpha0, beta0)

# Simulate event times
simdat = simulate(n, model)

# Administrative censoring.
cens = 10
status = simdat .< cens
simdat = min.(simdat, cens)

# Model fit from dataframe interface.
df = DataFrame(time=simdat, status=status, x1=des[:,1], x2=des[:,2], z1=des_t[:,1], z2=des_t[:,2])
model = fit(GeneralHazard{PowerGeneralizedWeibull},
    @formula(Surv(time, status) ~ x1 + x2),
    @formula(Surv(time, status) ~ z1 + z2),
    df)

General hazard model (PowerGeneralizedWeibull baseline)
  n: 100, events: 92
  log-likelihood: -71.241, AIC: 156.48, BIC: 174.72
  baseline: PowerGeneralizedWeibull(0.12554, 1.7538, 3.9287)
  coefficients:
    time-scale:
      z1  0.49282
      z2  1.1708
    hazard-level:
      x1  -0.50864
      x2  0.73921

Comparing the fitted values against the truth:

result = DataFrame(
    Parameter = ["θ₁", "θ₂", "θ₃", "α₁", "α₂","β₁", "β₂"],
    True      = vcat(theta0, alpha0, beta0),
    Fitted    = vcat(params(model.baseline)..., model.α, model.β)
)

Of course, increasing the number of observations would increase the quality of the fitted values. You can also use "subset" models (PH, AH, AFT) through the convenient constructors as follows:

PGW-PH model

model = ProportionalHazard(zeros(n), trues(n),
    PowerGeneralizedWeibull(theta0...),
    des, zeros(n,0),  # X2 is empty for PH
    zeros(0), beta0
)

# Simulate event times and censor them
simdat = simulate(n, model)
cens = 10
status = simdat .< cens
simdat = min.(simdat, cens)


# Build the model and fit it:
df = DataFrame(time=simdat, status=status, x1=des[:,1], x2=des[:,2])
model = fit(ProportionalHazard{PowerGeneralizedWeibull},
    @formula(Surv(time, status) ~ x1 + x2), df)

Proportional hazard model (PowerGeneralizedWeibull baseline)
  n: 100, events: 96
  log-likelihood: -61.941, AIC: 133.88, BIC: 146.91
  baseline: PowerGeneralizedWeibull(0.14499, 1.6597, 3.4394)
  coefficients:
    x1  -0.57227
    x2  0.73552

result = DataFrame(
    Parameter = ["θ₁", "θ₂", "θ₃", "β₁", "β₂"],
    True      = vcat(theta0, beta0),
    Fitted    = vcat(params(model.baseline)..., model.β)
)

PGW-AFT model

# Construct the model directly (no optimization)
model = AcceleratedFaillureTime(
    zeros(n), trues(n), PowerGeneralizedWeibull(theta0...),
    des, zeros(n,0),  # X2 is empty for AFT
    zeros(0), beta0
)

# Simulate event times
simdat = simulate(n, model)

# Censoring
cens = 10
status = simdat .< cens
simdat = min.(simdat, cens)

df = DataFrame(time=simdat, status=status, x1=des[:,1], x2=des[:,2])
model = fit(AcceleratedFaillureTime{PowerGeneralizedWeibull},
    @formula(Surv(time, status) ~ x1 + x2), df)

Accelerated failure time model (PowerGeneralizedWeibull baseline)
  n: 100, events: 97
  log-likelihood: -68.368, AIC: 146.74, BIC: 159.76
  baseline: PowerGeneralizedWeibull(0.091449, 2.2027, 5.8292)
  coefficients:
    x1  -0.54927
    x2  0.73101

result = DataFrame(
    Parameter = ["θ₁", "θ₂", "θ₃", "β₁", "β₂"],
    True      = vcat(theta0, beta0),
    Fitted    = vcat(params(model.baseline)..., model.β)
)

PGW-AH model

# Construct the model directly (no optimization)
model = AcceleratedHazard(zeros(n), trues(n),
    PowerGeneralizedWeibull(theta0...),
    zeros(n,0), des_t,  # X1 is empty for AH
    alpha0, zeros(0)
)

# Simulate event times
simdat = simulate(n, model)
cens = 10
status = simdat .< cens
simdat = min.(simdat, cens)

df = DataFrame(time=simdat, status=status, z1=des_t[:,1], z2=des_t[:,2])
model = fit(AcceleratedHazard{PowerGeneralizedWeibull},
    @formula(Surv(time, status) ~ z1 + z2), df)

Accelerated hazard model (PowerGeneralizedWeibull baseline)
  n: 100, events: 99
  log-likelihood: -53.545, AIC: 117.09, BIC: 130.12
  baseline: PowerGeneralizedWeibull(0.086929, 2.428, 5.7315)
  coefficients:
    z1  0.37483
    z2  0.91613

result = DataFrame(
    Parameter = ["θ₁", "θ₂", "θ₃", "α₁", "α₂"],
    True      = vcat(theta0, alpha0),
    Fitted    = vcat(params(model.baseline)..., model.α)
)