WeibullPH

Parametric proportional hazards with flexible monotone baseline hazard

from _api_doc_utils import *

1 Where it fits

Group: Survival / event-time models

WeibullPH generalizes the exponential PH model to

\[ h_i(t \mid x_i) = \lambda_0 \rho t^{\rho-1} \exp(x_i'\beta), \]

which nests increasing, decreasing, and constant baseline hazards depending on the shape parameter \(\rho\). Like ExponentialPH, it identifies absolute hazard and survival objects.

2 Likelihood and predictions

With \(\lambda=\exp(\alpha)\) and shape \(\rho=\exp(\gamma)\), the model is

\[ h(t\mid x) = \lambda\rho t^{\rho-1}\exp(x'\beta), \qquad H(t\mid x) = \lambda t^\rho\exp(x'\beta). \]

The right-censored log-likelihood is

\[ \ell(\alpha,\beta,\gamma) = \sum_i \left[ d_i\{\alpha+\gamma+x_i'\beta+(\rho-1)\log t_i\} -\exp(\alpha+x_i'\beta+\rho\log t_i) \right]. \]

Newton optimization uses analytic derivatives, a \(10^{-8}\) Hessian ridge, per-coordinate step clipping to \([-2,2]\), and backtracking. Initialization uses the exponential event-rate estimate for \(\alpha\), zero slopes, and \(\rho=1\). Convergence requires either \(\|\nabla\ell\|_\infty<\tau\) or an accepted likelihood change below \(\tau(1+|\ell|)\). Default prediction is \(\exp\{-H(t\mid x)\}\); the linear prediction is the time-specific log hazard.

3 Implementation walkthrough

WeibullPH reuses the exponential model’s native Newton driver but adds the log-shape coordinate analytically.

1. The parameter vector is \((\alpha,\beta',\gamma)'\) with \(\rho=e^\gamma\), which enforces a positive shape without constrained optimization. Initialization uses the exponential event-rate log hazard, zero slopes, and \(\gamma=0\) so \(\rho=1\).

2. For every row the derivative kernel computes \(L_i=\log t_i\), \(z_i=\exp(\alpha+x_i'\beta+\rho L_i)\), and the full censored likelihood contribution. The \((\alpha,\beta)\) score is \((d_i-z_i)(1,x_i')'\); the shape score is

   \[ d_i(1+\rho L_i)-z_i\rho L_i. \]

3. The Hessian’s \((\alpha,\beta)\) block is the negative weighted outer product from the exponential model. Cross-derivatives with \(\gamma\) use \(-z_i(1,x_i')'\rho L_i\), and the shape diagonal adds

   \[ d_i\rho L_i-z_i\{(\rho L_i)^2+\rho L_i\}. \]

   These are evaluated directly rather than by automatic or finite differentiation.

4. The optimizer subtracts \(10^{-8}\) from the log-likelihood Hessian diagonal, explicitly inverts it, clips every coordinate of \(H^{-1}g\) to \([-2,2]\), and tries up to 30 step halvings in \(\theta-H^{-1}g\). It accepts a finite nondecreasing likelihood and otherwise aborts.

5. Maximum-score and relative-objective stopping are checked exactly as in ExponentialPH; the iteration cap is not accepted as convergence. The final observed-information covariance is recomputed without the solve ridge.

6. Prediction reconstructs \(\rho=e^{\hat\gamma}\). The time-specific log hazard includes \(\log\rho+(\rho-1)\log t\), while cumulative hazard uses \(\exp(\alpha+X\beta)t^\rho\). Survival is its negative exponential.

Optimizing log shape avoids boundary handling, but large \(|\gamma|\) or \(|\log t|\) can magnify both \(z_i\) and its shape derivatives. Step clipping is the main safeguard; there is no exponential clipping or rescaling inside the implementation.

4 Inference

The covariance is inverse observed information for the parameter order

\[ (\alpha,\ \beta',\ \gamma)'. \]

It is model-based under a correctly specified Weibull PH likelihood, independent observations, and noninformative right censoring. The summary exposes the full covariance and shape but does not split out standard errors. There is no robust, cluster, bootstrap, or Wald method. Reaching max_iterations raises ValueError and leaves the estimator unfitted. A successful summary exposes converged, iterations, termination_reason, and objective, where objective is the negative log-likelihood.

5 Performance and numerical behavior

Dense Hessian construction costs \(O(np^2)\) and its Newton solve up to \(O(p^3)\) per iteration. Terms involving \(\exp(\gamma)\log t\) make the model more sensitive than ExponentialPH to extreme time scales; rescaling time can materially improve conditioning while changing the scale parameter in the expected way. Exponentials are not clipped. The implementation accepts only positive finite stop times and binary event indicators, with no delayed entry, weights, strata, or time-varying covariates.

6 Python API

Constructor: cm.WeibullPH

Use fit(x, time, event). predict_lin(x, time) returns the log hazard at a specific horizon, while predict_hazard, predict_cumulative_hazard, and predict_survival expose the absolute prediction surface. The default predict(x, time) returns survival probabilities. Times must be positive and finite, events must be binary, and x must be finite.

print(inspect.signature(cm.WeibullPH))

()

cls = cm.WeibullPH
display(HTML(html_table(["Public method"], public_methods(cls))))

Public method
fit(self, /, x, time, event, max_iterations=100, tolerance=1e-08)
predict(self, /, x, time)
predict_cumulative_hazard(self, /, x, time)
predict_hazard(self, /, x, time)
predict_lin(self, /, x, time)
predict_log_hazard(self, /, x, time)
predict_survival(self, /, x, time)
summary(self, /)
survival(self, /, x, time)

7 Minimal example

rng=np.random.default_rng(32)
x=rng.normal(size=(300,2)); rho=1.4; lam=0.03; u=rng.uniform(size=300); t_event=(( -np.log(u)) /(lam*np.exp(x@np.array([0.45,-0.25]))))**(1.0/rho); c=rng.exponential(35,size=300)
time=np.minimum(t_event,c); event=(t_event<=c).astype(float)
model=cm.WeibullPH(); model.fit(x,time,event)
fit=model.summary(); print({key: fit[key] for key in ['converged', 'iterations', 'termination_reason', 'objective']})
print(model.predict_lin(x[:3], time[:3]))
print(model.predict_hazard(x[:3], time[:3]))
print(model.predict_cumulative_hazard(x[:3], time[:3]))
print(model.predict(x[:3], time[:3]))

{'converged': True, 'iterations': 5, 'termination_reason': 'Relative objective tolerance reached', 'objective': 768.4966389774979}
[-1.9818025  -2.48230672 -2.03104686]
[0.13782059 0.08355028 0.1311981 ]
[0.94552302 0.49164741 1.19831413]
[0.38847634 0.61161798 0.30170242]

8 summary() contract

The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.

rng=np.random.default_rng(132); x=rng.normal(size=(140,2)); rho=1.3; lam=0.04; u=rng.uniform(size=140); te=(( -np.log(u)) /(lam*np.exp(x@np.array([0.35,-0.15]))))**(1.0/rho); c=rng.exponential(25,size=140); time=np.minimum(te,c); event=(te<=c).astype(float)
model=cm.WeibullPH(); model.fit(x,time,event)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))

summary() key	shape
log_scale_hazard	()
shape	()
coef	(2,)
hazard_ratio	(2,)
vcov	(4, 4)
log_likelihood	()
converged	()
iterations	()
termination_reason	()
objective	()