----------------------------------------------------------------------
This is the API documentation for the greenwood library.
----------------------------------------------------------------------


## The response

The Surv object, the spine of every analysis.


Surv(type: 'CensoringType', stop: 'Array', status: 'Array', start: 'Array | None' = None, lower: 'Array | None' = None, states: 'tuple[str, ...] | None' = None, weights: 'Array | None' = None) -> None

A validated time-to-event response for survival analysis.

`Surv` represents the outcome in survival models: a time at which each subject either
experienced an event (observed) or was censored (did not experience the event during
follow-up). `Surv` supports multiple censoring types:

- **Right-censored** (most common): The event time is at or after the recorded time.
  Use `Surv.right(time, event)`.
- **Left-censored**: The event time is before the recorded time.
  Use `Surv.left(time, event)`.
- **Counting-process** (left truncation, time-varying): Each subject enters the risk set
  at `start` and exits at `stop`. Use `Surv.counting(start, stop, event)`.
- **Interval-censored**: The event occurred within a time interval `[lower, upper)`.
  Use `Surv.interval(lower, upper)`.
- **Multi-state / competing risks**: Multiple mutually exclusive events.
  Use `Surv.multistate(time, event, states)`.

**Use the class methods** (`right`, `left`, `counting`, `interval`, `multistate`)
to construct `Surv` objects. They validate your input and set the censoring type
appropriately. As such, direct instantiation is not recommended.

Attributes
----------
type
    The `CensoringType` enum indicating the censoring mechanism.
stop
    Exit time (for interval censoring, the upper bound).
status
    Integer event code per observation: 0 = censored, 1+ = event code (for multi-state,
    codes >= 1 index into `states`).
start
    Entry time for the counting-process form (left truncation); `None` otherwise.
lower
    Lower bound for interval censoring; `None` otherwise.
states
    Event-state labels for multi-state/competing-risks endpoints; `None` for the
    single-event case.
weights
    Optional case weights (strictly positive); `None` if no weights provided.

Examples
--------
Here's an example of direct instantiation of `Surv`:

```{python}
import greenwood as gw
import numpy as np

y = gw.Surv(
    type=gw.CensoringType.RIGHT,
    stop=np.array([5, 6, 4, 9]),
    status=np.array([1, 0, 1, 0])
)
y
```

While this is fine, the preferred approach is to use the class method constructors for each
censoring type. They handle validation and conversion automatically.

Right-censored (the most common case): each subject has an exit time and an event
indicator.

```{python}
y = gw.Surv.right(time=[5, 6, 4, 9], event=[1, 0, 1, 0])
y
```

Counting-process form with left truncation (late entry):

```{python}
y = gw.Surv.counting(start=[0, 2, 1], stop=[5, 6, 4], event=[1, 0, 1])
y
```

Interval-censored (event known to occur in a time window):

```{python}
y = gw.Surv.interval(lower=[1, 3], upper=[3, 8])
y
```

Multi-state (competing risks, multiple mutually exclusive events):

```{python}
y = gw.Surv.multistate(time=[5, 6, 4], event=[1, 0, 2], states=("pcm", "death"))
y
```

CensoringType(*values)

The censoring flavor of a `Surv` response.

Examples
--------
The members correspond to the `Surv` constructors. A constructed response reports its
flavor through `Surv(...).type`, which is one of these values.

```{python}
from greenwood import CensoringType

list(CensoringType)
```


## Non-parametric estimators

Kaplan-Meier survival and Nelson-Aalen cumulative hazard.


KaplanMeier(*, conf_type: 'str' = 'log', conf_level: 'float' = 0.95) -> 'None'

Kaplan-Meier product-limit estimator of the survival function.

The Kaplan-Meier estimator is a non-parametric method to estimate the survival function
from right-censored data. It computes the survival probability at each observed event time
as the product of conditional survival probabilities, accounting for subjects still at risk.
This is the most widely used method for survival analysis and is the starting point for
comparing survival between groups or assessing model fit.

To use this estimator, call `fit()` with a right-censored `Surv` response (built with
`Surv.right()`). The estimator computes survival probabilities, standard errors, and
confidence intervals at each unique event time. Results can be accessed as aligned
arrays, exported to pandas/polars/pyarrow DataFrames, or queried through methods like
`median()`, `quantile()`, and `predict()`.

The implementation uses the product-limit formula

$$
S(t) = \prod_{t_i \le t} \frac{n_i - d_i}{n_i}
$$

where $n_i$ is the number at risk and $d_i$ is the number of events at time $t_i$.
Variance uses Greenwood's formula, and confidence intervals can be constructed on the
log, log-log, or identity scale.

Parameters
----------
conf_type
    Confidence-interval transform: `"log"` (default, as in R's `survfit`), `"plain"`,
    or `"log-log"`.
conf_level
    Confidence level for the interval (default 0.95).

Returns
-------
Fitted estimator
    Call `fit()` to produce a fitted estimator with cached results (`time_`, `surv_`,
    `std_error_`, `conf_low_`, `conf_high_`, `n_risk_`, `n_event_`, `n_censor_`),
    accessible as aligned arrays or exported to DataFrames.

Details
-------
Call `fit` with a `Surv` response. Results are exposed as aligned arrays (`time_`,
`survival_`, `std_error_`, `conf_low_`, `conf_high_`, `strata_`), as tidy frames via
`to_frame()` (optionally `format=`), and through `median`, `quantile`, and `predict`.

Examples
--------
Build a `Surv` response from the bundled `lung` dataset and fit the estimator. Printing
the fitted object reports the median survival and its confidence interval.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)
km
```

The full step function, one row per event time, is available with `to_frame`; pass
`format=` to choose the backend (here, Polars):

```{python}
km.to_frame(format="polars")
```

NelsonAalen(*, conf_type: 'str' = 'log', conf_level: 'float' = 0.95) -> 'None'

Nelson-Aalen estimator of the cumulative hazard.

The Nelson-Aalen estimator provides a non-parametric estimate of the cumulative hazard
function, which represents the total "accumulated risk" up to a given time. Unlike the
Kaplan-Meier estimator which models survival directly, this approach models the force of
mortality. The cumulative hazard at each event time is computed as a running sum of the
ratio of events to subjects at risk:

$$
H(t) = \sum_{t_i \le t} \frac{d_i}{n_i}
$$

This estimator is useful when you want to examine the hazard directly rather than survival
probabilities, and is often used as the basis for other analyses. You can convert the
cumulative hazard to a survival estimate via $S(t) = \exp(-H(t))$, though the Kaplan-Meier
estimator is typically preferred for direct survival estimation. Call `fit()` with a
right-censored `Surv` response to compute cumulative hazard at each event time.

The variance of the cumulative hazard estimate uses Aalen's formula:

$$
\mathrm{Var}(H(t)) = \sum_{t_i \le t} \frac{d_i}{n_i^2}
$$

Confidence intervals can be constructed on the plain or log scale, with the log scale
providing better coverage in the tails.

Parameters
----------
conf_type
    Confidence-interval transform: `"plain"` (default for Nelson-Aalen) or `"log"`.
conf_level
    Confidence level for the interval (default 0.95).

Returns
-------
Fitted estimator
    Call `fit()` to produce a fitted estimator with cached results (`time_`,
    `cumulative_hazard_`, `std_error_`, `conf_low_`, `conf_high_`, `n_risk_`, `n_event_`,
    `n_censor_`), accessible as aligned arrays or exported to DataFrames.

Details
-------
Call `fit()` with a `Surv` response. Results are exposed as aligned arrays, as tidy
frames via `to_frame()` (optionally `format=`), and through the `predict()`,
`quantile()`, and other methods.

Examples
--------
Build a `Surv` response from the bundled `lung` dataset and fit the estimator. Printing
the fitted object reports the counts and the maximum cumulative hazard reached.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
na = gw.NelsonAalen().fit(y)
na
```


## Regression

Cox proportional hazards and parametric AFT models.


CoxPH(*, ties: 'str' = 'efron', conf_level: 'float' = 0.95) -> 'None'

Cox proportional hazards model.

The Cox proportional hazards model is the most widely used regression method for survival
data. It models the hazard (instantaneous risk of an event) as a multiplicative function
of covariates: $h(t \mid x) = h_0(t) \exp(\beta^\top x)$. The model is semi-parametric:
the baseline hazard $h_0(t)$ is left unspecified (estimated non-parametrically), while
covariate effects are estimated parametrically through the log-hazard-ratio coefficients
$\beta$.

To use this model, call `fit()` with a right-censored or counting-process `Surv` response
and a design matrix of covariates (2-D array or DataFrame). The model automatically handles
stratification (via `by=` in fit), tied event times (via configurable tie-handling methods),
and can compute predictions, baseline hazards, and diagnostic residuals. Results include
coefficient estimates with confidence intervals, hazard ratios, standard errors, and
global significance tests.

The implementation uses maximum partial likelihood to estimate coefficients. Variance
estimates use the observed information matrix (Hessian). The model assumes proportional
hazards: the ratio of hazards between two subjects remains constant over time. This can
be checked using the `cox_zph()` method for formal tests or diagnostic plots.

Parameters
----------
ties
    Tie-handling method: `"efron"` (default, as in R) or `"breslow"`.
conf_level
    Confidence level for coefficient and hazard-ratio intervals (default 0.95).

Returns
-------
Fitted estimator
    Call `fit()` to produce a fitted estimator with cached results (`coef_`,
    `hazard_ratio_`, `std_error_`, `z_`, `p_value_`, `conf_low_`, `conf_high_`,
    `concordance_`, `lr_stat_`, `df_`), accessible as arrays or exported to DataFrames.

Details
-------
Call `fit(surv, covariates)` with a `Surv` response and a design (a 2-D array or a
dataframe of covariates). Rows with missing values are dropped (complete-case, as in R's
default `na.omit`). Results are exposed as arrays (`coef_`, `std_error_`, `hazard_ratio_`,
…) and as tidy frames via `to_frame()` (optionally `format=`) and `greenwood.tidy`.

Examples
--------
Build a `Surv` response from the bundled `lung` dataset and fit the model on `age` and
`sex`. Printing the fitted object reports the coefficient table and global tests in the
style of R's `summary.coxph`.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
cox
```

Hazard ratios (and their confidence limits) come from `tidy` with `exponentiate=True`;
pass `format=` to choose the backend (here, Polars):

```{python}
gw.tidy(cox, exponentiate=True, format="polars")
```

CoxNet(penalizer: 'float' = 0.1, l1_ratio: 'float' = 0.5, *, standardize: 'bool' = True, max_iter: 'int' = 1000, tol: 'float' = 1e-07) -> 'None'

Elastic-net penalized Cox proportional hazards model.

When the number of covariates is large relative to the sample size, unpenalized Cox models
may overfit or fail to converge due to multicollinearity. CoxNet addresses this by adding
a penalty term to the partial likelihood, which shrinks coefficients toward zero and can
perform automatic variable selection. The elastic-net penalty combines $L_1$ (lasso) and
$L_2$ (ridge) penalties: $\lambda(\alpha\|\beta\|_1 + \tfrac{1-\alpha}{2}\|\beta\|_2^2)$,
where the mixing parameter $\alpha$ controls the trade-off between sparsity and smoothness.

Fit the model with `fit()` supplying a right-censored or counting-process `Surv` response
and a design matrix of covariates. The algorithm uses FISTA (Fast Iterative Shrinkage-
Thresholding Algorithm) to optimize the penalized partial likelihood. By default, covariates
are standardized before penalizing (for fair comparison of penalties across features), but
coefficients are returned on the original scale. Ridge (`l1_ratio=0`) encourages small,
spread-out coefficients; lasso (`l1_ratio=1`) drives some coefficients exactly to zero.

The implementation follows the glmnet model for elastic-net regularization, using coordinate
descent-like optimization with soft-thresholding. Results include penalized coefficients,
standard errors, and indicators of which features were selected (non-zero coefficients).
Unlike unpenalized Cox, this model does not compute hazard ratios or perform formal
hypothesis tests on individual coefficients.

Parameters
----------
penalizer
    Overall penalty strength (`lambda`). `0` recovers the unpenalized Breslow Cox fit.
l1_ratio
    Elastic-net mixing in `[0, 1]`: `1` is lasso, `0` is ridge.
standardize
    Standardize covariates to unit variance before penalizing (default `True`, as in
    glmnet). Coefficients are returned on the original scale.
max_iter, tol
    Maximum FISTA iterations and the relative-change convergence tolerance.

Returns
-------
Fitted estimator
    Call `fit()` to produce a fitted estimator with cached results (`coef_`,
    `std_error_`, `n_features_in_`, `feature_names_in_`), accessible as arrays or
    exported to DataFrames.

Details
-------
Call `fit(surv, covariates)` with a right-censored or counting-process `Surv` response.
`covariates` may be a dataframe, a 2-D array, or a formula string with `data`. Stratified
penalized fits are not supported.

Examples
--------
Build a `Surv` response from the bundled `lung` dataset and fit a lasso (`l1_ratio=1.0`)
elastic-net Cox model over several covariates. The `ph.ecog`, `ph.karno`, and `wt.loss`
columns have missing values, which `CoxNet` drops automatically. Printing the fitted object
shows the penalized coefficients and how many were driven to zero.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cols = ["age", "sex", "ph.ecog", "ph.karno", "wt.loss"]
coxnet = gw.CoxNet(penalizer=0.05, l1_ratio=1.0).fit(y, lung[cols])
coxnet
```

ZPHResult(transform: 'str', per_term: 'dict[str, dict[str, float]]', global_test: 'dict[str, float]') -> None

Proportional-hazards test results (Grambsch-Therneau).

A key assumption of the Cox proportional hazards model is that the hazard ratio between
any two subjects is constant over time (hence "proportional"). When this assumption is
violated (for example, if a treatment effect diminishes over time) the Cox model may produce
biased estimates. The Grambsch-Therneau proportional hazards test checks this assumption
by testing whether scaled residuals are correlated with time.

`ZPHResult` holds the test results obtained from a fitted Cox model's `cox_zph()` method.
It provides both per-term tests (one for each covariate) and a global test (jointly across
all terms). Each test includes a chi-squared test statistic, degrees of freedom, and
p-value. Results can be printed, accessed via dictionary keys, or exported to pandas/polars/
pyarrow DataFrames for further analysis or visualization.

The test uses scaled Schoenfeld residuals, which have a known asymptotic distribution under
the proportional hazards assumption. Large chi-squared values or small p-values (typically
p < 0.05) suggest violation of the assumption. When the assumption is violated, stratified
analysis or time-dependent covariate models may be more appropriate.

Attributes
----------
transform
    The transformation applied to time when computing the test (e.g., identity, log, rank).
per_term
    Dictionary mapping each covariate name to `{chisq, df, p_value}` dict.
global_test
    Dictionary with `{chisq, df, p_value}` for the joint test across all terms.

Examples
--------
A `ZPHResult` comes from a fitted model's `cox_zph` method. Fit a Cox model to the
bundled `lung` dataset, run the proportional-hazards test, and print the result:

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
zph = cox.cox_zph()
zph
```

AFT(dist: 'str' = 'weibull', *, conf_level: 'float' = 0.95) -> 'None'

Parametric accelerated failure time model.

While the Cox proportional hazards model leaves the baseline hazard unspecified, AFT models
assume a fully parametric distribution for survival times and model how covariates
accelerate or decelerate the "clock" of failure. Specifically, $\log(T) = \mu + \beta^\top x
+ \sigma\varepsilon$, where $T$ is survival time, $\beta$ are log-time-scale coefficients,
$\sigma$ is a scale parameter, and $\varepsilon$ follows a specified error distribution
(e.g., extreme-value, logistic, normal). This means a unit increase in covariate $x$
multiplies survival time by $\exp(\beta)$.

AFT models are useful when you want explicit, interpretable survival time predictions or when
the parametric assumptions are reasonable. Unlike Cox models, they require choosing a
distributional family (Weibull, exponential, lognormal, or loglogistic). Call `fit()` with
a right-censored `Surv` response and a design matrix. The model automatically adds an
intercept and estimates coefficients (on the log-time scale), the scale parameter, and
standard errors via maximum likelihood.

The implementation uses numerical optimization (typically Newton-Raphson) to maximize the
likelihood. Coefficients on the log-time scale can be exponentiated to obtain time-
acceleration ratios: $\exp(\beta)$ is the multiplicative effect on median or mean survival.
The model also supports prediction of survival probabilities and quantiles at future times
given covariate values.

Parameters
----------
dist
    Error distribution: `"weibull"` (default), `"exponential"`, `"lognormal"`, or
    `"loglogistic"`.
conf_level
    Confidence level for coefficient intervals (default is `0.95`).

Returns
-------
Fitted estimator
    Call `fit()` to produce a fitted estimator with cached results (`coef_`,
    `scale_`, `std_error_`, `z_`, `p_value_`, `conf_low_`, `conf_high_`, `loglik_`,
    `aic_`, `bic_`), accessible as arrays or exported to DataFrames.

Details
-------
Call `fit(surv, covariates)` with a right-censored `Surv` response and a covariate design
(a 2-D array or a dataframe). An intercept is added automatically; rows with missing
covariates are dropped. Results are exposed as arrays (`coef_`, `scale_`, `std_error_`,
`z_`, `p_value_`) and as tidy frames via `to_frame()` (optionally `format=`) and
`greenwood.tidy`.

Examples
--------
Build a `Surv` response from the bundled `lung` dataset and fit a Weibull AFT model with
`age` and `sex` as covariates. Printing the fitted object reports the coefficients (on the
log-time scale), the scale, and the log-likelihood.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
aft = gw.AFT("weibull").fit(y, lung[["age", "sex"]])
aft
```

RoystonParmar(df: 'int' = 3, *, conf_level: 'float' = 0.95) -> 'None'

Royston-Parmar flexible parametric survival model (proportional hazards scale).

The Royston-Parmar model offers a middle ground between rigid parametric models (like
Weibull) and fully non-parametric methods (like Kaplan-Meier). It models the log baseline
cumulative hazard as a smooth spline function on the log-time scale, combined with
proportional-hazards covariate effects. This allows flexible baseline shapes while
maintaining interpretable proportional-hazards covariate coefficients. It's a key advantage
over fully parametric AFT models.

The model uses restricted cubic splines with a fixed number of degrees of freedom (controlled
by knots placed at quantiles of event times). A low df value (e.g., df=1) approaches a
Weibull fit; higher df values (e.g., df=3 or 4) provide greater flexibility. Call `fit()`
with a right-censored `Surv` response and a design matrix. The model reports spline and
covariate coefficients, fitted knot locations, log-likelihood, and supports predictions of
survival at specified times and covariate values.

The implementation uses maximum likelihood estimation with constraints that ensure monotone
increasing log cumulative hazard (valid hazard functions). The flexible baseline makes this
model useful when baseline hazard shape is unknown but important, yet you want interpretable
proportional-hazards effects of covariates. Results can be exported to tidy DataFrames or
accessed as coefficient arrays.

Parameters
----------
df
    Spline degrees of freedom: the number of spline terms beyond the intercept, equal to
    one more than the number of internal knots. `df=1` is a Weibull model; `df=3`
    (two internal knots) is a common flexible default.
conf_level
    Confidence level for coefficient intervals (default 0.95).

Returns
-------
Fitted estimator
    Call `fit()` to produce a fitted estimator with cached results (`coef_`,
    `std_error_`, `z_`, `p_value_`, `conf_low_`, `conf_high_`, `knots_`, `loglik_`,
    `aic_`, `bic_`), accessible as arrays or exported to DataFrames.

Details
-------
Call `fit(surv, covariates)` with a right-censored `Surv` response and a covariate design
(a dataframe, a 2-D array, or a formula string with `data`). Results are exposed as arrays
(`coef_`, `std_error_`, ...), the fitted `knots_`, and tidy frames via `to_frame()`
(optionally `format=`).

Examples
--------
Build a `Surv` response from the bundled `lung` dataset and fit a flexible model with three
spline degrees of freedom and `age` and `sex` as covariates. Printing the fitted object
reports the spline and covariate coefficients and the log-likelihood.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
rp = gw.RoystonParmar(df=3).fit(y, lung[["age", "sex"]])
rp
```


## Competing risks & multi-state

Cumulative incidence, the Fine-Gray model, and multi-state transition probabilities.


AalenJohansen(*, conf_level: 'float' = 0.95) -> 'None'

Aalen-Johansen estimator of cumulative incidence functions for competing risks.

In survival analysis, competing risks occur when subjects can experience multiple types
of events (e.g., progression to malignancy vs. death from other causes), and experiencing
one event precludes the others. The Aalen-Johansen estimator extends the Kaplan-Meier
approach to this setting by estimating the cumulative incidence function (CIF) for each
cause: the probability of experiencing that specific cause by time t, accounting for
competition from other causes.

Unlike naive estimates that ignore censoring or competing events, the Aalen-Johansen CIF
correctly accounts for both. It is computed using transition probabilities between states
via generalized Kaplan-Meier estimates. Call `fit()` with a multi-state `Surv` response
(built with `Surv.multistate()`) to obtain estimates for each competing cause. Results
are returned as tidy DataFrames with one row per combination of stratum, cause, and time.

The estimator uses the formula $\mathrm{CIF}_j(t) = \sum_{s \le t} \hat{S}(s^-) P_{0j}(s)$,
where $\hat{S}(s^-)$ is the probability of being event-free before time $s$, and
$P_{0j}(s)$ is the transition probability from censoring to cause $j$. Confidence intervals
use the Greenwood-style variance estimator on the complementary log-log scale for improved
coverage.

Parameters
----------
conf_level
    Confidence level for the (Wald) confidence intervals (default 0.95).

Returns
-------
Fitted estimator
    Call `fit()` to produce a fitted estimator with cached results (`states_`, and
    internal transition matrices), accessible as tidy DataFrames.

Details
-------
Call `fit(surv, by=...)` with a multi-state `Surv` response (built with
`Surv.multistate`, where `event` codes are 0 for censoring and `1..K` for the competing
causes). Results are tidy frames via `to_frame()` (optionally `format=`) with one row
per stratum, cause, and time.

Examples
--------
The bundled `mgus2` dataset follows monoclonal-gammopathy patients who may progress to
plasma-cell malignancy (`"pcm"`) or die first, a competing-risks setup. Build the
competing-risks response by combining the progression and death indicators into a single
cause code (0 censored, 1 progression, 2 death), then fit the estimator. Printing the
fitted object reports the final cumulative incidence for each cause.

```{python}
import greenwood as gw

