Starting PrEP usage

src/hivpy/column_names.py

@@ -88,20 +88,34 @@
 WHO4_OTHER = "who4_other"                       # Bool: True if other WHO4 disease occurs this timestep
 WHO4_OTHER_DIAGNOSED = "who4_other_diagnosed"   # Bool: True if other WHO4 disease diagnosed this timestep
 ADC = "AIDS_defining_condition"                 # Bool: presence of AIDS defining condition (any WHO4)
+
 R_PREP = "r_prep"                               # float: a semi-permanent random personal variable that determines whether someone is risk informed or suspects they are at risk enough to take PrEP
 PREP_ORAL_PREF = "prep_oral_pref"               # float: a value that determines and individual's preference for oral PrEP
 PREP_CAB_PREF = "prep_cab_pref"                 # float: a value that determines and individual's preference for injectable Cabotegravir PrEP
 PREP_LEN_PREF = "prep_len_pref"                 # float: a value that determines and individual's preference for injectable Lenacapavir PrEP
 PREP_VR_PREF = "prep_vr_pref"                   # float: a value that determines and individual's preference for vaginal ring PrEP
-PREP_PREF_RANKED = "prep_pref_ranked"           # None | prep.PrEPType(enum) list: a ranked list of PrEP type preferences
+PREP_ORAL_RANK = "prep_oral_rank"               # int: an individual's preference ranking for oral PrEP
+PREP_CAB_RANK = "prep_cab_rank"                 # int: an individual's preference ranking for injectable Cabotegravir PrEP
+PREP_LEN_RANK = "prep_len_rank"                 # int: an individual's preference ranking for injectable Lenacapavir PrEP
+PREP_VR_RANK = "prep_vr_rank"                   # int: an individual's preference ranking for vaginal ring PrEP
 PREP_ORAL_WILLING = "prep_oral_willing"         # Bool: True if an individual is willing to use oral PrEP
 PREP_CAB_WILLING = "prep_cab_willing"           # Bool: True if an individual is willing to use injectable Cabotegravir PrEP
 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
 PREP_ELIGIBLE = "prep_eligible"                 # Bool: True if an individual is eligible for PrEP usage
-PREP_TYPE = "prep_type"                         # None | prep.PrEPType(enum): Oral, Cabotegravir, Lenacapavir, or VaginalRing if PrEP is being used, o/w None (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
+FIRST_CAB_START_DATE = "first_cab_start_date"       # None | date: start date of first ever injectable Cab PrEP usage
+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
 

src/hivpy/data/prep.yaml

@@ -10,10 +10,20 @@ prob_risk_informed_prep: 0.05
 prob_greater_risk_informed_prep: 0.1
 prob_suspect_risk_prep: 0.5
 
+# prep preference beta
+prep_oral_pref_beta:
+  Value: [1.1, 1.3, 1.5]
+# effect of low vl on prep willingness
+vl_prevalence_prep_threshold:
+  Value: [0.005, 0.01]
+
 # starting prep
 rate_test_onprep_any: 1
 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?
 prob_prep_restart:
   Value: [0.05, 0.10, 0.20]

src/hivpy/prep.py

@@ -41,24 +41,36 @@ def __init__(self, **kwargs):
         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
 
-        self.prep_oral_pref_beta = rng.choice([1.1, 1.3, 1.5])
+        self.prep_oral_pref_beta = self.p_data.prep_oral_pref_beta.sample()
         self.prep_cab_pref_beta = self.prep_oral_pref_beta + 0.3
         self.prep_len_pref_beta = self.prep_cab_pref_beta
         self.prep_vr_pref_beta = self.prep_oral_pref_beta - 0.1
         self.vl_prevalence_affects_prep = rng.choice([True, False], p=[1/3, 2/3])
-        self.vl_prevalence_prep_threshold = rng.choice([0.005, 0.01])
+        self.vl_prevalence_prep_threshold = self.p_data.vl_prevalence_prep_threshold.sample()
 
         self.rate_test_onprep_any = self.p_data.rate_test_onprep_any
         self.prep_willing_threshold = self.p_data.prep_willing_threshold
         self.prob_test_prep_start = self.p_data.prob_test_prep_start.sample()
