Mercurial > repos > fubar > lifelines_km_cph_tool
changeset 2:dd5e65893cb8 draft default tip
add survival and collapsed life table outputs suggested by Wolfgang
| author | fubar |
|---|---|
| date | Thu, 10 Aug 2023 22:52:45 +0000 |
| parents | 232b874046a7 |
| children | |
| files | lifelines_tool/README.md lifelines_tool/agepartialrossi.png lifelines_tool/lifelineskmcph.xml lifelines_tool/parolepartialrossi.png lifelines_tool/plotlykm.py lifelines_tool/run_log.txt lifelines_tool/test-data/readme_sample |
| diffstat | 7 files changed, 297 insertions(+), 21 deletions(-) [+] |
line wrap: on
line diff
--- a/lifelines_tool/README.md Thu Aug 10 07:15:22 2023 +0000 +++ b/lifelines_tool/README.md Thu Aug 10 22:52:45 2023 +0000 @@ -1,28 +1,40 @@ # lifelines_tool - lifelines statistical package wrapped as a Galaxy tool. -## Galaxy tool to run failure time models using lifelines +A Galaxy tool for right censored failure time data. + +Provides Kaplan-Meier plots with confidence intervals, and optional Cox proportional hazards models -## Install to your Galaxy server from the toolshed - search for lifelines_km_cph_tool owned by fubar2 +Uses the [lifelines](https://lifelines.readthedocs.io/en/latest/index.html) package. + +### Install to your Galaxy server from the toolshed - search for lifelines_km_cph_tool owned by fubar2 -### More at https://lazarus.name/demo/ +More at https://lazarus.name/demo/ + +#### Using the Rossi sample recidivism data from lifelines: -#### Using the Rossi sample input data from lifelines, tool outputs include: +Runs Kaplan-Meier and generates a plot. Optional grouping variable. + +Plots show confidence intervals  -and - -and + +If 2 groups, runs a log-rank test for difference. +  - -Runs Kaplan-Meier and generates a plot. Optional grouping variable. -If 2 groups, runs a log-rank test for difference. -Plots show confidence intervals +If a comma separated list (for example: prio, age, race, mar, fin) of covariate column names is provided, +a Cox proportional hazards model is run, the assumption of proportionality is tested, and +recommendations made. + -If a list of covariate column names is provided, these are used in a -Cox Proportional Hazards model with tests for proportionality. +Also included are partial plots for each covariate like these from the Rossi recidivism lifelines sample data +used in the tool test. -Should work with any tabular data with the required columns - time and status for observations. + + + + +Uses pandas read_csv with tabular delimiters so should work with any tabular data with the required columns - time and status for observations. Issues to https://github.com/fubar2/lifelines_tool please. Autogenerated so pull requests are possibly meaningless but regeneration of a new version should work.
--- a/lifelines_tool/lifelineskmcph.xml Thu Aug 10 07:15:22 2023 +0000 +++ b/lifelines_tool/lifelineskmcph.xml Thu Aug 10 22:52:45 2023 +0000 @@ -1,6 +1,6 @@ <tool name="lifelineskmcph" id="lifelineskmcph" version="0.01"> <!--Source in git at: https://github.com/fubar2/galaxy_tf_overlay--> - <!--Created by toolfactory@galaxy.org at 10/08/2023 15:48:43 using the Galaxy Tool Factory.--> + <!--Created by toolfactory@galaxy.org at 10/08/2023 21:59:53 using the Galaxy Tool Factory.--> <description>Lifelines KM and optional Cox PH models</description> <requirements> <requirement version="1.5.3" type="package">pandas</requirement> @@ -104,11 +104,11 @@ else: colsok = (args.time in defaultcols) and (args.status in defaultcols) if colsok: - sys.stderr.write('replacing first row of data derived header %s with %s' % (testcols, defaultcols)) + print('Replacing first row of data derived header %s with %s' % (testcols, defaultcols)) df.columns = defaultcols else: sys.stderr.write('## CRITICAL USAGE ERROR (not a bug!): time %s and status %s do not match anything in the file header, supplied header or automatic default column names %s' % (args.time, args.status, defaultcols)) -print('## Lifelines tool starting.\nUsing data header =', df.columns, 'time column =', args.time, 'status column =', args.status) +print('## Lifelines tool\nInput data header =', df.columns, 'time column =', args.time, 'status column =', args.status) os.makedirs(args.image_dir, exist_ok=True) fig, ax = plt.subplots() if args.group > '': @@ -136,6 +136,24 @@ ax.set_title(args.title) fig.savefig(os.path.join(args.image_dir,'KM_%s.png' % args.title)) print('#### No grouping variable, so no log rank or other Kaplan-Meier statistical output is available') +survdf = lifelines.utils.survival_table_from_events(df[args.time], df[args.status]) +lifedf = lifelines.utils.survival_table_from_events(df[args.time], df[args.status], collapse=True) +print("#### Survival table using time %s and event %s" % (args.time, args.status)) +with pd.option_context('display.max_rows', None, + 'display.max_columns', None, + 'display.precision', 3, + ): + print(survdf) +print("#### Life table using time %s and event %s" % (args.time, args.status)) +with pd.option_context('display.max_rows', None, + 'display.max_columns', None, + 'display.precision', 3, + ): + print(lifedf) +outpath = os.path.join(args.image_dir,'survival_table.tabular') +survdf.to_csv(outpath, sep='\t') +outpath = os.path.join(args.image_dir,'life_table.tabular') +lifedf.to_csv(outpath, sep='\t') if len(args.cphcols) > 0: fig, ax = plt.subplots() ax.set_title('Cox-PH model: %s' % args.title) @@ -153,7 +171,7 @@ cph.fit(cphdf, duration_col=args.time, event_col=args.status) cph.print_summary() for i, cov in enumerate(colsdf.columns): - if ucolcounts[i] > 10: + if ucolcounts[i] > 10: # a hack - assume categories are sparse - if not imaginary quintiles will have to do v = pd.Series.tolist(cphdf[cov].quantile(QVALS)) vdt = df.dtypes[cov] if vdt == 'int64':
--- a/lifelines_tool/plotlykm.py Thu Aug 10 07:15:22 2023 +0000 +++ b/lifelines_tool/plotlykm.py Thu Aug 10 22:52:45 2023 +0000 @@ -98,6 +98,24 @@ ax.set_title(args.title) fig.savefig(os.path.join(args.image_dir,'KM_%s.png' % args.title)) print('#### No grouping variable, so no log rank or other Kaplan-Meier statistical output is available') +survdf = lifelines.utils.survival_table_from_events(df[args.time], df[args.status]) +lifedf = lifelines.utils.survival_table_from_events(df[args.time], df[args.status], collapse=True) +print("Survival table using time %s and event %s" % (args.time, args.status)) +with pd.option_context('display.max_rows', None, + 'display.max_columns', None, + 'display.precision', 3, + ): + print(survdf) +print("Life table using time %s and event %s" % (args.time, args.status)) +with pd.option_context('display.max_rows', None, + 'display.max_columns', None, + 'display.precision', 3, + ): + print(lifedf) +outpath = os.path.join(args.image_dir,'survival_table.tabular') +survdf.to_csv(outpath, sep='\t') +outpath = os.path.join(args.image_dir,'life_table.tabular') +lifedf.to_csv(outpath, sep='\t') if len(args.cphcols) > 0: fig, ax = plt.subplots() ax.set_title('Cox-PH model: %s' % args.title)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lifelines_tool/run_log.txt Thu Aug 10 22:52:45 2023 +0000 @@ -0,0 +1,168 @@ +## Lifelines tool +Input data header = Index(['Unnamed: 0', 'week', 'arrest', 'fin', 'age', 'race', 'wexp', 'mar', + 'paro', 'prio'], + dtype='object') time column = week status column = arrest +#### No grouping variable, so no log rank or other Kaplan-Meier statistical output is available +Survival table using time week and event arrest + removed observed censored entrance at_risk +event_at +0.0 0 0 0 432 432 +1.0 1 1 0 0 432 +2.0 1 1 0 0 431 +3.0 1 1 0 0 430 +4.0 1 1 0 0 429 +5.0 1 1 0 0 428 +6.0 1 1 0 0 427 +7.0 1 1 0 0 426 +8.0 5 5 0 0 425 +9.0 2 2 0 0 420 +10.0 1 1 0 0 418 +11.0 2 2 0 0 417 +12.0 2 2 0 0 415 +13.0 1 1 0 0 413 +14.0 3 3 0 0 412 +15.0 2 2 0 0 409 +16.0 2 2 0 0 407 +17.0 3 3 0 0 405 +18.0 3 3 0 0 402 +19.0 2 2 0 0 399 +20.0 5 5 0 0 397 +21.0 2 2 0 0 392 +22.0 1 1 0 0 390 +23.0 1 1 0 0 389 +24.0 4 4 0 0 388 +25.0 3 3 0 0 384 +26.0 3 3 0 0 381 +27.0 2 2 0 0 378 +28.0 2 2 0 0 376 +30.0 2 2 0 0 374 +31.0 1 1 0 0 372 +32.0 2 2 0 0 371 +33.0 2 2 0 0 369 +34.0 2 2 0 0 367 +35.0 4 4 0 0 365 +36.0 3 3 0 0 361 +37.0 4 4 0 0 358 +38.0 1 1 0 0 354 +39.0 2 2 0 0 353 +40.0 4 4 0 0 351 +42.0 2 2 0 0 347 +43.0 4 4 0 0 345 +44.0 2 2 0 0 341 +45.0 2 2 0 0 339 +46.0 4 4 0 0 337 +47.0 1 1 0 0 333 +48.0 2 2 0 0 332 +49.0 5 5 0 0 330 +50.0 3 3 0 0 325 +52.0 322 4 318 0 322 +Life table using time week and event arrest + removed observed censored at_risk +event_at +(-0.001, 13.844] 20 20 0 432 +(13.844, 27.687] 36 36 0 412 +(27.687, 41.531] 29 29 0 376 +(41.531, 55.374] 347 29 318 347 +### Lifelines test of Proportional Hazards results with prio, age, race, paro, mar, fin as covariates on test +<lifelines.CoxPHFitter: fitted with 432 total observations, 318 right-censored observations> + duration col = 'week' + event col = 'arrest' + baseline estimation = breslow + number of observations = 432 +number of events observed = 114 + partial log-likelihood = -659.00 + time fit was run = 2023-08-10 11:57:10 UTC + +--- + coef exp(coef) se(coef) coef lower 95% coef upper 95% exp(coef) lower 95% exp(coef) upper 95% +covariate +prio 0.10 1.10 0.03 0.04 0.15 1.04 1.16 +age -0.06 0.94 0.02 -0.10 -0.02 0.90 0.98 +race 0.32 1.38 0.31 -0.28 0.92 0.75 2.52 +paro -0.09 0.91 0.20 -0.47 0.29 0.62 1.34 +mar -0.48 0.62 0.38 -1.22 0.25 0.30 1.29 +fin -0.38 0.68 0.19 -0.75 -0.00 0.47 1.00 + + cmp to z p -log2(p) +covariate +prio 0.00 3.53 <0.005 11.26 +age 0.00 -2.95 <0.005 8.28 +race 0.00 1.04 0.30 1.75 +paro 0.00 -0.46 0.65 0.63 +mar 0.00 -1.28 0.20 2.32 +fin 0.00 -1.98 0.05 4.40 +--- +Concordance = 0.63 +Partial AIC = 1330.00 +log-likelihood ratio test = 32.77 on 6 df +-log2(p) of ll-ratio test = 16.39 + + + Bootstrapping lowess lines. May take a moment... + + + Bootstrapping lowess lines. May take a moment... + +The ``p_value_threshold`` is set at 0.01. Even under the null hypothesis of no violations, some +covariates will be below the threshold by chance. This is compounded when there are many covariates. +Similarly, when there are lots of observations, even minor deviances from the proportional hazard +assumption will be flagged. + +With that in mind, it's best to use a combination of statistical tests and visual tests to determine +the most serious violations. Produce visual plots using ``check_assumptions(..., show_plots=True)`` +and looking for non-constant lines. See link [A] below for a full example. + +<lifelines.StatisticalResult: proportional_hazard_test> + null_distribution = chi squared +degrees_of_freedom = 1 + model = <lifelines.CoxPHFitter: fitted with 432 total observations, 318 right-censored observations> + test_name = proportional_hazard_test + +--- + test_statistic p -log2(p) +age km 6.99 0.01 6.93 + rank 7.40 0.01 7.26 +fin km 0.02 0.90 0.15 + rank 0.01 0.91 0.13 +mar km 1.64 0.20 2.32 + rank 1.80 0.18 2.48 +paro km 0.06 0.81 0.31 + rank 0.07 0.79 0.34 +prio km 0.92 0.34 1.57 + rank 0.88 0.35 1.52 +race km 1.70 0.19 2.38 + rank 1.68 0.19 2.36 + + +1. Variable 'age' failed the non-proportional test: p-value is 0.0065. + + Advice 1: the functional form of the variable 'age' might be incorrect. That is, there may be +non-linear terms missing. The proportional hazard test used is very sensitive to incorrect +functional forms. See documentation in link [D] below on how to specify a functional form. + + Advice 2: try binning the variable 'age' using pd.cut, and then specify it in `strata=['age', +...]` in the call in `.fit`. See documentation in link [B] below. + + Advice 3: try adding an interaction term with your time variable. See documentation in link [C] +below. + + + Bootstrapping lowess lines. May take a moment... + + + Bootstrapping lowess lines. May take a moment... + + + Bootstrapping lowess lines. May take a moment... + + + Bootstrapping lowess lines. May take a moment... + + +--- +[A] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html +[B] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Bin-variable-and-stratify-on-it +[C] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Introduce-time-varying-covariates +[D] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Modify-the-functional-form +[E] https://lifelines.readthedocs.io/en/latest/jupyter_notebooks/Proportional%20hazard%20assumption.html#Stratification +
--- a/lifelines_tool/test-data/readme_sample Thu Aug 10 07:15:22 2023 +0000 +++ b/lifelines_tool/test-data/readme_sample Thu Aug 10 22:52:45 2023 +0000 @@ -1,5 +1,5 @@ -## Lifelines tool starting. -Using data header = Index(['Unnamed: 0', 'week', 'arrest', 'fin', 'age', 'race', 'wexp', 'mar', +## Lifelines tool +Input data header = Index(['Unnamed: 0', 'week', 'arrest', 'fin', 'age', 'race', 'wexp', 'mar', 'paro', 'prio'], dtype='object') time column = week status column = arrest Logrank test for race - 0 vs 1 @@ -14,6 +14,66 @@ --- test_statistic p -log2(p) 0.58 0.45 1.16 +#### Survival table using time week and event arrest + removed observed censored entrance at_risk +event_at +0.0 0 0 0 432 432 +1.0 1 1 0 0 432 +2.0 1 1 0 0 431 +3.0 1 1 0 0 430 +4.0 1 1 0 0 429 +5.0 1 1 0 0 428 +6.0 1 1 0 0 427 +7.0 1 1 0 0 426 +8.0 5 5 0 0 425 +9.0 2 2 0 0 420 +10.0 1 1 0 0 418 +11.0 2 2 0 0 417 +12.0 2 2 0 0 415 +13.0 1 1 0 0 413 +14.0 3 3 0 0 412 +15.0 2 2 0 0 409 +16.0 2 2 0 0 407 +17.0 3 3 0 0 405 +18.0 3 3 0 0 402 +19.0 2 2 0 0 399 +20.0 5 5 0 0 397 +21.0 2 2 0 0 392 +22.0 1 1 0 0 390 +23.0 1 1 0 0 389 +24.0 4 4 0 0 388 +25.0 3 3 0 0 384 +26.0 3 3 0 0 381 +27.0 2 2 0 0 378 +28.0 2 2 0 0 376 +30.0 2 2 0 0 374 +31.0 1 1 0 0 372 +32.0 2 2 0 0 371 +33.0 2 2 0 0 369 +34.0 2 2 0 0 367 +35.0 4 4 0 0 365 +36.0 3 3 0 0 361 +37.0 4 4 0 0 358 +38.0 1 1 0 0 354 +39.0 2 2 0 0 353 +40.0 4 4 0 0 351 +42.0 2 2 0 0 347 +43.0 4 4 0 0 345 +44.0 2 2 0 0 341 +45.0 2 2 0 0 339 +46.0 4 4 0 0 337 +47.0 1 1 0 0 333 +48.0 2 2 0 0 332 +49.0 5 5 0 0 330 +50.0 3 3 0 0 325 +52.0 322 4 318 0 322 +#### Life table using time week and event arrest + removed observed censored at_risk +event_at +(-0.001, 13.844] 20 20 0 432 +(13.844, 27.687] 36 36 0 412 +(27.687, 41.531] 29 29 0 376 +(41.531, 55.374] 347 29 318 347 ### Lifelines test of Proportional Hazards results with prio, age, race, paro, mar, fin as covariates on KM and CPH in lifelines test <lifelines.CoxPHFitter: fitted with 432 total observations, 318 right-censored observations> duration col = 'week' @@ -22,7 +82,7 @@ number of observations = 432 number of events observed = 114 partial log-likelihood = -659.00 - time fit was run = 2023-08-10 05:49:04 UTC + time fit was run = 2023-08-10 12:00:14 UTC --- coef exp(coef) se(coef) coef lower 95% coef upper 95% exp(coef) lower 95% exp(coef) upper 95%