mg = gw.load_dataset("mgus2", backend="pandas")
etime = mg["ptime"].where(mg["pstat"] == 1, mg["futime"])
cause = mg["pstat"].where(mg["pstat"] == 1, 2 * mg["death"])
cr = gw.Surv.multistate(etime, event=cause, states=("pcm", "death"))
aj = gw.AalenJohansen().fit(cr)
aj
```

FineGray(cause: 'Any', *, conf_level: 'float' = 0.95) -> 'None'

Fine-Gray subdistribution hazard model for a competing-risks endpoint.

The Fine-Gray model extends Cox regression to competing-risks settings where multiple
event types (e.g., disease progression vs. death) can occur and only the first is observed.
Rather than modeling the cause-specific hazard (as in standard Cox), it models the
subdistribution hazard: the rate at which subjects experience the target cause as if
competing events did not occur. This allows inference directly on the cumulative incidence
function, making it ideal for policy-relevant questions like "what is the effect of
treatment on my probability of experiencing event A?"

Technically, the Fine-Gray model uses a weighted Cox-like approach: subjects who experience
a competing event remain in the risk set but with decreasing inverse-probability-of-censoring
weights, reflecting their reduced ability to contribute information about the target cause.
Call `fit()` with a multi-state `Surv` response (built with `Surv.multistate()`) and specify
the target cause of interest. Coefficients, hazard ratios, and standard errors are computed
via weighted partial likelihood, with robust (clustered) standard errors accounting for the
weighting scheme.

The implementation automatically computes event weights and handles censoring. Standard errors
use the Lin-Wei robust (sandwich) estimator, validated against R's survival package. Unlike
the Aalen-Johansen non-parametric approach, Fine-Gray provides covariate-adjusted estimates
and supports prediction of cumulative incidence at specified times.

Parameters
----------
cause
    The target cause-of-interest label from the multi-state `Surv` response.
conf_level
    Confidence level for coefficient intervals (default 0.95).

Returns
-------
Fitted estimator
    Call `fit()` to produce a fitted estimator with cached results (`coef_`,
    `hazard_ratio_`, `std_error_`, `z_`, `p_value_`, `conf_low_`, `conf_high_`),
    accessible as arrays or exported to DataFrames.

Details
-------
Call `fit(surv, covariates)` with a multi-state `Surv` response (built with
`Surv.multistate()`) and specify the target `cause`. The model uses weighted Cox-like
optimization with robust standard errors. Results can be tidy frames via `to_frame()`
(optionally `format=`).

Examples
--------
Using the same competing-risks setup as `AalenJohansen`, model the cumulative incidence of
plasma-cell malignancy (`"pcm"`) as a function of age and sex. The `age` and `sex` columns
of the `mgus2` frame form the covariate design. Printing the fitted object reports the
subdistribution-hazard coefficient table.

```{python}
import greenwood as gw

mg = gw.load_dataset("mgus2", backend="pandas")
etime = mg["ptime"].where(mg["pstat"] == 1, mg["futime"])
cause = mg["pstat"].where(mg["pstat"] == 1, 2 * mg["death"])
cr = gw.Surv.multistate(etime, event=cause, states=("pcm", "death"))
fg = gw.FineGray("pcm").fit(cr, mg[["age", "sex"]])
fg
```

Passing `exponentiate=True` to `tidy` reports the subdistribution hazard ratios (with
their confidence limits) instead of the log-scale coefficients; pass `format=` to choose
the backend (here, Polars):

```{python}
gw.tidy(fg, exponentiate=True, format="polars")
```

MultiState()

Aalen-Johansen estimator of multi-state transition and occupancy probabilities.

Given counting-process intervals `(start, stop]` each labelled with the state occupied
(`state`) and the state transitioned to at `stop` (`event`, or a censoring marker), this
forms the Aalen-Johansen product $P(0, t) = \prod (I + dA(s))$ and reports the state
occupancy probabilities over time. Occupancy probabilities are validated to tolerance
against R's `survfit` multi-state `pstate`. (Competing risks and Kaplan-Meier are special
cases handled by `AalenJohansen` and `KaplanMeier`.)

Returns
-------
Fitted estimator
    Call `fit()` to produce a fitted estimator with cached results (`states_`,
    `time_`, `occupancy_`, and internal transition matrices), accessible as tidy
    DataFrames.

Examples
--------
The `mgus2` patients occupy three states in turn: `"mgus"` at entry, then possibly
`"pcm"` (plasma-cell malignancy), then `"death"`. Reshape the wide dataset into
counting-process intervals `(start, stop]`, one interval per state occupied, labelled with
the state entered next. A patient who progresses before dying contributes two intervals;
everyone else contributes one. Fitting reports the occupancy probability of each state
over time.

```{python}
import greenwood as gw

mg = gw.load_dataset("mgus2", backend="pandas")
start, stop, state, event = [], [], [], []
for i in range(len(mg)):
    pt, ft = mg["ptime"][i], mg["futime"][i]
    progressed, died = mg["pstat"][i] == 1, mg["death"][i] == 1
    if progressed and pt < ft:
        start += [0, pt]; stop += [pt, ft]; state += ["mgus", "pcm"]
        event += ["pcm", "death" if died else None]
    else:
        start += [0]; stop += [ft]; state += ["mgus"]
        event += ["death" if died else ("pcm" if progressed else None)]
rows = [(a, b, s, e) for a, b, s, e in zip(start, stop, state, event) if b > a]
start, stop, state, event = map(list, zip(*rows))
ms = gw.MultiState().fit(start, stop, state, event, states=("mgus", "pcm", "death"))
ms.to_frame(format="polars")
```


## Group comparisons

The log-rank test, trend tests for ordered groups, the G-rho (Fleming-Harrington) family, and restricted mean survival time (RMST) comparisons.


logrank_test(surv: 'Surv', group: 'Any', *, rho: 'float' = 0.0, gamma: 'float' = 0.0, strata: 'Any' = None) -> 'TestResult'

Compare survival across groups using the weighted log-rank (G-rho) test.

Tests whether survival curves differ significantly across two or more groups using a
chi-square test based on weighted event counts. The test is flexible: with default weights
(Fleming-Harrington rho=0, gamma=0), it gives equal weight to all event times (standard
log-rank). With rho=1, gamma=0 (Peto-Peto), it emphasizes early events where more subjects
are at risk. Other rho/gamma combinations allow custom emphasis on different phases of
follow-up.

The test compares observed vs. expected event counts under the null hypothesis of equal
survival. A large chi-square statistic indicates the groups differ; p-values are
interpreted as the probability of seeing such a statistic or larger if survival is truly
equal.

**Stratification**: Optionally stratify by a nuisance variable (e.g., site, gender) to
compute the test within each stratum, then combine results. This controls for confounding
while testing group differences.

Parameters
----------
surv
    A `Surv` response object representing censored survival times. Supports right-censored
    data (standard time-to-event) or counting-process format (interval-based data with
    entry/exit times). Constructed with `Surv.right()`, `Surv.counting()`, or
    `Surv.multistate()`.
group
    Group labels, one per observation. Can be a Narwhals series (Polars/Pandas), 1-D
    array, or Python sequence. Labels can be strings, integers, or other hashable types.
    Must have the same length as `surv`.
rho, gamma
    Fleming-Harrington weight exponents applied to the pooled Kaplan-Meier survival
    $S(t-)$ at each event time. The weight is $S(t-)^\rho \, (1-S(t-))^\gamma$.

    - `rho=0, gamma=0` (default): Standard log-rank test. Equal weight across all times.
    - `rho=1, gamma=0`: Peto-Peto (Wilcoxon) test. Emphasizes early events.
    - `rho=0, gamma=1`: Tarone-Ware. Alternative early-event emphasis.
    - Other (rho, gamma): Flexible emphasis. Higher values emphasize the chosen phase.

strata
    Optional stratifying factor, one per observation. Same length as `surv`. When
    provided, the test is computed separately within each stratum, then combined
    (stratified test). Use to control for confounding or variable that affects baseline
    hazard but not group differences. Example: stratify by site to account for
    site-specific differences in survival while testing an overall group effect.

Returns
-------
TestResult
    A result object with attributes:

    - `statistic`: Chi-square test statistic.
    - `df`: Degrees of freedom (number of groups minus one).
    - `p_value`: Upper-tail chi-square p-value. Small values indicate survival curves differ.
    - `method`: Description of the test (e.g., "Log-rank test", "Stratified log-rank test",
      "G-rho test (rho=1, gamma=0)").
    - `observed`: Dictionary mapping group labels to observed weighted event counts.
    - `expected`: Dictionary mapping group labels to expected event counts under null.

Details
-------
The log-rank test uses the hypergeometric variance for the chi-square statistic, matching
R's `survival::survdiff`. The pooled Kaplan-Meier survivor curve from all groups combined
is used to compute the Fleming-Harrington weights, ensuring the test is consistently
weighted regardless of group sample sizes.

Counting-process data (with entry times) are fully supported, allowing stratification and
left-truncation (delayed entry).

Examples
--------
Test whether survival differs between the two sexes in the bundled `lung` dataset:

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
result = gw.logrank_test(y, group=lung["sex"])
result
```

Extract individual components from the result:

```{python}
result.statistic  # Chi-square statistic
```

```{python}
result.p_value  # P-value for significance
```

```{python}
result.observed  # Observed event counts per group
```

Use the Peto-Peto (Wilcoxon) weighted test to emphasize differences in early survival:

```{python}
gw.logrank_test(y, group=lung["sex"], rho=1, gamma=0)
```

Run a stratified test to control for institution (if available in data):

```{python}
# gw.logrank_test(y, group=lung["sex"], strata=lung["institution"])
```

trend_test(surv: 'Surv', group: 'Any', *, scores: 'Array | None' = None, rho: 'float' = 0.0, gamma: 'float' = 0.0, strata: 'Any' = None) -> 'TestResult'

Test for linear trend across ordered groups using the log-rank test family.

A trend test is a one-degree-of-freedom test for whether survival changes linearly
across ordered groups (e.g., disease stages I, II, III, IV or dose levels: low,
medium, high). It's more powerful than the multi-degree-of-freedom log-rank test when
groups are naturally ordered.

The test assigns numeric scores to each group (default: 0, 1, 2, ... for sorted order)
and tests whether a linear relationship exists between group score and survival. It can
be combined with Fleming-Harrington weights and stratification, just like `logrank_test`.

**When to use trend tests**: Use when groups are naturally ordered and you want a
higher-power test for a linear trend, rather than testing all possible differences.
For exploratory analysis without assuming order, use `logrank_test` instead.

Parameters
----------
surv
    A `Surv` response object representing censored survival times. Supports right-censored
    data or counting-process format.
group
    Group labels (typically ordered categories like 0, 1, 2, 3 for dose levels or stages).
    Can be a Narwhals series, 1-D array, or Python sequence. Must have at least 2 groups.
    Must have the same length as `surv`.
scores
    Numeric scores assigned to each group to define the linear trend. If `None` (default),
    groups are sorted lexicographically and assigned scores 0, 1, 2, ..., (k-1) where k
    is the number of groups. Provide custom scores as a dictionary mapping group labels
    to numeric values (e.g., `{1: 0, 2: 1, 3: 2}` for stage labels 1,2,3) to use
    different scoring (e.g., unequal spacing). Scores can be any real numbers (including
    negative).
rho
    Fleming-Harrington weight exponent applied to the pooled survival probability.
    The default of `0` gives a standard trend test with equal weight across all times.
gamma
    Fleming-Harrington weight exponent applied to (1 - pooled survival probability).
    The default of `0`; combined with `rho=1` gives Peto-Peto (Wilcoxon) trend test
    emphasizing early events. See `logrank_test` for more details on Fleming-Harrington
    weighting.
strata
    Optional stratifying factor. When provided, the trend test is computed separately
    within each stratum, then combined (stratified trend test). Use to control for
    confounding while testing a linear trend.

Returns
-------
TestResult
    A result object with attributes:

    - `statistic`: chi-square test statistic (always 1 degree of freedom).
    - `df`: always `1` for trend tests.
    - `p_value`: upper-tail chi-square p-value.
    - `method`: description of the test, e.g., "Linear trend test", "Stratified linear
      trend test", "G-rho trend test (rho=1, gamma=0)".
    - `observed`: dictionary mapping each group to observed event count (unweighted).
    - `expected`: dictionary mapping each group to expected event count under null.

Details
-------
The test uses a linear contrast with assigned scores:

- $U = \sum_i \text{score}[i] \cdot (O[i] - E[i])$
- $V = \sum_i \text{score}[i]^2 \cdot \text{Var}[i]$
- $\chi^2 = U^2 / V \sim \chi^2(1)$

This is equivalent to fitting a Cox model with group encoded as the numeric score and
testing whether the coefficient is zero using a score test.

Examples
--------
Test for linear trend in survival across ordered cell types:

```{python}
import greenwood as gw

vet = gw.load_dataset("veteran", backend="polars")
y = gw.Surv.right(vet["time"], event=vet["status"])

# Default: cell types are sorted and assigned scores 0,1,2,3
result = gw.trend_test(y, group=vet["celltype"])
result
```

Use custom scores to give different weight to cell type severity:

```{python}
# Scores: adeno and large are low-risk (0,1), smallcell and squamous are high-risk (3,4)
scores = {"adeno": 0, "large": 1, "smallcell": 3, "squamous": 4}
gw.trend_test(y, group=vet["celltype"], scores=scores)
```

Use Peto-Peto weighting to emphasize early differences:

```{python}
# Peto-Peto: rho=1 gives more weight to early event times
gw.trend_test(y, group=vet["celltype"], rho=1, gamma=0)
```

Use Tarone-Ware weighting to emphasize late differences (gamma=1):

```{python}
# Tarone-Ware: gamma=1 gives more weight to late event times
gw.trend_test(y, group=vet["celltype"], rho=0, gamma=1)
```

Stratified by treatment to control for treatment effects:

```{python}
gw.trend_test(y, group=vet["celltype"], strata=vet["trt"])
```

pairwise_logrank_test(surv: 'Surv', group: 'Any', *, rho: 'float' = 0.0, gamma: 'float' = 0.0, strata: 'Any' = None, correction: 'str' = 'holm', format: 'str | None' = None) -> 'Any'

Pairwise log-rank tests for all group pairs with multiple-comparison correction.

Runs the log-rank test on every pair of groups, then adjusts p-values to control for
multiple testing. This answers the question: "Which pairs of groups have significantly
different survival?" when you have more than two groups.

After the global log-rank test (via `logrank_test`) indicates groups differ, this
pairwise test reveals which pairs are significantly different and by how much. P-values
are adjusted across all pairs using a chosen correction method to control the false
discovery rate or family-wise error rate.

**Typical workflow**: First run `logrank_test` to test overall group differences. If
significant, use `pairwise_logrank_test` to identify which pairs differ. The adjusted
p-values account for testing multiple pairs from the same data.

Parameters
----------
surv
    A `Surv` response object representing censored survival times. Supports right-censored
    data or counting-process format. Constructed with `Surv.right()`, `Surv.counting()`,
    or `Surv.multistate()`.
group
    Group labels, one per observation. Can be a Narwhals series, 1-D array, or Python
    sequence. Must have at least 3 unique levels (to create multiple pairs). Must have
    the same length as `surv`.
rho, gamma
    Fleming-Harrington weight exponents for the log-rank test (same as `logrank_test`).
    Default `(0, 0)` gives standard log-rank; `(1, 0)` gives Peto-Peto (emphasizes early
    events).
strata
    Optional stratifying factor. When provided, each pairwise test is stratified by this
    factor (computed within each stratum, then combined). Use to control for confounding.
correction
    Multiple-comparison adjustment method applied across all pairwise p-values:

    - `"holm"` (default): Controls family-wise error rate. Conservative; recommended for
      small numbers of pairs (fewer than ~10).
    - `"bh"`: Benjamini-Hochberg false-discovery rate. Less conservative; recommended for
      many pairs. Allows more false positives but focuses on their rate.
    - `"bonferroni"`: Bonferroni correction. Very conservative; adjusted $p = p \times m$,
      where $m$ is the number of pairs.
    - `"none"`: No adjustment. Use only if you're testing a single pre-planned pair
      (though use `logrank_test` directly in that case).
format
    Output format: `None` (default), `"pandas"`, `"polars"`, or `"pyarrow"`. When
    `None`, a backend is auto-detected (Polars, then Pandas, then PyArrow).

Returns
-------
pandas.DataFrame, polars.DataFrame, or pyarrow.Table
    One row per pair of groups with columns:

    - `group1`, `group2`: The pair of group labels being compared.
    - `statistic`: Chi-square test statistic for the pair.
    - `p_value`: Raw (unadjusted) log-rank p-value for the pair.
    - `p_adjusted`: Adjusted p-value after multiple-comparison correction. Use this for
      significance testing (e.g., p_adjusted < 0.05).

Details
-------
The number of pairs tested is $C(k, 2) = k(k-1)/2$, where $k$ is the number of groups. For
$k=3$, that's 3 pairs; for $k=5$, that's 10 pairs. Larger numbers of pairs can reduce power
per comparison (wider adjusted confidence intervals), so keep the number of groups
reasonable when possible.

The adjustment method affects stringency: Holm controls false discovery more strictly
(lower type-I error, higher type-II error), while Benjamini-Hochberg is more permissive
(higher type-I error rate overall, but controls the proportion of false discoveries).

Examples
--------
Test pairwise survival differences among the four cell types in the `veteran` dataset.
A global log-rank test first shows that cell types differ overall, but doesn't say which
pairs differ:

```{python}
import greenwood as gw

vet = gw.load_dataset("veteran", backend="polars")
y = gw.Surv.right(vet["time"], event=vet["status"])
gw.logrank_test(y, group=vet["celltype"])
```

The pairwise test compares all six pairs of cell types and returns a table with the
test statistic, raw p-value, and adjusted p-value for each pair. Pass `format=` to
choose the backend (here, Polars); use `p_adjusted` for significance testing:

```{python}
pairs = gw.pairwise_logrank_test(y, group=vet["celltype"], format="polars")
pairs
```

Filter to significant pairs (adjusted p-value < 0.05). Request `format="pandas"` here so
we can use boolean-mask filtering:

```{python}
pairs = gw.pairwise_logrank_test(y, group=vet["celltype"], format="pandas")
pairs[pairs["p_adjusted"] < 0.05]
```

Use the Peto-Peto (Wilcoxon) weighting to emphasize early survival differences:

```{python}
gw.pairwise_logrank_test(y, group=vet["celltype"], rho=1, format="polars")
```

Use Benjamini-Hochberg adjustment (less conservative) if you're interested in which pairs
show evidence of differences (false-discovery rate control rather than family-wise error):

```{python}
gw.pairwise_logrank_test(y, group=vet["celltype"], correction="bh", format="polars")
```

TestResult(statistic: 'float', df: 'int', p_value: 'float', method: 'str', observed: 'dict[Any, float]' = <factory>, expected: 'dict[Any, float]' = <factory>) -> None

The outcome of a log-rank group comparison test.

This class stores the results of `logrank_test` or `pairwise_logrank_test` in a
structured format. Access test statistics, significance (p-value), and per-group
observed vs. expected event counts.

Attributes
----------
statistic
    The chi-square test statistic. Larger values indicate stronger evidence against the
    null hypothesis of equal survival across groups.
df
    Degrees of freedom for the chi-square distribution (number of groups minus one for
    `logrank_test`, always 1 for pairwise tests).
p_value
    Upper-tail chi-square p-value. The probability of observing a chi-square statistic
    this large or larger under the null hypothesis of equal survival. Small p-values
    (typically p < 0.05) indicate significant differences between groups.
method
    Human-readable description of the test method and its configuration, e.g.,
    "Log-rank test", "Stratified log-rank test", "G-rho test (rho=1, gamma=0)".
observed
    Dictionary mapping each group label to its observed (actual) weighted event count.
    Useful for understanding which groups contribute more events.
expected
    Dictionary mapping each group label to its expected event count under the null
    hypothesis of equal survival. Comparison of observed vs. expected reveals which
    groups have more or fewer events than expected.

Details
-------
For a significant result (p_value < 0.05), examine the `observed` and `expected`
dictionaries to see which groups experienced more or fewer events than expected. Groups
with observed > expected have worse (shorter) survival; groups with observed < expected
have better (longer) survival.

Examples
--------
Run a log-rank test and examine results:

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
result = gw.logrank_test(y, group=lung["sex"])
result
```

Access individual components. The chi-square statistic:

```{python}
result.statistic
```

The p-value for significance:

```{python}
result.p_value
```

Observed event counts per group (actual events in data):

```{python}
result.observed
```

Expected event counts per group (under null hypothesis):

```{python}
result.expected
```

Test description:

```{python}
result.method
```

rmst_test(surv: 'Surv', tau: 'float', group: 'Any', *, estimand: 'str' = 'difference', strata: 'Any | None' = None, conf_level: 'float' = 0.95) -> 'RMSTResult'

Test for equality of RMST across two or more groups.

Compares restricted mean survival time (RMST) up to a fixed time tau across groups
using a z-test or t-test. Provides point estimate, standard error, confidence interval,
and p-value for the null hypothesis of equal RMST.

For two groups, this is equivalent to a z-test on the RMST difference (default) or
log-ratio (if estimand="ratio").

Parameters
----------
surv
    A right-censored `Surv` response (time-to-event data).
tau
    The restriction time, typically a clinically relevant horizon (e.g., 365, 1825).
group
    Group membership for each observation. Can be array-like or categorical variable.
estimand
    Type of estimand: `"difference"` (default, RMST1 - RMST2), `"ratio"` (RMST1 / RMST2),
    or `"percentage_difference"` ((RMST1 - RMST2) / RMST2 * 100).
strata
    (Optional) Stratification variable for stratified RMST comparison. If provided,
    RMST estimates are combined across strata before computing the test statistic.
conf_level
    Confidence level for confidence intervals (the default is `0.95` for 95% CI).

Returns
-------
RMSTResult
    A result object containing estimate, standard error, confidence interval,
    test statistic, and p-value.

Details
-------
For two groups (i=1,2), the RMST difference is:

$$
\Delta = \mathrm{RMST}_1(\tau) - \mathrm{RMST}_2(\tau)
$$

with standard error:

$$
\mathrm{SE}(\Delta) = \sqrt{\mathrm{SE}(\mathrm{RMST}_1)^2 + \mathrm{SE}(\mathrm{RMST}_2)^2}
$$

assuming independence. The z-statistic is $Z = \Delta / \mathrm{SE}(\Delta)$, with
two-tailed p-value from the standard normal.

For ratio estimand, the log-ratio variance uses the delta method:

$$
\mathrm{SE}(\log R) = \sqrt{\frac{\mathrm{SE}_1^2}{\mathrm{RMST}_1^2} +
\frac{\mathrm{SE}_2^2}{\mathrm{RMST}_2^2}}
$$

Examples
--------
Test RMST difference between two treatment groups:

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
result = gw.rmst_test(y, tau=365, group=lung["sex"])
result.estimate   # RMST difference
result.p_value    # significance
```

Using ratio estimand:

```{python}
result_ratio = gw.rmst_test(y, tau=365, group=lung["sex"], estimand="ratio")
```

rmst_diff(surv: 'Surv', tau: 'float', group: 'Any', *, strata: 'Any | None' = None, conf_level: 'float' = 0.95) -> 'Any'

Compute RMST difference between two groups with confidence interval.

Convenience function that calls `rmst_test()` with `estimand="difference"` and returns
a DataFrame with the comparison results.

Parameters
----------
surv
    A right-censored `Surv` response.
tau
    The restriction time.
group
    Group membership.
strata
    (Optional) Stratification variable.
