diff --git a/resources/healthsystem/ResourceFile_HealthSystem_parameters.csv b/resources/healthsystem/ResourceFile_HealthSystem_parameters.csv index c6bd6414e7..9822923092 100644 --- a/resources/healthsystem/ResourceFile_HealthSystem_parameters.csv +++ b/resources/healthsystem/ResourceFile_HealthSystem_parameters.csv @@ -14,6 +14,8 @@ equip_availability_postSwitch,default year_equip_availability_switch,2100 tclose_overwrite,0 tclose_days_offset_overwrite,7 +HR_scaling_by_year_and_officer_type_mode,no_historical_growth +HR_expansion_absorption_rate,1.0 HR_scaling_by_level_and_officer_type_mode,default year_HR_scaling_by_level_and_officer_type,2100 HR_scaling_by_district_mode,default diff --git a/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/custom.csv b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/custom.csv new file mode 100644 index 0000000000..5edeb2def5 --- /dev/null +++ b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/custom.csv @@ -0,0 +1,10 @@ +Officer_Category,2010,2019,2020,2021,2022,2023,2024,2025,2026,2027,2028,2029,2030 +Clinical,1,1,1,1,1,1,1,1,1,1,1,1,1 +DCSA,1,1,1,1,1,1,1,1,1,1,1,1,1 +Dental,1,1,1,1,1,1,1,1,1,1,1,1,1 +Laboratory,1,1,1,1,1,1,1,1,1,1,1,1,1 +Mental,1,1,1,1,1,1,1,1,1,1,1,1,1 +Nursing_and_Midwifery,1,1,1,1,1,1,1,1,1,1,1,1,1 +Nutrition,1,1,1,1,1,1,1,1,1,1,1,1,1 +Pharmacy,1,1,1,1,1,1,1,1,1,1,1,1,1 +Radiography,1,1,1,1,1,1,1,1,1,1,1,1,1 diff --git a/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/historical_cadre_mix.csv b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/historical_cadre_mix.csv new file mode 100644 index 0000000000..10fbad62d1 --- /dev/null +++ b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/historical_cadre_mix.csv @@ -0,0 +1,10 @@ +Officer_Category,2010,2019,2020,2021,2022,2023,2024,2025,2026,2027,2028,2029,2030 +Clinical,1,1,1.068268523,1.182359384,1.168541512,1.156685635,1.146388516,1,1,1,1,1,1 +DCSA,1,1,1,1,1,1,1,1,1,1,1,1,1 +Dental,1,1,1.019505292,1.052102681,1.048154718,1.044767324,1.04182529,1,1,1,1,1,1 +Laboratory,1,1,1.019505292,1.052102681,1.048154718,1.044767324,1.04182529,1,1,1,1,1,1 +Mental,1,1,1.019505292,1.052102681,1.048154718,1.044767324,1.04182529,1,1,1,1,1,1 +Nursing_and_Midwifery,1,1,1.039010585,1.104205363,1.096309435,1.089534649,1.08365058,1,1,1,1,1,1 +Nutrition,1,1,1.019505292,1.052102681,1.048154718,1.044767324,1.04182529,1,1,1,1,1,1 +Pharmacy,1,1,1.136537047,1.364718769,1.337083024,1.31337127,1.292777032,1,1,1,1,1,1 +Radiography,1,1,1.019505292,1.052102681,1.048154718,1.044767324,1.04182529,1,1,1,1,1,1 diff --git a/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/historical_uniform.csv b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/historical_uniform.csv new file mode 100644 index 0000000000..f312cdfc80 --- /dev/null +++ b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/historical_uniform.csv @@ -0,0 +1,10 @@ +Officer_Category,2010,2019,2020,2021,2022,2023,2024,2025,2026,2027,2028,2029,2030 +Clinical,1,1,1.027745477,1.076495291,1.076495291,1.076495291,1.076495291,1,1,1,1,1,1 +DCSA,1,1,1.027745477,1.076495291,1.076495291,1.076495291,1.076495291,1,1,1,1,1,1 +Dental,1,1,1.027745477,1.076495291,1.076495291,1.076495291,1.076495291,1,1,1,1,1,1 +Laboratory,1,1,1.027745477,1.076495291,1.076495291,1.076495291,1.076495291,1,1,1,1,1,1 +Mental,1,1,1.027745477,1.076495291,1.076495291,1.076495291,1.076495291,1,1,1,1,1,1 +Nursing_and_Midwifery,1,1,1.027745477,1.076495291,1.076495291,1.076495291,1.076495291,1,1,1,1,1,1 +Nutrition,1,1,1.027745477,1.076495291,1.076495291,1.076495291,1.076495291,1,1,1,1,1,1 +Pharmacy,1,1,1.027745477,1.076495291,1.076495291,1.076495291,1.076495291,1,1,1,1,1,1 +Radiography,1,1,1.027745477,1.076495291,1.076495291,1.076495291,1.076495291,1,1,1,1,1,1 diff --git a/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/no_historical_growth.csv b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/no_historical_growth.csv new file mode 100644 index 0000000000..5edeb2def5 --- /dev/null +++ b/resources/healthsystem/human_resources/scaling_capabilities/ResourceFile_HR_scaling_by_year_and_officer_type/no_historical_growth.csv @@ -0,0 +1,10 @@ +Officer_Category,2010,2019,2020,2021,2022,2023,2024,2025,2026,2027,2028,2029,2030 +Clinical,1,1,1,1,1,1,1,1,1,1,1,1,1 +DCSA,1,1,1,1,1,1,1,1,1,1,1,1,1 +Dental,1,1,1,1,1,1,1,1,1,1,1,1,1 +Laboratory,1,1,1,1,1,1,1,1,1,1,1,1,1 +Mental,1,1,1,1,1,1,1,1,1,1,1,1,1 +Nursing_and_Midwifery,1,1,1,1,1,1,1,1,1,1,1,1,1 +Nutrition,1,1,1,1,1,1,1,1,1,1,1,1,1 +Pharmacy,1,1,1,1,1,1,1,1,1,1,1,1,1 +Radiography,1,1,1,1,1,1,1,1,1,1,1,1,1 diff --git a/src/scripts/impact_of_historical_changes_in_hr/analysis_historical_changes_in_hr_extended.py b/src/scripts/impact_of_historical_changes_in_hr/analysis_historical_changes_in_hr_extended.py new file mode 100644 index 0000000000..1a333af5eb --- /dev/null +++ b/src/scripts/impact_of_historical_changes_in_hr/analysis_historical_changes_in_hr_extended.py @@ -0,0 +1,503 @@ +""" +Extract results of HRH staff counts, DALYs and Deaths across multiple historical HRH growth scenarios. +""" + +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 scripts.impact_of_historical_changes_in_hr.scenario_historical_changes_in_hr_extended import ( + HistoricalChangesInHRH, +) +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, the_target_period: Tuple[Date, Date] = None): + + TARGET_PERIOD = the_target_period + hrh_check_period = (Date(2018, 1, 1), Date(2026, 1, 1)) + + 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.""" + e = HistoricalChangesInHRH() + return tuple(e._scenarios.keys()) + + def get_num_deaths_by_year_cause(_df): + """Return total number of Deaths (total within the TARGET_PERIOD)""" + _df = _df.loc[pd.to_datetime(_df.date).between(*TARGET_PERIOD)] + _df['year'] = _df['date'].dt.year + _df = _df.groupby(['year', 'cause']) \ + .agg({'person_id': 'count'}) \ + .rename(columns={'person_id': 'num_deaths'})['num_deaths'] + + return _df + + 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 get_total_num_dalys_by_label_all_causes(_df): + """Return the total number of DALYS in the TARGET_PERIOD cause label.""" + _df = _df \ + .loc[_df['year'].between(*[d.year for d in TARGET_PERIOD])] \ + .drop(columns=['date', 'age_range', 'sex']) \ + .groupby('year') \ + .sum() \ + .reset_index() \ + .melt(id_vars='year', var_name='cause', value_name='dalys') \ + .set_index(['year', 'cause'])['dalys'] + + return _df + + def get_staff_counts(_df): + _df['year'] = _df['date'].dt.year + _df = _df.loc[pd.to_datetime(_df['date']).between(*hrh_check_period), ['year', 'GenericClinic'] + ].set_index('year').rename(columns={'GenericClinic': 'facility_officer'}) + _df_staff = _df['facility_officer'].apply(pd.Series).stack().reset_index() + _df_staff.columns = ['year', 'facility_officer', 'staff_count'] + _df_staff[['facility_id', 'officer_type']] = _df_staff['facility_officer'].str.extract( + r'FacilityID_(\d+)_Officer_(.*)' + ) + _df_staff['facility_id'] = _df_staff['facility_id'].astype(int) + + main_cadres = ['Clinical', 'Nursing_and_Midwifery', 'Pharmacy', 'DCSA'] + _df_staff.loc[ + ~_df_staff["officer_type"].isin(main_cadres), + "officer_type" + ] = "Other" + + _df_staff = _df_staff.groupby(['year', 'officer_type'])['staff_count'].sum() + return _df_staff + + 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_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 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 get_num_treatments_by_year_treatment(_df): + """Return the number of treatments by short treatment id and year (within the TARGET_PERIOD)""" + _df['year'] = _df['date'].dt.year + _df = _df.loc[pd.to_datetime(_df.date).between(*TARGET_PERIOD), ['year', 'TREATMENT_ID']].set_index('year') + _df = _df['TREATMENT_ID'].apply(pd.Series) + _df.columns = _df.columns.map(lambda x: x.split('_')[0] + "*") + _df = _df.T.groupby(level=0).sum().T + _df = _df.stack() + _df.index = _df.index.set_names(["year", "treatment_type"]) + _df.name = "count" + return _df + + 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(), + zorder=100, + ) + 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, _ = 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() + + # Check HRH staff counts + hcw_count = (extract_results( + results_folder, + module="tlo.methods.healthsystem.summary", + key="number_of_hcw_staff", + custom_generate_series=get_staff_counts, + do_scaling=False + )).pipe(set_param_names_as_column_index_level_0) + hcw_count = hcw_count.loc[:, hcw_count.columns.get_level_values("run") == 0] + hcw_count.columns = hcw_count.columns.get_level_values(0) + + hcw_count = ( + hcw_count.stack() + .reset_index(name='value') + .rename(columns={'draw': 'scenario'}) + ) + + hcw_count = hcw_count.sort_values(['officer_type', 'scenario', 'year']) + + hcw_count['scale_factor'] = ( + hcw_count['value'] / + hcw_count.groupby(['officer_type', 'scenario'])['value'].shift(1) + ) + + from matplotlib.lines import Line2D + # marker for each officer type + markers = { + 'Clinical': 'o', + 'Nursing_and_Midwifery': '*', + 'Pharmacy': '^', + 'DCSA': 'd', + 'Other': 'X', + } + marker_sizes = { + 'Clinical': 6, + 'Nursing_and_Midwifery': 8, + 'Pharmacy': 6, + 'DCSA': 6, + 'Other': 6, + } + + # color for each scenario + cmap = plt.cm.get_cmap('tab20', len(param_names)) + colors = { + scenario: cmap(i) + for i, scenario in enumerate(param_names) + } + + name_of_plot = "Number of healthcare workers" + fig, ax = plt.subplots(figsize=(12, 5)) + + for officer in hcw_count['officer_type'].unique(): + for scenario in hcw_count['scenario'].unique(): + subset = hcw_count[ + (hcw_count['officer_type'] == officer) & + (hcw_count['scenario'] == scenario) + ] + + ax.plot( + subset['year'], + subset['value'], + linestyle="--", + marker=markers[officer], + markersize=marker_sizes[officer], + color=colors[scenario], + label=f'{officer} - {scenario}', + ) + + ax.set_xlabel('Year') + ax.set_ylabel('Number of staff') + ax.set_xticks(sorted(hcw_count['year'].unique())) + + officer_handles = [ + Line2D( + [0], + [0], + marker=markers[officer], + color='black', + linestyle='None', + markersize=8, + label=officer + ) + for officer in markers.keys() + ] + + legend_officer = ax.legend( + handles=officer_handles, + title='Officer type', + loc='center left', + bbox_to_anchor=(1.02, 0.65), + fontsize=8, + title_fontsize=9 + ) + ax.add_artist(legend_officer) + + scenario_handles = [ + Line2D( + [0], + [0], + color=colors[scenario], + linestyle='--', + linewidth=2, + label=scenario + ) + for scenario in colors.keys() + ] + + ax.legend( + handles=scenario_handles, + title='Scenario', + loc='center left', + bbox_to_anchor=(1.02, 0.25), + fontsize=8, + title_fontsize=9 + ) + + plt.tight_layout() + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + plt.show() + plt.close(fig) + + # Check total DALYs + # %% Define parameter names + counterfactual_scenario = 'Main Counterfactual/No growth + Lower bound settings' + actual_scenario = 'Main Actual' + + # 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, xticklabels_horizontal_and_wrapped=True, + put_labels_in_legend=True) + ax.set_title(name_of_plot) + ax.set_ylabel('(Millions)') + fig.tight_layout() + ax.axhline(num_deaths_summarized.loc[counterfactual_scenario, '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'DALYs, {target_period()}' + fig, ax = do_bar_plot_with_ci(num_dalys_summarized / 1e6, xticklabels_horizontal_and_wrapped=True, + put_labels_in_legend=True) + ax.set_title(name_of_plot) + ax.set_ylabel('(Millions)') + ax.axhline(num_dalys_summarized.loc[counterfactual_scenario, '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 Actual + num_deaths_averted = summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison_series( + num_deaths.loc[0], + comparison=actual_scenario) + ).T + ).iloc[0].unstack().reindex(param_names).drop([actual_scenario]) + + pc_deaths_averted = 100.0 * summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison_series( + num_deaths.loc[0], + comparison=actual_scenario, + scaled=True) + ).T + ).iloc[0].unstack().reindex(param_names).drop([actual_scenario]) + + num_dalys_averted = summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison_series( + num_dalys.loc[0], + comparison=actual_scenario) + ).T + ).iloc[0].unstack().reindex(param_names).drop([actual_scenario]) + + pc_dalys_averted = 100.0 * summarize( + -1.0 * + pd.DataFrame( + find_difference_relative_to_comparison_series( + num_dalys.loc[0], + comparison=actual_scenario, + scaled=True) + ).T + ).iloc[0].unstack().reindex(param_names).drop([actual_scenario]) + + # DEATHS + name_of_plot = f'Deaths Averted vs Historical growth (uniform), {target_period()}' + fig, ax = do_bar_plot_with_ci( + pc_deaths_averted, # num_deaths_averted + annotations=None, + put_labels_in_legend=True, + xticklabels_horizontal_and_wrapped=True, + ) + # annotation = ( + # f"{int(round(num_deaths_averted.loc[actual_scenario, 'mean'], -3))} ({int(round(num_deaths_averted.loc[actual_scenario, 'lower'], -3))} - {int(round(num_deaths_averted.loc[actual_scenario, 'upper'], -3))})\n" + # f"{round(pc_deaths_averted.loc[actual_scenario, 'mean'])} ({round(pc_deaths_averted.loc[actual_scenario, 'lower'], 1)} - {round(pc_deaths_averted.loc[actual_scenario, 'upper'], 1)})% of that in Counterfactual" + # ) + ax.set_title(f"{name_of_plot}") + ax.set_ylabel('Deaths Averted vs Historical growth (uniform)') + # fig.set_figwidth(5) + 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 Historical growth (uniform), {target_period()}' + fig, ax = do_bar_plot_with_ci( + pc_dalys_averted, # (num_dalys_averted / 1e6), + annotations=None, + put_labels_in_legend=True, + xticklabels_horizontal_and_wrapped=True, + ) + # annotation = ( + # f"{int(round(num_dalys_averted.loc[actual_scenario, 'mean'], -4))} ({int(round(num_dalys_averted.loc[actual_scenario, 'lower'], -4))} - {int(round(num_dalys_averted.loc[actual_scenario, 'upper'], -4))})\n" + # f"{round(pc_dalys_averted.loc[actual_scenario, 'mean'])} ({round(pc_dalys_averted.loc[actual_scenario, 'lower'], 1)} - {round(pc_dalys_averted.loc[actual_scenario, 'upper'], 1)})% of that in Counterfactual" + # ) + ax.set_title(f"{name_of_plot}") + ax.set_ylabel('DALYS Averted vs Historical growth (uniform)') + # fig.set_figwidth(5) + fig.tight_layout() + fig.savefig(make_graph_file_name(name_of_plot.replace(' ', '_').replace(',', ''))) + fig.show() + plt.close(fig) + + # Prepare Absolute Number of Deaths and DALYs for Izzy + num_deaths_by_year_cause = extract_results( + results_folder, + module='tlo.methods.demography', + key='death', + custom_generate_series=get_num_deaths_by_year_cause, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0).stack(['draw', 'run']).reset_index(name='num_deaths') + + num_dalys_by_year_cause = extract_results( + results_folder, + module="tlo.methods.healthburden", + key="dalys_stacked_by_age_and_time", + custom_generate_series=get_total_num_dalys_by_label_all_causes, + do_scaling=True, + ).pipe(set_param_names_as_column_index_level_0).stack(['draw', 'run']).reset_index(name='num_dalys') + + # And absolute Number of treatments upon analysis needs + num_treatments_by_year_treatment = extract_results( + results_folder, + module='tlo.methods.healthsystem.summary', + key='HSI_Event_non_blank_appt_footprint', + custom_generate_series=get_num_treatments_by_year_treatment, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0).stack(['draw', 'run']).reset_index(name='num_treatments') + + num_dalys_by_year_cause.to_csv(output_folder / 'num_dalys_by_year_cause (for Izzy).csv', index=False) + num_deaths_by_year_cause.to_csv(output_folder / 'num_deaths_by_year_cause (for Izzy).csv', index=False) + num_treatments_by_year_treatment.to_csv( + output_folder / 'num_treatments_by_year_treatment (for Izzy).csv', index=False + ) + hcw_count.to_csv(output_folder / 'hcw_count (for Izzy).csv', index=False) + + +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() + + # # Produce results for short-term analysis - 2020 - 2024 (incl.) + # apply( + # results_folder=args.results_folder, + # output_folder=args.results_folder, + # resourcefilepath=Path('./resources'), + # the_target_period=(Date(2020, 1, 1), Date(2024, 12, 31)) + # ) + # Produce results for only later period 2025-2030 (incl.) + apply( + results_folder=args.results_folder, + output_folder=args.results_folder, + resourcefilepath=Path('./resources'), + the_target_period=(Date(2020, 1, 1), Date(2030, 12, 31)) + ) diff --git a/src/scripts/impact_of_historical_changes_in_hr/scenario_historical_changes_in_hr_extended.py b/src/scripts/impact_of_historical_changes_in_hr/scenario_historical_changes_in_hr_extended.py new file mode 100644 index 0000000000..ccfc79aa78 --- /dev/null +++ b/src/scripts/impact_of_historical_changes_in_hr/scenario_historical_changes_in_hr_extended.py @@ -0,0 +1,259 @@ +"""This Scenario file runs the model under different assumptions for the historical changes in Human Resources for Health + +Run on the batch system using: +``` +tlo batch-submit src/scripts/impact_of_historical_changes_in_hr/scenario_historical_changes_in_hr_extended.py +``` + +Or locally using: +``` + +tlo scenario-run src/scripts/impact_of_historical_changes_in_hr/scenario_historical_changes_in_hr_extended.py + +``` + +""" + +from pathlib import Path +from typing import Dict + +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 HistoricalChangesInHRH(BaseScenario): + def __init__(self): + super().__init__() + self.seed = 0 + self.start_date = Date(2010, 1, 1) + self.end_date = Date(2031, 1, 1) # <-- End at the end of year 2030 + self.pop_size = 20_000 # <-- Local small run: 30; previous study run: 20_000; for publication: 100_000 + self._scenarios = self._get_scenarios() + self.number_of_draws = len(self._scenarios) + self.runs_per_draw = 5 # for publication: 10 + + def log_configuration(self): + return { + 'filename': 'historical_changes_in_hr_extended', + '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() + [ImprovedHealthSystemAndCareSeekingScenarioSwitcher()] + ) + + def draw_parameters(self, draw_number, rng): + if draw_number < self.number_of_draws: + return list(self._scenarios.values())[draw_number] + + def draw_name(self, draw_number) -> str: + """Store scenario name. + (This name can be retrieved by the plotting scripts to make the graphs be labelled nicely). + """ + if draw_number < self.number_of_draws: + return list(self._scenarios.keys())[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.""" + + return { + "Main Counterfactual/No growth + Lower bound settings": + self._common_baseline(), + + "Main Actual": + self._hrh_growth_baseline(), + + "Historical growth (cadre-mix)": + mix_scenarios( + self._common_baseline(), + { + "HealthSystem": { + "HR_scaling_by_year_and_officer_type_mode": "historical_cadre_mix", + } + } + ), + + "Historical growth + LCOA policy": + mix_scenarios( + self._hrh_growth_baseline(), + { + "HealthSystem": { + "policy_name": "LCOA_EHP", + } + } + ), + + "Historical growth + Consumables (better)": + mix_scenarios( + self._hrh_growth_baseline(), + { + "HealthSystem": { + "cons_availability_postSwitch": "scenario6", + } + } + ), + + "Historical growth + Consumables (perfect)": + mix_scenarios( + self._hrh_growth_baseline(), + { + "HealthSystem": { + "cons_availability_postSwitch": "all", + } + } + ), + + "Historical growth + Low absorption rate/Historical growth + Lower bound settings": + mix_scenarios( + self._hrh_growth_baseline(), + { + "HealthSystem": { + "HR_expansion_absorption_rate": 0.5, + } + } + ), + + "Historical growth + Max HS performance": + mix_scenarios( + self._hrh_growth_baseline(), + { + 'ImprovedHealthSystemAndCareSeekingScenarioSwitcher': { + 'max_healthsystem_function': [False, True], + } + } + ), + + "No growth + LCOA policy": + mix_scenarios( + self._common_baseline(), + { + "HealthSystem": { + "policy_name": "LCOA_EHP", + } + } + ), + + "No growth + Consumables (better)": + mix_scenarios( + self._common_baseline(), + { + "HealthSystem": { + "cons_availability_postSwitch": "scenario6", + } + } + ), + + "No growth + Consumables (perfect)": + mix_scenarios( + self._common_baseline(), + { + "HealthSystem": { + "cons_availability_postSwitch": "all", + } + } + ), + + "No growth + Max HS performance": + mix_scenarios( + self._common_baseline(), + { + 'ImprovedHealthSystemAndCareSeekingScenarioSwitcher': { + 'max_healthsystem_function': [False, True], + } + } + ), + + "No growth + Upper bound settings": + mix_scenarios( + self._common_baseline(), + { + 'ImprovedHealthSystemAndCareSeekingScenarioSwitcher': { + 'max_healthsystem_function': [False, False], + }, + "HealthSystem": { + "policy_name": "LCOA_EHP", + "cons_availability_postSwitch": "all", + } + } + ), + + "Historical growth + Upper bound settings": + mix_scenarios( + self._common_baseline(), + { + 'ImprovedHealthSystemAndCareSeekingScenarioSwitcher': { + 'max_healthsystem_function': [False, False], + }, + "HealthSystem": { + "policy_name": "LCOA_EHP", + "cons_availability_postSwitch": "all", + "HR_scaling_by_year_and_officer_type_mode": "historical_cadre_mix", + } + } + ), + + } + + def _hrh_growth_baseline(self) -> Dict: + return mix_scenarios( + self._common_baseline(), + { + "HealthSystem": { + "HR_scaling_by_year_and_officer_type_mode": "historical_uniform", + } + }, + ) + + def _common_baseline(self) -> Dict: + 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": 2020, # <-- transition happens at start of 2020 when HRH starts to grow + + # Normalize the behaviour of Mode 2 + "policy_name": "Naive", # -- *For the alternative scenario of efficient implementation of EHP, otherwise use 'naive'* -- + "tclose_overwrite": 1, + "tclose_days_offset_overwrite": 7, + + # Clarify the consumable availability + "cons_availability": "default", + "cons_availability_postSwitch": "default", + "year_cons_availability_switch": 2020, + + # Clarify the historical HRH growth mode between 2020-2024 + "yearly_HR_scaling_mode": 'no_scaling', + "HR_scaling_by_year_and_officer_type_mode": 'no_historical_growth', + + # Clarify the HRH expansion absorption rate + "HR_expansion_absorption_rate": 1.0, + + }, + # -- *For the alternative scenario of increased demand and improved clinician performance* -- + 'ImprovedHealthSystemAndCareSeekingScenarioSwitcher': { + 'max_healthcare_seeking': [False, False], # <-- switch from False to True mid-way + 'max_healthsystem_function': [False, False], + 'year_of_switch': 2020, + } + }, + ) + + +if __name__ == '__main__': + from tlo.cli import scenario_run + scenario_run([__file__]) diff --git a/src/tlo/methods/healthsystem.py b/src/tlo/methods/healthsystem.py index f8b0f55a03..ea7aa280ba 100644 --- a/src/tlo/methods/healthsystem.py +++ b/src/tlo/methods/healthsystem.py @@ -308,6 +308,33 @@ class HealthSystem(Module): "(factors informed by survey data); and, `custom` (user can freely set these factors as " "parameters in the analysis).", ), + "HR_scaling_by_year_and_officer_type_table": Parameter( + Types.DICT, + "Factors by which capabilities of medical officer types will be scaled at the start of " + "the years considered. This is the imported from an Excel workbook: keys are the worksheet " + "names and values are the worksheets in the format of pd.DataFrames. Additional scenarios can " + "be added by adding worksheets to this workbook: the value of " + "`HR_scaling_by_year_and_officer_type_mode` indicates which sheet is used.", + ), + "HR_scaling_by_year_and_officer_type_mode": Parameter( + Types.STRING, + "Mode of HR scaling considered at the start of the simulation. This corresponds to the name" + "of the worksheet in `ResourceFile_HR_scaling_by_year_and_officer_type.xlsx` that should be used. " + "Options are: " + "`no_historical_growth` (capabilities are scaled by a factor of 1); " + "`historical_uniform` (factors across the cadres are the same, informated by the historical growth rate " + "in total HRH);" + "`historical_cadre_mix` (factors across the cadres are different, indicated by both the overall historical " + "increase rate and the gap_allocation HRH expansion scenario by She et al 2026); and," + "`custom` (user can freely set these factors as parameters in the analysis). " + "Note that if use historical_ modes in this HRH scaling, use no_scaling mode in DynamicHRHScaling " + "to avoid duplicates.", + ), + "HR_expansion_absorption_rate": Parameter( + Types.REAL, + "Rate at which HRH expansion is absorbed into the system, expressed as a fraction of the total HRH " + "expansion per year. A value of 1 meaning 100% of the expansion is absorbed each year.", + ), "HR_scaling_by_district_table": Parameter( Types.DICT, "Factors by which daily capabilities in different districts will be" @@ -694,6 +721,22 @@ def read_consumables(filename): f"{self.parameters['HR_scaling_by_level_and_officer_type_mode']}" ) + self.parameters["HR_scaling_by_year_and_officer_type_table"]: Dict = read_csv_files( + path_to_resourcefiles_for_healthsystem + / "human_resources" + / "scaling_capabilities" + / "ResourceFile_HR_scaling_by_year_and_officer_type", + files=None, # all sheets read in + ) + # Ensure the mode of HR scaling to be considered in included in the tables loaded + assert ( + self.parameters["HR_scaling_by_year_and_officer_type_mode"] + in self.parameters["HR_scaling_by_year_and_officer_type_table"] + ), ( + f"Value of `HR_scaling_by_year_and_officer_type_mode` not recognised: " + f"{self.parameters['HR_scaling_by_year_and_officer_type_mode']}" + ) + self.parameters["HR_scaling_by_district_table"]: Dict = read_csv_files( path_to_resourcefiles_for_healthsystem / "human_resources" @@ -915,6 +958,9 @@ def initialise_simulation(self, sim): # whilst the actual scaling will only take effect from 2011 onwards. sim.schedule_event(DynamicRescalingHRCapabilities(self), Date(sim.date)) + # Schedule recurring event which will rescale daily capabilities annually. + sim.schedule_event(RescalingHRCapabilities_ByYearAndOfficerAndAbsorption(self), Date(sim.date)) + # Schedule the logger to occur at the start of every year sim.schedule_event(HealthSystemLogger(self), Date(sim.date.year, 1, 1)) @@ -2984,6 +3030,45 @@ def apply(self, population): ] +class RescalingHRCapabilities_ByYearAndOfficerAndAbsorption(RegularEvent, PopulationScopeEventMixin): + """This event exists to scale the daily capabilities, with a factor for each Officer Type at each specified year.""" + + def __init__(self, module): + # set the frequency + super().__init__(module, frequency=DateOffset(years=1)) + + # Get the set of scaling_factors that are specified by the 'HR_scaling_by_level_and_officer_type_mode' + # assumption + self.HR_scaling_by_year_and_officer_type_factor = self.module.parameters[ + "HR_scaling_by_year_and_officer_type_table" + ][self.module.parameters["HR_scaling_by_year_and_officer_type_mode"]].set_index("Officer_Category") + + self.HR_scaling_by_year_and_officer_type_factor.columns = ( + self.HR_scaling_by_year_and_officer_type_factor.columns.astype(int) + ) + + self.years_for_scaling = np.array(self.HR_scaling_by_year_and_officer_type_factor.columns) + + def apply(self, population): + + most_recent_year_for_scaling = self.years_for_scaling[self.years_for_scaling <= self.sim.date.year].max() + + pattern = r"FacilityID_(\w+)_Officer_(\w+)" + + for clinic, clinic_cl in self.module._daily_capabilities.items(): + for officer in clinic_cl.keys(): + matches = re.match(pattern, officer) + # Extract officer type + officer_type = matches.group(2) + # Update capabilities by scaling factor and absorption rate + sf = self.HR_scaling_by_year_and_officer_type_factor.at[ + officer_type, most_recent_year_for_scaling + ] + ar = self.module.parameters["HR_expansion_absorption_rate"] if sf >= 1 else 1 + sf_ar_adjusted = 1 + (sf-1) * ar + self.module._daily_capabilities[clinic][officer] *= sf_ar_adjusted + + class RescaleHRCapabilities_ByDistrict(Event, PopulationScopeEventMixin): """This event exists to scale the daily capabilities, with a factor for each district."""