+        # probability of starting prep in people who are eligible, willing,
+        # and tested for HIV according to base rate of testing
+        self.prob_base_prep_start = self.p_data.prob_base_prep_start.sample()
+        # FIXME: add 4-year scale up for these probabilities
+        self.prob_oral_prep_start = self.prob_base_prep_start
+        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()
 
     def init_prep_variables(self, pop: Population):
         pop.init_variable(col.PREP_ORAL_PREF, 0)
         pop.init_variable(col.PREP_CAB_PREF, 0)
         pop.init_variable(col.PREP_LEN_PREF, 0)
         pop.init_variable(col.PREP_VR_PREF, 0)
-        pop.init_variable(col.PREP_PREF_RANKED, None)
+        pop.init_variable(col.PREP_ORAL_RANK, 0)
+        pop.init_variable(col.PREP_CAB_RANK, 0)
+        pop.init_variable(col.PREP_LEN_RANK, 0)
+        pop.init_variable(col.PREP_VR_RANK, 0)
         pop.init_variable(col.PREP_ORAL_WILLING, False)
         pop.init_variable(col.PREP_CAB_WILLING, False)
         pop.init_variable(col.PREP_LEN_WILLING, False)
@@ -67,7 +79,17 @@ def init_prep_variables(self, pop: Population):
         pop.init_variable(col.R_PREP, 1.0)
         pop.init_variable(col.PREP_ELIGIBLE, False)
         pop.init_variable(col.PREP_TYPE, None)
+        pop.init_variable(col.EVER_PREP, False)
+        pop.init_variable(col.FIRST_ORAL_START_DATE, None)
+        pop.init_variable(col.FIRST_CAB_START_DATE, None)
+        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.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.LTP_HIV_STATUS, False)
         pop.init_variable(col.LTP_HIV_DIAGNOSED, False)
         pop.init_variable(col.LTP_ON_ART, False)
@@ -108,6 +130,39 @@ def get_suspect_risk_pop(self, pop: Population):
                                    COND(col.LTP_HIV_STATUS, op.eq, True),
                                    COND(col.R_PREP, op.lt, self.prob_suspect_risk_prep)))
 
+    def get_presumed_hiv_neg_pop(self, pop: Population):
+        """
+        Return the sub-population that has been tested and is HIV positive but
+        received a false negative result.
+        """
+        false_neg_pop = pop.get_sub_pop(AND(COND(col.EVER_TESTED, op.eq, True),
+                                            COND(col.HIV_DIAGNOSED, op.eq, False),
+                                            COND(col.HIV_STATUS, op.eq, True)))
+
+        # general test sensitivity
+        eff_test_sens = pop.hiv_diagnosis.test_sens_general
+        if not pop.hiv_diagnosis.init_prep_inj_na:
+            # infected up to 3 months ago
+            recently_infected_pop = pop.get_sub_pop_intersection(
+                pop.get_sub_pop(COND(col.DATE_HIV_INFECTION, op.ge, pop.date - timedelta(months=3))), false_neg_pop)
+
+            # expand sensitivity into a list
+            eff_test_sens = [pop.hiv_diagnosis.test_sens_general] * len(false_neg_pop)
+            false_neg_list = list(false_neg_pop)
+            # find indices in false_neg_pop that correspond to people belonging to recently_infected_pop
+            common_i = [false_neg_list.index(i) for i in false_neg_list if i in recently_infected_pop]
+
+            # FIXME: is there a better way to do this?
+            for i in common_i:
+                # lower test sensitivity used to mimic more people starting prep when they have hiv
+                eff_test_sens[i] = pop.hiv_diagnosis.test_sens_primary_ab
+
+        # false negative outcomes
+        r = rng.uniform(size=len(false_neg_pop))
+        mask = r > eff_test_sens
+
+        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):
         """
         Set preference values for a specific type of PrEP and determine willingness.
@@ -162,10 +217,11 @@ def prep_willingness(self, pop: Population):
             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)
-            pop.set_present_variable(col.PREP_PREF_RANKED, pref_ranks)
-
-        # FIXME: do we need to keep track of everyone's highest PrEP preference here?
-        # having the actual ranking may be more useful depending on availability
+            # 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
@@ -184,11 +240,16 @@ def prep_willingness(self, pop: Population):
 
     def calc_prep_pref_ranks(self, oral_pref, cab_pref, len_pref, vr_pref):
         """