conf_level
    Confidence level for intervals.

Returns
-------
DataFrame or dict
    Comparison results in tabular format.

pairwise_rmst_test(surv: 'Surv', tau: 'float', group: 'Any', *, estimand: 'str' = 'difference', strata: 'Any | None' = None, correction: 'str' = 'holm', conf_level: 'float' = 0.95, format: 'str | None' = None) -> 'Any'

Pairwise RMST tests for all group pairs with multiple-comparison correction.

Compares RMST between all pairs of groups, with optional multiple-comparison adjustment.
This answers the question: "Which pairs of groups have significantly different RMST?"
when you have more than two groups.

Parameters
----------
surv
    A right-censored `Surv` response (time-to-event data).
tau
    The restriction time for RMST calculation.
group
    Group labels, one per observation. Can be array-like or categorical variable.
    Must have at least 2 unique levels to create pairs.
estimand
    Type of estimand: `"difference"` (default), `"ratio"`, or `"percentage_difference"`.
strata
    (Optional) Stratification variable. Each pairwise test is stratified by this factor.
correction
    Multiple-comparison adjustment: `"holm"` (default), `"bh"`, `"bonferroni"`, or `"none"`.
conf_level
    Confidence level for intervals (the default is `0.95`).
format
    Output format: None (auto-detect), `"pandas"`, `"polars"`, or `"pyarrow"`.

Returns
-------
DataFrame
    One row per pair of groups with columns for group1, group2, RMST estimates,
    difference/ratio, confidence interval, test statistic, p-value, and adjusted p-value.

Examples
--------
Compare RMST across multiple groups with pairwise comparisons:

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
# If there are multiple groups, e.g., by stage:
# result = gw.pairwise_rmst_test(y, tau=365, group=lung["stage"])
```

RMSTResult(estimate: 'float', se: 'float', lower_ci: 'float', upper_ci: 'float', statistic: 'float', p_value: 'float', method: 'str', group1: 'Any', group2: 'Any', rmst1: 'float', se1: 'float', rmst2: 'float', se2: 'float', tau: 'float', estimand: 'str' = 'difference', stratified: 'bool' = False, conf_level: 'float' = 0.95) -> None

Results of an RMST comparison test or difference calculation.

This class stores the results of RMST group comparisons in a structured format,
including point estimates, confidence intervals, and hypothesis test statistics.

Attributes
----------
estimate
    The point estimate of RMST difference, ratio, or percentage difference between groups.
lower_ci
    Lower bound of the confidence interval for the estimate.
upper_ci
    Upper bound of the confidence interval for the estimate.
se
    Standard error of the estimate.
statistic
    Test statistic (z-score for Wald test) for the null hypothesis of no difference.
p_value
    Two-tailed p-value for the hypothesis test. Small values (typically < 0.05)
    indicate significant differences between groups.
method
    Human-readable description of the comparison method, e.g.,
    `"RMST difference (tau=365)"`.
group1
    Label of the first group (minuend in difference).
group2
    Label of the second group (subtrahend in difference).
rmst1
    RMST estimate for group 1 at tau.
se1
    Standard error of RMST for group 1.
rmst2
    RMST estimate for group 2 at tau.
se2
    Standard error of RMST for group 2.
tau
    The restriction time tau used in the RMST calculation.
estimand
    The type of estimand: `"difference"`, `"ratio"`, or `"percentage_difference"`.
stratified
    Whether this is a stratified comparison (True) or pooled (False).
conf_level
    Confidence level used for interval estimation (the default is `0.95`).

Examples
--------
Compare RMST between two groups:

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

result = gw.rmst_test(y, tau=365, group=lung["sex"])
result
```

Access individual components:

```{python}
result.estimate  # difference between groups
result.p_value   # significance
result.se        # standard error
```

logrank_n_events(hazard_ratio: 'float', *, power: 'float' = 0.8, alpha: 'float' = 0.05, allocation: 'float' = 0.5, sides: 'int' = 2) -> 'int'

Number of events needed for the log-rank test to reach a target power.

Computes the minimum number of events required for a log-rank test to achieve a specified
power when detecting a given hazard ratio in a two-group survival study. This is the
foundational calculation in trial planning.

Under the proportional hazards assumption (Schoenfeld, 1981), power depends only on the
total number of events, not on follow-up duration, censoring distribution, or sample size
per se. This makes this calculation practical as it tells you "how many events do you need?"
and you plan your study to observe that many events through appropriate enrollment and
follow-up.

Use this when you want to know the required event count; use `logrank_sample_size` to
convert events to total enrollment given an expected event probability.

Parameters
----------
hazard_ratio
    The hazard ratio to detect (group 2 versus group 1). Can be < 1 (protective effect,
    better survival) or > 1 (harmful effect, worse survival). The result is symmetric:
    HR=0.5 and HR=2.0 require the same number of events. Typical ranges:

    - 0.5: 50% hazard reduction (strong effect)
    - 0.67: 33% hazard reduction (moderate effect)
    - 0.8: 20% hazard reduction (small effect)

power
    Target statistical power (default 0.8, i.e., 80%). The probability of detecting the
    effect if it truly exists. Standard choice is 0.8; use 0.9 or 0.95 for higher
    confidence (requires more events).
alpha
    Significance level or Type-I error rate (default 0.05). The probability of rejecting
    the null hypothesis if it's true. Use 0.05 for the conventional two-sided test with
    p < 0.05 threshold.
allocation
    Fraction of subjects in one group (default 0.5, a balanced design). For example:

    - 0.5: Equal allocation (balanced, minimizes total events needed)
    - 0.33 or 0.67: 1:2 allocation (unbalanced)
    - 0.25 or 0.75: 1:3 allocation (more unbalanced)

    Unbalanced allocation increases the events needed; balanced (0.5) is optimal for a
    given total sample size.

sides
    1 (one-sided) or 2 (two-sided, default). A two-sided test checks for differences in
    either direction (one-sided checks only one direction). One-sided tests require fewer
    events for the same power but are appropriate only when direction is known in advance.

Returns
-------
int
    Minimum number of events (rounded up), required across all groups combined.

Details
-------
**Schoenfeld's formula**: The exact number of events is

$$
d = \frac{(z_{1 - \alpha/\mathrm{sides}} + z_{\mathrm{power}})^2}
{p \, (1 - p) \, (\ln \mathrm{HR})^2},
$$

where:

- $z_{1 - \alpha/\mathrm{sides}}$: critical value for significance level and sides
- $z_{\mathrm{power}}$: critical value for desired power
- $p$: allocation fraction (default 0.5)
- $\ln(\mathrm{HR})$: natural log of hazard ratio

Balanced allocation ($p = 0.5$) minimizes the denominator and thus minimizes events needed.

**Practical use**: Once you know the required event count, estimate sample size by dividing
by the expected event rate: `n_subjects = ceil(n_events / prob_event)`. Then design your
study (enrollment, duration, follow-up) to observe that many events.

Examples
--------
For a trial detecting a hazard ratio of 0.5 (50% hazard reduction) with conventional 80%
power and two-sided significance level 0.05, how many events are needed?

```{python}
import greenwood as gw

gw.logrank_n_events(hazard_ratio=0.5)
```

Repeat for higher power (90%) to see the increase:

```{python}
gw.logrank_n_events(hazard_ratio=0.5, power=0.9)
```

Explore sensitivity to effect size. Smaller effects (HR closer to 1) require more events:

```{python}
for hr in [0.5, 0.6, 0.7, 0.8]:
    events = gw.logrank_n_events(hazard_ratio=hr)
    print(f"HR {hr}: {events} events needed")
```

Use `sides=1` for one-sided testing (higher power, fewer events, but assumes direction):

```{python}
gw.logrank_n_events(hazard_ratio=0.5, sides=1)
```

Use `allocation=0.33` for unbalanced assignment (e.g., 1 treated per 2 control):

```{python}
gw.logrank_n_events(hazard_ratio=0.5, allocation=0.33)
```

After determining event count, convert to sample size using expected event probability.
For instance, if 40% of subjects are expected to have events, divide by 0.4:

```{python}
events = gw.logrank_n_events(hazard_ratio=0.5)
subjects = gw.logrank_sample_size(hazard_ratio=0.5, prob_event=0.4)
print(f"Events needed: {events}")
print(f"Subjects needed: {subjects}")
print(f"Ratio: {subjects / events:.1f}")
```

logrank_power(hazard_ratio: 'float', n_events: 'float', *, alpha: 'float' = 0.05, allocation: 'float' = 0.5, sides: 'int' = 2) -> 'float'

Power of the log-rank test given the number of observed events.

Computes the statistical power to detect a specified hazard ratio using the log-rank test,
given a fixed number of observed events in a two-group survival study. This is the inverse
calculation of `logrank_n_events`: instead of finding the events needed for a target power,
this function finds the power achieved with a given number of events.

Power depends on three factors:

1. **Number of events** (`n_events`): More events → higher power
2. **Effect size** (`hazard_ratio`): Larger effects (HR far from 1.0) → higher power
3. **Significance level** (`alpha`): More stringent (smaller alpha) → lower power

Under the proportional hazards assumption, power depends only on the total event count, not
the follow-up duration, censoring rate, or sample size separately. This makes it a practical
tool for updating power calculations as events accumulate during a trial.

Parameters
----------
hazard_ratio
    The hazard ratio to detect (group 2 versus group 1). Can be < 1 (group 2 has lower
    hazard/better survival) or > 1 (group 2 has higher hazard/worse survival). The result
    is symmetric: HR=0.5 and HR=2.0 give the same power.
n_events
    Total number of observed events. Must be positive. Power increases with more events;
    even small trials can have high power if many events occur.
alpha
    Significance level (Type-I error rate, default 0.05). The probability of rejecting the
    null hypothesis when it's true. Use `alpha=0.05` for two-sided tests with p < 0.05
    threshold.
allocation
    Fraction of subjects in one group (default 0.5, a balanced design). For unbalanced
    designs (e.g., 0.3, 0.7), power decreases; balanced allocation minimizes the total
    sample size needed for a target power.
sides
    1 (one-sided) or 2 (two-sided, default). One-sided tests have higher power but test
    directional hypotheses only. Two-sided tests are standard but require more events for
    the same power.

Returns
-------
float
    Statistical power, a value between 0 and 1. Power of 0.8 (80%) is conventional in many
    fields; higher power (0.9, 0.95) requires more events.

Details
-------
**Schoenfeld's formula**: Under proportional hazards, the power of the log-rank test is

$$
\mathrm{Power} = \Phi\!\left(\sqrt{d \, p \, (1-p)} \; |\ln(\mathrm{HR})|
- z_{1-\alpha/\mathrm{sides}}\right),
$$

where $d$ is the number of events, $p$ is the allocation fraction, and $\Phi$ is the
cumulative normal distribution function. This formula is exact under proportional hazards
and asymptotically valid for finite samples.

**Practical use**: During a running trial, as events accumulate, you can use this function
to assess interim power. If interim power is low despite many events, the effect size may
be smaller than anticipated.

Examples
--------
A study expects to observe 60 events over its follow-up period. What power does it have
to detect a hazard ratio of 0.5 (50% hazard reduction)?

```{python}
import greenwood as gw

gw.logrank_power(hazard_ratio=0.5, n_events=60)
```

This power (~0.9) is typical for a well-designed trial. Lower power suggests more events
are needed, or the effect size is smaller than assumed. Compute power for different
effect sizes to understand study sensitivity:

```{python}
for hr in [0.5, 0.6, 0.7, 0.8]:
    power = gw.logrank_power(hazard_ratio=hr, n_events=60)
    print(f"HR {hr}: power = {power:.2%}")
```

Use `sides=1` for a one-sided test (higher power, but assumes direction is known):

```{python}
gw.logrank_power(hazard_ratio=0.5, n_events=60, sides=1)
```

Use unbalanced allocation if one group is larger. Power decreases with imbalance:

```{python}
gw.logrank_power(hazard_ratio=0.5, n_events=60, allocation=0.3)
```

logrank_sample_size(hazard_ratio: 'float', prob_event: 'float', *, power: 'float' = 0.8, alpha: 'float' = 0.05, allocation: 'float' = 0.5, sides: 'int' = 2) -> 'int'

Total sample size needed for the log-rank test to reach a target power.

Computes the number of subjects required to observe enough events for a log-rank test to
achieve a target power when detecting a specified hazard ratio. This function combines two
calculations:

1. First, it computes the required number of events using `logrank_n_events` (Schoenfeld's
   formula under proportional hazards).
2. Then, it converts events to subjects using the expected probability that a subject
   experiences the event during follow-up.

**Workflow**: Start here to plan study size. You provide the expected effect size (hazard
ratio), the fraction of subjects expected to have the event (based on baseline hazard,
accrual, and follow-up time), and your desired power. The result is the total enrollment
needed.

Parameters
----------
hazard_ratio
    The hazard ratio to detect (group 2 versus group 1). Smaller HR (e.g., 0.5 = 50%
    hazard reduction) requires fewer subjects for a given power; larger HR (e.g., 0.8 =
    20% hazard reduction) requires more.
prob_event
    Probability (or fraction) that a subject experiences the event (death, hospitalization,
    etc.) during the study. Range: (0, 1]. Typical values depend on the condition and
    follow-up duration:

    - Rare disease (annual incidence 1%): 0.01-0.05 per year of follow-up
    - Common condition (annual incidence 20%): 0.2 per year
    - Study design: shorter follow-up → lower prob_event; longer follow-up → higher

    This is usually estimated from historical data, Kaplan-Meier curves, or clinical
    judgment. Use sensitivity analysis (try 0.3, 0.4, 0.5) if uncertain.

power
    Target statistical power (default 0.8, i.e., 80%). Interpretation: the probability
    that the study detects the effect if it truly exists. Common choices:

    - 0.80 (80%): Conventional, implies 20% Type-II error rate
    - 0.90 (90%): Higher confidence, requires more subjects
    - 0.95 (95%): Stringent, requires many more subjects

alpha
    Significance level (Type-I error rate, default 0.05). Probability of rejecting the
    null hypothesis if it's true. Use 0.05 for two-sided tests with p < 0.05 threshold;
    use 0.01 for stricter control or 0.10 for more exploratory studies.
allocation
    Fraction of subjects in one group (default 0.5, balanced design). For unbalanced
    allocation (e.g., control-to-treatment ratio 2:1), use `allocation=0.33`. Unbalanced
    allocation increases total sample size needed; use only if required by design or
    logistics.
sides
    1 (one-sided) or 2 (two-sided, default). One-sided tests have higher power and require
    fewer subjects but test directional hypotheses only. Two-sided tests are standard but
    require more enrollment.

Returns
-------
int
    Total number of subjects (enrollment) needed, rounded up. This is the total across all
    groups.

Details
-------
**Relationship between events and subjects**:

    n_subjects = ceil(n_events / prob_event)

More subjects are needed when:

- `prob_event` is low (most subjects censored before event): e.g., 50% power, 50 events
  needed, but if only 25% get the event, you need 200 subjects.
- The effect size is smaller: smaller HR requires more events
- Power is higher: 90% power requires more events than 80%
- Design is unbalanced: 3:1 allocation needs more subjects than 1:1

**Estimating prob_event**: Use Kaplan-Meier curves from historical data or prior studies,
or calculate from baseline rates:
$\text{prob\_event} \approx 1 - \exp(-\lambda_0 \times t_{\text{follow-up}})$.

**Sensitivity analysis**: If prob_event is uncertain, compute sample size for a range of
values (e.g., 0.3 to 0.5) to understand robustness.

Examples
--------
A trial aims to detect a hazard ratio of 0.5 (50% hazard reduction) with 90% power. Based
on historical data, about 40% of subjects are expected to have the event during follow-up.
How many subjects must be enrolled?

```{python}
import greenwood as gw

gw.logrank_sample_size(hazard_ratio=0.5, prob_event=0.4, power=0.9)
```

This sample size (~350) is much larger than the event count from `logrank_n_events`
(~140) because most subjects will be censored before the event occurs.

Perform sensitivity analysis for uncertain event probability. How does sample size change
if only 30% or 50% of subjects have events?

```{python}
for prob in [0.3, 0.4, 0.5]:
    n = gw.logrank_sample_size(hazard_ratio=0.5, prob_event=prob, power=0.9)
    print(f"prob_event={prob}: n={n} subjects, {int(prob * n)} events")
```

Compare sample size for different effect sizes (smaller HR → fewer subjects):

```{python}
for hr in [0.5, 0.6, 0.7]:
    n = gw.logrank_sample_size(hazard_ratio=hr, prob_event=0.4, power=0.9)
    print(f"HR {hr}: n={n} subjects")
```

Use higher power (0.95) if you want to be very confident the effect is detected:

```{python}
gw.logrank_sample_size(hazard_ratio=0.5, prob_event=0.4, power=0.95)
```

Use unbalanced allocation (e.g., 2:1 control:treatment) if logistics require it:

```{python}
gw.logrank_sample_size(hazard_ratio=0.5, prob_event=0.4, power=0.9, allocation=1/3)
```


## Prediction performance

Concordance and the IPCW Brier score.


concordance_index(surv: 'Surv', risk: 'Any') -> 'float'

Harrell's concordance index: discrimination of risk scores against observed survival.

Computes Harrell's C-statistic, a measure of how well a risk score discriminates between
subjects who experience early events and those who survive longer. The concordance index
compares all comparable pairs of subjects: those with an observed event are compared to
those still under observation at the same time or later. Higher risk should correspond to
earlier failure; if predictions match reality better than chance, the index exceeds 0.5.

**Interpretation**:

- 0.5: Random discrimination (no better than a coin flip)
- 0.6-0.7: Moderate discrimination (clinically useful)
- 0.7-0.8: Strong discrimination
- 0.8+: Excellent discrimination (rare in practice)

**Practical use**: Validates a model's ability to rank subjects by risk. After fitting a
Cox model or other survival model, compute the concordance index of its predictions to
assess out-of-sample discrimination. A model with high concordance index generalizes well
to ranking future subjects by risk.

Parameters
----------
surv
    A right-censored `Surv` response (time-to-event data).
risk
    Risk score for each subject, one per observation. Can be a 1-D array, Pandas/Polars
    series, or Python sequence. Higher values indicate higher risk (earlier expected
    failure). Examples: Cox model linear predictor, predicted log-hazard, or predicted
    cumulative incidence.

Returns
-------
float
    Concordance index between 0 and 1. Values above 0.5 indicate the model discriminates
    better than random; below 0.5 indicates worse-than-random discrimination (possibly
    inverted risk scale).

Details
-------
**Pair comparison rule**:

- A subject with an observed event at time t is compared to all subjects still under
  observation at time t or later (including censored subjects at exactly t).
- Pairs are concordant if the subject with the event has higher risk than the subject
  without.
- Pairs with tied risk scores are counted as 0.5 (half-concordant).
- Pairs with the same event time are excluded.

**Censoring handling**: Censored subjects are handled through the comparable pairs
definition. A censored subject at time t can only be compared to subjects with events at
times strictly greater than t. This avoids artificial inflation of concordance from
censored subjects and matches R's `survival::concordance`.

**Relationship to other metrics**: The concordance index is related to the rank correlation
between risk and observed survival. It's invariant to monotonic transformations of risk
(e.g., $\exp(\text{lp})$ vs. $\text{lp}$ both give the same concordance).

Examples
--------
Fit a Cox model on the `lung` dataset and evaluate its discrimination. Higher linear
predictor values should correspond to earlier deaths.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

lp = cox.predict(type="lp")
c_index = gw.concordance_index(y, lp)
c_index
```

A concordance index of ~0.6 indicates moderate discrimination. Compare with baseline
(naive model assuming all subjects are at equal risk):

```{python}
import numpy as np
baseline_c = gw.concordance_index(y, np.zeros(len(y)))
print(f"Baseline: {baseline_c:.3f}")
print(f"Cox model: {c_index:.3f}")
print(f"Improvement: {c_index - baseline_c:.3f}")
```

Evaluate on hold-out test data to assess generalization (use Cox model fit on training
data, predict on test data):

```{python}
# cox = CoxPH().fit(y_train, x_train)
# lp_test = cox.predict(x_test, type="lp")
# c_test = gw.concordance_index(y_test, lp_test)
```

brier_score(surv: 'Surv', survival_prob: 'Any', times: 'Any') -> 'Array'

IPCW (Graf) Brier score of predicted survival probabilities at specified times.

Measures calibration and accuracy of predicted survival probabilities at fixed time points
using inverse-probability-of-censoring-weighted (IPCW) averaging. The Brier score is the
mean squared difference between predicted and observed outcomes, weighted so censored
subjects contribute honestly without bias.

**Interpretation**:

- Ranges from 0 (perfect predictions) to 1 (worst possible).
- Lower is better. A Brier score of 0.25 means, on average, predictions are off by
  $\pm 0.5$ in terms of squared deviation.
- **Null model baseline**: A model predicting 50% survival at every time has Brier score
  $\approx 0.25$. Compare your model to this baseline to assess practical improvement.
- Score typically increases with time (harder to predict farther into the future).

**Practical use**: After fitting a survival model (Cox, parametric, flexible), evaluate
calibration at important clinical horizons (e.g., 1-year, 5-year survival). Compute Brier
scores at multiple times, then use `integrated_brier_score()` for a single summary.

Parameters
----------
surv
    A right-censored `Surv` response (time-to-event data).
survival_prob
    Predicted survival probabilities, shape `(n_subjects, n_times)`. Each entry is a
    predicted probability of surviving beyond the corresponding time. Must be between 0
    and 1. Example: columns from `CoxPH.predict(type="survival", times=...)` (excluding
    the `time` column), transposed so rows are times and columns are subjects.
times
    Evaluation times where Brier scores are computed. 1-D array-like. Must have length
    equal to the second dimension of `survival_prob`.

Returns
-------
ndarray
    Brier score at each time, shape `(len(times),)`. Lower is better.

Details
-------
**Graf (IPCW) Brier score**: The unbiased Brier score under censoring is

$$
BS(t) = E\left[(S(t) - \hat{S}(t))^2 \cdot \text{weights}\right],
$$

where:

- $S(t)$ is the true survival status at time $t$ (1 if alive, 0 if dead)
- $\hat{S}(t)$ is the predicted survival probability
- $\text{weights}$ are inverse-probability-of-censoring: inverse of the censoring survival
  function

Mathematically:

$$
BS(t) = \frac{1}{n} \sum_{\text{dead at } t} \frac{(\hat{S}_i(t))^2}{G(t_i^-)}
        + \frac{1}{n} \sum_{\text{alive at } t} \frac{(1 - \hat{S}_i(t))^2}{G(t)}
$$

where $G(u)$ is the Kaplan-Meier estimate of the censoring distribution (probability of not
being censored).

**Advantages over MSE**: IPCW weighting makes the Brier score unbiased under censoring,
unlike naive MSE which would be biased (censored subjects look artificially "correct").

**Time dependence**: Brier scores typically increase with time in chronic-disease settings
(longer prediction horizons are harder); they may decrease in acute-illness settings
(events cluster early, predictions stabilize).

Examples
--------
Fit a Cox model on the `lung` dataset and evaluate its calibration (how accurately does
it predict survival?) at a few horizons.

```{python}
import greenwood as gw
import numpy as np

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

