diff --git a/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_level_and_officer_type.xlsx b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_level_and_officer_type.xlsx index 3d804bbc77..1dc4b6ea7e 100644 --- a/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_level_and_officer_type.xlsx +++ b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_level_and_officer_type.xlsx @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:af86c2c2af5c291c18c5d481681d6d316526b81806c8c8e898517e850160e6fd -size 12465 +oid sha256:80651d157772a292bf9617c86e2616d8165a20385ada6d85a5244aca9c55aa0c +size 21938 diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/analysis_hss_elements.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/analysis_hss_elements.py new file mode 100644 index 0000000000..76708f7c25 --- /dev/null +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/analysis_hss_elements.py @@ -0,0 +1,272 @@ +"""Produce plots to show the impact each the healthcare system (overall health impact) when running under different +scenarios (scenario_impact_of_healthsystem.py)""" + +import argparse +import textwrap +from pathlib import Path +from typing import Tuple + +import numpy as np +import pandas as pd +from matplotlib import pyplot as plt + +from tlo import Date +from tlo.analysis.utils import extract_results, make_age_grp_lookup, summarize + + +def apply(results_folder: Path, output_folder: Path, resourcefilepath: Path = None): + """Produce standard set of plots describing the effect of each TREATMENT_ID. + - We estimate the epidemiological impact as the EXTRA deaths that would occur if that treatment did not occur. + - We estimate the draw on healthcare system resources as the FEWER appointments when that treatment does not occur. + """ + + TARGET_PERIOD = (Date(2020, 1, 1), Date(2030, 12, 31)) + + # Definitions of general helper functions + make_graph_file_name = lambda stub: output_folder / f"{stub.replace('*', '_star_')}.png" # noqa: E731 + + _, age_grp_lookup = make_age_grp_lookup() + + def target_period() -> str: + """Returns the target period as a string of the form YYYY-YYYY""" + return "-".join(str(t.year) for t in TARGET_PERIOD) + + def get_parameter_names_from_scenario_file() -> Tuple[str]: + """Get the tuple of names of the scenarios from `Scenario` class used to create the results.""" + from scripts.comparison_of_horizontal_and_vertical_programs.scenario_hss_elements import ( + HSSElements, + ) + e = HSSElements() + return tuple(e._scenarios.keys()) + + def get_num_deaths(_df): + """Return total number of Deaths (total within the TARGET_PERIOD)""" + return pd.Series(data=len(_df.loc[pd.to_datetime(_df.date).between(*TARGET_PERIOD)])) + + def get_num_dalys(_df): + """Return total number of DALYS (Stacked) by label (total within the TARGET_PERIOD). + Throw error if not a record for every year in the TARGET PERIOD (to guard against inadvertently using + results from runs that crashed mid-way through the simulation. + """ + years_needed = [i.year for i in TARGET_PERIOD] + assert set(_df.year.unique()).issuperset(years_needed), "Some years are not recorded." + return pd.Series( + data=_df + .loc[_df.year.between(*years_needed)] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum().sum() + ) + + def set_param_names_as_column_index_level_0(_df): + """Set the columns index (level 0) as the param_names.""" + ordered_param_names_no_prefix = {i: x for i, x in enumerate(param_names)} + names_of_cols_level0 = [ordered_param_names_no_prefix.get(col) for col in _df.columns.levels[0]] + assert len(names_of_cols_level0) == len(_df.columns.levels[0]) + _df.columns = _df.columns.set_levels(names_of_cols_level0, level=0) + return _df + + def find_difference_relative_to_comparison(_ser: pd.Series, + comparison: str, + scaled: bool = False, + drop_comparison: bool = True, + ): + """Find the difference in the values in a pd.Series with a multi-index, between the draws (level 0) + within the runs (level 1), relative to where draw = `comparison`. + The comparison is `X - COMPARISON`.""" + return _ser \ + .unstack(level=0) \ + .apply(lambda x: (x - x[comparison]) / (x[comparison] if scaled else 1.0), axis=1) \ + .drop(columns=([comparison] if drop_comparison else [])) \ + .stack() + + def do_bar_plot_with_ci(_df, annotations=None, xticklabels_horizontal_and_wrapped=False, put_labels_in_legend=True): + """Make a vertical bar plot for each row of _df, using the columns to identify the height of the bar and the + extent of the error bar.""" + + substitute_labels = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' + + yerr = np.array([ + (_df['mean'] - _df['lower']).values, + (_df['upper'] - _df['mean']).values, + ]) + + xticks = {(i + 0.5): k for i, k in enumerate(_df.index)} + + # Define colormap (used only with option `put_labels_in_legend=True`) + cmap = plt.get_cmap("tab20") + rescale = lambda y: (y - np.min(y)) / (np.max(y) - np.min(y)) # noqa: E731 + colors = list(map(cmap, rescale(np.array(list(xticks.keys()))))) if put_labels_in_legend else None + + fig, ax = plt.subplots(figsize=(10, 5)) + ax.bar( + xticks.keys(), + _df['mean'].values, + yerr=yerr, + alpha=0.8, + ecolor='black', + color=colors, + capsize=10, + label=xticks.values() + ) + if annotations: + for xpos, ypos, text in zip(xticks.keys(), _df['upper'].values, annotations): + ax.text(xpos, ypos*1.15, text, horizontalalignment='center', rotation='vertical', fontsize='x-small') + ax.set_xticks(list(xticks.keys())) + + if put_labels_in_legend: + # Update xticks label with substitute labels + # Insert legend with updated labels that shows correspondence between substitute label and original label + xtick_values = [letter for letter, label in zip(substitute_labels, xticks.values())] + xtick_legend = [f'{letter}: {label}' for letter, label in zip(substitute_labels, xticks.values())] + h, legs = ax.get_legend_handles_labels() + ax.legend(h, xtick_legend, loc='center left', fontsize='small', bbox_to_anchor=(1, 0.5)) + ax.set_xticklabels(list(xtick_values)) + else: + if not xticklabels_horizontal_and_wrapped: + # xticklabels will be vertical and not wrapped + ax.set_xticklabels(list(xticks.values()), rotation=90) + else: + wrapped_labs = ["\n".join(textwrap.wrap(_lab, 20)) for _lab in xticks.values()] + ax.set_xticklabels(wrapped_labs) + + ax.grid(axis="y") + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + fig.tight_layout() + + return fig, ax + + # %% Define parameter names + param_names = get_parameter_names_from_scenario_file() + + # %% Quantify the health gains associated with all interventions combined. + + # Absolute Number of Deaths and DALYs + num_deaths = extract_results( + results_folder, + module='tlo.methods.demography', + key='death', + custom_generate_series=get_num_deaths, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + + num_dalys = extract_results( + results_folder, + module='tlo.methods.healthburden', + key='dalys_stacked', + custom_generate_series=get_num_dalys, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + + # %% Charts of total numbers of deaths / DALYS + num_dalys_summarized = summarize(num_dalys).loc[0].unstack().reindex(param_names) + num_deaths_summarized = summarize(num_deaths).loc[0].unstack().reindex(param_names) + + name_of_plot = f'Deaths, {target_period()}' + fig, ax = do_bar_plot_with_ci(num_deaths_summarized / 1e6) + ax.set_title(name_of_plot) + ax.set_ylabel('(Millions)') + fig.tight_layout() + ax.axhline(num_deaths_summarized.loc['Baseline', 'mean']/1e6, color='black', alpha=0.5) + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + name_of_plot = f'All Scenarios: DALYs, {target_period()}' + fig, ax = do_bar_plot_with_ci(num_dalys_summarized / 1e6) + ax.set_title(name_of_plot) + ax.set_ylabel('(Millions)') + ax.axhline(num_dalys_summarized.loc['Baseline', 'mean']/1e6, color='black', alpha=0.5) + fig.tight_layout() + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + + # %% Deaths and DALYS averted relative to Status Quo + num_deaths_averted = summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison( + num_deaths.loc[0], + comparison='Baseline') + ).T + ).iloc[0].unstack().reindex(param_names).drop(['Baseline']) + + pc_deaths_averted = 100.0 * summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison( + num_deaths.loc[0], + comparison='Baseline', + scaled=True) + ).T + ).iloc[0].unstack().reindex(param_names).drop(['Baseline']) + + num_dalys_averted = summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison( + num_dalys.loc[0], + comparison='Baseline') + ).T + ).iloc[0].unstack().reindex(param_names).drop(['Baseline']) + + pc_dalys_averted = 100.0 * summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison( + num_dalys.loc[0], + comparison='Baseline', + scaled=True) + ).T + ).iloc[0].unstack().reindex(param_names).drop(['Baseline']) + + # DEATHS + name_of_plot = f'Additional Deaths Averted vs Baseline, {target_period()}' + fig, ax = do_bar_plot_with_ci( + num_deaths_averted.clip(lower=0.0), + annotations=[ + f"{round(row['mean'], 0)} ({round(row['lower'], 1)}-{round(row['upper'], 1)}) %" + for _, row in pc_deaths_averted.clip(lower=0.0).iterrows() + ] + ) + ax.set_title(name_of_plot) + ax.set_ylabel('Additional Deaths Averted') + fig.tight_layout() + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + # DALYS + name_of_plot = f'Additional DALYs Averted vs Baseline, {target_period()}' + fig, ax = do_bar_plot_with_ci( + (num_dalys_averted / 1e6).clip(lower=0.0), + annotations=[ + f"{round(row['mean'])} ({round(row['lower'], 1)}-{round(row['upper'], 1)}) %" + for _, row in pc_dalys_averted.clip(lower=0.0).iterrows() + ] + ) + ax.set_title(name_of_plot) + ax.set_ylabel('Additional DALYS Averted \n(Millions)') + fig.tight_layout() + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + # todo: Neaten graphs + # todo: Graph showing difference broken down by disease (this can be cribbed from the calcs about wealth from the + # third set of analyses in the overview paper). + # todo: other metrics of health + # todo: other graphs, broken down by age/sex (this can also be cribbed from overview paper stuff) + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("results_folder", type=Path) # outputs/horizontal_and_vertical_programs-2024-05-16 + args = parser.parse_args() + + apply( + results_folder=args.results_folder, + output_folder=args.results_folder, + resourcefilepath=Path('./resources') + ) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/analysis_vertical_programs_with_and_without_hss.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/analysis_vertical_programs_with_and_without_hss.py new file mode 100644 index 0000000000..f0dd083d97 --- /dev/null +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/analysis_vertical_programs_with_and_without_hss.py @@ -0,0 +1,363 @@ +"""Produce plots to show the impact each the healthcare system (overall health impact) when running under different +scenarios (scenario_impact_of_healthsystem.py)""" + +import argparse +import textwrap +from pathlib import Path +from typing import Tuple + +import numpy as np +import pandas as pd +from matplotlib import pyplot as plt + +from tlo import Date +from tlo.analysis.utils import extract_results, make_age_grp_lookup, summarize + + +def apply(results_folder: Path, output_folder: Path, resourcefilepath: Path = None): + """Produce standard set of plots describing the effect of each TREATMENT_ID. + - We estimate the epidemiological impact as the EXTRA deaths that would occur if that treatment did not occur. + - We estimate the draw on healthcare system resources as the FEWER appointments when that treatment does not occur. + """ + + TARGET_PERIOD = (Date(2020, 1, 1), Date(2030, 12, 31)) + + # Definitions of general helper functions + make_graph_file_name = lambda stub: output_folder / f"{stub.replace('*', '_star_')}.png" # noqa: E731 + + _, age_grp_lookup = make_age_grp_lookup() + + def target_period() -> str: + """Returns the target period as a string of the form YYYY-YYYY""" + return "-".join(str(t.year) for t in TARGET_PERIOD) + + def get_parameter_names_from_scenario_file() -> Tuple[str]: + """Get the tuple of names of the scenarios from `Scenario` class used to create the results.""" + from scripts.comparison_of_horizontal_and_vertical_programs.scenario_vertical_programs_with_and_without_hss import ( + HTMWithAndWithoutHSS, + ) + e = HTMWithAndWithoutHSS() + return tuple(e._scenarios.keys()) + + def get_num_deaths(_df): + """Return total number of Deaths (total within the TARGET_PERIOD)""" + return pd.Series(data=len(_df.loc[pd.to_datetime(_df.date).between(*TARGET_PERIOD)])) + + def get_num_dalys(_df): + """Return total number of DALYS (Stacked) by label (total within the TARGET_PERIOD). + Throw error if not a record for every year in the TARGET PERIOD (to guard against inadvertently using + results from runs that crashed mid-way through the simulation. + """ + years_needed = [i.year for i in TARGET_PERIOD] + assert set(_df.year.unique()).issuperset(years_needed), "Some years are not recorded." + return pd.Series( + data=_df + .loc[_df.year.between(*years_needed)] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum().sum() + ) + + def set_param_names_as_column_index_level_0(_df): + """Set the columns index (level 0) as the param_names.""" + ordered_param_names_no_prefix = {i: x for i, x in enumerate(param_names)} + names_of_cols_level0 = [ordered_param_names_no_prefix.get(col) for col in _df.columns.levels[0]] + assert len(names_of_cols_level0) == len(_df.columns.levels[0]) + _df.columns = _df.columns.set_levels(names_of_cols_level0, level=0) + return _df + + def find_difference_relative_to_comparison_series( + _ser: pd.Series, + comparison: str, + scaled: bool = False, + drop_comparison: bool = True, + ): + """Find the difference in the values in a pd.Series with a multi-index, between the draws (level 0) + within the runs (level 1), relative to where draw = `comparison`. + The comparison is `X - COMPARISON`.""" + return _ser \ + .unstack(level=0) \ + .apply(lambda x: (x - x[comparison]) / (x[comparison] if scaled else 1.0), axis=1) \ + .drop(columns=([comparison] if drop_comparison else [])) \ + .stack() + + def find_difference_relative_to_comparison_series_dataframe(_df: pd.DataFrame, **kwargs): + """Apply `find_difference_relative_to_comparison_series` to each row in a dataframe""" + return pd.concat({ + _idx: find_difference_relative_to_comparison_series(row, **kwargs) + for _idx, row in _df.iterrows() + }, axis=1).T + + def do_bar_plot_with_ci(_df, annotations=None, xticklabels_horizontal_and_wrapped=False, put_labels_in_legend=True): + """Make a vertical bar plot for each row of _df, using the columns to identify the height of the bar and the + extent of the error bar.""" + + substitute_labels = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' + + yerr = np.array([ + (_df['mean'] - _df['lower']).values, + (_df['upper'] - _df['mean']).values, + ]) + + xticks = {(i + 0.5): k for i, k in enumerate(_df.index)} + + # Define colormap (used only with option `put_labels_in_legend=True`) + cmap = plt.get_cmap("tab20") + rescale = lambda y: (y - np.min(y)) / (np.max(y) - np.min(y)) # noqa: E731 + colors = list(map(cmap, rescale(np.array(list(xticks.keys()))))) if put_labels_in_legend else None + + fig, ax = plt.subplots(figsize=(10, 5)) + ax.bar( + xticks.keys(), + _df['mean'].values, + yerr=yerr, + alpha=0.8, + ecolor='black', + color=colors, + capsize=10, + label=xticks.values() + ) + if annotations: + for xpos, ypos, text in zip(xticks.keys(), _df['upper'].values, annotations): + ax.text(xpos, ypos*1.15, text, horizontalalignment='center', rotation='vertical', fontsize='x-small') + ax.set_xticks(list(xticks.keys())) + + if put_labels_in_legend: + # Update xticks label with substitute labels + # Insert legend with updated labels that shows correspondence between substitute label and original label + xtick_values = [letter for letter, label in zip(substitute_labels, xticks.values())] + xtick_legend = [f'{letter}: {label}' for letter, label in zip(substitute_labels, xticks.values())] + h, legs = ax.get_legend_handles_labels() + ax.legend(h, xtick_legend, loc='center left', fontsize='small', bbox_to_anchor=(1, 0.5)) + ax.set_xticklabels(list(xtick_values)) + else: + if not xticklabels_horizontal_and_wrapped: + # xticklabels will be vertical and not wrapped + ax.set_xticklabels(list(xticks.values()), rotation=90) + else: + wrapped_labs = ["\n".join(textwrap.wrap(_lab, 20)) for _lab in xticks.values()] + ax.set_xticklabels(wrapped_labs) + + ax.grid(axis="y") + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + fig.tight_layout() + + return fig, ax + + # %% Define parameter names + param_names = get_parameter_names_from_scenario_file() + + # %% Quantify the health gains associated with all interventions combined. + + # Absolute Number of Deaths and DALYs + num_deaths = extract_results( + results_folder, + module='tlo.methods.demography', + key='death', + custom_generate_series=get_num_deaths, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + + num_dalys = extract_results( + results_folder, + module='tlo.methods.healthburden', + key='dalys_stacked', + custom_generate_series=get_num_dalys, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + + # %% Charts of total numbers of deaths / DALYS + num_dalys_summarized = summarize(num_dalys).loc[0].unstack().reindex(param_names) + num_deaths_summarized = summarize(num_deaths).loc[0].unstack().reindex(param_names) + + name_of_plot = f'Deaths, {target_period()}' + fig, ax = do_bar_plot_with_ci(num_deaths_summarized / 1e6) + ax.set_title(name_of_plot) + ax.set_ylabel('(Millions)') + fig.tight_layout() + ax.axhline(num_deaths_summarized.loc['Baseline', 'mean']/1e6, color='black', alpha=0.5) + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + name_of_plot = f'All Scenarios: DALYs, {target_period()}' + fig, ax = do_bar_plot_with_ci(num_dalys_summarized / 1e6) + ax.set_title(name_of_plot) + ax.set_ylabel('(Millions)') + ax.axhline(num_dalys_summarized.loc['Baseline', 'mean']/1e6, color='black', alpha=0.5) + fig.tight_layout() + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + + # %% Deaths and DALYS averted relative to Status Quo + num_deaths_averted = summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison_series( + num_deaths.loc[0], + comparison='Baseline') + ).T + ).iloc[0].unstack().reindex(param_names).drop(['Baseline']) + + pc_deaths_averted = 100.0 * summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison_series( + num_deaths.loc[0], + comparison='Baseline', + scaled=True) + ).T + ).iloc[0].unstack().reindex(param_names).drop(['Baseline']) + + num_dalys_averted = summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison_series( + num_dalys.loc[0], + comparison='Baseline') + ).T + ).iloc[0].unstack().reindex(param_names).drop(['Baseline']) + + pc_dalys_averted = 100.0 * summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison_series( + num_dalys.loc[0], + comparison='Baseline', + scaled=True) + ).T + ).iloc[0].unstack().reindex(param_names).drop(['Baseline']) + + # DEATHS + name_of_plot = f'Additional Deaths Averted vs Baseline, {target_period()}' + fig, ax = do_bar_plot_with_ci( + num_deaths_averted.clip(lower=0.0), + annotations=[ + f"{round(row['mean'], 0)} ({round(row['lower'], 1)}-{round(row['upper'], 1)}) %" + for _, row in pc_deaths_averted.clip(lower=0.0).iterrows() + ] + ) + ax.set_title(name_of_plot) + ax.set_ylabel('Additional Deaths Averted vs Baseline') + fig.tight_layout() + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + # DALYS + name_of_plot = f'DALYs Averted vs Baseline, {target_period()}' + fig, ax = do_bar_plot_with_ci( + (num_dalys_averted / 1e6).clip(lower=0.0), + annotations=[ + f"{round(row['mean'])} ({round(row['lower'], 1)}-{round(row['upper'], 1)}) %" + for _, row in pc_dalys_averted.clip(lower=0.0).iterrows() + ] + ) + ax.set_title(name_of_plot) + ax.set_ylabel('Additional DALYS Averted vs Baseline \n(Millions)') + fig.tight_layout() + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + + # %% DALYS averted relative to Baseline - broken down by major cause (HIV, TB, MALARIA) + + def get_total_num_dalys_by_label(_df): + """Return the total number of DALYS in the TARGET_PERIOD by wealth and cause label.""" + y = _df \ + .loc[_df['year'].between(*[d.year for d in TARGET_PERIOD])] \ + .drop(columns=['date', 'year', 'li_wealth']) \ + .sum(axis=0) + + # define course cause mapper for HIV, TB, MALARIA and OTHER + causes = { + 'AIDS': 'HIV/AIDS', + 'TB (non-AIDS)': 'TB', + 'Malaria': 'Malaria', + '': 'Other', # defined in order to use this dict to determine ordering of the causes in output + } + causes_relabels = y.index.map(causes).fillna('Other') + + return y.groupby(by=causes_relabels).sum()[list(causes.values())] + + total_num_dalys_by_label_results = extract_results( + results_folder, + module="tlo.methods.healthburden", + key="dalys_by_wealth_stacked_by_age_and_time", + custom_generate_series=get_total_num_dalys_by_label, + do_scaling=True, + ).pipe(set_param_names_as_column_index_level_0) + + total_num_dalys_by_label_results_averted_vs_baseline = summarize( + -1.0 * find_difference_relative_to_comparison_series_dataframe( + total_num_dalys_by_label_results, + comparison='Baseline' + ), + only_mean=True + ) + + # Check that when we sum across the causes, we get the same total as calculated when we didn't split by cause. + assert ( + (total_num_dalys_by_label_results_averted_vs_baseline.sum(axis=0).sort_index() + - num_dalys_averted['mean'].sort_index() + ) < 1e-6 + ).all() + + # Make a separate plot for the scale-up of each program/programs + plots = { + 'HIV programs': [ + 'HIV Programs Scale-up WITHOUT HSS PACKAGE', + 'HIV Programs Scale-up WITH HSS PACKAGE', + ], + 'TB programs': [ + 'TB Programs Scale-up WITHOUT HSS PACKAGE', + 'TB Programs Scale-up WITH HSS PACKAGE', + ], + 'Malaria programs': [ + 'Malaria Programs Scale-up WITHOUT HSS PACKAGE', + 'Malaria Programs Scale-up WITH HSS PACKAGE', + ], + 'All programs': [ + 'FULL HSS PACKAGE', + 'HIV/Tb/Malaria Programs Scale-up WITHOUT HSS PACKAGE', + 'HIV/Tb/Malaria Programs Scale-up WITH HSS PACKAGE', + ] + } + + for plot_name, scenario_names in plots.items(): + name_of_plot = f'{plot_name}' + fig, ax = plt.subplots() + total_num_dalys_by_label_results_averted_vs_baseline[scenario_names].T.plot.bar( + stacked=True, + ax=ax, + rot=0, + alpha=0.75 + ) + ax.set_ylim([0, 10e7]) + ax.set_title(name_of_plot) + ax.set_ylabel(f'DALYs Averted vs Baseline, {target_period()}\n(Millions)') + wrapped_labs = ["\n".join(textwrap.wrap(_lab.get_text(), 20)) for _lab in ax.get_xticklabels()] + ax.set_xticklabels(wrapped_labs) + fig.tight_layout() + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + # todo: Neaten graphs + # todo: other metrics of health + # todo: other graphs, broken down by age/sex (this can also be cribbed from overview paper stuff) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("results_folder", type=Path) # outputs/horizontal_and_vertical_programs-2024-05-16 + args = parser.parse_args() + + apply( + results_folder=args.results_folder, + output_folder=args.results_folder, + resourcefilepath=Path('./resources') + ) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/mini_analysis_for_testing/mini_version_scenario.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/mini_analysis_for_testing/mini_version_scenario.py new file mode 100644 index 0000000000..a991c019ba --- /dev/null +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/mini_analysis_for_testing/mini_version_scenario.py @@ -0,0 +1,109 @@ +"""This Scenario file is intended to help with debugging the scale-up of HIV. Tb and Malaria services, per issue #1413. + +Changes to the main analysis: + +* We're running this in MODE 1 and we're only looking. +* We're capturing the logged output from HIV, Tb and malaria +* We're limiting it to few scenarios: baseline + the scale-up of all HTM programs (no HealthSystem scale-up) + +""" + +from pathlib import Path +from typing import Dict + +from scripts.comparison_of_horizontal_and_vertical_programs.scenario_definitions import ( + ScenarioDefinitions, +) +from tlo import Date, logging +from tlo.analysis.utils import get_parameters_for_status_quo, mix_scenarios +from tlo.methods.fullmodel import fullmodel +from tlo.methods.scenario_switcher import ImprovedHealthSystemAndCareSeekingScenarioSwitcher +from tlo.scenario import BaseScenario + + +class MiniRunHTMWithAndWithoutHSS(BaseScenario): + def __init__(self): + super().__init__() + self.seed = 0 + self.start_date = Date(2010, 1, 1) + self.end_date = Date(2031, 1, 1) + self.pop_size = 100_000 + self._scenarios = self._get_scenarios() + self.number_of_draws = len(self._scenarios) + self.runs_per_draw = 1 + + def log_configuration(self): + return { + 'filename': 'mini_htm_with_and_without_hss', + 'directory': Path('./outputs'), + 'custom_levels': { + '*': logging.WARNING, + 'tlo.methods.demography': logging.INFO, + 'tlo.methods.demography.detail': logging.WARNING, + 'tlo.methods.healthburden': logging.INFO, + 'tlo.methods.healthsystem': logging.WARNING, + 'tlo.methods.healthsystem.summary': logging.INFO, + 'tlo.methods.hiv': logging.INFO, + 'tlo.methods.tb': logging.INFO, + 'tlo.methods.malaria': logging.INFO, + } + } + + def modules(self): + return ( + fullmodel(resourcefilepath=self.resources) + + [ImprovedHealthSystemAndCareSeekingScenarioSwitcher(resourcefilepath=self.resources)] + ) + + def draw_parameters(self, draw_number, rng): + if draw_number < len(self._scenarios): + return list(self._scenarios.values())[draw_number] + + def _get_scenarios(self) -> Dict[str, Dict]: + """Return the Dict with values for the parameters that are changed, keyed by a name for the scenario.""" + # Load helper class containing the definitions of the elements of all the scenarios + scenario_definitions = ScenarioDefinitions() + + return { + "Baseline": + self._baseline(), + + # - - - HIV & TB & MALARIA SCALE-UP WITHOUT HSS PACKAGE- - - + "HIV/Tb/Malaria Programs Scale-up WITHOUT HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._hiv_scaleup(), + scenario_definitions._tb_scaleup(), + scenario_definitions._malaria_scaleup(), + ), + } + + def _baseline(self): + self.YEAR_OF_CHANGE_FOR_HSS = 2019 + + return mix_scenarios( + get_parameters_for_status_quo(), + { + "HealthSystem": { + "mode_appt_constraints": 1, # <-- Mode 1 prior to change to preserve calibration + "mode_appt_constraints_postSwitch": 1, # <-- ***** NO CHANGE --- STAYING IN MODE 1 + "scale_to_effective_capabilities": False, # <-- irrelevant, as not changing mode + "year_mode_switch": 2100, # <-- irrelevant as not changing modes + + # Baseline scenario is with absence of HCW + 'year_HR_scaling_by_level_and_officer_type': self.YEAR_OF_CHANGE_FOR_HSS, + 'HR_scaling_by_level_and_officer_type_mode': 'with_absence', + + # Normalize the behaviour of Mode 2 (irrelevant as in Mode 1) + "policy_name": "Naive", + "tclose_overwrite": 1, + "tclose_days_offset_overwrite": 7, + } + }, + ) + + +if __name__ == '__main__': + from tlo.cli import scenario_run + + scenario_run([__file__]) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_definitions.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_definitions.py new file mode 100644 index 0000000000..6465b57f49 --- /dev/null +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_definitions.py @@ -0,0 +1,90 @@ +"""The file contains all the definitions of scenarios used the Horizontal and Vertical Program Impact Analyses""" +from typing import Dict + +from tlo.analysis.utils import get_parameters_for_status_quo, mix_scenarios + + +class ScenarioDefinitions: + + @property + def YEAR_OF_CHANGE_FOR_HSS(self) -> int: + """Year in which Health Systems Strengthening changes are made.""" + return 2019 # <-- baseline year of Human Resources for Health is 2018, and this is consistent with calibration + # during 2015-2019 period. + + + @property + def YEAR_OF_CHANGE_FOR_HTM(self) -> int: + """Year in which HIV, TB, Malaria scale-up changes are made.""" + return 2019 # todo <-- what is the natural year of scale-up? Should this be the same as the when the HSS + # changes happen? + + def _baseline(self) -> Dict: + """Return the Dict with values for the parameter changes that define the baseline scenario. """ + return mix_scenarios( + get_parameters_for_status_quo(), + { + "HealthSystem": { + "mode_appt_constraints": 1, # <-- Mode 1 prior to change to preserve calibration + "mode_appt_constraints_postSwitch": 2, # <-- Mode 2 post-change to show effects of HRH + "scale_to_effective_capabilities": True, + # <-- Transition into Mode2 with the effective capabilities in HRH 'revealed' in Mode 1 + "year_mode_switch": self.YEAR_OF_CHANGE_FOR_HSS, + + # Baseline scenario is with absence of HCW + 'year_HR_scaling_by_level_and_officer_type': self.YEAR_OF_CHANGE_FOR_HSS, + 'HR_scaling_by_level_and_officer_type_mode': 'with_absence', + # todo <-- Do we want the first part of the run be with_abscence too...? (Although that will mean + # that there is actually greater capacity if we do the rescaling) + + # Normalize the behaviour of Mode 2 + "policy_name": "Naive", + "tclose_overwrite": 1, + "tclose_days_offset_overwrite": 7, + } + }, + ) + + def _hss_package(self) -> Dict: + """The parameters for the Health System Strengthening Package""" + return { + 'ImprovedHealthSystemAndCareSeekingScenarioSwitcher': { + 'max_healthsystem_function': [False, True], # <-- switch from False to True mid-way + 'max_healthcare_seeking': [False, True], # <-- switch from False to True mid-way + 'year_of_switch': self.YEAR_OF_CHANGE_FOR_HSS + }, + 'HealthSystem': { + 'year_cons_availability_switch': self.YEAR_OF_CHANGE_FOR_HSS, + 'cons_availability_postSwitch': 'all', + 'yearly_HR_scaling_mode': 'GDP_growth_fHE_case5', + 'year_HR_scaling_by_level_and_officer_type': self.YEAR_OF_CHANGE_FOR_HSS, + 'HR_scaling_by_level_and_officer_type_mode': 'no_absence_&_x2_fac0+1', + } + } + + def _hiv_scaleup(self) -> Dict: + """The parameters for the scale-up of the HIV program""" + return { + "Hiv": { + 'do_scaleup': True, + 'scaleup_start_year': self.YEAR_OF_CHANGE_FOR_HTM, + } + } + + def _tb_scaleup(self) -> Dict: + """The parameters for the scale-up of the TB program""" + return { + "Tb": { + 'do_scaleup': True, + 'scaleup_start_year': self.YEAR_OF_CHANGE_FOR_HTM, + } + } + + def _malaria_scaleup(self) -> Dict: + """The parameters for the scale-up of the Malaria program""" + return { + 'Malaria': { + 'do_scaleup': True, + 'scaleup_start_year': self.YEAR_OF_CHANGE_FOR_HTM, + } + } diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_hss_elements.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_hss_elements.py new file mode 100644 index 0000000000..e1aebec8f6 --- /dev/null +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_hss_elements.py @@ -0,0 +1,242 @@ +"""This Scenario file run the model under different assumptions for the HealthSystem and Vertical Program Scale-up + +Run on the batch system using: +``` +tlo batch-submit + src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_hss_elements.py +``` + +""" + +from pathlib import Path +from typing import Dict + +# from scripts.comparison_of_horizontal_and_vertical_programs.scenario_definitions import ( +# ScenarioDefinitions, +# ) +from tlo import Date, logging +from tlo.analysis.utils import get_parameters_for_status_quo, mix_scenarios +from tlo.methods.fullmodel import fullmodel +from tlo.methods.scenario_switcher import ImprovedHealthSystemAndCareSeekingScenarioSwitcher +from tlo.scenario import BaseScenario + + +class HSSElements(BaseScenario): + def __init__(self): + super().__init__() + self.seed = 0 + self.start_date = Date(2010, 1, 1) + self.end_date = Date(2031, 1, 1) + self.pop_size = 100_000 + self._scenarios = self._get_scenarios() + self.number_of_draws = len(self._scenarios) + self.runs_per_draw = 3 # <--- todo: N.B. Very small number of repeated run, to be efficient for now + + def log_configuration(self): + return { + 'filename': 'hss_elements', + 'directory': Path('./outputs'), + 'custom_levels': { + '*': logging.WARNING, + 'tlo.methods.demography': logging.INFO, + 'tlo.methods.demography.detail': logging.WARNING, + 'tlo.methods.healthburden': logging.INFO, + 'tlo.methods.healthsystem': logging.WARNING, + 'tlo.methods.healthsystem.summary': logging.INFO, + } + } + + def modules(self): + return ( + fullmodel(resourcefilepath=self.resources) + + [ImprovedHealthSystemAndCareSeekingScenarioSwitcher(resourcefilepath=self.resources)] + ) + + def draw_parameters(self, draw_number, rng): + if draw_number < len(self._scenarios): + return list(self._scenarios.values())[draw_number] + + def _get_scenarios(self) -> Dict[str, Dict]: + """Return the Dict with values for the parameters that are changed, keyed by a name for the scenario.""" + # todo - decide on final definition of scenarios and the scenario package + # todo - refactorize to use the ScenariosDefinitions helperclass, which will make sure that this script and + # 'scenario_vertical_programs)_with_and_without_hss.py' are synchronised (e.g. baseline and HSS pkg scenarios) + + self.YEAR_OF_CHANGE = 2019 + # <-- baseline year of Human Resources for Health is 2018, and this is consistent with calibration during + # 2015-2019 period. + + return { + "Baseline": self._baseline(), + + # *************************** + # HEALTH SYSTEM STRENGTHENING + # *************************** + + # - - - Human Resource for Health - - - + + "Reduced Absence": + mix_scenarios( + self._baseline(), + { + 'HealthSystem': { + 'year_HR_scaling_by_level_and_officer_type': self.YEAR_OF_CHANGE, + 'HR_scaling_by_level_and_officer_type_mode': 'no_absence', + } + } + ), + + "Reduced Absence + Double Capacity at Primary Care": + mix_scenarios( + self._baseline(), + { + 'HealthSystem': { + 'year_HR_scaling_by_level_and_officer_type': self.YEAR_OF_CHANGE, + 'HR_scaling_by_level_and_officer_type_mode': 'no_absence_&_x2_fac0+1', + } + } + ), + + "HRH Keeps Pace with Population Growth": + mix_scenarios( + self._baseline(), + { + 'HealthSystem': { + 'yearly_HR_scaling_mode': 'scaling_by_population_growth', + # This is in-line with population growth _after 2018_ (baseline year for HRH) + } + } + ), + + "HRH Increases at GDP Growth": + mix_scenarios( + self._baseline(), + { + 'HealthSystem': { + 'yearly_HR_scaling_mode': 'GDP_growth', + # This is GDP growth after 2018 (baseline year for HRH) + } + } + ), + + "HRH Increases above GDP Growth": + mix_scenarios( + self._baseline(), + { + 'HealthSystem': { + 'yearly_HR_scaling_mode': 'GDP_growth_fHE_case5', + # This is above-GDP growth after 2018 (baseline year for HRH) + } + } + ), + + + # - - - Quality of Care - - - + "Perfect Clinical Practice": + mix_scenarios( + self._baseline(), + { + 'ImprovedHealthSystemAndCareSeekingScenarioSwitcher': { + 'max_healthsystem_function': [False, True], # <-- switch from False to True mid-way + 'year_of_switch': self.YEAR_OF_CHANGE, + } + }, + ), + + "Perfect Healthcare Seeking": + mix_scenarios( + get_parameters_for_status_quo(), + { + 'ImprovedHealthSystemAndCareSeekingScenarioSwitcher': { + 'max_healthcare_seeking': [False, True], + 'year_of_switch': self.YEAR_OF_CHANGE, + } + }, + ), + + # - - - Supply Chains - - - + "Perfect Availability of Vital Items": + mix_scenarios( + self._baseline(), + { + 'HealthSystem': { + 'year_cons_availability_switch': self.YEAR_OF_CHANGE, + 'cons_availability_postSwitch': 'all_vital_available', + } + } + ), + + "Perfect Availability of Medicines": + mix_scenarios( + self._baseline(), + { + 'HealthSystem': { + 'year_cons_availability_switch': self.YEAR_OF_CHANGE, + 'cons_availability_postSwitch': 'all_medicines_available', + } + } + ), + + "Perfect Availability of All Consumables": + mix_scenarios( + self._baseline(), + { + 'HealthSystem': { + 'year_cons_availability_switch': self.YEAR_OF_CHANGE, + 'cons_availability_postSwitch': 'all', + } + } + ), + + # - - - FULL PACKAGE OF HEALTH SYSTEM STRENGTHENING - - - + "FULL PACKAGE": + mix_scenarios( + self._baseline(), + { + 'ImprovedHealthSystemAndCareSeekingScenarioSwitcher': { + 'max_healthsystem_function': [False, True], # <-- switch from False to True mid-way + 'max_healthcare_seeking': [False, True], # <-- switch from False to True mid-way + 'year_of_switch': self.YEAR_OF_CHANGE + }, + 'HealthSystem': { + 'year_cons_availability_switch': self.YEAR_OF_CHANGE, + 'cons_availability_postSwitch': 'all', + 'yearly_HR_scaling_mode': 'GDP_growth_fHE_case5', + 'year_HR_scaling_by_level_and_officer_type': self.YEAR_OF_CHANGE, + 'HR_scaling_by_level_and_officer_type_mode': 'no_absence_&_x2_fac0+1', + } + }, + ), + + } + + def _baseline(self) -> Dict: + """Return the Dict with values for the parameter changes that define the baseline scenario. """ + return mix_scenarios( + get_parameters_for_status_quo(), + { + "HealthSystem": { + "mode_appt_constraints": 1, # <-- Mode 1 prior to change to preserve calibration + "mode_appt_constraints_postSwitch": 2, # <-- Mode 2 post-change to show effects of HRH + "scale_to_effective_capabilities": True, # <-- Transition into Mode2 with the effective + # capabilities in HRH 'revealed' in Mode 1 + "year_mode_switch": self.YEAR_OF_CHANGE, + + # Baseline scenario is with absence of HCW + 'year_HR_scaling_by_level_and_officer_type': self.YEAR_OF_CHANGE, + 'HR_scaling_by_level_and_officer_type_mode': 'with_absence', + # todo <-- Do we want the first part of the run be with_abscence too...? (Although that will mean + # that there is actually greater capacity if we do the rescaling) + + # Normalize the behaviour of Mode 2 + "policy_name": "Naive", + "tclose_overwrite": 1, + "tclose_days_offset_overwrite": 7, + } + }, + ) + +if __name__ == '__main__': + from tlo.cli import scenario_run + + scenario_run([__file__]) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_vertical_programs_with_and_without_hss.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_vertical_programs_with_and_without_hss.py new file mode 100644 index 0000000000..ce7f664fe9 --- /dev/null +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_vertical_programs_with_and_without_hss.py @@ -0,0 +1,144 @@ +"""This Scenario file run the model under different assumptions for the HealthSystem and Vertical Program Scale-up + +Run on the batch system using: +``` +tlo batch-submit + src/scripts/comparison_of_horizontal_and_vertical_programs/scenario_vertical_programs_with_and_without_hss.py +``` + +""" + +from pathlib import Path +from typing import Dict + +from scripts.comparison_of_horizontal_and_vertical_programs.scenario_definitions import ( + ScenarioDefinitions, +) +from tlo import Date, logging +from tlo.analysis.utils import mix_scenarios +from tlo.methods.fullmodel import fullmodel +from tlo.methods.scenario_switcher import ImprovedHealthSystemAndCareSeekingScenarioSwitcher +from tlo.scenario import BaseScenario + + +class HTMWithAndWithoutHSS(BaseScenario): + def __init__(self): + super().__init__() + self.seed = 0 + self.start_date = Date(2010, 1, 1) + self.end_date = Date(2031, 1, 1) + self.pop_size = 100_000 + self._scenarios = self._get_scenarios() + self.number_of_draws = len(self._scenarios) + self.runs_per_draw = 3 # <--- todo: N.B. Very small number of repeated run, to be efficient for now + + def log_configuration(self): + return { + 'filename': 'htm_with_and_without_hss', + 'directory': Path('./outputs'), + 'custom_levels': { + '*': logging.WARNING, + 'tlo.methods.demography': logging.INFO, + 'tlo.methods.demography.detail': logging.WARNING, + 'tlo.methods.healthburden': logging.INFO, + 'tlo.methods.healthsystem': logging.WARNING, + 'tlo.methods.healthsystem.summary': logging.INFO, + } + } + + def modules(self): + return ( + fullmodel(resourcefilepath=self.resources) + + [ImprovedHealthSystemAndCareSeekingScenarioSwitcher(resourcefilepath=self.resources)] + ) + + def draw_parameters(self, draw_number, rng): + if draw_number < len(self._scenarios): + return list(self._scenarios.values())[draw_number] + + def _get_scenarios(self) -> Dict[str, Dict]: + """Return the Dict with values for the parameters that are changed, keyed by a name for the scenario.""" + # Load helper class containing the definitions of the elements of all the scenarios + scenario_definitions = ScenarioDefinitions() + + return { + "Baseline": + scenario_definitions._baseline(), + + # - - - FULL PACKAGE OF HEALTH SYSTEM STRENGTHENING - - - + "FULL HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._hss_package(), + ), + + # ************************************************** + # VERTICAL PROGRAMS WITH AND WITHOUT THE HSS PACKAGE + # ************************************************** + + # - - - HIV SCALE-UP WITHOUT HSS PACKAGE- - - + "HIV Programs Scale-up WITHOUT HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._hiv_scaleup(), + ), + # - - - HIV SCALE-UP *WITH* HSS PACKAGE- - - + "HIV Programs Scale-up WITH HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._hiv_scaleup(), + scenario_definitions._hss_package(), + ), + + # - - - TB SCALE-UP WITHOUT HSS PACKAGE- - - + "TB Programs Scale-up WITHOUT HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._tb_scaleup(), + ), + # - - - TB SCALE-UP *WITH* HSS PACKAGE- - - + "TB Programs Scale-up WITH HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._tb_scaleup(), + scenario_definitions._hss_package(), + ), + + # - - - MALARIA SCALE-UP WITHOUT HSS PACKAGE- - - + "Malaria Programs Scale-up WITHOUT HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._malaria_scaleup(), + ), + # - - - MALARIA SCALE-UP *WITH* HSS PACKAGE- - - + "Malaria Programs Scale-up WITH HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._malaria_scaleup(), + scenario_definitions._hss_package(), + ), + + # - - - HIV & TB & MALARIA SCALE-UP WITHOUT HSS PACKAGE- - - + "HIV/Tb/Malaria Programs Scale-up WITHOUT HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._hiv_scaleup(), + scenario_definitions._tb_scaleup(), + scenario_definitions._malaria_scaleup(), + ), + # - - - HIV & TB & MALARIA SCALE-UP *WITH* HSS PACKAGE- - - + "HIV/Tb/Malaria Programs Scale-up WITH HSS PACKAGE": + mix_scenarios( + scenario_definitions._baseline(), + scenario_definitions._hiv_scaleup(), + scenario_definitions._tb_scaleup(), + scenario_definitions._malaria_scaleup(), + scenario_definitions._hss_package(), + ), + } + + +if __name__ == '__main__': + from tlo.cli import scenario_run + + scenario_run([__file__])