From fe027349a47501e9b792679205407d3aac156567 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Fri, 26 Sep 2025 12:17:30 +0100 Subject: [PATCH 01/14] Changes to accommodate revised scenarios - Note that Xpert costs no longer need to be estimated separately. --- .../roi_analysis_horizontal_vs_vertical.py | 281 ++++++------------ 1 file changed, 92 insertions(+), 189 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index 8d43c24b2f..a74b2660c4 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -51,13 +51,13 @@ # Create folders to store results resourcefilepath = Path("./resources") outputfilepath = Path('./outputs/t.mangal@imperial.ac.uk') -main_figurespath = Path('./outputs/horizontal_v_vertical') +main_figurespath = Path('./outputs/horizontal_v_vertical_firstrevision') if not os.path.exists(main_figurespath): os.makedirs(main_figurespath) # Load result files # ------------------------------------------------------------------------------------------------------------------ -results_folder = get_scenario_outputs('htm_and_hss_runs-2025-01-16T135243Z.py', outputfilepath)[0] +results_folder = get_scenario_outputs('htm_and_hss_runs-2025-09-16T141811Z.py', outputfilepath)[0] # Check can read results from draw=0, run=0 log = load_pickled_dataframes(results_folder, 0, 0) # look at one log (so can decide what to extract) @@ -83,24 +83,54 @@ # Scenarios # Full list of scenarios used in the manuscript -all_manuscript_scenarios = {0: "Baseline", 1: "Pessimistic HRH Scale-up", 2: "Historical HRH Scale-up", +main_manuscript_scenarios = {0: "Baseline", + 1: "Pessimistic HRH Scale-up", 2: "Historical HRH Scale-up", 3: "Optimistic HRH Scale-up", - 4: "Primary Healthcare Workforce Scale-up", 5: "Consumables Increased to 75th Percentile", - 6: "Consumables Increased to HIV levels", 7: "Consumables Increased to EPI Levels", - 8: "HSS Expansion Package", - 9: "HIV Program Scale-up Without HSS Expansion", - 17: "HIV Program Scale-up With HSS Expansion Package", - 18: "TB Program Scale-up Without HSS Expansion", - 26: "TB Program Scale-up With HSS Expansion Package", - 27: "Malaria Program Scale-up Without HSS Expansion", - 35: "Malaria Program Scale-up With HSS Expansion Package", - 36: "HTM Programs Scale-up Without HSS Expansion", - 44: "HTM Programs Scale-up With HSS Expansion Package"} -# 39: "HTM Program Scale-up With HRH Scale-up (6%)", 41: "HTM Program Scale-up With Consumables at 75th Percentile", + 4: "Consumables Increased to 75th Percentile", + 5: "Consumables Increased to HIV levels", 6: "Consumables Increased to EPI Levels", + 7: "HSS Expansion Package", + 8: "HIV Program Scale-up Without HSS Expansion", + 15: "HIV Program Scale-up With HSS Expansion Package", + 16: "TB Program Scale-up Without HSS Expansion", + 23: "TB Program Scale-up With HSS Expansion Package", + 24: "Malaria Program Scale-up Without HSS Expansion", + 31: "Malaria Program Scale-up With HSS Expansion Package", + 32: "HTM Programs Scale-up Without HSS Expansion", + 39: "HTM Programs Scale-up With HSS Expansion Package"} + +frontier_scenarios = {40:"HIV + TB Scale-up Without HSS Expansion", + 41:"HIV + TB Scale-up With Historical HRH Scale-up", + 42:"HIV + TB Scale-up With Consumables Increased to 75th Percentile", + 43:"HIV + TB Scale-up With HSS Expansion", + 44:"HIV + Malaria Scale-up Without HSS Expansion", + 45:"HIV + Malaria Scale-up With Historical HRH Scale-up", + 46:"HIV + Malaria Scale-up With Consumables Increased to 75th Percentile", + 47:"HIV + Malaria Scale-up With HSS Expansion", + 48:"TB + Malaria Scale-up Without HSS Expansion", + 49:"TB + Malaria Scale-up With Historical HRH Scale-up", + 50:"TB + Malaria Scale-up With Consumables Increased to 75th Percentile", + 51:"TB + Malaria Scale-up With HSS Expansion"} + +all_manuscript_scenarios = {**main_manuscript_scenarios, **frontier_scenarios} + +all_manuscript_scenarios_reverse = {v: k for k, v in all_manuscript_scenarios.items()} + +baseline = all_manuscript_scenarios_reverse.get("Baseline") +horizontal_hss = all_manuscript_scenarios_reverse.get("HSS Expansion Package") +vertical_hiv = all_manuscript_scenarios_reverse.get("HIV Program Scale-up Without HSS Expansion") +vertical_tb = all_manuscript_scenarios_reverse.get("TB Program Scale-up Without HSS Expansion") +vertical_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up Without HSS Expansion") +vertical_htm = all_manuscript_scenarios_reverse.get("HTM Programs Scale-up Without HSS Expansion") +diagonal_hiv = all_manuscript_scenarios_reverse.get("HIV Program Scale-up Without HSS Expansion") +diagonal_tb = all_manuscript_scenarios_reverse.get("TB Program Scale-up Without HSS Expansion") +diagonal_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up Without HSS Expansion") +diagonal_htm = all_manuscript_scenarios_reverse.get("HTM Programs Scale-up With HSS Expansion Package") # Use letters instead of full scenario name for figures all_manuscript_scenarios_substitutedict = {0: "0", 1: "A", 2: "B", 3: "C", 4: "D", 5: "E", 6: "F", 7: "G", 8: "H", - 9: "I", 17: "J", 18: "K", 26: "L", 27: "M", 35: "N", 36: "O", 44: "P"} + 15: "I", 16: "J", 23: "K", 24: "L", 31: "M", 32: "N", 39: "O", + 40: "P", 41: "Q", 42: "R", 43: "S", 44: "T", 45: "U", 46: "V", 47: "W", 48: "X", + 49: "Y", 50: "Z", 51: "2"} # Function to adjust color brightness (lighten/darken) @@ -118,7 +148,6 @@ def adjust_color(hex_color, factor=0.5): "Pessimistic HRH Scale-up": adjust_color('#9e0142', 0.5), "Historical HRH Scale-up": adjust_color('#9e0142', 0.5), "Optimistic HRH Scale-up": adjust_color('#9e0142', 0.5), - "Primary Healthcare Workforce Scale-up": adjust_color('#9e0142', 0.5), # Consumables scenarios "Consumables Increased to 75th Percentile": adjust_color('#9e0142', 0.5), "Consumables Increased to HIV levels": adjust_color('#9e0142', 0.5), @@ -468,139 +497,13 @@ def melt_and_label_malaria_scaleup_cost(_df, label): _total_input_costs=input_costs, _relevant_period_for_costing=relevant_period_for_costing) - - def estimate_xpert_costs(_results_folder, _relevant_period_for_costing): - # Load health spending projections - unit_costs = load_unit_cost_assumptions(resourcefilepath) - unit_price_consumable = unit_costs["consumables"] - - # Read parameters for consumables costs - # Load consumables cost data - unit_price_consumable = unit_price_consumable[unit_price_consumable['Item_Code'].notna()] - - # Add cost of Xpert consumables which was missed in the current analysis - def get_counts_of_items_requested(_df): - counts_of_used = defaultdict(lambda: defaultdict(int)) - counts_of_available = defaultdict(lambda: defaultdict(int)) - counts_of_not_available = defaultdict(lambda: defaultdict(int)) - - for _, row in _df.iterrows(): - date = row['date'] - for item, num in row['Item_Used'].items(): - counts_of_used[date][item] += num - for item, num in row['Item_Available'].items(): - counts_of_available[date][item] += num - for item, num in row['Item_NotAvailable'].items(): - counts_of_not_available[date][item] += num - - used_df = pd.DataFrame(counts_of_used).fillna(0).astype(int).stack().rename('Used') - available_df = pd.DataFrame(counts_of_available).fillna(0).astype(int).stack().rename('Available') - not_available_df = pd.DataFrame(counts_of_not_available).fillna(0).astype(int).stack().rename( - 'Not_Available') - - # Combine the two dataframes into one series with MultiIndex (date, item, availability_status) - combined_df = pd.concat([used_df, available_df, not_available_df], axis=1).fillna(0).astype(int) - - # Convert to a pd.Series, as expected by the custom_generate_series function - return combined_df.stack() - - cons_req = extract_results( - _results_folder, - module='tlo.methods.healthsystem.summary', - key='Consumables', - custom_generate_series=get_counts_of_items_requested, - do_scaling=True) - keep_xpert = cons_req.index.get_level_values(0) == '187' - keep_instances_logged_as_not_available = cons_req.index.get_level_values(2) == 'Not_Available' - cons_req = cons_req[keep_xpert & keep_instances_logged_as_not_available] - cons_req = cons_req.reset_index() - - # Keep only relevant draws - # Filter columns based on keys from all_manuscript_scenarios - col_subset = [col for col in cons_req.columns if - ((col[0] in all_manuscript_scenarios.keys()) | (col[0] == 'level_1'))] - # Keep only the relevant columns - cons_req = cons_req[col_subset] - - def transform_cons_requested_for_costing(_df, date_column): - _df['year'] = pd.to_datetime(_df[date_column]).dt.year - - # Validate that all necessary years are in the DataFrame - if not set(_relevant_period_for_costing).issubset(_df['year'].unique()): - raise ValueError("Some years are not recorded in the dataset.") - - # Filter for relevant years and return the total population as a Series - return _df.loc[ - _df['year'].between(min(_relevant_period_for_costing), max(_relevant_period_for_costing))].drop( - columns=date_column).set_index( - 'year') - - xpert_cost_per_cartridge = unit_price_consumable[unit_price_consumable.Item_Code == 187][ - 'Price_per_unit'] - xpert_availability_adjustment = 0.77 - # note that as per 2018 Openlmis mean availability was 0.5, 0.77, 0.809, 0.861, 0.744 during Jul, Sep - Dec. - # The median value is 0.77 - for the rest of the months availability was close to 0 - - xpert_dispensed_cost = (transform_cons_requested_for_costing(_df=cons_req, - date_column=('level_1', '')) - * xpert_availability_adjustment - * xpert_cost_per_cartridge.iloc[0]) - draws_with_positive_xpert_costs = ( - input_costs[(input_costs.cost_subgroup == 'Xpert') & (input_costs.cost > 0)].groupby('draw')[ - 'cost'].sum().reset_index()['draw'] - .unique()).tolist() - - def melt_and_label_xpert_cost(_df): - multi_index = pd.MultiIndex.from_tuples(_df.columns) - _df.columns = multi_index - - # reshape dataframe and assign 'draw' and 'run' as the correct column headers - melted_df = pd.melt(_df.reset_index(), id_vars=['year']).rename( - columns={'variable_0': 'draw', 'variable_1': 'run'}) - # For draws where the costing is already correct, set additional costs to 0 - melted_df.loc[melted_df.draw.isin(draws_with_positive_xpert_costs), 'value'] = 0 - # Replace item_code with consumable_name_tlo - melted_df['cost_category'] = 'medical consumables' - melted_df['cost_subgroup'] = 'Xpert' - melted_df['Facility_Level'] = 'all' - melted_df = melted_df.rename(columns={'value': 'cost'}) - - # Replicate and estimate cost of consumables stocked and supply chain costs - df_with_all_cost_subcategories = pd.concat([melted_df] * 3, axis=0, ignore_index=True) - # Define cost subcategory values - cost_categories = ['cost_of_consumables_dispensed', 'cost_of_excess_consumables_stocked', 'supply_chain'] - # Assign values to the new 'cost_subcategory' column - df_with_all_cost_subcategories['cost_subcategory'] = np.tile(cost_categories, len(melted_df)) - # The excess stock ratio of Xpert as per 2018 LMIS data is 0.125833 - df_with_all_cost_subcategories.loc[df_with_all_cost_subcategories[ - 'cost_subcategory'] == 'cost_of_excess_consumables_stocked', 'cost'] *= 0.125833 - # Supply chain costs are 0.12938884672119721 of the cost of dispensed + stocked - df_with_all_cost_subcategories.loc[ - df_with_all_cost_subcategories['cost_subcategory'] == 'supply_chain', 'cost'] *= ( - 0.12938884672119721 * (1 + 0.125833)) - return df_with_all_cost_subcategories - - xpert_total_cost = melt_and_label_xpert_cost(xpert_dispensed_cost) - xpert_total_cost.to_csv('./outputs/horizontal_v_vertical/xpert_cost.csv') - xpert_total_cost_discounted = apply_discounting_to_cost_data(xpert_total_cost, - _discount_rate=discount_rate_cost, - _initial_year=_relevant_period_for_costing[0]) - return xpert_total_cost_discounted - - - print("Appending Xpert costs") - xpert_total_cost = estimate_xpert_costs(_results_folder=results_folder, - _relevant_period_for_costing=relevant_period_for_costing) - - # Update input costs to include Xpert costs - input_costs = pd.concat([input_costs, xpert_total_cost], ignore_index=True) input_costs = input_costs.groupby(['draw', 'run', 'year', 'cost_subcategory', 'Facility_Level', 'cost_subgroup', 'cost_category'])['cost'].sum().reset_index() # Keep costs for relevant draws input_costs = input_costs[input_costs['draw'].isin(list(all_manuscript_scenarios.keys()))] # Extract input_costs for browsing - # input_costs.groupby(['draw', 'run', 'cost_category', 'cost_subcategory', 'cost_subgroup','year'])['cost'].sum().to_csv(figurespath / 'cost_detailed.csv') + #input_costs.groupby(['draw', 'run', 'cost_category', 'cost_subcategory', 'cost_subgroup','year'])['cost'].sum().to_csv(figurespath / 'cost_detailed.csv') # %% # Return on Invesment analysis @@ -631,7 +534,7 @@ def find_difference_relative_to_comparison(_ser: pd.Series, comparison=0) # sets the comparator to draw 0 which is the Actual scenario ).T.iloc[0].unstack()).T incremental_scenario_cost_subset_for_figure = incremental_scenario_cost[ - incremental_scenario_cost.index.get_level_values('draw').isin(list(all_manuscript_scenarios.keys()))] + incremental_scenario_cost.index.get_level_values('draw').isin(list(main_manuscript_scenarios.keys()))] incremental_scenario_cost_summarised = summarize_cost_data(incremental_scenario_cost_subset_for_figure, _metric=chosen_metric) @@ -677,7 +580,7 @@ def get_num_dalys(_df): # Plot DALYs num_dalys_averted_subset_for_figure = num_dalys_averted[ - num_dalys_averted.index.get_level_values('draw').isin(list(all_manuscript_scenarios.keys()))] + num_dalys_averted.index.get_level_values('draw').isin(list(main_manuscript_scenarios.keys()))] num_dalys_averted_summarised = summarize_cost_data(num_dalys_averted_subset_for_figure, _metric=chosen_metric) name_of_plot = f'Incremental DALYs averted compared to baseline {relevant_period_for_costing[0]}-{relevant_period_for_costing[1]}' fig, ax = do_standard_bar_plot_with_ci( @@ -712,7 +615,7 @@ def get_monetary_value_of_incremental_health(_num_dalys_averted, _chosen_value_o dominated_scenarios = icers_summarized['median'] < 0 icers_summarized[dominated_scenarios] = np.nan # The dominanted scenarios are assigned as ICER of NaN icers_summarized_subset_for_figure = icers_summarized[ - icers_summarized.index.get_level_values('draw').isin(list(all_manuscript_scenarios.keys()))] + icers_summarized.index.get_level_values('draw').isin(list(main_manuscript_scenarios.keys()))] # Plot ICERs annotations_icers = [] @@ -780,9 +683,9 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, } # Assign scenario groups - horizontal_scenarios = [1, 2, 3, 4, 5, 6, 7, 8] - vertical_scenarios = [9, 18, 27, 36] - diagonal_scenarios = [17, 26, 35, 44] + horizontal_scenarios = [1, 2, 3, 4, 5, 6, 7] + vertical_scenarios = [8, 16, 24, 32] + diagonal_scenarios = [15, 23, 31, 39] darker_scenarios = {'HSS Expansion Package', 'HTM Programs Scale-up Without HSS Expansion', 'HTM Programs Scale-up With HSS Expansion Package'} @@ -915,7 +818,7 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, do_incremental_cost_and_health_plot(incremental_cost_df=incremental_scenario_cost_summarised, incremental_dalys_df=num_dalys_averted_summarised, - draws_with_icer_labels=[8, 36, 44], + draws_with_icer_labels=[7, 32, 39], figname='cea_plane.png') # 4. Return on Investment @@ -1098,8 +1001,8 @@ def compute_breakeven_with_ci( monetary_benefit_df=get_monetary_value_of_incremental_health(num_dalys_averted, vsly), incremental_cost_df=incremental_scenario_cost, implementation_cost_df=0 * incremental_scenario_cost, - draw_target=44, - draw_compare=36 + draw_target=diagonal_htm, + draw_compare=vertical_htm ) # Breakeven implementation cost for joint HTM diagonal scenario (V joint HTM vertical scenario); # Implementation cost = 138% of Input costs @@ -1108,8 +1011,8 @@ def compute_breakeven_with_ci( monetary_benefit_df=get_monetary_value_of_incremental_health(num_dalys_averted, vsly), incremental_cost_df=incremental_scenario_cost, implementation_cost_df=implementation_cost_upper_limit, - draw_target=44, - draw_compare=36 + draw_target=diagonal_htm, + draw_compare=vertical_htm ) # Breakeven implementation cost for joint HTM diagonal scenario (V horizontal scenario); # Implementation cost = 0 @@ -1118,8 +1021,8 @@ def compute_breakeven_with_ci( monetary_benefit_df=get_monetary_value_of_incremental_health(num_dalys_averted, vsly), incremental_cost_df=incremental_scenario_cost, implementation_cost_df=0 * incremental_scenario_cost, - draw_target=44, - draw_compare=8 + draw_target=diagonal_htm, + draw_compare=horizontal_hss ) # Breakeven implementation cost for joint HTM diagonal scenario (V horizontal scenario); # Implementation cost = 138% of Input costs @@ -1128,8 +1031,8 @@ def compute_breakeven_with_ci( monetary_benefit_df=get_monetary_value_of_incremental_health(num_dalys_averted, vsly), incremental_cost_df=incremental_scenario_cost, implementation_cost_df=implementation_cost_upper_limit, - draw_target=44, - draw_compare=8 + draw_target=diagonal_htm, + draw_compare=horizontal_hss ) if (vsly_fig_suffixes[i] == 'MAIN') & (rates["discounting_scenario"] == 'MAIN (0.03,0.03)'): # The breakeven @@ -1137,7 +1040,7 @@ def compute_breakeven_with_ci( # Generate data for the representation of breakeven implementation costs on the ROI plots additional_horizontal_lines_for_interpretation = [ { - 'y_value': roi_at_0_implementation_cost_dict[36][chosen_metric], + 'y_value': roi_at_0_implementation_cost_dict[vertical_htm][chosen_metric], 'x_start': 0, 'x_end': breakeven_lower_v_htm[0], # where the horizontal line intersects the ROI curve of the diagonal strategy @@ -1146,8 +1049,8 @@ def compute_breakeven_with_ci( 'scenario_label': 'Diagonal versus Vertical (at ASC = 0)' }, { - 'y_value': roi_at_upper_limit_implementation_cost_dict[36][chosen_metric], - 'x_start': implementation_cost_upper_limit_dict[36][chosen_metric] / 1e6, + 'y_value': roi_at_upper_limit_implementation_cost_dict[vertical_htm][chosen_metric], + 'x_start': implementation_cost_upper_limit_dict[vertical_htm][chosen_metric] / 1e6, 'x_end': breakeven_upper_v_htm[0], # where the horizontal line intersects the ROI curve of the diagonal strategy 'label': f"b = ${breakeven_upper_v_htm[0]:.0f}M [${breakeven_upper_v_htm[1]:.0f}M–${breakeven_upper_v_htm[2]:.0f}M]", @@ -1156,7 +1059,7 @@ def compute_breakeven_with_ci( 'y_label_offset': 0.25 }, { - 'y_value': roi_at_0_implementation_cost_dict[8][chosen_metric], + 'y_value': roi_at_0_implementation_cost_dict[horizontal_hss][chosen_metric], 'x_start': 0, 'x_end': breakeven_lower_v_hss[3], # where the horizontal line intersects the ROI curve of the diagonal strategy @@ -1166,8 +1069,8 @@ def compute_breakeven_with_ci( 'y_label_offset': -0.25 }, { - 'y_value': roi_at_upper_limit_implementation_cost_dict[8][chosen_metric], - 'x_start': implementation_cost_upper_limit_dict[8][chosen_metric] / 1e6, + 'y_value': roi_at_upper_limit_implementation_cost_dict[horizontal_hss][chosen_metric], + 'x_start': implementation_cost_upper_limit_dict[horizontal_hss][chosen_metric] / 1e6, 'x_end': breakeven_upper_v_hss[3], # where the horizontal line intersects the ROI curve of the diagonal strategy 'label': f"d = ${breakeven_upper_v_hss[0]:.0f}M [${breakeven_upper_v_hss[1]:.0f}M–${breakeven_upper_v_hss[2]:.0f}M]", @@ -1181,12 +1084,12 @@ def compute_breakeven_with_ci( additional_horizontal_lines_for_interpretation = None legend_switch_for_main_roi_plot = True - draw_colors = {8: '#9e0142', 36: '#fdae61', 44: '#66c2a5'} + draw_colors = {horizontal_hss: '#9e0142', vertical_htm: '#fdae61', diagonal_htm: '#66c2a5'} generate_multiple_scenarios_roi_plot( _monetary_value_of_incremental_health=get_monetary_value_of_incremental_health(num_dalys_averted, _chosen_value_of_life_year=vsly), _incremental_input_cost=incremental_scenario_cost, - _draws=[8, 36, 44], + _draws=[horizontal_hss, vertical_htm, diagonal_htm], _scenario_dict=all_manuscript_scenarios, _outputfilepath=figurespath, _value_of_life_suffix=vsly_fig_suffixes[i], @@ -1205,15 +1108,15 @@ def compute_breakeven_with_ci( f"the diagonal approach provided a higher ROI up to an even higher threshold of " f"${breakeven_upper_v_htm[0]: .2f}[${breakeven_upper_v_htm[1]: .2f} - ${breakeven_upper_v_htm[2]: .2f}] million " f"incremental above service level costs in comparison with the vertical approach, or equal to " - f"{(breakeven_upper_v_htm[0] * 1e6) / incremental_scenario_cost_summarised[44][chosen_metric] * 100: .2f}% of its own incremental service level cost") + f"{(breakeven_upper_v_htm[0] * 1e6) / incremental_scenario_cost_summarised[diagonal_htm][chosen_metric] * 100: .2f}% of its own incremental service level cost") # HIV scenarios with and without HSS - draw_colors = {9: '#fdae61', 17: '#66c2a5'} + draw_colors = {vertical_hiv: '#fdae61', diagonal_hiv: '#66c2a5'} generate_multiple_scenarios_roi_plot( _monetary_value_of_incremental_health=get_monetary_value_of_incremental_health(num_dalys_averted, _chosen_value_of_life_year=chosen_value_of_statistical_life), _incremental_input_cost=incremental_scenario_cost, - _draws=[9, 17], + _draws=[vertical_hiv, diagonal_hiv], _scenario_dict=all_manuscript_scenarios, _metric=chosen_metric, _outputfilepath=figurespath, @@ -1222,12 +1125,12 @@ def compute_breakeven_with_ci( _draw_colors=draw_colors) # TB scenarios with and without HSS - draw_colors = {18: '#fdae61', 26: '#66c2a5'} + draw_colors = {vertical_tb: '#fdae61', diagonal_tb: '#66c2a5'} generate_multiple_scenarios_roi_plot( _monetary_value_of_incremental_health=get_monetary_value_of_incremental_health(num_dalys_averted, _chosen_value_of_life_year=chosen_value_of_statistical_life), _incremental_input_cost=incremental_scenario_cost, - _draws=[18, 26], + _draws=[vertical_tb, diagonal_tb], _scenario_dict=all_manuscript_scenarios, _metric=chosen_metric, _outputfilepath=figurespath, @@ -1237,12 +1140,12 @@ def compute_breakeven_with_ci( _y_axis_lim=30) # Malaria scenarios with and without HSS - draw_colors = {27: '#fdae61', 35: '#66c2a5'} + draw_colors = {vertical_malaria: '#fdae61', diagonal_malaria: '#66c2a5'} generate_multiple_scenarios_roi_plot( _monetary_value_of_incremental_health=get_monetary_value_of_incremental_health(num_dalys_averted, _chosen_value_of_life_year=chosen_value_of_statistical_life), _incremental_input_cost=incremental_scenario_cost, - _draws=[27, 35], + _draws=[vertical_malaria, diagonal_malaria], _scenario_dict=all_manuscript_scenarios, _metric=chosen_metric, _outputfilepath=figurespath, @@ -1304,7 +1207,7 @@ def compute_breakeven_with_ci( lower=0.0) # monetary value - change in costs max_ability_to_pay_for_implementation_subset_for_figure = max_ability_to_pay_for_implementation[ max_ability_to_pay_for_implementation.index.get_level_values('draw').isin( - list(all_manuscript_scenarios.keys()))] + list(main_manuscript_scenarios.keys()))] max_ability_to_pay_for_implementation_summarized = summarize_cost_data( max_ability_to_pay_for_implementation_subset_for_figure, _metric=chosen_metric) @@ -1485,25 +1388,25 @@ def remove_discounting(_df, _discount_rate=0, _year=None): # Baseline do_line_plot_of_cost(_df=input_costs_for_plot_summarized, _cost_category='all', - _year=list_of_relevant_years_for_costing, _draws=[0], + _year=list_of_relevant_years_for_costing, _draws=[baseline], disaggregate_by='cost_category', _outputfilepath=figurespath) # HSS alone do_line_plot_of_cost(_df=input_costs_for_plot_summarized, _cost_category='all', - _year=list_of_relevant_years_for_costing, _draws=[8], + _year=list_of_relevant_years_for_costing, _draws=[horizontal_hss], disaggregate_by='cost_category', _outputfilepath=figurespath) # HTM without HSS do_line_plot_of_cost(_df=input_costs_for_plot_summarized, _cost_category='all', - _year=list_of_relevant_years_for_costing, _draws=[36], + _year=list_of_relevant_years_for_costing, _draws=[vertical_htm], disaggregate_by='cost_category', _outputfilepath=figurespath) # HTM with HSS do_line_plot_of_cost(_df=input_costs_for_plot_summarized, _cost_category='all', - _year=list_of_relevant_years_for_costing, _draws=[44], + _year=list_of_relevant_years_for_costing, _draws=[diagonal_htm], disaggregate_by='cost_category', _outputfilepath=figurespath) @@ -1511,17 +1414,17 @@ def remove_discounting(_df, _discount_rate=0, _year=None): # -------------------------- # ICER results icer_result = convert_results_to_dict(icers_summarized) -hr_scenario_with_lowest_icer = min([1, 2, 3, 4], key=lambda k: icer_result[k][chosen_metric]) -hr_scenario_with_highest_icer = max([1, 2, 3, 4], key=lambda k: icer_result[k][chosen_metric]) -cons_scenario_with_lowest_icer = min([5, 6, 7], key=lambda k: icer_result[k][chosen_metric]) -cons_scenario_with_highest_icer = max([5, 6, 7], key=lambda k: icer_result[k][chosen_metric]) +hr_scenario_with_lowest_icer = min([1, 2, 3], key=lambda k: icer_result[k][chosen_metric]) +hr_scenario_with_highest_icer = max([1, 2, 3], key=lambda k: icer_result[k][chosen_metric]) +cons_scenario_with_lowest_icer = min([4, 5, 6], key=lambda k: icer_result[k][chosen_metric]) +cons_scenario_with_highest_icer = max([4, 5, 6], key=lambda k: icer_result[k][chosen_metric]) # Get DALYs averted with vertical strategy as comparator num_dalys_averted_v_vertical = (-1.0 * pd.DataFrame( find_difference_relative_to_comparison( num_dalys.loc[0], - comparison=36) # sets the comparator to 0 which is the Actual scenario + comparison=vertical_htm) # sets the comparator to 0 which is the Actual scenario ).T.iloc[0].unstack(level='run')) num_dalys_averted_v_vertical = num_dalys_averted_v_vertical[ num_dalys_averted_v_vertical.index.get_level_values('draw').isin( @@ -1531,13 +1434,13 @@ def remove_discounting(_df, _discount_rate=0, _year=None): dalys_averted_v_vertical_result = convert_results_to_dict(num_dalys_averted_v_vertical_summarised) print(f"The ICER of vertical strategy relative to the baseline scenario was " - f"${icer_result[36][chosen_metric]:.2f} [${icer_result[36]['lower']:.2f} - ${icer_result[36]['upper']:.2f}] " + f"${icer_result[vertical_htm][chosen_metric]:.2f} [${icer_result[vertical_htm]['lower']:.2f} - ${icer_result[vertical_htm]['upper']:.2f}] " f"per DALY averted, assuming no additional implementation costs. The ICER of the diagonal HTM with HSS strategy " f"relative to the baseline scenario was " - f"${icer_result[44][chosen_metric]:.2f} [${icer_result[44]['lower']:.2f} - ${icer_result[44]['upper']:.2f}], " + f"${icer_result[diagonal_htm][chosen_metric]:.2f} [${icer_result[diagonal_htm]['lower']:.2f} - ${icer_result[diagonal_htm]['upper']:.2f}], " f"demonstrating that the diagonal strategy was more cost-effective than the vertical one. While the horizontal " f"strategy, HSS expansion, averted " - f"{dalys_averted_v_vertical_result[8][chosen_metric] / 1e6:.2f} [{dalys_averted_v_vertical_result[8]['lower'] / 1e6:.2f} - {dalys_averted_v_vertical_result[8]['upper'] / 1e6:.2f}] " + f"{dalys_averted_v_vertical_result[horizontal_hss][chosen_metric] / 1e6:.2f} [{dalys_averted_v_vertical_result[horizontal_hss]['lower'] / 1e6:.2f} - {dalys_averted_v_vertical_result[horizontal_hss]['upper'] / 1e6:.2f}] " f"million more DALYs than HTM expansion, it did so at a higher cost per DALY averted " - f"(${icer_result[8][chosen_metric]:.2f} [${icer_result[8]['lower']:.2f} - ${icer_result[8]['upper']:.2f}] versus " - f"${icer_result[36][chosen_metric]:.2f} [${icer_result[36]['lower']:.2f} - ${icer_result[36]['upper']:.2f}]), meaning it was less cost-effective.") + f"(${icer_result[horizontal_hss][chosen_metric]:.2f} [${icer_result[horizontal_hss]['lower']:.2f} - ${icer_result[horizontal_hss]['upper']:.2f}] versus " + f"${icer_result[vertical_htm][chosen_metric]:.2f} [${icer_result[vertical_htm]['lower']:.2f} - ${icer_result[vertical_htm]['upper']:.2f}]), meaning it was less cost-effective.") From ce1a216e8130f20dbdd670ffab8b64ee1e0015e8 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Fri, 26 Sep 2025 19:27:19 +0100 Subject: [PATCH 02/14] first take - Frontier plot --- .../roi_analysis_horizontal_vs_vertical.py | 364 ++++++++++++------ 1 file changed, 251 insertions(+), 113 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index a74b2660c4..58928acff6 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -16,6 +16,7 @@ import matplotlib.colors as mcolors import matplotlib.patches as mpatches import matplotlib.pyplot as plt +from scipy.spatial import ConvexHull # for frontier plot import numpy as np import pandas as pd import seaborn as sns @@ -98,7 +99,31 @@ 32: "HTM Programs Scale-up Without HSS Expansion", 39: "HTM Programs Scale-up With HSS Expansion Package"} -frontier_scenarios = {40:"HIV + TB Scale-up Without HSS Expansion", +frontier_scenarios = {9:"HIV Program Scale-up With Pessimistic HRH Scale-up", + 10:"HIV Program Scale-up With Historical HRH Scale-up", + 11:"HIV Program Scale-up With Historical HRH Scale-up", + 12:"HIV Program Scale-up With Consumables Increased to 75th Percentile", + 13:"HIV Program Scale-up With Consumables Increased to HIV levels", + 14:"HIV Program Scale-up With Consumables Increased to EPI Levels", + 17:"TB Program Scale-up With Pessimistic HRH Scale-up", + 18:"TB Program Scale-up With Historical HRH Scale-up", + 19:"TB Program Scale-up With Optimistic HRH Scale-up", + 20:"TB Program Scale-up With Consumables Increased to 75th Percentile", + 21:"TB Program Scale-up With Consumables Increased to HIV levels", + 22:"TB Program Scale-up With Consumables Increased to EPI Levels", + 25:"Malaria Program Scale-up With Pessimistic HRH Scale-up", + 26:"Malaria Program Scale-up With Historical HRH Scale-up", + 27:"Malaria Program Scale-up With Historical HRH Scale-up", + 28:"Malaria Program Scale-up With Consumables Increased to 75th Percentile", + 29:"Malaria Program Scale-up With Consumables Increased to HIV levels", + 30:"Malaria Program Scale-up With Consumables Increased to EPI Levels", + 33:"HTM Scale-up With Pessimistic HRH Scale-up", + 34:"HTM Scale-up With Historical HRH Scale-up", + 35:"HTM Scale-up With Optimistic HRH Scale-up", + 36:"HTM Scale-up With Consumables Increased to 75th Percentile", + 37:"HTM Scale-up With Consumables Increased to HIV levels", + 38:"HTM Scale-up With Consumables Increased to EPI Levels", + 40:"HIV + TB Scale-up Without HSS Expansion", 41:"HIV + TB Scale-up With Historical HRH Scale-up", 42:"HIV + TB Scale-up With Consumables Increased to 75th Percentile", 43:"HIV + TB Scale-up With HSS Expansion", @@ -128,9 +153,7 @@ # Use letters instead of full scenario name for figures all_manuscript_scenarios_substitutedict = {0: "0", 1: "A", 2: "B", 3: "C", 4: "D", 5: "E", 6: "F", 7: "G", 8: "H", - 15: "I", 16: "J", 23: "K", 24: "L", 31: "M", 32: "N", 39: "O", - 40: "P", 41: "Q", 42: "R", 43: "S", 44: "T", 45: "U", 46: "V", 47: "W", 48: "X", - 49: "Y", 50: "Z", 51: "2"} + 15: "I", 16: "J", 23: "K", 24: "L", 31: "M", 32: "N", 39: "O"} # Function to adjust color brightness (lighten/darken) @@ -267,6 +290,83 @@ def do_standard_bar_plot_with_ci(_df, set_colors=None, annotations=None, return fig, ax +# Thi sis a helper function for the upper_convex_frontier_from_dict +def prune_inward_tail(frontier_cd, frontier_draws, tol=1e-9): + """Remove any final points that violate the upper-hull orientation or reduce DALYs.""" + keep = [] + for i, (c, v) in enumerate(frontier_cd): + while len(keep) >= 2: + o = frontier_cd[keep[-2]] + a = frontier_cd[keep[-1]] + # same orientation test we used to build the hull + z = (a[0] - o[0]) * (v - a[1]) - (a[1] - o[1]) * (c - a[0]) + if z > tol: # concave upward -> drop 'a' + keep.pop() + else: + break + # also enforce non-decreasing DALYs with increasing cost + if not keep or v >= frontier_cd[keep[-1], 1] - tol: + keep.append(i) + return frontier_cd[keep], [frontier_draws[i] for i in keep] + +# The function below helps construct a frontier on a cost-effectiveness plane +def upper_convex_frontier_from_dict(points_by_draw, include_collinear=False, tol=1e-12): + """ + points_by_draw: dict[draw_id] -> (cost, dalys), both incremental vs the *same baseline*. + Returns: + draws_on_frontier (ascending cost), frontier_cd (Nx2 array of (cost, dalys)). + Removes strictly and extendedly dominated points (upper convex hull). + """ + items = sorted(((int(d), float(c), float(v)) for d, (c, v) in points_by_draw.items()), + key=lambda t: (t[1], -t[2])) # Sort primarily by cost ascending (t[1]), + # and if two costs are equal, sort the one with higher DALYs first (-t[2]). + + # collapse exact cost ties -> keep the one with max effect + # The monotone-chain hull expects strictly increasing x (cost). + # If two rows have (nearly) the same cost, keep only the one with greater DALYs. + dedup = [] + for t in items: + if dedup and abs(t[1] - dedup[-1][1]) <= tol: + if t[2] > dedup[-1][2] + tol: + dedup[-1] = t + else: + dedup.append(t) + + # Calculate the z component or 2D cross product to make sure that the next point involves a right turn + # (i.e. lower ICER than other subsequent options from the previous point) + def cross(o, a, b): + # cross of OA x AB in (x=cost, y=dalys) + return (a[0]-o[0])*(b[1]-a[1]) - (a[1]-o[1])*(b[0]-a[0]) + # z > 0: the path o→a→b turns left (counter-clockwise). + # z < 0: it turns right (clockwise) + # z ≈ 0: the three points are collinear. + + hull = [] # (cost, dalys, draw) + for d, c, v in dedup: + p = (c, v, d) + while len(hull) >= 2: + o = hull[-2]; a = hull[-1] # These the the two previous points on the hull + z = cross((o[0], o[1]), (a[0], a[1]), (c, v)) + # For the UPPER hull we keep clockwise turns (z < 0). + # If z > 0 (concave) or (collinear and we don't include collinear), pop. + if z > tol or (not include_collinear and abs(z) <= tol): + hull.pop() + else: + break + hull.append(p) + + draws_on_frontier = [d for (_, _, d) in hull] + frontier_cd = np.array([(c, v) for (c, v, _) in hull], dtype=float) + + frontier_cd, draws_on_frontier = prune_inward_tail(frontier_cd, draws_on_frontier, tol=1e-9) + + return draws_on_frontier, frontier_cd + +def incremental_icers(frontier_cd): + """Adjacent-segment ICERs along the convex frontier.""" + dc = np.diff(frontier_cd[:, 0]) + dv = np.diff(frontier_cd[:, 1]) + return dc, dv, dc / dv # Define alternative discount rate sets as a list of dictionaries alternative_discount_rates = [ @@ -674,130 +774,146 @@ def get_monetary_value_of_incremental_health(_num_dalys_averted, _chosen_value_o def do_incremental_cost_and_health_plot(incremental_cost_df, incremental_dalys_df, figname, - draws_with_icer_labels: Optional[list[int]] = None): - # Define colors - color_map_scatter = { - "Horizontal": "#9e0142", - "Vertical": "#fdae61", - "Diagonal": "#66c2a5" - } - - # Assign scenario groups - horizontal_scenarios = [1, 2, 3, 4, 5, 6, 7] - vertical_scenarios = [8, 16, 24, 32] - diagonal_scenarios = [15, 23, 31, 39] - darker_scenarios = {'HSS Expansion Package', 'HTM Programs Scale-up Without HSS Expansion', - 'HTM Programs Scale-up With HSS Expansion Package'} - - # Reorder DataFrames + draws_with_icer_labels: Optional[list[int]] = None, + horizontal_scenarios: Optional[list[int]] = None, + vertical_scenarios: Optional[list[int]] = None, + diagonal_scenarios: Optional[list[int]] = None, + scenario_dict: Optional[dict[int, str]] = None, + # NEW frontier options + overlay_frontier: bool = True, + plot_icers_on_frontier: bool = False, + icer_lookup: Optional[dict[int, str]] = None, + # needs to be provided if plot_icers_on_frontier = True + frontier_color: str = "black", + frontier_linewidth: float = 2.0): + + # Define colors for high level grouping of scenarios + color_map_scatter = {"HSS Investments": "#9e0142", "HTM-focussed Investments": "#fdae61", + "Combined Investments": "#66c2a5"} + horizontal_scenarios = horizontal_scenarios + vertical_scenarios = vertical_scenarios + diagonal_scenarios = diagonal_scenarios + # Draws to highlight, if any + darker_draws = draws_with_icer_labels + darker_set = set(darker_draws) + + # Reorder new_order = horizontal_scenarios + vertical_scenarios + diagonal_scenarios incremental_dalys_df = incremental_dalys_df.loc[new_order] incremental_cost_df = incremental_cost_df.loc[new_order] - # Generate letter labels - scenario_labels = {draw: letter for draw, letter in zip(incremental_dalys_df.index, string.ascii_uppercase)} + # NUMERIC labels 1..N (based on the plotting order above) + scenario_labels = {draw: str(i + 1) for i, draw in enumerate(incremental_dalys_df.index)} - # Extract values for axes + # Extract stats for axes x_median = incremental_dalys_df['median'] / 1e6 x_lower = incremental_dalys_df['lower'] / 1e6 x_upper = incremental_dalys_df['upper'] / 1e6 - y_median = incremental_cost_df['median'] y_lower = incremental_cost_df['lower'] y_upper = incremental_cost_df['upper'] - # Create scatter plot - plt.figure(figsize=(12, 7)) + fig, ax = plt.subplots(figsize=(16, 9)) + texts = [] - # Loop through each draw to plot points with assigned colors + # Scatter + error bars + numeric tiles for draw in incremental_dalys_df.index: - scenario_name = all_manuscript_scenarios.get(draw, f"Scenario {draw}") - - # Assign colors based on scenario category if draw in horizontal_scenarios: - color = color_map_scatter["Horizontal"] + color = color_map_scatter["HSS Investments"] elif draw in vertical_scenarios: - color = color_map_scatter["Vertical"] + color = color_map_scatter["HTM-focussed Investments"] elif draw in diagonal_scenarios: - color = color_map_scatter["Diagonal"] + color = color_map_scatter["Combined Investments"] else: - color = "gray" # Default fallback color - - is_darker = scenario_name in darker_scenarios + color = "gray" - # Scatter plot with error bars - plt.scatter( - x_median[draw], y_median[draw], color=color, - edgecolor='black' if is_darker else 'none', - linewidth=3.5 if is_darker else 0, - label=scenario_labels[draw] - ) + is_darker = draw in darker_set - plt.errorbar( - x_median[draw], y_median[draw], - xerr=[[x_median[draw] - x_lower[draw]], [x_upper[draw] - x_median[draw]]], - yerr=[[y_median[draw] - y_lower[draw]], [y_upper[draw] - y_median[draw]]], - fmt='o', color=color, alpha=0.6, capsize=5, elinewidth=1 - ) + # Plot the central values + ax.scatter(x_median[draw], y_median[draw], + s=36, color=color, + edgecolor='black' if is_darker else 'none', + linewidth=2.0 if is_darker else 0.0, zorder=3) - # Add letter labels above the dots - # Show the draw label (e.g., "A") above the dot - texts = [] - # Add text with no manual offset + # Plot the error bars + ax.errorbar(x_median[draw], y_median[draw], + xerr=[[x_median[draw] - x_lower[draw]], [x_upper[draw] - x_median[draw]]], + yerr=[[y_median[draw] - y_lower[draw]], [y_upper[draw] - y_median[draw]]], + fmt='none', color=color, alpha=0.6, capsize=4, elinewidth=1, zorder=2) + # Add number labels above the dots texts.append( - plt.text( - x_median[draw], - y_median[draw] + 3e7, # adjust upward offset - scenario_labels[draw], - fontsize=8, - ha='center', - va='bottom', - color='white', - weight='bold', - bbox=dict(facecolor=color, edgecolor='none', boxstyle='round,pad=0.2', alpha=0.9) - )) + ax.text(x_median[draw], y_median[draw] + 1e7, + scenario_labels[draw], + fontsize=8, ha='center', va='bottom', color='white', weight='bold', + bbox=dict(facecolor=color, edgecolor='none', boxstyle='round,pad=0.2', alpha=0.9), + zorder=4) + ) # If this draw should show ICER, add it below the dot - if draws_with_icer_labels: - icer_label_offset = 1.2 - if draw in draws_with_icer_labels: - icer_text = icer_lookup.get(draw, "") - texts.append(plt.text( - x_median[draw] + icer_label_offset, - y_median[draw] + 25e7, - f"$ICER_{{{scenario_labels[draw]}}} =$\n{icer_text}", - fontsize=7.5, - weight='bold', - ha='center', - va='top', - color='black' - )) - adjust_text(texts, arrowprops=dict(arrowstyle='->')) - - # Manually create legend - legend_labels = [] - dummy_handles = [] + if draws_with_icer_labels and draw in draws_with_icer_labels: + icer_text = icer_lookup.get(draw, "") + texts.append(ax.text( + x_median[draw] + 0.5, y_median[draw] + 0.1, + f"$ICER_{{{scenario_labels[draw]}}} =$\n{icer_text}", + fontsize=7.5, weight='bold', ha='center', va='top', color='black', zorder=4 + )) + + if texts: + adjust_text(texts, arrowprops=dict(arrowstyle='->')) + + # --- Frontier overlay (upper convex hull in (cost, dalys) vs baseline) --- + if overlay_frontier: + # Build dict in (cost, dalys) using *medians vs baseline* + points_by_draw = { + int(d): (float(y_median[d]), float(x_median[d] * 1e6)) # (cost, dalys) + for d in incremental_dalys_df.index + } + # add baseline origin because we're dealing with incremental costs + points_by_draw.setdefault(0, (0.0, 0.0)) + + frontier_draws, frontier_cd = upper_convex_frontier_from_dict(points_by_draw, include_collinear=False) + + # convert for plotting on scaled axes (x=DALYs millions, y=cost billions) + fx = frontier_cd[:, 1] / 1e6 + fy = frontier_cd[:, 0] + (frontier_line,) = ax.plot(fx, fy, '-', lw=frontier_linewidth, + color=frontier_color, label='Frontier (median)', zorder=5) + ax.scatter(fx, fy, s=24, facecolors='white', edgecolors=frontier_color, + linewidths=1.5, zorder=6) + + if plot_icers_on_frontier and len(frontier_cd) >= 2: + dc, dv, icers = incremental_icers(frontier_cd) + for i in range(len(icers)): + mx = (fx[i] + fx[i + 1]) / 2 + 1.5 + my = (fy[i] + fy[i + 1]) / 2 - 0.5 + ax.text(mx, my, f"ICER compared to\n previous scenario = ${icers[i]:.2f}", fontsize=7, ha='center', + va='bottom', color='black', + bbox=dict( + facecolor='white', # background color + edgecolor=frontier_color, + boxstyle='round,pad=0.2', + linewidth=0.8, + alpha=0.9 + ), + zorder=10) + + # --- Legend using numeric labels --- + legend_labels, dummy_handles = [], [] scenario_categories = { - "Horizontal Approach Scenarios": horizontal_scenarios, - "Vertical Approach Scenarios": vertical_scenarios, - "Diagonal Approach Scenarios": diagonal_scenarios + "HSS Investments": horizontal_scenarios, + "HTM-focussed Investments": vertical_scenarios, + "Combined Investments": diagonal_scenarios } - for category, draws in scenario_categories.items(): - legend_labels.append(category) # Add category header - dummy_handles.append(mpatches.Patch(color="white", label="")) # Invisible spacing patch - + legend_labels.append(category) + dummy_handles.append(mpatches.Patch(color="white", label="")) for draw in draws: - scenario_name = all_manuscript_scenarios.get(draw, f"Scenario {draw}") - color = ( - color_map_scatter["Horizontal"] if draw in horizontal_scenarios else - color_map_scatter["Vertical"] if draw in vertical_scenarios else - color_map_scatter["Diagonal"] if draw in diagonal_scenarios else - "gray" - ) - is_darker = scenario_name in darker_scenarios - + scenario_name = scenario_dict.get(draw, f"Scenario {draw}") + color = (color_map_scatter["HSS Investments"] if draw in horizontal_scenarios else + color_map_scatter["HTM-focussed Investments"] if draw in vertical_scenarios else + color_map_scatter["Combined Investments"] if draw in diagonal_scenarios else "gray") + is_darker = draw in darker_set legend_labels.append(f"{scenario_labels[draw]}: {scenario_name}") dummy_handles.append( plt.Line2D([0], [0], marker='o', color='w', @@ -805,21 +921,43 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, markersize=8, markeredgewidth=2 if is_darker else 0) ) - plt.legend(dummy_handles, legend_labels, loc='upper left', fontsize=8, title="Scenario Key", frameon=True) - - # Labels and title - plt.xlabel("DALYs Averted, millions") - plt.ylabel("Incremental Scenario Cost, billions (USD)") - plt.grid(True) - - plt.savefig(figurespath / figname) - plt.close() - - - do_incremental_cost_and_health_plot(incremental_cost_df=incremental_scenario_cost_summarised, - incremental_dalys_df=num_dalys_averted_summarised, + ax.legend(dummy_handles, legend_labels, bbox_to_anchor=(1.1, 1), loc='upper left', fontsize=8, + title="Scenario Key", frameon=True) + + ax.set_xlabel("DALYs Averted, millions") + ax.set_ylabel("Incremental Scenario Cost, billions (USD)") + ax.grid(True) + fig.tight_layout() + plt.savefig(figurespath / figname, dpi=300) + plt.close(fig) + + + do_incremental_cost_and_health_plot(incremental_cost_df=summarize_cost_data(incremental_scenario_cost, chosen_metric), + incremental_dalys_df=summarize_cost_data(num_dalys_averted, chosen_metric), + draws_with_icer_labels=[1], + horizontal_scenarios=[2, 4, 7], + vertical_scenarios=[8, 16, 24, 32, 40, 44, 48], + diagonal_scenarios=[10, 12, 15, + 18, 20, 23, + 26, 28, 31, + 34, 36, 39, + 41, 42, 43, 45, 46, 47, 49, 50, 51], + plot_icers_on_frontier=True, + scenario_dict=all_manuscript_scenarios, + figname='cea_plane_frontier.png', + icer_lookup=icer_lookup, + overlay_frontier=True) + + do_incremental_cost_and_health_plot(incremental_cost_df=summarize_cost_data(incremental_scenario_cost, chosen_metric), + incremental_dalys_df=summarize_cost_data(num_dalys_averted, chosen_metric), draws_with_icer_labels=[7, 32, 39], - figname='cea_plane.png') + horizontal_scenarios=[1, 2, 3, 4, 5, 6, 7], + vertical_scenarios=[8, 16, 24, 32], + diagonal_scenarios=[15, 23, 31, 39], + scenario_dict=all_manuscript_scenarios, + figname='cea_plane.png', + icer_lookup=icer_lookup, + overlay_frontier=False) # 4. Return on Investment # ---------------------------------------------------- From 1a502d09dfac4ff2e722abc40becbaf7b0d8cff1 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Sun, 28 Sep 2025 15:51:33 +0100 Subject: [PATCH 03/14] corrections to Scenario dictionary --- .../roi_analysis_horizontal_vs_vertical.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index 58928acff6..dafe951e82 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -101,7 +101,7 @@ frontier_scenarios = {9:"HIV Program Scale-up With Pessimistic HRH Scale-up", 10:"HIV Program Scale-up With Historical HRH Scale-up", - 11:"HIV Program Scale-up With Historical HRH Scale-up", + 11:"HIV Program Scale-up With Optimistic HRH Scale-up", 12:"HIV Program Scale-up With Consumables Increased to 75th Percentile", 13:"HIV Program Scale-up With Consumables Increased to HIV levels", 14:"HIV Program Scale-up With Consumables Increased to EPI Levels", @@ -113,7 +113,7 @@ 22:"TB Program Scale-up With Consumables Increased to EPI Levels", 25:"Malaria Program Scale-up With Pessimistic HRH Scale-up", 26:"Malaria Program Scale-up With Historical HRH Scale-up", - 27:"Malaria Program Scale-up With Historical HRH Scale-up", + 27:"Malaria Program Scale-up With Optimistic HRH Scale-up", 28:"Malaria Program Scale-up With Consumables Increased to 75th Percentile", 29:"Malaria Program Scale-up With Consumables Increased to HIV levels", 30:"Malaria Program Scale-up With Consumables Increased to EPI Levels", @@ -146,9 +146,9 @@ vertical_tb = all_manuscript_scenarios_reverse.get("TB Program Scale-up Without HSS Expansion") vertical_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up Without HSS Expansion") vertical_htm = all_manuscript_scenarios_reverse.get("HTM Programs Scale-up Without HSS Expansion") -diagonal_hiv = all_manuscript_scenarios_reverse.get("HIV Program Scale-up Without HSS Expansion") -diagonal_tb = all_manuscript_scenarios_reverse.get("TB Program Scale-up Without HSS Expansion") -diagonal_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up Without HSS Expansion") +diagonal_hiv = all_manuscript_scenarios_reverse.get("HIV Program Scale-up With HSS Expansion Package") +diagonal_tb = all_manuscript_scenarios_reverse.get("TB Program Scale-up With HSS Expansion Package") +diagonal_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up With HSS Expansion Package") diagonal_htm = all_manuscript_scenarios_reverse.get("HTM Programs Scale-up With HSS Expansion Package") # Use letters instead of full scenario name for figures From a05d98244b53384fbf40e6b3157ccb9b03f4fc3e Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Sun, 28 Sep 2025 22:55:46 +0100 Subject: [PATCH 04/14] improve legend positioning in Cost effectiveness plane plots --- .../roi_analysis_horizontal_vs_vertical.py | 26 +++++++++++++------ 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index dafe951e82..b6cd3cd9af 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -785,7 +785,8 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, icer_lookup: Optional[dict[int, str]] = None, # needs to be provided if plot_icers_on_frontier = True frontier_color: str = "black", - frontier_linewidth: float = 2.0): + frontier_linewidth: float = 2.0, + legend_position: Literal['inset', 'offset'] = "offset"): # Define colors for high level grouping of scenarios color_map_scatter = {"HSS Investments": "#9e0142", "HTM-focussed Investments": "#fdae61", @@ -885,10 +886,10 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, if plot_icers_on_frontier and len(frontier_cd) >= 2: dc, dv, icers = incremental_icers(frontier_cd) for i in range(len(icers)): - mx = (fx[i] + fx[i + 1]) / 2 + 1.5 - my = (fy[i] + fy[i + 1]) / 2 - 0.5 - ax.text(mx, my, f"ICER compared to\n previous scenario = ${icers[i]:.2f}", fontsize=7, ha='center', - va='bottom', color='black', + mx = (fx[i] + fx[i + 1]) / 2 + my = (fy[i] + fy[i + 1]) / 2 + 0.7 + ax.text(mx, my, f"ICER = ${icers[i]:.2f} \nper DALY averted", fontsize=7, ha='center', va='bottom', + color='black', bbox=dict( facecolor='white', # background color edgecolor=frontier_color, @@ -921,7 +922,14 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, markersize=8, markeredgewidth=2 if is_darker else 0) ) - ax.legend(dummy_handles, legend_labels, bbox_to_anchor=(1.1, 1), loc='upper left', fontsize=8, + if (legend_position == "inset"): + bbox_to_anchor = (0, 1) + elif (legend_position == "offset"): + bbox_to_anchor = (1.1, 1) + else: + bbox_to_anchor = (0, 1) + + ax.legend(dummy_handles, legend_labels, bbox_to_anchor=bbox_to_anchor, loc='upper left', fontsize=8, title="Scenario Key", frameon=True) ax.set_xlabel("DALYs Averted, millions") @@ -946,7 +954,8 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, scenario_dict=all_manuscript_scenarios, figname='cea_plane_frontier.png', icer_lookup=icer_lookup, - overlay_frontier=True) + overlay_frontier=True, + legend_position="offset") do_incremental_cost_and_health_plot(incremental_cost_df=summarize_cost_data(incremental_scenario_cost, chosen_metric), incremental_dalys_df=summarize_cost_data(num_dalys_averted, chosen_metric), @@ -957,7 +966,8 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, scenario_dict=all_manuscript_scenarios, figname='cea_plane.png', icer_lookup=icer_lookup, - overlay_frontier=False) + overlay_frontier=False, + legend_position="inset") # 4. Return on Investment # ---------------------------------------------------- From 2a467dc40115f2e2257ebccdc42844df71971e16 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Mon, 29 Sep 2025 17:34:46 +0100 Subject: [PATCH 05/14] remove horizontal lines from ROI plots and formatting changes to CEA Plane --- .../roi_analysis_horizontal_vs_vertical.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index b6cd3cd9af..3119733cbf 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -846,8 +846,8 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, texts.append( ax.text(x_median[draw], y_median[draw] + 1e7, scenario_labels[draw], - fontsize=8, ha='center', va='bottom', color='white', weight='bold', - bbox=dict(facecolor=color, edgecolor='none', boxstyle='round,pad=0.2', alpha=0.9), + fontsize=12, ha='center', va='bottom', color='white', weight='bold', + bbox=dict(facecolor=color, edgecolor='none', boxstyle='round,pad=0.3', alpha=0.9), zorder=4) ) @@ -857,7 +857,7 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, texts.append(ax.text( x_median[draw] + 0.5, y_median[draw] + 0.1, f"$ICER_{{{scenario_labels[draw]}}} =$\n{icer_text}", - fontsize=7.5, weight='bold', ha='center', va='top', color='black', zorder=4 + fontsize=12, weight='bold', ha='center', va='top', color='black', zorder=4 )) if texts: @@ -929,7 +929,7 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, else: bbox_to_anchor = (0, 1) - ax.legend(dummy_handles, legend_labels, bbox_to_anchor=bbox_to_anchor, loc='upper left', fontsize=8, + ax.legend(dummy_handles, legend_labels, bbox_to_anchor=bbox_to_anchor, loc='upper left', fontsize=9.2, title="Scenario Key", frameon=True) ax.set_xlabel("DALYs Averted, millions") @@ -1244,7 +1244,7 @@ def compute_breakeven_with_ci( _metric=chosen_metric, _year_suffix=f' ({str(relevant_period_for_costing[0])} - {str(relevant_period_for_costing[1])})', _projected_health_spending=projected_health_spending_baseline, - _additional_horizontal_lines_for_interpretation=additional_horizontal_lines_for_interpretation, + _additional_horizontal_lines_for_interpretation=None, _draw_colors=draw_colors, show_title_and_legend=legend_switch_for_main_roi_plot) From 39bf7536bb78c183f80fdbc8ca0e6e925a6551af Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Mon, 29 Sep 2025 18:00:06 +0100 Subject: [PATCH 06/14] remove horizontal lines from ROI plots and formatting changes to CEA Plane + Add ROI bar plots at two implementation cost assumptions --- .../roi_analysis_horizontal_vs_vertical.py | 59 +++++++++++++++++-- 1 file changed, 55 insertions(+), 4 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index 3119733cbf..cfecba1eac 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -223,8 +223,15 @@ def do_standard_bar_plot_with_ci(_df, set_colors=None, annotations=None, xticks = {(i + 0.5): k for i, k in enumerate(_df.index)} - if set_colors: - colors = [color_map.get(series, 'grey') for series in _df.index] + if set_colors is not None: + # dict mapping -> use index keys; list/tuple/Series -> use as-is + if isinstance(set_colors, dict): + colors = [set_colors.get(k, 'grey') for k in _df.index] + # Optional debug: + # missing = [k for k in _df.index if k not in set_colors] + # if missing: print("No color for:", missing) + else: + colors = list(set_colors) else: cmap = sns.color_palette('Spectral', as_cmap=True) rescale = lambda y: (y - np.min(y)) / (np.max(y) - np.min(y)) # noqa: E731 @@ -1078,8 +1085,8 @@ def convert_results_to_dict(_df): - implementation_cost_upper_limit) roi_at_upper_limit_implementation_cost = benefit_at_0_implementation_cost.div( abs(incremental_scenario_cost + implementation_cost_upper_limit)) - roi_at_upper_limit_implementation_cost_dict = convert_results_to_dict( - summarize_cost_data(roi_at_upper_limit_implementation_cost, _metric=chosen_metric)) + roi_at_upper_limit_implementation_cost_summarized = summarize_cost_data(roi_at_upper_limit_implementation_cost, _metric=chosen_metric) + roi_at_upper_limit_implementation_cost_dict = convert_results_to_dict(roi_at_upper_limit_implementation_cost_summarized) # Create a function to generate threshold (or maximum) implementation costs which the scneario with the higher ROI @@ -1348,6 +1355,50 @@ def compute_breakeven_with_ci( roi_table_small.to_csv(figurespath / f'tabulated_roi_for_manuscript_{roi_table_label[i]}.csv', index=False) i += 1 + # ROI Bar plots + chosen_draws = [horizontal_hss, vertical_htm, diagonal_htm] + roi_at_0_implementation_cost_for_figure = roi_at_0_implementation_cost_summarized[roi_at_0_implementation_cost_summarized.index.isin(chosen_draws)] + roi_at_upper_limit_implementation_cost_for_figure = roi_at_upper_limit_implementation_cost_summarized[roi_at_upper_limit_implementation_cost_summarized.index.isin(chosen_draws)] + draw_colors = {horizontal_hss: '#9e0142', vertical_htm: '#fdae61', diagonal_htm: '#66c2a5'} + + name_of_plot = f'ROI assuming zero implementation costs {relevant_period_for_costing[0]}-{relevant_period_for_costing[1]}' + fig, ax = do_standard_bar_plot_with_ci( + (roi_at_0_implementation_cost_for_figure), + annotations=[ + f"{row['median']:.2f} ({row['lower'] :.2f}- {row['upper']:.2f})" + for _, row in roi_at_0_implementation_cost_for_figure.iterrows() + ], + xticklabels_horizontal_and_wrapped=False, + put_labels_in_legend=True, + offset=0.2, + set_colors = draw_colors + ) + ax.set_title(name_of_plot) + ax.set_ylabel('Return on Investment') + ax.set_ylim(bottom=0) + fig.tight_layout() + fig.savefig(figurespath / name_of_plot.replace(' ', '_').replace(',', '')) + plt.close(fig) + + name_of_plot = f'ROI assuming non-zero implementation costs {relevant_period_for_costing[0]}-{relevant_period_for_costing[1]}' + fig, ax = do_standard_bar_plot_with_ci( + (roi_at_upper_limit_implementation_cost_for_figure), + annotations=[ + f"{row['median']:.2f} ({row['lower'] :.2f}- {row['upper']:.2f})" + for _, row in roi_at_upper_limit_implementation_cost_for_figure.iterrows() + ], + xticklabels_horizontal_and_wrapped=False, + put_labels_in_legend=True, + offset=0.2, + set_colors = draw_colors + ) + ax.set_title(name_of_plot) + ax.set_ylabel('Return on Investment') + ax.set_ylim(bottom=0) + fig.tight_layout() + fig.savefig(figurespath / name_of_plot.replace(' ', '_').replace(',', '')) + plt.close(fig) + # 5. Plot Maximum ability-to-pay at CET # ---------------------------------------------------- max_ability_to_pay_for_implementation = (get_monetary_value_of_incremental_health(num_dalys_averted, From 72ae25b172c019e61884de0aca082bd4920e20f8 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Wed, 1 Oct 2025 18:56:08 +0100 Subject: [PATCH 07/14] move all the functions up + fix calculation of roi at upper bound ASC + formatting fixes to plots --- .../roi_analysis_horizontal_vs_vertical.py | 1165 ++++++++--------- 1 file changed, 573 insertions(+), 592 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index cfecba1eac..deda3879cb 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -7,18 +7,17 @@ """ import datetime import os -import string import textwrap -from collections import defaultdict from pathlib import Path -from typing import Optional +from typing import Optional, Literal, Callable, Iterable +from dataclasses import dataclass import matplotlib.colors as mcolors import matplotlib.patches as mpatches import matplotlib.pyplot as plt -from scipy.spatial import ConvexHull # for frontier plot import numpy as np import pandas as pd +from pandas import DataFrame, Series import seaborn as sns from adjustText import adjust_text # For the CEA plane figure to avoid overlaps in data labels @@ -56,105 +55,86 @@ if not os.path.exists(main_figurespath): os.makedirs(main_figurespath) -# Load result files -# ------------------------------------------------------------------------------------------------------------------ -results_folder = get_scenario_outputs('htm_and_hss_runs-2025-09-16T141811Z.py', outputfilepath)[0] - -# Check can read results from draw=0, run=0 -log = load_pickled_dataframes(results_folder, 0, 0) # look at one log (so can decide what to extract) -params = extract_params(results_folder) -info = get_scenario_info(results_folder) - -# Declare default parameters for cost analysis -# ------------------------------------------------------------------------------------------------------------------ -# Population scaling factor for malaria scale-up projections -population_scaling_factor = log['tlo.methods.demography']['scaling_factor']['scaling_factor'].iloc[0] -# Load the list of districts and their IDs -district_dict = pd.read_csv(resourcefilepath / 'demography' / 'ResourceFile_Population_2010.csv')[ - ['District_Num', 'District']].drop_duplicates() -district_dict = dict(zip(district_dict['District_Num'], district_dict['District'])) +# % +# Define necessary functions +# Convert results to dictionary to write text extracts for manuscript +def convert_results_to_dict(_df, _metric: str = 'median'): + draws = _df.index.to_list() + values = { + draw: { + chosen_metric: _df.loc[_df.index.get_level_values('draw') == draw, _metric].iloc[0], + "lower": _df.loc[_df.index.get_level_values('draw') == draw, 'lower'].iloc[0], + "upper": _df.loc[_df.index.get_level_values('draw') == draw, 'upper'].iloc[0] + } + for draw in draws + } + return values + +# Function to get incremental values +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() -# Period relevant for costing -TARGET_PERIOD = (Date(2025, 1, 1), Date(2035, 12, 31)) # This is the period that is costed -relevant_period_for_costing = [i.year for i in TARGET_PERIOD] -list_of_relevant_years_for_costing = list(range(relevant_period_for_costing[0], relevant_period_for_costing[1] + 1)) -# Choose central metric used - mean or median -chosen_metric = 'median' +def calculate_detailed_incremental_costs(_df, comparison_draw=0): + """ + Calculate incremental costs relative to a specified comparison draw. + The df needs to have the columns 'draw' and 'cost' and this will do a 'run' by run comparison for each disaggregated + cost component -# Scenarios -# Full list of scenarios used in the manuscript -main_manuscript_scenarios = {0: "Baseline", - 1: "Pessimistic HRH Scale-up", 2: "Historical HRH Scale-up", - 3: "Optimistic HRH Scale-up", - 4: "Consumables Increased to 75th Percentile", - 5: "Consumables Increased to HIV levels", 6: "Consumables Increased to EPI Levels", - 7: "HSS Expansion Package", - 8: "HIV Program Scale-up Without HSS Expansion", - 15: "HIV Program Scale-up With HSS Expansion Package", - 16: "TB Program Scale-up Without HSS Expansion", - 23: "TB Program Scale-up With HSS Expansion Package", - 24: "Malaria Program Scale-up Without HSS Expansion", - 31: "Malaria Program Scale-up With HSS Expansion Package", - 32: "HTM Programs Scale-up Without HSS Expansion", - 39: "HTM Programs Scale-up With HSS Expansion Package"} + Parameters: + - _df (pd.DataFrame): The input DataFrame with cost data. + - comparison_draw (int): The draw to use as the baseline for comparison. -frontier_scenarios = {9:"HIV Program Scale-up With Pessimistic HRH Scale-up", - 10:"HIV Program Scale-up With Historical HRH Scale-up", - 11:"HIV Program Scale-up With Optimistic HRH Scale-up", - 12:"HIV Program Scale-up With Consumables Increased to 75th Percentile", - 13:"HIV Program Scale-up With Consumables Increased to HIV levels", - 14:"HIV Program Scale-up With Consumables Increased to EPI Levels", - 17:"TB Program Scale-up With Pessimistic HRH Scale-up", - 18:"TB Program Scale-up With Historical HRH Scale-up", - 19:"TB Program Scale-up With Optimistic HRH Scale-up", - 20:"TB Program Scale-up With Consumables Increased to 75th Percentile", - 21:"TB Program Scale-up With Consumables Increased to HIV levels", - 22:"TB Program Scale-up With Consumables Increased to EPI Levels", - 25:"Malaria Program Scale-up With Pessimistic HRH Scale-up", - 26:"Malaria Program Scale-up With Historical HRH Scale-up", - 27:"Malaria Program Scale-up With Optimistic HRH Scale-up", - 28:"Malaria Program Scale-up With Consumables Increased to 75th Percentile", - 29:"Malaria Program Scale-up With Consumables Increased to HIV levels", - 30:"Malaria Program Scale-up With Consumables Increased to EPI Levels", - 33:"HTM Scale-up With Pessimistic HRH Scale-up", - 34:"HTM Scale-up With Historical HRH Scale-up", - 35:"HTM Scale-up With Optimistic HRH Scale-up", - 36:"HTM Scale-up With Consumables Increased to 75th Percentile", - 37:"HTM Scale-up With Consumables Increased to HIV levels", - 38:"HTM Scale-up With Consumables Increased to EPI Levels", - 40:"HIV + TB Scale-up Without HSS Expansion", - 41:"HIV + TB Scale-up With Historical HRH Scale-up", - 42:"HIV + TB Scale-up With Consumables Increased to 75th Percentile", - 43:"HIV + TB Scale-up With HSS Expansion", - 44:"HIV + Malaria Scale-up Without HSS Expansion", - 45:"HIV + Malaria Scale-up With Historical HRH Scale-up", - 46:"HIV + Malaria Scale-up With Consumables Increased to 75th Percentile", - 47:"HIV + Malaria Scale-up With HSS Expansion", - 48:"TB + Malaria Scale-up Without HSS Expansion", - 49:"TB + Malaria Scale-up With Historical HRH Scale-up", - 50:"TB + Malaria Scale-up With Consumables Increased to 75th Percentile", - 51:"TB + Malaria Scale-up With HSS Expansion"} + Returns: + - pd.DataFrame: DataFrame with an additional column 'incremental_cost'. + """ + # Identify all grouping columns (everything except 'draw' and 'cost') + group_cols = [col for col in _df.columns if col not in ['draw', 'cost']] -all_manuscript_scenarios = {**main_manuscript_scenarios, **frontier_scenarios} + # Extract baseline costs for the specified comparison draw + baseline_costs = _df[_df['draw'] == comparison_draw][group_cols + ['cost']] + baseline_costs = baseline_costs.rename(columns={'cost': 'baseline_cost'}) -all_manuscript_scenarios_reverse = {v: k for k, v in all_manuscript_scenarios.items()} + # Merge baseline costs with the original dataframe + merged_df = _df.merge(baseline_costs, on=group_cols, how='left') -baseline = all_manuscript_scenarios_reverse.get("Baseline") -horizontal_hss = all_manuscript_scenarios_reverse.get("HSS Expansion Package") -vertical_hiv = all_manuscript_scenarios_reverse.get("HIV Program Scale-up Without HSS Expansion") -vertical_tb = all_manuscript_scenarios_reverse.get("TB Program Scale-up Without HSS Expansion") -vertical_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up Without HSS Expansion") -vertical_htm = all_manuscript_scenarios_reverse.get("HTM Programs Scale-up Without HSS Expansion") -diagonal_hiv = all_manuscript_scenarios_reverse.get("HIV Program Scale-up With HSS Expansion Package") -diagonal_tb = all_manuscript_scenarios_reverse.get("TB Program Scale-up With HSS Expansion Package") -diagonal_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up With HSS Expansion Package") -diagonal_htm = all_manuscript_scenarios_reverse.get("HTM Programs Scale-up With HSS Expansion Package") + # Calculate incremental costs + merged_df['incremental_cost'] = merged_df['cost'] - merged_df['baseline_cost'] + merged_df.drop(columns=['baseline_cost', 'cost'], inplace=True) # Drop the baseline cost column + return merged_df[merged_df['draw'] != comparison_draw] -# Use letters instead of full scenario name for figures -all_manuscript_scenarios_substitutedict = {0: "0", 1: "A", 2: "B", 3: "C", 4: "D", 5: "E", 6: "F", 7: "G", 8: "H", - 15: "I", 16: "J", 23: "K", 24: "L", 31: "M", 32: "N", 39: "O"} +def summarise_detailed_costs_df(_df:pd.DataFrame, _metric:str = 'median', column_to_summarise='cost'): + # Summarize values + agg_func = np.mean if _metric == 'mean' else np.median + groupby_cols = [col for col in _df.columns if col not in ['run', column_to_summarise]] + _df = pd.concat( + { + chosen_metric: _df.groupby(by=groupby_cols, sort=False)[column_to_summarise].agg(agg_func), + 'lower': _df.groupby(by=groupby_cols, sort=False)[column_to_summarise].quantile(0.025), + 'upper': _df.groupby(by=groupby_cols, sort=False)[column_to_summarise].quantile(0.975), + }, + axis=1 + ) + summarised_df = pd.melt( + _df.reset_index(), + id_vars=groupby_cols, # Columns to keep + value_vars=[_metric, 'lower', 'upper'], # Columns to unpivot + var_name='stat', # New column name for the 'sub-category' of cost + value_name=column_to_summarise + ) + return summarised_df # Function to adjust color brightness (lighten/darken) def adjust_color(hex_color, factor=0.5): @@ -162,52 +142,8 @@ def adjust_color(hex_color, factor=0.5): adjusted_rgb = [(1 - factor) * c + factor * 1 for c in rgb] # Lighten return mcolors.to_hex(adjusted_rgb) - -# Generate color map -color_map = { - # Baseline (single color) - 'Baseline': 'black', - # HR scenarios - "Pessimistic HRH Scale-up": adjust_color('#9e0142', 0.5), - "Historical HRH Scale-up": adjust_color('#9e0142', 0.5), - "Optimistic HRH Scale-up": adjust_color('#9e0142', 0.5), - # Consumables scenarios - "Consumables Increased to 75th Percentile": adjust_color('#9e0142', 0.5), - "Consumables Increased to HIV levels": adjust_color('#9e0142', 0.5), - "Consumables Increased to EPI Levels": adjust_color('#9e0142', 0.5), - # HIV scenarios - "HIV Program Scale-up Without HSS Expansion": adjust_color('#fdae61', 0.5), - "HIV Program Scale-up With HSS Expansion Package": adjust_color('#66c2a5', 0.5), - # TB scenarios - "TB Program Scale-up Without HSS Expansion": adjust_color('#fdae61', 0.5), - "TB Program Scale-up With HSS Expansion Package": adjust_color('#66c2a5', 0.5), - # Malaria scenarios - "Malaria Program Scale-up Without HSS Expansion": adjust_color('#fdae61', 0.5), - "Malaria Program Scale-up With HSS Expansion Package": adjust_color('#66c2a5', 0.5), - # HSS scenarios - "HSS Expansion Package": adjust_color('#9e0142', 0.1), # Darker - # HTM scenarios - "HTM Programs Scale-up Without HSS Expansion": adjust_color('#fdae61', 0.1), - "HTM Programs Scale-up With HSS Expansion Package": adjust_color('#66c2a5', 0.1), -} - -# Cost-effectiveness threshold -chosen_cet = 191.4304166 # This is based on the estimate from Lomas et al (2023)- $160.595987085533 in 2019 USD coverted to 2023 USD -# based on Ochalek et al (2018) - the paper provided the value $61 in 2016 USD terms, this value is $77.4 in 2023 USD terms -chosen_value_of_statistical_life = 834 # This is based on Munthali et al (2020) National Planning Commission Report on -# "Medium and long-term impacts of a moderate lockdown (social restrictions) in response to the COVID-19 pandemic in Malawi" -chosen_value_of_statistical_life_upper = 2427.31 # upper bound estimated using Robinson et al method (income elasticity of VSL = 1) -chosen_value_of_statistical_life_lower = 425.96 # lower bound estimated using Robinson et al method (income elasticity of VSL = 1.5) -# lomas_consumption_value_of_health = 257.472 # this value is for 2025 (converted to 2023 USD) -# and assumed income elasticity of consumption value of health to be 1. - -# Above service level costs as a percentage of service level costs (This is used for the interpretation of ROI results) -# As per Opuni et al (2023), service level costs were 42% of the total cost --> above service level costs were 58% of total or 138% of sevice level costs -above_service_level_cost_proportion = 1.38 - - # Define a function to create bar plots -def do_standard_bar_plot_with_ci(_df, set_colors=None, annotations=None, +def do_standard_bar_plot_with_ci(_df: pd.DataFrame, set_colors=None, annotations=None, xticklabels_horizontal_and_wrapped=False, put_labels_in_legend=True, offset=1e6): @@ -297,7 +233,7 @@ def do_standard_bar_plot_with_ci(_df, set_colors=None, annotations=None, return fig, ax -# Thi sis a helper function for the upper_convex_frontier_from_dict +# This is a helper function for the upper_convex_frontier_from_dict def prune_inward_tail(frontier_cd, frontier_draws, tol=1e-9): """Remove any final points that violate the upper-hull orientation or reduce DALYs.""" keep = [] @@ -331,7 +267,7 @@ def upper_convex_frontier_from_dict(points_by_draw, include_collinear=False, tol # collapse exact cost ties -> keep the one with max effect # The monotone-chain hull expects strictly increasing x (cost). # If two rows have (nearly) the same cost, keep only the one with greater DALYs. - dedup = [] + dedup = [] # Deduplicated (cost) pairs of cost and DALYs for t in items: if dedup and abs(t[1] - dedup[-1][1]) <= tol: if t[2] > dedup[-1][2] + tol: @@ -355,7 +291,7 @@ def cross(o, a, b): o = hull[-2]; a = hull[-1] # These the the two previous points on the hull z = cross((o[0], o[1]), (a[0], a[1]), (c, v)) # For the UPPER hull we keep clockwise turns (z < 0). - # If z > 0 (concave) or (collinear and we don't include collinear), pop. + # If z > 0 (concave) or (we don't include collinear and it is collinear), pop. if z > tol or (not include_collinear and abs(z) <= tol): hull.pop() else: @@ -375,6 +311,441 @@ def incremental_icers(frontier_cd): dv = np.diff(frontier_cd[:, 1]) return dc, dv, dc / dv +# Create a function to generate threshold (or maximum) implementation costs which the scneario with the higher ROI +# can incur over and above those incurred by the scenario with the lower ROI while still delivering a higher ROI +def compute_breakeven_with_ci( + roi_df: pd.DataFrame, + monetary_benefit_df: pd.DataFrame, + incremental_cost_df: pd.DataFrame, + implementation_cost_df: pd.DataFrame, + draw_target: int, + draw_compare: int +): + """ + Compute breakeven implementation costs across runs and return median and 95% CI, i.e. how much additional cost + the target draw can incur until its ROI drops below the comparator draw + + Parameters: + ----------- + roi_df : pd.DataFrame + ROI DataFrame with index as (draw) and columns as runs. + monetary_benefit_df : pd.DataFrame + Monetised health benefits with same structure as roi_df. + incremental_cost_df : pd.DataFrame + Incremental input cost with same structure. + + draw_target : int + The draw index for the target scenario. + + draw_compare : int + The draw index for the scenario used to compute the threshold (e.g., vertical or HSS). + + Returns: + -------- + str + A label string in the format "$XXXM [$XX–$XXM]" + """ + + values = [] + for run in roi_df.columns: + roi_val = roi_df.loc[draw_compare, run] + health_benefit = monetary_benefit_df.loc[draw_target, run] + cost = incremental_cost_df.loc[draw_target, run] + + # Breakeven formula + val = (health_benefit - cost * (roi_val + 1)) / (roi_val + 1) + val = val - implementation_cost_df.loc[draw_compare, run] + values.append(val) + + values = np.array(values) / 1e6 # Convert to millions + median, lower, upper = np.median(values), np.percentile(values, 2.5), np.percentile(values, 97.5) + + # estimate median point on the ROI curves for the figure + roi_summary = convert_results_to_dict(summarize_cost_data(roi_df, chosen_metric), _metric = chosen_metric) + health_summary = convert_results_to_dict(summarize_cost_data(monetary_benefit_df, chosen_metric), _metric = chosen_metric) + cost_summary = convert_results_to_dict(summarize_cost_data(incremental_cost_df, chosen_metric), _metric = chosen_metric) + breakeven_central_value = (health_summary[draw_target][chosen_metric] - cost_summary[draw_target][ + chosen_metric] * (roi_summary[draw_compare][chosen_metric] + 1)) / ( + roi_summary[draw_compare][chosen_metric] + 1) + breakeven_central_value = breakeven_central_value / 1e6 + + return median, lower, upper, breakeven_central_value + +def do_incremental_cost_and_health_plot(incremental_cost_df: pd.DataFrame, + incremental_dalys_df: pd.DataFrame, + figname: str, + draws_with_icer_labels: Optional[list[int]] = None, + horizontal_scenarios: Optional[list[int]] = None, + vertical_scenarios: Optional[list[int]] = None, + diagonal_scenarios: Optional[list[int]] = None, + scenario_dict: Optional[dict[int, str]] = None, + # NEW frontier options + overlay_frontier: bool = True, + plot_icers_on_frontier: bool = False, + icer_lookup: Optional[dict[int, str]] = None, + # needs to be provided if plot_icers_on_frontier = True + frontier_color: str = "black", + frontier_linewidth: float = 2.0, + legend_position: Literal['inset', 'offset'] = "offset"): + + # Define colors for high level grouping of scenarios + color_map_scatter = {"HSS Investments": "#9e0142", "HTM-focussed Investments": "#fdae61", + "Combined Investments": "#66c2a5"} + horizontal_scenarios = horizontal_scenarios + vertical_scenarios = vertical_scenarios + diagonal_scenarios = diagonal_scenarios + # Draws to highlight, if any + darker_draws = draws_with_icer_labels + darker_set = set(darker_draws) + + # Reorder + new_order = horizontal_scenarios + vertical_scenarios + diagonal_scenarios + incremental_dalys_df = incremental_dalys_df.loc[new_order] + incremental_cost_df = incremental_cost_df.loc[new_order] + + # NUMERIC labels 1..N (based on the plotting order above) + scenario_labels = {draw: str(i + 1) for i, draw in enumerate(incremental_dalys_df.index)} + + # Extract stats for axes + x_median = incremental_dalys_df['median'] / 1e6 + x_lower = incremental_dalys_df['lower'] / 1e6 + x_upper = incremental_dalys_df['upper'] / 1e6 + y_median = incremental_cost_df['median'] + y_lower = incremental_cost_df['lower'] + y_upper = incremental_cost_df['upper'] + + fig, ax = plt.subplots(figsize=(16, 9)) + texts = [] + + # Scatter + error bars + numeric tiles + for draw in incremental_dalys_df.index: + if draw in horizontal_scenarios: + color = color_map_scatter["HSS Investments"] + elif draw in vertical_scenarios: + color = color_map_scatter["HTM-focussed Investments"] + elif draw in diagonal_scenarios: + color = color_map_scatter["Combined Investments"] + else: + color = "gray" + + is_darker = draw in darker_set + + # Plot the central values + ax.scatter(x_median[draw], y_median[draw], + s=36, color=color, + edgecolor='black' if is_darker else 'none', + linewidth=2.0 if is_darker else 0.0, zorder=3) + + # Plot the error bars + ax.errorbar(x_median[draw], y_median[draw], + xerr=[[x_median[draw] - x_lower[draw]], [x_upper[draw] - x_median[draw]]], + yerr=[[y_median[draw] - y_lower[draw]], [y_upper[draw] - y_median[draw]]], + fmt='none', color=color, alpha=0.6, capsize=4, elinewidth=1, zorder=2) + + # Add number labels above the dots + texts.append( + ax.text(x_median[draw], y_median[draw] + 1e7, + scenario_labels[draw], + fontsize=12, ha='center', va='bottom', color='white', weight='bold', + bbox=dict(facecolor=color, edgecolor='none', boxstyle='round,pad=0.3', alpha=0.9), + zorder=4) + ) + + # If this draw should show ICER, add it below the dot + if draws_with_icer_labels and draw in draws_with_icer_labels: + icer_text = icer_lookup.get(draw, "") + texts.append(ax.text( + x_median[draw] + 0.5, y_median[draw] + 0.1, + f"$ICER_{{{scenario_labels[draw]}}} =$\n{icer_text}", + fontsize=12, weight='bold', ha='center', va='top', color='black', zorder=4 + )) + + if texts: + adjust_text(texts, arrowprops=dict(arrowstyle='->')) + + # --- Frontier overlay (upper convex hull in (cost, dalys) vs baseline) --- + if overlay_frontier: + # Build dict in (cost, dalys) using *medians vs baseline* + points_by_draw = { + int(d): (float(y_median[d]), float(x_median[d] * 1e6)) # (cost, dalys) + for d in incremental_dalys_df.index + } + # add baseline origin because we're dealing with incremental costs + points_by_draw.setdefault(0, (0.0, 0.0)) + + frontier_draws, frontier_cd = upper_convex_frontier_from_dict(points_by_draw, include_collinear=False) + + # convert for plotting on scaled axes (x=DALYs millions, y=cost billions) + fx = frontier_cd[:, 1] / 1e6 + fy = frontier_cd[:, 0] + (frontier_line,) = ax.plot(fx, fy, '-', lw=frontier_linewidth, + color=frontier_color, label='Frontier (median)', zorder=3) + ax.scatter(fx, fy, s=24, facecolors='white', edgecolors=frontier_color, + linewidths=1.5, zorder=4) + + if plot_icers_on_frontier and len(frontier_cd) >= 2: + dc, dv, icers = incremental_icers(frontier_cd) + for i in range(len(icers)): + mx = (fx[i] + fx[i + 1]) / 2 + my = (fy[i] + fy[i + 1]) / 2 + 0.7 + ax.text(mx, my, f"ICER = ${icers[i]:.2f} \nper DALY averted", fontsize=7, ha='center', va='bottom', + color='black', + bbox=dict( + facecolor='white', # background color + edgecolor=frontier_color, + boxstyle='round,pad=0.2', + linewidth=0.8, + alpha=0.9 + ), + zorder=8) + + # --- Legend using numeric labels --- + legend_labels, dummy_handles = [], [] + scenario_categories = { + "HSS Investments": horizontal_scenarios, + "HTM-focussed Investments": vertical_scenarios, + "Combined Investments": diagonal_scenarios + } + for category, draws in scenario_categories.items(): + legend_labels.append(category) + dummy_handles.append(mpatches.Patch(color="white", label="")) + for draw in draws: + scenario_name = scenario_dict.get(draw, f"Scenario {draw}") + color = (color_map_scatter["HSS Investments"] if draw in horizontal_scenarios else + color_map_scatter["HTM-focussed Investments"] if draw in vertical_scenarios else + color_map_scatter["Combined Investments"] if draw in diagonal_scenarios else "gray") + is_darker = draw in darker_set + legend_labels.append(f"{scenario_labels[draw]}: {scenario_name}") + dummy_handles.append( + plt.Line2D([0], [0], marker='o', color='w', + markerfacecolor=color, markeredgecolor='black' if is_darker else 'none', + markersize=8, markeredgewidth=2 if is_darker else 0) + ) + + if (legend_position == "inset"): + bbox_to_anchor = (0, 1) + elif (legend_position == "offset"): + bbox_to_anchor = (1.1, 1) + else: + bbox_to_anchor = (0, 1) + + ax.legend(dummy_handles, legend_labels, bbox_to_anchor=bbox_to_anchor, loc='upper left', fontsize=9.2, + title="Scenario Key", frameon=True) + + ax.set_xlabel("DALYs Averted, millions", fontsize=12) + ax.set_ylabel("Incremental Scenario Cost, billions (USD)", fontsize=12) + ax.grid(True) + fig.tight_layout() + plt.savefig(figurespath / figname, dpi=300) + plt.close(fig) + +def get_manuscript_ready_table_of_projected_health_spending(_relevant_period_for_costing): + def get_total_population_by_year(_df): + years_needed = _relevant_period_for_costing # Malaria scale-up period years + _df['year'] = pd.to_datetime(_df['date']).dt.year + + # Validate that all necessary years are in the DataFrame + if not set(years_needed).issubset(_df['year'].unique()): + raise ValueError("Some years are not recorded in the dataset.") + + # Filter for relevant years and return the total population as a Series + return \ + _df.loc[_df['year'].between(min(years_needed), max(years_needed)), ['year', 'total']].set_index('year')[ + 'total'] + + # Get total population by year + total_population_by_year = extract_results( + results_folder, + module='tlo.methods.demography', + key='population', + custom_generate_series=get_total_population_by_year, + do_scaling=True + ).unstack().reset_index().rename(columns={0: 'population'}) + total_population_summary = total_population_by_year[total_population_by_year.draw == 0].groupby("year")[ + "population"].agg( + population="median" + ).reset_index() + unit_costs = load_unit_cost_assumptions(resourcefilepath) + health_spending_per_capita = unit_costs["health_spending_projections"] + health_spending_per_capita = health_spending_per_capita[health_spending_per_capita.year.isin( + list(range(_relevant_period_for_costing[0], _relevant_period_for_costing[1] + 1)))] + health_spending_per_capita = health_spending_per_capita[['year', 'total_mean']].apply( + pd.to_numeric, errors='coerce') + health_spending_per_capita_table = health_spending_per_capita.merge(total_population_summary, on="year", + how="left", validate="1:1") + health_spending_per_capita_table["total_health_spending"] = health_spending_per_capita_table['total_mean'] * \ + health_spending_per_capita_table['population'] + return health_spending_per_capita_table + + +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 = relevant_period_for_costing + assert set(_df.year.unique()).issuperset(years_needed), "Some years are not recorded." + _df = _df.loc[_df.year.between(*years_needed)].drop(columns=['date', 'sex', 'age_range']).groupby( + 'year').sum().sum(axis=1) + + # Initial year and discount rate + initial_year = min(_df.index.unique()) + + # Calculate the discounted values + discounted_values = _df / (1 + discount_rate_health) ** (_df.index - initial_year) + + return pd.Series(discounted_values.sum()) + +# The monetary value of the health benefit is delta health times CET (negative values are set to 0) +def get_monetary_value_of_incremental_health(_num_dalys_averted, _chosen_value_of_life_year): + monetary_value_of_incremental_health = (_num_dalys_averted * _chosen_value_of_life_year).clip(lower=0.0) + return monetary_value_of_incremental_health + +# %% CONFIG +# Load result files +# ------------------------------------------------------------------------------------------------------------------ +results_folder = get_scenario_outputs('htm_and_hss_runs-2025-09-16T141811Z.py', outputfilepath)[0] + +# Check can read results from draw=0, run=0 +log = load_pickled_dataframes(results_folder, 0, 0) # look at one log (so can decide what to extract) +params = extract_params(results_folder) +info = get_scenario_info(results_folder) + +# Declare default parameters for cost analysis +# ------------------------------------------------------------------------------------------------------------------ +# Population scaling factor for malaria scale-up projections +population_scaling_factor = log['tlo.methods.demography']['scaling_factor']['scaling_factor'].iloc[0] +# Load the list of districts and their IDs +district_dict = pd.read_csv(resourcefilepath / 'demography' / 'ResourceFile_Population_2010.csv')[ + ['District_Num', 'District']].drop_duplicates() +district_dict = dict(zip(district_dict['District_Num'], district_dict['District'])) + +# Period relevant for costing +TARGET_PERIOD = (Date(2025, 1, 1), Date(2035, 12, 31)) # This is the period that is costed +relevant_period_for_costing = [i.year for i in TARGET_PERIOD] +list_of_relevant_years_for_costing = list(range(relevant_period_for_costing[0], relevant_period_for_costing[1] + 1)) + +# Choose central metric used - mean or median +chosen_metric = 'median' + +# Scenarios +# Full list of scenarios used in the manuscript +main_manuscript_scenarios = {0: "Baseline", + 1: "Pessimistic HRH Scale-up", 2: "Historical HRH Scale-up", + 3: "Optimistic HRH Scale-up", + 4: "Consumables Increased to 75th Percentile", + 5: "Consumables Increased to HIV levels", 6: "Consumables Increased to EPI Levels", + 7: "HSS Expansion Package", + 8: "HIV Program Scale-up Without HSS Expansion", + 15: "HIV Program Scale-up With HSS Expansion Package", + 16: "TB Program Scale-up Without HSS Expansion", + 23: "TB Program Scale-up With HSS Expansion Package", + 24: "Malaria Program Scale-up Without HSS Expansion", + 31: "Malaria Program Scale-up With HSS Expansion Package", + 32: "HTM Programs Scale-up Without HSS Expansion", + 39: "HTM Programs Scale-up With HSS Expansion Package"} + +frontier_scenarios = {9:"HIV Program Scale-up With Pessimistic HRH Scale-up", + 10:"HIV Program Scale-up With Historical HRH Scale-up", + 11:"HIV Program Scale-up With Optimistic HRH Scale-up", + 12:"HIV Program Scale-up With Consumables Increased to 75th Percentile", + 13:"HIV Program Scale-up With Consumables Increased to HIV levels", + 14:"HIV Program Scale-up With Consumables Increased to EPI Levels", + 17:"TB Program Scale-up With Pessimistic HRH Scale-up", + 18:"TB Program Scale-up With Historical HRH Scale-up", + 19:"TB Program Scale-up With Optimistic HRH Scale-up", + 20:"TB Program Scale-up With Consumables Increased to 75th Percentile", + 21:"TB Program Scale-up With Consumables Increased to HIV levels", + 22:"TB Program Scale-up With Consumables Increased to EPI Levels", + 25:"Malaria Program Scale-up With Pessimistic HRH Scale-up", + 26:"Malaria Program Scale-up With Historical HRH Scale-up", + 27:"Malaria Program Scale-up With Optimistic HRH Scale-up", + 28:"Malaria Program Scale-up With Consumables Increased to 75th Percentile", + 29:"Malaria Program Scale-up With Consumables Increased to HIV levels", + 30:"Malaria Program Scale-up With Consumables Increased to EPI Levels", + 33:"HTM Scale-up With Pessimistic HRH Scale-up", + 34:"HTM Scale-up With Historical HRH Scale-up", + 35:"HTM Scale-up With Optimistic HRH Scale-up", + 36:"HTM Scale-up With Consumables Increased to 75th Percentile", + 37:"HTM Scale-up With Consumables Increased to HIV levels", + 38:"HTM Scale-up With Consumables Increased to EPI Levels", + 40:"HIV + TB Scale-up Without HSS Expansion", + 41:"HIV + TB Scale-up With Historical HRH Scale-up", + 42:"HIV + TB Scale-up With Consumables Increased to 75th Percentile", + 43:"HIV + TB Scale-up With HSS Expansion", + 44:"HIV + Malaria Scale-up Without HSS Expansion", + 45:"HIV + Malaria Scale-up With Historical HRH Scale-up", + 46:"HIV + Malaria Scale-up With Consumables Increased to 75th Percentile", + 47:"HIV + Malaria Scale-up With HSS Expansion", + 48:"TB + Malaria Scale-up Without HSS Expansion", + 49:"TB + Malaria Scale-up With Historical HRH Scale-up", + 50:"TB + Malaria Scale-up With Consumables Increased to 75th Percentile", + 51:"TB + Malaria Scale-up With HSS Expansion"} + +all_manuscript_scenarios = {**main_manuscript_scenarios, **frontier_scenarios} + +all_manuscript_scenarios_reverse = {v: k for k, v in all_manuscript_scenarios.items()} + + # Identify scenarios to highlight in the analysis +baseline = all_manuscript_scenarios_reverse.get("Baseline") +horizontal_hss = all_manuscript_scenarios_reverse.get("HSS Expansion Package") +vertical_hiv = all_manuscript_scenarios_reverse.get("HIV Program Scale-up Without HSS Expansion") +vertical_tb = all_manuscript_scenarios_reverse.get("TB Program Scale-up Without HSS Expansion") +vertical_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up Without HSS Expansion") +vertical_htm = all_manuscript_scenarios_reverse.get("HTM Programs Scale-up Without HSS Expansion") +diagonal_hiv = all_manuscript_scenarios_reverse.get("HIV Program Scale-up With HSS Expansion Package") +diagonal_tb = all_manuscript_scenarios_reverse.get("TB Program Scale-up With HSS Expansion Package") +diagonal_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up With HSS Expansion Package") +diagonal_htm = all_manuscript_scenarios_reverse.get("HTM Programs Scale-up With HSS Expansion Package") + +# Use letters instead of full scenario name for figures +all_manuscript_scenarios_substitutedict = {0: "0", 1: "A", 2: "B", 3: "C", 4: "D", 5: "E", 6: "F", 7: "G", 8: "H", + 15: "I", 16: "J", 23: "K", 24: "L", 31: "M", 32: "N", 39: "O"} + +# Generate color map +color_map = { + # Baseline (single color) + 'Baseline': 'black', + # HR scenarios + "Pessimistic HRH Scale-up": adjust_color('#9e0142', 0.5), + "Historical HRH Scale-up": adjust_color('#9e0142', 0.5), + "Optimistic HRH Scale-up": adjust_color('#9e0142', 0.5), + # Consumables scenarios + "Consumables Increased to 75th Percentile": adjust_color('#9e0142', 0.5), + "Consumables Increased to HIV levels": adjust_color('#9e0142', 0.5), + "Consumables Increased to EPI Levels": adjust_color('#9e0142', 0.5), + # HIV scenarios + "HIV Program Scale-up Without HSS Expansion": adjust_color('#fdae61', 0.5), + "HIV Program Scale-up With HSS Expansion Package": adjust_color('#66c2a5', 0.5), + # TB scenarios + "TB Program Scale-up Without HSS Expansion": adjust_color('#fdae61', 0.5), + "TB Program Scale-up With HSS Expansion Package": adjust_color('#66c2a5', 0.5), + # Malaria scenarios + "Malaria Program Scale-up Without HSS Expansion": adjust_color('#fdae61', 0.5), + "Malaria Program Scale-up With HSS Expansion Package": adjust_color('#66c2a5', 0.5), + # HSS scenarios + "HSS Expansion Package": adjust_color('#9e0142', 0.1), # Darker + # HTM scenarios + "HTM Programs Scale-up Without HSS Expansion": adjust_color('#fdae61', 0.1), + "HTM Programs Scale-up With HSS Expansion Package": adjust_color('#66c2a5', 0.1), +} +plt.rcParams['xtick.labelsize'] = 12 +plt.rcParams['ytick.labelsize'] = 12 + +# Cost-effectiveness threshold +chosen_cet = 191.4304166 # This is based on the estimate from Lomas et al (2023)- $160.595987085533 in 2019 USD coverted to 2023 USD +# based on Ochalek et al (2018) - the paper provided the value $61 in 2016 USD terms, this value is $77.4 in 2023 USD terms +chosen_value_of_statistical_life = 834 # This is based on Munthali et al (2020) National Planning Commission Report on +# "Medium and long-term impacts of a moderate lockdown (social restrictions) in response to the COVID-19 pandemic in Malawi" +chosen_value_of_statistical_life_upper = 2427.31 # upper bound estimated using Robinson et al method (income elasticity of VSL = 1) +chosen_value_of_statistical_life_lower = 425.96 # lower bound estimated using Robinson et al method (income elasticity of VSL = 1.5) +# lomas_consumption_value_of_health = 257.472 # this value is for 2025 (converted to 2023 USD) +# and assumed income elasticity of consumption value of health to be 1. + +# Above service level costs as a percentage of service level costs (This is used for the interpretation of ROI results) +# As per Opuni et al (2023), service level costs were 42% of the total cost --> above service level costs were 58% of total or 138% of sevice level costs +above_service_level_cost_proportion = 1.38 + # Define alternative discount rate sets as a list of dictionaries alternative_discount_rates = [ {"discount_rate_cost": 0.03, "discount_rate_health": 0, "discounting_scenario": 'WHO-CHOICE (0.03,0)'}, @@ -392,7 +763,6 @@ def incremental_icers(frontier_cd): print(f"NOW RUNNING SCENARIO {rates['discounting_scenario']}...") - # %% # Estimate standard input costs of scenario # ----------------------------------------------------------------------------------------------------------------------- input_costs = estimate_input_cost_of_scenarios(results_folder, resourcefilepath, @@ -617,56 +987,20 @@ def melt_and_label_malaria_scaleup_cost(_df, label): # 1. Calculate incremental cost # ----------------------------------------------------------------------------------------------------------------------- total_input_cost = input_costs.groupby(['draw', 'run'])['cost'].sum() - total_input_cost_summarized = summarize_cost_data(total_input_cost.unstack(level='run'), _metric=chosen_metric) - - - 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() - + #total_input_cost_summarized = summarize_cost_data(total_input_cost.unstack(level='run'), _metric=chosen_metric) incremental_scenario_cost = (pd.DataFrame( find_difference_relative_to_comparison( total_input_cost, comparison=0) # sets the comparator to draw 0 which is the Actual scenario ).T.iloc[0].unstack()).T + incremental_scenario_cost_subset_for_figure = incremental_scenario_cost[ incremental_scenario_cost.index.get_level_values('draw').isin(list(main_manuscript_scenarios.keys()))] - incremental_scenario_cost_summarised = summarize_cost_data(incremental_scenario_cost_subset_for_figure, - _metric=chosen_metric) # 2. Monetary value of health impact # ----------------------------------------------------------------------------------------------------------------------- - 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 = relevant_period_for_costing - assert set(_df.year.unique()).issuperset(years_needed), "Some years are not recorded." - _df = _df.loc[_df.year.between(*years_needed)].drop(columns=['date', 'sex', 'age_range']).groupby( - 'year').sum().sum(axis=1) - - # Initial year and discount rate - initial_year = min(_df.index.unique()) - - # Calculate the discounted values - discounted_values = _df / (1 + discount_rate_health) ** (_df.index - initial_year) - - return pd.Series(discounted_values.sum()) - - num_dalys = extract_results( results_folder, module='tlo.methods.healthburden', @@ -682,8 +1016,6 @@ def get_num_dalys(_df): num_dalys.loc[0], comparison=0) # sets the comparator to 0 which is the Actual scenario ).T.iloc[0].unstack(level='run')) - num_dalys_averted = num_dalys_averted[num_dalys_averted.index.get_level_values('draw').isin( - list(all_manuscript_scenarios.keys()))] # keep only relevant draws # Plot DALYs num_dalys_averted_subset_for_figure = num_dalys_averted[ @@ -698,7 +1030,7 @@ def get_num_dalys(_df): ], xticklabels_horizontal_and_wrapped=False, put_labels_in_legend=True, - offset=2, + offset=0.5, ) ax.set_title(name_of_plot) ax.set_ylabel('DALYs \n(Millions)') @@ -707,20 +1039,13 @@ def get_num_dalys(_df): fig.savefig(figurespath / name_of_plot.replace(' ', '_').replace(',', '')) plt.close(fig) - - # The monetary value of the health benefit is delta health times CET (negative values are set to 0) - def get_monetary_value_of_incremental_health(_num_dalys_averted, _chosen_value_of_life_year): - monetary_value_of_incremental_health = (_num_dalys_averted * _chosen_value_of_life_year).clip(lower=0.0) - return monetary_value_of_incremental_health - - # 3. Estimate and plot ICERs # ---------------------------------------------------- icers = incremental_scenario_cost.div(num_dalys_averted) # Element-wise division - # icers = icers.mask(num_dalys_averted < 0) + icers_summarized = summarize_cost_data(icers, _metric=chosen_metric) dominated_scenarios = icers_summarized['median'] < 0 - icers_summarized[dominated_scenarios] = np.nan # The dominanted scenarios are assigned as ICER of NaN + icers_summarized[dominated_scenarios] = np.nan # The dominated scenarios are assigned as ICER of NaN icers_summarized_subset_for_figure = icers_summarized[ icers_summarized.index.get_level_values('draw').isin(list(main_manuscript_scenarios.keys()))] @@ -752,9 +1077,7 @@ def get_monetary_value_of_incremental_health(_num_dalys_averted, _chosen_value_o plt.close(fig) # Extract ICERs into a table for manuscript - icers_summarized_subset_for_table = icers_summarized[ - icers_summarized.index.get_level_values('draw').isin(list(all_manuscript_scenarios.keys()))] - icers_summarized_subset_for_table = icers_summarized_subset_for_table.reset_index() + icers_summarized_subset_for_table = icers_summarized_subset_for_figure.reset_index() icers_summarized_subset_for_table['scenario'] = icers_summarized_subset_for_table['draw'].map( all_manuscript_scenarios) icers_summarized_subset_for_table['ICER (2023 USD)'] = icers_summarized_subset_for_table.apply( @@ -777,176 +1100,7 @@ def get_monetary_value_of_incremental_health(_num_dalys_averted, _chosen_value_o # Plot incremental health and cost in a scatterplot - # ---------------------------------------------------- - def do_incremental_cost_and_health_plot(incremental_cost_df, - incremental_dalys_df, - figname, - draws_with_icer_labels: Optional[list[int]] = None, - horizontal_scenarios: Optional[list[int]] = None, - vertical_scenarios: Optional[list[int]] = None, - diagonal_scenarios: Optional[list[int]] = None, - scenario_dict: Optional[dict[int, str]] = None, - # NEW frontier options - overlay_frontier: bool = True, - plot_icers_on_frontier: bool = False, - icer_lookup: Optional[dict[int, str]] = None, - # needs to be provided if plot_icers_on_frontier = True - frontier_color: str = "black", - frontier_linewidth: float = 2.0, - legend_position: Literal['inset', 'offset'] = "offset"): - - # Define colors for high level grouping of scenarios - color_map_scatter = {"HSS Investments": "#9e0142", "HTM-focussed Investments": "#fdae61", - "Combined Investments": "#66c2a5"} - horizontal_scenarios = horizontal_scenarios - vertical_scenarios = vertical_scenarios - diagonal_scenarios = diagonal_scenarios - # Draws to highlight, if any - darker_draws = draws_with_icer_labels - darker_set = set(darker_draws) - - # Reorder - new_order = horizontal_scenarios + vertical_scenarios + diagonal_scenarios - incremental_dalys_df = incremental_dalys_df.loc[new_order] - incremental_cost_df = incremental_cost_df.loc[new_order] - - # NUMERIC labels 1..N (based on the plotting order above) - scenario_labels = {draw: str(i + 1) for i, draw in enumerate(incremental_dalys_df.index)} - - # Extract stats for axes - x_median = incremental_dalys_df['median'] / 1e6 - x_lower = incremental_dalys_df['lower'] / 1e6 - x_upper = incremental_dalys_df['upper'] / 1e6 - y_median = incremental_cost_df['median'] - y_lower = incremental_cost_df['lower'] - y_upper = incremental_cost_df['upper'] - - fig, ax = plt.subplots(figsize=(16, 9)) - texts = [] - - # Scatter + error bars + numeric tiles - for draw in incremental_dalys_df.index: - if draw in horizontal_scenarios: - color = color_map_scatter["HSS Investments"] - elif draw in vertical_scenarios: - color = color_map_scatter["HTM-focussed Investments"] - elif draw in diagonal_scenarios: - color = color_map_scatter["Combined Investments"] - else: - color = "gray" - - is_darker = draw in darker_set - - # Plot the central values - ax.scatter(x_median[draw], y_median[draw], - s=36, color=color, - edgecolor='black' if is_darker else 'none', - linewidth=2.0 if is_darker else 0.0, zorder=3) - - # Plot the error bars - ax.errorbar(x_median[draw], y_median[draw], - xerr=[[x_median[draw] - x_lower[draw]], [x_upper[draw] - x_median[draw]]], - yerr=[[y_median[draw] - y_lower[draw]], [y_upper[draw] - y_median[draw]]], - fmt='none', color=color, alpha=0.6, capsize=4, elinewidth=1, zorder=2) - - # Add number labels above the dots - texts.append( - ax.text(x_median[draw], y_median[draw] + 1e7, - scenario_labels[draw], - fontsize=12, ha='center', va='bottom', color='white', weight='bold', - bbox=dict(facecolor=color, edgecolor='none', boxstyle='round,pad=0.3', alpha=0.9), - zorder=4) - ) - - # If this draw should show ICER, add it below the dot - if draws_with_icer_labels and draw in draws_with_icer_labels: - icer_text = icer_lookup.get(draw, "") - texts.append(ax.text( - x_median[draw] + 0.5, y_median[draw] + 0.1, - f"$ICER_{{{scenario_labels[draw]}}} =$\n{icer_text}", - fontsize=12, weight='bold', ha='center', va='top', color='black', zorder=4 - )) - - if texts: - adjust_text(texts, arrowprops=dict(arrowstyle='->')) - - # --- Frontier overlay (upper convex hull in (cost, dalys) vs baseline) --- - if overlay_frontier: - # Build dict in (cost, dalys) using *medians vs baseline* - points_by_draw = { - int(d): (float(y_median[d]), float(x_median[d] * 1e6)) # (cost, dalys) - for d in incremental_dalys_df.index - } - # add baseline origin because we're dealing with incremental costs - points_by_draw.setdefault(0, (0.0, 0.0)) - - frontier_draws, frontier_cd = upper_convex_frontier_from_dict(points_by_draw, include_collinear=False) - - # convert for plotting on scaled axes (x=DALYs millions, y=cost billions) - fx = frontier_cd[:, 1] / 1e6 - fy = frontier_cd[:, 0] - (frontier_line,) = ax.plot(fx, fy, '-', lw=frontier_linewidth, - color=frontier_color, label='Frontier (median)', zorder=5) - ax.scatter(fx, fy, s=24, facecolors='white', edgecolors=frontier_color, - linewidths=1.5, zorder=6) - - if plot_icers_on_frontier and len(frontier_cd) >= 2: - dc, dv, icers = incremental_icers(frontier_cd) - for i in range(len(icers)): - mx = (fx[i] + fx[i + 1]) / 2 - my = (fy[i] + fy[i + 1]) / 2 + 0.7 - ax.text(mx, my, f"ICER = ${icers[i]:.2f} \nper DALY averted", fontsize=7, ha='center', va='bottom', - color='black', - bbox=dict( - facecolor='white', # background color - edgecolor=frontier_color, - boxstyle='round,pad=0.2', - linewidth=0.8, - alpha=0.9 - ), - zorder=10) - - # --- Legend using numeric labels --- - legend_labels, dummy_handles = [], [] - scenario_categories = { - "HSS Investments": horizontal_scenarios, - "HTM-focussed Investments": vertical_scenarios, - "Combined Investments": diagonal_scenarios - } - for category, draws in scenario_categories.items(): - legend_labels.append(category) - dummy_handles.append(mpatches.Patch(color="white", label="")) - for draw in draws: - scenario_name = scenario_dict.get(draw, f"Scenario {draw}") - color = (color_map_scatter["HSS Investments"] if draw in horizontal_scenarios else - color_map_scatter["HTM-focussed Investments"] if draw in vertical_scenarios else - color_map_scatter["Combined Investments"] if draw in diagonal_scenarios else "gray") - is_darker = draw in darker_set - legend_labels.append(f"{scenario_labels[draw]}: {scenario_name}") - dummy_handles.append( - plt.Line2D([0], [0], marker='o', color='w', - markerfacecolor=color, markeredgecolor='black' if is_darker else 'none', - markersize=8, markeredgewidth=2 if is_darker else 0) - ) - - if (legend_position == "inset"): - bbox_to_anchor = (0, 1) - elif (legend_position == "offset"): - bbox_to_anchor = (1.1, 1) - else: - bbox_to_anchor = (0, 1) - - ax.legend(dummy_handles, legend_labels, bbox_to_anchor=bbox_to_anchor, loc='upper left', fontsize=9.2, - title="Scenario Key", frameon=True) - - ax.set_xlabel("DALYs Averted, millions") - ax.set_ylabel("Incremental Scenario Cost, billions (USD)") - ax.grid(True) - fig.tight_layout() - plt.savefig(figurespath / figname, dpi=300) - plt.close(fig) - - + # ----------------------------------------------------= do_incremental_cost_and_health_plot(incremental_cost_df=summarize_cost_data(incremental_scenario_cost, chosen_metric), incremental_dalys_df=summarize_cost_data(num_dalys_averted, chosen_metric), draws_with_icer_labels=[1], @@ -990,71 +1144,12 @@ def do_incremental_cost_and_health_plot(incremental_cost_df, # Extract projected health spending table for appendix - def get_manuscript_ready_table_of_projected_health_spending(_relevant_period_for_costing): - def get_total_population_by_year(_df): - years_needed = _relevant_period_for_costing # Malaria scale-up period years - _df['year'] = pd.to_datetime(_df['date']).dt.year - - # Validate that all necessary years are in the DataFrame - if not set(years_needed).issubset(_df['year'].unique()): - raise ValueError("Some years are not recorded in the dataset.") - - # Filter for relevant years and return the total population as a Series - return \ - _df.loc[_df['year'].between(min(years_needed), max(years_needed)), ['year', 'total']].set_index('year')[ - 'total'] - - # Get total population by year - total_population_by_year = extract_results( - results_folder, - module='tlo.methods.demography', - key='population', - custom_generate_series=get_total_population_by_year, - do_scaling=True - ).unstack().reset_index().rename(columns={0: 'population'}) - total_population_summary = total_population_by_year[total_population_by_year.draw == 0].groupby("year")[ - "population"].agg( - population="median" - ).reset_index() - unit_costs = load_unit_cost_assumptions(resourcefilepath) - health_spending_per_capita = unit_costs["health_spending_projections"] - health_spending_per_capita = health_spending_per_capita[health_spending_per_capita.year.isin( - list(range(_relevant_period_for_costing[0], _relevant_period_for_costing[1] + 1)))] - health_spending_per_capita = health_spending_per_capita[['year', 'total_mean']].apply( - pd.to_numeric, errors='coerce') - health_spending_per_capita_table = health_spending_per_capita.merge(total_population_summary, on="year", - how="left", validate="1:1") - health_spending_per_capita_table["total_health_spending"] = health_spending_per_capita_table['total_mean'] * \ - health_spending_per_capita_table['population'] - return health_spending_per_capita_table - - health_spending_per_capita_table = get_manuscript_ready_table_of_projected_health_spending( _relevant_period_for_costing=relevant_period_for_costing) health_spending_per_capita_table.to_csv(figurespath / 'projected_health_spending.csv', index=False) # Combined ROI plot of relevant scenarios # ROI plot comparing HSS alone, HTM without HSS, and HTM with HSS - if rates["discounting_scenario"] == 'MAIN (0.03,0.03)': - legend_switch_for_main_roi_plot = False # Don't show legend for main plot - else: - legend_switch_for_main_roi_plot = True - - - # Convert results to dictionary to write text extracts for manuscript - def convert_results_to_dict(_df): - draws = _df.index.to_list() - values = { - draw: { - chosen_metric: _df.loc[_df.index.get_level_values('draw') == draw, chosen_metric].iloc[0], - "lower": _df.loc[_df.index.get_level_values('draw') == draw, 'lower'].iloc[0], - "upper": _df.loc[_df.index.get_level_values('draw') == draw, 'upper'].iloc[0] - } - for draw in draws - } - return values - - # Do ROI plots for different VSLY values i = 0 vsly_fig_suffixes = ['LOWER', 'UPPER', 'MAIN'] @@ -1066,13 +1161,7 @@ def convert_results_to_dict(_df): roi_at_0_implementation_cost = benefit_at_0_implementation_cost.div(abs(incremental_scenario_cost)) roi_at_0_implementation_cost_summarized = summarize_cost_data(roi_at_0_implementation_cost, _metric=chosen_metric) - roi_at_0_implementation_cost_dict = convert_results_to_dict(roi_at_0_implementation_cost_summarized) - - # health_benefit_summarised = convert_results_to_dict(summarize_cost_data( - # get_monetary_value_of_incremental_health(num_dalys_averted, vsly), - # _metric=chosen_metric)) - # incremental_scenario_cost_summarised = convert_results_to_dict( - # summarize_cost_data(incremental_scenario_cost, _metric=chosen_metric)) + roi_at_0_implementation_cost_dict = convert_results_to_dict(roi_at_0_implementation_cost_summarized, _metric = chosen_metric) # ROI at implementation costs = 138% of input costs implementation_cost_upper_limit = incremental_scenario_cost * above_service_level_cost_proportion @@ -1083,75 +1172,15 @@ def convert_results_to_dict(_df): vsly) - incremental_scenario_cost - implementation_cost_upper_limit) - roi_at_upper_limit_implementation_cost = benefit_at_0_implementation_cost.div( + roi_at_upper_limit_implementation_cost = benefit_at_upper_limit_implementation_cost.div( abs(incremental_scenario_cost + implementation_cost_upper_limit)) roi_at_upper_limit_implementation_cost_summarized = summarize_cost_data(roi_at_upper_limit_implementation_cost, _metric=chosen_metric) - roi_at_upper_limit_implementation_cost_dict = convert_results_to_dict(roi_at_upper_limit_implementation_cost_summarized) - - - # Create a function to generate threshold (or maximum) implementation costs which the scneario with the higher ROI - # can incur over and above those incurred by the scenario with the lower ROI while still delivering a higher ROI - def compute_breakeven_with_ci( - roi_df: pd.DataFrame, - monetary_benefit_df: pd.DataFrame, - incremental_cost_df: pd.DataFrame, - implementation_cost_df: pd.DataFrame, - draw_target: int, - draw_compare: int - ): - """ - Compute breakeven implementation costs across runs and return median and 95% CI. - - Parameters: - ----------- - roi_df : pd.DataFrame - ROI DataFrame with index as (draw) and columns as runs. - monetary_benefit_df : pd.DataFrame - Monetised health benefits with same structure as roi_df. - incremental_cost_df : pd.DataFrame - Incremental input cost with same structure. - - draw_target : int - The draw index for the target scenario. - - draw_compare : int - The draw index for the scenario used to compute the threshold (e.g., vertical or HSS). - - Returns: - -------- - str - A label string in the format "$XXXM [$XX–$XXM]" - """ - - values = [] - for run in roi_df.columns: - roi_val = roi_df.loc[draw_compare, run] - health_benefit = monetary_benefit_df.loc[draw_target, run] - cost = incremental_cost_df.loc[draw_target, run] - - # Breakeven formula - val = (health_benefit - cost * (roi_val + 1)) / (roi_val + 1) - val = val - implementation_cost_df.loc[draw_compare, run] - values.append(val) - - values = np.array(values) / 1e6 # Convert to millions - median, lower, upper = np.median(values), np.percentile(values, 2.5), np.percentile(values, 97.5) - - # estimate median point on the ROI curves for the figure - roi_summary = convert_results_to_dict(summarize_cost_data(roi_df, chosen_metric)) - health_summary = convert_results_to_dict(summarize_cost_data(monetary_benefit_df, chosen_metric)) - cost_summary = convert_results_to_dict(summarize_cost_data(incremental_cost_df, chosen_metric)) - breakeven_central_value = (health_summary[draw_target][chosen_metric] - cost_summary[draw_target][ - chosen_metric] * (roi_summary[draw_compare][chosen_metric] + 1)) / ( - roi_summary[draw_compare][chosen_metric] + 1) - breakeven_central_value = breakeven_central_value / 1e6 - - return median, lower, upper, breakeven_central_value - + roi_at_upper_limit_implementation_cost_dict = convert_results_to_dict(roi_at_upper_limit_implementation_cost_summarized, _metric = chosen_metric) # Breakeven implementation cost for joint HTM diagonal scenario (V joint HTM vertical scenario); + draw_colors = {horizontal_hss: '#9e0142', vertical_htm: '#fdae61', diagonal_htm: '#66c2a5'} # Implementation cost = 0 - breakeven_lower_v_htm = compute_breakeven_with_ci( + breakeven_lower_htmhss_v_htm = compute_breakeven_with_ci( roi_df=roi_at_0_implementation_cost, monetary_benefit_df=get_monetary_value_of_incremental_health(num_dalys_averted, vsly), incremental_cost_df=incremental_scenario_cost, @@ -1161,7 +1190,7 @@ def compute_breakeven_with_ci( ) # Breakeven implementation cost for joint HTM diagonal scenario (V joint HTM vertical scenario); # Implementation cost = 138% of Input costs - breakeven_upper_v_htm = compute_breakeven_with_ci( + breakeven_upper_htmhss_v_htm = compute_breakeven_with_ci( roi_df=roi_at_upper_limit_implementation_cost, monetary_benefit_df=get_monetary_value_of_incremental_health(num_dalys_averted, vsly), incremental_cost_df=incremental_scenario_cost, @@ -1169,27 +1198,28 @@ def compute_breakeven_with_ci( draw_target=diagonal_htm, draw_compare=vertical_htm ) - # Breakeven implementation cost for joint HTM diagonal scenario (V horizontal scenario); + # Breakeven implementation cost for horizontal scenario (V vertical HTM); # Implementation cost = 0 - breakeven_lower_v_hss = compute_breakeven_with_ci( + breakeven_lower_hss_v_htm = compute_breakeven_with_ci( roi_df=roi_at_0_implementation_cost, monetary_benefit_df=get_monetary_value_of_incremental_health(num_dalys_averted, vsly), incremental_cost_df=incremental_scenario_cost, implementation_cost_df=0 * incremental_scenario_cost, - draw_target=diagonal_htm, - draw_compare=horizontal_hss + draw_target=horizontal_hss, + draw_compare=vertical_htm ) # Breakeven implementation cost for joint HTM diagonal scenario (V horizontal scenario); # Implementation cost = 138% of Input costs - breakeven_upper_v_hss = compute_breakeven_with_ci( + breakeven_upper_hss_v_htm = compute_breakeven_with_ci( roi_df=roi_at_upper_limit_implementation_cost, monetary_benefit_df=get_monetary_value_of_incremental_health(num_dalys_averted, vsly), incremental_cost_df=incremental_scenario_cost, implementation_cost_df=implementation_cost_upper_limit, - draw_target=diagonal_htm, - draw_compare=horizontal_hss + draw_target=horizontal_hss, + draw_compare=vertical_htm ) + # For the MAIN plot generate a version of the ROI plots with breakeven implementation costs if (vsly_fig_suffixes[i] == 'MAIN') & (rates["discounting_scenario"] == 'MAIN (0.03,0.03)'): # The breakeven # ASC costs are only plotted in th main results plot # Generate data for the representation of breakeven implementation costs on the ROI plots @@ -1197,49 +1227,62 @@ def compute_breakeven_with_ci( { 'y_value': roi_at_0_implementation_cost_dict[vertical_htm][chosen_metric], 'x_start': 0, - 'x_end': breakeven_lower_v_htm[0], + 'x_end': breakeven_lower_htmhss_v_htm[0], # where the horizontal line intersects the ROI curve of the diagonal strategy - 'label': f"a = ${breakeven_lower_v_htm[0]:.0f}M [${breakeven_lower_v_htm[1]:.0f}M–${breakeven_lower_v_htm[2]:.0f}M]", - 'color': '#fdae61', + 'label': f"a = ${breakeven_lower_htmhss_v_htm[0]:.0f}M [${breakeven_lower_htmhss_v_htm[1]:.0f}M–${breakeven_lower_htmhss_v_htm[2]:.0f}M]", + 'color': '#66c2a5', 'scenario_label': 'Diagonal versus Vertical (at ASC = 0)' }, { 'y_value': roi_at_upper_limit_implementation_cost_dict[vertical_htm][chosen_metric], 'x_start': implementation_cost_upper_limit_dict[vertical_htm][chosen_metric] / 1e6, - 'x_end': breakeven_upper_v_htm[0], + 'x_end': breakeven_upper_htmhss_v_htm[0], # where the horizontal line intersects the ROI curve of the diagonal strategy - 'label': f"b = ${breakeven_upper_v_htm[0]:.0f}M [${breakeven_upper_v_htm[1]:.0f}M–${breakeven_upper_v_htm[2]:.0f}M]", - 'color': '#fdae61', + 'label': f"b = ${breakeven_upper_htmhss_v_htm[0]:.0f}M [${breakeven_upper_htmhss_v_htm[1]:.0f}M–${breakeven_upper_htmhss_v_htm[2]:.0f}M]", + 'color': '#66c2a5', 'scenario_label': 'Diagonal versus Vertical (at ASC = 138% X SC)', 'y_label_offset': 0.25 }, { - 'y_value': roi_at_0_implementation_cost_dict[horizontal_hss][chosen_metric], + 'y_value': roi_at_0_implementation_cost_dict[vertical_htm][chosen_metric], 'x_start': 0, - 'x_end': breakeven_lower_v_hss[3], + 'x_end': breakeven_lower_hss_v_htm[3], # where the horizontal line intersects the ROI curve of the diagonal strategy - 'label': f"c = ${breakeven_lower_v_hss[0]:.0f}M [${breakeven_lower_v_hss[1]:.0f}M–${breakeven_lower_v_hss[2]:.0f}M]", + 'label': f"c = ${breakeven_lower_hss_v_htm[0]:.0f}M [${breakeven_lower_hss_v_htm[1]:.0f}M–${breakeven_lower_hss_v_htm[2]:.0f}M]", 'color': '#9e0142', 'scenario_label': 'Diagonal versus Horizontal (at ASC = 0)', 'y_label_offset': -0.25 }, { - 'y_value': roi_at_upper_limit_implementation_cost_dict[horizontal_hss][chosen_metric], - 'x_start': implementation_cost_upper_limit_dict[horizontal_hss][chosen_metric] / 1e6, - 'x_end': breakeven_upper_v_hss[3], + 'y_value': roi_at_upper_limit_implementation_cost_dict[vertical_htm][chosen_metric], + 'x_start': implementation_cost_upper_limit_dict[vertical_htm][chosen_metric] / 1e6, + 'x_end': breakeven_upper_hss_v_htm[3], # where the horizontal line intersects the ROI curve of the diagonal strategy - 'label': f"d = ${breakeven_upper_v_hss[0]:.0f}M [${breakeven_upper_v_hss[1]:.0f}M–${breakeven_upper_v_hss[2]:.0f}M]", + 'label': f"d = ${breakeven_upper_hss_v_htm[0]:.0f}M [${breakeven_upper_hss_v_htm[1]:.0f}M–${breakeven_upper_hss_v_htm[2]:.0f}M]", 'color': '#9e0142', 'scenario_label': 'Diagonal versus Horizontal (at ASC = 138% X SC)', 'y_label_offset': -0.25 } ] - legend_switch_for_main_roi_plot = False + + generate_multiple_scenarios_roi_plot( + _monetary_value_of_incremental_health=get_monetary_value_of_incremental_health(num_dalys_averted, + _chosen_value_of_life_year=vsly), + _incremental_input_cost=incremental_scenario_cost, + _draws=[horizontal_hss, vertical_htm, diagonal_htm], + _scenario_dict=all_manuscript_scenarios, + _outputfilepath=figurespath, + _value_of_life_suffix=vsly_fig_suffixes[i], + _metric=chosen_metric, + _year_suffix=f' ({str(relevant_period_for_costing[0])} - {str(relevant_period_for_costing[1])}) with breakeven implementation costs', + _projected_health_spending=projected_health_spending_baseline, + _additional_horizontal_lines_for_interpretation=additional_horizontal_lines_for_interpretation, + _draw_colors=draw_colors, + show_title_and_legend=False) + else: additional_horizontal_lines_for_interpretation = None - legend_switch_for_main_roi_plot = True - draw_colors = {horizontal_hss: '#9e0142', vertical_htm: '#fdae61', diagonal_htm: '#66c2a5'} generate_multiple_scenarios_roi_plot( _monetary_value_of_incremental_health=get_monetary_value_of_incremental_health(num_dalys_averted, _chosen_value_of_life_year=vsly), @@ -1253,18 +1296,10 @@ def compute_breakeven_with_ci( _projected_health_spending=projected_health_spending_baseline, _additional_horizontal_lines_for_interpretation=None, _draw_colors=draw_colors, - show_title_and_legend=legend_switch_for_main_roi_plot) + show_title_and_legend=True) i = i + 1 - incremental_scenario_cost_summarised = convert_results_to_dict(incremental_scenario_cost_summarised) - print(f"Under an alternative assumption that the vertical approach incurs incremental above service level costs " - f"equal to 138% of its incremental service level cost (based on estimates from Opuni et al (2023)), " - f"the diagonal approach provided a higher ROI up to an even higher threshold of " - f"${breakeven_upper_v_htm[0]: .2f}[${breakeven_upper_v_htm[1]: .2f} - ${breakeven_upper_v_htm[2]: .2f}] million " - f"incremental above service level costs in comparison with the vertical approach, or equal to " - f"{(breakeven_upper_v_htm[0] * 1e6) / incremental_scenario_cost_summarised[diagonal_htm][chosen_metric] * 100: .2f}% of its own incremental service level cost") - # HIV scenarios with and without HSS draw_colors = {vertical_hiv: '#fdae61', diagonal_hiv: '#66c2a5'} generate_multiple_scenarios_roi_plot( @@ -1345,7 +1380,7 @@ def compute_breakeven_with_ci( _chosen_value_of_life_year=vsly), _incremental_input_cost=incremental_scenario_cost, _non_zero_implementation_cost_proportion=above_service_level_cost_proportion, - _draws=list(all_manuscript_scenarios.keys()), + _draws=list(main_manuscript_scenarios.keys()), _metric='median') roi_table_small['scenario'] = roi_table_small['draw'].map(all_manuscript_scenarios) # Drop 'draw' and move 'scenario' to the first column @@ -1453,8 +1488,7 @@ def compute_breakeven_with_ci( # 6. Plot costs # ---------------------------------------------------- # First summarize all input costs - - input_costs_subset_for_figure = input_costs[input_costs['draw'].isin(list(all_manuscript_scenarios.keys()))] + input_costs_subset_for_figure = input_costs[input_costs['draw'].isin(list(main_manuscript_scenarios.keys()))] agg_funcs = { chosen_metric: ('cost', chosen_metric), 'lower': ('cost', lambda x: x.quantile(0.025)), @@ -1498,61 +1532,9 @@ def compute_breakeven_with_ci( # %% # The following figures and results are generated only for the main (0.03 discount for both health and cost) scenario # ---------------------------------------------------------- -# Extract summary of incremental costs for appendix - -def calculate_detailed_incremental_costs(_df, comparison_draw=0): - """ - Calculate incremental costs relative to a specified comparison draw. - - Parameters: - - _df (pd.DataFrame): The input DataFrame with cost data. - - comparison_draw (int): The draw to use as the baseline for comparison. - - Returns: - - pd.DataFrame: DataFrame with an additional column 'incremental_cost'. - """ - # Identify all grouping columns (everything except 'draw' and 'cost') - group_cols = [col for col in _df.columns if col not in ['draw', 'cost']] - - # Extract baseline costs for the specified comparison draw - baseline_costs = _df[_df['draw'] == comparison_draw][group_cols + ['cost']] - baseline_costs = baseline_costs.rename(columns={'cost': 'baseline_cost'}) - - # Merge baseline costs with the original dataframe - merged_df = _df.merge(baseline_costs, on=group_cols, how='left') - - # Calculate incremental costs - merged_df['incremental_cost'] = merged_df['cost'] - merged_df['baseline_cost'] - merged_df.drop(columns=['baseline_cost', 'cost'], inplace=True) # Drop the baseline cost column - return merged_df[merged_df['draw'] != comparison_draw] - - incremental_costs_df = calculate_detailed_incremental_costs(input_costs_subset_for_figure, comparison_draw=0) - -def summarise_detailed_costs_df(_df, _metric, column_to_summarise='cost'): - # Summarize values - agg_func = np.mean if _metric == 'mean' else np.median - groupby_cols = [col for col in _df.columns if col not in ['run', column_to_summarise]] - _df = pd.concat( - { - chosen_metric: _df.groupby(by=groupby_cols, sort=False)[column_to_summarise].agg(agg_func), - 'lower': _df.groupby(by=groupby_cols, sort=False)[column_to_summarise].quantile(0.025), - 'upper': _df.groupby(by=groupby_cols, sort=False)[column_to_summarise].quantile(0.975), - }, - axis=1 - ) - - summarised_df = pd.melt( - _df.reset_index(), - id_vars=groupby_cols, # Columns to keep - value_vars=[_metric, 'lower', 'upper'], # Columns to unpivot - var_name='stat', # New column name for the 'sub-category' of cost - value_name=column_to_summarise - ) - return summarised_df - - +# Extract summary of incremental costs for appendix incremental_costs_summary = summarise_detailed_costs_df(_df=incremental_costs_df, _metric=chosen_metric, column_to_summarise='incremental_cost') incremental_costs_summary = \ @@ -1612,7 +1594,7 @@ def remove_discounting(_df, _discount_rate=0, _year=None): # Extracts for manuscript # -------------------------- # ICER results -icer_result = convert_results_to_dict(icers_summarized) +icer_result = convert_results_to_dict(icers_summarized, _metric = chosen_metric) hr_scenario_with_lowest_icer = min([1, 2, 3], key=lambda k: icer_result[k][chosen_metric]) hr_scenario_with_highest_icer = max([1, 2, 3], key=lambda k: icer_result[k][chosen_metric]) cons_scenario_with_lowest_icer = min([4, 5, 6], key=lambda k: icer_result[k][chosen_metric]) @@ -1630,16 +1612,15 @@ def remove_discounting(_df, _discount_rate=0, _year=None): list(all_manuscript_scenarios.keys()))] # keep only relevant draws # Summarize and convert to dictionary num_dalys_averted_v_vertical_summarised = summarize_cost_data(num_dalys_averted_v_vertical, _metric=chosen_metric) -dalys_averted_v_vertical_result = convert_results_to_dict(num_dalys_averted_v_vertical_summarised) +dalys_averted_v_vertical_result = convert_results_to_dict(num_dalys_averted_v_vertical_summarised, _metric = chosen_metric) print(f"The ICER of vertical strategy relative to the baseline scenario was " f"${icer_result[vertical_htm][chosen_metric]:.2f} [${icer_result[vertical_htm]['lower']:.2f} - ${icer_result[vertical_htm]['upper']:.2f}] " f"per DALY averted, assuming no additional implementation costs. The ICER of the diagonal HTM with HSS strategy " f"relative to the baseline scenario was " f"${icer_result[diagonal_htm][chosen_metric]:.2f} [${icer_result[diagonal_htm]['lower']:.2f} - ${icer_result[diagonal_htm]['upper']:.2f}], " - f"demonstrating that the diagonal strategy was more cost-effective than the vertical one. While the horizontal " - f"strategy, HSS expansion, averted " + f"The horizontal strategy, HSS expansion, averted " f"{dalys_averted_v_vertical_result[horizontal_hss][chosen_metric] / 1e6:.2f} [{dalys_averted_v_vertical_result[horizontal_hss]['lower'] / 1e6:.2f} - {dalys_averted_v_vertical_result[horizontal_hss]['upper'] / 1e6:.2f}] " - f"million more DALYs than HTM expansion, it did so at a higher cost per DALY averted " - f"(${icer_result[horizontal_hss][chosen_metric]:.2f} [${icer_result[horizontal_hss]['lower']:.2f} - ${icer_result[horizontal_hss]['upper']:.2f}] versus " - f"${icer_result[vertical_htm][chosen_metric]:.2f} [${icer_result[vertical_htm]['lower']:.2f} - ${icer_result[vertical_htm]['upper']:.2f}]), meaning it was less cost-effective.") + f"million more DALYs than HTM expansion and it did so at a lower cost " + f"(${convert_results_to_dict(incremental_scenario_cost_summarized)[horizontal_hss][chosen_metric] / 1e6:.2f}m [${convert_results_to_dict(incremental_scenario_cost_summarized)[horizontal_hss]['lower'] / 1e6:.2f}m - ${convert_results_to_dict(incremental_scenario_cost_summarized)[horizontal_hss]['upper'] / 1e6:.2f}m] versus " + f"${convert_results_to_dict(incremental_scenario_cost_summarized)[vertical_htm][chosen_metric] / 1e6:.2f}m [${convert_results_to_dict(incremental_scenario_cost_summarized)[vertical_htm]['lower'] / 1e6:.2f}m - ${convert_results_to_dict(incremental_scenario_cost_summarized)[vertical_htm]['upper'] / 1e6:.2f}m]), meaning it was dominant.") From 603d0be7ff2fa1d8d53482f603be0819dc2ba80b Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Wed, 1 Oct 2025 19:12:04 +0100 Subject: [PATCH 08/14] update the order in which scenarios appear in extracted csvs --- .../roi_analysis_horizontal_vs_vertical.py | 23 +++++++++++++++---- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index deda3879cb..93d3fa2ef3 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -601,6 +601,15 @@ def get_monetary_value_of_incremental_health(_num_dalys_averted, _chosen_value_o monetary_value_of_incremental_health = (_num_dalys_averted * _chosen_value_of_life_year).clip(lower=0.0) return monetary_value_of_incremental_health +def reorder_rows_for_csv(_df, desired_order): + if desired_order is None: + return _df + else: + present = [s for s in desired_order if s in _df['scenario'].unique()] + _df['scenario'] = pd.Categorical(_df['scenario'], categories=present, ordered=True) + df = _df.sort_values('scenario') + return df + # %% CONFIG # Load result files # ------------------------------------------------------------------------------------------------------------------ @@ -698,6 +707,10 @@ def get_monetary_value_of_incremental_health(_num_dalys_averted, _chosen_value_o diagonal_malaria = all_manuscript_scenarios_reverse.get("Malaria Program Scale-up With HSS Expansion Package") diagonal_htm = all_manuscript_scenarios_reverse.get("HTM Programs Scale-up With HSS Expansion Package") +# Desired order of scenarios for tables +ordered_ids = [1,2,3,4,5,6,7,8,16,24,32,15,23,31,39] +desired_order = [all_manuscript_scenarios[i] for i in ordered_ids] + # Use letters instead of full scenario name for figures all_manuscript_scenarios_substitutedict = {0: "0", 1: "A", 2: "B", 3: "C", 4: "D", 5: "E", 6: "F", 7: "G", 8: "H", 15: "I", 16: "J", 23: "K", 24: "L", 31: "M", 32: "N", 39: "O"} @@ -1089,7 +1102,7 @@ def melt_and_label_malaria_scaleup_cost(_df, label): ), axis=1 ) - icers_summarized_subset_for_table[['scenario', 'ICER (2023 USD)']].to_csv(figurespath / 'tabulated_icers.csv', + reorder_rows_for_csv(icers_summarized_subset_for_table[['scenario', 'ICER (2023 USD)']], desired_order).to_csv(figurespath / 'tabulated_icers.csv', index=False) # Create a lookup from draw to formatted ICER string @@ -1352,7 +1365,7 @@ def melt_and_label_malaria_scaleup_cost(_df, label): _monetary_value_of_incremental_health=get_monetary_value_of_incremental_health(num_dalys_averted, _chosen_value_of_life_year=vsly), _incremental_input_cost=incremental_scenario_cost, - _draws=list(all_manuscript_scenarios.keys()), + _draws=list(main_manuscript_scenarios.keys()), _metric='median') # Extract ROIs into a table for manuscript @@ -1371,8 +1384,8 @@ def melt_and_label_malaria_scaleup_cost(_df, label): # Reset index to move 'implementation_cost' to columns and reshape for the manuscript roi_pivot_table = roi_pivot_table['ROI at VSLY = $834'].unstack(level='implementation_cost') roi_pivot_table.columns.name = None # Remove multi-index column name - roi_pivot_table.index.name = 'Scenario' - roi_pivot_table.reset_index().drop(columns='draw').to_csv( + roi_pivot_table.index.name = 'scenario' + reorder_rows_for_csv(roi_pivot_table.reset_index(), desired_order).reset_index().drop(columns='draw').to_csv( figurespath / f'tabulated_roi_for_all_implementation_costs_{roi_table_label[i]}.csv', index=False) roi_table_small = extract_roi_at_specific_implementation_costs( @@ -1387,7 +1400,7 @@ def melt_and_label_malaria_scaleup_cost(_df, label): roi_table_small = roi_table_small.drop(columns='draw') cols = ['scenario'] + [col for col in roi_table_small.columns if col != 'scenario'] roi_table_small = roi_table_small[cols] - roi_table_small.to_csv(figurespath / f'tabulated_roi_for_manuscript_{roi_table_label[i]}.csv', index=False) + reorder_rows_for_csv(roi_table_small.reset_index(), desired_order).to_csv(figurespath / f'tabulated_roi_for_manuscript_{roi_table_label[i]}.csv', index=False) i += 1 # ROI Bar plots From b175cb6ab2b1db480930188b1b7c4e2bb280c192 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Wed, 1 Oct 2025 19:12:52 +0100 Subject: [PATCH 09/14] remove superfluous imports --- .../roi_analysis_horizontal_vs_vertical.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index 93d3fa2ef3..20cb706563 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -9,15 +9,13 @@ import os import textwrap from pathlib import Path -from typing import Optional, Literal, Callable, Iterable -from dataclasses import dataclass +from typing import Optional, Literal import matplotlib.colors as mcolors import matplotlib.patches as mpatches import matplotlib.pyplot as plt import numpy as np import pandas as pd -from pandas import DataFrame, Series import seaborn as sns from adjustText import adjust_text # For the CEA plane figure to avoid overlaps in data labels From 766a4529170a72cfb04e8fbaa94d65c51a743e0f Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Fri, 14 Nov 2025 19:37:28 +0000 Subject: [PATCH 10/14] new results folder --- .../roi_analysis_horizontal_vs_vertical.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index 20cb706563..8d0fe6d7a4 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -611,7 +611,7 @@ def reorder_rows_for_csv(_df, desired_order): # %% CONFIG # Load result files # ------------------------------------------------------------------------------------------------------------------ -results_folder = get_scenario_outputs('htm_and_hss_runs-2025-09-16T141811Z.py', outputfilepath)[0] +results_folder = get_scenario_outputs('htm_and_hss_runs-2025-10-14T084418Z.py', outputfilepath)[0] # Check can read results from draw=0, run=0 log = load_pickled_dataframes(results_folder, 0, 0) # look at one log (so can decide what to extract) @@ -1082,7 +1082,7 @@ def melt_and_label_malaria_scaleup_cost(_df, label): ) ax.set_title(name_of_plot) ax.set_ylabel('ICERs \n($/DALY averted)') - ax.set_ylim(bottom=0) + ax.set_ylim(bottom=0, top = icers_summarized_subset_for_figure['upper'].max() * 1.1) fig.tight_layout() fig.savefig(figurespath / name_of_plot.replace(' ', '_').replace(',', '')) plt.close(fig) From 94c8ae91fd451e0440ff04d0bb4df4e7e445e239 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Wed, 19 Nov 2025 21:24:22 +0000 Subject: [PATCH 11/14] Add consumables plot + change title of ROI plots --- .../roi_analysis_horizontal_vs_vertical.py | 64 ++++++++++++++++++- 1 file changed, 62 insertions(+), 2 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index 8d0fe6d7a4..5cd0706517 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -1407,7 +1407,7 @@ def melt_and_label_malaria_scaleup_cost(_df, label): roi_at_upper_limit_implementation_cost_for_figure = roi_at_upper_limit_implementation_cost_summarized[roi_at_upper_limit_implementation_cost_summarized.index.isin(chosen_draws)] draw_colors = {horizontal_hss: '#9e0142', vertical_htm: '#fdae61', diagonal_htm: '#66c2a5'} - name_of_plot = f'ROI assuming zero implementation costs {relevant_period_for_costing[0]}-{relevant_period_for_costing[1]}' + name_of_plot = f'ROI assuming zero above service-level costs' fig, ax = do_standard_bar_plot_with_ci( (roi_at_0_implementation_cost_for_figure), annotations=[ @@ -1426,7 +1426,7 @@ def melt_and_label_malaria_scaleup_cost(_df, label): fig.savefig(figurespath / name_of_plot.replace(' ', '_').replace(',', '')) plt.close(fig) - name_of_plot = f'ROI assuming non-zero implementation costs {relevant_period_for_costing[0]}-{relevant_period_for_costing[1]}' + name_of_plot = f'ROI assuming non-zero above service-level costs' fig, ax = do_standard_bar_plot_with_ci( (roi_at_upper_limit_implementation_cost_for_figure), annotations=[ @@ -1635,3 +1635,63 @@ def remove_discounting(_df, _discount_rate=0, _year=None): f"million more DALYs than HTM expansion and it did so at a lower cost " f"(${convert_results_to_dict(incremental_scenario_cost_summarized)[horizontal_hss][chosen_metric] / 1e6:.2f}m [${convert_results_to_dict(incremental_scenario_cost_summarized)[horizontal_hss]['lower'] / 1e6:.2f}m - ${convert_results_to_dict(incremental_scenario_cost_summarized)[horizontal_hss]['upper'] / 1e6:.2f}m] versus " f"${convert_results_to_dict(incremental_scenario_cost_summarized)[vertical_htm][chosen_metric] / 1e6:.2f}m [${convert_results_to_dict(incremental_scenario_cost_summarized)[vertical_htm]['lower'] / 1e6:.2f}m - ${convert_results_to_dict(incremental_scenario_cost_summarized)[vertical_htm]['upper'] / 1e6:.2f}m]), meaning it was dominant.") + +# Generate consumable availability plot +#--------------------------------------- +# Load availability data +path_for_cons_resourcefiles = resourcefilepath / "healthsystem/consumables" +full_df_with_scenario = pd.read_csv(path_for_cons_resourcefiles / "ResourceFile_Consumables_availability_small.csv") + +# Import item_category +program_item_mapping = pd.read_csv(path_for_cons_resourcefiles / 'ResourceFile_Consumables_Item_Designations.csv')[['Item_Code', 'item_category']] +program_item_mapping = program_item_mapping.rename(columns ={'Item_Code': 'item_code'})[program_item_mapping.item_category.notna()] + +# Get TLO Facility_ID for each district and facility level +mfl = pd.read_csv(resourcefilepath / "healthsystem" / "organisation" / "ResourceFile_Master_Facilities_List.csv") +districts = set(pd.read_csv(resourcefilepath / 'demography' / 'ResourceFile_Population_2010.csv')['District']) +fac_levels = {'0', '1a', '1b', '2', '3', '4'} + +df_for_plots = full_df_with_scenario.merge(mfl[['Facility_ID', 'Facility_Level']], on = 'Facility_ID', how = 'left', validate = "m:1") +df_for_plots = df_for_plots.merge(program_item_mapping, on = 'item_code', how = 'left', validate = "m:1") + +# Choose scenarios to plot +scenario_list = [6,10,11] +chosen_availability_columns = ['available_prop'] + [f'available_prop_scenario{i}' for i in + scenario_list] +scenario_names_dict = {'available_prop': 'Actual', + 'available_prop_scenario6': '75th percentile\n facility', + 'available_prop_scenario10': 'HIV supply \n chain', 'available_prop_scenario11': 'EPI supply \n chain'} +# recreate the chosen columns list based on the mapping above +chosen_availability_columns = [scenario_names_dict[col] for col in chosen_availability_columns] +df_for_plots = df_for_plots.rename(columns = scenario_names_dict) + +for level in ['1a', '1b']: + # Generate a heatmap + # Pivot the DataFrame + aggregated_df = df_for_plots.groupby(['item_category', 'Facility_Level'])[chosen_availability_columns].mean().reset_index() + aggregated_df = aggregated_df[aggregated_df.Facility_Level.isin([level])] + heatmap_data = aggregated_df.set_index('item_category').drop(columns = 'Facility_Level') + + # Calculate the aggregate row and column + aggregate_col= df_for_plots.loc[df_for_plots.Facility_Level.isin([level]), chosen_availability_columns].mean() + #overall_aggregate = aggregate_col.mean() + + # Add aggregate row and column + #heatmap_data['Average'] = aggregate_row + #aggregate_col['Average'] = overall_aggregate + heatmap_data.loc['Average'] = aggregate_col + + # Generate the heatmap + sns.set(font_scale=0.8) + plt.figure(figsize=(10, 8)) + sns.heatmap(heatmap_data, annot=True, cmap='RdYlGn', cbar_kws={'label': 'Proportion of days on which consumable is available'}) + + # Customize the plot + plt.title(f'Facility Level {level}') + plt.xlabel('Scenarios') + plt.ylabel('Disease/Public health \n program') + plt.xticks(rotation=90, fontsize=8) + plt.yticks(rotation=0, fontsize=8) + + plt.savefig(figurespath /f'consumable_availability_heatmap_{level}.png', dpi=300, bbox_inches='tight') + plt.close() From 54b818a94b10ade021816b53b0370c3c05d8b606 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Thu, 20 Nov 2025 08:39:57 +0000 Subject: [PATCH 12/14] formatting edits to consumables plot --- .../roi_analysis_horizontal_vs_vertical.py | 46 +++++++++++-------- 1 file changed, 27 insertions(+), 19 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index 5cd0706517..bc96254911 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -1643,55 +1643,63 @@ def remove_discounting(_df, _discount_rate=0, _year=None): full_df_with_scenario = pd.read_csv(path_for_cons_resourcefiles / "ResourceFile_Consumables_availability_small.csv") # Import item_category -program_item_mapping = pd.read_csv(path_for_cons_resourcefiles / 'ResourceFile_Consumables_Item_Designations.csv')[['Item_Code', 'item_category']] -program_item_mapping = program_item_mapping.rename(columns ={'Item_Code': 'item_code'})[program_item_mapping.item_category.notna()] +program_item_mapping = pd.read_csv(path_for_cons_resourcefiles / 'ResourceFile_Consumables_Item_Designations.csv')[ + ['Item_Code', 'item_category']] +program_item_mapping = program_item_mapping.rename(columns={'Item_Code': 'item_code'})[ + program_item_mapping.item_category.notna()] # Get TLO Facility_ID for each district and facility level mfl = pd.read_csv(resourcefilepath / "healthsystem" / "organisation" / "ResourceFile_Master_Facilities_List.csv") districts = set(pd.read_csv(resourcefilepath / 'demography' / 'ResourceFile_Population_2010.csv')['District']) fac_levels = {'0', '1a', '1b', '2', '3', '4'} -df_for_plots = full_df_with_scenario.merge(mfl[['Facility_ID', 'Facility_Level']], on = 'Facility_ID', how = 'left', validate = "m:1") -df_for_plots = df_for_plots.merge(program_item_mapping, on = 'item_code', how = 'left', validate = "m:1") +df_for_plots = full_df_with_scenario.merge(mfl[['Facility_ID', 'Facility_Level']], on='Facility_ID', how='left', + validate="m:1") +df_for_plots = df_for_plots.merge(program_item_mapping, on='item_code', how='left', validate="m:1") # Choose scenarios to plot -scenario_list = [6,10,11] +scenario_list = [6, 10, 11] chosen_availability_columns = ['available_prop'] + [f'available_prop_scenario{i}' for i in - scenario_list] + scenario_list] scenario_names_dict = {'available_prop': 'Actual', - 'available_prop_scenario6': '75th percentile\n facility', - 'available_prop_scenario10': 'HIV supply \n chain', 'available_prop_scenario11': 'EPI supply \n chain'} + 'available_prop_scenario6': 'Consumables increased\n to 75th percentile', + 'available_prop_scenario10': 'Consumables increased\n to HIV level', + 'available_prop_scenario11': 'Consumables increased\n to EPI level'} # recreate the chosen columns list based on the mapping above chosen_availability_columns = [scenario_names_dict[col] for col in chosen_availability_columns] -df_for_plots = df_for_plots.rename(columns = scenario_names_dict) +df_for_plots = df_for_plots.rename(columns=scenario_names_dict) for level in ['1a', '1b']: # Generate a heatmap # Pivot the DataFrame - aggregated_df = df_for_plots.groupby(['item_category', 'Facility_Level'])[chosen_availability_columns].mean().reset_index() + aggregated_df = df_for_plots.groupby(['item_category', 'Facility_Level'])[ + chosen_availability_columns].mean().reset_index() aggregated_df = aggregated_df[aggregated_df.Facility_Level.isin([level])] - heatmap_data = aggregated_df.set_index('item_category').drop(columns = 'Facility_Level') + heatmap_data = aggregated_df.set_index('item_category').drop(columns='Facility_Level') # Calculate the aggregate row and column - aggregate_col= df_for_plots.loc[df_for_plots.Facility_Level.isin([level]), chosen_availability_columns].mean() - #overall_aggregate = aggregate_col.mean() + aggregate_col = df_for_plots.loc[df_for_plots.Facility_Level.isin([level]), chosen_availability_columns].mean() + # overall_aggregate = aggregate_col.mean() # Add aggregate row and column - #heatmap_data['Average'] = aggregate_row - #aggregate_col['Average'] = overall_aggregate + # heatmap_data['Average'] = aggregate_row + # aggregate_col['Average'] = overall_aggregate heatmap_data.loc['Average'] = aggregate_col # Generate the heatmap sns.set(font_scale=0.8) plt.figure(figsize=(10, 8)) - sns.heatmap(heatmap_data, annot=True, cmap='RdYlGn', cbar_kws={'label': 'Proportion of days on which consumable is available'}) + sns.heatmap(heatmap_data, annot=True, cmap='RdYlGn', + cbar_kws={'label': 'Proportion of days on which consumable is available'}, vmin=0, vmax=1, + annot_kws={"size": 12}) # Customize the plot plt.title(f'Facility Level {level}') plt.xlabel('Scenarios') plt.ylabel('Disease/Public health \n program') - plt.xticks(rotation=90, fontsize=8) - plt.yticks(rotation=0, fontsize=8) + plt.xticks(rotation=0, fontsize=10) + plt.yticks(rotation=0, fontsize=10) - plt.savefig(figurespath /f'consumable_availability_heatmap_{level}.png', dpi=300, bbox_inches='tight') + plt.savefig(figurespath / f'consumable_availability_heatmap_{level}.png', dpi=500, bbox_inches='tight') plt.close() + From 60c65237280cf237e04e88e062359b0087ce6362 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Thu, 20 Nov 2025 22:28:15 +0000 Subject: [PATCH 13/14] fix issues with negative costs/health impact when estimating ICERs and ROIs --- .../roi_analysis_horizontal_vs_vertical.py | 58 ++++++++++--------- 1 file changed, 32 insertions(+), 26 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index bc96254911..1f79561a37 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -773,6 +773,7 @@ def reorder_rows_for_csv(_df, desired_order): os.makedirs(figurespath) print(f"NOW RUNNING SCENARIO {rates['discounting_scenario']}...") + plt.style.use('default') # reset plot format # Estimate standard input costs of scenario # ----------------------------------------------------------------------------------------------------------------------- @@ -1052,14 +1053,41 @@ def melt_and_label_malaria_scaleup_cost(_df, label): # 3. Estimate and plot ICERs # ---------------------------------------------------- - icers = incremental_scenario_cost.div(num_dalys_averted) # Element-wise division + invalid = (num_dalys_averted < 0) & (incremental_scenario_cost > 0) + # TODO this is a quick and dirty workaround to calculate the ICER for HIV without HSS which has a negative health impact for only 1 out of 5 runs + # Might need to fix the above + icers = (incremental_scenario_cost[~invalid] / + num_dalys_averted[~invalid]) # Element-wise division icers_summarized = summarize_cost_data(icers, _metric=chosen_metric) - dominated_scenarios = icers_summarized['median'] < 0 - icers_summarized[dominated_scenarios] = np.nan # The dominated scenarios are assigned as ICER of NaN + dominant_scenarios = (summarize_cost_data(incremental_scenario_cost,_metric=chosen_metric)['median'] < 0) & (summarize_cost_data(num_dalys_averted,_metric=chosen_metric)['median'] > 0) + dominated_scenarios = (summarize_cost_data(incremental_scenario_cost,_metric=chosen_metric)['median'] > 0) & (summarize_cost_data(num_dalys_averted,_metric=chosen_metric)['median'] < 0) + icers_summarized[(dominated_scenarios|dominant_scenarios)] = np.nan # The dominated and dominant scenarios are assigned as ICER of NaN icers_summarized_subset_for_figure = icers_summarized[ icers_summarized.index.get_level_values('draw').isin(list(main_manuscript_scenarios.keys()))] + # Extract ICERs into a table for manuscript + icers_summarized_subset_for_table = icers_summarized_subset_for_figure.reset_index() + icers_summarized_subset_for_table['scenario'] = icers_summarized_subset_for_table['draw'].map( + all_manuscript_scenarios) + icers_summarized_subset_for_table['ICER (2023 USD)'] = icers_summarized_subset_for_table.apply( + lambda + row: "Dominated" # if incremental health is negative, the scenario is dominated - reporting a negative ICER can be confusing + if row['median'] < 0 else ( + f"${row['median']:.2f} [${row['lower']:.2f} - ${row['upper']:.2f}]" + if not any(pd.isna(row[['median', 'lower', 'upper']])) else "Dominant" # TODO add the Dominated case + ), + axis=1 + ) + reorder_rows_for_csv(icers_summarized_subset_for_table[['scenario', 'ICER (2023 USD)']], desired_order).to_csv(figurespath / 'tabulated_icers.csv', + index=False) + + # Create a lookup from draw to formatted ICER string + icer_lookup = dict(zip( + icers_summarized_subset_for_table['draw'], + icers_summarized_subset_for_table['ICER (2023 USD)'] + )) + # Plot ICERs annotations_icers = [] for _, row in icers_summarized_subset_for_figure.iterrows(): @@ -1087,28 +1115,6 @@ def melt_and_label_malaria_scaleup_cost(_df, label): fig.savefig(figurespath / name_of_plot.replace(' ', '_').replace(',', '')) plt.close(fig) - # Extract ICERs into a table for manuscript - icers_summarized_subset_for_table = icers_summarized_subset_for_figure.reset_index() - icers_summarized_subset_for_table['scenario'] = icers_summarized_subset_for_table['draw'].map( - all_manuscript_scenarios) - icers_summarized_subset_for_table['ICER (2023 USD)'] = icers_summarized_subset_for_table.apply( - lambda - row: "Dominated" # if incremental health is negative, the scenario is dominated - reporting a negative ICER can be confusing - if row['median'] < 0 else ( - f"${row['median']:.2f} [${row['lower']:.2f} - ${row['upper']:.2f}]" - if not any(pd.isna(row[['median', 'lower', 'upper']])) else "Dominated" - ), - axis=1 - ) - reorder_rows_for_csv(icers_summarized_subset_for_table[['scenario', 'ICER (2023 USD)']], desired_order).to_csv(figurespath / 'tabulated_icers.csv', - index=False) - - # Create a lookup from draw to formatted ICER string - icer_lookup = dict(zip( - icers_summarized_subset_for_table['draw'], - icers_summarized_subset_for_table['ICER (2023 USD)'] - )) - # Plot incremental health and cost in a scatterplot # ----------------------------------------------------= @@ -1175,7 +1181,7 @@ def melt_and_label_malaria_scaleup_cost(_df, label): roi_at_0_implementation_cost_dict = convert_results_to_dict(roi_at_0_implementation_cost_summarized, _metric = chosen_metric) # ROI at implementation costs = 138% of input costs - implementation_cost_upper_limit = incremental_scenario_cost * above_service_level_cost_proportion + implementation_cost_upper_limit = np.maximum(incremental_scenario_cost * above_service_level_cost_proportion,0) implementation_cost_upper_limit_dict = convert_results_to_dict( summarize_cost_data(implementation_cost_upper_limit, _metric=chosen_metric)) From d9b2560baf15738e724169f9ac83c6fdef13b731 Mon Sep 17 00:00:00 2001 From: sakshimohan Date: Tue, 25 Nov 2025 10:51:59 +0000 Subject: [PATCH 14/14] edit roi bar plots for inset figures --- .../roi_analysis_horizontal_vs_vertical.py | 26 ++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py index 1f79561a37..c72ad13062 100644 --- a/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py +++ b/src/scripts/comparison_of_horizontal_and_vertical_programs/economic_analysis_for_manuscript/roi_analysis_horizontal_vs_vertical.py @@ -1413,6 +1413,11 @@ def melt_and_label_malaria_scaleup_cost(_df, label): roi_at_upper_limit_implementation_cost_for_figure = roi_at_upper_limit_implementation_cost_summarized[roi_at_upper_limit_implementation_cost_summarized.index.isin(chosen_draws)] draw_colors = {horizontal_hss: '#9e0142', vertical_htm: '#fdae61', diagonal_htm: '#66c2a5'} + full_xticklabels = [ + textwrap.fill(all_manuscript_scenarios[key], width=15) + for key in roi_at_0_implementation_cost_for_figure.index + ] + name_of_plot = f'ROI assuming zero above service-level costs' fig, ax = do_standard_bar_plot_with_ci( (roi_at_0_implementation_cost_for_figure), @@ -1420,13 +1425,21 @@ def melt_and_label_malaria_scaleup_cost(_df, label): f"{row['median']:.2f} ({row['lower'] :.2f}- {row['upper']:.2f})" for _, row in roi_at_0_implementation_cost_for_figure.iterrows() ], - xticklabels_horizontal_and_wrapped=False, + xticklabels_horizontal_and_wrapped=True, put_labels_in_legend=True, offset=0.2, set_colors = draw_colors ) - ax.set_title(name_of_plot) + ax.set_title(name_of_plot, fontsize = 12) + for text in ax.texts: # annotation font size + text.set_fontsize(14) ax.set_ylabel('Return on Investment') + ax.tick_params(axis='both', labelsize=14) # tick label font size + n = len(full_xticklabels) + # --- Extract bar positions (this fixes your misalignment) --- + bar_centers = [patch.get_x() + patch.get_width() / 2 + 0.25 for patch in ax.patches] + ax.set_xticks(bar_centers) + ax.set_xticklabels(full_xticklabels, rotation=0, ha='right', fontsize=12) ax.set_ylim(bottom=0) fig.tight_layout() fig.savefig(figurespath / name_of_plot.replace(' ', '_').replace(',', '')) @@ -1444,9 +1457,16 @@ def melt_and_label_malaria_scaleup_cost(_df, label): offset=0.2, set_colors = draw_colors ) - ax.set_title(name_of_plot) + ax.set_title(name_of_plot, fontsize = 12) + for text in ax.texts: # annotation font size + text.set_fontsize(14) ax.set_ylabel('Return on Investment') ax.set_ylim(bottom=0) + ax.tick_params(axis='both', labelsize=14) # tick label font size + n = len(full_xticklabels) + bar_centers = [patch.get_x() + patch.get_width() / 2 + 0.25 for patch in ax.patches] + ax.set_xticks(bar_centers) + ax.set_xticklabels(full_xticklabels, rotation=0, ha='right', fontsize=12) fig.tight_layout() fig.savefig(figurespath / name_of_plot.replace(' ', '_').replace(',', '')) plt.close(fig)