times = [180, 365, 540]
surv_pred = cox.predict(lung[["age", "sex"]], type="survival", times=times, format="pandas")
probs = surv_pred.iloc[:, 1:].to_numpy().T
brier = gw.brier_score(y, probs, times)
brier
```

Brier scores at three time points. Scores typically increase over time. Compare to a
null model (all subjects at 50% survival) to assess improvement:

```{python}
null_probs = np.full_like(probs, 0.5)
null_brier = gw.brier_score(y, null_probs, times)
print(f"Null model Brier: {null_brier}")
print(f"Cox model Brier: {brier}")
print(f"Improvement: {null_brier - brier}")
```

Summarize Brier scores across times into a single summary via time-averaged Brier score:

```{python}
ibs = gw.integrated_brier_score(y, probs, times)
print(f"Integrated Brier Score: {ibs:.3f}")
```

integrated_brier_score(surv: 'Surv', survival_prob: 'Any', times: 'Any') -> 'float'

Integrated (time-averaged) Brier score across multiple time points.

Summarizes Brier scores computed at multiple time horizons into a single summary metric
via trapezoidal integration. This provides an overall calibration measure that doesn't
emphasize any particular time point.

**Use this to**: Reduce multiple Brier scores (one per time point) to a single number
for model comparison. A single IBS score makes it easier to compare two models or report
a single "calibration quality" metric.

**Interpretation**: Same scale as Brier score (0 = perfect, 1 = worst). Values of
0.15-0.25 are typical for reasonable survival models; values > 0.30 suggest poor
calibration.

Parameters
----------
surv
    A right-censored `Surv` response.
survival_prob
    Predicted survival probabilities, shape `(n_subjects, n_times)`.
times
    Evaluation times (must be at least 2 to define an interval). The IBS is computed as
    the area under the Brier-score curve from times[0] to times[-1], normalized by the
    time span.

Returns
-------
float
    Integrated Brier score (time-averaged). Lower is better.

Details
-------
**Computation**: The integrated Brier score is

$$
IBS = \frac{1}{t_{\max} - t_{\min}} \int BS(t) \, dt
$$

Using trapezoidal rule to approximate the integral. This ensures the summary score
balances contributions from all times without emphasizing early or late horizons.

**Time scale sensitivity**: IBS weights contribution proportionally to time intervals.
If you have many time points clustered early (e.g., 10 times in [0, 100] and 1 time at
[1000]), early times dominate. Use evenly-spaced time points for balanced assessment.

Examples
--------
Fit a Cox model and compute its integrated Brier score over a range of times:

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

times = [180, 365, 540]
surv_pred = cox.predict(lung[["age", "sex"]], type="survival", times=times, format="pandas")
probs = surv_pred.iloc[:, 1:].to_numpy().T
ibs = gw.integrated_brier_score(y, probs, times)
ibs
```

Compare two models via their integrated Brier scores. Lower is better:

```{python}
# cox2 = CoxPH().fit(y, lung[["age", "sex", "ph.ecog"]])  # More covariates
# surv_pred2 = cox2.predict(...)
# ibs2 = gw.integrated_brier_score(y, probs2, times)
# print(f"Model 1 IBS: {ibs:.3f}")
# print(f"Model 2 IBS: {ibs2:.3f}")
# print(f"Better model: {'Model 2' if ibs2 < ibs else 'Model 1'}")
```

Compute integrated Brier score over a wide range of times to get an overall calibration
assessment:

```{python}
times_wide = list(range(100, 700, 50))
surv_pred_wide = cox.predict(
    lung[["age", "sex"]], type="survival", times=times_wide, format="pandas"
)
probs_wide = surv_pred_wide.iloc[:, 1:].to_numpy().T
ibs_wide = gw.integrated_brier_score(y, probs_wide, times_wide)
print(f"IBS over {len(times_wide)} time points: {ibs_wide:.3f}")
```

cross_validate(model: 'Any', surv: 'Surv', covariates: 'Any', *, data: 'Any' = None, k: 'int' = 5, metric: 'str' = 'concordance', times: 'Any' = None, seed: 'int | None' = None) -> 'dict[str, Any]'

Evaluate a survival model's out-of-sample performance using k-fold cross-validation.

Provides an honest, unbiased estimate of model performance by splitting data into folds,
fitting on training folds, and evaluating on held-out test folds. This avoids overfitting
bias that occurs when fitting and scoring on the same data.

**Why cross-validate?** Fitting and scoring on the training data gives overly optimistic
performance estimates. A model may fit the training data well due to overfitting, not
true predictive ability. Cross-validation repeatedly fits on different training splits
and evaluates on held-out data, simulating performance on new subjects.

**Metrics**:

- `"concordance"` (default): Harrell's C-statistic on the test fold. Higher is better
  (0.5 = random, 1.0 = perfect). Requires CoxPH, CoxNet, or AFT model.
- `"brier"`: Integrated IPCW Brier score over specified times. Lower is better
  (0 = perfect calibration, 1 = worst). Requires explicit `times=` parameter.

Parameters
----------
model
    An unfitted estimator instance (e.g., `CoxPH()`, `CoxNet()`, `AFT("weibull")`).
    A fresh copy is fit on each training fold, leaving the passed object unchanged.
    Supported: CoxPH, CoxNet, AFT (for concordance) and any of those (for Brier).
surv
    A `Surv` response (time-to-event data). Can be right-censored or counting-process.
    Weights in the response are carried through the cross-validation.
covariates
    Covariates/predictors for the model. Can be:

    - A 2-D array or pandas/Polars DataFrame with one row per subject
    - A formula string (as in `CoxPH.fit()`), evaluated against `data`

data
    If `covariates` is a formula string, the data frame to evaluate it against.
k
    Number of folds (default 5). Each fold serves as test data once; subjects are split
    randomly and evenly across folds. Typical choices: 5 or 10.
metric
    Performance metric for evaluation:

    - `"concordance"` (default): Harrell's C-statistic. Requires CoxPH, CoxNet, or AFT.
    - `"brier"`: Integrated inverse-probability-of-censoring-weighted (IPCW) Brier
      score. Requires `times=` with at least 2 time points.

times
    For `metric="brier"`, evaluation time points (1-D array-like, length $\ge 2$). The Brier
    score is computed at each time, then integrated (time-averaged). Example:
    `times=[365, 730, 1095]` for 1, 2, 3-year predictions.
seed
    Random seed for fold shuffling, ensures reproducibility. If `None`, results may vary
    between runs. Use a fixed seed for consistent comparisons.

Returns
-------
dict
    Dictionary with keys:

    - `"metric"`: Metric name used (`"concordance"` or `"brier"`).
    - `"k"`: Number of folds.
    - `"scores"`: List of per-fold scores (one per fold).
    - `"mean"`: Mean score across folds (primary summary).
    - `"std"`: Standard deviation of scores (variability estimate).

    For concordance, higher mean is better. For Brier, lower mean is better.

Details
-------
**How folds work**: Subjects are randomly shuffled and split into k roughly equal-sized
groups. On iteration i, fold i is held out for testing, while the other k-1 folds are
combined for training. This repeats k times until each fold has served as test data once.

**Completeness**: Subjects with missing covariates are dropped before folding. This
ensures all folds use the same cleaned data, avoiding alignment issues.

**AFT model note**: For AFT, concordance uses the negated linear predictor (since in AFT,
larger lp means longer survival, opposite to Cox). This is handled automatically.

**Reproducibility**: Set `seed=` to ensure the same folds are used across runs. This is
important for comparing different models or reporting consistent results.

Examples
--------
Evaluate a Cox model with 5-fold cross-validation using concordance:

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
result = gw.cross_validate(
    gw.CoxPH(), y, lung[["age", "sex"]], k=5, metric="concordance", seed=1
)
result
```

Access individual components. The mean concordance across folds:

```{python}
result["mean"]
```

Per-fold scores (variability check):

```{python}
result["scores"]
```

Standard deviation (estimate of generalization uncertainty):

```{python}
result["std"]
```

Use Brier score (calibration) instead of concordance (discrimination):

```{python}
result_brier = gw.cross_validate(
    gw.CoxPH(), y, lung[["age", "sex"]], k=5,
    metric="brier", times=[180, 365, 540], seed=1
)
result_brier
```

Compare two models via cross-validation. Model with higher mean concordance (or lower
mean Brier) generalizes better:

```{python}
# simple_model = gw.CoxPH()
# complex_model = gw.CoxPH()
# simple_cv = gw.cross_validate(simple_model, y, lung[["age"]], seed=1)
# complex_cv = gw.cross_validate(complex_model, y, lung[["age", "sex", "ph.ecog"]], seed=1)
# print(f"Simple model C-index: {simple_cv['mean']:.3f} ± {simple_cv['std']:.3f}")
# print(f"Complex model C-index: {complex_cv['mean']:.3f} ± {complex_cv['std']:.3f}")
```


## Visualization

plotnine survival curves and aligned numbers-at-risk tables.


plot_survival(km: 'KaplanMeier', *, conf_int: 'bool' = True, censor_marks: 'bool' = True, risk_table: 'bool' = False, times: 'Any' = None, xlab: 'str' = 'Time', ylab: 'str' = 'Survival probability', width: 'int' = 500, height: 'int' = 300) -> 'Any'

Plot Kaplan-Meier survival curve(s) with Altair.

Renders one or more Kaplan-Meier survival curves as an interactive Vega-Lite chart. Each
curve shows the proportion of subjects surviving (event-free) over time as a
right-continuous step function, with an optional shaded confidence band and censoring
marks. Stratified fits produce one colored curve per group with a legend.

The result is a composable Altair object: layer, facet, or restyle it with Altair's API,
and it renders interactively (tooltips and zoom) in notebooks and browsers. Pass
`risk_table=True` to stack an aligned numbers-at-risk table beneath the curve (the x-axis
scale is shared). Requires Altair (install with `pip install greenwood[altair]`).

Parameters
----------
km
    A fitted `KaplanMeier` object, unstratified (single curve) or stratified.
conf_int
    If `True` (default), draw the point-wise confidence band as a shaded step area.
censor_marks
    If `True` (default), mark censoring times with `+` symbols on the curve.
risk_table
    If `True`, return an `alt.VConcatChart` stacking the curve over an aligned
    numbers-at-risk table. If `False` (default), return only the curve.
times
    Query times for the numbers-at-risk table (used only if `risk_table=True`). Defaults
    to six evenly spaced, rounded times from 0 to the maximum observed follow-up time.
xlab, ylab
    Axis labels (defaults `"Time"` and `"Survival probability"`).
width, height
    Curve dimensions in pixels (defaults 500 x 300).

Returns
-------
An `alt.LayerChart` (or an `alt.VConcatChart` combining the curve and table if
`risk_table=True`).

Examples
--------
```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y, by=lung["sex"])

gw.viz.altair.plot_survival(km, risk_table=True)
```

risk_table(km: 'KaplanMeier', times: 'Any' = None) -> 'Any'

Return the numbers-at-risk table as a standalone Altair chart.

Examples
--------
```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y, by=lung["sex"])

gw.viz.altair.risk_table(km, times=[0, 250, 500, 750, 1000])
```


## Core kernel

The risk-set / event-table tabulation shared by KM, log-rank, and Cox.


EventTable(time: 'Array', n_risk: 'Array', n_event: 'Array', n_censor: 'Array', strata: 'Array | None' = None) -> None

Per-time risk-set tabulation (optionally within strata).

Every array is aligned row-wise. When `strata` is not `None`, rows are grouped by
stratum (each stratum's times are ascending). Counts are weighted when case weights
are supplied, so they may be floats.

Examples
--------
An `EventTable` is produced by `event_table`. Build one from the bundled `lung`
dataset and view it as a Polars frame with `to_frame`.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
et = gw.event_table(y)
et.to_frame(format="polars")
```

event_table(surv: 'Surv', *, group: 'Any' = None, weights: 'Any' = None) -> 'EventTable'

Tabulate the event history: risk sets and events at each observed time.

Creates a structured summary of the survival data at each unique event time. The table
shows how many subjects were at risk, how many experienced events, and how many were
censored at each time point. This is the foundational data structure used by
non-parametric estimators (Kaplan-Meier, Nelson-Aalen) and the log-rank test.

**Uses**:

- **Verification**: Inspect risk sets to understand data structure and check for
  censoring patterns.
- **Manual calculations**: Compute survival estimates, cumulative event rates, or
  other summaries directly from the risk-set counts.
- **Understanding censoring**: See censoring patterns by time and stratification.
- **Reporting**: Present summary tables in publications (common in clinical trials).

Parameters
----------
surv
    A `Surv` response (time-to-event data). Supports right-censored or counting-process
    format. Weighted responses are supported; weights are incorporated into risk-set
    counts.
group
    Optional grouping variable for stratification, one value per subject. Can be a
    Pandas/Polars series, 1-D array, or Python sequence. When provided, the table is
    split into blocks with a `strata` column, one group per block. Groups appear in
    order of first appearance in the data.
weights
    Optional case weights. Can be a 1-D array or series. If `None` (default), uses
    weights from the `surv` response if present, otherwise treats all subjects as
    weight 1.

Returns
-------
EventTable
    A structured result with the following attributes:

    - `time`: Unique event times (ascending, per stratum if grouped).
    - `n_risk`: Weighted number of subjects at risk (alive/uncensored and under follow-up)
      at each time.
    - `n_event`: Weighted number of events at each time.
    - `n_censor`: Weighted number of censored subjects at each time.
    - `strata` (if grouped): Stratum label for each row.

    Access columns via `.to_frame()` (optionally `format=`), or iterate directly.

Details
-------
**Risk-set definition**: At time $t$, subjects "at risk" are those with:

- Entry time $\le t$ (for counting-process data)
- Exit time $> t$ (not yet having an event or censoring)

For right-censored data, entry is always 0, so the condition simplifies.

**Censoring and events at the same time**: Subjects censored at time t are handled
carefully. A subject censored at exactly time t is at risk for any event at t (following
convention in survival analysis). The table counts them in `n_risk` at time t, but
removes them from `n_risk` at times > t.

**Stratification order**: If grouped, strata appear in the order of first appearance in
the input data, not alphabetically. This allows meaningful orderings (e.g., control,
then treatment).

**Weights**: If weights are provided, all counts (n_risk, n_event, n_censor) are sums
of weights, not subject counts. This handles case weights or frequency weights.

Examples
--------
View the basic event table for the `lung` dataset: how many subjects are at risk,
events, and censoring at each time.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
et = gw.event_table(y)
et.to_frame(format="polars").head(10)
```

The first row shows the first event time: how many subjects were at risk, how many
experienced an event, and how many were censored. Note that `n_risk` decreases over
time as subjects leave the risk set (events or censoring).

Stratify by sex to see event patterns for each group separately:

```{python}
et_sex = gw.event_table(y, group=lung["sex"])
et_sex.to_frame(format="polars").head(15)
```

Now each unique time appears twice (once per stratum) with a `strata` column indicating
the group. This is useful for inspecting whether event rates and censoring patterns
differ by group.

Compute a manual survival estimate from risk-set counts. The survival probability at
time $t$ is the product of `(1 - n_event / n_risk)` over all times $\le t$:

```{python}
import numpy as np
et = gw.event_table(y)
df = et.to_frame(format="pandas")
# Kaplan-Meier survival at each time
df["surv"] = np.cumprod(1 - df["n_event"] / df["n_risk"])
df[["time", "n_risk", "n_event", "surv"]].head(10)
```

This manual calculation matches the Kaplan-Meier estimate from `KaplanMeier().fit()`.


----------------------------------------------------------------------
This is the User Guide documentation for the package.
----------------------------------------------------------------------


## Getting started

### Get started

Greenwood turns time-to-event data into estimates, tests, models, and figures. It is built
on [Narwhals](https://narwhals-dev.github.io/narwhals/) for backend-agnostic compute,
validated against R's `survival` package, and visualized with
[plotnine](https://plotnine.org/).

## The pipeline

Every analysis starts from a `Surv` response and flows through a shared kernel:

1. **`Surv` response**: captures the outcome and its censoring (right, left, interval, or
the counting-process form with left truncation), plus optional weights and multi-state
endpoints. It is validated eagerly, so downstream code can trust it.
2. **Risk-set / event-table kernel**: at each unique event time, how many are at risk, how
many had events, how many were censored. This single tabulation underlies Kaplan-Meier,
the log-rank test, and Cox.
3. **Estimators and models**: Kaplan-Meier and Nelson-Aalen, group tests, Cox regression,
parametric models, competing risks, and multi-state models.
4. **Outputs**: tidy frames, plotnine figures, and publication tables through Great Tables
and Great Summaries.

Every statistic in Greenwood is validated to tolerance against R's `survival` package, so
results match the reference implementation the field already trusts.

## How this guide is organized

The guide builds up from data to models. If you are new to survival analysis, read it in
order; if you know what you need, jump straight to the relevant page.

- **Foundations**: [Survival data and the Surv object](03-survival-data.qmd) explains
  censoring and how to represent it, and [Data sources and formats](04-data-sources.qmd)
  covers loading data from any backend.
- **Description and comparison**: [Kaplan-Meier](05-kaplan-meier.qmd),
  [Comparing groups](06-comparing-groups.qmd), and
  [Visualizing survival](07-visualization.qmd).
- **Regression**: [Cox regression](08-cox-regression.qmd),
  [Cox model diagnostics](09-cox-diagnostics.qmd), and
  [Parametric models](10-parametric-models.qmd).
- **Multiple event types**: [Competing risks](11-competing-risks.qmd) and
  [Multi-state models](12-multi-state.qmd).
- **Evaluation**: [Prediction performance](13-prediction-performance.qmd).

## Next steps

- New to the library? Start with [Installation](01-installation.qmd), then the
  [Quick start](02-quick-start.qmd) for a fast tour.
- New to survival analysis? Begin with
  [Survival data and the Surv object](03-survival-data.qmd).


### Installation

Greenwood targets Python 3.10+ and is not yet on PyPI. Once released:

```bash
pip install greenwood
```

## Optional extras

The core depends only on Narwhals, NumPy, and SciPy. Backend- and feature-specific pieces
are optional extras:

| Extra | Adds |
|---|---|
| `pd` | pandas |
| `pl` | Polars |
| `formula` | formulaic (the `Surv(...) ~ x` formula API) |
| `altair` | Altair (interactive, Narwhals-native visualization) |
| `plotnine` | plotnine (grammar-of-graphics visualization) |
| `tables` | Great Tables (publication tables) |
| `fast` | numba (accelerated kernels) |
| `all` | pd, pl, formula, altair, plotnine, tables |

```bash
pip install "greenwood[all]"
```

## From source

```bash
git clone https://github.com/rich-iannone/greenwood.git
cd greenwood
make install   # pip install -e ".[dev]"
make check     # ruff + pyright + pytest
```

The version is derived from git tags by `setuptools_scm`; no version string is committed.


### Quick start

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

This is a short tour that goes from raw data to a survival curve, a plot, and a fitted model.
Each step is only a few lines. The rest of the guide explains every piece in depth, and links
are collected at the end.

## Represent the data

Survival data pairs a follow-up time with an indicator of whether the event was actually
observed or the subject was censored (still event-free when we last saw them). Greenwood keeps
the two together in a `Surv` object. We use the bundled `lung` dataset, the survival times of
patients with advanced lung cancer.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

The summary reports 228 patients, of whom 165 died during the study; the rest were censored.
One detail to note: this dataset codes `status` as 1 for censored and 2 for dead, so we turn
it into a plain event indicator with `status == 2`.

## Estimate the survival curve

The Kaplan-Meier estimator turns those times into a survival curve, the probability of
surviving past each point in time. Printing the fitted estimator reports the median survival.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

km = gw.KaplanMeier().fit(y)
km
```

The median is 310 days: half of the patients are still alive at 310 days. The two columns
after it give a 95 percent confidence interval for that median, 285 to 363 days.

## Plot the curve

`plot_survival` draws the curve as an interactive chart, with a shaded confidence band and a
notch wherever a patient was censored.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

km = gw.KaplanMeier().fit(y)

gw.plot_survival(km)
```

The curve starts at 1.0 and steps down at each death. It does not reach 0 because some
patients were still alive at the end of follow-up, so their survival time is only known to be
at least that long.

## Model a covariate effect

To ask how a characteristic changes the risk of death, fit a Cox proportional hazards model.
Here we include `age` and `sex`. Printing the model shows a coefficient table.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
cox
```

The `exp(coef)` column holds hazard ratios. For `sex` (coded 1 for male, 2 for female) the
hazard ratio is about 0.60, which means women had roughly 40 percent lower risk of death than
men over the study, holding age fixed. A ratio below 1 is lower risk, above 1 is higher risk.

## Next steps

That is the core loop: represent the data, estimate a curve, visualize it, and model an
effect. Each part of the guide takes one of these steps further.

- [Survival data and the Surv object](03-survival-data.qmd) and
  [Data sources and formats](04-data-sources.qmd) cover the input side in full.
- [Estimating survival with Kaplan-Meier](05-kaplan-meier.qmd),
  [Comparing groups](06-comparing-groups.qmd), and
  [Visualizing survival](07-visualization.qmd) cover description and comparison.
- [Cox regression](08-cox-regression.qmd) and
  [Parametric survival models](10-parametric-models.qmd) cover modeling, and
  [Prediction performance](13-prediction-performance.qmd) covers evaluation.



## Foundations

### Survival data and the Surv object

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

Survival analysis studies the time until an event happens: death, relapse, machine failure,
customer churn, or any other well-defined endpoint. What makes this kind of data special,
and what separates survival analysis from ordinary regression, is that we usually cannot
observe the event time for everyone. Some subjects are still event-free when the study ends,
some drop out early, and some enter late. Handling these partially observed times correctly
is the entire point of the field, and it starts with how you represent the data.

This page explains the structure of survival data and introduces the `Surv` object, the
response type that every estimator and model in Greenwood consumes.

## Why survival data is different

Imagine a study that follows patients for two years and records the time until relapse. At
the end of the study, three situations are possible for any given patient.

The first is a fully observed event: the patient relapsed at a known time, say 8 months. The
second is right censoring: the patient was relapse-free at their last visit, say at 20
months, so all we know is that their true relapse time is greater than 20 months. The third
is a variation on censoring caused by how subjects enter and leave the study.

Right censoring is by far the most common situation, and it is why we cannot simply average
the observed times or run a standard regression on them. A censored time of 20 months is not
the same as an event at 20 months, and treating it as one would badly bias every estimate.
Survival methods are built to use the partial information in a censored observation without
pretending it is a complete one.

