Calibrate Use of Anti-Epilepsy Medicines (#1511)

* log proportion of all epi on ASM

* introduce and use prob_start_anti_epilep_when_seizures_detected_in_generic_first_appt

* add plot of prop_ever_elipsey_on_antiepilepsy_meds

* make script faster and show progress bar

* update output to be infrequent/frequent seizures

* squash! update output to be infrequent/frequent seizures

* soft-code max_num_of_failed_attempts_before_defaulting and apply to initiation

* update plot to be red

* updated parameters

* add prob_start_anti_epilep_when_seizures_detected_in_generic_first_appt --> 1.0 for "max_healthsystem_function" scenario

* put script back up to fullsize

* update parameter file for final time

* update word document

* tidy up script and update all plots in write-up

src/scripts/epilepsy_analyses/analysis_epilepsy.py

@@ -28,7 +28,7 @@
 
 start_date = Date(2010, 1, 1)
 end_date = Date(2020,  1, 1)
-popsize = 200000
+popsize = 100_000
 
 # Establish the simulation object
 log_config = {
@@ -40,10 +40,11 @@
         'tlo.methods.demography': logging.INFO,
         'tlo.methods.healthsystem': logging.WARNING,
         'tlo.methods.healthburden': logging.WARNING,
+        'tlo.methods.population': logging.INFO,
     }
 }
 
-sim = Simulation(start_date=start_date, seed=0, log_config=log_config)
+sim = Simulation(start_date=start_date, seed=0, log_config=log_config, show_progress_bar=True)
 
 # make a dataframe that contains the switches for which interventions are allowed or not allowed
 # during this run. NB. These must use the exact 'registered strings' that the disease modules allow
@@ -125,7 +126,8 @@
 )
 n_seiz_stat_1_3.plot()
 plt.title('Number with epilepsy (past or current)')
-plt.ylim(0, 800000)
+plt.gca().set_ylim(bottom=0)
+plt.ylabel("Number (not scaled)")
 plt.tight_layout()
 plt.show()
 
@@ -135,11 +137,25 @@
 )
 n_seiz_stat_2_3.plot()
 plt.title('Number with epilepsy (infrequent or frequent seizures)')
-plt.ylim(0, 300000)
+plt.gca().set_ylim(bottom=0)
+plt.ylabel("Number (not scaled)")
 plt.tight_layout()
 plt.show()
 plt.clf()
 
+
+prop_antiepilep_seiz_infreq_or_freq = pd.Series(
+    output['tlo.methods.epilepsy']['epilepsy_logging']['prop_freq_or_infreq_seiz_on_antiep'].values,
+    index=output['tlo.methods.epilepsy']['epilepsy_logging']['date']
+)
+prop_antiepilep_seiz_infreq_or_freq.plot(color='r')
+plt.title('Proportion on antiepileptics\namongst people that have infrequent or frequent epileptic seizures')
+plt.ylim(0, 1)
+plt.tight_layout()
+plt.show()
+plt.clf()
+
+
 prop_antiepilep_seiz_stat_1 = pd.Series(
     output['tlo.methods.epilepsy']['epilepsy_logging']['prop_antiepilep_seiz_stat_1'].values,
     index=output['tlo.methods.epilepsy']['epilepsy_logging']['date']
@@ -179,7 +195,8 @@
 )
 n_epi_death.plot()
 plt.title('Number of deaths from epilepsy')
-plt.ylim(0, 50)
+plt.gca().set_ylim(bottom=0)
+plt.ylabel("Number (not scaled)")
 plt.tight_layout()
 plt.show()
 plt.clf()
@@ -190,11 +207,21 @@
 )
 n_antiep.plot()
 plt.title('Number of people on antiepileptics')
-plt.ylim(0, 50000)
+plt.gca().set_ylim(bottom=0)
+plt.ylabel("Number (not scaled)")
 plt.tight_layout()
 plt.show()
 plt.clf()
 
