diff --git a/pyproject.toml b/pyproject.toml index 6069aa3c1f..6f7b92fcc9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,7 +35,7 @@ dependencies = [ "azure-storage-file-share", # For saving and loading simulation state "dill", - "sdv", + "sdv==1.24.1", ] description = "Thanzi la Onse Epidemiology Model" dynamic = ["version"] diff --git a/resources/ResourceFile_RTI/RTI_emulator.pkl b/resources/ResourceFile_RTI/RTI_emulator.pkl new file mode 100644 index 0000000000..5312c280f5 --- /dev/null +++ b/resources/ResourceFile_RTI/RTI_emulator.pkl @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5bc5d14a555d635888a7fb760158493c79969281ed574aa871d1296e1bdb4ec9 +size 35333224 diff --git a/src/scripts/healthsystem/impact_of_actual_vs_funded/analysis_extract_data.py b/src/scripts/healthsystem/impact_of_actual_vs_funded/analysis_extract_data.py new file mode 100644 index 0000000000..4f2f0b8f35 --- /dev/null +++ b/src/scripts/healthsystem/impact_of_actual_vs_funded/analysis_extract_data.py @@ -0,0 +1,317 @@ +"""Produce plots to show the health impact (deaths, dalys) each the healthcare system (overall health impact) when +running under different MODES and POLICIES (scenario_impact_of_actual_vs_funded.py)""" + +# short tclose -> ideal case +# long tclose -> status quo +import argparse +from pathlib import Path +from typing import Tuple + +import numpy as np +import pandas as pd +from matplotlib import pyplot as plt +from matplotlib.colors import Normalize +from scipy.optimize import curve_fit + +from tlo import Date +from tlo.analysis.utils import extract_results, summarize +from tlo.analysis.life_expectancy import get_life_expectancy_estimates + + +# Range of years considered +min_year = 2010 +max_year = 2040 + + +def apply(results_folder: Path, output_folder: Path, resourcefilepath: Path = None, ): + """Produce standard set of plots describing the effect of each TREATMENT_ID. + - We estimate the epidemiological impact as the EXTRA deaths that would occur if that treatment did not occur. + - We estimate the draw on healthcare system resources as the FEWER appointments when that treatment does not occur. + """ + + TARGET_PERIOD = (Date(min_year, 1, 1), Date(max_year, 1, 1)) + + # Definitions of general helper functions + make_graph_file_name = lambda stub: output_folder / f"{stub.replace('*', '_star_')}.png" # noqa: E731 + + def target_period() -> str: + """Returns the target period as a string of the form YYYY-YYYY""" + return "-".join(str(t.year) for t in TARGET_PERIOD) + + def get_parameter_names_from_scenario_file() -> Tuple[str]: + """Get the tuple of names of the scenarios from `Scenario` class used to create the results.""" + from scripts.healthsystem.impact_of_actual_vs_funded.scenario_impact_of_actual_vs_funded import ( + ImpactOfHealthSystemMode, + ) + e = ImpactOfHealthSystemMode() + return tuple(e._scenarios.keys()) + + def get_num_deaths(_df): + """Return total number of Deaths (total within the TARGET_PERIOD) + """ + return pd.Series(data=len(_df.loc[pd.to_datetime(_df.date).between(*TARGET_PERIOD)])) + + def get_num_dalys(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year.between(*[i.year for i in TARGET_PERIOD])] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum().sum() + ) + + def get_num_dalys_by_cause(_df): + """Return number of DALYs by cause by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year.between(*[i.year for i in TARGET_PERIOD])] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum() + ) + + def set_param_names_as_column_index_level_0(_df): + """Set the columns index (level 0) as the param_names.""" + ordered_param_names_no_prefix = {i: x for i, x in enumerate(param_names)} + names_of_cols_level0 = [ordered_param_names_no_prefix.get(col) for col in _df.columns.levels[0]] + assert len(names_of_cols_level0) == len(_df.columns.levels[0]) + _df.columns = _df.columns.set_levels(names_of_cols_level0, level=0) + return _df + + def find_difference_relative_to_comparison(_ser: pd.Series, + comparison: str, + scaled: bool = False, + drop_comparison: bool = True, + ): + """Find the difference in the values in a pd.Series with a multi-index, between the draws (level 0) + within the runs (level 1), relative to where draw = `comparison`. + The comparison is `X - COMPARISON`.""" + return _ser \ + .unstack(level=0) \ + .apply(lambda x: (x - x[comparison]) / (x[comparison] if scaled else 1.0), axis=1) \ + .drop(columns=([comparison] if drop_comparison else [])) \ + .stack() + + + def get_counts_of_hsi_by_treatment_id(_df): + """Get the counts of the short TREATMENT_IDs occurring""" + _counts_by_treatment_id = _df \ + .loc[pd.to_datetime(_df['date']).between(*TARGET_PERIOD), 'TREATMENT_ID'] \ + .apply(pd.Series) \ + .sum() \ + .astype(int) + return _counts_by_treatment_id.groupby(level=0).sum() + + year_target = 2023 + def get_counts_of_hsi_by_treatment_id_by_year(_df): + """Get the counts of the short TREATMENT_IDs occurring""" + _counts_by_treatment_id = _df \ + .loc[pd.to_datetime(_df['date']).dt.year ==year_target, 'TREATMENT_ID'] \ + .apply(pd.Series) \ + .sum() \ + .astype(int) + return _counts_by_treatment_id.groupby(level=0).sum() + + def get_counts_of_hsi_by_short_treatment_id(_df): + """Get the counts of the short TREATMENT_IDs occurring (shortened, up to first underscore)""" + _counts_by_treatment_id = get_counts_of_hsi_by_treatment_id(_df) + _short_treatment_id = _counts_by_treatment_id.index.map(lambda x: x.split('_')[0] + "*") + return _counts_by_treatment_id.groupby(by=_short_treatment_id).sum() + + def get_counts_of_hsi_by_short_treatment_id_by_year(_df): + """Get the counts of the short TREATMENT_IDs occurring (shortened, up to first underscore)""" + _counts_by_treatment_id = get_counts_of_hsi_by_treatment_id_by_year(_df) + _short_treatment_id = _counts_by_treatment_id.index.map(lambda x: x.split('_')[0] + "*") + return _counts_by_treatment_id.groupby(by=_short_treatment_id).sum() + + + # Obtain parameter names for this scenario file + param_names = get_parameter_names_from_scenario_file() + print(param_names) + + # ================================================================================================ + # TIME EVOLUTION OF TOTAL DALYs + # Plot DALYs averted compared to the ``No Policy'' policy + + year_target = 2023 # This global variable will be passed to custom function + def get_num_dalys_by_year(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year == year_target] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum().sum() + ) + + ALL = {} + # Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred + # are consistent across different policies + this_min_year = 2010 + for year in range(this_min_year, max_year+1): + year_target = year + num_dalys_by_year = extract_results( + results_folder, + module='tlo.methods.healthburden', + key='dalys_stacked', + custom_generate_series=get_num_dalys_by_year, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = num_dalys_by_year + # Concatenate the DataFrames into a single DataFrame + concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + concatenated_df.index = concatenated_df.index.set_names(['date', 'index_original']) + concatenated_df = concatenated_df.reset_index(level='index_original',drop=True) + dalys_by_year = concatenated_df + print(dalys_by_year) + dalys_by_year.to_csv('ConvertedOutputs/Total_DALYs_with_time.csv', index=True) + + # ================================================================================================ + # Print population under each scenario + pop_model = extract_results(results_folder, + module="tlo.methods.demography", + key="population", + column="total", + index="date", + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + + pop_model.index = pop_model.index.year + pop_model = pop_model[(pop_model.index >= this_min_year) & (pop_model.index <= max_year)] + print(pop_model) + assert dalys_by_year.index.equals(pop_model.index) + assert all(dalys_by_year.columns == pop_model.columns) + pop_model.to_csv('ConvertedOutputs/Population_with_time.csv', index=True) + + # ================================================================================================ + # DALYs BROKEN DOWN BY CAUSES AND YEAR + # DALYs by cause per year + # %% Quantify the health losses associated with all interventions combined. + + year_target = 2023 # This global variable will be passed to custom function + def get_num_dalys_by_year_and_cause(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year == year_target] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum() + ) + + ALL = {} + # Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred + # are consistent across different policies + this_min_year = 2010 + for year in range(this_min_year, max_year+1): + year_target = year + num_dalys_by_year = extract_results( + results_folder, + module='tlo.methods.healthburden', + key='dalys_stacked', + custom_generate_series=get_num_dalys_by_year_and_cause, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = num_dalys_by_year #summarize(num_dalys_by_year) + + # Concatenate the DataFrames into a single DataFrame + concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + + concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) + + df_total = concatenated_df + df_total.to_csv('ConvertedOutputs/DALYS_by_cause_with_time.csv', index=True) + + ALL = {} + # Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred + # are consistent across different policies + for year in range(min_year, max_year+1): + year_target = year + + hsi_delivered_by_year = extract_results( + results_folder, + module='tlo.methods.healthsystem.summary', + key='HSI_Event', + custom_generate_series=get_counts_of_hsi_by_short_treatment_id_by_year, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = hsi_delivered_by_year + + # Concatenate the DataFrames into a single DataFrame + concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) + HSI_ran_by_year = concatenated_df + + del ALL + + ALL = {} + # Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred + # are consistent across different policies + for year in range(min_year, max_year+1): + year_target = year + + hsi_not_delivered_by_year = extract_results( + results_folder, + module='tlo.methods.healthsystem.summary', + key='Never_ran_HSI_Event', + custom_generate_series=get_counts_of_hsi_by_short_treatment_id_by_year, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = hsi_not_delivered_by_year + + # Concatenate the DataFrames into a single DataFrame + concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) + HSI_never_ran_by_year = concatenated_df + + HSI_never_ran_by_year = HSI_never_ran_by_year.fillna(0) #clean_df( + HSI_ran_by_year = HSI_ran_by_year.fillna(0) + HSI_total_by_year = HSI_ran_by_year.add(HSI_never_ran_by_year, fill_value=0) + HSI_ran_by_year.to_csv('ConvertedOutputs/HSIs_ran_by_area_with_time.csv', index=True) + HSI_never_ran_by_year.to_csv('ConvertedOutputs/HSIs_never_ran_by_area_with_time.csv', index=True) + print(HSI_ran_by_year) + print(HSI_never_ran_by_year) + print(HSI_total_by_year) + +if __name__ == "__main__": + rfp = Path('resources') + + parser = argparse.ArgumentParser( + description="Produce plots to show the impact each set of treatments", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument( + "--output-path", + help=( + "Directory to write outputs to. If not specified (set to None) outputs " + "will be written to value of --results-path argument." + ), + type=Path, + default=None, + required=False, + ) + parser.add_argument( + "--resources-path", + help="Directory containing resource files", + type=Path, + default=Path('resources'), + required=False, + ) + parser.add_argument( + "--results-path", + type=Path, + help=( + "Directory containing results from running " + "src/scripts/healthsystem/impact_of_actual_vs_funded/scenario_impact_of_actual_vs_funded.py " + ), + default=None, + required=False + ) + args = parser.parse_args() + assert args.results_path is not None + results_path = args.results_path + + output_path = results_path if args.output_path is None else args.output_path + + apply( + results_folder=results_path, + output_folder=output_path, + resourcefilepath=args.resources_path + ) diff --git a/src/scripts/healthsystem/impact_of_actual_vs_funded/scenario_impact_of_actual_vs_funded.py b/src/scripts/healthsystem/impact_of_actual_vs_funded/scenario_impact_of_actual_vs_funded.py new file mode 100644 index 0000000000..8c5cb7d932 --- /dev/null +++ b/src/scripts/healthsystem/impact_of_actual_vs_funded/scenario_impact_of_actual_vs_funded.py @@ -0,0 +1,117 @@ +"""This Scenario file run the model under different assumptions for HR capabilities expansion in order to estimate the +impact that is achieved under each. + +Run on the batch system using: +``` +tlo batch-submit src/scripts/healthsystem/impact_of_policy/scenario_impact_actual_vs_funded.py +``` + +or locally using: +``` +tlo scenario-run src/scripts/healthsystem/impact_of_policy/scenario_impact_actual_vs_funded.py +``` + +""" +from pathlib import Path +from typing import Dict + +import pandas as pd + +from tlo import Date, logging +from tlo.analysis.utils import get_parameters_for_status_quo, mix_scenarios +from tlo.methods.fullmodel import fullmodel +from tlo.methods.scenario_switcher import ImprovedHealthSystemAndCareSeekingScenarioSwitcher +from tlo.scenario import BaseScenario + + +class ImpactOfHealthSystemMode(BaseScenario): + def __init__(self): + super().__init__() + self.seed = 0 + self.start_date = Date(2010, 1, 1) + self.end_date = self.start_date + pd.DateOffset(years=31) + self.pop_size = 100_000 + self._scenarios = self._get_scenarios() + self.number_of_draws = len(self._scenarios) + self.runs_per_draw = 10 + + def log_configuration(self): + return { + 'filename': 'effect_of_actual_vs_funded', + 'directory': Path('./outputs'), # <- (specified only for local running) + 'custom_levels': { + '*': logging.WARNING, + 'tlo.methods.demography': logging.INFO, + 'tlo.methods.demography.detail': logging.WARNING, + 'tlo.methods.healthburden': logging.INFO, + 'tlo.methods.healthsystem.summary': logging.INFO, + } + } + + def modules(self): + return ( + fullmodel(resourcefilepath=self.resources) + + [ImprovedHealthSystemAndCareSeekingScenarioSwitcher(resourcefilepath=self.resources)] + ) + + def draw_parameters(self, draw_number, rng): + if draw_number < self.number_of_draws: + return list(self._scenarios.values())[draw_number] + else: + return + + def _get_scenarios(self) -> Dict[str, Dict]: + """Return the Dict with values for the parameters that are changed, keyed by a name for the scenario. + """ + + self.YEAR_OF_CHANGE = 2024 + + return { + + "No growth actual": + mix_scenarios( + self._baseline(), + { + "HealthSystem": { + "use_funded_or_actual_staffing_postSwitch": "actual", + }, + } + ), + + "No growth funded": + mix_scenarios( + self._baseline(), + { + "HealthSystem": { + "use_funded_or_actual_staffing_postSwitch": "funded", + }, + } + ), + + } + + def _baseline(self) -> Dict: + """Return the Dict with values for the parameter changes that define the baseline scenario. """ + return mix_scenarios( + get_parameters_for_status_quo(), + { + "HealthSystem": { + "mode_appt_constraints": 1, # <-- Mode 1 prior to change to preserve calibration + "mode_appt_constraints_postSwitch": 2, # <-- Mode 2 post-change to show effects of HRH + "year_mode_switch": self.YEAR_OF_CHANGE, + "use_funded_or_actual_staffing": 'actual', + "year_use_funded_or_actual_staffing_switch": self.YEAR_OF_CHANGE, + "scale_to_effective_capabilities": False, + "policy_name": "Naive", + "tclose_overwrite": 1, + "tclose_days_offset_overwrite": 7, + "cons_availability": "default", + } + }, + ) + + +if __name__ == '__main__': + from tlo.cli import scenario_run + + scenario_run([__file__]) diff --git a/src/scripts/rti_emulator/analysis_extract_data.py b/src/scripts/rti_emulator/analysis_extract_data.py new file mode 100644 index 0000000000..16ba64c492 --- /dev/null +++ b/src/scripts/rti_emulator/analysis_extract_data.py @@ -0,0 +1,569 @@ +"""Produce plots to show the health impact (deaths, dalys) each the healthcare system (overall health impact) when +running under different MODES and POLICIES (scenario_test_rti_emulator)""" + +# short tclose -> ideal case +# long tclose -> status quo +import argparse +from pathlib import Path +from typing import Tuple + +import numpy as np +import pandas as pd +from matplotlib import pyplot as plt +from matplotlib.colors import Normalize +from scipy.optimize import curve_fit + +from tlo import Date +from tlo.analysis.utils import extract_results, summarize +from tlo.analysis.life_expectancy import get_life_expectancy_estimates +from tlo.analysis.utils import ( + CAUSE_OF_DEATH_OR_DALY_LABEL_TO_COLOR_MAP, + extract_results, + format_gbd, + get_color_cause_of_death_or_daly_label, + load_pickled_dataframes, + make_age_grp_lookup, + make_age_grp_types, + make_calendar_period_lookup, + make_calendar_period_type, + order_of_cause_of_death_or_daly_label, + plot_clustered_stacked, + summarize, +) +tag = 'new_normal_All' + +#tag = 'new_emulated_with_conditionality_All' +#outputs/test_rti_emulator-2025-02-20T092114Z + +#tag = 'new_normal_All' +#outputs/test_rti_emulator-2025-02-10T095938Z + +#tag = 'new_emulated_with_conditionality_All_x250_other' +#outputs/test_rti_emulator-2025-02-19T153503Z + +#tag = 'new_normal_All_x250_other' +#outputs/test_rti_emulator-2025-02-10T171625Z +""" +""" +#tag = 'normal_Nothing' +#tag = 'normal' + +# Range of years considered +min_year = 2010 +max_year = 2019 + +appt_dict = {'Under5OPD': 'OPD', + 'Over5OPD': 'OPD', + 'AntenatalFirst': 'AntenatalTotal', + 'ANCSubsequent': 'AntenatalTotal', + 'NormalDelivery': 'Delivery', + 'CompDelivery': 'Delivery', + 'EstMedCom': 'EstAdult', + 'EstNonCom': 'EstAdult', + 'VCTPositive': 'VCTTests', + 'VCTNegative': 'VCTTests', + 'DentAccidEmerg': 'DentalAll', + 'DentSurg': 'DentalAll', + 'DentU5': 'DentalAll', + 'DentO5': 'DentalAll', + 'MentOPD': 'MentalAll', + 'MentClinic': 'MentalAll' + } + +def apply(results_folder: Path, output_folder: Path, resourcefilepath: Path = None, ): + """Produce standard set of plots describing the effect of each TREATMENT_ID. + - We estimate the epidemiological impact as the EXTRA deaths that would occur if that treatment did not occur. + - We estimate the draw on healthcare system resources as the FEWER appointments when that treatment does not occur. + """ + + TARGET_PERIOD = (Date(min_year, 1, 1), Date(max_year, 1, 1)) + + # Definitions of general helper functions + make_graph_file_name = lambda stub: output_folder / f"{stub.replace('*', '_star_')}.png" # noqa: E731 + + def target_period() -> str: + """Returns the target period as a string of the form YYYY-YYYY""" + return "-".join(str(t.year) for t in TARGET_PERIOD) + + def get_parameter_names_from_scenario_file() -> Tuple[str]: + """Get the tuple of names of the scenarios from `Scenario` class used to create the results.""" + from scripts.rti_emulator.scenario_test_rti_emulator import ( + GenerateDataChains, + ) + e = GenerateDataChains() + return tuple(e._scenarios.keys()) + + def get_num_deaths(_df): + """Return total number of Deaths (total within the TARGET_PERIOD) + """ + return pd.Series(data=len(_df.loc[pd.to_datetime(_df.date).between(*TARGET_PERIOD)])) + + def get_num_dalys(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year.between(*[i.year for i in TARGET_PERIOD])] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum().sum() + ) + + def get_num_dalys_by_cause(_df): + """Return number of DALYs by cause by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year.between(*[i.year for i in TARGET_PERIOD])] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum() + ) + + def set_param_names_as_column_index_level_0(_df): + """Set the columns index (level 0) as the param_names.""" + ordered_param_names_no_prefix = {i: x for i, x in enumerate(param_names)} + names_of_cols_level0 = [ordered_param_names_no_prefix.get(col) for col in _df.columns.levels[0]] + assert len(names_of_cols_level0) == len(_df.columns.levels[0]) + _df.columns = _df.columns.set_levels(names_of_cols_level0, level=0) + return _df + + def find_difference_relative_to_comparison(_ser: pd.Series, + comparison: str, + scaled: bool = False, + drop_comparison: bool = True, + ): + """Find the difference in the values in a pd.Series with a multi-index, between the draws (level 0) + within the runs (level 1), relative to where draw = `comparison`. + The comparison is `X - COMPARISON`.""" + return _ser \ + .unstack(level=0) \ + .apply(lambda x: (x - x[comparison]) / (x[comparison] if scaled else 1.0), axis=1) \ + .drop(columns=([comparison] if drop_comparison else [])) \ + .stack() + + + def get_counts_of_hsi_by_treatment_id(_df): + """Get the counts of the short TREATMENT_IDs occurring""" + _counts_by_treatment_id = _df \ + .loc[pd.to_datetime(_df['date']).between(*TARGET_PERIOD), 'TREATMENT_ID'] \ + .apply(pd.Series) \ + .sum() \ + .astype(int) + return _counts_by_treatment_id.groupby(level=0).sum() + + year_target = 2023 + def get_counts_of_hsi_by_treatment_id_by_year(_df): + """Get the counts of the short TREATMENT_IDs occurring""" + _counts_by_treatment_id = _df \ + .loc[pd.to_datetime(_df['date']).dt.year ==year_target, 'TREATMENT_ID'] \ + .apply(pd.Series) \ + .sum() \ + .astype(int) + return _counts_by_treatment_id.groupby(level=0).sum() + + def get_counts_of_hsi_by_short_treatment_id(_df): + """Get the counts of the short TREATMENT_IDs occurring (shortened, up to first underscore)""" + _counts_by_treatment_id = get_counts_of_hsi_by_treatment_id(_df) + _short_treatment_id = _counts_by_treatment_id.index.map(lambda x: x.split('_')[0] + "*") + return _counts_by_treatment_id.groupby(by=_short_treatment_id).sum() + + def get_counts_of_hsi_by_short_treatment_id_by_year(_df): + """Get the counts of the short TREATMENT_IDs occurring (shortened, up to first underscore)""" + _counts_by_treatment_id = get_counts_of_hsi_by_treatment_id_by_year(_df) + _short_treatment_id = _counts_by_treatment_id.index.map(lambda x: x.split('_')[0] + "*") + return _counts_by_treatment_id.groupby(by=_short_treatment_id).sum() + + def get_annual_num_appts_by_level(results_folder: Path) -> pd.DataFrame: + """Return pd.DataFrame gives the (mean) simulated annual number of appointments of each type at each level.""" + + def get_counts_of_appts(_df): + """Get the mean number of appointments of each type being used each year at each level. + Need to rename appts to match standardized categories from the DHIS2 data.""" + + def unpack_nested_dict_in_series(_raw: pd.Series): + return pd.concat( + { + idx: pd.DataFrame.from_dict(mydict) for idx, mydict in _raw.items() + } + ).unstack().fillna(0.0).astype(int) + + return _df \ + .loc[pd.to_datetime(_df['date']).between(*TARGET_PERIOD), 'Number_By_Appt_Type_Code_And_Level'] \ + .pipe(unpack_nested_dict_in_series) \ + .rename(columns=appt_dict, level=1) \ + .rename(columns={'1b': '2'}, level=0) \ + .groupby(level=[0, 1], axis=1).sum() \ + .mean(axis=0) # mean over each year (row) + + return extract_results( + results_folder, + module='tlo.methods.healthsystem.summary', + key='HSI_Event', + custom_generate_series=get_counts_of_appts, + do_scaling=True + ) + + + def get_annual_num_appts_by_level_with_confidence_interval(results_folder: Path) -> pd.DataFrame: + """Return pd.DataFrame gives the (mean) simulated annual number of appointments of each type at each level, + with 95% confidence interval.""" + #param_names = get_parameter_names_from_scenario_file() + def get_counts_of_appts(_df): + """Get the mean number of appointments of each type being used each year at each level. + Need to rename appts to match standardized categories from the DHIS2 data.""" + + def unpack_nested_dict_in_series(_raw: pd.Series): + return pd.concat( + { + idx: pd.DataFrame.from_dict(mydict) for idx, mydict in _raw.iteritems() + } + ).unstack().fillna(0.0).astype(int) + + return _df \ + .loc[pd.to_datetime(_df['date']).between(*TARGET_PERIOD), 'Number_By_Appt_Type_Code_And_Level'] \ + .pipe(unpack_nested_dict_in_series) \ + .rename(columns=appt_dict, level=1) \ + .rename(columns={'1b': '2'}, level=0) \ + .groupby(level=[0, 1], axis=1).sum() \ + .mean(axis=0) # mean over each year (row) + + return summarize( + extract_results( + results_folder, + module='tlo.methods.healthsystem.summary', + key='HSI_Event', + custom_generate_series=get_counts_of_appts, + do_scaling=True + ), + only_mean=False, + collapse_columns=True, + ).unstack().astype(int) + + + model = get_annual_num_appts_by_level(results_folder=results_folder) + model_summarised = summarize(model,only_mean=False,collapse_columns=True).unstack().astype(int) + model.to_csv('ConvertedOutputs/Emulator_Files/Total_Appt_Footprint_' + tag + '.csv', index=True) + model_summarised.to_csv('ConvertedOutputs/Emulator_Files/Total_Appt_Footprint_summarised_' + tag + '.csv', index=True) + #exit(-1) + # Obtain parameter names for this scenario file + param_names = get_parameter_names_from_scenario_file() + print(param_names) + + # ================================================================================================ + # TIME EVOLUTION OF TOTAL DALYs + # Plot DALYs averted compared to the ``No Policy'' policy + + year_target = 2023 # This global variable will be passed to custom function + def get_num_dalys_by_year(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year == year_target] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum().sum() + ) + + ALL = {} + # Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred + # are consistent across different policies + this_min_year = 2010 + for year in range(this_min_year, max_year+1): + year_target = year + num_dalys_by_year = extract_results( + results_folder, + module='tlo.methods.healthburden', + key='dalys_stacked', + custom_generate_series=get_num_dalys_by_year, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = num_dalys_by_year + # Concatenate the DataFrames into a single DataFrame + concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + concatenated_df.index = concatenated_df.index.set_names(['date', 'index_original']) + concatenated_df = concatenated_df.reset_index(level='index_original',drop=True) + dalys_by_year = concatenated_df + print(dalys_by_year) + dalys_by_year_summarise = summarize(dalys_by_year).unstack().astype(int) + + dalys_by_year.to_csv('ConvertedOutputs/Emulator_Files/Total_DALYs_with_time_' + tag + '.csv', index=True) + dalys_by_year_summarise.to_csv('ConvertedOutputs/Emulator_Files/Total_DALYs_with_time_summarised_' + tag + '.csv', index=True) + + # ================================================================================================ + # Print population under each scenario + pop_model = extract_results(results_folder, + module="tlo.methods.demography", + key="population", + column="total", + index="date", + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + + pop_model.index = pop_model.index.year + pop_model = pop_model[(pop_model.index >= this_min_year) & (pop_model.index <= max_year)] + print(pop_model) + assert dalys_by_year.index.equals(pop_model.index) + assert all(dalys_by_year.columns == pop_model.columns) + pop_model.to_csv('ConvertedOutputs/Emulator_Files/Population_with_time_' + tag + '.csv', index=True) + + # ================================================================================================ + # DEATHSs BROKEN DOWN BY CAUSES + # %% Quantify the health losses associated with all interventions combined. + + def get_num_deaths_by_cause_label(_df): + """Return total number of Deaths by label (total by age-group within the TARGET_PERIOD) + """ + return _df \ + .loc[pd.to_datetime(_df.date).between(*TARGET_PERIOD)] \ + .groupby(_df['label']) \ + .size() + + num_deaths_by_cause_label = summarize( + extract_results( + results_folder, + module='tlo.methods.demography', + key='death', + custom_generate_series=get_num_deaths_by_cause_label, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ) + print("Total number of deaths") + print(num_deaths_by_cause_label) + + year_target = 2023 # This global variable will be passed to custom function + def get_num_deaths_by_year_and_cause(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + newTARGET_PERIOD = (Date(year_target, 1, 1), Date(year_target, 12, 31)) + return _df \ + .loc[pd.to_datetime(_df.date).between(*newTARGET_PERIOD)] \ + .groupby(_df['label']) \ + .size() + #.loc[pd.to_datetime(_df.date).dt.year == year_target] \ + + ALL = {} + # Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred + # are consistent across different policies + this_min_year = 2010 + for year in range(this_min_year, max_year+1): + year_target = year + num_deaths_by_year_and_cause = extract_results( + results_folder, + module='tlo.methods.demography', + key='death', + custom_generate_series=get_num_deaths_by_year_and_cause, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = num_deaths_by_year_and_cause #summarize(num_dalys_by_year) + + # Concatenate the DataFrames into a single DataFrame + concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + + concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) + + df_total = concatenated_df + print(df_total) + print(df_total.groupby('cause').cumsum()) + print(summarize(df_total.groupby('cause').sum())) + df_total.to_csv('ConvertedOutputs/Emulator_Files/Deaths_by_cause_with_time_' + tag + '.csv', index=True) + df_total_summarise = summarize(df_total).unstack().astype(int) + df_total_summarise.to_csv('ConvertedOutputs/Emulator_Files/Deaths_by_cause_with_time_summarised_' + tag + '.csv', index=True) + + + # ================================================================================================ + # DALYs BROKEN DOWN BY CAUSES AND YEAR + # DALYs by cause per year + # %% Quantify the health losses associated with all interventions combined. + + year_target = 2023 # This global variable will be passed to custom function + def get_num_dalys_by_year_and_cause(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year == year_target] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum() + ) + + ALL = {} + # Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred + # are consistent across different policies + this_min_year = 2010 + for year in range(this_min_year, max_year+1): + year_target = year + num_dalys_by_year = extract_results( + results_folder, + module='tlo.methods.healthburden', + key='dalys_stacked', + custom_generate_series=get_num_dalys_by_year_and_cause, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = num_dalys_by_year #summarize(num_dalys_by_year) + + # Concatenate the DataFrames into a single DataFrame + concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + + concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) + + df_total = concatenated_df + print(df_total) + df_total.to_csv('ConvertedOutputs/Emulator_Files/DALYS_by_cause_with_time_' + tag + '.csv', index=True) + df_total_summarise = summarize(df_total).unstack().astype(int) + print(df_total_summarise) + df_total_summarise.to_csv('ConvertedOutputs/Emulator_Files/DALYS_by_cause_with_time_summarised_' + tag + '.csv', index=True) + + + + ALL = {} + # Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred + # are consistent across different policies + for year in range(min_year, max_year+1): + year_target = year + + hsi_delivered_by_year = extract_results( + results_folder, + module='tlo.methods.healthsystem.summary', + key='HSI_Event', + custom_generate_series=get_counts_of_hsi_by_short_treatment_id_by_year, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = hsi_delivered_by_year + + # Concatenate the DataFrames into a single DataFrame + concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) + HSI_ran_by_year = concatenated_df + + del ALL + + ALL = {} + # Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred + # are consistent across different policies + for year in range(min_year, max_year+1): + year_target = year + + hsi_not_delivered_by_year = extract_results( + results_folder, + module='tlo.methods.healthsystem.summary', + key='Never_ran_HSI_Event', + custom_generate_series=get_counts_of_hsi_by_short_treatment_id_by_year, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = hsi_not_delivered_by_year + + # Concatenate the DataFrames into a single DataFrame + concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) + HSI_never_ran_by_year = concatenated_df + + HSI_never_ran_by_year = HSI_never_ran_by_year.fillna(0) #clean_df( + HSI_ran_by_year = HSI_ran_by_year.fillna(0) + HSI_total_by_year = HSI_ran_by_year.add(HSI_never_ran_by_year, fill_value=0) + HSI_ran_by_year.to_csv('ConvertedOutputs/Emulator_Files/HSIs_ran_by_area_with_time_' + tag + '.csv', index=True) + HSI_never_ran_by_year.to_csv('ConvertedOutputs/Emulator_Files/HSIs_never_ran_by_area_with_time_' + tag + '.csv', index=True) + print(HSI_ran_by_year) + print(HSI_never_ran_by_year) + print(HSI_total_by_year) + + def get_dalys_by_period_sex_agegrp_label(df): + """Sum the dalys by period, sex, age-group and label""" + _, calperiodlookup = make_calendar_period_lookup() + + df['age_grp'] = df['age_range'].astype(make_age_grp_types()) + df["period"] = df["year"].map(calperiodlookup).astype(make_calendar_period_type()) + df = df.drop(columns=['date', 'age_range', 'year']) + df = df.groupby(by=["period", "sex", "age_grp"]).sum().stack() + df.index = df.index.set_names('label', level=3) + return df + + results = extract_results( + results_folder, + module="tlo.methods.healthburden", + key="dalys_stacked_by_age_and_time", # <-- for DALYS stacked by age and time + custom_generate_series=get_dalys_by_period_sex_agegrp_label, + do_scaling=True + ) + + # divide by five to give the average number of deaths per year within the five year period: + results = results.div(5.0) + results_to_store = (results.loc['2010-2014'] + results.loc['2015-2019'])/2 + results_to_store.to_csv('ConvertedOutputs/Emulator_Files/DALYs_by_sex_age_' + tag + '.csv', index=True) + results_to_store = summarize(results_to_store) + results_to_store.to_csv('ConvertedOutputs/Emulator_Files/DALYs_by_sex_age_summarised_' + tag + '.csv', index=True) + + + def get_total_num_dalys_by_wealth_and_label(_df): + """Return the total number of DALYS in the TARGET_PERIOD by wealth and cause label.""" + wealth_cats = {5: '0-19%', 4: '20-39%', 3: '40-59%', 2: '60-79%', 1: '80-100%'} + + return _df \ + .loc[_df['year'].between(*[d.year for d in TARGET_PERIOD])] \ + .drop(columns=['date', 'year']) \ + .assign( + li_wealth=lambda x: x['li_wealth'].map(wealth_cats).astype( + pd.CategoricalDtype(wealth_cats.values(), ordered=True) + ) + ).melt(id_vars=['li_wealth'], var_name='label') \ + .groupby(by=['li_wealth', 'label'])['value'] \ + .sum() + + total_num_dalys_by_wealth_and_label = extract_results( + results_folder, + module="tlo.methods.healthburden", + key="dalys_by_wealth_stacked_by_age_and_time", + custom_generate_series=get_total_num_dalys_by_wealth_and_label, + do_scaling=True, + ) + + total_num_dalys_by_wealth_and_label.to_csv('ConvertedOutputs/Emulator_Files/DALYs_by_wealth_' + tag + '.csv', index=True) + + total_num_dalys_by_wealth_and_label = summarize(total_num_dalys_by_wealth_and_label, + collapse_columns=True, + only_mean=False, + ).unstack() + + total_num_dalys_by_wealth_and_label.to_csv('ConvertedOutputs/Emulator_Files/DALYs_by_wealth_summarised_' + tag + '.csv', index=True) + + + + +if __name__ == "__main__": + rfp = Path('resources') + + parser = argparse.ArgumentParser( + description="Produce plots to show the impact each set of treatments", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument( + "--output-path", + help=( + "Directory to write outputs to. If not specified (set to None) outputs " + "will be written to value of --results-path argument." + ), + type=Path, + default=None, + required=False, + ) + parser.add_argument( + "--resources-path", + help="Directory containing resource files", + type=Path, + default=Path('resources'), + required=False, + ) + parser.add_argument( + "--results-path", + type=Path, + help=( + "Directory containing results from running " + "src/scripts/rti_emulator/scenario_test_rti_emulator " + ), + default=None, + required=False + ) + args = parser.parse_args() + assert args.results_path is not None + results_path = args.results_path + + output_path = results_path if args.output_path is None else args.output_path + + apply( + results_folder=results_path, + output_folder=output_path, + resourcefilepath=args.resources_path + ) diff --git a/src/scripts/rti_emulator/final_analysis_extract_data.py b/src/scripts/rti_emulator/final_analysis_extract_data.py new file mode 100644 index 0000000000..ca702699a4 --- /dev/null +++ b/src/scripts/rti_emulator/final_analysis_extract_data.py @@ -0,0 +1,826 @@ +"""Produce plots to show the health impact (deaths, dalys) each the healthcare system (overall health impact) when +running under different MODES and POLICIES (scenario_impact_of_capabilities_expansion_scaling.py)""" + +# short tclose -> ideal case +# long tclose -> status quo +import argparse +from pathlib import Path +from typing import Tuple + +import pandas as pd +from tlo import Date +from tlo.analysis.utils import extract_results, make_calendar_period_type, make_calendar_period_lookup, make_age_grp_lookup, make_age_grp_types +import matplotlib.pyplot as plt + +# Archive of outputs +# Batch based on /Users/mm2908/Desktop/EmuIBM/Save_With_WellPerforming/emulators/latest_CTGANSynthesizer_epochs500_dsF_batch_size500_num_k_folds10_Nsubsample10000_InAndOutC_test_k_folding_UniformEncoder_CTGANtest3_repeat_seed42_k_fold0.pkl +# 'standard_RTI': 'outputs/test_rti_emulator-2025-08-12T205454Z' +# 'emulated_RTI': 'outputs/test_rti_emulator-2025-08-13T080302Z' +# 'standard_RTI_x250': 'outputs/test_rti_emulator-2025-09-06T055808Z' +# 'emulated_RTI_x250': 'outputs/test_rti_emulator-2025-09-05T085159Z' +# 'no_RTI': 'outputs/test_rti_emulator-2025-08-13T143342Z' + + +# Latest CTGAN runs + number of runs per scenario increased from 10 -> 50 +outputs = { + 'standard_RTI': {'results_folder' : Path('outputs/test_rti_emulator-2025-08-12T205454Z'), 'data': {}}, + 'emulated_RTI': {'results_folder' : Path('outputs/test_rti_emulator-2025-08-13T080302Z'), 'data' : {}}, + 'standard_RTI_x250': {'results_folder' : Path('outputs/test_rti_emulator-2025-09-06T055808Z'), 'data': {}}, + 'emulated_RTI_x250': {'results_folder' : Path('outputs/test_rti_emulator-2025-09-05T085159Z'), 'data' : {}}, + 'no_RTI': {'results_folder' : Path('outputs/test_rti_emulator-2025-08-13T143342Z'), 'data' : {}}, + } + + +def apply(results_folder: Path) -> Tuple: + + tag = results_folder.name + + def get_parameter_names_from_scenario_file() -> Tuple[str]: + """Get the tuple of names of the scenarios from `Scenario` class used to create the results.""" + from scripts.healthsystem.impact_of_const_capabilities_expansion.scenario_impact_of_capabilities_expansion_scaling import ( + ImpactOfHealthSystemMode, + ) + e = ImpactOfHealthSystemMode() + return tuple(e._scenarios.keys()) + + def set_param_names_as_column_index_level_0(_df): + """Set the columns index (level 0) as the param_names.""" + ordered_param_names_no_prefix = {i: x for i, x in enumerate(param_names)} + names_of_cols_level0 = [ordered_param_names_no_prefix.get(col) for col in _df.columns.levels[0]] + assert len(names_of_cols_level0) == len(_df.columns.levels[0]) + _df.columns = _df.columns.set_levels(names_of_cols_level0, level=0) + return _df + + def get_counts_of_death_by_year(df): + df["year"] = df["date"].dt.year + result = df.groupby(by=["year", "label"])["person_id"].count() + return result + + def get_num_dalys_by_year(_df): + drop_cols = ["date", "sex", "age_range"] + result = _df.drop(columns=drop_cols).groupby("year", as_index=True).sum().stack() + result.index.set_names(["year", "label"], inplace=True) + return result + + def get_dalys_by_period_sex_agegrp_label(df): + """Sum the dalys by period, sex, age-group and label""" + _, calperiodlookup = make_calendar_period_lookup() + + df['age_grp'] = df['age_range'].astype(make_age_grp_types()) + df["period"] = df["year"].map(calperiodlookup).astype(make_calendar_period_type()) + df = df.drop(columns=['date', 'age_range', 'year']) + df = df.groupby(by=["period", "sex", "age_grp"]).sum().stack() + df.index = df.index.set_names('label', level=3) + return df + + def rename_index_levels(df): + rename_dict = {} + for level in df.index.names: + if level == 'date': + rename_dict['date'] = 'year' + elif level == 'label': + rename_dict['label'] = 'cause' + + if rename_dict: + df = df.rename_axis(index=rename_dict) + + if "year" in df.index.names: + years = df.index.get_level_values("year") + mask = (years >= 2010) & (years < 2020) + df = df[mask] + return df + + param_names = get_parameter_names_from_scenario_file() + + num_deaths_by_year_and_cause = rename_index_levels(extract_results( + results_folder, + module="tlo.methods.demography", + key="death", + custom_generate_series=get_counts_of_death_by_year, + do_scaling=True + )) + + num_dalys_by_year_and_cause = rename_index_levels(extract_results( + results_folder, + module="tlo.methods.healthburden", + key="dalys_stacked", + custom_generate_series=get_num_dalys_by_year, + do_scaling=True + )) + + num_individuals = extract_results(results_folder, + module="tlo.methods.demography", + key="population", + column="total", + index="date", + do_scaling=True + ) + num_individuals.index = num_individuals.index.year + num_individuals = rename_index_levels(num_individuals) + + cfr = num_deaths_by_year_and_cause.div(num_individuals, level="year") + + dalys_by_sex_and_age_time = extract_results( + results_folder, + module="tlo.methods.healthburden", + key="dalys_stacked_by_age_and_time", # <-- for DALYS stacked by age and time + custom_generate_series=get_dalys_by_period_sex_agegrp_label, + do_scaling=True + ) + # divide by five to give the average number of deaths per year within the five year period: + dalys_by_sex_and_age_time = dalys_by_sex_and_age_time.div(5.0) + dalys_by_sex_and_age = (dalys_by_sex_and_age_time.loc['2010-2014'] + dalys_by_sex_and_age_time.loc['2015-2019'])/2 + + # Collect all data for this output + data = {'deaths' : num_deaths_by_year_and_cause, 'dalys' : num_dalys_by_year_and_cause, 'pop' : num_individuals, 'cfr' : cfr, 'dalys_by_sex_and_age' : dalys_by_sex_and_age} + + return data + + +# Extract and calculate summary statistics +def compute_summary_stats(df): + df_mean = df.mean(axis=1) + df[('0','mean')] = df_mean + df_lower = df.quantile((1-0.682)/2.0,axis=1) + df_upper = df.quantile((1-0.682)/2.0+0.682,axis=1) + df[('0','mean')] = df_mean + df[('0','lower')] = df_lower + df[('0','upper')] = df_upper + return df + + +def compare_and_make_plots_age_sex_distr(outputs, first_scenario, second_scenario, third_scenario, fourth_scenario, target, compare_to_other_x250=True): + # Load from outputs + df_map = { + "normal": outputs[first_scenario]['data'][target], + "emulated": outputs[second_scenario]['data'][target], + "normal_x250": outputs[third_scenario]['data'][target], + "emulated_x250": outputs[fourth_scenario]['data'][target], + } + labels = [first_scenario, second_scenario, third_scenario, fourth_scenario] + + # Focus cause + cause = 'Transport Injuries' + for k, df in df_map.items(): + df_map[k] = compute_summary_stats(df.loc[(slice(None), slice(None), cause)]) + + # Age group mapping + age_grps = df_map["normal"].index.get_level_values('age_grp').unique() + age_grp_mapping = {age: i for i, age in enumerate(age_grps)} + + # Init p-values + pvalues = pd.DataFrame(columns=['p-value ks', 'p-value t', 'p-value w'], + index=pd.MultiIndex.from_tuples([], names=['case', 'sex', 'age_grp'])) + """ + # ---- Helper: compute p-values + def run_tests(case, df1, df2): + for sex in ['F', 'M']: + for age in age_grps: + data1, data2 = df1.loc[(sex, age)].values, df2.loc[(sex, age)].values + p_ks = ks_2samp(data1, data2).pvalue + p_t = stats.ttest_ind(data1, data2).pvalue + p_w = stats.wilcoxon(data1, data2).pvalue + pvalues.loc[(case, sex, age), :] = [p_ks, p_t, p_w] + + run_tests("other_x1", df_map["normal"], df_map["emulated"]) + if compare_to_other_x250: + run_tests("other_x250", df_map["normal_x250"], df_map["emulated_x250"]) + """ + + # ---- Plot + fig, axes = plt.subplots(1, 2, figsize=(14, 7), sharey=True) + for ax in axes: ax.set_ylim(-10, 120000) + + for i, sex in enumerate(['F', 'M']): + x_pos = [age_grp_mapping[a] for a in age_grps] + + def plot_band_and_mean(df, color, marker, label, alpha_fill, alpha_line, ls): + mean = df.loc[(sex, slice(None)), ('0', 'mean')].values + low = df.loc[(sex, slice(None)), ('0', 'lower')].values + up = df.loc[(sex, slice(None)), ('0', 'upper')].values + axes[i].fill_between(x_pos, low, up, color=color, alpha=alpha_fill) + axes[i].plot(x_pos, mean, marker=marker, color=color, alpha=alpha_line, linestyle=ls, label=label) + return mean + + # Plot main (normal vs emulated) + ls, fill, line = ('--', 0.1, 0.4) if compare_to_other_x250 else ('-', 0.3, 1.0) + m_norm = plot_band_and_mean(df_map["normal"], "blue", "o", labels[0] if i==0 else None, fill, line, ls) + m_emul = plot_band_and_mean(df_map["emulated"], "orange", "s", labels[1] if i==0 else None, fill, line, ls) + + # Optional comparison to x250 + if compare_to_other_x250: + m_norm_x250 = plot_band_and_mean(df_map["normal_x250"], "blue", "o", labels[2] if i==0 else None, 0.3, 1.0, '-') + m_emul_x250 = plot_band_and_mean(df_map["emulated_x250"], "orange", "s", labels[3] if i==0 else None, 0.3, 1.0, '-') + + """ + # Annotate KS values + y_base = 42000 if sex=="M" else 30000 + for case, offset in [("other_x1", 7000), ("other_x250", 0)] if compare_to_other_x250 else [("other_x1", 7000)]: + for age in age_grps: + axes[i].text(age_grp_mapping[age], y_base+offset, + f"{pvalues.loc[(case, sex, age), 'p-value ks']:.2f}", + fontsize=8, ha="center", rotation=90) + axes[i].text(age_grp_mapping['0-4'], y_base+offset+3000, + "KS test p-value" + (" (incr. other mortality)" if case=="other_x250" else ""), + fontsize=8, ha="left") + """ + # Formatting + axes[i].set_title(sex) + axes[i].set_xlabel("Age Group") + axes[i].set_ylabel("Averaged Yearly DALYs") + axes[i].set_xticks(x_pos) + axes[i].set_xticklabels(age_grps, rotation=45) + axes[i].legend() + + plt.tight_layout() + fname = "plots/final_DALYs_Breakdown_by_age_grp_DC.png" if compare_to_other_x250 else "plots/final_DALYs_Breakdown_by_age_grp.png" + plt.savefig(fname) + plt.close() + +def compare_old(outputs, first_scenario, second_scenario, third_scenario, fourth_scenario, target): + + # Load output files, which include standard other mortality and other mortality rates increased by x250 + df_af_normal = pd.read_csv('DALYs_by_sex_age_' + first_data_set + '.csv', header=[0, 1], index_col=[0,1,2]) + df_af_emulated = pd.read_csv('DALYs_by_sex_age_' + second_data_set + '.csv', header=[0, 1], index_col=[0,1,2]) + df_af_normal_other_x250 = pd.read_csv('DALYs_by_sex_age_' + first_data_set + '_x250_other.csv', header=[0, 1], index_col=[0,1,2]) + df_af_emulated_other_x250 = pd.read_csv('DALYs_by_sex_age_'+ second_data_set + '_x250_other.csv', header=[0, 1], index_col=[0,1,2]) + + # This is the cause that will be considered in the plots. Could also make loop over this + cause = 'Transport Injuries' + df_normal = df_af_normal.loc[(slice(None), slice(None), cause)] + df_normal_other_x250 = df_af_normal_other_x250.loc[(slice(None), slice(None), cause)] + df_emulated = df_af_emulated.loc[(slice(None), slice(None), cause)] + df_emulated_other_x250 = df_af_emulated_other_x250.loc[(slice(None), slice(None), cause)] + + # Map 'age_grp' to numeric values for easier plotting + age_grp_mapping = {age: i for i, age in enumerate(df_normal.index.get_level_values('age_grp').unique())} + + df_normal = compute_summary_stats(df_normal) + df_emulated = compute_summary_stats(df_emulated) + df_normal_other_x250 = compute_summary_stats(df_normal_other_x250) + df_emulated_other_x250 = compute_summary_stats(df_emulated_other_x250) + + # Store p-values for all x-values + index = pd.MultiIndex.from_tuples([], names=['case', 'sex', 'age_grp']) + columns = ['p-value ks', 'p-value t', 'p-value w'] + pvalues = pd.DataFrame(index=index, columns=columns) + + for sex in ['F','M']: + for age_group in df_normal_other_x250.index.get_level_values('age_grp').unique(): + + # Compute statistical tests for x1 other mortality + data_normal = df_normal_other_x250.loc[(sex,age_group)].values + data_emulat = df_emulated_other_x250.loc[(sex,age_group)].values + + statistic, p_value_ks = ks_2samp(data_normal, data_emulat) + t_stat, p_value_t = stats.ttest_ind(data_normal, data_emulat) + res_w = stats.wilcoxon(data_normal,data_emulat) + p_value_w = res_w.pvalue + pvalues.loc[('other_x250',sex,age_group),:] = [p_value_ks, p_value_t, p_value_w] + + # Compute statistical tests for x250 other mortality + data_normal = df_normal.loc[(sex,age_group)].values + data_emulat = df_emulated.loc[(sex,age_group)].values + + statistic, p_value_ks = ks_2samp(data_normal, data_emulat) + t_stat, p_value_t = stats.ttest_ind(data_normal, data_emulat) + res_w = stats.wilcoxon(data_normal,data_emulat) + p_value_w = res_w.pvalue + pvalues.loc[('other_x1',sex,age_group),:] = [p_value_ks, p_value_t, p_value_w] + + + # Create a figure for the subplots (one for each sex) + fig, axes = plt.subplots(1, 2, figsize=(14, 7), sharey=True) + for ax in axes: + ax.set_ylim(-1000, 55000) + + # Plot for each sex + for sex in ['F','M']: + if sex == 'F': + i = 0 + else: + i = 1 + + # Extract the age groups from the dataframe + age_groups = df_normal.index.get_level_values('age_grp').unique() + + # Create a numeric x_pos based on age groups for easier plotting + x_pos = [age_grp_mapping[age] for age in age_groups] + + # Prepare the values for plotting + normal_means = df_normal.loc[(sex, slice(None)),('0','mean')].values + emulated_means = df_emulated.loc[(sex, slice(None)),('0','mean')].values + normal_lower = df_normal.loc[(sex, slice(None)),('0','lower')].values + normal_upper = df_normal.loc[(sex, slice(None)),('0','upper')].values + emulated_lower = df_emulated.loc[(sex, slice(None)),('0','lower')].values + emulated_upper = df_emulated.loc[(sex, slice(None)),('0','upper')].values + + if compare_to_other_x250: + normal_means_other_x250 = df_normal_other_x250.loc[(sex, slice(None)),('0','mean')].values + emulated_means_other_x250 = df_emulated_other_x250.loc[(sex, slice(None)),('0','mean')].values + normal_lower_other_x250 = df_normal_other_x250.loc[(sex, slice(None)),('0','lower')].values + normal_upper_other_x250 = df_normal_other_x250.loc[(sex, slice(None)),('0','upper')].values + emulated_lower_other_x250 = df_emulated_other_x250.loc[(sex, slice(None)),('0','lower')].values + emulated_upper_other_x250 = df_emulated_other_x250.loc[(sex, slice(None)),('0','upper')].values + fill_main_comparison = 0.1 + line_main_compariosn = 0.4 + main_linestyle = '--' + print(main_linestyle) + else: + fill_main_comparison = 0.3 + line_main_compariosn = 1.0 + main_linestyle = '-' + + # Calculate % error: ((emulated - original) / original) * 100 + percent_error = 100 * (emulated_means - normal_means) / normal_means + if compare_to_other_x250: + percent_error_other_x250 = 100 * (emulated_means_other_x250 - normal_means_other_x250) / normal_means_other_x250 + + # Plot the shaded area for normal data (F or M) + axes[i].fill_between(x_pos, normal_lower, normal_upper, + color='blue', alpha=fill_main_comparison) + if compare_to_other_x250: + axes[i].fill_between(x_pos, normal_lower_other_x250, normal_upper_other_x250, + color='blue', alpha=0.3) + + # Plot the shaded area for emulated data (F or M) + axes[i].fill_between(x_pos, emulated_lower, emulated_upper, + color='orange', alpha=fill_main_comparison) + if compare_to_other_x250: + axes[i].fill_between(x_pos, emulated_lower_other_x250, emulated_upper_other_x250, + color='orange', alpha=0.3) + + # Plot the mean values for normal and emulated data + axes[i].plot(x_pos, normal_means, 'o-', color='blue', alpha=line_main_compariosn, linestyle=main_linestyle, label=labels[first_data_set] if i == 0 else "") # Normal mean + axes[i].plot(x_pos, emulated_means, 's-', color='orange', alpha=line_main_compariosn, linestyle=main_linestyle, label=labels[second_data_set] if i == 0 else "") # Emulated mean + if compare_to_other_x250: + axes[i].plot(x_pos, normal_means_other_x250, 'o-', color='blue', label=f'Original (incr. other mortality)' if i == 0 else "") # Normal mean + axes[i].plot(x_pos, emulated_means_other_x250, 's-', color='orange', label=f'Emulated (incr. other mortality)' if i == 0 else "") # Emulated mean + + # Include KS test p-values in the plot + if sex == 'M': + ks_yloc_other_x250 = 42000 + else: + ks_yloc_other_x250 = 30000 + ks_yloc_other_x1 = ks_yloc_other_x250 + 7000 + + for age_grp in df_normal.index.get_level_values('age_grp').unique(): + ks_value = pvalues.loc[('other_x1',sex, age_grp), 'p-value ks'] + axes[i].text( + x=age_grp_mapping[age_grp], + y=ks_yloc_other_x1, + s=f"{ks_value:.2f}", + fontsize=8, + ha='center', + rotation=90, + #bbox=dict(boxstyle="round,pad=0.1", facecolor="white", alpha=0.7) + ) + + #ks_value = pvalues.loc[(date, cause), 'p-value ks'] + leg_ypos = ks_yloc_other_x1+3000 + axes[i].text( + x=age_grp_mapping['0-4'], + y=leg_ypos, + s=f"KS test p-value", + fontsize=8, + ha='left', + rotation=0, + #bbox=dict(boxstyle="round,pad=0.2", facecolor="white", alpha=0.7) + ) + if compare_to_other_x250: + for age_grp in df_normal.index.get_level_values('age_grp').unique(): + ks_value = pvalues.loc[('other_x250',sex, age_grp), 'p-value ks'] + axes[i].text( + x=age_grp_mapping[age_grp], + y=ks_yloc_other_x250, + s=f"{ks_value:.2f}", + fontsize=8, + ha='center', + rotation=90, + #bbox=dict(boxstyle="round,pad=0.2", facecolor="white", alpha=0.7) + ) + + leg_ypos = ks_yloc_other_x250+3000 + axes[i].text( + x=age_grp_mapping['0-4'], + y=leg_ypos, + s=f"KS test p-value (incr. other mortality)", + fontsize=8, + ha='left', + rotation=0, + #bbox=dict(boxstyle="round,pad=0.2", facecolor="white", alpha=0.7) + ) + + # Set titles and labels + axes[i].set_title(f'{sex}') + axes[i].set_xlabel('Age Group') + axes[i].set_ylabel('Averaged Yearly DALYs') + + # Set x-ticks and labels + axes[i].set_xticks(x_pos) + axes[i].set_xticklabels(age_groups, rotation=45) + + # Add a legend + axes[i].legend() + + # Adjust layout for better spacing + plt.tight_layout() + + # Show the plot + if compare_to_other_x250: + plt.savefig('plots/final_DALYs_Breakdown_by_age_grp_DC.png') + else: + plt.savefig('plots/final_DALYs_Breakdown_by_age_grp.png') + + +def compare_and_plot(outputs, first_scenario, second_scenario, target, factor=None, ylabel=None): + + if factor is None: + df_first = outputs[first_scenario]['data'][target] # formally df + df_second = outputs[second_scenario]['data'][target] # formally df_emulated + else: + df_first = outputs[first_scenario]['data'][target]*factor # formally df + df_second = outputs[second_scenario]['data'][target]*factor # formally df_emulated + + # Define label names + labels = [first_scenario, second_scenario] + + # Iterate over causes + causes = df_first.index.get_level_values('cause').unique() + for cause in causes: + # Extract relevant data for both scenarios + data = {} + for label_name, dataset in zip(labels, [df_first, df_second]): + df_cause = dataset.loc[dataset.index.get_level_values('cause') == cause, ('0', ['lower', 'mean', 'upper'])] + data[label_name] = { + 'lower': df_cause[('0', 'lower')], + 'mean': df_cause[('0', 'mean')], + 'upper': df_cause[('0', 'upper')] + } + + years = data[labels[0]]['mean'].index.get_level_values('year') + + # % error + percent_error = 100 * (data[labels[1]]['mean'].values - data[labels[0]]['mean'].values) / data[labels[0]]['mean'].values + + # Create figure with two panels + fig, (ax_main, ax_error) = plt.subplots(2, 1, figsize=(10, 8), + gridspec_kw={'height_ratios': [3, 1]}, sharex=True) + + colors = ['blue', 'orange'] + markers = ['o', 'x'] + for i, label_name in enumerate(labels): + ax_main.plot(years, data[label_name]['mean'], label=label_name, color=colors[i], linestyle='-', marker=markers[i]) + ax_main.fill_between(years, data[label_name]['lower'], data[label_name]['upper'], color=colors[i], alpha=0.2) + + if ylabel is None: + ax_main.set_ylabel(f'{target}') + else: + ax_main.set_ylabel(f'{ylabel}') + ax_main.set_title(f'{target} due to {cause}') + ax_main.legend(loc='lower right') + ax_main.grid(True) + ax_main.set_ylim(0, 1.2e6) + + + # % error panel + print("For cause ", cause, " and target ", target, " MEAN PERC ERROR +++++++", percent_error.mean()) + ax_error.plot(years, percent_error, color='red', marker='o', linestyle='-', label='Mean') + ax_error.axhline(0, color='black', linestyle='--', linewidth=0.8) + ax_error.set_ylabel('% Error') + ax_error.set_xlabel('Year') + ax_error.legend(loc='upper right') + ax_error.grid(True) + + plt.tight_layout() + plt.savefig(f'plots/final_{target}_due_to_{cause}.png') + plt.close() + + +# First collect all data +for key in outputs.keys(): + outputs[key]['data'] = apply(outputs[key]['results_folder']) + for data_type in outputs[key]['data'].keys(): + outputs[key]['data'][data_type] = compute_summary_stats(outputs[key]['data'][data_type]) + + +compare_and_plot(outputs, 'standard_RTI', 'emulated_RTI', 'deaths', None, 'Deaths') +compare_and_plot(outputs, 'standard_RTI', 'emulated_RTI', 'dalys', None, 'DALYs') +compare_and_plot(outputs, 'standard_RTI', 'emulated_RTI', 'cfr', 1000, 'Crude mortality rate (/year/1000 individuals)') +compare_and_make_plots_age_sex_distr(outputs, 'standard_RTI', 'emulated_RTI', 'standard_RTI_x250', 'emulated_RTI_x250', 'dalys_by_sex_and_age', compare_to_other_x250=True) + + +exit(-1) + + + + +TARGET_PERIOD = (Date(min_year, 1, 1), Date(max_year, 1, 1)) + +# Definitions of general helper functions +lambda stub: output_folder / f"{stub.replace('*', '_star_')}.png" # noqa: E731 + +def target_period() -> str: + """Returns the target period as a string of the form YYYY-YYYY""" + return "-".join(str(t.year) for t in TARGET_PERIOD) + + + +def get_num_deaths(_df): + """Return total number of Deaths (total within the TARGET_PERIOD) + """ + return pd.Series(data=len(_df.loc[pd.to_datetime(_df.date).between(*TARGET_PERIOD)])) + +def get_num_dalys(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year.between(*[i.year for i in TARGET_PERIOD])] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum().sum() + ) + +def get_num_dalys_by_cause(_df): + """Return number of DALYs by cause by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year.between(*[i.year for i in TARGET_PERIOD])] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum() + ) + +def set_param_names_as_column_index_level_0(_df): + """Set the columns index (level 0) as the param_names.""" + ordered_param_names_no_prefix = {i: x for i, x in enumerate(param_names)} + names_of_cols_level0 = [ordered_param_names_no_prefix.get(col) for col in _df.columns.levels[0]] + assert len(names_of_cols_level0) == len(_df.columns.levels[0]) + _df.columns = _df.columns.set_levels(names_of_cols_level0, level=0) + return _df + +def find_difference_relative_to_comparison(_ser: pd.Series, + comparison: str, + scaled: bool = False, + drop_comparison: bool = True, + ): + """Find the difference in the values in a pd.Series with a multi-index, between the draws (level 0) + within the runs (level 1), relative to where draw = `comparison`. + The comparison is `X - COMPARISON`.""" + return _ser \ + .unstack(level=0) \ + .apply(lambda x: (x - x[comparison]) / (x[comparison] if scaled else 1.0), axis=1) \ + .drop(columns=([comparison] if drop_comparison else [])) \ + .stack() + + +def get_counts_of_hsi_by_treatment_id(_df): + """Get the counts of the short TREATMENT_IDs occurring""" + _counts_by_treatment_id = _df \ + .loc[pd.to_datetime(_df['date']).between(*TARGET_PERIOD), 'TREATMENT_ID'] \ + .apply(pd.Series) \ + .sum() \ + .astype(int) + return _counts_by_treatment_id.groupby(level=0).sum() + +year_target = 2023 +def get_counts_of_hsi_by_treatment_id_by_year(_df): + """Get the counts of the short TREATMENT_IDs occurring""" + _counts_by_treatment_id = _df \ + .loc[pd.to_datetime(_df['date']).dt.year ==year_target, 'TREATMENT_ID'] \ + .apply(pd.Series) \ + .sum() \ + .astype(int) + return _counts_by_treatment_id.groupby(level=0).sum() + +def get_counts_of_hsi_by_short_treatment_id(_df): + """Get the counts of the short TREATMENT_IDs occurring (shortened, up to first underscore)""" + _counts_by_treatment_id = get_counts_of_hsi_by_treatment_id(_df) + _short_treatment_id = _counts_by_treatment_id.index.map(lambda x: x.split('_')[0] + "*") + return _counts_by_treatment_id.groupby(by=_short_treatment_id).sum() + +def get_counts_of_hsi_by_short_treatment_id_by_year(_df): + """Get the counts of the short TREATMENT_IDs occurring (shortened, up to first underscore)""" + _counts_by_treatment_id = get_counts_of_hsi_by_treatment_id_by_year(_df) + _short_treatment_id = _counts_by_treatment_id.index.map(lambda x: x.split('_')[0] + "*") + return _counts_by_treatment_id.groupby(by=_short_treatment_id).sum() + + +# Obtain parameter names for this scenario file +param_names = get_parameter_names_from_scenario_file() +print(param_names) + + +def get_counts_of_death_by_year(df): + """Aggregate the model outputs into five-year periods for age and time""" + #_, agegrplookup = make_age_grp_lookup() + #_, calperiodlookup = make_calendar_period_lookup() + df["year"] = df["date"].dt.year + #df["age_grp"] = df["age"].map(agegrplookup).astype(make_age_grp_types()) + #df["period"] = df["year"].map(calperiodlookup).astype(make_calendar_period_type()) + return df.groupby(by=["year", "label"])["person_id"].count() + + +deaths_by_year_and_cause = extract_results( + results_folder, + module="tlo.methods.demography", + key="death", + custom_generate_series=get_counts_of_death_by_year, + do_scaling=True +).pipe(set_param_names_as_column_index_level_0) + + +print(deaths_by_year_and_cause) +print(list(deaths_by_year_and_cause.index)) +deaths_by_year_and_cause.to_csv('ConvertedOutputs/Deaths_by_cause_with_time.csv', index=True) +exit(-1) + +# ================================================================================================ +# TIME EVOLUTION OF TOTAL DALYs +# Plot DALYs averted compared to the ``No Policy'' policy + +year_target = 2023 # This global variable will be passed to custom function +def get_num_dalys_in_year(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year == year_target] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum().sum() + ) + +ALL = {} +# Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred +# are consistent across different policies +this_min_year = 2010 +for year in range(this_min_year, max_year+1): + year_target = year + num_dalys_by_year = extract_results( + results_folder, + module='tlo.methods.healthburden', + key='dalys_stacked', + custom_generate_series=get_num_dalys_in_year, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = num_dalys_by_year + print(num_dalys_by_year.columns) +# Concatenate the DataFrames into a single DataFrame +concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) +concatenated_df.index = concatenated_df.index.set_names(['date', 'index_original']) +concatenated_df = concatenated_df.reset_index(level='index_original',drop=True) +dalys_by_year = concatenated_df +print(dalys_by_year) +dalys_by_year.to_csv('ConvertedOutputs/Total_DALYs_with_time.csv', index=True) + +# ================================================================================================ +# Print population under each scenario +pop_model = extract_results(results_folder, + module="tlo.methods.demography", + key="population", + column="total", + index="date", + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + +pop_model.index = pop_model.index.year +pop_model = pop_model[(pop_model.index >= this_min_year) & (pop_model.index <= max_year)] +print(pop_model) +assert dalys_by_year.index.equals(pop_model.index) +assert all(dalys_by_year.columns == pop_model.columns) +pop_model.to_csv('ConvertedOutputs/Population_with_time.csv', index=True) + +# ================================================================================================ +# DALYs BROKEN DOWN BY CAUSES AND YEAR +# DALYs by cause per year +# %% Quantify the health losses associated with all interventions combined. + +year_target = 2023 # This global variable will be passed to custom function +def get_num_dalys_by_year_and_cause(_df): + """Return total number of DALYs (Stacked) by label (total within the TARGET_PERIOD)""" + return pd.Series( + data=_df + .loc[_df.year == year_target] + .drop(columns=['date', 'sex', 'age_range', 'year']) + .sum() + ) + +ALL = {} +# Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred +# are consistent across different policies +this_min_year = 2010 +for year in range(this_min_year, max_year+1): + year_target = year + num_dalys_by_year = extract_results( + results_folder, + module='tlo.methods.healthburden', + key='dalys_stacked', + custom_generate_series=get_num_dalys_by_year_and_cause, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = num_dalys_by_year #summarize(num_dalys_by_year) + +# Concatenate the DataFrames into a single DataFrame +concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) + +concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) + +df_total = concatenated_df +df_total.to_csv('ConvertedOutputs/DALYS_by_cause_with_time.csv', index=True) + +ALL = {} +# Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred +# are consistent across different policies +for year in range(min_year, max_year+1): + year_target = year + + hsi_delivered_by_year = extract_results( + results_folder, + module='tlo.methods.healthsystem.summary', + key='HSI_Event', + custom_generate_series=get_counts_of_hsi_by_short_treatment_id_by_year, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = hsi_delivered_by_year + +# Concatenate the DataFrames into a single DataFrame +concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) +concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) +HSI_ran_by_year = concatenated_df + +del ALL + +ALL = {} +# Plot time trend show year prior transition as well to emphasise that until that point DALYs incurred +# are consistent across different policies +for year in range(min_year, max_year+1): + year_target = year + + hsi_not_delivered_by_year = extract_results( + results_folder, + module='tlo.methods.healthsystem.summary', + key='Never_ran_HSI_Event', + custom_generate_series=get_counts_of_hsi_by_short_treatment_id_by_year, + do_scaling=True + ).pipe(set_param_names_as_column_index_level_0) + ALL[year_target] = hsi_not_delivered_by_year + +# Concatenate the DataFrames into a single DataFrame +concatenated_df = pd.concat(ALL.values(), keys=ALL.keys()) +concatenated_df.index = concatenated_df.index.set_names(['date', 'cause']) +HSI_never_ran_by_year = concatenated_df + +HSI_never_ran_by_year = HSI_never_ran_by_year.fillna(0) #clean_df( +HSI_ran_by_year = HSI_ran_by_year.fillna(0) +HSI_total_by_year = HSI_ran_by_year.add(HSI_never_ran_by_year, fill_value=0) +HSI_ran_by_year.to_csv('ConvertedOutputs/HSIs_ran_by_area_with_time.csv', index=True) +HSI_never_ran_by_year.to_csv('ConvertedOutputs/HSIs_never_ran_by_area_with_time.csv', index=True) +print(HSI_ran_by_year) +print(HSI_never_ran_by_year) +print(HSI_total_by_year) +exit(-1) + +""" +if __name__ == "__main__": + rfp = Path('resources') + + parser = argparse.ArgumentParser( + description="Produce plots to show the impact each set of treatments", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument( + "--output-path", + help=( + "Directory to write outputs to. If not specified (set to None) outputs " + "will be written to value of --results-path argument." + ), + type=Path, + default=None, + required=False, + ) + parser.add_argument( + "--resources-path", + help="Directory containing resource files", + type=Path, + default=Path('resources'), + required=False, + ) + parser.add_argument( + "--results-path", + type=Path, + help=( + "Directory containing results from running " + "src/scripts/healthsystem/impact_of_const_capabilities_expansion/scenario_impact_of_capabilities_expansion_scaling.py " + ), + default=None, + required=False + ) + args = parser.parse_args() + assert args.results_path is not None + results_path = args.results_path + + output_path = results_path if args.output_path is None else args.output_path + + apply( + results_folder=results_path, + output_folder=output_path, + resourcefilepath=args.resources_path + ) +""" diff --git a/src/scripts/rti_emulator/scenario_test_rti_emulator.py b/src/scripts/rti_emulator/scenario_test_rti_emulator.py new file mode 100644 index 0000000000..c60c682322 --- /dev/null +++ b/src/scripts/rti_emulator/scenario_test_rti_emulator.py @@ -0,0 +1,145 @@ +"""This Scenario file run the model to test the use of the RTI emulator + +Run on the batch system using: +``` +tlo batch-submit + src/scripts/analysis_data_generation/scenario_test_rti_emulator.py +``` + +or locally using: +``` + tlo scenario-run src/scripts/rti_emulator/scenario_test_rti_emulator.py +``` + +""" +from pathlib import Path +from typing import Dict + +import pandas as pd + +from tlo import Date, logging +from tlo.analysis.utils import get_parameters_for_status_quo, mix_scenarios +from tlo.methods.fullmodel import fullmodel +from tlo.methods.scenario_switcher import ImprovedHealthSystemAndCareSeekingScenarioSwitcher +from tlo.scenario import BaseScenario +from tlo.methods import ( + alri, + cardio_metabolic_disorders, + care_of_women_during_pregnancy, + simplified_births, + contraception, + demography, + depression, + diarrhoea, + enhanced_lifestyle, + epi, + healthburden, + healthseekingbehaviour, + healthsystem, + hiv, + rti, + labour, + malaria, + newborn_outcomes, + postnatal_supervisor, + pregnancy_supervisor, + stunting, + symptommanager, + tb, + wasting, +) + +class GenerateDataChains(BaseScenario): + def __init__(self): + super().__init__() + self.seed = 0 + self.start_date = Date(2010, 1, 1) + self.end_date = self.start_date + pd.DateOffset(years=10) + self.pop_size = 50_000 + self._scenarios = self._get_scenarios() + self.number_of_draws = len(self._scenarios) + self.runs_per_draw = 50 + + def log_configuration(self): + return { + 'filename': 'test_rti_emulator', + 'directory': Path('./outputs'), # <- (specified only for local running) + 'custom_levels': { + '*': logging.WARNING, + 'tlo.methods.demography': logging.INFO, + 'tlo.methods.events': logging.INFO, + 'tlo.methods.demography.detail': logging.WARNING, + 'tlo.methods.healthburden': logging.INFO, + 'tlo.methods.healthsystem.summary': logging.INFO, + } + } + + def modules(self): + # MODIFY + # Here instead of running full module + return [demography.Demography(), #resourcefilepath=self.resources), + enhanced_lifestyle.Lifestyle(),#resourcefilepath=self.resources), + healthburden.HealthBurden(),#resourcefilepath=self.resources), + symptommanager.SymptomManager(), #(resourcefilepath=self.resources, spurious_symptoms=False), + rti.RTI(), #resourcefilepath=self.resources), + healthseekingbehaviour.HealthSeekingBehaviour(), #(resourcefilepath=self.resources), + simplified_births.SimplifiedBirths(), #(resourcefilepath=self.resources), + healthsystem.HealthSystem(#resourcefilepath=self.resources, + mode_appt_constraints=1, + cons_availability='all')] + + # return ( + # fullmodel(resourcefilepath=self.resources) + # + [ImprovedHealthSystemAndCareSeekingScenarioSwitcher(resourcefilepath=self.resources)] + # ) + + def draw_parameters(self, draw_number, rng): + if draw_number < self.number_of_draws: + return list(self._scenarios.values())[draw_number] + else: + return + + # case 1: gfHE = -0.030, factor = 1.01074 + # case 2: gfHE = -0.020, factor = 1.02116 + # case 3: gfHE = -0.015, factor = 1.02637 + # case 4: gfHE = 0.015, factor = 1.05763 + # case 5: gfHE = 0.020, factor = 1.06284 + # case 6: gfHE = 0.030, factor = 1.07326 + + def _get_scenarios(self) -> Dict[str, Dict]: + """Return the Dict with values for the parameters that are changed, keyed by a name for the scenario. + """ + + return { + + # =========== STATUS QUO ============ + "Baseline": + mix_scenarios( + self._baseline(), + #{ + # "HealthSystem": { + # "yearly_HR_scaling_mode": "no_scaling", + # }, + #} + ), + + } + + def _baseline(self) -> Dict: + """Return the Dict with values for the parameter changes that define the baseline scenario. """ + return mix_scenarios( + get_parameters_for_status_quo(), + { + "HealthSystem": { + "mode_appt_constraints": 1, + "use_funded_or_actual_staffing": "actual", + "cons_availability": "all", + "Service_Availability": [] + } + }, + ) + +if __name__ == '__main__': + from tlo.cli import scenario_run + + scenario_run([__file__]) diff --git a/src/tlo/methods/emulated_rti.py b/src/tlo/methods/emulated_rti.py new file mode 100644 index 0000000000..4491c047cc --- /dev/null +++ b/src/tlo/methods/emulated_rti.py @@ -0,0 +1,1452 @@ +""" +Road traffic injury module. + +""" +from __future__ import annotations + +from pathlib import Path +from typing import TYPE_CHECKING, List, Optional +from collections import Counter + +import numpy as np +import pandas as pd + +from tlo import DateOffset, Module, Parameter, Property, Types, logging, Date + +from tlo.events import Event, IndividualScopeEventMixin, PopulationScopeEventMixin, RegularEvent +from tlo.lm import LinearModel, LinearModelType, Predictor +from tlo.methods import Metadata +from tlo.methods.causes import Cause +from tlo.methods.hsi_event import HSI_Event +from tlo.methods.hsi_generic_first_appts import GenericFirstAppointmentsMixin +from tlo.methods.symptommanager import Symptom +from tlo.util import read_csv_files + +from sdv.single_table import CTGANSynthesizer +from sdv.single_table import GaussianCopulaSynthesizer +from sdv.single_table import TVAESynthesizer +from sdv.sampling import Condition + +if TYPE_CHECKING: + from tlo.methods.hsi_generic_first_appts import HSIEventScheduler + from tlo.population import IndividualProperties + +# --------------------------------------------------------------------------------------------------------- +# MODULE DEFINITIONS +# --------------------------------------------------------------------------------------------------------- + +logger = logging.getLogger(__name__) +logger.setLevel(logging.DEBUG) + + +class EmulatedRTI(Module, GenericFirstAppointmentsMixin): + """ + The road traffic injuries module for the TLO model, handling all injuries related to road traffic accidents. + """ + + def __init__(self, name=None): + # NB. Parameters passed to the module can be inserted in the __init__ definition. + super().__init__(name) + + INIT_DEPENDENCIES = {"SymptomManager", + "HealthBurden"} + + ADDITIONAL_DEPENDENCIES = { + 'Demography', + 'Lifestyle', + 'HealthSystem', + } + + # ================================================================================ + # EMULATOR PARAMETERS + # Counters tracking use of HealthSystem by RTI module under use of emulator + HS_Use_Type = [ + 'Level2_AccidentsandEmerg', 'Level2_DiagRadio', 'Level2_EPI', + 'Level2_IPAdmission', 'Level2_InpatientDays', 'Level2_MajorSurg', + 'Level2_MinorSurg', 'Level2_Over5OPD', 'Level2_Tomography', 'Level2_Under5OPD' + ] + + # Initialize the counter with all items set to 0 + HS_Use_by_RTI = Counter({col: 0 for col in HS_Use_Type}) + + Rti_Services = ['Rti_AcutePainManagement','Rti_BurnManagement','Rti_FractureCast','Rti_Imaging','Rti_MajorSurgeries','Rti_MedicalIntervention','Rti_MinorSurgeries','Rti_OpenFractureTreatment','Rti_ShockTreatment','Rti_Suture','Rti_TetanusVaccine'] + + HS_conditions = {} + + RTI_emulator = None + + # Module parameters + PARAMETERS = { + 'main_polling_frequency': Parameter( + Types.INT, + 'Frequency in months for RTI polling events that determine new injuries' + ), + 'base_rate_injrti': Parameter( + Types.REAL, + 'Base rate of RTI per year', + ), + 'rr_injrti_age04': Parameter( + Types.REAL, + 'risk ratio of RTI in age 0-4 compared to base rate of RTI' + ), + 'rr_injrti_age59': Parameter( + Types.REAL, + 'risk ratio of RTI in age 5-9 compared to base rate of RTI' + ), + 'rr_injrti_age1017': Parameter( + Types.REAL, + 'risk ratio of RTI in age 10-17 compared to base rate of RTI' + ), + 'rr_injrti_age1829': Parameter( + Types.REAL, + 'risk ratio of RTI in age 18-29 compared to base rate of RTI', + ), + 'rr_injrti_age3039': Parameter( + Types.REAL, + 'risk ratio of RTI in age 30-39 compared to base rate of RTI', + ), + 'rr_injrti_age4049': Parameter( + Types.REAL, + 'risk ratio of RTI in age 40-49 compared to base rate of RTI', + ), + 'rr_injrti_age5059': Parameter( + Types.REAL, + 'risk ratio of RTI in age 50-59 compared to base rate of RTI', + ), + 'rr_injrti_age6069': Parameter( + Types.REAL, + 'risk ratio of RTI in age 60-69 compared to base rate of RTI', + ), + 'rr_injrti_age7079': Parameter( + Types.REAL, + 'risk ratio of RTI in age 70-79 compared to base rate of RTI', + ), + 'rr_injrti_male': Parameter( + Types.REAL, + 'risk ratio of RTI when male compared to females', + ), + 'rr_injrti_excessalcohol': Parameter( + Types.REAL, + 'risk ratio of RTI in those that consume excess alcohol compared to those who do not' + ), + 'imm_death_proportion_rti': Parameter( + Types.REAL, + 'Proportion of those involved in an RTI that die at site of accident or die before seeking medical ' + 'intervention' + ), + 'prob_bleeding_leads_to_shock': Parameter( + Types.REAL, + 'The proportion of those with heavily bleeding injuries who go into shock' + ), + 'prob_death_iss_less_than_9': Parameter( + Types.REAL, + 'Proportion of people who pass away in the following month after medical treatment for injuries with an ISS' + 'score less than or equal to 9' + ), + 'prob_death_iss_10_15': Parameter( + Types.REAL, + 'Proportion of people who pass away in the following month after medical treatment for injuries with an ISS' + 'score from 10 to 15' + ), + 'prob_death_iss_16_24': Parameter( + Types.REAL, + 'Proportion of people who pass away in the following month after medical treatment for injuries with an ISS' + 'score from 16 to 24' + ), + 'prob_death_iss_25_35': Parameter( + Types.REAL, + 'Proportion of people who pass away in the following month after medical treatment for injuries with an ISS' + 'score from 25 to 34' + ), + 'prob_death_iss_35_plus': Parameter( + Types.REAL, + 'Proportion of people who pass away in the following month after medical treatment for injuries with an ISS' + 'score 35 and above' + ), + 'prob_perm_disability_with_treatment_severe_TBI': Parameter( + Types.REAL, + 'probability that someone with a treated severe TBI is permanently disabled' + ), + 'prob_death_TBI_SCI_no_treatment': Parameter( + Types.REAL, + 'probability that someone with a spinal cord injury will die without treatment' + ), + 'prop_death_burns_no_treatment': Parameter( + Types.REAL, + 'probability that someone with a burn injury will die without treatment' + ), + 'prob_death_fractures_no_treatment': Parameter( + Types.REAL, + 'probability that someone with a fracture injury will die without treatment' + ), + 'prob_TBI_require_craniotomy': Parameter( + Types.REAL, + 'probability that someone with a traumatic brain injury will require a craniotomy surgery' + ), + 'prob_exploratory_laparotomy': Parameter( + Types.REAL, + 'probability that someone with an internal organ injury will require a exploratory_laparotomy' + ), + 'prob_depressed_skull_fracture': Parameter( + Types.REAL, + 'Probability that a skull fracture will be depressed and therefore require surgery' + ), + 'prob_mild_burns': Parameter( + Types.REAL, + 'Probability that a burn within a region will result in < 10% total body surface area' + ), + 'prob_dislocation_requires_surgery': Parameter( + Types.REAL, + 'Probability that a dislocation will require surgery to relocate the joint.' + ), + 'number_of_injured_body_regions_distribution': Parameter( + Types.LIST, + 'The distribution of number of injured AIS body regions, used to decide how many injuries a person has' + ), + 'injury_location_distribution': Parameter( + Types.LIST, + 'The distribution of where injuries are located in the body, based on the AIS body region definition' + ), + # Length of stay + 'mean_los_ISS_less_than_4': Parameter( + Types.REAL, + 'Mean length of stay for someone with an ISS score < 4' + ), + 'sd_los_ISS_less_than_4': Parameter( + Types.REAL, + 'Standard deviation in length of stay for someone with an ISS score < 4' + ), + 'mean_los_ISS_4_to_8': Parameter( + Types.REAL, + 'Mean length of stay for someone with an ISS score between 4 and 8' + ), + 'sd_los_ISS_4_to_8': Parameter( + Types.REAL, + 'Standard deviation in length of stay for someone with an ISS score between 4 and 8' + ), + 'mean_los_ISS_9_to_15': Parameter( + Types.REAL, + 'Mean length of stay for someone with an ISS score between 9 and 15' + ), + 'sd_los_ISS_9_to_15': Parameter( + Types.REAL, + 'Standard deviation in length of stay for someone with an ISS score between 9 and 15' + ), + 'mean_los_ISS_16_to_24': Parameter( + Types.REAL, + 'Mean length of stay for someone with an ISS score between 16 and 24' + ), + 'sd_los_ISS_16_to_24': Parameter( + Types.REAL, + 'Standard deviation in length of stay for someone with an ISS score between 16 and 24' + ), + 'mean_los_ISS_more_than_25': Parameter( + Types.REAL, + 'Mean length of stay for someone with an ISS score between 16 and 24' + ), + 'sd_los_ISS_more_that_25': Parameter( + Types.REAL, + 'Standard deviation in length of stay for someone with an ISS score between 16 and 24' + ), + # DALY weights + 'daly_wt_unspecified_skull_fracture': Parameter( + Types.REAL, + 'daly_wt_unspecified_skull_fracture - code 1674' + ), + 'daly_wt_basilar_skull_fracture': Parameter( + Types.REAL, + 'daly_wt_basilar_skull_fracture - code 1675' + ), + 'daly_wt_epidural_hematoma': Parameter( + Types.REAL, + 'daly_wt_epidural_hematoma - code 1676' + ), + 'daly_wt_subdural_hematoma': Parameter( + Types.REAL, + 'daly_wt_subdural_hematoma - code 1677' + ), + 'daly_wt_subarachnoid_hematoma': Parameter( + Types.REAL, + 'daly_wt_subarachnoid_hematoma - code 1678' + ), + 'daly_wt_brain_contusion': Parameter( + Types.REAL, + 'daly_wt_brain_contusion - code 1679' + ), + 'daly_wt_intraventricular_haemorrhage': Parameter( + Types.REAL, + 'daly_wt_intraventricular_haemorrhage - code 1680' + ), + 'daly_wt_diffuse_axonal_injury': Parameter( + Types.REAL, + 'daly_wt_diffuse_axonal_injury - code 1681' + ), + 'daly_wt_subgaleal_hematoma': Parameter( + Types.REAL, + 'daly_wt_subgaleal_hematoma - code 1682' + ), + 'daly_wt_midline_shift': Parameter( + Types.REAL, + 'daly_wt_midline_shift - code 1683' + ), + 'daly_wt_facial_fracture': Parameter( + Types.REAL, + 'daly_wt_facial_fracture - code 1684' + ), + 'daly_wt_facial_soft_tissue_injury': Parameter( + Types.REAL, + 'daly_wt_facial_soft_tissue_injury - code 1685' + ), + 'daly_wt_eye_injury': Parameter( + Types.REAL, + 'daly_wt_eye_injury - code 1686' + ), + 'daly_wt_neck_soft_tissue_injury': Parameter( + Types.REAL, + 'daly_wt_neck_soft_tissue_injury - code 1687' + ), + 'daly_wt_neck_internal_bleeding': Parameter( + Types.REAL, + 'daly_wt_neck_internal_bleeding - code 1688' + ), + 'daly_wt_neck_dislocation': Parameter( + Types.REAL, + 'daly_wt_neck_dislocation - code 1689' + ), + 'daly_wt_chest_wall_bruises_hematoma': Parameter( + Types.REAL, + 'daly_wt_chest_wall_bruises_hematoma - code 1690' + ), + 'daly_wt_hemothorax': Parameter( + Types.REAL, + 'daly_wt_hemothorax - code 1691' + ), + 'daly_wt_lung_contusion': Parameter( + Types.REAL, + 'daly_wt_lung_contusion - code 1692' + ), + 'daly_wt_diaphragm_rupture': Parameter( + Types.REAL, + 'daly_wt_diaphragm_rupture - code 1693' + ), + 'daly_wt_rib_fracture': Parameter( + Types.REAL, + 'daly_wt_rib_fracture - code 1694' + ), + 'daly_wt_flail_chest': Parameter( + Types.REAL, + 'daly_wt_flail_chest - code 1695' + ), + 'daly_wt_chest_wall_laceration': Parameter( + Types.REAL, + 'daly_wt_chest_wall_laceration - code 1696' + ), + 'daly_wt_closed_pneumothorax': Parameter( + Types.REAL, + 'daly_wt_closed_pneumothorax - code 1697' + ), + 'daly_wt_open_pneumothorax': Parameter( + Types.REAL, + 'daly_wt_open_pneumothorax - code 1698' + ), + 'daly_wt_surgical_emphysema': Parameter( + Types.REAL, + 'daly_wt_surgical_emphysema aka subcuteal emphysema - code 1699' + ), + 'daly_wt_abd_internal_organ_injury': Parameter( + Types.REAL, + 'daly_wt_abd_internal_organ_injury - code 1700' + ), + 'daly_wt_spinal_cord_lesion_neck_with_treatment': Parameter( + Types.REAL, + 'daly_wt_spinal_cord_lesion_neck_with_treatment - code 1701' + ), + 'daly_wt_spinal_cord_lesion_neck_without_treatment': Parameter( + Types.REAL, + 'daly_wt_spinal_cord_lesion_neck_without_treatment - code 1702' + ), + 'daly_wt_spinal_cord_lesion_below_neck_with_treatment': Parameter( + Types.REAL, + 'daly_wt_spinal_cord_lesion_below_neck_with_treatment - code 1703' + ), + 'daly_wt_spinal_cord_lesion_below_neck_without_treatment': Parameter( + Types.REAL, + 'daly_wt_spinal_cord_lesion_below_neck_without_treatment - code 1704' + ), + 'daly_wt_vertebrae_fracture': Parameter( + Types.REAL, + 'daly_wt_vertebrae_fracture - code 1705' + ), + 'daly_wt_clavicle_scapula_humerus_fracture': Parameter( + Types.REAL, + 'daly_wt_clavicle_scapula_humerus_fracture - code 1706' + ), + 'daly_wt_hand_wrist_fracture_with_treatment': Parameter( + Types.REAL, + 'daly_wt_hand_wrist_fracture_with_treatment - code 1707' + ), + 'daly_wt_hand_wrist_fracture_without_treatment': Parameter( + Types.REAL, + 'daly_wt_hand_wrist_fracture_without_treatment - code 1708' + ), + 'daly_wt_radius_ulna_fracture_short_term_with_without_treatment': Parameter( + Types.REAL, + 'daly_wt_radius_ulna_fracture_short_term_with_without_treatment - code 1709' + ), + 'daly_wt_radius_ulna_fracture_long_term_without_treatment': Parameter( + Types.REAL, + 'daly_wt_radius_ulna_fracture_long_term_without_treatment - code 1710' + ), + 'daly_wt_dislocated_shoulder': Parameter( + Types.REAL, + 'daly_wt_dislocated_shoulder - code 1711' + ), + 'daly_wt_amputated_finger': Parameter( + Types.REAL, + 'daly_wt_amputated_finger - code 1712' + ), + 'daly_wt_amputated_thumb': Parameter( + Types.REAL, + 'daly_wt_amputated_thumb - code 1713' + ), + 'daly_wt_unilateral_arm_amputation_with_treatment': Parameter( + Types.REAL, + 'daly_wt_unilateral_arm_amputation_with_treatment - code 1714' + ), + 'daly_wt_unilateral_arm_amputation_without_treatment': Parameter( + Types.REAL, + 'daly_wt_unilateral_arm_amputation_without_treatment - code 1715' + ), + 'daly_wt_bilateral_arm_amputation_with_treatment': Parameter( + Types.REAL, + 'daly_wt_bilateral_arm_amputation_with_treatment - code 1716' + ), + 'daly_wt_bilateral_arm_amputation_without_treatment': Parameter( + Types.REAL, + 'daly_wt_bilateral_arm_amputation_without_treatment - code 1717' + ), + 'daly_wt_foot_fracture_short_term_with_without_treatment': Parameter( + Types.REAL, + 'daly_wt_foot_fracture_short_term_with_without_treatment - code 1718' + ), + 'daly_wt_foot_fracture_long_term_without_treatment': Parameter( + Types.REAL, + 'daly_wt_foot_fracture_long_term_without_treatment - code 1719' + ), + 'daly_wt_patella_tibia_fibula_fracture_with_treatment': Parameter( + Types.REAL, + 'daly_wt_patella_tibia_fibula_fracture_with_treatment - code 1720' + ), + 'daly_wt_patella_tibia_fibula_fracture_without_treatment': Parameter( + Types.REAL, + 'daly_wt_patella_tibia_fibula_fracture_without_treatment - code 1721' + ), + 'daly_wt_hip_fracture_short_term_with_without_treatment': Parameter( + Types.REAL, + 'daly_wt_hip_fracture_short_term_with_without_treatment - code 1722' + ), + 'daly_wt_hip_fracture_long_term_with_treatment': Parameter( + Types.REAL, + 'daly_wt_hip_fracture_long_term_with_treatment - code 1723' + ), + 'daly_wt_hip_fracture_long_term_without_treatment': Parameter( + Types.REAL, + 'daly_wt_hip_fracture_long_term_without_treatment - code 1724' + ), + 'daly_wt_pelvis_fracture_short_term': Parameter( + Types.REAL, + 'daly_wt_pelvis_fracture_short_term - code 1725' + ), + 'daly_wt_pelvis_fracture_long_term': Parameter( + Types.REAL, + 'daly_wt_pelvis_fracture_long_term - code 1726' + ), + 'daly_wt_femur_fracture_short_term': Parameter( + Types.REAL, + 'daly_wt_femur_fracture_short_term - code 1727' + ), + 'daly_wt_femur_fracture_long_term_without_treatment': Parameter( + Types.REAL, + 'daly_wt_femur_fracture_long_term_without_treatment - code 1728' + ), + 'daly_wt_dislocated_hip': Parameter( + Types.REAL, + 'daly_wt_dislocated_hip - code 1729' + ), + 'daly_wt_dislocated_knee': Parameter( + Types.REAL, + 'daly_wt_dislocated_knee - code 1730' + ), + 'daly_wt_amputated_toes': Parameter( + Types.REAL, + 'daly_wt_amputated_toes - code 1731' + ), + 'daly_wt_unilateral_lower_limb_amputation_with_treatment': Parameter( + Types.REAL, + 'daly_wt_unilateral_lower_limb_amputation_with_treatment - code 1732' + ), + 'daly_wt_unilateral_lower_limb_amputation_without_treatment': Parameter( + Types.REAL, + 'daly_wt_unilateral_lower_limb_amputation_without_treatment - code 1733' + ), + 'daly_wt_bilateral_lower_limb_amputation_with_treatment': Parameter( + Types.REAL, + 'daly_wt_bilateral_lower_limb_amputation_with_treatment - code 1734' + ), + 'daly_wt_bilateral_lower_limb_amputation_without_treatment': Parameter( + Types.REAL, + 'daly_wt_bilateral_lower_limb_amputation_without_treatment - code 1735' + ), + 'rt_emergency_care_ISS_score_cut_off': Parameter( + Types.INT, + 'A parameter to determine which level of injury severity corresponds to the emergency health care seeking ' + 'symptom and which to the non-emergency generic injury symptom' + ), + 'prob_death_non_serious': Parameter( + Types.REAL, + 'A parameter to determine the probability of death for non serious condition' + ), 'prob_death_MAIS1': Parameter( + Types.REAL, + 'A parameter to determine the probability of death without medical intervention with a military AIS' + 'score of 1' + ), + 'prob_death_MAIS2': Parameter( + Types.REAL, + 'A parameter to determine the probability of death without medical intervention with a military AIS' + 'score of 2' + ), + 'prob_death_MAIS3': Parameter( + Types.REAL, + 'A parameter to determine the probability of death without medical intervention with a military AIS' + 'score of 3' + ), + 'prob_death_MAIS4': Parameter( + Types.REAL, + 'A parameter to determine the probability of death without medical intervention with a military AIS' + 'score of 4' + ), + 'prob_death_MAIS5': Parameter( + Types.REAL, + 'A parameter to determine the probability of death without medical intervention with a military AIS' + 'score of 5' + ), + 'prob_death_MAIS6': Parameter( + Types.REAL, + 'A parameter to determine the probability of death without medical intervention with a military AIS' + 'score of 6' + ), + 'femur_fracture_skeletal_traction_mean_los': Parameter( + Types.INT, + 'The mean length of stay for a person with a femur fracture being treated with skeletal traction' + ), + 'other_skeletal_traction_los': Parameter( + Types.INT, + 'The mean length of stay for a person with a non-femur fracture being treated with skeletal traction' + ), + 'prob_foot_frac_require_cast': Parameter( + Types.REAL, + 'The probability that a person with a foot fracture will be treated with a plaster cast' + ), + 'prob_foot_frac_require_maj_surg': Parameter( + Types.REAL, + 'The probability that a person with a foot fracture will be treated with a major surgery' + ), + 'prob_foot_frac_require_min_surg': Parameter( + Types.REAL, + 'The probability that a person with a foot fracture will be treated with a major surgery' + ), + 'prob_foot_frac_require_amp': Parameter( + Types.REAL, + 'The probability that a person with a foot fracture will be treated with amputation via a major surgery' + ), + 'prob_tib_fib_frac_require_cast': Parameter( + Types.REAL, + 'The probability that a person with a tibia/fibula fracture will be treated with a plaster cast' + ), + 'prob_tib_fib_frac_require_maj_surg': Parameter( + Types.REAL, + 'The probability that a person with a tibia/fibula fracture will be treated with a major surgery' + ), + 'prob_tib_fib_frac_require_min_surg': Parameter( + Types.REAL, + 'The probability that a person with a tibia/fibula fracture will be treated with a minor surgery' + ), + 'prob_tib_fib_frac_require_amp': Parameter( + Types.REAL, + 'The probability that a person with a tibia/fibula fracture will be treated with an amputation via major ' + 'surgery' + ), + 'prob_tib_fib_frac_require_traction': Parameter( + Types.REAL, + 'The probability that a person with a tibia/fibula fracture will be treated with skeletal traction' + ), + 'prob_femural_fracture_require_major_surgery': Parameter( + Types.REAL, + 'The probability that a person with a femur fracture will be treated with major surgery' + ), + 'prob_femural_fracture_require_minor_surgery': Parameter( + Types.REAL, + 'The probability that a person with a femur fracture will be treated with minor surgery' + ), + 'prob_femural_fracture_require_cast': Parameter( + Types.REAL, + 'The probability that a person with a femur fracture will be treated with a plaster cast' + ), + 'prob_femural_fracture_require_amputation': Parameter( + Types.REAL, + 'The probability that a person with a femur fracture will be treated with amputation via major surgery' + ), + 'prob_femural_fracture_require_traction': Parameter( + Types.REAL, + 'The probability that a person with a femur fracture will be treated with skeletal traction' + ), + 'prob_pelvis_fracture_traction': Parameter( + Types.REAL, + 'The probability that a person with a pelvis fracture will be treated with skeletal traction' + ), + 'prob_pelvis_frac_major_surgery': Parameter( + Types.REAL, + 'The probability that a person with a pelvis fracture will be treated with major surgery' + ), + 'prob_pelvis_frac_minor_surgery': Parameter( + Types.REAL, + 'The probability that a person with a pelvis fracture will be treated with minor surgery' + ), + 'prob_pelvis_frac_cast': Parameter( + Types.REAL, + 'The probability that a person with a pelvis fracture will be treated with a cast' + ), + 'prob_dis_hip_require_maj_surg': Parameter( + Types.REAL, + 'The probability that a person with a dislocated hip will be treated with a major surgery' + ), + 'prob_dis_hip_require_cast': Parameter( + Types.REAL, + 'The probability that a person with a dislocated hip will be treated with a plaster cast' + ), + 'prob_hip_dis_require_traction': Parameter( + Types.REAL, + 'The probability that a person with a dislocated hip will be treated with skeletal traction' + ), + 'hdu_cut_off_iss_score': Parameter( + Types.INT, + 'The ISS score used as a criteria to admit patients to the HDU/ICU units' + ), + 'mean_icu_days': Parameter( + Types.REAL, + 'The mean length of stay in the ICUfor those without TBI' + ), + 'sd_icu_days': Parameter( + Types.REAL, + 'The standard deviation in length of stay in the ICU for those without TBI' + ), + 'mean_tbi_icu_days': Parameter( + Types.REAL, + 'The mean length of stay in the ICU for those with TBI' + ), + 'sd_tbi_icu_days': Parameter( + Types.REAL, + 'The standard deviation in length of stay in the ICU for those with TBI' + ), + 'prob_open_fracture_contaminated': Parameter( + Types.REAL, + 'The probability that an open fracture will be contaminated' + ), + 'allowed_interventions': Parameter( + Types.LIST, + 'List of additional interventions that can be included when performing model analysis' + ), + 'head_prob_112': Parameter( + Types.REAL, + "The probability that this person's head injury is a skull fracture" + ), + 'head_prob_113': Parameter( + Types.REAL, + "The probability that this person's head injury is a basilar skull fracture" + ), + 'head_prob_133a': Parameter( + Types.REAL, + "The probability that this person's head injury is a Subarachnoid hematoma" + ), + 'head_prob_133b': Parameter( + Types.REAL, + "The probability that this person's head injury is a Brain contusion" + ), + 'head_prob_133c': Parameter( + Types.REAL, + "The probability that this person's head injury is an Intraventricular haemorrhage" + ), + 'head_prob_133d': Parameter( + Types.REAL, + "The probability that this person's head injury is a Subgaleal hematoma" + ), + 'head_prob_134a': Parameter( + Types.REAL, + "The probability that this person's head injury is an Epidural hematoma" + ), + 'head_prob_134b': Parameter( + Types.REAL, + "The probability that this person's head injury is a Subdural hematoma" + ), + 'head_prob_135': Parameter( + Types.REAL, + "The probability that this person's head injury is a Diffuse axonal injury/midline shift" + ), + 'head_prob_1101': Parameter( + Types.REAL, + "The probability that this person's head injury is a laceration" + ), + 'head_prob_1114': Parameter( + Types.REAL, + "The probability that this person's head injury is a burn" + ), + 'face_prob_211': Parameter( + Types.REAL, + "The probability that this person's face injury is a Facial fracture (nasal/unspecified)" + ), + 'face_prob_212': Parameter( + Types.REAL, + "The probability that this person's face injury is a Facial fracture (mandible/zygomatic)" + ), + 'face_prob_241': Parameter( + Types.REAL, + "The probability that this person's face injury is a soft tissue injury" + ), + 'face_prob_2101': Parameter( + Types.REAL, + "The probability that this person's face injury is a laceration" + ), + 'face_prob_2114': Parameter( + Types.REAL, + "The probability that this person's face injury is a burn" + ), + 'face_prob_291': Parameter( + Types.REAL, + "The probability that this person's face injury is an eye injury" + ), + 'neck_prob_3101': Parameter( + Types.REAL, + "The probability that this person's neck injury is a laceration" + ), + 'neck_prob_3113': Parameter( + Types.REAL, + "The probability that this person's neck injury is a burn" + ), + 'neck_prob_342': Parameter( + Types.REAL, + "The probability that this person's neck injury is a Soft tissue injury in neck (vertebral artery " + "laceration)" + ), + 'neck_prob_343': Parameter( + Types.REAL, + "The probability that this person's neck injury is a Soft tissue injury in neck (pharynx contusion)" + ), + 'neck_prob_361': Parameter( + Types.REAL, + "The probability that this person's neck injury is a Sternomastoid m. hemorrhage/ Hemorrhage, " + "supraclavicular triangle/Hemorrhage, posterior triangle/Anterior vertebral vessel hemorrhage/ Neck muscle " + "hemorrhage" + ), + 'neck_prob_363': Parameter( + Types.REAL, + "The probability that this person's neck injury is a Hematoma in carotid sheath/Carotid sheath hemorrhage" + ), + 'neck_prob_322': Parameter( + Types.REAL, + "The probability that this person's neck injury is an Atlanto-occipital subluxation" + ), + 'neck_prob_323': Parameter( + Types.REAL, + "The probability that this person's neck injury is an Atlanto-axial subluxation" + ), + 'thorax_prob_4101': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a laceration" + ), + 'thorax_prob_4113': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a burn" + ), + 'thorax_prob_461': Parameter( + Types.REAL, + "The probability that this person's thorax injury is Chest wall bruises/haematoma" + ), + 'thorax_prob_463': Parameter( + Types.REAL, + "The probability that this person's thorax injury is Haemothorax" + ), + 'thorax_prob_453a': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a Lung contusion" + ), + 'thorax_prob_453b': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a Diaphragm rupture" + ), + 'thorax_prob_412': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a rib fracture" + ), + 'thorax_prob_414': Parameter( + Types.REAL, + "The probability that this person's thorax injury is flail chest" + ), + 'thorax_prob_441': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a Chest wall lacerations/avulsions" + ), + 'thorax_prob_442': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a Surgical emphysema" + ), + 'thorax_prob_443': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a Closed pneumothorax/ open pneumothorax" + ), + 'abdomen_prob_5101': Parameter( + Types.REAL, + "The probability that this person's abdomen injury is a laceration" + ), + 'abdomen_prob_5113': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a burn" + ), + 'abdomen_prob_552': Parameter( + Types.REAL, + "The probability that this person's thorax injury is a skull fracture" + ), + 'abdomen_prob_553': Parameter( + Types.REAL, + "The probability that this person's thorax injury is an Injury to stomach/intestines/colon" + ), + 'abdomen_prob_554': Parameter( + Types.REAL, + "The probability that this person's thorax injury is an Injury to spleen/Urinary bladder/Liver/Urethra/" + "Diaphragm" + ), + 'spine_prob_612': Parameter( + Types.REAL, + "The probability that this person's spine injury is a vertabrae fracture" + ), + 'spine_prob_673a': Parameter( + Types.REAL, + "The probability that this person's spine injury is a Spinal cord injury at neck level" + ), + 'spine_prob_673b': Parameter( + Types.REAL, + "The probability that this person's spine injury is a Spinal cord injury below neck level" + ), + 'spine_prob_674a': Parameter( + Types.REAL, + "The probability that this person's spine injury is a Spinal cord injury at neck level" + ), + 'spine_prob_674b': Parameter( + Types.REAL, + "The probability that this person's spine injury is a Spinal cord injury below neck level" + ), + 'spine_prob_675a': Parameter( + Types.REAL, + "The probability that this person's spine injury is a Spinal cord injury at neck level" + ), + 'spine_prob_675b': Parameter( + Types.REAL, + "The probability that this person's spine injury is a Spinal cord injury below neck level" + ), + 'spine_prob_676': Parameter( + Types.REAL, + "The probability that this person's spine injury is a Spinal cord injury at neck level" + ), + 'upper_ex_prob_7101': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is a laceration" + ), + 'upper_ex_prob_7113': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is a burn" + ), + 'upper_ex_prob_712a': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is a Fracture to Clavicle, scapula, humerus" + ), + 'upper_ex_prob_712b': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is a Fracture to Hand/wrist" + ), + 'upper_ex_prob_712c': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is a Fracture to Radius/ulna" + ), + 'upper_ex_prob_722': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is a dislocated shoulder" + ), + 'upper_ex_prob_782a': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is an Amputated finger" + ), + 'upper_ex_prob_782b': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is a Unilateral arm amputation" + ), + 'upper_ex_prob_782c': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is a Thumb amputation" + ), + 'upper_ex_prob_783': Parameter( + Types.REAL, + "The probability that this person's upper extremity injury is a bilateral arm amputation" + ), + 'lower_ex_prob_8101': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a laceration" + ), + 'lower_ex_prob_8113': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a burn" + ), + 'lower_ex_prob_811': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a foot fracture" + ), + 'lower_ex_prob_813do': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is an open foot fracture" + ), + 'lower_ex_prob_812': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a Fracture to patella, tibia, fibula, ankle" + ), + 'lower_ex_prob_813eo': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is an open Fracture to patella, tibia, fibula, " + "ankle" + ), + 'lower_ex_prob_813a': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a Hip fracture" + ), + 'lower_ex_prob_813b': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a Pelvis fracture" + ), + 'lower_ex_prob_813bo': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is an open Pelvis fracture" + ), + 'lower_ex_prob_813c': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a Femur fracture" + ), + 'lower_ex_prob_813co': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is an open Femur fracture" + ), + 'lower_ex_prob_822a': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a Dislocated hip" + ), + 'lower_ex_prob_822b': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a Dislocated knee" + ), + 'lower_ex_prob_882': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is an Amputation of toes" + ), + 'lower_ex_prob_883': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a Unilateral leg amputation" + ), + 'lower_ex_prob_884': Parameter( + Types.REAL, + "The probability that this person's lower extremity injury is a Bilateral leg amputation " + ), + 'blocked_interventions': Parameter( + Types.LIST, + "A list of interventions that are blocked in a simulation" + ), + 'unavailable_treatment_mortality_mais_cutoff': Parameter( + Types.INT, + "A cut-off score above which an injury will result in additional mortality if the person has " + "sought healthcare and not received it." + ), + 'consider_death_no_treatment_ISS_cut_off': Parameter( + Types.INT, + "A cut-off score above which an injuries will be considered severe enough to cause mortality in those who" + "have not sought care." + ), + 'maximum_number_of_times_HSI_events_should_run': Parameter( + Types.INT, + "limit on the number of times an HSI event can run" + ), + 'use_RTI_emulator': Parameter( + Types.BOOL, + "Replace module with RTI emulator, valid if running in mode 1 with actual consumable availability" + ), + 'hsi_schedule_window_days': Parameter( + Types.INT, + 'Number of days window to schedule HSI appointment' + ), + + 'rti_check_death_no_med_event_frequency_days': Parameter( + Types.INT, + 'Frequency in days for RTI check death no med event' + ), + 'rti_recovery_event_frequency_days': Parameter( + Types.INT, + 'Frequency in days for RTI recovery event' + ), + 'incidence_rate_frequency': Parameter( + Types.INT, + 'Number of months per year for rate conversion calculations' + ), + 'incidence_rate_per_population': Parameter( + Types.INT, + 'Population base for incidence rate calculations' + ), + 'days_to_death_without_treatment': Parameter( + Types.INT, + 'Number of days until death for untreated RTI patients' + ), + 'max_treatment_duration_days': Parameter( + Types.INT, + 'Maximum number of days for RTI treatment duration' + ), + 'intervention_incidence_reduction_factor': Parameter( + Types.REAL, + 'Factor by which interventions reduce RTI incidence ' + '(applied when reduce_incidence in allowed_interventions)' + ), + 'laceration_recovery_days': Parameter( + Types.INT, + 'Number of days for recovery assessment period' + ), + 'hsi_opening_delay_days': Parameter( + Types.INT, + 'Number of days delay before HSI appointments can be scheduled' + ), + 'main_polling_initialisation_delay_months': Parameter( + Types.INT, + 'Delay in months at initialisation for first main polling event' + ), + 'rti_recovery_initialisation_delay_months': Parameter( + Types.INT, + 'Delay in months at initialisation for rti recovery event' + ), + 'rti_check_death_no_med_initialisation_delay_months': Parameter( + Types.INT, + 'Delay in months at initialisation for rti check death no med event' + ), + 'non_permanent_tbi_recovery_months': Parameter( + Types.INT, + 'Recovery duration in months for tbi if non permanent; sets date for date_to_remove_daly_column' + ), + 'rti_fracture_cast_recovery_weeks': Parameter( + Types.INT, + 'Recovery duration in weeks for fracture cast; sets date for date_to_remove_daly_column' + ), + 'rti_open_fracture_recovery_months': Parameter( + Types.INT, + 'Recovery duration in months for open fracture; sets date for date_to_remove_daly_column' + ), + 'rti_burn_recovery_weeks': Parameter( + Types.INT, + 'Recovery duration in weeks for rti burn; sets date for date_to_remove_daly_column' + ), + + } + + PROPERTIES = { + 'rt_disability': Property(Types.REAL, 'disability weight for current month'), + 'rt_disability_permanent': Property(Types.REAL, 'disability weight incurred permanently'), + 'rt_date_inj': Property(Types.DATE, 'date of latest injury'), + 'rt_road_traffic_inc': Property(Types.BOOL, 'involved in a road traffic injury'), + } + + # Declare Metadata + METADATA = { + Metadata.DISEASE_MODULE, # Disease modules: Any disease module should carry this label. + Metadata.USES_SYMPTOMMANAGER, # The 'Symptom Manager' recognises modules with this label. + Metadata.USES_HEALTHSYSTEM, # The 'HealthSystem' recognises modules with this label. + Metadata.USES_HEALTHBURDEN # The 'HealthBurden' module recognises modules with this label. + } + + # Declare Causes of Death + CAUSES_OF_DEATH = { + 'RTI_death_without_med': Cause(gbd_causes='Road injuries', label='Transport Injuries'), + 'RTI_death_with_med': Cause(gbd_causes='Road injuries', label='Transport Injuries'), + 'RTI_unavailable_med': Cause(gbd_causes='Road injuries', label='Transport Injuries'), + 'RTI_imm_death': Cause(gbd_causes='Road injuries', label='Transport Injuries'), + 'RTI_death_shock': Cause(gbd_causes='Road injuries', label='Transport Injuries'), + } + + # Declare Causes of Death and Disability + CAUSES_OF_DISABILITY = { + 'RTI': Cause(gbd_causes='Road injuries', label='Transport Injuries') + } + + def read_parameters(self, resourcefilepath: Optional[Path] = None): + """ Reads the parameters used in the RTI module""" + p = self.parameters + + dfd = read_csv_files(resourcefilepath / 'ResourceFile_RTI', files='parameter_values') + self.load_parameters_from_dataframe(dfd) + + # Load emulator + self.RTI_emulator = CTGANSynthesizer.load(filepath= resourcefilepath / 'ResourceFile_RTI/RTI_emulator.pkl') + + def initialise_population(self, population): + """Sets up the default properties used in the RTI module and applies them to the dataframe. The default state + for the RTI module is that people haven't been involved in a road traffic accident and are therefor alive and + healthy.""" + df = population.props + df.loc[df.is_alive, 'rt_disability'] = 0 # default: no DALY + df.loc[df.is_alive, 'rt_disability_permanent'] = 0 # default: no DALY + df.loc[df.is_alive, 'rt_road_traffic_inc'] = False + df.loc[df.is_alive, 'rt_date_inj'] = pd.NaT + + def initialise_simulation(self, sim): + """At the start of the simulation we schedule a logging event, which records the relevant information + regarding road traffic injuries in the last month. + + Afterwards, we schedule three RTI events, the first is the main RTI event which takes parts + of the population and assigns them to be involved in road traffic injuries and providing they survived will + begin the interaction with the healthcare system. This event runs monthly. + + The second is the begin scheduling the RTI recovery event, which looks at those in the population who have been + injured in a road traffic accident, checking every day whether enough time has passed for their injuries to have + healed. When the injury has healed the associated daly weight is removed. + + The final event is one which checks if this person has not sought sought care or been given care, if they + haven't then it asks whether they should die away from their injuries + """ + p = self.parameters + # Begin modelling road traffic injuries + + sim.schedule_event(RTIPollingEvent(self), sim.date + + DateOffset(months=p['main_polling_initialisation_delay_months'])) + + # If all services are included, set everything to True + if sim.modules['HealthSystem'].service_availability == ['*']: + for i in sim.modules['EmulatedRTI'].Rti_Services: + sim.modules['EmulatedRTI'].HS_conditions[i] = True + else: + for i in sim.modules['EmulatedRTI'].Rti_Services: + if (i + '_*') in sim.modules['HealthSystem'].service_availability: + sim.modules['EmulatedRTI'].HS_conditions[i] = True + else: + sim.modules['EmulatedRTI'].HS_conditions[i] = False + + # Begin logging the RTI events + sim.schedule_event(RTI_Logging_Event(self), sim.date + DateOffset(months=1)) + + def on_birth(self, mother_id, child_id): + """ + When a person is born this function sets up the default properties for the road traffic injuries module + :param mother_id: The mother + :param child_id: The newborn + :return: n/a + """ + df.at[child_id, 'rt_disability'] = 0 # default: no DALY + df.at[child_id, 'rt_disability_permanent'] = 0 # default: no DALY + df.at[child_id, 'rt_road_traffic_inc'] = False + df.at[child_id, 'rt_date_inj'] = pd.NaT + + def report_daly_values(self): + # This must send back a pd.Series or pd.DataFrame that reports on the average daly-weights that have been + # experienced by persons in the previous month. Only rows for alive-persons must be returned. + # The names of the series of columns is taken to be the label of the cause of this disability. + # It will be recorded by the healthburden module as _. + logger.debug(key='rti_general_message', data='This is RTI reporting my daly values') + df = self.sim.population.props + disability_series_for_alive_persons = df.loc[df.is_alive, "rt_disability"] + return disability_series_for_alive_persons + +class RTIPollingEvent(RegularEvent, PopulationScopeEventMixin): + """The regular RTI event which handles all the initial RTI related changes to the dataframe. It can be thought of + as the actual road traffic accident occurring. Specifically the event decides who is involved in a road traffic + accident every month (via the linear model helper class). The emulator then determines the outcome of that traffic accident. + """ + + def __init__(self, module): + """Schedule to take place + """ + super().__init__(module, frequency=DateOffset(months=module.parameters['main_polling_frequency'])) + p = module.parameters + # Parameters which transition the model between states + self.base_1m_prob_rti = (p['base_rate_injrti'] / 12) + if 'reduce_incidence' in p['allowed_interventions']: + self.base_1m_prob_rti = self.base_1m_prob_rti * p['intervention_incidence_reduction_factor'] + self.rr_injrti_age04 = p['rr_injrti_age04'] + self.rr_injrti_age59 = p['rr_injrti_age59'] + self.rr_injrti_age1017 = p['rr_injrti_age1017'] + self.rr_injrti_age1829 = p['rr_injrti_age1829'] + self.rr_injrti_age3039 = p['rr_injrti_age3039'] + self.rr_injrti_age4049 = p['rr_injrti_age4049'] + self.rr_injrti_age5059 = p['rr_injrti_age5059'] + self.rr_injrti_age6069 = p['rr_injrti_age6069'] + self.rr_injrti_age7079 = p['rr_injrti_age7079'] + self.rr_injrti_male = p['rr_injrti_male'] + self.rr_injrti_excessalcohol = p['rr_injrti_excessalcohol'] + self.imm_death_proportion_rti = p['imm_death_proportion_rti'] + self.prob_bleeding_leads_to_shock = p['prob_bleeding_leads_to_shock'] + self.rt_emergency_care_ISS_score_cut_off = p['rt_emergency_care_ISS_score_cut_off'] + + def apply(self, population): + """Apply this event to the population. + + :param population: the current population + """ + df = population.props + now = self.sim.date + + rt_current_non_ind = df.index[df.is_alive & ~df.rt_road_traffic_inc] + + # ========= Update for people currently not involved in a RTI, make some involved in a RTI event ============== + # Use linear model helper class + eq = LinearModel(LinearModelType.MULTIPLICATIVE, + self.base_1m_prob_rti, + Predictor('sex').when('M', self.rr_injrti_male), + Predictor( + 'age_years', + conditions_are_mutually_exclusive=True + ) + .when('.between(0,4)', self.rr_injrti_age04) + .when('.between(5,9)', self.rr_injrti_age59) + .when('.between(10,17)', self.rr_injrti_age1017) + .when('.between(18,29)', self.rr_injrti_age1829) + .when('.between(30,39)', self.rr_injrti_age3039) + .when('.between(40,49)', self.rr_injrti_age4049) + .when('.between(50,59)', self.rr_injrti_age5059) + .when('.between(60,69)', self.rr_injrti_age6069) + .when('.between(70,79)', self.rr_injrti_age7079), + Predictor('li_ex_alc').when(True, self.rr_injrti_excessalcohol) + ) + pred = eq.predict(df.loc[rt_current_non_ind]) + random_draw_in_rti = self.module.rng.random_sample(size=len(rt_current_non_ind)) + selected_for_rti = rt_current_non_ind[pred > random_draw_in_rti] + + # Update to say they have been involved in a rti + df.loc[selected_for_rti, 'rt_road_traffic_inc'] = True + # Set the date that people were injured to now + df.loc[selected_for_rti, 'rt_date_inj'] = now + + if len(selected_for_rti)>0: + # This is where we want to replace normal course of events for RTI with emulator. + + # First, sample outcomes for individuals which were selected_for_rti. + # For now, don't consider properties of individual when sampling outcome. All we care about is the number of samples. + condition_for_Rti = Condition( + num_rows=len(selected_for_rti), + column_values=self.sim.modules['EmulatedRTI'].HS_conditions + ) + NN_model = self.sim.modules['EmulatedRTI'].RTI_emulator.sample_from_conditions( + conditions=[condition_for_Rti], + ) + + # HS USAGE + # Get the total number of different types of appts that will be accessed as a result of this polling event and add to rolling count. + for column in self.sim.modules['EmulatedRTI'].HS_Use_Type: + self.sim.modules['EmulatedRTI'].HS_Use_by_RTI[column] += NN_model[column].sum() # Sum all values in the column + + # Change current properties of the individual and schedule resolution of event. + count = 0 + for person_id in selected_for_rti: + + # These properties are determined by the NN sampling + is_alive_after_RTI = NN_model.loc[count,'is_alive_after_RTI'] + + # Why does this require an int wrapper to work with DateOffset? + duration_days = int(NN_model.loc[count,'duration_days']) + + # Individual experiences an immediate death + if is_alive_after_RTI is False and duration_days == 0: + + # For each person selected to experience pre-hospital mortality, schedule an InstantaneosDeath event + self.sim.modules['Demography'].do_death(individual_id=individual_id, cause="RTI_imm_death", + originating_module=self.module) + + # Else individual doesn't immediately die, therefore schedule resolution + else: + # Set disability to what will be the average over duration of the episodef + df.loc[person_id,'rt_disability'] = NN_model.loc[count,'rt_disability_average'] + + # Make sure this person is not 'discoverable' by polling event next month. + #df.loc[person_id,'rt_inj_severity'] = 'mild' # instead of "none", but actually we don't know how severe it is + + # Schedule resolution + if is_alive_after_RTI: + # Store permanent disability incurred now to be accessed when Recovery Event is invoked. + df.loc[person_id,'rt_disability_permanent'] = NN_model.loc[count,'rt_disability_permanent'] + self.sim.schedule_event(RTI_NNResolution_Recovery_Event(self.module, person_id), df.loc[person_id, 'rt_date_inj'] + DateOffset(days=duration_days)) + else: + self.sim.schedule_event(RTI_NNResolution_Death_Event(self.module, person_id), df.loc[person_id, 'rt_date_inj'] + DateOffset(days=duration_days)) + + count += 1 + +class RTI_NNResolution_Death_Event(Event, IndividualScopeEventMixin): + """This is an individual-level event that determines the end of the incindent for individual via death""" + def __init__(self, module, person_id): + super().__init__(module, person_id=person_id) + def apply(self, person_id): + # For now invoking RTI_death_with_med, although technically we don't know whether individual accessed HSI or not. + # How finely do we want to really resolve RTI deaths? + self.sim.modules['Demography'].do_death(individual_id=person_id, cause="RTI_death_with_med", + originating_module=self.module) + +class RTI_NNResolution_Recovery_Event(Event, IndividualScopeEventMixin): + """This is an individual-level event that determines the end of the incindent for individual via death""" + def __init__(self, module, person_id): + super().__init__(module, person_id=person_id) + def apply(self, person_id): + df = self.sim.population.props + # Updat rt disability with long term outcome of accident + df.loc[person_id,'rt_disability'] = df.loc[person_id,'rt_disability_permanent'] + # Ensure that this person will be 'eligible' for rti injury again + df.loc[person_id,'rt_road_traffic_inc'] = False + +# --------------------------------------------------------------------------------------------------------- +# LOGGING EVENTS +# +# Put the logging events here. There should be a regular logger outputting current states of the +# population. There may also be a loggig event that is driven by particular events. +# --------------------------------------------------------------------------------------------------------- + + +class RTI_Logging_Event(RegularEvent, PopulationScopeEventMixin): + + def __init__(self, module): + """Produce a summary of the numbers of people with respect to the action of this module. + This is a regular event that can output current states of people or cumulative events since last logging event. + """ + + # run this event every month + self.repeat = 1 + super().__init__(module, frequency=DateOffset(months=self.repeat)) + assert isinstance(module, EmulatedRTI) + # Create variables used to store simulation data in + # Number of injured body region data + self.tot1inj = 0 + self.tot2inj = 0 + self.tot3inj = 0 + self.tot4inj = 0 + self.tot5inj = 0 + self.tot6inj = 0 + self.tot7inj = 0 + self.tot8inj = 0 + # Injury category data + self.totfracnumber = 0 + self.totdisnumber = 0 + self.tottbi = 0 + self.totsoft = 0 + self.totintorg = 0 + self.totintbled = 0 + self.totsci = 0 + self.totamp = 0 + self.toteye = 0 + self.totextlac = 0 + self.totburns = 0 + # Injury location on body data + self.totAIS1 = 0 + self.totAIS2 = 0 + self.totAIS3 = 0 + self.totAIS4 = 0 + self.totAIS5 = 0 + self.totAIS6 = 0 + self.totAIS7 = 0 + self.totAIS8 = 0 + # Injury severity data + self.totmild = 0 + self.totsevere = 0 + # More model progression data + self.totinjured = 0 + self.deathonscene = 0 + self.soughtmedcare = 0 + self.deathaftermed = 0 + self.deathwithoutmed = 0 + self.permdis = 0 + self.ISSscore = [] + self.severe_pain = 0 + self.moderate_pain = 0 + self.mild_pain = 0 + # Create variables for averages over time in the model + self.numerator = 0 + self.denominator = 0 + self.death_inc_numerator = 0 + self.death_in_denominator = 0 + self.fracdenominator = 0 + # Create variables to measure where certain injuries are located on the body + self.fracdist = [0, 0, 0, 0, 0, 0, 0, 0] + self.openwounddist = [0, 0, 0, 0, 0, 0, 0, 0] + self.burndist = [0, 0, 0, 0, 0, 0, 0, 0] + + def apply(self, population): + # Make some summary statistics + # Get the dataframe and isolate the important information + df = population.props + population_with_injuries = df.loc[df.rt_road_traffic_inc] + # ==================================== Incidence ============================================================== + # How many were involved in a RTI + n_in_RTI = int(df.rt_road_traffic_inc.sum()) + children_in_RTI = len(df.loc[df.rt_road_traffic_inc & (df['age_years'] < 19)]) + children_alive = len(df.loc[df['age_years'] < 19]) + self.numerator += n_in_RTI + self.totinjured += n_in_RTI + # How many were disabled + n_perm_disabled = int((df.is_alive & df.rt_disability_permanent>0).sum()) + # self.permdis += n_perm_disabled + n_alive = int(df.is_alive.sum()) + self.denominator += (n_alive - n_in_RTI) * (1 / 12) + diedfromrtiidx = df.index[df.cause_of_death == 'RTI'] + if (n_alive - n_in_RTI) > 0: + inc_rti = (n_in_RTI / ((n_alive - n_in_RTI) * (1 / 12))) * 100000 + else: + inc_rti = 0.0 + if (children_alive - children_in_RTI) > 0: + inc_rti_in_children = (children_in_RTI / ((children_alive - children_in_RTI) * (1 / 12))) * 100000 + else: + inc_rti_in_children = 0.0 + if (n_alive - len(diedfromrtiidx)) > 0: + inc_rti_death = (len(diedfromrtiidx) / ((n_alive - len(diedfromrtiidx)) * (1 / 12))) * 100000 + else: + inc_rti_death = 0.0 + + dict_to_output = { + 'number involved in a rti': n_in_RTI, + 'incidence of rti per 100,000': inc_rti, + 'incidence of rti per 100,000 in children': inc_rti_in_children, + 'incidence of rti death per 100,000': inc_rti_death, + 'number alive': n_alive, + 'number permanently disabled': n_perm_disabled, + } + logger.info(key='summary_1m', + data=dict_to_output, + description='Summary of the rti injuries in the last month') + # =========================== Get population demographics of those with RTIs ================================== + columnsOfInterest = ['sex', 'age_years', 'li_ex_alc'] + injuredDemographics = df.loc[df.rt_road_traffic_inc] + + injuredDemographics = injuredDemographics.loc[:, columnsOfInterest] + try: + percent_related_to_alcohol = len(injuredDemographics.loc[injuredDemographics.li_ex_alc]) / \ + len(injuredDemographics) + except ZeroDivisionError: + percent_related_to_alcohol = 0.0 + injured_demography_summary = { + 'males_in_rti': injuredDemographics['sex'].value_counts()['M'], + 'females_in_rti': injuredDemographics['sex'].value_counts()['F'], + 'age': injuredDemographics['age_years'].values.tolist(), + 'male_age': injuredDemographics.loc[injuredDemographics['sex'] == 'M', 'age_years'].values.tolist(), + 'female_age': injuredDemographics.loc[injuredDemographics['sex'] == 'F', 'age_years'].values.tolist(), + 'percent_related_to_alcohol': percent_related_to_alcohol, + } + logger.info(key='rti_demography', + data=injured_demography_summary, + description='Demographics of those in rti') + + # =================================== Flows through the model ================================================== + dict_to_output = {'total_injured': self.totinjured, + 'total_permanently_disabled': n_perm_disabled} + logger.info(key='model_progression', + data=dict_to_output, + description='Flows through the rti module') diff --git a/src/tlo/methods/fullmodel.py b/src/tlo/methods/fullmodel.py index 3f0c79434e..a30fac0baa 100644 --- a/src/tlo/methods/fullmodel.py +++ b/src/tlo/methods/fullmodel.py @@ -13,6 +13,7 @@ demography, depression, diarrhoea, + emulated_rti, enhanced_lifestyle, epi, epilepsy, @@ -41,6 +42,7 @@ def fullmodel( use_simplified_births: bool = False, + use_emulated_RTI = False, module_kwargs: Optional[Dict[str, Dict]] = {}, ) -> List[Module]: """Return a list of modules that should be registered in a run of the full model. @@ -48,6 +50,7 @@ def fullmodel( :param resourcefilepath: Path to root of directory containing resource files. :param use_simplified_births: Whether to use ``SimplifiedBirths`` module in place of full pregnancy related modules. + :param use_emulated_RTI: Option to use emulated version of RTI module :param module_kwargs: Dictionary mapping from module class names to dictionaries of keyword argument names and values to set for the module. If ``None`` (the default), the default values for all module keyword arguments are used other @@ -89,6 +92,12 @@ def fullmodel( postnatal_supervisor.PostnatalSupervisor, ] ), + *( + [emulated_rti.EmulatedRTI] if use_emulated_RTI else + [ + rti.RTI + ] + ), # Conditions of Early Childhood alri.Alri, diarrhoea.Diarrhoea, @@ -110,8 +119,6 @@ def fullmodel( prostate_cancer.ProstateCancer, # - Cardio-metabolic Disorders cardio_metabolic_disorders.CardioMetabolicDisorders, - # - Injuries - rti.RTI, # - Other Non-Communicable Conditions copd.Copd, depression.Depression, diff --git a/src/tlo/methods/healthsystem.py b/src/tlo/methods/healthsystem.py index b77ec86722..5ab8b6dc1d 100644 --- a/src/tlo/methods/healthsystem.py +++ b/src/tlo/methods/healthsystem.py @@ -2013,6 +2013,34 @@ def on_end_of_year(self) -> None: and self.parameters['scale_to_effective_capabilities'] ): self._rescale_capabilities_to_capture_effective_capability() + + # Add emulated appts to real ones in HS summary counters before logging the latter. + # Notes: + # 1) The emulated count is defined in the disease module itself, rather than in the HS, but + # is reset at the HS level after the emulated appts have been added to the total count. + # I *think* this is correct, but we may wish to discuss further. + # 2) Currently this is RTI-specific: in the future, this would have to be generalised to + # any emulated module + + # To total appt count + if 'EmulatedRTI' in self.sim.modules: + for key, value in self.sim.modules['EmulatedRTI'].HS_Use_by_RTI.items(): + # Extract the category part (ignoring the level) + _, category = key.split('_', 1) + self._summary_counter._appts[category] += value + + # To appt countbroken down by level + for key, value in self.sim.modules['EmulatedRTI'].HS_Use_by_RTI.items(): + # Split key into level and category + level, category = key.split('_', 1) + level_number = level.replace("Level", "") + + if level_number in self._summary_counter._appts_by_level: # Ensure the level exists in the second dictionary + self._summary_counter._appts_by_level[level_number][category] += value + + # Reset emulator counter. + self.sim.modules['EmulatedRTI'].HS_Use_by_RTI = Counter({col: 0 for col in self.sim.modules['EmulatedRTI'].HS_Use_Type}) + self._summary_counter.write_to_log_and_reset_counters() self.consumables.on_end_of_year() self.bed_days.on_end_of_year()