::: {#fig-censoring}
<svg viewBox="0 0 680 210" xmlns="http://www.w3.org/2000/svg" style="max-width:680px;width:100%;height:auto;font-family:sans-serif;color:inherit">
<text x="125" y="59" text-anchor="end" font-size="13" fill="currentColor">Subject A</text>
<line x1="140" y1="55" x2="280" y2="55" stroke="currentColor" stroke-width="3"/>
<circle cx="280" cy="55" r="6" fill="#4c78a8"/>
<text x="292" y="59" font-size="13" fill="currentColor">Event (t = 8)</text>
<text x="125" y="109" text-anchor="end" font-size="13" fill="currentColor">Subject B</text>
<line x1="140" y1="105" x2="490" y2="105" stroke="currentColor" stroke-width="3"/>
<circle cx="490" cy="105" r="6" fill="none" stroke="currentColor" stroke-width="2"/>
<text x="502" y="109" font-size="13" fill="currentColor">Censored (t = 20)</text>
<line x1="140" y1="150" x2="560" y2="150" stroke="currentColor" stroke-width="1"/>
<g font-size="11" fill="currentColor" text-anchor="middle">
<line x1="140" y1="150" x2="140" y2="156" stroke="currentColor"/><text x="140" y="170">0</text>
<line x1="210" y1="150" x2="210" y2="156" stroke="currentColor"/><text x="210" y="170">4</text>
<line x1="280" y1="150" x2="280" y2="156" stroke="currentColor"/><text x="280" y="170">8</text>
<line x1="350" y1="150" x2="350" y2="156" stroke="currentColor"/><text x="350" y="170">12</text>
<line x1="420" y1="150" x2="420" y2="156" stroke="currentColor"/><text x="420" y="170">16</text>
<line x1="490" y1="150" x2="490" y2="156" stroke="currentColor"/><text x="490" y="170">20</text>
<line x1="560" y1="150" x2="560" y2="156" stroke="currentColor"/><text x="560" y="170">24</text>
</g>
<text x="350" y="195" text-anchor="middle" font-size="12" fill="currentColor">Months of follow-up</text>
</svg>

Two subjects observed over a study window. Subject A reaches the event at month 8. Subject B is still event-free at month 20, so their true event time is only known to exceed 20: this is right censoring.
:::

Greenwood also supports less common but important patterns: left censoring (the event
happened before a known time), interval censoring (the event happened between two known
times), and left truncation or late entry (a subject only becomes observable after some
delay). All of these are expressed through the same response object.

::: {.callout-note}
## The key idea
A survival observation carries two pieces of information: a time, and a status that says
what the time means (an event, or a censoring point). You must always keep them together.
Greenwood does this for you with the `Surv` object.
:::

## The Surv object

`Surv` is the response you pass to estimators and models. It bundles the times with their
status codes and validates them, so that downstream code can rely on a clean, consistent
representation. You build one with a constructor that names the censoring type.

The most common constructor is `gw.Surv.right`, for right-censored data. It takes the observed
times and an event indicator, where a truthy value (or `1`) means the event occurred and a
falsy value (or `0`) means the observation was censored.

```{python}
import greenwood as gw

# Four subjects: events at 5 and 4, censored at 6 and 9.
y = gw.Surv.right(time=[5, 6, 4, 9], event=[1, 0, 1, 0])
y
```

The `repr` gives a quick summary: the censoring type, the number of observations, and how
many of them were events. You can also reach the derived quantities directly.

```{python}
print("observations:", y.n)
print("events:", y.n_events)
print("censored:", y.n_censored)
```

### Other censoring types

When subjects enter the risk set after time zero, or when a covariate changes during
follow-up, you use the counting-process form. Each observation is an interval `(start,
stop]` with an event indicator at the stop time. This is also how left truncation is
expressed.

```{python}
gw.Surv.counting(start=[0, 2, 1], stop=[5, 6, 4], event=[1, 0, 1])
```

Left censoring, where all you know is that the event had already happened by the time of
observation, uses `gw.Surv.left`. It takes the same `(time, event)` arguments as `gw.Surv.right`,
but the time is an upper bound on the unknown event time rather than a lower one.

```{python}
gw.Surv.left(time=[5, 6, 4], event=[1, 0, 1])
```

Interval censoring, where the event is known only to fall between two times, uses
`gw.Surv.interval`. Use `numpy.inf` for the upper bound to mark right censoring.

```{python}
import numpy as np

gw.Surv.interval(lower=[1, 2, 3], upper=[2, np.inf, 5])
```

Competing risks and multi-state endpoints, where more than one kind of event can occur, use
`gw.Surv.multistate`. The event codes are `0` for censoring and `1, 2, ...` for the competing
causes, and you name the states.

```{python}
gw.Surv.multistate(time=[5, 6, 7, 8], event=[1, 2, 0, 1], states=("relapse", "death"))
```

## Building a response from a data frame

In practice your data lives in a data frame. Greenwood is dataframe-agnostic through
Narwhals, so you can pass pandas or Polars columns directly. The bundled `lung` dataset,
from the North Central Cancer Treatment Group, is a good example.

```{python}
lung = gw.load_dataset("lung", backend="polars")
lung
```

There is a subtlety here that trips up almost everyone at least once. In this dataset, and
in several others that come from R, the `status` column is coded `1` for censored and `2`
for dead, rather than the `0`/`1` convention. Greenwood does not guess at this coding, and
it will raise an error if you pass `1`/`2` values directly. You convert the column
explicitly instead.

```{python}
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

::: {.callout-warning}
## Check your event coding
Greenwood requires the event indicator to be boolean or `0`/`1`, where `1` means the event
occurred. If your data uses another convention, such as `1` for censored and `2` for the
event, convert it first with a comparison like `status == 2`. This deliberate strictness
prevents a silent and serious error where censoring and events are swapped.
:::

## Validation and reproducibility

The `Surv` constructor validates its inputs immediately. Times must be finite and
non-negative, the arrays must have matching lengths, and for the counting-process form each
start must be strictly less than its stop. When something is wrong, you get a clear error at
construction time rather than a confusing failure deep inside an estimator.

```{python}
try:
    gw.Surv.right([5, -1, 3], [1, 1, 1])
except ValueError as error:
    print(error)
```

Another common data quality issue is when a subject's entry time equals their exit time
(zero follow-up time). This indicates a data problem: if a subject has no time at risk, they
cannot experience an event. Greenwood catches this early:

```{python}
#| eval: false
try:
    # Subject enters at time 1, exits at time 1 (follow-up = 0)
    gw.Surv.counting(start=[1, 0], stop=[1, 2], event=[1, 0])
except ValueError as error:
    print(error)
```

Output:

```
Each `start` must be strictly less than its `stop`.
```

If you encounter this error, check your data:

- **Most common fix**: Exclude subjects with `start >= stop` (they have zero or negative follow-up time)
- **Data quality check**: Why do some subjects have identical entry and exit times? Is this a data entry error?
- **If intentional**: If the event truly occurred instantaneously, add a small epsilon to the stop time: `stop = start + 1e-10`

A response also serializes to and from JSON without loss, which is useful for saving an
analysis input or moving it between systems. Calling `to_json()` returns a plain string, so you
can write it to a file or send it over the wire. Here is the serialized form of the `lung`
response, truncated to its first stretch so you can see the shape of it.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

text = y.to_json(indent=None)
text[:200]
```

That string carries the times and status codes together, which is exactly what a `Surv`
object holds. Reading it back with `from_json` reconstructs an equivalent response, and we can
confirm nothing was lost by comparing the JSON of the original and the restored object.

```{python}
restored = gw.Surv.from_json(text)
print("round-trips exactly:", restored.to_json() == y.to_json())
```

For working in memory rather than as text, `to_dict` and `from_dict` perform the same
round-trip with a plain Python dictionary, and `to_frame()` (optionally with a `format=` of
`"polars"`, `"pandas"`, or `"pyarrow"`) gives a tidy, one-row-per-subject view that is
convenient for inspection or export.

```{python}
y.to_frame(format="polars")
```

The dictionary form is the same structure `to_json` serializes, so `gw.Surv.from_dict` rebuilds
an equivalent response from it.

```{python}
gw.Surv.from_dict(gw.Surv.right([5, 6, 4], event=[1, 0, 1]).to_dict())
```

## The risk-set table

Under every non-parametric estimate is one tabulation: at each distinct event time, how many
subjects are still at risk, how many have the event, and how many are censored. This risk-set
table is what Kaplan-Meier, the log-rank test, and Cox all build on, and you can compute it
directly with `event_table`. It matches R's `survfit` tabulation.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

gw.event_table(y).to_frame(format="polars")
```

Each row is a distinct time. The `n_risk` column counts the subjects under observation just
before that time, `n_event` the events at it, and `n_censor` the censorings. Pass `group=` to
tabulate within strata, which adds a `strata` column.

```{python}
gw.event_table(y, group=lung["sex"]).to_frame(format="polars")
```

You rarely need this table directly, but it is the shared foundation the estimators are built
on, and it is handy when you want to check counts by hand or drive a custom calculation.

## Next steps

You now know how to represent survival data and how to avoid the most common coding mistake.
From here you can bring in your own data or start estimating.

- [Data sources and formats](04-data-sources.qmd) shows how to load data from pandas, Polars,
  and other backends, and lists the datasets bundled with the package.
- [Estimating survival with Kaplan-Meier](05-kaplan-meier.qmd) shows how to estimate and
  summarize the survival curve.
- [Comparing groups](06-comparing-groups.qmd) covers the log-rank family of tests.
- The [Quick start](02-quick-start.qmd) is a fast tour of everything if you prefer to skim
  first.


### Data sources and formats

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

Real survival data arrives in many shapes. It might be a pandas data frame you built from a
CSV, a Polars data frame from a Parquet file, a table pulled straight from a database, or one
of the datasets bundled with Greenwood for learning and testing. One of the design goals of
this package is that none of this should matter: you bring the data in whatever form you
already have, and Greenwood works with it directly. This page explains how, and lists the
datasets that ship with the package.

## Backend-agnostic by design

Greenwood does not have its own data frame type and does not force you to convert your data
into one. Instead it is built on [Narwhals](https://narwhals-dev.github.io/narwhals/), a
lightweight compatibility layer that lets a library speak to many data frame backends through
one interface. When you hand Greenwood a data frame or a column, it reads the values through
Narwhals and runs its numerical work on the underlying arrays. Your original object is never
mutated and never copied into a foreign format.

In practice this means you can pass data from any of the backends Narwhals supports. Three of
these are covered directly by the Greenwood test suite, which checks that they all produce
identical estimates on the same data.

- **pandas**, the most widely used data frame library, and the default for the bundled
  datasets.
- **Polars**, a fast, multi-threaded, Arrow-based data frame library, well suited to large
  datasets.
- **PyArrow** tables, the in-memory Arrow format that many file readers and databases produce.

Other backends that Narwhals supports, such as Modin, cuDF, and DuckDB, use the same interface
and are expected to work the same way, though the test suite pins parity on the three above.

::: {.callout-note}
## What "backend-agnostic" buys you
You do not have to rewrite your analysis when you switch data frame libraries. Code written
for a pandas workflow runs unchanged when the data arrives as Polars, because Greenwood reads
both through the same interface. This also means results do not depend on the backend: the
same data produces the same estimates regardless of where the frame came from. The test suite
verifies this for pandas, Polars, and PyArrow on identical inputs.
:::

## Passing columns and passing frames

There are two ways to feed data to Greenwood, and both are always available.

The first is to pass individual columns to a `Surv` constructor. This is how you build the
response. A column can be a pandas `Series`, a Polars `Series`, a NumPy array, or an ordinary
Python list. Greenwood treats them all the same.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

The second is to pass a whole data frame of covariates to a model. The model reads the columns
it needs, converts non-numeric columns to indicator variables, and drops rows with missing
values, just as R's modeling functions do. We fit the model in its own cell, then set it aside
and look at its results through the tidy functions below.

```{python}
cox = gw.CoxPH().fit(y, lung[["age", "sex", "ph.ecog"]])
```

Printing a fitted model shows a readable summary, but for programmatic use you read its
results through tidy output frames. The `tidy` function returns one row per model term, and
passing `exponentiate=True` reports hazard ratios rather than raw coefficients.

```{python}
gw.tidy(cox, exponentiate=True, format="polars")
```

## Choosing a backend

The `load_dataset` function is itself dataframe-agnostic. With its default it returns a Polars
data frame when Polars is installed, falls back to pandas when only pandas is available, and
raises a clear error if neither is. When you want a specific type, name it with the `backend`
argument.

```{python}
lung_pd = gw.load_dataset("lung", backend="pandas")
lung_pl = gw.load_dataset("lung", backend="polars")

type(lung_pd).__module__.split(".")[0], type(lung_pl).__module__.split(".")[0]
```

The backend makes no difference to any result. The same model fit on the pandas frame and on
the Polars frame returns the same coefficients, a property the test suite checks on every
release. To show this, we first build a response from the Polars columns. The resulting `Surv`
object looks the same as one built from pandas: it records the same times and events, because
Greenwood read both through Narwhals.

```{python}
y_pl = gw.Surv.right(lung_pl["time"], event=(lung_pl["status"] == 2))
y_pl
```

Fitting a Cox model on the Polars-backed response and covariates, then viewing the result with
`to_frame()`, gives the same coefficient table you would get from the Polars frame.

```{python}
gw.CoxPH().fit(y_pl, lung_pl[["age", "sex"]]).to_frame(format="polars")
```

::: {.callout-tip}
## Reading from files and databases
Greenwood does not read files itself, and it does not need to. Use your data frame library's
reader (`pandas.read_csv`, `polars.read_parquet`, a database query, and so on) to get a frame,
then pass its columns to `Surv` exactly as above. Anything that produces a supported data
frame is a valid source.
:::

## Tidy output frames

Every estimator and model returns its results as tidy data frames, through the `to_frame()`
method and through the `tidy` and `glance` functions. Pass `format=` (`"pandas"`, `"polars"`,
or `"pyarrow"`) to choose the backend; when you leave it off, Greenwood auto-detects one
(Polars if installed, otherwise pandas, otherwise PyArrow). This makes the results easy to
inspect, filter, join, or pass on to a plotting or table library.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

cox = gw.CoxPH().fit(y, lung[["age", "sex", "ph.ecog"]])

gw.glance(cox, format="polars")
```

The tidy and glance functions follow the same conventions as R's broom package, so a fitted
model always yields a per-term table (`tidy`) and a one-row model summary (`glance`). This is
the seam that lets Greenwood feed publication tables in Great Summaries and Great Tables.

## Datasets bundled with Greenwood

For learning, testing, and reproducible examples, Greenwood ships several classic survival
datasets drawn from R's `survival` package. Load any of them by name with
`gw.load_dataset`, and list them with `gw.available_datasets()`.

```{python}
gw.available_datasets()
```

Each one targets a different situation. The `lung` dataset (228 subjects), from the North
Central Cancer Treatment Group, is the standard right-censored example used throughout this
guide. `veteran` (137 subjects) is a randomized lung cancer trial that includes a cell-type
covariate, which makes it a good fit for group comparisons. `ovarian` (26 subjects) is a small
trial that is convenient for quick checks and teaching. `pbc` (418 subjects) comes from the
Mayo Clinic trial in primary biliary cholangitis and carries many covariates for regression.
`colon` (1858 subjects) records adjuvant therapy for colon cancer and has both recurrence and
death endpoints. Finally, `mgus2` (1384 subjects) follows patients with monoclonal gammopathy
and provides competing risks of progression and death, which the competing-risks and
multi-state chapters build on.

::: {.callout-warning}
## Event coding in the bundled data
Several of these datasets, `lung` among them, code the status column as `1` for censored and
`2` for the event, following R's convention rather than the `0`/`1` convention. Convert the
column explicitly when you build the response, as in the examples above with
`event=(lung["status"] == 2)`. See [Survival data](03-survival-data.qmd) for why Greenwood
does not guess at this coding.
:::

## Next steps

You can now bring survival data into Greenwood from whatever source you have, and you know
what the bundled datasets contain.

- [Estimating survival with Kaplan-Meier](05-kaplan-meier.qmd) puts this data to work with
  the first and most important estimator.
- [Survival data and the Surv object](03-survival-data.qmd) covers the response type in more
  detail if you skipped ahead.
- The [Quick start](02-quick-start.qmd) is a fast end-to-end tour if you prefer to skim first.



## Estimation

### Estimating survival with Kaplan-Meier

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

The Kaplan-Meier estimator is the workhorse of survival analysis. It answers the most basic
question you can ask of time-to-event data: what fraction of subjects are still event-free
at each point in time? The result is the survival curve, a step function that starts at 1.0
and drops at each observed event. This page shows how to estimate it, quantify its
uncertainty, summarize it with medians and restricted means, and compare subgroups.

Every example here uses the bundled `lung` dataset and the `Surv` response introduced in
[Survival data and the Surv object](03-survival-data.qmd).

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

The response `y` is a `Surv` object holding one entry per subject: the follow-up time and
whether that time ended in an event or a censoring. Displaying it shows a compact summary
where a trailing `+` marks a censored observation. Everything on this page is estimated from
this single object.

## Estimating the survival function

You fit the estimator by calling `fit` with a response. The `KaplanMeier` object follows the
familiar fit-then-inspect pattern, so we fit it first and set it aside.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

km = gw.KaplanMeier().fit(y)
km
```

The fitted `km` object holds the whole curve internally. The most direct way to read it is
`to_frame()`, which lays the estimate out one row per observed time. We show the first
few rows.

```{python}
km.to_frame(format="polars")
```

Each row corresponds to a distinct observed time. The columns tell a complete story:
`n_risk` is the number of subjects still under observation just before that time,
`n_event` and `n_censor` count the events and censorings that occurred at it, `estimate` is
the Kaplan-Meier survival probability, and `conf_low` and `conf_high` bound it. The estimate
only steps down at times where events occur; censorings reduce the risk set but do not move
the curve.

The survival probability at any set of times is available through `predict`, which
evaluates the step function. It returns a plain array with one probability per requested
time, in the same order you asked for them.

```{python}
km.predict([180, 365, 730])
```

::: {.callout-note}
## Reading the survival probability
A survival probability of 0.53 at 365 days means the estimator predicts that 53 percent of
subjects remain event-free one year in. Because the curve is a step function, the value
holds constant between event times.
:::

## Confidence intervals

The confidence band around the curve comes from Greenwood's variance formula. Greenwood
supports three transforms for the interval, chosen with `conf_type`. The default is `"log"`,
which matches the default in R's survival package. The `"log-log"` transform keeps the
limits inside the valid range from 0 to 1 more reliably at the tails, and `"plain"` gives a
symmetric interval on the probability scale.

You set the transform when you construct the estimator, so we fit a second `KaplanMeier` with
`conf_type="log-log"` and keep it separate from the default fit above.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

km_loglog = gw.KaplanMeier(conf_type="log-log").fit(y)
km_loglog
```

Pulling out just the estimate and its interval limits lets you compare this fit against the
default one. We select those columns to keep the table narrow.

```{python}
km_loglog.to_frame(format="polars")[["time", "estimate", "conf_low", "conf_high"]]
```

The point estimates are identical across transforms; only the interval limits differ. Choose
the transform before fitting, and report which one you used.

### Customizing the confidence level

By default, Greenwood computes 95% confidence intervals. You can adjust this with the
`conf_level` parameter, useful for regulatory submissions or sensitivity analyses that require
different coverage levels (e.g., 90% or 99%).

```{python}
km_90 = gw.KaplanMeier(conf_level=0.90).fit(y)
km_99 = gw.KaplanMeier(conf_level=0.99).fit(y)

# Compare intervals at a few time points
km_90.to_frame(format="polars")[["time", "estimate", "conf_low", "conf_high"]].head(3)
```

Note that narrower confidence levels (0.90) produce tighter bands, while wider levels (0.99)
give more conservative coverage. The confidence level affects the median and quantile intervals
as well.

### Comparing confidence interval transforms

The three confidence interval transforms give the same point estimates; they differ only in how
they construct the interval. Here is a comparison to help you choose:

```{python}
import polars as pl

# Fit KM with each transform
km_plain = gw.KaplanMeier(conf_type="plain").fit(y)
km_log = gw.KaplanMeier(conf_type="log").fit(y)
km_loglog = gw.KaplanMeier(conf_type="log-log").fit(y)

# Compare interval widths at selected quantile times
results = []
for km, name in [(km_plain, "plain"), (km_log, "log"), (km_loglog, "log-log")]:
    df = km.to_frame(format="polars")
    # Select every 30th row to show interval differences across follow-up
    df_subset = df[::30].select(["time", "estimate", "conf_low", "conf_high"])
    df_subset = df_subset.with_columns(pl.lit(name).alias("transform"))
    results.append(df_subset)

pl.concat(results)
```

**Transform summary:**

- log** (default): symmetric on the log scale; matches R's survival package. Good all-purpose choice.
- log-log: bounds intervals more reliably at the tails; best when survival is near 0 or 1.
- plain: symmetric on probability scale; intuitive but can exceed [0,1] at extreme tails.

Choose the transform before fitting and report which you used.

## Median survival and quantiles

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)

km.median(ci=True)
```

The result is the median together with its confidence limits, obtained by inverting the
confidence band. Other quantiles work the same way through `quantile`.

```{python}
km.quantile(0.25)
```

::: {.callout-tip}
## When the median is not defined
If the survival curve never falls to 0.5 within the observed follow-up, the median is not
estimable and Greenwood returns `nan`. This is not an error; it is an honest statement that
the data do not reach the median.
:::

## Restricted mean survival time

The mean survival time is usually not estimable from censored data, because the tail of the
curve is unknown. The restricted mean survival time, or RMST, solves this by measuring the
area under the survival curve up to a chosen horizon `tau`. It has a direct interpretation
as the average event-free time over the first `tau` units, and it is a good summary when
survival curves cross, which makes hazard ratios hard to interpret.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)

km.rmst(365, ci=True)
```

This reports the average number of event-free days during the first year, with a confidence
interval. Choose `tau` based on a clinically or practically meaningful horizon, and use the
same `tau` when comparing groups.

## Restricted mean residual life

A related question is how much additional time a subject can expect, given that they have
already survived to some landmark time `s`. The restricted mean residual life answers this: it
is the area under the conditional survival curve from `s` to `tau`, divided by the survival at
`s`. In other words, it is the restricted mean survival measured from `s` onward, among
subjects still event-free at `s`. Setting `s=0` recovers the restricted mean survival time.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier().fit(y)

km.rmrl(180, 730, ci=True)
```

Here the value is the expected number of additional event-free days between 180 and 730 days,
for a subject known to be alive at 180 days. This is useful for updating a prognosis partway
through follow-up, and it complements the conditional survival curves shown in
[Cox model diagnostics](09-cox-diagnostics.qmd).

## Comparing subgroups

Passing a grouping variable to `by` fits a separate curve within each stratum, which is the
starting point for any subgroup comparison. The tidy output gains a `strata` column.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

km_sex = gw.KaplanMeier().fit(y, by=lung["sex"])
km_sex
```

The grouped fit `km_sex` holds one curve per level of `sex`. Its tidy output gains a `strata`
column identifying which curve each row belongs to. We show the first couple of rows within
each stratum so you can see both curves side by side without printing every time point.

```{python}
km_sex.to_frame(format="polars").group_by("strata").head(2)
```

Per-stratum summaries follow naturally. `median` and `rmst` return a dictionary keyed by
stratum when the fit is grouped.