+(n_antiep / popsize).plot()
+plt.title('Proportion of of people (whole population) on antiepileptics')
+plt.gca().set_ylim(bottom=0)
+plt.ylabel("Number (not scaled)")
+plt.tight_layout()
+plt.show()
+plt.clf()
+
+
 epi_death_rate = pd.Series(
     output['tlo.methods.epilepsy']['epilepsy_logging']['epi_death_rate'].values,
     index=output['tlo.methods.epilepsy']['epilepsy_logging']['date']
@@ -233,8 +260,7 @@
 for _row, period in enumerate(('2010-2014', '2015-2019')):
     ax = axs[_row]
     comparison.loc[(period, slice(None), slice(None), CAUSE_NAME)]\
-              .droplevel([0, 1, 3])\
-              .groupby(axis=0, level=0)\
+              .groupby(axis=0, level=1)\
               .sum()\
               .plot(use_index=True, ax=ax)
     ax.set_ylabel('Deaths per year')

src/tlo/methods/epilepsy.py

@@ -100,6 +100,15 @@ def __init__(self, name=None, resourcefilepath=None):
         'daly_wt_epilepsy_seizure_free': Parameter(
             Types.REAL, 'disability weight for less severe epilepsy' 'controlled phase - code 862'
         ),
+        'prob_start_anti_epilep_when_seizures_detected_in_generic_first_appt': Parameter(
+            Types.REAL, 'probability that someone who has had a seizure is started on anti-epileptics. This is '
+                        'calibrated to induce the correct proportion of persons with epilepsy currently receiving '
+                        'anti-epileptics.'
+        ),
+        'max_num_of_failed_attempts_before_defaulting': Parameter(
+            Types.INT, 'maximum number of time an HSI can be repeated if the relevant essential consumables are not '
+                       'available.'
+        ),
     }
 
     """
@@ -406,8 +415,14 @@ def do_at_generic_first_appt_emergency(
         **kwargs,
     ) -> None:
         if "seizures" in symptoms:
-            event = HSI_Epilepsy_Start_Anti_Epileptic(person_id=person_id, module=self)
-            schedule_hsi_event(event, priority=0, topen=self.sim.date)
+            # Determine if treatment will start - depends on probability of prescribing, which is calibrated to
+            # induce the right proportion of persons with epilepsy receiving treatment.
+
+            prob_start = self.parameters['prob_start_anti_epilep_when_seizures_detected_in_generic_first_appt']
+
+            if self.rng.random_sample() < prob_start:
+                event = HSI_Epilepsy_Start_Anti_Epileptic(person_id=person_id, module=self)
+                schedule_hsi_event(event, priority=0, topen=self.sim.date)
 
 
 class EpilepsyEvent(RegularEvent, PopulationScopeEventMixin):
@@ -576,12 +591,17 @@ def apply(self, population):
 
         cum_deaths = (~df.is_alive).sum()
 
+        # Proportion of those with infrequent or frequent seizures currently on anti-epileptics
+        prop_freq_or_infreq_seiz_on_antiep = status_groups[2:].ep_antiep.sum() / status_groups[2:].is_alive.sum() \
+            if status_groups[2:].is_alive.sum() > 0 else 0
+
         logger.info(key='epilepsy_logging',
                     data={
                         'prop_seiz_stat_0': status_groups['prop_seiz_stats'].iloc[0],
                         'prop_seiz_stat_1': status_groups['prop_seiz_stats'].iloc[1],
                         'prop_seiz_stat_2': status_groups['prop_seiz_stats'].iloc[2],
                         'prop_seiz_stat_3': status_groups['prop_seiz_stats'].iloc[3],
+                        'prop_freq_or_infreq_seiz_on_antiep': prop_freq_or_infreq_seiz_on_antiep,
                         'prop_antiepilep_seiz_stat_0': status_groups['prop_seiz_stat_on_anti_ep'].iloc[0],
                         'prop_antiepilep_seiz_stat_1': status_groups['prop_seiz_stat_on_anti_ep'].iloc[1],
                         'prop_antiepilep_seiz_stat_2': status_groups['prop_seiz_stat_on_anti_ep'].iloc[2],
@@ -608,6 +628,9 @@ def __init__(self, module, person_id):
         self.EXPECTED_APPT_FOOTPRINT = self.make_appt_footprint({'Over5OPD': 1})
         self.ACCEPTED_FACILITY_LEVEL = '1b'
 
+        self._MAX_NUMBER_OF_FAILED_ATTEMPTS_BEFORE_DEFAULTING = module.parameters['max_num_of_failed_attempts_before_defaulting']
+        self._counter_of_failed_attempts_due_to_unavailable_medicines = 0
+
     def apply(self, person_id, squeeze_factor):
         df = self.sim.population.props
         hs = self.sim.modules["HealthSystem"]
@@ -639,8 +662,12 @@ def apply(self, person_id, squeeze_factor):
                 priority=0
             )
 
-        else:
+        elif (
+            self._counter_of_failed_attempts_due_to_unavailable_medicines
+            < self._MAX_NUMBER_OF_FAILED_ATTEMPTS_BEFORE_DEFAULTING
+        ):
             # If no medicine is available, run this HSI again next month
+            self._counter_of_failed_attempts_due_to_unavailable_medicines += 1
             self.module.sim.modules['HealthSystem'].schedule_hsi_event(hsi_event=self,
                                                                        topen=self.sim.date + pd.DateOffset(months=1),
                                                                        tclose=None,
@@ -652,7 +679,7 @@ class HSI_Epilepsy_Follow_Up(HSI_Event, IndividualScopeEventMixin):
     def __init__(self, module, person_id):
         super().__init__(module, person_id=person_id)
 
-        self._MAX_NUMBER_OF_FAILED_ATTEMPTS_BEFORE_DEFAULTING = 2
+        self._MAX_NUMBER_OF_FAILED_ATTEMPTS_BEFORE_DEFAULTING = module.parameters['max_num_of_failed_attempts_before_defaulting']
         self._DEFAULT_APPT_FOOTPRINT = self.make_appt_footprint({'Over5OPD': 1})
         self._REPEATED_APPT_FOOTPRINT = self.make_appt_footprint({'PharmDispensing': 1})