mixle.inference.survival module

Time-to-event (survival) estimators and hazard regression.

Survival analysis models the time until an event in the presence of right-censoring (subjects who leave the study before the event – their time is a lower bound, not a missing value). The toolkit here covers the estimators and the regression layer:

  • kaplan_meier() / nelson_aalen() – nonparametric survival and cumulative-hazard curves with Greenwood / Poisson variance and POINTWISE confidence intervals (per time point; no simultaneous band is implemented).

  • cox_ph() – the Cox proportional-hazards regression: how covariates multiply the hazard, estimated from the partial likelihood (Efron or Breslow tie handling), with stratification and time-varying covariates (counting-process start, stop input), Breslow baseline hazard, and the concordance index.

  • discrete_time_hazard() (+ to_person_period()) – discrete-time hazard models fit as a binary GLM on the person-period array (logit or complementary-log-log), supporting offsets and fixed effects through the design matrix.

  • aalen_johansen() – competing-risks cumulative incidence functions (cause-specific).

  • aalen_additive() – Aalen’s additive-hazards regression (cumulative covariate effects).

  • frailty_cox() – shared gamma-frailty Cox for clustered survival (random effect per group), fit by EM.

Event indicators are 1 for an observed event and 0 for right-censoring (for competing risks, an integer cause label with 0 = censored).

Assumptions shared by everything here, stated once: subjects (or, for the frailty model, groups) are INDEPENDENT, and censoring is NON-INFORMATIVE (independent of the event time given the covariates). Informative censoring – subjects leaving because they are about to fail – biases every estimator in this module, and nothing in the data can detect it; it is a design claim the caller must own.

kaplan_meier(time, event=None, *, ci_level=0.95)[source]

Kaplan–Meier product-limit estimate of the survival function S(t).

Parameters:

  • time (ndarray) – (n,) observed times (event or censoring).

  • event (ndarray | None) – (n,) 1 = event, 0 = right-censored (defaults to all events).

  • ci_level (float) – confidence level for the POINTWISE log–log interval at each time point.

Returns:

{'time', 'survival', 'se', 'ci_low', 'ci_high', 'at_risk', 'n_events', 'median'}.

Return type:

dict[str, ndarray]

ci_low/ci_high are PER-TIME-POINT (pointwise) intervals, not a simultaneous confidence band: each is a MARGINAL, asymptotic (Greenwood log–log) interval for S(t) at its own t with probability ci_level under independent censoring, but the probability that the true curve stays inside the envelope at EVERY time point is far lower – measured simultaneous coverage of the pointwise 95% envelope is 0.60 / 0.51 / 0.43 at n = 50 / 100 / 200 (it FALLS with n, because more event times mean more chances to escape). Reading the plotted envelope as “the curve lives in here” is a simultaneous claim these intervals do not support; a Hall–Wellner band is the standard tool for that and is not implemented here.

nelson_aalen(time, event=None)[source]

Nelson–Aalen estimate of the cumulative hazard H(t) = sum d_i / Y_i.

Returns:

{'time', 'cumhaz', 'se'} with the Poisson-type standard error of the cumulative hazard.

Parameters:
Return type:

dict[str, ndarray]

class CoxResult(coef, se, cov, loglik, baseline_time, baseline_cumhaz, concordance, n_iter, converged, baseline_by_stratum=<factory>)[source]

Bases: object

Fitted Cox proportional-hazards model.

Parameters:
coef

(p,) log-hazard-ratio coefficients.

Type:

numpy.ndarray

se

(p,) standard errors (inverse observed information).

Type:

numpy.ndarray

cov

(p, p) covariance.

Type:

numpy.ndarray

loglik

maximised partial log-likelihood.

Type:

float

baseline_time / baseline_cumhaz

Breslow baseline cumulative hazard for an unstratified fit. These arrays are empty for a multi-stratum fit.

baseline_by_stratum

mapping from each stratum label to its own time and cumhaz arrays.

Type:

dict[object, dict[str, numpy.ndarray]]

concordance

Harrell’s C-index, evaluated only within strata. This is None for counting-process input unless subject identifiers are supplied.

Type:

float | None

n_iter

Newton iterations.

Type:

int

converged

whether Newton’s method reached tol with a finite, well-identified fit. False means either it hit max_iter without reaching tol, a step produced a non-finite iterate (kept the last finite one instead), or the final standard errors are enormous – the classic symptom of (quasi-)complete separation, where the partial likelihood has no finite maximizer and Newton’s steps can shrink below tol while the coefficients themselves keep drifting toward infinity.