```{python}
km_sex.median()
```

To test whether the difference between groups is statistically significant, rather than just
describing it, use the log-rank family covered in
[Comparing groups](06-comparing-groups.qmd).

### Handling case weights

When your data contains frequency weights (e.g., summarized counts rather than individual
records), or case weights for design-based analysis, pass them via the `weights` parameter.

```{python}
# Example with case weights: create sample weights for demonstration
import numpy as np

weights = np.random.uniform(0.5, 2.0, len(y))

km_weighted = gw.KaplanMeier().fit(y, weights=weights)
km_weighted
```

The estimates adjust for the weights using the weighted Kaplan-Meier formula. Confidence
intervals and other summaries (median, RMST) also account for the weights automatically.

## Nelson-Aalen cumulative hazard

Closely related to survival is the cumulative hazard, the accumulated risk of the event over
time. The Nelson-Aalen estimator provides it directly, and it is often preferred for
diagnostics because it is roughly linear when the hazard is constant.

The estimator follows the same fit-then-inspect pattern as `KaplanMeier`, so we fit it in its
own step.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

na = gw.NelsonAalen().fit(y)
na
```

Reading `na` with `to_frame()` gives the same table layout as before, but now the
`estimate` column holds the cumulative hazard rather than the survival probability. Unlike
survival, it starts at 0 and climbs as events accumulate.

```{python}
na.to_frame(format="polars")
```

## Next steps

You can now estimate survival curves, summarize them, and split them by group.

- [Comparing groups](06-comparing-groups.qmd) tests whether survival differs between
  strata with the log-rank family.
- [Visualizing survival](07-visualization.qmd) turns these curves into publication-quality
  figures with plotnine, including numbers-at-risk tables.
- [Cox regression](08-cox-regression.qmd) moves from describing groups to modeling the
  effect of continuous and multiple covariates.


### Comparing groups

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

Once you have survival curves for two or more groups, the natural question is whether the
differences between them are real or could be explained by chance. The log-rank test and its
relatives answer this. They compare the observed number of events in each group against the
number expected if the groups shared a common survival experience, accumulated across every
event time. This page explains the test, its weighted variants, and how to read the results.

We continue with the `lung` dataset, comparing survival between the two sexes.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

The response `y` is a `Surv` object pairing each subject's follow-up time with an event
indicator, where a trailing `+` marks a censored observation. The grouping variable
`lung["sex"]` splits these same subjects into the two groups we are about to compare.

## The log-rank test

The log-rank test is the standard method for comparing survival across groups. At each event
time it forms a small contingency table of who was at risk and who had an event in each
group, computes the events expected under the null hypothesis of no difference, and combines
the observed-minus-expected discrepancies into a single chi-square statistic.

```{python}
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

result = gw.logrank_test(y, group=lung["sex"])
```

The call returns a `TestResult` object, which we store in `result`. Displaying it prints a
readable summary of the whole test at once.

```{python}
result
```

The result carries the test statistic, its degrees of freedom (one less than the number of
groups), and the p-value. A small p-value is evidence that the survival curves differ. You
can also inspect the weighted observed and expected event counts per group, which show the
direction of the difference. These live on the `result` object as the `observed` and
`expected` attributes, one entry per group.

```{python}
print("observed:", result.observed)
print("expected:", result.expected)
```

Here one group has more observed events than expected and the other fewer, which is the
signature of a survival difference. The test works for any number of groups; with three or
more it produces an overall test on the corresponding degrees of freedom.

::: {.callout-note}
## What the log-rank test does and does not tell you

The test tells you whether the curves differ overall. It does not tell you by how much, at
what times, or in which direction beyond the observed-versus-expected counts. Pair it with a
Kaplan-Meier plot and with effect measures such as a hazard ratio from
[Cox regression](08-cox-regression.qmd) or a restricted-mean-survival difference.
:::

## Weighted tests: the G-rho family

The standard log-rank test weights every event time equally, which makes it most sensitive
to differences in the tails of the curves where events accumulate. Sometimes you care more
about early differences. The Fleming-Harrington, or G-rho, family generalizes the log-rank
test by weighting each event time by a power of the pooled survival probability.

The `rho` parameter controls this. With `rho=0` you get the ordinary log-rank test. With
`rho=1` you get the Peto-Peto test, a Wilcoxon-type test that gives more weight to early
event times and is therefore more sensitive to early separation of the curves.

This returns the same kind of `TestResult` as before, so you can read it the same way. Only
the weighting has changed, so compare its statistic and p-value against the unweighted test
above to see how much the emphasis on early event times matters here.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

gw.logrank_test(y, group=lung["sex"], rho=1)
```

Choosing between them is a matter of what alternative you want to detect, and the choice
should be made before looking at the results rather than by picking whichever gives the
smaller p-value.

::: {.callout-warning}
## Choose the weighting in advance

Running several weighted tests and reporting the most significant one inflates the false
positive rate. Decide on the test, and its `rho`, based on subject-matter reasoning before
you see the data.
:::

::: {.callout-tip}
## Robust handling of edge cases

Greenwood validates input data upfront, so weighted tests handle edge cases robustly.
Times can be zero (events at baseline are valid in survival analysis), and weighted tests
like Fleming-Harrington work correctly without silent data loss or dimension mismatches.
This robustness extends to all censoring patterns and follows the validation rules
documented in the [Validation and reproducibility](03-survival-data.qmd#validation-and-reproducibility)
section.
:::

## Comparing more than two groups

The log-rank test extends to any number of groups, giving one overall test. We switch to the
`veteran` dataset, which records a cell type with four levels.

```{python}
vet = gw.load_dataset("veteran", backend="polars")

y_vet = gw.Surv.right(vet["time"], event=vet["status"])

gw.logrank_test(y_vet, group=vet["celltype"])
```

A small p-value here says the four cell types do not all share the same survival, but it does
not say which ones differ. For that, test each pair and correct for the fact that you are
running several tests at once. `pairwise_logrank_test` does both.

```{python}
gw.pairwise_logrank_test(y_vet, group=vet["celltype"], format="polars")
```

Each row is one pair, with its own statistic, the raw `p_value`, and a `p_adjusted` value
that accounts for the multiple comparisons. The default correction is Holm's; pass
`correction="bh"` for Benjamini-Hochberg, `"bonferroni"`, or `"none"`. Read significance from
the adjusted column, not the raw one.

## Restricted mean survival time (RMST) comparisons

While the log-rank test tells you whether groups differ, it does not tell you by how much.
The restricted mean survival time (RMST) is an intuitive effect measure: the expected
survival time over a fixed time window `tau`, computed as the area under the survival curve.
For example, 1-year RMST is the average survival time during the first year. Unlike hazard
ratios, RMST is on the survival time scale, making it easy to interpret and explain to
non-technical audiences.

### Two-group RMST comparison

Compare RMST between two groups with confidence intervals and hypothesis tests using
`rmst_test()`. Here we compare 1-year RMST between the two sexes in the `lung` dataset.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Compare 1-year RMST between groups
result = gw.rmst_test(y, tau=365, group=lung["sex"])
result
```

The result shows the RMST for each group, their difference, standard error, 95% confidence
interval, and a two-sided p-value testing equality. A positive difference means group 1 has
higher RMST (longer average survival).

You can extract specific components:

```{python}
print(f"RMST difference: {result.estimate:.1f} days")
print(f"95% CI: [{result.lower_ci:.1f}, {result.upper_ci:.1f}]")
print(f"P-value: {result.p_value:.4f}")
```

### Alternative estimands

By default, `rmst_test()` returns the RMST difference. You can also compute the ratio or
percentage difference using the `estimand` parameter.

**Ratio of RMSTs** (group 1 / group 2):

```{python}
result_ratio = gw.rmst_test(y, tau=365, group=lung["sex"], estimand="ratio")
result_ratio.estimate  # Ratio
```

A ratio of `0.8` means group 1 has 80% of the RMST of group 2 (20% shorter average survival).

**Percentage difference** ((group 1 - group 2) / group 2 × 100):

```{python}
result_pct = gw.rmst_test(y, tau=365, group=lung["sex"], estimand="percentage_difference")
result_pct.estimate  # Percentage
```

All three estimands use the delta method to construct standard errors and confidence intervals,
and support flexible confidence levels via the `conf_level=` parameter.

### Pairwise RMST comparisons

When you have three or more groups, use `pairwise_rmst_test()` to compare all pairs with
multiple-comparison adjustment. Here we test which cell types in the `veteran` dataset differ
in 1-year RMST.

```{python}
vet = gw.load_dataset("veteran", backend="polars")

y_vet = gw.Surv.right(vet["time"], event=vet["status"])

# Pairwise RMST differences with Holm adjustment
pairs = gw.pairwise_rmst_test(y_vet, tau=365, group=vet["celltype"])
pairs
```

Each row compares one pair of groups, showing the RMST estimates, difference, and both raw
and adjusted p-values. The adjustment method (default: Holm) controls for multiple testing.
Pass `correction="bh"` for Benjamini-Hochberg or `"bonferroni"` for Bonferroni correction.

### Choosing the restriction time tau

The choice of `tau` is critical: RMST at `tau=365` only includes follow-up through 1 year,
while `tau=1825` extends to 5 years. Choose `tau` based on clinical relevance before looking
at the data. Common choices are clinically meaningful milestones (1 year, 5 years, etc.).

```{python}
# Compare 5-year RMST instead
result_5y = gw.rmst_test(y, tau=1825, group=lung["sex"])
result_5y.estimate
```

::: {.callout-note}
## RMST vs. log-rank: complementary tools

- log-rank test: answers whether curves differ overall (yes/no).
- RMST comparison: answers by how much they differ and how long the difference persists.

Use both: run the log-rank test for significance, then estimate RMST difference for the
magnitude and direction of the effect. They often agree, but RMST is easier to communicate
to stakeholders unfamiliar with hazards.
:::

## Stratified tests

Sometimes you want to compare groups while controlling for a nuisance variable, so that the
comparison is made within each level of that variable and then combined. A stratified
log-rank does this. Here we compare the two treatments while stratifying by cell type, which
removes any cell-type imbalance from the comparison.

```{python}
vet = gw.load_dataset("veteran", backend="polars")

y_vet = gw.Surv.right(vet["time"], event=vet["status"])

gw.logrank_test(y_vet, group=vet["trt"], strata=vet["celltype"])
```

The statistic is pooled across the strata, matching R's `survdiff` with a `strata()` term.
Stratification is the test-level counterpart to a stratified Cox model, where each stratum
has its own baseline hazard.

## Testing for trend in ordered groups

When your groups are naturally ordered—disease stages (I, II, III, IV), dose levels (low,
medium, high), or educational attainment (primary, secondary, tertiary), a general log-rank
test asks the right question but it doesn't use the information about ordering. A trend test is
more powerful because it looks specifically for a linear relationship across the ordered levels.

The trend test compares a linear contrast across groups to see if survival changes systematically
with group ordering. It produces a one-degree-of-freedom chi-square test, which is more powerful
than the multi-degree-of-freedom log-rank test when groups are truly ordered.

### Default scoring

By default, groups are sorted and assigned scores 0, 1, 2, ..., (*k* - 1). Here we test whether
survival differs across cell types in the `veteran` dataset, treating them as an ordered grouping.

```{python}
vet = gw.load_dataset("veteran", backend="polars")

y = gw.Surv.right(vet["time"], event=vet["status"])

# Test linear trend across ordered cell types
gw.trend_test(y, group=vet["celltype"])
```

The result is a `TestResult` with `df=1`, indicating a one-degree-of-freedom test. The
`observed` and `expected` event counts per group reveal which groups contribute most to the
trend, and the p-value tells you if the trend is statistically significant.

### Custom scoring

If the ordering is not uniform—for instance, some groups have very different prognoses—you
can specify custom scores. Scores must be numeric and reflect your belief about the ordering.

```{python}
vet = gw.load_dataset("veteran", backend="polars")

y = gw.Surv.right(vet["time"], event=vet["status"])

# Custom scores: adeno and large get weight 1, smallcell and squamous get weight 3
scores = {"adeno": 1, "large": 1, "smallcell": 3, "squamous": 3}
gw.trend_test(y, group=vet["celltype"], scores=scores)
```

Different scores can change the test statistic, reflecting your scientific hypothesis about
how the groups are ordered.

### Weighted trend tests

Like the log-rank test, the trend test can use Fleming-Harrington weights to emphasize early
or late differences. With `rho=1`, the Peto-Peto variant emphasizes early event times, which
is useful if you suspect early separation of the curves.

```{python}
vet = gw.load_dataset("veteran", backend="polars")

y = gw.Surv.right(vet["time"], event=vet["status"])

# Peto-Peto trend test: weights early events more heavily
gw.trend_test(y, group=vet["celltype"], rho=1)
```

### Stratified trend tests

You can also stratify a trend test to control for confounders while still testing for linear
trend across ordered groups. The stratified test computes the trend separately in each stratum
and then combines the results.

```{python}
vet = gw.load_dataset("veteran", backend="polars")

y = gw.Surv.right(vet["time"], event=vet["status"])

# Test trend in celltype while controlling for treatment
gw.trend_test(y, group=vet["celltype"], strata=vet["trt"])
```

::: {.callout-note}
## When to use trend tests vs. pairwise tests

- use `trend_test()` when groups are ordered and you care about detecting a linear trend.
  It's more powerful for this specific hypothesis.
- use `pairwise_logrank_test()` when you want to know which specific pairs of unordered
  groups differ. It tests all pairwise differences, not a linear trend.
- use `logrank_test()` for an overall test of differences without assuming any structure
  on the groups.
:::

## Study design: sample size and power

The log-rank test also answers a planning question asked before any data are collected: how
large must a study be? Under proportional hazards, its power depends on the data only through
the number of events, so Greenwood provides Schoenfeld's calculators. To detect a given hazard
ratio at a target power, `logrank_n_events` returns the events required.

```{python}
gw.logrank_n_events(hazard_ratio=0.5, power=0.9, alpha=0.05)
```

Running it the other way, `logrank_power` gives the power a planned study would have for a
given number of events, which is useful for checking whether an existing dataset can support a
comparison.

```{python}
gw.logrank_power(hazard_ratio=0.5, n_events=60)
```

Events are not the same as subjects. To turn the required events into a sample size, divide by
the probability that a subject has the event during the study, which you estimate from the
expected survival and the planned follow-up. The `logrank_sample_size()` function does this.

```{python}
gw.logrank_sample_size(hazard_ratio=0.5, prob_event=0.4, power=0.9)
```

::: {.callout-note}
## What the calculators assume

These formulas assume proportional hazards, a two-group comparison, and (by default) a
balanced, two-sided test at `alpha=0.05`. Use `allocation` for an unbalanced design and
`sides=1` for a one-sided test. A balanced design needs the fewest events.
:::

## From tests to models

A group comparison test is the right tool when you have one categorical grouping and want a
yes-or-no answer about a difference. When you need to adjust for other variables, handle a
continuous predictor, or estimate the size of an effect, you move to regression. The
log-rank test has a close connection to the Cox model: the score test from a Cox model with
a single group indicator is essentially the log-rank test, which is why the two usually
agree.

## Next steps

You can now describe and formally compare survival across groups.

- [Visualizing survival](07-visualization.qmd) shows how to present grouped curves clearly.
- [Cox regression](08-cox-regression.qmd) estimates hazard ratios and adjusts for multiple
  covariates.
- Revisit [Kaplan-Meier](05-kaplan-meier.qmd) for restricted-mean-survival summaries, which
  give an interpretable effect size to accompany a significant test.


### Visualizing survival

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

A survival curve is far easier to understand as a picture than as a table of numbers. By
default Greenwood draws curves with [Altair](https://altair-viz.github.io/), which produces
interactive charts (tooltips and zoom) and is Narwhals-native: the numbers flow straight from
your Polars fit into the chart with no Pandas or Matplotlib in the way, so a pure
Polars/Narwhals workflow stays pure. This page covers the survival curve plot, confidence
bands, censoring marks, and the numbers-at-risk table. A [plotnine](https://plotnine.org/)
backend, for a grammar-of-graphics workflow, is described near the end.

Altair is an optional dependency. If you installed Greenwood with the `altair` extra, or with
`all`, you already have it. We use the `lung` dataset throughout. First we build the response.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

We then fit a Kaplan-Meier estimator stratified by sex, one curve per group, using the
log-log confidence interval transform. Printing the fitted estimator summarizes each stratum.

```{python}
km = gw.KaplanMeier(conf_type="log-log").fit(y, by=lung["sex"])
km
```

Everything below draws from this fitted `km` object. Before plotting, it also helps to see
the numbers the curve is built from. Calling `to_frame()` on the fitted estimator returns
the step-function estimates in tidy form, one row per stratum and time, with the survival
probability and its confidence limits. The preview shows the first and last rows along with
the full dimensions of the table.

```{python}
km.to_frame(format="polars")
```

Each row is a step in one of the curves. The `time` and `estimate` columns give the point at
which the survival probability drops and its value after the drop, while `conf_low` and
`conf_high` give the band that the plot will shade. The `strata` column identifies the group,
so there is one set of steps per sex, and `n_risk`, `n_event`, and `n_censor` record how many
subjects were at risk, failed, or were censored at each time. The plot turns exactly these
numbers into a picture.

## The survival curve plot

The main entry point is `plot_survival`. Given a fitted `KaplanMeier`, it draws the step
curve for each stratum, a shaded confidence band, and a small notch at each censoring time.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
km = gw.KaplanMeier(conf_type="log-log").fit(y, by=lung["sex"])

gw.plot_survival(km)
```

Each element carries meaning. The steps show the estimated survival probability, the band
shows its uncertainty, and the censoring notches — angled so they sit on the curve without
obscuring it — show where subjects left the risk set without an event. Because the chart is
interactive, you can hover a curve to read the survival probability at a time, and drag or
scroll to zoom. Curves that separate and stay apart suggest a real group difference, which you
would confirm with a log-rank test from [Comparing groups](06-comparing-groups.qmd).

You can turn individual elements off. Removing the band and the censoring marks produces a
cleaner curve for a slide or a small figure.

```{python}
gw.plot_survival(km, conf_int=False, censor_marks=False)
```

## Adding a numbers-at-risk table

A survival plot is much easier to trust when it is accompanied by the number of subjects
still at risk at regular intervals. As the risk set shrinks, the curve becomes less certain,
and the table makes that visible. Pass `risk_table=True` to stack an aligned table beneath
the curve.

```{python}
gw.plot_survival(km, risk_table=True)
```

The table shares the horizontal axis with the curve above it, so the counts line up with the
times they describe — and because they share one x scale, panning or zooming the curve moves
the table with it. You control the times shown with the `times` argument.

```{python}
gw.plot_survival(km, risk_table=True, times=[0, 250, 500, 750, 1000])
```

If you only need the numbers, `risk_table_data` returns them as a tidy frame, which you can
render however you like, including through Great Tables. It is backend-neutral, so it works
whether or not you have a plotting library installed.

```{python}
gw.viz.risk_table_data(km, times=[0, 250, 500, 750, 1000], format="polars")
```

To draw just the aligned table as its own chart, without a curve above it, use `risk_table`.

```{python}
gw.risk_table(km, times=[0, 250, 500, 750, 1000])
```

## Customizing the plot

`plot_survival` returns an ordinary Altair chart, so you refine it with Altair's API. Chaining
`.properties` sets the title and dimensions, and you can keep layering or restyling from there.

```{python}
gw.plot_survival(km).properties(
    title="Survival by sex", width=640, height=360
)
```

::: {.callout-note}
## The plot is a chart object, not an image

`plot_survival` returns an Altair chart (an `alt.LayerChart`, or an `alt.VConcatChart` when
`risk_table=True`). It renders interactively in notebooks and browsers, and you keep composing
it with Altair's API. Nothing about the figure is locked down.
:::

::: {.callout-tip}
## Saving figures

Call `save` on the returned chart. Use `chart.save("survival.html")` for the interactive
version, or `chart.save("survival.png")` (also `.svg`, `.pdf`) for a static image — static
export uses `vl-convert`, which ships with the `altair` extra.
:::

## The plotnine backend

If you prefer a grammar-of-graphics workflow, or you already build figures with ggplot2-style
layers, Greenwood also draws the same curves with [plotnine](https://plotnine.org/) under
`greenwood.viz.plotnine`. It is an optional dependency, installed with the `plotnine` extra
(or with `all`). The plotnine backend is built on Matplotlib and Pandas, so it is not part of
a Pandas-free workflow, but it produces static, publication-ready figures you compose with the
`+` operator.

The API mirrors the default one: `plot_survival` takes the same fitted `KaplanMeier` and the
same `conf_int`, `censor_marks`, `risk_table`, and `times` arguments, but returns a plotnine
`ggplot` (or a plotnine composition when `risk_table=True`).

```{python}
from plotnine import labs, theme_bw

(
    gw.viz.plotnine.plot_survival(km, risk_table=True)
    + theme_bw()
    + labs(title="Survival by sex", x="Days since enrollment", color="Sex")
)
```

Because the result is a plotnine object, you extend it with any plotnine layer, scale, or
theme, and you save it with its `save` method, for example
`plot.save("survival.png", width=8, height=6, dpi=300)`.

## Next steps

You can now produce and customize publication-quality survival figures.

- [Cox regression](08-cox-regression.qmd) introduces models whose coefficients you will
  later visualize as forest plots.
- [Competing risks](11-competing-risks.qmd) produces cumulative incidence curves, which are
  plotted the same way.
- Revisit [Comparing groups](06-comparing-groups.qmd) to pair a grouped figure with a formal
  test.



## Regression

### Cox regression

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

The Kaplan-Meier estimator and the log-rank test describe and compare whole groups. To
measure the effect of a continuous predictor, or to adjust for several variables at once, you
need a regression model. The Cox proportional hazards model is the most widely used tool for
this. It models the hazard, the instantaneous risk of the event among those still at risk,
and expresses how covariates multiply that risk without requiring you to specify the shape of
the baseline hazard. This page covers fitting the model, reading hazard ratios, and the
tests it reports.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

The response `y` is a `Surv` object that pairs each subject's follow-up time with an event
indicator, marking censored times with a `+`. It records 165 events among 228 subjects, and it
is the outcome the Cox model below regresses on the covariates.

## The proportional hazards idea

The Cox model assumes that each covariate multiplies the hazard by a constant factor that
does not change over time. If being in a treatment group halves the hazard at one month, the
model assumes it halves the hazard at every month. That multiplicative factor is the hazard
ratio, and estimating it is the point of the model. The word "proportional" refers to this
constant-factor assumption, which you should check after fitting; see
[Cox model diagnostics](09-cox-diagnostics.qmd).

The great convenience of the Cox model is that it estimates these hazard ratios without
assuming any particular form for how the baseline risk changes over time. This is why it is
called semiparametric, and why it is the default choice for most analyses.

## Fitting a model

You fit the model with a `Surv` response and a set of covariates. Covariates can be a pandas
or Polars data frame; numeric columns are used directly and non-numeric columns are turned
into indicator variables automatically. Rows with missing covariate values are dropped, as
in a standard complete-case analysis.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex", "ph.ecog"]])
```

The call returns a fitted `CoxPH` estimator, held in `cox`. Printing it gives a compact
summary modeled on R's `coxph` output: the coefficient table, the sample size and event
count, and the overall likelihood ratio test.

```{python}
cox
```

The printed summary is meant for reading. When you want the numbers as data, to filter, join,
or plot them, ask for the coefficient table with `to_frame()`, which returns a tidy frame.

```{python}
cox.to_frame(format="polars")
```

The coefficient table reports, for each covariate, the estimated log hazard ratio
(`estimate`), its standard error, a Wald z-statistic and p-value, and a confidence interval.
The coefficients are on the log scale, which is convenient for the arithmetic of the model
but not for interpretation.

### Handling ties in event times

When multiple subjects experience events at the same time, the Cox model needs a method to
handle these "ties". Greenwood supports two tie-handling methods, controlled by the `ties`
parameter in the constructor:

- **"efron"** (default): Efron's method. Recommended; matches R's survival package default.
- **"breslow"**: Breslow's method. Computationally simpler but less accurate when ties are common.

The choice usually has minimal impact on results unless ties are very common. When in doubt,
stick with the default (Efron).

```{python}
import pandas as pd

# Compare methods
cox_efron = gw.CoxPH(ties="efron").fit(y, lung[["age", "sex"]])
cox_breslow = gw.CoxPH(ties="breslow").fit(y, lung[["age", "sex"]])

# Coefficients are nearly identical unless ties are very common
pd.DataFrame({
    "Efron": cox_efron.to_frame(format="polars")["estimate"],
    "Breslow": cox_breslow.to_frame(format="polars")["estimate"],
})
```

## Hazard ratios and their interpretation

To interpret the model you exponentiate the coefficients, which turns log hazard ratios into
hazard ratios. Greenwood does this for you through the tidy layer.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex", "ph.ecog"]])

gw.tidy(cox, exponentiate=True, format="polars")
```

A hazard ratio above 1 means the covariate increases the hazard, and a value below 1 means
it decreases it. For a continuous covariate such as `age`, the hazard ratio is the
multiplicative change in hazard per one-unit increase; a hazard ratio of 1.02 for age means
roughly a 2 percent higher hazard for each additional year. For the `sex` indicator, coded 1
and 2 in this dataset, a hazard ratio below 1 means the higher-coded group has lower risk.

::: {.callout-note}
## Hazard ratio, not risk ratio

A hazard ratio compares instantaneous rates among those still at risk, not the probability of
the event over the whole study. The two are related but not identical, and the difference
matters when events are common. Report hazard ratios as such, and consider an absolute
measure like a survival difference or a restricted mean difference alongside them.
:::

## Model-fit statistics

The `glance` view gives one-row summary statistics for the whole model, including the
log-likelihood, the AIC, and the likelihood-ratio test of the null hypothesis that all
coefficients are zero.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex", "ph.ecog"]])

gw.glance(cox, format="polars")
```

The model reports three classical global tests, all of which assess whether the covariates
jointly improve the fit: the likelihood-ratio test, the Wald test, and the score test. They
usually agree closely, and the likelihood-ratio test is generally preferred.

```{python}
print("likelihood ratio:", round(cox.lr_stat_, 3))
print("Wald:", round(cox.wald_stat_, 3))
print("score:", round(cox.score_stat_, 3))
```

## Handling tied event times

When two or more events occur at exactly the same recorded time, the partial likelihood must
account for the tie. Greenwood defaults to the Efron approximation, which is accurate and is
also the default in R. The Breslow approximation is available and is faster but slightly less
accurate when ties are common.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

gw.CoxPH(ties="breslow").fit(y, lung[["age", "sex"]]).to_frame(format="polars")[["term", "estimate"]]
```

Unless you have a specific reason to match another tool's Breslow output, the Efron default
is the right choice.

::: {.callout-tip}
## Flexibility in tie-handling methods

Greenwood's support for both Efron and Breslow tie-handling methods gives you flexibility
that is useful in several situations:

- matching textbook examples: if you're learning from a textbook that uses SAS (which
defaults to Breslow), you can use `ties="breslow"` to get results that match exactly.
- comparing across tools: if you're validating results across R, SAS, and Python, you
can switch the tie method in Greenwood to match any tool's output with a single parameter
change, rather than rewriting your analysis in a different language.
- sensitivity analysis: you can easily compare whether your results change substantially
between the two methods, which is a good diagnostic when ties are very common.

Most analyses will stick with the default Efron method (which matches R), but this flexibility
means you're never locked into one approach.
:::

## Time-varying covariates

Some covariates change during follow-up: a treatment that starts partway through, a lab value
that is remeasured, or a status that switches once. The Cox model handles these through the
counting-process form, where each subject contributes one row per interval over which their
covariates are constant. Each row records the interval `(start, stop]`, whether the event
happened at its end, and the covariate values that held during it.

We build a small illustrative dataset where a treatment switches on for some subjects. The
response is a `gw.Surv.counting` object rather than `gw.Surv.right`.

```{python}
intervals = pd.DataFrame(
    {
        "subject": [1, 1, 2, 3, 3, 4, 5, 5],
        "start":   [0, 4, 0, 0, 5, 0, 0, 3],
        "stop":    [4, 10, 7, 5, 14, 9, 3, 12],
        "event":   [0, 1, 1, 0, 0, 1, 0, 1],
        "treated": [0, 1, 0, 0, 1, 0, 0, 1],
    }
)

intervals
```

Subject 1, for example, is untreated over `(0, 4]` and treated over `(4, 10]`, with the event
at day 10. We pass the interval endpoints to `gw.Surv.counting` and fit as usual.

```{python}
y_tv = gw.Surv.counting(start=intervals["start"], stop=intervals["stop"], event=intervals["event"])

gw.CoxPH().fit(y_tv, intervals[["treated"]]).to_frame(format="polars")[["term", "estimate", "p_value"]]
```

The risk set at each event time correctly includes only the intervals that span it, using each
subject's covariate values as of that moment. This is exactly R's start-stop `coxph`, and
Greenwood matches it to tolerance.

::: {.callout-note}
## One row per interval, not per subject

The only change from an ordinary Cox fit is the data layout: a subject with a covariate that
changes `k` times contributes `k + 1` rows, and the event indicator is 1 only on the interval
where the event occurred. Left truncation and delayed entry use the same counting-process
response.
:::

::: {.callout-important}
## Subject timelines must start at `0`

Each subject needs their own timeline starting at 0 (subject-relative time), not calendar time.

**Correct** (subject-relative): Each subject's earliest interval starts at 0

- Subject 1: (0, 4], (4, 10] (personal follow-up time)
- Subject 2: (0, 7] (personal follow-up time)

**Incorrect** (calendar time): Subjects enter the study at different dates

- Subject 1: (2024-01-01, 2024-01-04], (2024-01-04, 2024-01-10]
- Subject 2: (2024-06-15, 2024-06-22] (different calendar dates)

If your data uses calendar time (e.g., from a wide dataset where subjects enroll on different dates),
subtract each subject's entry date from their start/stop times.

**Pandas**:

```{python}
# Example data in calendar time (subjects enroll at different dates)
df = pd.DataFrame({
    'subject': [1, 1, 2, 2, 3, 3],
    'start': [0, 10, 365, 375, 730, 740],      # Different enrollment dates
    'stop': [10, 25, 375, 390, 740, 755],
    'event': [0, 1, 0, 1, 0, 1]
})

# Convert calendar time to subject-relative time (pandas)
df['entry_date'] = df.groupby('subject')['start'].transform('min')
df['start_relative'] = df['start'] - df['entry_date']
df['stop_relative'] = df['stop'] - df['entry_date']

# Drop helper column
df = df.drop(columns=['entry_date'])
df
```

**Polars**:

```{python}
import polars as pl

# Example data in calendar time (subjects enroll at different dates)
df = pl.DataFrame({
    'subject': [1, 1, 2, 2, 3, 3],
    'start': [0, 10, 365, 375, 730, 740],      # Different enrollment dates
    'stop': [10, 25, 375, 390, 740, 755],
    'event': [0, 1, 0, 1, 0, 1]
})

# Convert calendar time to subject-relative time (polars)
df = df.with_columns([
    pl.col('start').min().over('subject').alias('entry_date')
]).with_columns([
    (pl.col('start') - pl.col('entry_date')).alias('start_relative'),
    (pl.col('stop') - pl.col('entry_date')).alias('stop_relative')
]).drop('entry_date')

df
```

Now you can use `start_relative` and `stop_relative` in `Surv.counting()`.

**Note**: Greenwood will warn if it detects counting-process data where subjects don't start at 0,
as this usually indicates an error in data preparation.

This matters because risk sets are formed by subjects at risk at each event time. Using calendar
time instead of subject-relative time creates an imbalance where early times have many subjects
at risk and late times have few, distorting the baseline hazard and model coefficients.
:::

## Robust and clustered standard errors

By default, the Cox model uses the model-based (Fisher information) standard errors. When the
model may be misspecified, or when observations are not independent (e.g., repeated measures,
family data), use the sandwich (robust) variance estimator by setting `robust=True`:

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

cox_robust = gw.CoxPH().fit(y, lung[["age", "sex"]], robust=True)

gw.tidy(cox_robust, format="polars")
```

The standard errors adjust upward (more conservative) if there is overdispersion, or
downward if there is underdispersion. For clustered data (e.g., family members, multiple
measurements per subject), pass the cluster variable to ensure standard errors account for
within-cluster dependence:

```{python}
# Simulated example: suppose subjects cluster by family
lung_with_cluster = lung.with_columns(
    family=((pl.int_range(len(lung)) % 10) + 1)
)

cox_cluster = gw.CoxPH().fit(y, lung_with_cluster[["age", "sex"]],
                              robust=True, cluster=lung_with_cluster["family"])

gw.tidy(cox_cluster, format="polars")
```

Clustered sandwich standard errors are wider than model-based ones, reflecting the loss of
information from within-cluster correlation.

Instead of selecting columns yourself, you can describe the model with a formula, passing the
right-hand side as a string and the data frame as `data`. This is convenient for categorical
variables, interactions, and transformations, and it mirrors the notation used in R. The
formula support uses [formulaic](https://matthewwardrop.github.io/formulaic/), installed with
the `formula` extra.

```{python}
gw.CoxPH().fit(y, "age + sex", data=lung).to_frame(format="polars")[["term", "estimate"]]
```

Categorical columns are expanded into indicator terms automatically, and interactions are
written with `*` (which expands to the main effects plus their product) or `:` (the product
alone). A term like `C(ph.ecog)` forces a column to be treated as categorical.

```{python}
gw.CoxPH().fit(y, "age * sex", data=lung).to_frame(format="polars")[["term", "estimate"]]
```

Rows with a missing value in any formula term are dropped, the same complete-case rule used
when you pass columns directly. The formula interface is also available on `AFT`.

## Penalized regression

When you have many covariates, or they are collinear, an unpenalized fit can overfit or become
unstable. `CoxNet` fits an elastic-net penalized Cox model: it shrinks the coefficients and,
for the lasso, sets some of them to exactly zero, which selects variables. The `penalizer`
argument sets the overall strength, and `l1_ratio` mixes the lasso (`1.0`) and ridge (`0.0`)
penalties.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

cols = ["age", "sex", "ph.ecog", "ph.karno", "wt.loss"]

lasso = gw.CoxNet(penalizer=0.05, l1_ratio=1.0).fit(y, lung[cols])
lasso
```

The printed summary reports how many coefficients survived: the lasso has driven the weakest
ones to exactly zero, leaving a smaller model. A pure ridge penalty (`l1_ratio=0.0`) instead
shrinks every coefficient toward zero without removing any. Because penalized estimates are
biased on purpose, `CoxNet` reports coefficients for prediction and selection but not
p-values, and with `penalizer=0` it reduces to the ordinary Cox fit.

The penalty strength is usually chosen by cross-validation. `cross_validate` accepts a
`CoxNet` directly, so you can compare candidate values by held-out concordance.

```{python}
for lam in [0.02, 0.05, 0.1]:
    cv = gw.cross_validate(gw.CoxNet(penalizer=lam, l1_ratio=1.0), y, lung[cols], k=5, seed=1)
    print(f"penalizer={lam}: CV concordance = {cv['mean']:.4f}")
```

::: {.callout-note}
## Standardization

`CoxNet` standardizes covariates to unit variance before applying the penalty, as glmnet
does, so that the penalty treats variables on comparable scales. Coefficients are returned on
the original scale, so you read and use them exactly as you would from `CoxPH`.
:::

## Next steps

You can now fit a Cox model, read hazard ratios, and assess overall fit.

- [Cox model diagnostics](09-cox-diagnostics.qmd) checks the proportional hazards
assumption, computes residuals, predicts survival curves, and adds robust variance and
stratification.
- [Parametric survival models](10-parametric-models.qmd) offers an alternative modeling
approach when you are willing to specify the shape of the survival distribution.
- [Prediction performance](13-prediction-performance.qmd) evaluates how well a fitted model
discriminates and calibrates.


### Cox model diagnostics

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

Fitting a Cox model is only half the work. Before you trust its hazard ratios, you should
check that its central assumption holds, look at residuals for influential or poorly fit
observations, and understand what the model predicts. This page covers the proportional
hazards test, residuals, baseline hazard and survival prediction, stratification, robust
standard errors, and the concordance index. Together these turn a fitted model into a
defensible one.

We begin with the `lung` outcomes. The response `y` is a `Surv` object pairing each follow-up
time with an event indicator, and it is the target every model on this page is fit against.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

The printed response records 165 events among 228 subjects, with a `+` marking each censored
observation. We fit a Cox model to it using `age` and `sex` as covariates. This is the model
we diagnose throughout the page, so we fit it once here and reuse it below.

```{python}
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
```

## Checking the proportional hazards assumption

The Cox model assumes each covariate's effect on the hazard is constant over time. When this
fails, for example if a treatment helps early but not late, the reported hazard ratio is a
misleading average. The Grambsch-Therneau test checks the assumption by looking for a trend
in the scaled Schoenfeld residuals against time. A small p-value is evidence that the
assumption is violated for that covariate.

Calling `cox_zph` returns a `ZPHResult`, an object that bundles the test statistics together.
Displaying it directly gives a compact summary of the test.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

zph = cox.cox_zph()
zph
```

For the full per-covariate breakdown, call `to_frame()` on that result. This is the table
you will read most often.

```{python}
zph.to_frame(format="polars")
```

The table has one row per covariate plus a `GLOBAL` row that tests all covariates jointly.
Large p-values, as here, are reassuring: they give no evidence against proportional hazards.
The test uses a time transform, which defaults to `"identity"`; `"log"` is also available.

::: {.callout-warning}
## A non-significant test is not proof

Failing to reject the proportional hazards assumption does not prove it holds, especially in
small samples. Complement the test by plotting scaled Schoenfeld residuals or by fitting
time-stratified models when you have reason to suspect a time-varying effect.
:::

## Residuals

Residuals reveal how individual observations relate to the fitted model. Martingale
residuals, one per subject, are the difference between the observed number of events and the
number the model expected; large negative values flag subjects who lived much longer than
predicted. Schoenfeld residuals, one per event, underlie the proportional hazards test and
are useful for spotting time trends in a covariate's effect.

The martingale residuals come back as a NumPy array with one entry per subject. There are 228
of them, so we look at the first five rather than the whole array.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

cox.residuals("martingale")[:5]
```

The Schoenfeld residuals come back as a DataFrame instead, with one row per event and one
column per covariate.

```{python}
cox.residuals("schoenfeld", format="polars")
```

## Baseline hazard and predicted survival

Although the Cox model does not assume a shape for the baseline hazard, it can estimate one
after fitting. The baseline cumulative hazard, and the survival curve derived from it,
describe a reference subject.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

cox.baseline_hazard(format="polars")
```

More useful in practice is predicting the survival curve for specific covariate values.
First we build a small frame of the subjects we want predictions for: a 50-year-old with
`sex` 1 and a 70-year-old with `sex` 2. Showing it confirms the covariate values before we
predict from them.

```{python}
import pandas as pd

newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 2]})
newdata
```

Given that new data and a set of times, `predict` returns the estimated survival probability
for each subject.

```{python}
cox.predict(newdata, type="survival", times=[180, 365, 730], format="polars")
```

Each column is a subject from `newdata`, and each row is a requested time. The `type="lp"`
and `type="risk"` options instead return the linear predictor and the relative risk.

### Confidence bands and conditional survival

Pass `ci=True` to add a pointwise confidence band around each predicted curve. Every subject
then gains `_lower` and `_upper` columns, built from the standard error of the cumulative
hazard and matching R's `survfit` (which combines baseline-hazard and coefficient
uncertainty).

```{python}
cox.predict(newdata, type="survival", times=[180, 365, 730], ci=True, format="polars")
```

You can also predict conditional on a subject having already survived some time, which is
useful for updating a prognosis partway through follow-up. With `conditional_after=180`, the
returned probabilities are $P(T > t \mid T > 180) = S(t) / S(180)$.

```{python}
cox.predict(newdata, type="survival", times=[365, 730], conditional_after=180, format="polars")
```

## Stratification

Sometimes a variable violates the proportional hazards assumption but is not of direct
interest, such as enrolling hospital. Stratification lets each level of that variable have
its own baseline hazard while sharing the coefficients of the covariates you do care about.
You pass the stratifying variable to `strata`. Here we stratify by `sex` while keeping `age`
and `ph.ecog` as covariates.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

stratified = gw.CoxPH().fit(y, lung[["age", "ph.ecog"]], strata=lung["sex"])
```