-        Calculates PrEP preference rankings based on all preference values.
+        Returns PrEP preference rankings based on all preference values.
         """
+        ranks = [0, 0, 0, 0]
         prefs = [oral_pref, cab_pref, len_pref, vr_pref]
+        # reverse sort preference values (position indicates rank, value indicates prep type)
         sorted_pref_indices = sorted(range(len(prefs)), key=lambda x: prefs[x], reverse=True)
-        return [[PrEPType(i) for i in sorted_pref_indices]]
+        # 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]
 
     def prep_eligibility(self, pop: Population):
         """
@@ -356,3 +417,131 @@ 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):
+        """
+        Update people starting PrEP for the first time after testing to start PrEP.
+        """
+        # only start if specific prep type has been introduced
+        if pop.date >= self.date_prep_intro[prep_type]:
+            # tested explicitly to start prep
+            starting_prep_pop = pop.get_sub_pop_intersection(
+                prep_eligible_pop, pop.get_sub_pop(COND(prep_tested_col, op.eq, True)))
+
+            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.LAST_PREP_START_DATE, pop.date, starting_prep_pop)
+                pop.set_present_variable(first_start_col, pop.date, starting_prep_pop)
+
+    def general_start_prep(self, pop: Population, prep_eligible_pop):
+        """
+        Update people starting PrEP for the first time without specifically testing to start PrEP.
+        """
+        # not tested explicitly to start any prep
+        starting_prep_pop = pop.get_sub_pop_intersection(
+           prep_eligible_pop, pop.get_sub_pop(AND(COND(col.PREP_ORAL_TESTED, op.eq, False),
+                                                  COND(col.PREP_CAB_TESTED, op.eq, False),
+                                                  COND(col.PREP_LEN_TESTED, op.eq, False),
+                                                  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)
+
+            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)
+
+            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))))
+
+            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 calc_willing_start_prep(self, oral_pref, cab_pref, len_pref, vr_pref,
+                                oral_willing, cab_willing, len_willing, vr_willing, 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:
+            starting = r < self.prob_oral_prep_start
+        elif PrEPType(starting_prep) is PrEPType.Cabotegravir:
+            starting = r < self.prob_cab_prep_start
+        elif PrEPType(starting_prep) is PrEPType.Lenacapavir:
+            starting = r < self.prob_len_prep_start
+        elif PrEPType(starting_prep) is PrEPType.VaginalRing:
+            starting = r < self.prob_vr_prep_start
+        prep = [starting_prep if s else None for s in starting]
+
+        return prep
+
+    def start_prep(self, pop: Population):
+        """
+        Update PrEP usage for people starting PrEP 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_union(
+                pop.get_sub_pop(COND(col.HIV_STATUS, op.eq, False)), self.get_presumed_hiv_neg_pop(pop)))
+
+        # starting oral prep after testing
+        self.tested_start_prep(
+            pop, starting_prep_pop, PrEPType.Oral, col.PREP_ORAL_TESTED, col.FIRST_ORAL_START_DATE)
+        # 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)
+        # 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)
+        # starting vr prep after testing
+        self.tested_start_prep(
+            pop, starting_prep_pop, PrEPType.VaginalRing, col.PREP_VR_TESTED, col.FIRST_VR_START_DATE)
+
+        # not tested explicitly to start prep
+        self.general_start_prep(pop, starting_prep_pop)
+
+    def prep_usage(self, pop: Population):
+        """
+        Update PrEP usage for people starting, restarting, switching, and stopping PrEP.
+        """
+        # starting prep for the first time
+        self.start_prep(pop)