restricted_mean_survival_time(..., return_variance=True) returns variance of the truncated RV, not the sampling variance — silently produces wrong p-values

[sebasj13]

Summary

lifelines.utils.restricted_mean_survival_time(model, t, return_variance=True) returns a "variance" that is not the sampling variance of the RMST estimator. It returns

E[min(T,t)²] − E[min(T,t)]²

i.e. the variance of the truncated survival-time random variable — a property of the survival distribution, not of the estimator. The sampling variance needed for confidence intervals or a 2-sample z-test follows the Greenwood-based formula in Klein & Moeschberger (2003, eq. 4.4.4), as implemented in survRM2::rmst2.

The docstring just says "variance of the RMST", so anyone using return_variance=True for inference will silently get a grossly inflated standard error and conservative p-values.

Source

lifelines/utils/__init__.py, lines ~248–253 on master:

mean = _expected_value_of_survival_up_to_t(model_or_survival_function, t)
if return_variance:
    sq = _expected_value_of_survival_squared_up_to_t(model_or_survival_function, t)
    return (mean, sq - mean**2)

sq - mean**2 is the variance of the truncated RV, not the sampling variance of the estimator.

Minimal reproduction (uses only the bundled waltons dataset)

import numpy as np
from lifelines import KaplanMeierFitter
from lifelines.datasets import load_waltons
from lifelines.utils import restricted_mean_survival_time

df = load_waltons()
ix = df["group"] == "miR-137"
kmf = KaplanMeierFitter().fit(df.loc[ix]["T"], df.loc[ix]["E"])

mean, var = restricted_mean_survival_time(kmf, t=10, return_variance=True)
print(f"RMST = {mean:.4f}")
print(f"reported 'variance' = {var:.4f}, implied SE = {np.sqrt(var):.4f}")

Output:

RMST = 9.7941
reported 'variance' = 0.5164, implied SE = 0.7186

Reference values from survRM2::rmst2 (and from PR #1526's own test fixture for this exact case): RMST = 9.794, SE = 0.123. The lifelines variance is ~34× too large; the implied SE is ~6× too large.

Effect on inference

Two-sample z-test of RMST_A − RMST_B on waltons at τ=10:

variance source	p-value
lifelines current return_variance=True	0.81
Greenwood (survRM2 / Klein-Moeschberger)	0.15

For real survival data on a months-or-years scale, the inflation is large enough to hide clear differences between visibly divergent KM curves.

Suggested resolution

PR #1526 already implements the Greenwood-based variance and a difference_in_restricted_mean_survival_time 2-sample test; reviving / merging that would resolve this. If the existing return_variance keyword is kept for backward compatibility, the docstring should at minimum spell out that the returned value is not suitable for hypothesis testing or confidence intervals.

Happy to help with either path if useful.