diff --git a/src/hivpy/column_names.py b/src/hivpy/column_names.py index a51af3b..ceec1d6 100644 --- a/src/hivpy/column_names.py +++ b/src/hivpy/column_names.py @@ -119,7 +119,12 @@ PREP_LEN_WILLING = "prep_len_willing" # Bool: True if an individual is willing to use injectable Lenacapavir PrEP PREP_VR_WILLING = "prep_vr_willing" # Bool: True if an individual is willing to use vaginal ring PrEP PREP_ANY_WILLING = "prep_any_willing" # Bool: True if an individual is willing to use at least one type of PrEP +FAVOURED_PREP_TYPE = "favoured_prep_type" # None | prep.PrEPType(enum): PrEP type with highest preference an individual is willing to take that is also currently available, o/w None PREP_ELIGIBLE = "prep_eligible" # Bool: True if an individual is eligible for PrEP usage +PREP_ORAL_TESTED = "prep_oral_tested" # Bool: True if an individual has tested explicitly to start oral PrEP (DUMMY) +PREP_CAB_TESTED = "prep_cab_tested" # Bool: True if an individual has tested explicitly to start injectable Cab PrEP (DUMMY) +PREP_LEN_TESTED = "prep_len_tested" # Bool: True if an individual has tested explicitly to start injectable Len PrEP (DUMMY) +PREP_VR_TESTED = "prep_vr_tested" # Bool: True if an individual has tested explicitly to start vaginal ring PrEP (DUMMY) PREP_TYPE = "prep_type" # None | prep.PrEPType(enum): Oral, Cabotegravir, Lenacapavir, or VaginalRing if PrEP is being used, o/w None EVER_PREP = "ever_prep" # Bool: True if an individual has ever been on PrEP FIRST_ORAL_START_DATE = "first_oral_start_date" # None | date: start date of first ever oral PrEP usage @@ -127,14 +132,16 @@ FIRST_LEN_START_DATE = "first_len_start_date" # None | date: start date of first ever injectable Len PrEP usage FIRST_VR_START_DATE = "first_vr_start_date" # None | date: start date of first ever vaginal ring PrEP usage LAST_PREP_START_DATE = "last_prep_start_date" # None | date: start date of most recent PrEP usage -PREP_JUST_STARTED = "prep_just_started" # Bool: True if PrEP usage began this time step (DUMMY) -PREP_ORAL_TESTED = "prep_oral_tested" # Bool: True if an individual has tested explicitly to start oral PrEP (DUMMY) -PREP_CAB_TESTED = "prep_cab_tested" # Bool: True if an individual has tested explicitly to start injectable Cab PrEP (DUMMY) -PREP_LEN_TESTED = "prep_len_tested" # Bool: True if an individual has tested explicitly to start injectable Len PrEP (DUMMY) -PREP_VR_TESTED = "prep_vr_tested" # Bool: True if an individual has tested explicitly to start vaginal ring PrEP (DUMMY) - -ART_ADHERENCE = "art_adherence" # DUMMY +LAST_PREP_STOP_DATE = "last_prep_stop_date" # None | date: stop date of most recent PrEP usage +PREP_JUST_STARTED = "prep_just_started" # Bool: True if PrEP usage began this time step +CONT_ON_PREP = "cont_on_prep" # None | timedelta: total length of continuous usage of current PrEP based on user intention (breaks due to ineligibility do not count against continuity) +CONT_ACTIVE_ON_PREP = "cont_active_on_prep" # None | timedelta: actual total length of continuous usage of current PrEP +CUMULATIVE_PREP_ORAL = "cumulative_prep_oral" # None | timedelta: total length of cumulative oral PrEP usage +CUMULATIVE_PREP_CAB = "cumulative_prep_cab" # None | timedelta: total length of cumulative injectable Cab PrEP usage +CUMULATIVE_PREP_LEN = "cumulative_prep_len" # None | timedelta: total length of cumulative injectable Len PrEP usage +CUMULATIVE_PREP_VR = "cumulative_prep_vr" # None | timedelta: total length of cumulative injectable vaginal ring PrEP usage +ART_ADHERENCE = "art_adherence" # None | Float: percentage of ART intake that is adhered to TA_MUTATION = "tam" # X_MUTATION: drug resistance cols (TODO: all 24 currently DUMMIED) M184_MUTATION = "m184m" K65_MUTATION = "k65m" diff --git a/src/hivpy/data/prep.yaml b/src/hivpy/data/prep.yaml index 18427e5..4a4c18b 100644 --- a/src/hivpy/data/prep.yaml +++ b/src/hivpy/data/prep.yaml @@ -23,7 +23,11 @@ prep_willing_threshold: 0.2 prob_test_prep_start: Value: [0.25, 0.50, 0.75] prob_base_prep_start: - Value: [0.05, 0.10, 0.20] -# FIXME: do we need this to be separate from prob_base_prep_start? + Value: [0.10, 0.30] prob_prep_restart: Value: [0.05, 0.10, 0.20] +prob_base_prep_stop: + Value: [0.05, 0.15, 0.30] +prob_base_prep_stop_nonuniform: + Value: [0.05, 0.15, 0.30] + Probability: [0.8, 0.1, 0.1] diff --git a/src/hivpy/population.py b/src/hivpy/population.py index f1a5f66..1fabd8e 100644 --- a/src/hivpy/population.py +++ b/src/hivpy/population.py @@ -271,6 +271,17 @@ def general_func(g): # Use Dummy column to in order to enable transform method and avoid any risks to data return df.groupby(param_list, dropna=dropna)["Dummy"].transform(general_func) + def col_apply(self, param_list, func, sub_pop=None): + """ + Applies a function to specific columns in the dataframe. + """ + if sub_pop is not None: + df = self.data.loc[sub_pop] + else: + df = self.data + # Lambda function used to extract param column contents + return df.apply(lambda x: func(*[x[p] for p in param_list]), axis=1) + def evolve(self, time_step: timedelta): """ Advance the population by one time step. @@ -287,7 +298,7 @@ def evolve(self, time_step: timedelta): if self.HIV_introduced: self.hiv_status.set_primary_infection(self) self.hiv_status.set_viral_load_groups(self) - self.prep.prep_willingness(self) + self.prep.prep_propensity(self) self.prep.prep_eligibility(self) if self.circumcision.vmmc_disrup_covid: @@ -307,7 +318,7 @@ def evolve(self, time_step: timedelta): if (n_deaths and self.apply_death): self.drop_from_population(HIV_deaths) self.hiv_diagnosis.update_HIV_diagnosis(self) - self.prep.prep_usage(self) + self.prep.prep_usage(self, time_step) # Some population cleanup self.pregnancy.reset_anc_at_birth(self) diff --git a/src/hivpy/prep.py b/src/hivpy/prep.py index 2ebf029..dac38f3 100644 --- a/src/hivpy/prep.py +++ b/src/hivpy/prep.py @@ -37,6 +37,7 @@ def __init__(self, **kwargs): date(self.p_data.date_prep_cab_intro), date(self.p_data.date_prep_len_intro), date(self.p_data.date_prep_vr_intro)] + self.cab_available = True self.prob_risk_informed_prep = self.p_data.prob_risk_informed_prep self.prob_greater_risk_informed_prep = self.p_data.prob_greater_risk_informed_prep self.prob_suspect_risk_prep = self.p_data.prob_suspect_risk_prep @@ -59,8 +60,12 @@ def __init__(self, **kwargs): self.prob_cab_prep_start = self.prob_base_prep_start self.prob_len_prep_start = self.prob_base_prep_start self.prob_vr_prep_start = self.prob_base_prep_start - # FIXME: values are the same as for prob_base_prep_start, can we just re-sample that? self.prob_prep_restart = self.p_data.prob_prep_restart.sample() + # FIXME: stop probabilities dependent on time step length + self.prob_oral_prep_stop = self.p_data.prob_base_prep_stop.sample() + self.prob_cab_prep_stop = self.p_data.prob_base_prep_stop.sample() + self.prob_len_prep_stop = self.prob_cab_prep_stop + self.prob_vr_prep_stop = self.p_data.prob_base_prep_stop_nonuniform.sample() def init_prep_variables(self, pop: Population): pop.init_variable(col.PREP_ORAL_PREF, 0) @@ -76,8 +81,13 @@ def init_prep_variables(self, pop: Population): pop.init_variable(col.PREP_LEN_WILLING, False) pop.init_variable(col.PREP_VR_WILLING, False) pop.init_variable(col.PREP_ANY_WILLING, False) + pop.init_variable(col.FAVOURED_PREP_TYPE, None) pop.init_variable(col.R_PREP, 1.0) pop.init_variable(col.PREP_ELIGIBLE, False) + pop.init_variable(col.PREP_ORAL_TESTED, False) + pop.init_variable(col.PREP_CAB_TESTED, False) + pop.init_variable(col.PREP_LEN_TESTED, False) + pop.init_variable(col.PREP_VR_TESTED, False) pop.init_variable(col.PREP_TYPE, None) pop.init_variable(col.EVER_PREP, False) pop.init_variable(col.FIRST_ORAL_START_DATE, None) @@ -85,15 +95,31 @@ def init_prep_variables(self, pop: Population): pop.init_variable(col.FIRST_LEN_START_DATE, None) pop.init_variable(col.FIRST_VR_START_DATE, None) pop.init_variable(col.LAST_PREP_START_DATE, None) + pop.init_variable(col.LAST_PREP_STOP_DATE, None) pop.init_variable(col.PREP_JUST_STARTED, False) - pop.init_variable(col.PREP_ORAL_TESTED, False) - pop.init_variable(col.PREP_CAB_TESTED, False) - pop.init_variable(col.PREP_LEN_TESTED, False) - pop.init_variable(col.PREP_VR_TESTED, False) + pop.init_variable(col.CONT_ON_PREP, timedelta(months=0)) + pop.init_variable(col.CONT_ACTIVE_ON_PREP, timedelta(months=0)) + pop.init_variable(col.CUMULATIVE_PREP_ORAL, timedelta(months=0)) + pop.init_variable(col.CUMULATIVE_PREP_CAB, timedelta(months=0)) + pop.init_variable(col.CUMULATIVE_PREP_LEN, timedelta(months=0)) + pop.init_variable(col.CUMULATIVE_PREP_VR, timedelta(months=0)) pop.init_variable(col.LTP_HIV_STATUS, False) pop.init_variable(col.LTP_HIV_DIAGNOSED, False) pop.init_variable(col.LTP_ON_ART, False) + # FIXME: should this function be in another module? + def get_vl_prevalence(self, pop: Population): + """ + Return the prevalence of people between 15 and 50 years old with a viral load of over 3.0. + Affects willingness to take PrEP. + """ + gen_pop = len(pop.get_sub_pop([(col.AGE, op.ge, 15), (col.AGE, op.lt, 50)])) + # find prevalence of people with a viral load of over 3.0 + return (len(pop.get_sub_pop([(col.VIRAL_LOAD, op.ge, 3.0), + (col.AGE, op.ge, 15), + (col.AGE, op.lt, 50)])) / gen_pop + if gen_pop > 0 else 0) + def reroll_r_prep(self, pop: Population): """ Reroll the r_prep value for each individual that was ineligible for PrEP last time step. @@ -165,9 +191,25 @@ def get_presumed_hiv_neg_pop(self, pop: Population): return pop.apply_bool_mask(mask, false_neg_pop) - def set_prep_preference(self, pop: Population, date_intro, pref_beta, pref_col, willing_col, sub_pop_mod=None): + def prep_preference(self, pop: Population): """ - Set preference values for a specific type of PrEP and determine willingness. + Determine PrEP preferences for all PrEP types. + """ + # oral prep pref + self.set_prep_preference(pop, self.date_prep_intro[PrEPType.Oral], + self.prep_oral_pref_beta, col.PREP_ORAL_PREF) + # injectable prep pref + self.set_prep_preference(pop, self.date_prep_intro[PrEPType.Cabotegravir], + self.prep_cab_pref_beta, col.PREP_CAB_PREF) + self.set_prep_preference(pop, self.date_prep_intro[PrEPType.Lenacapavir], + self.prep_len_pref_beta, col.PREP_LEN_PREF) + # vr prep pref (women only) + self.set_prep_preference(pop, self.date_prep_intro[PrEPType.VaginalRing], self.prep_vr_pref_beta, + col.PREP_VR_PREF, sub_pop_mod=pop.get_sub_pop([(col.SEX, op.eq, SexType.Female)])) + + def set_prep_preference(self, pop: Population, date_intro, pref_beta, pref_col, sub_pop_mod=None): + """ + Set preference values for a specific type of PrEP. """ if pop.date >= date_intro: # find those who turned 15 this time step @@ -183,55 +225,12 @@ def set_prep_preference(self, pop: Population, date_intro, pref_beta, pref_col, # random preference beta distribution pref = rng.beta(pref_beta, 5, size=len(sub_pop)) pop.set_present_variable(pref_col, pref, sub_pop) - # determine willingness by comparing to threshold - willingness = pref > self.prep_willing_threshold - pop.set_present_variable(willing_col, willingness, sub_pop) - pop.set_present_variable(col.PREP_ANY_WILLING, True, pop.apply_bool_mask(willingness, sub_pop)) def prep_willingness(self, pop: Population): """ - Determine which individuals are willing to take PrEP, as well as their PrEP preferences. + Determine PrEP willingness for all PrEP types. """ - # initial preference values - init_prefs = pop.data[[col.PREP_ORAL_PREF, col.PREP_CAB_PREF, col.PREP_LEN_PREF, col.PREP_VR_PREF]] - # oral prep pref + willingness - self.set_prep_preference(pop, self.date_prep_intro[PrEPType.Oral], self.prep_oral_pref_beta, - col.PREP_ORAL_PREF, col.PREP_ORAL_WILLING) - # injectable prep pref + willingness - # FIXME: should Cab be controlled by an availability flag instead of introduction date? - self.set_prep_preference(pop, self.date_prep_intro[PrEPType.Cabotegravir], self.prep_cab_pref_beta, - col.PREP_CAB_PREF, col.PREP_CAB_WILLING) - self.set_prep_preference(pop, self.date_prep_intro[PrEPType.Lenacapavir], self.prep_len_pref_beta, - col.PREP_LEN_PREF, col.PREP_LEN_WILLING) - # vr prep pref + willingness (women only) - self.set_prep_preference(pop, self.date_prep_intro[PrEPType.VaginalRing], self.prep_vr_pref_beta, - col.PREP_VR_PREF, col.PREP_VR_WILLING, - sub_pop_mod=pop.get_sub_pop([(col.SEX, op.eq, SexType.Female)])) - - # new preference values - new_prefs = pop.data[[col.PREP_ORAL_PREF, col.PREP_CAB_PREF, col.PREP_LEN_PREF, col.PREP_VR_PREF]] - # find people whose preference has changed this time step - changed_pref_pop = new_prefs.compare(init_prefs).index - - if len(changed_pref_pop) > 0: - # get ranking outcomes - # FIXME: not sure if transform group is the best way to do this, but it works for now - pref_ranks = pop.transform_group([col.PREP_ORAL_PREF, col.PREP_CAB_PREF, - col.PREP_LEN_PREF, col.PREP_VR_PREF], - self.calc_prep_pref_ranks, sub_pop=changed_pref_pop, use_size=False) - # set ranks for each prep type - pop.set_present_variable(col.PREP_ORAL_RANK, [i[0] for i in pref_ranks], changed_pref_pop) - pop.set_present_variable(col.PREP_CAB_RANK, [i[1] for i in pref_ranks], changed_pref_pop) - pop.set_present_variable(col.PREP_LEN_RANK, [i[2] for i in pref_ranks], changed_pref_pop) - pop.set_present_variable(col.PREP_VR_RANK, [i[3] for i in pref_ranks], changed_pref_pop) - - gen_pop = len(pop.get_sub_pop([(col.AGE, op.ge, 15), (col.AGE, op.lt, 50)])) - # find prevalence of people with a viral load of over 1000 - vl_prevalence = (len(pop.get_sub_pop([(col.VIRAL_LOAD, op.gt, 1000), - (col.AGE, op.ge, 15), - (col.AGE, op.lt, 50)])) / gen_pop - if gen_pop > 0 else 0) - + vl_prevalence = self.get_vl_prevalence(pop) # there's a chance nobody is willing to take PrEP if unsuppressed viral load prevalence is too low if self.vl_prevalence_affects_prep and vl_prevalence < self.vl_prevalence_prep_threshold: pop.set_present_variable(col.PREP_ORAL_WILLING, False) @@ -239,6 +238,35 @@ def prep_willingness(self, pop: Population): pop.set_present_variable(col.PREP_LEN_WILLING, False) pop.set_present_variable(col.PREP_VR_WILLING, False) pop.set_present_variable(col.PREP_ANY_WILLING, False) + # otherwise set willingness as normal + else: + self.set_prep_willingness(pop, col.PREP_ORAL_PREF, col.PREP_ORAL_WILLING) + self.set_prep_willingness(pop, col.PREP_CAB_PREF, col.PREP_CAB_WILLING) + self.set_prep_willingness(pop, col.PREP_LEN_PREF, col.PREP_LEN_WILLING) + self.set_prep_willingness(pop, col.PREP_VR_PREF, col.PREP_VR_WILLING) + + def set_prep_willingness(self, pop: Population, pref_col, willing_col): + """ + Set willingness values for a specific type of PrEP. + """ + # determine willingness by comparing to threshold + willingness = pop.get_variable(pref_col) > self.prep_willing_threshold + pop.set_present_variable(willing_col, willingness) + pop.set_present_variable(col.PREP_ANY_WILLING, True, pop.apply_bool_mask(willingness)) + + def prep_pref_ranks(self, pop: Population, sub_pop=None): + """ + Rank PrEP preferences. + """ + # get ranking outcomes + pref_ranks = pop.col_apply([col.PREP_ORAL_PREF, col.PREP_CAB_PREF, + col.PREP_LEN_PREF, col.PREP_VR_PREF], + self.calc_prep_pref_ranks, sub_pop=sub_pop) + # set ranks for each prep type + pop.set_present_variable(col.PREP_ORAL_RANK, [i[0] for i in pref_ranks], sub_pop) + pop.set_present_variable(col.PREP_CAB_RANK, [i[1] for i in pref_ranks], sub_pop) + pop.set_present_variable(col.PREP_LEN_RANK, [i[2] for i in pref_ranks], sub_pop) + pop.set_present_variable(col.PREP_VR_RANK, [i[3] for i in pref_ranks], sub_pop) def calc_prep_pref_ranks(self, oral_pref, cab_pref, len_pref, vr_pref): """ @@ -251,7 +279,65 @@ def calc_prep_pref_ranks(self, oral_pref, cab_pref, len_pref, vr_pref): # assign rank per prep type (position indicates prep type, value indicates rank) for i in range(len(prefs)): ranks[sorted_pref_indices[i]] = i+1 - return [ranks] + return ranks + + def favoured_prep(self, pop: Population, sub_pop=None): + """ + Determine favoured PrEP type using preference ranks. Favoured PrEP is the type of PrEP an individual + is willing to take with the highest preference value that is also currently available. + """ + # FIXME: can we pass the date to transform_group in a better way? + self.date = pop.date + # find prep type with highest preference an individual is willing to take that is also currently available + favoured_prep = pop.transform_group([col.PREP_ORAL_RANK, col.PREP_CAB_RANK, + col.PREP_LEN_RANK, col.PREP_VR_RANK, + col.PREP_ORAL_WILLING, col.PREP_CAB_WILLING, + col.PREP_LEN_WILLING, col.PREP_VR_WILLING], + self.calc_favoured_prep, sub_pop=sub_pop, use_size=False) + pop.set_present_variable(col.FAVOURED_PREP_TYPE, favoured_prep, sub_pop) + + def calc_favoured_prep(self, oral_rank, cab_rank, len_rank, vr_rank, + oral_willing, cab_willing, len_willing, vr_willing): + """ + Returns favoured PrEP type based on willingness, preference rank and availability. + """ + # group pref ranks and willingness + prefs = [oral_rank, cab_rank, len_rank, vr_rank] + willing = [oral_willing, cab_willing, len_willing, vr_willing] + # zip prep type and willingness together and sort by pref rank + sorted_zipped = sorted(enumerate(willing), key=lambda x: prefs[x[0]]) + sorted_dict = dict(sorted_zipped) + + favoured_prep = None + # find prep type someone is willing to take with the highest pref that is currently available + for prep_type in sorted_dict: + willing = sorted_dict[prep_type] + if self.date >= self.date_prep_intro[prep_type] and willing: + if PrEPType(prep_type) is not PrEPType.Cabotegravir or self.cab_available: + favoured_prep = prep_type + break + + return favoured_prep + + def prep_propensity(self, pop: Population): + """ + Determine PrEP preference values, willingness to take PrEP, PrEP preference ranks, and favoured PrEP type. + """ + # store initial preference values + init_prefs = pop.data[[col.PREP_ORAL_PREF, col.PREP_CAB_PREF, col.PREP_LEN_PREF, col.PREP_VR_PREF]] + # set preference values + self.prep_preference(pop) + # set willingness values + self.prep_willingness(pop) + # get new preference values + new_prefs = pop.data[[col.PREP_ORAL_PREF, col.PREP_CAB_PREF, col.PREP_LEN_PREF, col.PREP_VR_PREF]] + # find people whose preference has changed this time step + changed_pref_pop = new_prefs.compare(init_prefs).index + if len(changed_pref_pop) > 0: + # update preference ranks + self.prep_pref_ranks(pop, changed_pref_pop) + # update favoured prep + self.favoured_prep(pop, changed_pref_pop) def prep_eligibility(self, pop: Population): """ @@ -420,7 +506,8 @@ def prep_eligibility(self, pop: Population): if len(prep_eligible_pop) > 0: pop.set_present_variable(col.PREP_ELIGIBLE, True, prep_eligible_pop) - def tested_start_prep(self, pop: Population, prep_eligible_pop, prep_type, prep_tested_col, first_start_col): + def tested_start_prep(self, pop: Population, prep_eligible_pop, prep_type, + prep_tested_col, first_start_col, time_step): """ Update people starting PrEP for the first time after testing to start PrEP. """ @@ -433,10 +520,17 @@ def tested_start_prep(self, pop: Population, prep_eligible_pop, prep_type, prep_ if len(starting_prep_pop) > 0: pop.set_present_variable(col.PREP_TYPE, prep_type, starting_prep_pop) pop.set_present_variable(col.EVER_PREP, True, starting_prep_pop) + pop.set_present_variable(col.PREP_JUST_STARTED, True, starting_prep_pop) + # set start dates pop.set_present_variable(col.LAST_PREP_START_DATE, pop.date, starting_prep_pop) pop.set_present_variable(first_start_col, pop.date, starting_prep_pop) + # set continuous use + pop.set_present_variable(col.CONT_ON_PREP, time_step, starting_prep_pop) + pop.set_present_variable(col.CONT_ACTIVE_ON_PREP, time_step, starting_prep_pop) + # increment cumulative use + self.set_all_prep_cumulative(pop, starting_prep_pop, time_step) - def general_start_prep(self, pop: Population, prep_eligible_pop): + def general_start_prep(self, pop: Population, prep_eligible_pop, time_step): """ Update people starting PrEP for the first time without specifically testing to start PrEP. """ @@ -448,103 +542,236 @@ def general_start_prep(self, pop: Population, prep_eligible_pop): COND(col.PREP_VR_TESTED, op.eq, False)))) if len(starting_prep_pop) > 0: - # FIXME: can we pass the date to transform_group in a better way? - self.date = pop.date # starting prep outcomes - prep_types = pop.transform_group([col.PREP_ORAL_RANK, col.PREP_CAB_RANK, - col.PREP_LEN_RANK, col.PREP_VR_RANK, - col.PREP_ORAL_WILLING, col.PREP_CAB_WILLING, - col.PREP_LEN_WILLING, col.PREP_VR_WILLING], - self.calc_willing_start_prep, sub_pop=starting_prep_pop) - + prep_types = pop.transform_group([col.FAVOURED_PREP_TYPE], self.calc_starting_prep, + sub_pop=starting_prep_pop, dropna=True) pop.set_present_variable(col.PREP_TYPE, prep_types, starting_prep_pop) pop.set_present_variable(col.EVER_PREP, True, starting_prep_pop) - pop.set_present_variable(col.LAST_PREP_START_DATE, pop.date, starting_prep_pop) + pop.set_present_variable(col.PREP_JUST_STARTED, True, starting_prep_pop) + # set start dates + self.set_all_prep_start_dates(pop, starting_prep_pop) + # set continuous use + pop.set_present_variable(col.CONT_ON_PREP, time_step, starting_prep_pop) + pop.set_present_variable(col.CONT_ACTIVE_ON_PREP, time_step, starting_prep_pop) + # increment cumulative use + self.set_all_prep_cumulative(pop, starting_prep_pop, time_step) + + def set_all_prep_start_dates(self, pop: Population, starting_prep_pop): + """ + Set current PrEP start date and all first PrEP start date columns. + """ + pop.set_present_variable(col.LAST_PREP_START_DATE, pop.date, starting_prep_pop) + self.set_prep_first_start_date(pop, starting_prep_pop, PrEPType.Oral, col.FIRST_ORAL_START_DATE) + self.set_prep_first_start_date(pop, starting_prep_pop, PrEPType.Cabotegravir, col.FIRST_CAB_START_DATE) + self.set_prep_first_start_date(pop, starting_prep_pop, PrEPType.Lenacapavir, col.FIRST_LEN_START_DATE) + self.set_prep_first_start_date(pop, starting_prep_pop, PrEPType.VaginalRing, col.FIRST_VR_START_DATE) - def set_prep_start_date(pop: Population, starting_prep_pop, prep_type, start_date_col): - """ - Set a specific start date column for the population starting a corresponding PrEP type. - """ - pop.set_present_variable(start_date_col, pop.date, - pop.get_sub_pop_intersection( - starting_prep_pop, - pop.get_sub_pop(COND(col.PREP_TYPE, op.eq, prep_type)))) + def set_prep_first_start_date(self, pop: Population, starting_prep_pop, prep_type, first_start_date_col): + """ + Set a specific start date column for the population starting a corresponding + PrEP type for the first time. Only overwrites if first start date is None. + """ + pop.set_present_variable( + first_start_date_col, pop.date, + pop.get_sub_pop_intersection( + starting_prep_pop, pop.get_sub_pop(AND(COND(col.PREP_TYPE, op.eq, prep_type), + COND(first_start_date_col, op.eq, None))))) - set_prep_start_date(pop, starting_prep_pop, PrEPType.Oral, col.FIRST_ORAL_START_DATE) - set_prep_start_date(pop, starting_prep_pop, PrEPType.Cabotegravir, col.FIRST_CAB_START_DATE) - set_prep_start_date(pop, starting_prep_pop, PrEPType.Lenacapavir, col.FIRST_LEN_START_DATE) - set_prep_start_date(pop, starting_prep_pop, PrEPType.VaginalRing, col.FIRST_VR_START_DATE) + def set_all_prep_cumulative(self, pop: Population, using_prep_pop, time_step): + """ + Increment all cumulative PrEP usage columns for active PrEP users. + """ + self.set_prep_cumulative(pop, using_prep_pop, PrEPType.Oral, col.CUMULATIVE_PREP_ORAL, time_step) + self.set_prep_cumulative(pop, using_prep_pop, PrEPType.Cabotegravir, col.CUMULATIVE_PREP_CAB, time_step) + self.set_prep_cumulative(pop, using_prep_pop, PrEPType.Lenacapavir, col.CUMULATIVE_PREP_LEN, time_step) + self.set_prep_cumulative(pop, using_prep_pop, PrEPType.VaginalRing, col.CUMULATIVE_PREP_VR, time_step) - def calc_willing_start_prep(self, oral_pref, cab_pref, len_pref, vr_pref, - oral_willing, cab_willing, len_willing, vr_willing, size): + def set_prep_cumulative(self, pop: Population, using_prep_pop, prep_type, cumulative_col, time_step): + """ + Increment a specific cumulative PrEP usage column for the population using a + corresponding PrEP type. Applies to those starting, continuing, or switching PrEP. + """ + prep_cont = pop.get_variable(cumulative_col) + time_step + pop.set_present_variable( + cumulative_col, prep_cont, + pop.get_sub_pop_intersection( + using_prep_pop, pop.get_sub_pop(COND(col.PREP_TYPE, op.eq, prep_type)))) + + def calc_starting_prep(self, favoured_prep, size): """ Returns PrEP types for people starting PrEP for the first time without explicitly testing to start PrEP. Individual preferences and availability are taken into account. """ - # group pref ranks and willingness - prefs = [oral_pref, cab_pref, len_pref, vr_pref] - willing = [oral_willing, cab_willing, len_willing, vr_willing] - # zip prep type and willingness together and sort by pref rank - sorted_zipped = sorted(enumerate(willing), key=lambda x: prefs[x[0]]) - sorted_dict = dict(sorted_zipped) - - starting_prep = None - # find prep type someone is willing to take with the highest pref that is currently available - for prep_type in sorted_dict: - willing = sorted_dict[prep_type] - if self.date >= self.date_prep_intro[prep_type] and willing: - starting_prep = prep_type - break - # outcomes r = rng.uniform(size=size) - if PrEPType(starting_prep) is PrEPType.Oral: + if PrEPType(favoured_prep) is PrEPType.Oral: starting = r < self.prob_oral_prep_start - elif PrEPType(starting_prep) is PrEPType.Cabotegravir: + elif PrEPType(favoured_prep) is PrEPType.Cabotegravir: starting = r < self.prob_cab_prep_start - elif PrEPType(starting_prep) is PrEPType.Lenacapavir: + elif PrEPType(favoured_prep) is PrEPType.Lenacapavir: starting = r < self.prob_len_prep_start - elif PrEPType(starting_prep) is PrEPType.VaginalRing: + elif PrEPType(favoured_prep) is PrEPType.VaginalRing: starting = r < self.prob_vr_prep_start - prep = [starting_prep if s else None for s in starting] + prep = [favoured_prep if s else None for s in starting] return prep - def start_prep(self, pop: Population): + def start_prep(self, pop: Population, time_step): """ Update PrEP usage for people starting PrEP for the first time. """ + # clear just_started flag + pop.set_present_variable(col.PREP_JUST_STARTED, False, + pop.get_sub_pop([(col.LAST_PREP_START_DATE, op.ne, pop.date)])) + # find people eligible to start for the first time eligible = pop.get_sub_pop([(col.HARD_REACH, op.eq, False), (col.HIV_DIAGNOSED, op.eq, False), (col.PREP_ELIGIBLE, op.eq, True), (col.PREP_ANY_WILLING, op.eq, True), (col.EVER_PREP, op.eq, False), (col.LAST_TEST_DATE, op.eq, pop.date)]) - # factor in both true and false negatives in hiv status - starting_prep_pop = pop.get_sub_pop_intersection( - eligible, pop.get_sub_pop(OR(COND(col.HIV_STATUS, op.eq, False), - AND(COND(col.HIV_STATUS, op.eq, True), - COND(col.HIV_DIAGNOSED, op.eq, False))))) # starting oral prep after testing self.tested_start_prep( - pop, starting_prep_pop, PrEPType.Oral, col.PREP_ORAL_TESTED, col.FIRST_ORAL_START_DATE) + pop, eligible, PrEPType.Oral, col.PREP_ORAL_TESTED, col.FIRST_ORAL_START_DATE, time_step) # starting injectable cab prep after testing self.tested_start_prep( - pop, starting_prep_pop, PrEPType.Cabotegravir, col.PREP_CAB_TESTED, col.FIRST_CAB_START_DATE) + pop, eligible, PrEPType.Cabotegravir, col.PREP_CAB_TESTED, col.FIRST_CAB_START_DATE, time_step) # starting injectable len prep after testing self.tested_start_prep( - pop, starting_prep_pop, PrEPType.Lenacapavir, col.PREP_LEN_TESTED, col.FIRST_LEN_START_DATE) + pop, eligible, PrEPType.Lenacapavir, col.PREP_LEN_TESTED, col.FIRST_LEN_START_DATE, time_step) # starting vr prep after testing self.tested_start_prep( - pop, starting_prep_pop, PrEPType.VaginalRing, col.PREP_VR_TESTED, col.FIRST_VR_START_DATE) - + pop, eligible, PrEPType.VaginalRing, col.PREP_VR_TESTED, col.FIRST_VR_START_DATE, time_step) # not tested explicitly to start prep - self.general_start_prep(pop, starting_prep_pop) + self.general_start_prep(pop, eligible, time_step) + + def continue_prep(self, pop: Population, time_step): + """ + Update PrEP usage for people continuing PrEP. + """ + # people who have used prep before but not yet started this time step + eligible = pop.get_sub_pop(AND(COND(col.PREP_ELIGIBLE, op.eq, True), + COND(col.EVER_PREP, op.eq, True), + COND(col.PREP_JUST_STARTED, op.eq, False), + COND(col.LAST_PREP_STOP_DATE, op.eq, None), + OR(COND(col.LAST_TEST_DATE, op.ne, pop.date), + AND(COND(col.LAST_TEST_DATE, op.eq, pop.date), + COND(col.HIV_DIAGNOSED, op.eq, False))))) + + if len(eligible) > 0: + # continuous prep outcomes + prep_types = pop.transform_group([col.PREP_TYPE, col.FAVOURED_PREP_TYPE], + self.calc_current_prep, sub_pop=eligible, dropna=True) + # find various sub-populations + # people who are continuing current prep + continuing_prep_mask = pop.get_variable(col.PREP_TYPE, eligible) == prep_types + continuing_prep_pop = pop.apply_bool_mask(continuing_prep_mask, eligible) + # people who are switching prep + switching_prep_mask = (pop.get_variable(col.PREP_TYPE, eligible) != prep_types) & prep_types.notnull() + switching_prep_pop = pop.apply_bool_mask(switching_prep_mask, eligible) + # people who are either continuing or switching prep + using_prep_pop = pop.apply_bool_mask(prep_types.notnull(), eligible) + # people who are stopping prep + stopping_prep_pop = pop.apply_bool_mask(prep_types.isnull(), eligible) + + if len(continuing_prep_pop) > 0: + prep_cont = pop.get_variable(col.CONT_ON_PREP, eligible) + time_step + prep_active_cont = pop.get_variable(col.CONT_ACTIVE_ON_PREP, eligible) + time_step + # increment continuous use + pop.set_present_variable(col.CONT_ON_PREP, prep_cont, continuing_prep_pop) + pop.set_present_variable(col.CONT_ACTIVE_ON_PREP, prep_active_cont, continuing_prep_pop) + + if len(switching_prep_pop) > 0: + # set new prep types + pop.set_present_variable(col.PREP_TYPE, prep_types, switching_prep_pop) + # set start dates + self.set_all_prep_start_dates(pop, switching_prep_pop) + # reset continuous use + pop.set_present_variable(col.CONT_ON_PREP, time_step, switching_prep_pop) + pop.set_present_variable(col.CONT_ACTIVE_ON_PREP, time_step, switching_prep_pop) + + if len(using_prep_pop) > 0: + # increment cumulative use + self.set_all_prep_cumulative(pop, using_prep_pop, time_step) + + if len(stopping_prep_pop) > 0: + # stop continuous use + pop.set_present_variable(col.CONT_ON_PREP, timedelta(months=0), stopping_prep_pop) + pop.set_present_variable(col.CONT_ACTIVE_ON_PREP, timedelta(months=0), stopping_prep_pop) + # set stop date + pop.set_present_variable(col.LAST_PREP_STOP_DATE, pop.date, stopping_prep_pop) + + def calc_current_prep(self, prep_type, favoured_prep, size): + """ + Returns PrEP types for people continuing PrEP. + Individual preferences and availability are taken into account. + """ + # outcomes + r = rng.uniform(size=size) + if PrEPType(prep_type) is PrEPType.Oral: + continuing = r < (1 - self.prob_oral_prep_stop) + elif PrEPType(prep_type) is PrEPType.Cabotegravir: + continuing = r < (1 - self.prob_cab_prep_stop) + elif PrEPType(prep_type) is PrEPType.Lenacapavir: + continuing = r < (1 - self.prob_len_prep_stop) + elif PrEPType(prep_type) is PrEPType.VaginalRing: + continuing = r < (1 - self.prob_vr_prep_stop) + prep = [favoured_prep if c else None for c in continuing] + + return prep + + def restart_prep(self, pop: Population, time_step): + """ + Update PrEP usage for people restarting PrEP. + """ + # people who have used prep before and previously stopped using it + eligible = pop.get_sub_pop(AND(COND(col.HIV_DIAGNOSED, op.eq, False), + COND(col.PREP_ELIGIBLE, op.eq, True), + COND(col.EVER_PREP, op.eq, True), + COND(col.LAST_PREP_STOP_DATE, op.lt, pop.date), + COND(col.LAST_TEST_DATE, op.eq, pop.date))) + + if len(eligible) > 0: + # starting prep outcomes + prep_types = pop.transform_group([col.FAVOURED_PREP_TYPE], self.calc_restarting_prep, + sub_pop=eligible, dropna=True) + # people who are restarting prep + restarting_prep_pop = pop.apply_bool_mask(prep_types.notnull(), eligible) + + if len(restarting_prep_pop) > 0: + # set prep types + pop.set_present_variable(col.PREP_TYPE, prep_types, restarting_prep_pop) + pop.set_present_variable(col.PREP_JUST_STARTED, True, restarting_prep_pop) + # set start dates + self.set_all_prep_start_dates(pop, restarting_prep_pop) + # set continuous use + pop.set_present_variable(col.CONT_ON_PREP, time_step, restarting_prep_pop) + pop.set_present_variable(col.CONT_ACTIVE_ON_PREP, time_step, restarting_prep_pop) + # increment cumulative use + self.set_all_prep_cumulative(pop, restarting_prep_pop, time_step) + # unset stop date + pop.set_present_variable(col.LAST_PREP_STOP_DATE, None, restarting_prep_pop) + + def calc_restarting_prep(self, favoured_prep, size): + """ + Returns PrEP types for people restarting PrEP. + Individual preferences and availability are taken into account. + """ + # outcomes + r = rng.uniform(size=size) + restarting = r < self.prob_prep_restart + prep = [favoured_prep if r else None for r in restarting] + + return prep - def prep_usage(self, pop: Population): + def prep_usage(self, pop: Population, time_step): """ Update PrEP usage for people starting, continuing, switching, restarting, and stopping PrEP. """ # starting prep for the first time - self.start_prep(pop) + self.start_prep(pop, time_step) + # continuing prep + self.continue_prep(pop, time_step) + # restarting prep + self.restart_prep(pop, time_step) diff --git a/src/hivpy/prep_data.py b/src/hivpy/prep_data.py index b3a0df3..4699e5c 100644 --- a/src/hivpy/prep_data.py +++ b/src/hivpy/prep_data.py @@ -29,6 +29,8 @@ def __init__(self, filename): self.prob_test_prep_start = self._get_discrete_dist("prob_test_prep_start") self.prob_base_prep_start = self._get_discrete_dist("prob_base_prep_start") self.prob_prep_restart = self._get_discrete_dist("prob_prep_restart") + self.prob_base_prep_stop = self._get_discrete_dist("prob_base_prep_stop") + self.prob_base_prep_stop_nonuniform = self._get_discrete_dist("prob_base_prep_stop_nonuniform") except KeyError as ke: print(ke.args) diff --git a/src/tests/test_prep.py b/src/tests/test_prep.py index 9a02ae6..b8a2d7f 100644 --- a/src/tests/test_prep.py +++ b/src/tests/test_prep.py @@ -9,27 +9,110 @@ from hivpy.prep import PrEPType -def reset_prep_willingness_cols(pop: Population): +def reset_prep_propensity_cols(pop: Population): pop.set_present_variable(col.PREP_ORAL_PREF, 0) pop.set_present_variable(col.PREP_CAB_PREF, 0) pop.set_present_variable(col.PREP_LEN_PREF, 0) pop.set_present_variable(col.PREP_VR_PREF, 0) - pop.set_present_variable(col.PREP_ORAL_RANK, 0) - pop.set_present_variable(col.PREP_CAB_RANK, 0) - pop.set_present_variable(col.PREP_LEN_RANK, 0) - pop.set_present_variable(col.PREP_VR_RANK, 0) pop.set_present_variable(col.PREP_ORAL_WILLING, False) pop.set_present_variable(col.PREP_CAB_WILLING, False) pop.set_present_variable(col.PREP_LEN_WILLING, False) pop.set_present_variable(col.PREP_VR_WILLING, False) pop.set_present_variable(col.PREP_ANY_WILLING, False) + pop.set_present_variable(col.PREP_ORAL_RANK, 0) + pop.set_present_variable(col.PREP_CAB_RANK, 0) + pop.set_present_variable(col.PREP_LEN_RANK, 0) + pop.set_present_variable(col.PREP_VR_RANK, 0) + pop.set_present_variable(col.FAVOURED_PREP_TYPE, None) + + +def test_at_risk_pop(): + N = 100 + pop = Population(size=N, start_date=date(2020, 1, 1)) + # at_risk = num_stp >= 1 OR (ltp_diag AND not ltp_on_art) + pop.data[col.NUM_PARTNERS] = [0, 1] * (N // 2) + pop.data[col.LTP_HIV_DIAGNOSED] = [True, False, False, False] * (N // 4) + pop.data[col.LTP_ON_ART] = False + # 3/4 of people fulfill one of the conditions for being at risk + assert len(pop.prep.get_at_risk_pop(pop)) == N//4 * 3 + + +def test_risk_informed_pop(): + N = 1000 + pop = Population(size=N, start_date=date(2020, 1, 1)) + # risk_informed = ltp AND not ltp_on_art AND r < prob_risk_informed_prep + pop.data[col.LONG_TERM_PARTNER] = True + pop.data[col.LTP_ON_ART] = False + pop.data[col.LTP_HIV_STATUS] = False + pop.prep.reroll_r_prep(pop) + pop.prep.prob_risk_informed_prep = 0.1 + + # get stats + no_risk_informed = len(pop.prep.get_risk_informed_pop(pop, pop.prep.prob_risk_informed_prep)) + mean = N * pop.prep.prob_risk_informed_prep + stdev = sqrt(mean * (1 - pop.prep.prob_risk_informed_prep)) + # expecting ~10% of the population to be risk informed + assert mean - 3 * stdev <= no_risk_informed <= mean + 3 * stdev + + pop.data[col.LONG_TERM_PARTNER] = False + no_risk_informed = len(pop.prep.get_risk_informed_pop(pop, pop.prep.prob_risk_informed_prep)) + assert no_risk_informed == 0 + + +def test_suspect_risk_pop(): + N = 1000 + pop = Population(size=N, start_date=date(2020, 1, 1)) + # suspect_risk = ltp AND not ltp_on_art AND ltp_infected AND r < prob_suspect_risk_prep + pop.data[col.LONG_TERM_PARTNER] = True + pop.data[col.LTP_ON_ART] = [True, False] * (N//2) + pop.data[col.LTP_HIV_STATUS] = True + pop.prep.reroll_r_prep(pop) + pop.prep.prob_suspect_risk_prep = 0.5 + + # get stats + no_suspect_risk = len(pop.prep.get_suspect_risk_pop(pop)) + mean = (N/2) * pop.prep.prob_suspect_risk_prep + stdev = sqrt(mean * (1 - pop.prep.prob_suspect_risk_prep)) + # expecting ~50% of the population with partners NOT on ART to suspect they are at risk + assert mean - 3 * stdev <= no_suspect_risk <= mean + 3 * stdev -def test_prep_willingness(): +def test_presumed_hiv_neg_pop(): + N = 1000 + pop = Population(size=N, start_date=date(2020, 1, 1)) + pop.data[col.EVER_TESTED] = True + pop.data[col.HIV_DIAGNOSED] = True + pop.data[col.HIV_STATUS] = True + pop.data[col.DATE_HIV_INFECTION] = date(2019, 12, 1) + pop.hiv_diagnosis.init_prep_inj_na = True + pop.hiv_diagnosis.test_sens_general = 0.8 + pop.hiv_diagnosis.test_sens_primary_ab = 0.5 + + # nobody should be false negative because everyone is diagnosed + assert len(pop.prep.get_presumed_hiv_neg_pop(pop)) == 0 + + pop.data[col.HIV_DIAGNOSED] = False + # get stats (general test sensitivity) + no_presumed_hiv_neg = len(pop.prep.get_presumed_hiv_neg_pop(pop)) + mean = N * (1 - pop.hiv_diagnosis.test_sens_general) + stdev = sqrt(mean * pop.hiv_diagnosis.test_sens_general) + # expecting ~20% of the population to be false negative + assert mean - 3 * stdev <= no_presumed_hiv_neg <= mean + 3 * stdev + + pop.hiv_diagnosis.init_prep_inj_na = False + # get stats (primary test sensitivity) + no_presumed_hiv_neg = len(pop.prep.get_presumed_hiv_neg_pop(pop)) + mean = N * (1 - pop.hiv_diagnosis.test_sens_primary_ab) + stdev = sqrt(mean * pop.hiv_diagnosis.test_sens_primary_ab) + # expecting ~50% of the population to be false negative + assert mean - 3 * stdev <= no_presumed_hiv_neg <= mean + 3 * stdev + + +def test_prep_propensity(): N = 100 pop = Population(size=N, start_date=date(1999, 1, 1)) pop.data[col.AGE] = [10, 20] * (N // 2) - pop.data[col.VIRAL_LOAD] = 10000 + pop.data[col.VIRAL_LOAD] = 5.0 # all prep types have different intro dates pop.prep.date_prep_intro = [date(2000), date(3000), date(4000), date(5000)] # adjust chances of higher preference @@ -42,7 +125,7 @@ def test_prep_willingness(): pop.prep.vl_prevalence_prep_threshold = 0.5 # no willingness before prep intro date - pop.prep.prep_willingness(pop) + pop.prep.prep_propensity(pop) assert sum(pop.data[col.PREP_ANY_WILLING]) == 0 # find sub-pops @@ -50,7 +133,7 @@ def test_prep_willingness(): over_15s = pop.get_sub_pop([(col.AGE, op.ge, 15)]) # willingness calculated for under and over 15s pop.date = date(2000, 1, 1) - pop.prep.prep_willingness(pop) + pop.prep.prep_propensity(pop) # no willingness established for under 15s assert sum(pop.get_variable(col.PREP_ORAL_WILLING, under_15s)) == 0 @@ -64,9 +147,9 @@ def test_prep_willingness(): assert all(pop.get_variable(col.PREP_ORAL_RANK, under_15s) == 0) assert all(pop.get_variable(col.PREP_ORAL_RANK, over_15s) == 1) - reset_prep_willingness_cols(pop) + reset_prep_propensity_cols(pop) pop.date = date(3000, 1, 1) - pop.prep.prep_willingness(pop) + pop.prep.prep_propensity(pop) # no willingness established for under 15s assert sum(pop.get_variable(col.PREP_CAB_WILLING, under_15s)) == 0 @@ -80,9 +163,9 @@ def test_prep_willingness(): assert all(pop.get_variable(col.PREP_CAB_RANK, under_15s) == 0) assert all(pop.get_variable(col.PREP_CAB_RANK, over_15s) == 1) - reset_prep_willingness_cols(pop) + reset_prep_propensity_cols(pop) pop.date = date(4000, 1, 1) - pop.prep.prep_willingness(pop) + pop.prep.prep_propensity(pop) # no willingness established for under 15s assert sum(pop.get_variable(col.PREP_LEN_WILLING, under_15s)) == 0 @@ -96,9 +179,9 @@ def test_prep_willingness(): assert all(pop.get_variable(col.PREP_LEN_RANK, under_15s) == 0) assert all(pop.get_variable(col.PREP_LEN_RANK, over_15s) == 1) - reset_prep_willingness_cols(pop) + reset_prep_propensity_cols(pop) pop.date = date(5000, 1, 1) - pop.prep.prep_willingness(pop) + pop.prep.prep_propensity(pop) # no willingness established for under 15s assert sum(pop.get_variable(col.PREP_VR_WILLING, under_15s)) == 0 @@ -118,37 +201,37 @@ def test_prep_willingness(): over_15s, pop.get_sub_pop([(col.SEX, op.eq, SexType.Female)]))) == 1) # willingness calculated for those who turned 15 this time step - reset_prep_willingness_cols(pop) + reset_prep_propensity_cols(pop) pop.date = date(2020, 1, 1) pop.data[col.AGE] = 15 - pop.prep.prep_willingness(pop) + pop.prep.prep_propensity(pop) # some oral willingness established assert sum(pop.data[col.PREP_ORAL_WILLING]) > 0 assert sum(pop.data[col.PREP_CAB_WILLING]) == 0 assert sum(pop.data[col.PREP_LEN_WILLING]) == 0 assert sum(pop.data[col.PREP_VR_WILLING]) == 0 - reset_prep_willingness_cols(pop) + reset_prep_propensity_cols(pop) pop.date = date(3020, 1, 1) - pop.prep.prep_willingness(pop) + pop.prep.prep_propensity(pop) # some oral + cab willingness established assert sum(pop.data[col.PREP_ORAL_WILLING]) > 0 assert sum(pop.data[col.PREP_CAB_WILLING]) > 0 assert sum(pop.data[col.PREP_LEN_WILLING]) == 0 assert sum(pop.data[col.PREP_VR_WILLING]) == 0 - reset_prep_willingness_cols(pop) + reset_prep_propensity_cols(pop) pop.date = date(4020, 1, 1) - pop.prep.prep_willingness(pop) + pop.prep.prep_propensity(pop) # some oral + cab + len willingness established assert sum(pop.data[col.PREP_ORAL_WILLING]) > 0 assert sum(pop.data[col.PREP_CAB_WILLING]) > 0 assert sum(pop.data[col.PREP_LEN_WILLING]) > 0 assert sum(pop.data[col.PREP_VR_WILLING]) == 0 - reset_prep_willingness_cols(pop) + reset_prep_propensity_cols(pop) pop.date = date(5020, 1, 1) - pop.prep.prep_willingness(pop) + pop.prep.prep_propensity(pop) # some oral + cab + len + vr willingness established assert sum(pop.data[col.PREP_ORAL_WILLING]) > 0 assert sum(pop.data[col.PREP_CAB_WILLING]) > 0 @@ -156,94 +239,12 @@ def test_prep_willingness(): assert sum(pop.data[col.PREP_VR_WILLING]) > 0 # reset willingness with low viral load prevalence - pop.data[col.VIRAL_LOAD] = 100 - pop.prep.prep_willingness(pop) + pop.data[col.VIRAL_LOAD] = 2.0 + pop.prep.prep_propensity(pop) # no willingness remains assert sum(pop.data[col.PREP_ANY_WILLING]) == 0 -def test_at_risk_pop(): - N = 100 - pop = Population(size=N, start_date=date(2020, 1, 1)) - # at_risk = num_stp >= 1 OR (ltp_diag AND not ltp_on_art) - pop.data[col.NUM_PARTNERS] = [0, 1] * (N // 2) - pop.data[col.LTP_HIV_DIAGNOSED] = [True, False, False, False] * (N // 4) - pop.data[col.LTP_ON_ART] = False - # 3/4 of people fulfill one of the conditions for being at risk - assert len(pop.prep.get_at_risk_pop(pop)) == N//4 * 3 - - -def test_risk_informed_pop(): - N = 1000 - pop = Population(size=N, start_date=date(2020, 1, 1)) - # risk_informed = ltp AND not ltp_on_art AND r < prob_risk_informed_prep - pop.data[col.LONG_TERM_PARTNER] = True - pop.data[col.LTP_ON_ART] = False - pop.data[col.LTP_HIV_STATUS] = False - pop.prep.reroll_r_prep(pop) - pop.prep.prob_risk_informed_prep = 0.1 - - # get stats - no_risk_informed = len(pop.prep.get_risk_informed_pop(pop, pop.prep.prob_risk_informed_prep)) - mean = N * pop.prep.prob_risk_informed_prep - stdev = sqrt(mean * (1 - pop.prep.prob_risk_informed_prep)) - # expecting ~10% of the population to be risk informed - assert mean - 3 * stdev <= no_risk_informed <= mean + 3 * stdev - - pop.data[col.LONG_TERM_PARTNER] = False - no_risk_informed = len(pop.prep.get_risk_informed_pop(pop, pop.prep.prob_risk_informed_prep)) - assert no_risk_informed == 0 - - -def test_suspect_risk_pop(): - N = 1000 - pop = Population(size=N, start_date=date(2020, 1, 1)) - # suspect_risk = ltp AND not ltp_on_art AND ltp_infected AND r < prob_suspect_risk_prep - pop.data[col.LONG_TERM_PARTNER] = True - pop.data[col.LTP_ON_ART] = [True, False] * (N//2) - pop.data[col.LTP_HIV_STATUS] = True - pop.prep.reroll_r_prep(pop) - pop.prep.prob_suspect_risk_prep = 0.5 - - # get stats - no_suspect_risk = len(pop.prep.get_suspect_risk_pop(pop)) - mean = (N/2) * pop.prep.prob_suspect_risk_prep - stdev = sqrt(mean * (1 - pop.prep.prob_suspect_risk_prep)) - # expecting ~50% of the population with partners NOT on ART to suspect they are at risk - assert mean - 3 * stdev <= no_suspect_risk <= mean + 3 * stdev - - -def test_presumed_hiv_neg_pop(): - N = 1000 - pop = Population(size=N, start_date=date(2020, 1, 1)) - pop.data[col.EVER_TESTED] = True - pop.data[col.HIV_DIAGNOSED] = True - pop.data[col.HIV_STATUS] = True - pop.data[col.DATE_HIV_INFECTION] = date(2019, 12, 1) - pop.hiv_diagnosis.init_prep_inj_na = True - pop.hiv_diagnosis.test_sens_general = 0.8 - pop.hiv_diagnosis.test_sens_primary_ab = 0.5 - - # nobody should be false negative because everyone is diagnosed - assert len(pop.prep.get_presumed_hiv_neg_pop(pop)) == 0 - - pop.data[col.HIV_DIAGNOSED] = False - # get stats (general test sensitivity) - no_presumed_hiv_neg = len(pop.prep.get_presumed_hiv_neg_pop(pop)) - mean = N * (1 - pop.hiv_diagnosis.test_sens_general) - stdev = sqrt(mean * pop.hiv_diagnosis.test_sens_general) - # expecting ~20% of the population to be false negative - assert mean - 3 * stdev <= no_presumed_hiv_neg <= mean + 3 * stdev - - pop.hiv_diagnosis.init_prep_inj_na = False - # get stats (primary test sensitivity) - no_presumed_hiv_neg = len(pop.prep.get_presumed_hiv_neg_pop(pop)) - mean = N * (1 - pop.hiv_diagnosis.test_sens_primary_ab) - stdev = sqrt(mean * pop.hiv_diagnosis.test_sens_primary_ab) - # expecting ~50% of the population to be false negative - assert mean - 3 * stdev <= no_presumed_hiv_neg <= mean + 3 * stdev - - @pytest.mark.parametrize("prep_strategy", [i for i in range(1, 17)]) def test_prep_ineligible(prep_strategy): N = 1000 @@ -586,6 +587,7 @@ def test_prep_eligibility_all(): def test_starting_prep(): N = 1000 + time_step = timedelta(months=1) pop = Population(size=N, start_date=date(5000, 1, 1)) pop.prep.date_prep_intro = [date(2000), date(3000), date(4000), date(5000)] pop.data[col.HARD_REACH] = False @@ -595,18 +597,31 @@ def test_starting_prep(): pop.data[col.PREP_ANY_WILLING] = True pop.data[col.EVER_PREP] = False pop.data[col.LAST_TEST_DATE] = pop.date + pop.data[col.CONT_ON_PREP] = None + pop.data[col.CONT_ACTIVE_ON_PREP] = None + pop.data[col.CUMULATIVE_PREP_ORAL] = timedelta(months=0) + pop.data[col.CUMULATIVE_PREP_CAB] = timedelta(months=0) + pop.data[col.CUMULATIVE_PREP_LEN] = timedelta(months=0) + pop.data[col.CUMULATIVE_PREP_VR] = timedelta(months=0) # tested explicitly to start prep pop.data[col.PREP_ORAL_TESTED] = [True, False, False, False] * (N // 4) pop.data[col.PREP_CAB_TESTED] = [False, True, False, False] * (N // 4) pop.data[col.PREP_LEN_TESTED] = [False, False, True, False] * (N // 4) pop.data[col.PREP_VR_TESTED] = [False, False, False, True] * (N // 4) - pop.prep.start_prep(pop) + pop.prep.start_prep(pop, time_step) # prep types spread evenly among population assert sum(pop.data[col.PREP_TYPE] == PrEPType.Oral) == N/4 assert sum(pop.data[col.PREP_TYPE] == PrEPType.Cabotegravir) == N/4 assert sum(pop.data[col.PREP_TYPE] == PrEPType.Lenacapavir) == N/4 assert sum(pop.data[col.PREP_TYPE] == PrEPType.VaginalRing) == N/4 + # check continuous and cumulative prep usage + assert sum(pop.data[col.CONT_ON_PREP] == time_step) == N + assert sum(pop.data[col.CONT_ACTIVE_ON_PREP] == time_step) == N + assert sum(pop.data[col.CUMULATIVE_PREP_ORAL] == time_step) == N/4 + assert sum(pop.data[col.CUMULATIVE_PREP_CAB] == time_step) == N/4 + assert sum(pop.data[col.CUMULATIVE_PREP_LEN] == time_step) == N/4 + assert sum(pop.data[col.CUMULATIVE_PREP_VR] == time_step) == N/4 pop.data[col.PREP_TYPE] = None pop.data[col.EVER_PREP] = [True, False] * (N // 2) @@ -615,7 +630,7 @@ def test_starting_prep(): pop.data[col.FIRST_LEN_START_DATE] = None pop.data[col.FIRST_VR_START_DATE] = None pop.data[col.LAST_PREP_START_DATE] = None - pop.prep.start_prep(pop) + pop.prep.start_prep(pop, time_step) # only 50% eligible to start prep for the first time assert sum(pop.data[col.PREP_TYPE].isnull()) == N/2 @@ -649,7 +664,8 @@ def test_starting_prep(): pop.prep.prob_len_prep_start = 0.7 pop.prep.prob_vr_prep_start = 0.6 - pop.prep.start_prep(pop) + pop.prep.favoured_prep(pop, None) + pop.prep.start_prep(pop, time_step) # test oral prep type start probability no_on_oral = sum(pop.data[col.PREP_TYPE] == PrEPType.Oral) mean = N/4 * pop.prep.prob_oral_prep_start @@ -673,25 +689,39 @@ def test_starting_prep(): pop.data[col.PREP_TYPE] = None pop.data[col.EVER_PREP] = False + pop.data[col.CONT_ON_PREP] = None + pop.data[col.CONT_ACTIVE_ON_PREP] = None + pop.data[col.CUMULATIVE_PREP_ORAL] = timedelta(months=0) + pop.data[col.CUMULATIVE_PREP_CAB] = timedelta(months=0) + pop.data[col.CUMULATIVE_PREP_LEN] = timedelta(months=0) + pop.data[col.CUMULATIVE_PREP_VR] = timedelta(months=0) # 100% chance to start prep pop.prep.prob_oral_prep_start = 1 pop.prep.prob_cab_prep_start = 1 pop.prep.prob_len_prep_start = 1 pop.prep.prob_vr_prep_start = 1 - pop.prep.start_prep(pop) + pop.prep.start_prep(pop, time_step) # everyone starts their most preferred prep type assert all((pop.data[col.PREP_TYPE] == PrEPType.Oral) == (pop.data[col.PREP_ORAL_RANK] == 1)) assert all((pop.data[col.PREP_TYPE] == PrEPType.Cabotegravir) == (pop.data[col.PREP_CAB_RANK] == 1)) assert all((pop.data[col.PREP_TYPE] == PrEPType.Lenacapavir) == (pop.data[col.PREP_LEN_RANK] == 1)) assert all((pop.data[col.PREP_TYPE] == PrEPType.VaginalRing) == (pop.data[col.PREP_VR_RANK] == 1)) + # check continuous and cumulative prep usage + assert sum(pop.data[col.CONT_ON_PREP] == time_step) == N + assert sum(pop.data[col.CONT_ACTIVE_ON_PREP] == time_step) == N + assert all((pop.data[col.CUMULATIVE_PREP_ORAL] == time_step) == (pop.data[col.PREP_ORAL_RANK] == 1)) + assert all((pop.data[col.CUMULATIVE_PREP_CAB] == time_step) == (pop.data[col.PREP_CAB_RANK] == 1)) + assert all((pop.data[col.CUMULATIVE_PREP_LEN] == time_step) == (pop.data[col.PREP_LEN_RANK] == 1)) + assert all((pop.data[col.CUMULATIVE_PREP_VR] == time_step) == (pop.data[col.PREP_VR_RANK] == 1)) pop.data[col.PREP_TYPE] = None pop.data[col.EVER_PREP] = False # nobody is willing to take oral or cab pop.data[col.PREP_ORAL_WILLING] = False pop.data[col.PREP_CAB_WILLING] = False - pop.prep.start_prep(pop) + pop.prep.favoured_prep(pop, None) + pop.prep.start_prep(pop, time_step) # everyone is either on len or vr assert sum(pop.data[col.PREP_TYPE] == PrEPType.Lenacapavir) == N * 0.75 @@ -704,7 +734,137 @@ def test_starting_prep(): pop.data[col.PREP_TYPE] = None pop.data[col.EVER_PREP] = False pop.prep.date_prep_intro = [date(2000), date(3000), date(4000), date(6000)] - pop.prep.start_prep(pop) + pop.prep.favoured_prep(pop, None) + pop.prep.start_prep(pop, time_step) # everyone is on len because vr is not yet available assert all(pop.data[col.PREP_TYPE] == PrEPType.Lenacapavir) assert all(pop.data[col.FIRST_LEN_START_DATE] == pop.date) + + +def test_continuing_prep(): + N = 100 + time_step = timedelta(months=1) + pop = Population(size=N, start_date=date(5000, 1, 1)) + pop.prep.date_prep_intro = [date(2000), date(3000), date(4000), date(5000)] + pop.data[col.HIV_DIAGNOSED] = False + pop.data[col.PREP_ELIGIBLE] = True + pop.data[col.EVER_PREP] = True + pop.data[col.LAST_PREP_STOP_DATE] = None + pop.data[col.PREP_JUST_STARTED] = False + pop.data[col.CONT_ON_PREP] = timedelta(months=3) + pop.data[col.CONT_ACTIVE_ON_PREP] = timedelta(months=2) + pop.data[col.CUMULATIVE_PREP_ORAL] = time_step + pop.data[col.CUMULATIVE_PREP_CAB] = time_step + pop.data[col.CUMULATIVE_PREP_LEN] = time_step + pop.data[col.CUMULATIVE_PREP_VR] = time_step + pop.data[col.LAST_TEST_DATE] = pop.date - timedelta(months=3) + # prep types spread evenly among population + pop.data[col.PREP_TYPE] = [PrEPType.Oral, PrEPType.Cabotegravir, + PrEPType.Lenacapavir, PrEPType.VaginalRing] * (N // 4) + # everyone is taking their favoured prep type + pop.data[col.FAVOURED_PREP_TYPE] = [PrEPType.Oral, PrEPType.Cabotegravir, + PrEPType.Lenacapavir, PrEPType.VaginalRing] * (N // 4) + # 10% chance to stop prep + prob_base_prep_stop = 0.1 + pop.prep.prob_oral_prep_stop = prob_base_prep_stop + pop.prep.prob_cab_prep_stop = prob_base_prep_stop + pop.prep.prob_len_prep_stop = prob_base_prep_stop + pop.prep.prob_vr_prep_stop = prob_base_prep_stop + + pop.prep.continue_prep(pop, time_step) + # expecting 90% of people to continue prep + no_on_prep = sum(pop.data[col.CONT_ACTIVE_ON_PREP] == timedelta(months=3)) + mean = N * (1 - prob_base_prep_stop) + stdev = sqrt(mean * prob_base_prep_stop) + assert mean - 3 * stdev <= no_on_prep <= mean + 3 * stdev + # expecting 10% of people to stop prep + no_off_prep = sum(pop.data[col.LAST_PREP_STOP_DATE] == pop.date) + mean = N * prob_base_prep_stop + stdev = sqrt(mean * (1 - prob_base_prep_stop)) + assert mean - 3 * stdev <= no_off_prep <= mean + 3 * stdev + + # check cumulative prep usage (those that stopped should not be incremented) + assert all(((pop.data[col.CUMULATIVE_PREP_ORAL] == timedelta(months=2)) == + (pop.data[col.PREP_TYPE] == PrEPType.Oral)) | + ((pop.data[col.CUMULATIVE_PREP_ORAL] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE] == pop.date))) + assert all(((pop.data[col.CUMULATIVE_PREP_CAB] == timedelta(months=2)) == + (pop.data[col.PREP_TYPE] == PrEPType.Cabotegravir)) | + ((pop.data[col.CUMULATIVE_PREP_CAB] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE] == pop.date))) + assert all(((pop.data[col.CUMULATIVE_PREP_LEN] == timedelta(months=2)) == + (pop.data[col.PREP_TYPE] == PrEPType.Lenacapavir)) | + ((pop.data[col.CUMULATIVE_PREP_LEN] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE] == pop.date))) + assert all(((pop.data[col.CUMULATIVE_PREP_VR] == timedelta(months=2)) == + (pop.data[col.PREP_TYPE] == PrEPType.VaginalRing)) | + ((pop.data[col.CUMULATIVE_PREP_VR] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE] == pop.date))) + + pop.data[col.LAST_PREP_STOP_DATE] = None + # nobody is taking their favoured prep type + pop.data[col.FAVOURED_PREP_TYPE] = [PrEPType.VaginalRing, PrEPType.Lenacapavir, + PrEPType.Cabotegravir, PrEPType.Oral] * (N // 4) + + pop.prep.continue_prep(pop, time_step) + # expecting 90% of people to switch prep + no_on_prep = sum(pop.data[col.CONT_ACTIVE_ON_PREP] == timedelta(months=1)) + mean = N * (1 - prob_base_prep_stop) + stdev = sqrt(mean * prob_base_prep_stop) + assert mean - 3 * stdev <= no_on_prep <= mean + 3 * stdev + # expecting 10% of people to stop prep + no_off_prep = sum(pop.data[col.LAST_PREP_STOP_DATE] == pop.date) + mean = N * prob_base_prep_stop + stdev = sqrt(mean * (1 - prob_base_prep_stop)) + assert mean - 3 * stdev <= no_off_prep <= mean + 3 * stdev + + # check cumulative prep usage (those that stopped should not be incremented) + assert all(((pop.data[col.CUMULATIVE_PREP_ORAL] == timedelta(months=2)) == + (pop.data[col.PREP_TYPE] == PrEPType.Oral)) | + ((pop.data[col.CUMULATIVE_PREP_ORAL] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE] == pop.date))) + assert all(((pop.data[col.CUMULATIVE_PREP_CAB] == timedelta(months=2)) == + (pop.data[col.PREP_TYPE] == PrEPType.Cabotegravir)) | + ((pop.data[col.CUMULATIVE_PREP_CAB] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE] == pop.date))) + assert all(((pop.data[col.CUMULATIVE_PREP_LEN] == timedelta(months=2)) == + (pop.data[col.PREP_TYPE] == PrEPType.Lenacapavir)) | + ((pop.data[col.CUMULATIVE_PREP_LEN] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE] == pop.date))) + assert all(((pop.data[col.CUMULATIVE_PREP_VR] == timedelta(months=2)) == + (pop.data[col.PREP_TYPE] == PrEPType.VaginalRing)) | + ((pop.data[col.CUMULATIVE_PREP_VR] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE] == pop.date))) + + +def test_restarting_prep(): + N = 100 + time_step = timedelta(months=1) + pop = Population(size=N, start_date=date(5000, 1, 1)) + pop.prep.date_prep_intro = [date(2000), date(3000), date(4000), date(5000)] + pop.data[col.HIV_DIAGNOSED] = False + pop.data[col.PREP_ELIGIBLE] = True + pop.data[col.EVER_PREP] = True + pop.data[col.LAST_PREP_STOP_DATE] = pop.date - time_step + pop.data[col.PREP_JUST_STARTED] = False + pop.data[col.CONT_ON_PREP] = timedelta(months=0) + pop.data[col.CONT_ACTIVE_ON_PREP] = timedelta(months=0) + pop.data[col.LAST_TEST_DATE] = pop.date + pop.data[col.PREP_TYPE] = None + # everyone is taking their favoured prep type + pop.data[col.FAVOURED_PREP_TYPE] = [PrEPType.Oral, PrEPType.Cabotegravir, + PrEPType.Lenacapavir, PrEPType.VaginalRing] * (N // 4) + # 50% chance to restart prep + pop.prep.prob_prep_restart = 0.5 + + pop.prep.restart_prep(pop, time_step) + # expecting 90% of people to restart prep + no_on_prep = sum(pop.data[col.LAST_PREP_STOP_DATE].isnull()) + mean = N * pop.prep.prob_prep_restart + stdev = sqrt(mean * (1 - pop.prep.prob_prep_restart)) + assert mean - 3 * stdev <= no_on_prep <= mean + 3 * stdev + # check continuous prep usage + assert all((pop.data[col.CONT_ON_PREP] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE].isnull())) + assert all((pop.data[col.CONT_ACTIVE_ON_PREP] == time_step) == + (pop.data[col.LAST_PREP_STOP_DATE].isnull()))