As with any fitted model, we can print `stratified` for a summary, or read its coefficient
table as data with `to_frame()`.

```{python}
stratified.to_frame(format="polars")
```

Notice that the stratifying variable, `sex`, does not appear as a coefficient. Its effect is
absorbed into the separate baselines, which is exactly the point of stratification.

## Robust and clustered standard errors

When observations are not independent, for example when the same subject contributes several
rows or when subjects are grouped into clusters, the model-based standard errors are too
small. The robust, or sandwich, variance corrects for this. Set `robust=True` for the
Lin-Wei sandwich estimator, and pass `cluster` to group correlated observations.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

robust = gw.CoxPH().fit(y, lung[["age", "sex"]], robust=True)

robust.to_frame(format="polars")[["term", "std_error"]]
```

The coefficients are unchanged; only the standard errors, and therefore the p-values and
intervals, are adjusted. The model-based standard errors remain available as
`naive_std_error_` for comparison.

## The concordance index

The concordance index, or C-statistic, measures how well the model's risk ordering agrees
with the observed order of events. It is the probability that, for a random comparable pair
of subjects, the one who failed first had the higher predicted risk. A value of 0.5 is no
better than chance and 1.0 is perfect discrimination.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

cox.concordance()
```

## Next steps

You can now validate a Cox model and use it to predict.

- [Prediction performance](13-prediction-performance.qmd) extends discrimination and adds
the Brier score for calibration, including for external risk scores.
- [Parametric survival models](10-parametric-models.qmd) provide an alternative when the
proportional hazards assumption is untenable.
- Revisit [Cox regression](08-cox-regression.qmd) for the modeling basics if any of this is
unfamiliar.


### Parametric survival models

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

The Cox model deliberately avoids specifying the shape of the baseline hazard. Parametric
models take the opposite approach: they assume the survival times follow a particular
distribution, such as Weibull or log-normal. In exchange for that assumption you gain a fully
specified model that can extrapolate beyond the observed follow-up, produce smooth survival
and hazard curves, and sometimes fit more efficiently. Greenwood provides these as
accelerated failure time models, which describe how covariates stretch or compress the time
scale. This page shows how to fit them and how to read their coefficients.

We work from the `lung` outcomes. The response `y` is a `Surv` object that pairs each
follow-up time with an event indicator, and it is what every model on this page is fit
against.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

The printed response shows 165 events among 228 subjects, with a `+` marking each censored
time.

## The accelerated failure time model

An accelerated failure time model, or AFT model, works on the logarithm of survival time. It
says that covariates act by multiplying the time scale: a covariate might make time pass
twice as fast, halving survival, or twice as slowly, doubling it. This is a different and
often more intuitive framing than the Cox model's multiplication of hazards. The model is
fit by maximum likelihood, and it includes an intercept because it estimates the actual
location of the survival distribution, not just relative effects.

You fit one by naming a distribution. The Weibull is the default and the most common choice.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

aft = gw.AFT("weibull").fit(y, lung[["age", "sex"]])
aft
```

Printing the fitted model gives a summary in the style of R's `survreg`: the coefficient
table, the scale parameter, the sample size, and the log-likelihood. For the coefficients as
data, pass the model to `gw.tidy`, which returns a tidy frame.

```{python}
gw.tidy(aft, format="polars")
```

The coefficients are on the log-time scale. A positive coefficient lengthens survival time,
and a negative coefficient shortens it. This is the opposite direction from a Cox hazard
ratio, where a positive coefficient means higher risk and shorter survival, so take care when
comparing the two.

## Choosing a distribution

Greenwood supports four parametric distributions, each giving a different hazard shape.

- exponential: Constant hazard (h(t) = λ). Simplest; often too restrictive.
- Weibull: monotonic hazard (increasing or decreasing). Default; most flexible while remaining parsimonious.
- log-normal: hazard rises then falls; many peak early. Good for disease incidence.
- log-logistic: similar to log-normal but with heavier tails.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

for dist in ("exponential", "weibull", "lognormal", "loglogistic"):
    model = gw.AFT(dist).fit(y, lung[["age", "sex"]])
    print(f"{dist:20s} loglik={model.loglik_:10.3f} scale={model.scale_:.4f}")
```

To compare and choose among distributions, use the AIC (available from `glance`), where a
lower value indicates a better fit:

```{python}
import pandas as pd

results = []
for dist in ("exponential", "weibull", "lognormal", "loglogistic"):
    model = gw.AFT(dist).fit(y, lung[["age", "sex"]])
    glance_result = gw.glance(model, format="pandas")
    glance_result["distribution"] = dist
    results.append(glance_result)

comparison = pd.concat(results)[["distribution", "loglik", "aic"]]

print(comparison.sort_values("aic"))
```

The scale parameter describes the spread of the distribution. For the exponential, it is fixed at 1 (constant hazard). For Weibull, values > 1 indicate increasing hazard, < 1 indicate decreasing. Other distributions use scale differently; see the model's docstring for details.

## Model comparison and selection

Each call to `gw.glance` returns a one-row DataFrame of model-level summaries for a
single fit. We can compare all four distributions systematically:

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

comparisons = []
for dist in ("exponential", "weibull", "lognormal", "loglogistic"):
    model = gw.AFT(dist).fit(y, lung[["age", "sex"]])
    result = gw.glance(model, format="pandas")
    result["distribution"] = dist
    comparisons.append(result)

comparison = pd.concat(comparisons)[["distribution", "loglik", "aic"]]