Type:

bool

hazard_ratios()[source]

Return exponentiated Cox coefficients.

Return type:

ndarray

z_values()[source]

Return Wald z statistics for Cox coefficients.

Return type:

ndarray

p_values()[source]

Return two-sided normal-approximation p-values for Cox coefficients.

Return type:

ndarray

cox_ph(x, time, event, *, start=None, strata=None, subject=None, ties='efron', max_iter=50, tol=1e-9)[source]

Cox proportional-hazards regression by Newton–Raphson on the partial likelihood.

The hazard is h(t | x) = h0(t) exp(x' beta); only the ordering of event times enters, so the baseline h0 is left unspecified (semi-parametric). Time-varying covariates are supported through the counting-process form: pass start so each row is an at-risk interval (start, stop] (a subject contributes several rows), and the risk set at an event time is every interval covering it.

Parameters:

  • x (ndarray) – (n, p) covariates (no intercept – it is absorbed into the baseline).

  • time (ndarray) – (n,) event/censoring times (the interval stop times).

  • event (ndarray) – (n,) 1 = event, 0 = censored.

  • start (ndarray | None) – optional (n,) interval start times for time-varying covariates / left truncation.

  • strata (ndarray | None) – optional (n,) labels; each stratum gets its own baseline hazard (coefficients shared).

  • subject (ndarray | None) – optional (n,) subject identifiers. Required to compute concordance for counting-process input with multiple rows per subject.

  • ties (str) – "efron" (default, more accurate) or "breslow" tie handling.

  • max_iter (int) – Newton controls.

  • tol (float) – Newton controls.

Returns:

A CoxResult.

Return type:

CoxResult

What the standard errors assume – read this before using counting-process input for clustered or recurrent events. The SEs are the inverse observed information of the partial likelihood, which is valid when SUBJECTS are independent. The (start, stop] interface happily accepts multiple rows per subject (that is its purpose: time-varying covariates for one subject) and also LOOKS like it accepts recurrent events or cluster members as rows – but rows from the same subject/cluster are dependent, the information-based SEs do not know that, and they understate: on a simulated clustered design the mean reported SE was 0.68 of the true sampling SD. For clustered data use frailty_cox() (models the dependence) or bootstrap over SUBJECTS/CLUSTERS – resampling rows would repeat the same mistake. Proportional hazards itself (covariate effects constant over time) is assumed, not checked: no diagnostic ships in this module, so inspect Schoenfeld-style residuals or stratified log(-log) survival plots externally before reading coef as a time-constant log-hazard ratio.

to_person_period(time, event, covariates=None)[source]

Expand right-censored durations into a person-period (long) array for discrete-time models.

Each subject contributes one row per discrete period they were at risk; the binary outcome is 1 in the period the event occurred and 0 otherwise. Integer time is the number of periods observed.

Returns:

{'period', 'outcome', 'subject', 'covariates'} (covariates repeated per period if given).

Parameters:
Return type:

dict[str, ndarray]

discrete_time_hazard(x, outcome, *, link='cloglog', offset=None)[source]

Discrete-time hazard model: a binary GLM on the person-period array.

Fit on the long-format data from to_person_period() (the design x typically holds period indicators / a time trend plus covariates). cloglog gives the grouped-proportional-hazards (interval-censored Cox) interpretation; logit gives the proportional-odds hazard.

Returns:

a mixle.inference.glm.GLMResult (binomial family with the chosen link).

Parameters:
aalen_johansen(time, event, *, causes=None)[source]

Aalen–Johansen cumulative incidence functions for competing risks.

With several mutually exclusive event types, the cause-specific CIF F_k(t) is the probability of failing from cause k by time t accounting for the competing causes (it is not 1 - KM on the cause, which overstates incidence).

Parameters:

  • time (ndarray) – (n,) event/censoring times.

  • event (ndarray) – (n,) integer cause label, 0 = censored, 1..K = causes.

  • causes (ndarray | None) – optional explicit list of cause labels; inferred from event if None.

Returns:

{'time', 'cif': {cause: array}, 'overall_survival'}.

Return type:

dict

aalen_additive(x, time, event, *, intercept=True)[source]

Aalen’s additive-hazards regression: cumulative regression functions B(t).

Models h(t | x) = b0(t) + sum_j x_j b_j(t) with time-varying additive effects. At each event time the increment dB is the least-squares solution over the risk set; the cumulative B(t) (returned) has interpretable slopes – a rising B_j means covariate j adds hazard.

The estimator is only defined while the risk-set design has full column rank, and it degrades BEFORE it breaks: as the risk set shrinks toward p subjects the least-squares increments are dominated by noise, and the accumulated tail of B(t) is numerical garbage that LOOKS like a dramatic late effect (measured: max abs(B) growing 1.4 -> 19 over the last few event times of a null design). Estimation therefore STOPS at the first event time whose risk-set design is rank-deficient or numerically ill-conditioned – the standard Aalen practice – instead of silently substituting a minimum-norm solution and accumulating it. truncated_at reports where (None if the full curve was estimable); the curve up to that point is unaffected by the truncation.

Returns:

{'time', 'cum_coef', 'truncated_at'} where cum_coef is (n_estimated_times, p[+1]) cumulative coefficients (the first column is the baseline when intercept is True) and truncated_at is the first event time excluded because the risk-set design lost full rank (None when no truncation occurred).

Parameters:
Return type:

dict

class FrailtyCoxResult(coef, se, theta, frailties, groups, frailty_variance=<factory>, frailty_log_mean=<factory>, n_iter=0, converged=False, ties='breslow', se_method='complete-data')[source]

Bases: object

Shared gamma-frailty Cox result.

Parameters:
coef / se

fixed-effect log-hazard-ratios and standard errors (see se_method for what the SEs account for – the default plug-in SEs understate).

theta

estimated frailty variance. A POINT estimate with no uncertainty attached: no SE or test accompanies it, the null theta = 0 sits on the parameter boundary (a likelihood ratio there is a 50:50 chi2_0/chi2_1 mixture, not implemented), and the EM estimate never actually reaches 0 under a true theta of 0 (measured mean 0.046 on null simulations) – so a small positive theta is what NO clustering looks like, not weak evidence of some.

Type:

float

frailties

posterior mean random effect per group.

Type:

numpy.ndarray

frailty_variance

posterior variance of each group random effect.

Type:

numpy.ndarray

frailty_log_mean

posterior expectation of log(w_g).

Type:

numpy.ndarray

groups

group labels aligned to frailties.

Type:

numpy.ndarray

n_iter

EM iterations.

Type:

int

converged

whether the EM convergence criterion was met.

Type:

bool

se_method

"complete-data" (plug-in information, understates) or "jackknife" (delete-one-group, accounts for frailty/theta estimation).

Type:

str

frailty_cox(x, time, event, groups, *, max_iter=50, tol=1e-5, ties='breslow', se_method='complete-data')[source]

Shared gamma-frailty Cox model for clustered survival, by EM.

Subjects in the same group share an unobserved frailty w_g ~ Gamma(1/theta, 1/theta) (mean 1, variance theta) that multiplies the hazard, capturing within-group correlation. The E-step retains E[w_g], Var(w_g), and E[log(w_g)] under the conjugate gamma posterior. The coefficient/baseline M-step uses E[w_g] in the integrated-hazard term, while the dispersion M-step maximises the expected gamma-prior log likelihood using both E[w_g] and E[log(w_g)].

On theta: it is returned as a point estimate ONLY – see the honesty notes on FrailtyCoxResult before reading a small value as evidence for (or against) clustering.

Parameters:

  • se_method (str) – how the fixed-effect SEs are computed. "complete-data" (default) inverts the partial-likelihood information with the posterior-mean frailties and theta PLUGGED IN as if known – it prices in none of their estimation uncertainty and understates: measured against the true sampling SD at theta = 1, se/SD was 0.79-0.88 and Wald-95 coverage 0.87-0.91 (clean at theta = 0). "jackknife" deletes one GROUP at a time and refits, pricing in frailty and theta estimation; it is honest but costs G extra EM fits, and its own validity is asymptotic in the NUMBER OF GROUPS (a handful of groups gives a noisy SE from few leave-outs).

  • x (ndarray)

  • time (ndarray)

  • event (ndarray)

  • groups (ndarray)

  • max_iter (int)

  • tol (float)

  • ties (str)

Returns:

A FrailtyCoxResult (se_method records which SEs you got).

Return type:

FrailtyCoxResult