print(comparison.sort_values("aic"))
```

Compare the models by reading the AIC column and preferring the lowest value. A difference of
10 in AIC is substantial; differences of 2-3 are subtle and don't warrant changing models.

**When to use each distribution:**

- exponential: only if you're confident hazard is constant (rare).
- Weibull: safe default that fits many data types well.
- log-normal/log-logistic: When you expect a hazard peak early in follow-up.
- Gompertz: for aging-related or mortality data where hazard accelerates exponentially.

::: {.callout-note}
## Parametric or Cox?

Choose a parametric model when you need to extrapolate, want smooth hazard or quantile
predictions, or have a distributional form suggested by theory. Choose the Cox model when
you want to avoid distributional assumptions and care mainly about relative effects. Both are
valid; they answer slightly different questions.
:::

::: {.callout-tip}
## Interpreting the acceleration factor

Exponentiating an AFT coefficient gives a time-acceleration factor. A factor of 1.5 for a
covariate means subjects with a one-unit-higher value survive 1.5 times as long on average,
holding other covariates fixed.
:::

## Predicting survival and quantiles

Because a parametric model specifies the whole distribution, it can predict smoothly at any
time and extrapolate beyond the observed follow-up. We predict for two new subjects, a
50-year-old and a 70-year-old, both with `sex` coded 1.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
aft = gw.AFT("weibull").fit(y, lung[["age", "sex"]])

newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 1]})

aft.predict(newdata, type="quantile", p=[0.25, 0.5, 0.75], format="polars")
```

With `type="quantile"` the model returns survival-time quantiles: each column is a subject and
each row a failure probability. The middle row (`p = 0.5`) is the predicted median survival
time, and the older subject's is shorter, as expected. These quantiles match R's `survreg`.

For the survival curve itself, use `type="survival"` with the times you want.

```{python}
aft.predict(newdata, type="survival", times=[180, 365, 730], format="polars")
```

Each column is a subject and each row a requested time, giving the estimated probability of
surviving past that time.

## Flexible parametric (Royston-Parmar) models

The AFT distributions impose a fixed hazard shape. When none of them fits well but you still
want the smooth, extrapolatable curves of a parametric model, a Royston-Parmar model is a good
middle ground. It models the log cumulative hazard as a restricted cubic spline in log time,
so the baseline shape is estimated from the data rather than assumed. The flexibility is set
by `df`, the number of spline degrees of freedom.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

rp = gw.RoystonParmar(df=3).fit(y, lung[["age", "sex"]])

rp.to_frame(format="polars")
```

The `gamma` terms are the spline coefficients for the baseline log cumulative hazard, and the
named terms are the covariate effects on that scale.

### Choosing the degrees of freedom

The `df` parameter controls spline flexibility. `df=1` is exactly a Weibull model (no spline);
higher values let the hazard adapt to the data shape. More flexibility fits the observed data
better but risks overfitting and unstable extrapolation.

```{python}
# Compare different df values using loglik
results = []
for df in (1, 2, 3, 4, 5):
    rp_df = gw.RoystonParmar(df=df).fit(y, lung[["age", "sex"]])
    results.append({"df": df, "loglik": rp_df.loglik_})

comparison = pd.DataFrame(results)

print(comparison)
```

Notice that loglik generally improves up to `df=3`, after which gains diminish, suggesting that
additional flexibility is not justified by the data. Two or three degrees of freedom are typical
defaults. Choose based on the trade-off between fit and complexity:

- **df=1**: equivalent to Weibull; simplest.
- **df=2-3**: common choices; balance parsimony and fit.
- **df=4+**: use only if AIC clearly favors it and validation shows good extrapolation.

Prediction works as for the AFT model, with survival, hazard, or cumulative hazard curves.

```{python}
newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 1]})

rp.predict(newdata, type="survival", times=[180, 365, 730], format="polars")
```

## Next steps

You can now fit, compare, and predict from parametric survival models.

- [Cox regression](08-cox-regression.qmd) is the semiparametric alternative, useful when you
prefer not to assume a distribution.
- [Prediction performance](13-prediction-performance.qmd) evaluates any fitted model,
parametric or not.
- [Competing risks](11-competing-risks.qmd) handles the case of more than one type of event.



## Multiple events

### Competing risks

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

So far every method has assumed a single kind of event. Often that is not realistic. A
patient with a precancerous condition might progress to cancer, but they might also die of
something else first, and death from another cause prevents the progression from ever being
observed. These are competing risks: mutually exclusive event types where the occurrence of
one removes the possibility of the others. Analyzing them with ordinary survival methods, one
event type at a time, gives biased and even nonsensical results. This page explains why, and
shows the two standard tools Greenwood provides.

```{mermaid}
%%| label: fig-competing
%%| fig-cap: "A competing-risks model: from an initial state, a subject can move to exactly one of several absorbing event states."
flowchart LR
    A[Event-free] --> B[Progression]
    A --> C[Death]
```

We use the bundled `mgus2` dataset, which follows patients with monoclonal gammopathy. The
competing events are progression to a plasma-cell malignancy and death without progression.
We first build the competing-risks endpoint from two intermediate arrays: a single event time
and a cause code that is 0 for censored, 1 for progression, and 2 for death.

```{python}
import numpy as np
import greenwood as gw

mg = gw.load_dataset("mgus2", backend="polars")

etime = np.where(mg["pstat"] == 1, mg["ptime"], mg["futime"])
cause = np.where(mg["pstat"] == 1, 1, 2 * mg["death"])   # 0 censor, 1 pcm, 2 death
```

The `etime` array picks each subject's event time: the progression time `ptime` when the
subject progressed (`pstat == 1`), and otherwise the follow-up time `futime`. Because a
subject can experience only the first of the competing events, one time per subject is all we
need. Here are the first several values, in months.

```{python}
etime[:8]
```

The `cause` array records which event that time refers to, using the coding noted in the
comment above. Where a subject progressed it is 1 (`pcm`); otherwise it is `2 * death`, which
is 2 when the subject died without progressing and 0 when the subject was censored. Reading
`etime` and `cause` position by position gives each subject's (time, cause) pair.

```{python}
cause[:8]
```

We pass both arrays to `gw.Surv.multistate`, mapping cause code 1 to `pcm` and code 2 to `death`,
with code 0 left as censoring. The resulting response object shows the event times alongside
their decoded states.

```{python}
y = gw.Surv.multistate(etime, event=cause, states=("pcm", "death"))
y
```

## Cumulative incidence, not one-minus-survival

The intuitive but wrong approach is to treat one cause as the event, censor the other, run
Kaplan-Meier, and report one minus the survival curve as the probability of that cause. This
overestimates the probability, sometimes badly, because it treats subjects who died of the
other cause as if they could still progress. The correct quantity is the cumulative incidence
function, which gives the probability of experiencing a specific cause by a given time while
accounting for the competing cause.

The Aalen-Johansen estimator computes it. You fit it to the multi-state response, and it
returns a curve for every cause.

```{python}
mg = gw.load_dataset("mgus2", backend="polars")

etime = np.where(mg["pstat"] == 1, mg["ptime"], mg["futime"])
cause = np.where(mg["pstat"] == 1, 1, 2 * mg["death"])   # 0 censor, 1 pcm, 2 death
y = gw.Surv.multistate(etime, event=cause, states=("pcm", "death"))

aj = gw.AalenJohansen().fit(y)
aj
```

The fitted estimator holds a cumulative incidence curve for each cause. To inspect it we call
`to_frame()` and take the first couple of rows within each cause, so both curves are
visible side by side.

```{python}
aj.to_frame(format="polars").group_by("cause").head(2)
```

Each cause has its own cumulative incidence estimate, standard error, and confidence
interval at every event time. Because the causes are competing, their cumulative incidences
plus the event-free probability sum to one at every time, which is exactly the accounting
that one-minus-Kaplan-Meier gets wrong.

::: {.callout-warning}
## Do not use one-minus-Kaplan-Meier for a competing cause

Censoring the competing event and applying Kaplan-Meier treats those subjects as if they
remained at risk for the cause of interest, which they do not. Always use the cumulative
incidence function for competing-risks probabilities.
:::

You can stratify by a grouping variable, just as with Kaplan-Meier, to compare cumulative
incidence across groups.

```{python}
aj_sex = gw.AalenJohansen().fit(y, by=mg["sex"])
aj_sex
```

The stratified fit carries a separate set of curves for each level of `sex`. Here we pull the
`pcm` cause and keep the last row within each stratum, which is the cumulative incidence of
progression at the final event time in that group.

```{python}
import polars as pl

aj_sex.to_frame(format="polars").filter(pl.col("cause") == "pcm").group_by("strata").tail(1)
```

## Modeling with the Fine-Gray model

The cumulative incidence function describes groups. To model how covariates affect the
incidence of a specific cause, use the Fine-Gray model. It is a regression on the
subdistribution hazard, constructed so that its coefficients describe the effect of
covariates on the cumulative incidence itself, which is usually what you want to report.

```{python}
mg = gw.load_dataset("mgus2", backend="polars")

etime = np.where(mg["pstat"] == 1, mg["ptime"], mg["futime"])
cause = np.where(mg["pstat"] == 1, 1, 2 * mg["death"])   # 0 censor, 1 pcm, 2 death
y = gw.Surv.multistate(etime, event=cause, states=("pcm", "death"))

fg = gw.FineGray("pcm").fit(y, mg[["age", "sex"]])
fg
```

We name `pcm` as the target cause and fit against age and sex. To read the fitted model we
tidy it into a coefficient table, asking for exponentiated estimates so the numbers are on the
hazard-ratio scale rather than the log scale.

```{python}
gw.tidy(fg, exponentiate=True, format="polars")
```

The exponentiated coefficients are subdistribution hazard ratios. A value above 1 means the
covariate increases the cumulative incidence of the target cause, `pcm` here, over time. The
model reports clustered robust standard errors by default, which are appropriate given the
inverse-probability weighting it uses internally.

### Modeling multiple competing causes

In a competing-risks study with more than one cause of interest, you can fit separate
Fine-Gray models, one for each cause, to compare the covariate effects across causes.

```{python}
import pandas as pd

# Fit FineGray models for both competing causes
fg_pcm = gw.FineGray("pcm").fit(y, mg[["age", "sex"]])
fg_death = gw.FineGray("death").fit(y, mg[["age", "sex"]])

# Compare the hazard ratios side-by-side
pcm_tidy = gw.tidy(fg_pcm, exponentiate=True, format="pandas")
death_tidy = gw.tidy(fg_death, exponentiate=True, format="pandas")

pcm_tidy["cause"] = "PCM"
death_tidy["cause"] = "Death"

pd.concat([pcm_tidy, death_tidy])
```

This comparison reveals how covariates affect each cause differently. A covariate might
increase the incidence of progression to PCM while decreasing the incidence of death from
other causes (or vice versa) which both tell important clinical stories.

::: {.callout-note}
## Two kinds of hazard, two kinds of question

A cause-specific Cox model answers "what drives the rate of this event among those still at
risk?", while the Fine-Gray model answers "what drives the cumulative probability of this
event?". They can point in different directions for the same covariate, and both can be
correct. Decide which question you are asking before choosing the model.
:::

## Next steps

You can now estimate and model competing-risks data correctly.

- [Multi-state models](12-multi-state.qmd) generalize competing risks to settings where
subjects move between several states over time.
- [Visualizing survival](07-visualization.qmd) shows the grammar-of-graphics approach you
can reuse to plot cumulative incidence curves.
- Revisit [Survival data](03-survival-data.qmd) for the details of the multi-state response.


### Multi-state models

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

Competing risks describe a single transition from an initial state to one of several
end states. Multi-state models generalize this to processes where subjects move through
several states over time, possibly passing through intermediate ones. The classic example is
the illness-death model: a subject starts healthy, may become ill, and may die either
directly or after becoming ill. Greenwood estimates the transition and occupancy
probabilities for such models with the Aalen-Johansen estimator. This page shows how to
prepare multi-state data and interpret the results.

```{mermaid}
%%| label: fig-illness-death
%%| fig-cap: "The illness-death model. Subjects can move from health to illness, and reach death from either state."
flowchart LR
    H[Healthy] --> I[Ill]
    H --> D[Death]
    I --> D
```

## Preparing multi-state data

A multi-state analysis needs one row per period a subject spends in a state. Each row records
the interval `(start, stop]`, the state occupied during it, and the state transitioned to at
the end, or a censoring marker if none. A subject who becomes ill and then dies contributes
two rows: one for the healthy period ending in the illness transition, and one for the ill
period ending in death.

We build this from `mgus2`, treating progression as the intermediate "pcm" state and death
as the absorbing state.

```{python}
import greenwood as gw

mg = gw.load_dataset("mgus2", backend="polars")

start, stop, state, event = [], [], [], []
for i in range(len(mg)):
    pt, ft = mg["ptime"][i], mg["futime"][i]
    progressed, died = mg["pstat"][i] == 1, mg["death"][i] == 1
    if progressed and pt < ft:
        start += [0, pt]
        stop += [pt, ft]
        state += ["mgus", "pcm"]
        event += ["pcm", "death" if died else None]
    else:
        start += [0]
        stop += [ft]
        state += ["mgus"]
        event += ["death" if died else ("pcm" if progressed else None)]

# Keep only positive-length intervals.
rows = [(a, b, s, e) for a, b, s, e in zip(start, stop, state, event) if b > a]
start, stop, state, event = map(list, zip(*rows))
```

The loop appends to four parallel lists, so a single interval is spread across `start[i]`,
`stop[i]`, `state[i]`, and `event[i]` at the same index. After the loop we drop any zero-length
intervals and unpack the surviving rows back into the four lists. Counting the entries tells us
how many intervals we ended up with, which is more than the number of subjects because
progressing subjects contribute two rows each.

```{python}
len(start)
```

Reading the four lists position by position is easier if we zip them back together and look at
the first several intervals as `(start, stop, state, event)` tuples.

```{python}
list(zip(start, stop, state, event))[:6]
```

Each tuple is one interval `(start, stop]`: the subject occupies `state` from just after
`start` up to and including `stop`, and at `stop` transitions to `event`. A `None` in the
`event` position marks a censored interval, one that ends without a transition. A subject who
only ever holds `mgus` produces a single row, while a subject who progresses produces two
consecutive rows, the second of which has state `pcm`, exactly the illness-death pattern the
model needs.

## Estimating transition and occupancy probabilities

The `MultiState` estimator forms the Aalen-Johansen product of the per-time transition
matrices and reports, at each time, the probability of occupying each state. This is the
multi-state generalization of the survival curve.

```{python}
mg = gw.load_dataset("mgus2", backend="polars")

start, stop, state, event = [], [], [], []
for i in range(len(mg)):
    pt, ft = mg["ptime"][i], mg["futime"][i]
    progressed, died = mg["pstat"][i] == 1, mg["death"][i] == 1
    if progressed and pt < ft:
        start += [0, pt]
        stop += [pt, ft]
        state += ["mgus", "pcm"]
        event += ["pcm", "death" if died else None]
    else:
        start += [0]
        stop += [ft]
        state += ["mgus"]
        event += ["death" if died else ("pcm" if progressed else None)]

# Keep only positive-length intervals.
rows = [(a, b, s, e) for a, b, s, e in zip(start, stop, state, event) if b > a]
start, stop, state, event = map(list, zip(*rows))

ms = gw.MultiState().fit(start, stop, state, event, states=("mgus", "pcm", "death"))
```

The fitted estimator holds an occupancy curve for each of the three states. In its data frame
the columns after `time` give the occupancy probability of `mgus`, `pcm`, and `death`
respectively, and the preview reports the full number of time points.

```{python}
ms.to_frame(format="polars")
```

Each column after `time` is the probability that a subject is in that state at that time.
The probabilities sum to one across states, because every subject is always in exactly one
state. Early on, almost everyone is in the initial `mgus` state; over time, mass flows into
`pcm` and, predominantly, `death`.

You can read the occupancy probabilities at specific times with `predict`.

```{python}
ms.predict([60, 120, 240], format="polars")
```

The full transition-probability matrices, giving the probability of moving from each state to
each other state over the interval from time zero, are available as the `transition_`
attribute for users who need them.

::: {.callout-note}
## Occupancy versus transition probabilities

Occupancy probabilities answer "what is the chance a subject is in state k at time t?". The
transition matrix answers the more detailed "given a subject in state j at the start, what is
the chance they are in state k at time t?". The occupancy probabilities are the initial-state
distribution multiplied through the transition matrix.
:::

::: {.callout-tip}
## Competing risks is a special case

A competing-risks model is a multi-state model with one initial state and several absorbing
states and no intermediate transitions. If that is your situation, use
[`AalenJohansen`](11-competing-risks.qmd) directly, which also gives standard errors for the
cumulative incidence.
:::

## Next steps

You can now estimate occupancy probabilities for general multi-state processes.

- [Competing risks](11-competing-risks.qmd) covers the important special case with standard
errors and the Fine-Gray regression model.
- [Survival data](03-survival-data.qmd) explains the counting-process response that
underlies multi-state data.
- The [Quick start](02-quick-start.qmd) links every topic together if you want a map of the
whole package.



## Evaluation

### Prediction performance

```{python}
#| echo: false
#| output: false
import great_docs as gd

gd.enable_tbl_preview(n_head=8, n_tail=3)
```

A model that fits the data in front of it is not necessarily a model that predicts well. To
judge a survival model as a predictor, or to compare competing models, you need metrics that
account for censoring. Greenwood provides two complementary ones: the concordance index,
which measures how well a model ranks subjects by risk, and the Brier score, which measures
how accurate its probability predictions are. This page explains both and shows how to
compute them for a fitted model or for any external risk score.

We start from the `lung` outcomes and a Cox model fit to them. This is the model whose
predictions we will judge.

```{python}
import greenwood as gw

lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```

The response `y` records 165 events among 228 subjects. We fit a Cox model to it, then set the
model aside: everything below needs only `y` and the risk scores and survival probabilities
that `cox` produces.

```{python}
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
```

## Discrimination: the concordance index

Discrimination is the ability to rank subjects correctly: to give higher risk scores to
subjects who go on to have the event sooner. The concordance index, or C-statistic, measures
exactly this. It is the probability that, for a randomly chosen comparable pair of subjects,
the one who failed first had the higher risk score. A value of 0.5 is no better than a coin
flip, and 1.0 is perfect ordering.

You can compute it for any risk score, which makes it useful for evaluating a model on new
data or for comparing an external score. Here we use the Cox model's linear predictor. Asking
for `type="lp"` returns one number per subject, where a larger value means higher predicted
risk.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

risk = cox.predict(type="lp")
risk[:8]
```

::: {.callout-tip}
## Linear predictors are centered

The linear predictor returned by `predict(type="lp")` is always centered around the mean of
the training covariates. This ensures consistent, numerically stable predictions across all
variable types: continuous, categorical, and mixed. The centering is applied uniformly, so you
always get mathematically sound risk scores, even in edge cases like categorical-only models
with all positive associations.
:::

These are the scores the index will rank against the observed event order. Because higher risk
should correspond to earlier failure, we pass a score where larger means higher risk, such as
this linear predictor or its exponential.

```{python}
gw.concordance_index(y, risk)
```

::: {.callout-note}
## Discrimination is only part of the story

A model can rank subjects well yet produce poorly calibrated probabilities, or vice versa.
The concordance index says nothing about whether predicted probabilities are numerically
accurate. That is what the Brier score adds.
:::

## Accuracy: the Brier score

The Brier score is the mean squared error between predicted survival probabilities and actual
outcomes, evaluated at a fixed time. Lower is better. Because outcomes are censored, the
score uses inverse-probability-of-censoring weighting so that censored subjects contribute
correctly rather than being ignored or mishandled.

To compute it you supply predicted survival probabilities for each subject at each evaluation
time. Start by asking the Cox model for survival curves at three fixed times.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

times = [180, 365, 540]
surv = cox.predict(lung[["age", "sex"]], type="survival", times=times, format="polars")
surv.shape
```

The returned frame has a leading `time` column and one column per subject, so its shape is 3
rows (the requested times) by 229 columns (228 subjects plus the time column). The Brier score
wants the transpose of this: an array with one row per subject and one column per time. We
drop the `time` column and transpose to get there.

```{python}
probs = surv.drop("time").to_numpy().T
probs.shape
```

With the probabilities in the expected `(n_subjects, n_times)` shape, we can score them.

```{python}
gw.brier_score(y, probs, times)
```

The result is one Brier score per time. As a rule of thumb, a useful model should beat the
Brier score of a model that ignores covariates, and a score approaching 0.25 at a time where
about half the subjects have had the event is close to uninformative.

To summarize accuracy across a range of times in a single number, integrate the Brier score
over time.

```{python}
gw.integrated_brier_score(y, probs, times)
```

::: {.callout-tip}
## Evaluate on data the model did not see

Computing these metrics on the same data used to fit the model gives optimistic results.
For an honest assessment, evaluate on a held-out test set or use cross-validation, fitting on
one part of the data and scoring predictions on another.
:::

## Cross-validation

Rather than splitting the data by hand, `cross_validate` runs k-fold cross-validation for
you: it fits a fresh copy of the model on each set of training folds and scores its
predictions on the held-out fold. Pass an unfitted model, the response, and the covariates.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

gw.cross_validate(gw.CoxPH(), y, lung[["age", "sex"]], k=5, metric="concordance", seed=1)
```

The result reports the score on each fold plus their mean and standard deviation. The
cross-validated concordance is usually a little lower than the value computed on the training
data, which is exactly the optimism we are trying to remove. Set `metric="brier"` with a list
of `times` to cross-validate the integrated Brier score instead, and the same call works for
an `AFT` model.

## Calibration

Discrimination and the Brier score both summarize a model in a single number. Calibration
asks a more detailed question: at a chosen time, do the predicted survival probabilities match
what actually happens? `calibration` groups subjects into bins by their predicted survival at
a horizon, then compares the mean prediction in each bin against the observed survival, a
Kaplan-Meier estimate for that bin.

```{python}
lung = gw.load_dataset("lung", backend="polars")

y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])

surv_at_year = cox.predict(
    lung[["age", "sex"]], type="survival", times=[365.0], format="polars"
)

predicted = surv_at_year.slice(0, 1).drop("time").to_numpy().flatten()

gw.calibration(y, predicted, 365.0, n_bins=5, format="polars")
```

Each row is a bin, from the lowest predicted survival to the highest. A well-calibrated model
has `observed` close to `predicted` in every bin, with the observed values inside their
confidence limits. Systematic gaps, for example predictions consistently higher than observed,
signal miscalibration even when discrimination is good.

::: {.callout-tip}
## Assess calibration out of sample too

Like the other metrics, calibration is optimistic on the training data. For a fair picture,
compute it on held-out predictions, for instance the test-fold predictions from a
cross-validation split.
:::

## Next steps

You can now quantify how well a survival model discriminates and predicts.

- [Cox model diagnostics](09-cox-diagnostics.qmd) covers the assumption checks and residuals
that complement these performance metrics.
- [Parametric survival models](10-parametric-models.qmd) and
[Cox regression](08-cox-regression.qmd) produce the predictions these metrics evaluate.
- The [Quick start](02-quick-start.qmd) ties the whole workflow together.