HIV diagnosis for people in primary infection and general population

src/hivpy/column_names.py

@@ -58,8 +58,10 @@
 HIV_STATUS = "HIV_status"                       # bool: True if person is HIV positive, o/w False
 DATE_HIV_INFECTION = "date_HIV_infection"       # None | date: date of HIV infection if HIV+, o/w None
 IN_PRIMARY_INFECTION = "in_primary_infection"   # bool: True if a person contracted HIV within 3 months of the current date, o/w False
+HIV_INFECTION_GE6M = "HIV_infection_ge6m"       # bool: True is a person has been infected with HIV for 6 months or more (DUMMY)
 HIV_DIAGNOSED = "HIV_diagnosed"                 # bool: True if individual had a positive HIV test
 HIV_DIAGNOSIS_DATE = "HIV_Diagnosis_Date"       # None | datetime.date: date of HIV diagnosis (to nearest timestep) if HIV+, o/w None
+UNDER_CARE = "under_care"                       # bool: True if under care after a positive HIV diagnosis
 VIRAL_LOAD_GROUP = "viral_load_group"           # int: value 0-5 placing bounds on viral load for an HIV positive person
 VIRAL_LOAD = "viral_load"                       # float: viral load for HIV+ person
 CD4 = "cd4"                                     # None | float: CD4 count per cubic millimeter; set to None for people w/o HIV
@@ -79,6 +81,8 @@
 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)
+PREP_TYPE = "prep_type"                         # None | prep.PrEPType(enum): Oral, Cabotegravir, Lenacapavir, or VaginalRing if PrEP is being used, o/w None (DUMMY)
+PREP_JUST_STARTED = "prep_just_started"         # Bool: True if PrEP usage began this time step (DUMMY)
 
 ART_ADHERENCE = "art_adherence"                 # DUMMY
 

src/hivpy/hiv_status.py

@@ -6,6 +6,7 @@
     from .population import Population
 
 import operator as op
+from enum import IntEnum
 
 import numpy as np
 import pandas as pd
@@ -14,6 +15,14 @@
 
 from . import output
 from .common import COND, SexType, opposite_sex, rng, timedelta
+from .prep import PrEPType
+
+
+# Ab (default), PCR (RNA VL), Ag/Ab
+class HIVTestType(IntEnum):
+    Ab = 0
+    PCR = 1
+    AgAb = 2
 
 
 class HIVStatusModule:
@@ -82,14 +91,37 @@ def __init__(self):
         self.cm_mortality_factor = rng.choice([3, 5, 10])
         self.other_adc_mortality_factor = rng.choice([1.5, 2, 3])
 
+        # HIV diagnosis
+        self.hiv_test_type = HIVTestType.Ab
+        self.init_prep_inj_pcr = rng.choice([True, False], p=[0.5, 0.5])
+        self.prep_inj_pcr = rng.choice([True, False], p=[0.5, 0.5]) if self.init_prep_inj_pcr else 0
+        self.test_sens_general = 0.98
+        self.test_sens_primary_ab = rng.choice([0.5, 0.75])
+        self.test_sens_prep_inj_primary_ab = rng.choice([0, 0.1])
+        self.test_sens_prep_inj_3m_ab = rng.choice([0, 0.2])
+        self.test_sens_prep_inj_ge6m_ab = rng.choice([0.10, 0.25, 0.50])
+        self.tests_sens_prep_inj = rng.choice([0, 1, 2, 3])
+        self.test_sens_prep_inj_primary_pcr = [0.7, 0.5, 0.3, 0.2][self.tests_sens_prep_inj]
+        self.test_sens_prep_inj_3m_pcr = [0.85, 0.7, 0.5, 0.3][self.tests_sens_prep_inj]
+        self.test_sens_prep_inj_ge6m_pcr = [0.95, 0.8, 0.7, 0.5][self.tests_sens_prep_inj]
+        self.prob_loss_at_diag = rng.choice([0.02, 0.05, 0.15, 0.35, 0.50],
+                                            p=[0.60, 0.30, 0.05, 0.04, 0.01])
+        # FIXME: may be 2 or 3 if sw_art_disadv=1
+        self.sw_incr_prob_loss_at_diag = 1
+        self.higher_newp_less_engagement = rng.choice([True, False], p=[0.2, 0.8])
+        self.prob_loss_at_diag_adc_tb = rng.beta(5, 95)
+        self.prob_loss_at_diag_non_tb_who3 = rng.beta(15, 85)
+
     def init_HIV_variables(self, population: Population):
         population.init_variable(col.HIV_STATUS, False)
         population.init_variable(col.DATE_HIV_INFECTION, None)
         population.init_variable(col.IN_PRIMARY_INFECTION, False)
+        population.init_variable(col.HIV_INFECTION_GE6M, False)  # FIXME: DUMMY variable
         population.init_variable(col.CD4, 0.0)
         population.init_variable(col.MAX_CD4, 6.6 + rng.normal(0, 0.25, size=population.size))
         population.init_variable(col.HIV_DIAGNOSED, False)
         population.init_variable(col.HIV_DIAGNOSIS_DATE, None)
+        population.init_variable(col.UNDER_CARE, False)
         population.init_variable(col.VIRAL_LOAD_GROUP, None)
         population.init_variable(col.VIRAL_LOAD, 0.0)
         population.init_variable(col.X4_VIRUS, False)
@@ -106,6 +138,8 @@ def init_HIV_variables(self, population: Population):
         population.init_variable(col.SBI_DIAGNOSED, False)
         population.init_variable(col.WHO4_OTHER, False)
         population.init_variable(col.WHO4_OTHER_DIAGNOSED, False)
+        population.init_variable(col.PREP_TYPE, None)
+        population.init_variable(col.PREP_JUST_STARTED, False)
 
         self.init_resistance_mutations(population)
 
@@ -426,3 +460,177 @@ def disease_and_diagnosis(disease_col, diagnosis_col, disease_rate, diagnosis_pr
         HIV_deaths = r_death < prob_death
         self.output.record_HIV_deaths(pop, HIV_deaths)
         return HIV_deaths
+
+    def update_HIV_diagnosis(self, pop: Population):
+        """
+        Diagnose people that have been tested this time step. The default test type used
+        is Ab, but certain policy options make use of PCR (RNA VL) or Ag/Ab tests.
+        Accuracy depends on test sensitivity, PrEP usage, as well as CD4 count.
+        """
+        # tested population in primary infection
+        primary_pop = pop.get_sub_pop([(col.IN_PRIMARY_INFECTION, op.eq, True),
+                                       (col.LAST_TEST_DATE, op.eq, pop.date),
+                                       (col.HIV_DIAGNOSED, op.eq, False)])
+
+        if len(primary_pop) > 0:
+            # primary infection diagnosis outcomes
+            diagnosed = pop.transform_group([col.PREP_TYPE, col.PREP_JUST_STARTED],
+                                            self.calc_primary_diag_outcomes, sub_pop=primary_pop)
+            # set outcomes
+            pop.set_present_variable(col.HIV_DIAGNOSED, diagnosed, primary_pop)
+            pop.set_present_variable(col.HIV_DIAGNOSIS_DATE, pop.date,
+                                     sub_pop=pop.apply_bool_mask(diagnosed, primary_pop))
+
+            # some people lost at diagnosis
+            lost = pop.transform_group([col.SEX_WORKER], self.calc_primary_loss_at_diag,
+                                       sub_pop=pop.apply_bool_mask(diagnosed, primary_pop))
+            pop.set_present_variable(col.UNDER_CARE, ~lost, sub_pop=pop.apply_bool_mask(diagnosed, primary_pop))
+
+        # remaining tested general population
+        general_pop = pop.get_sub_pop([(col.IN_PRIMARY_INFECTION, op.eq, False),
+                                       (col.LAST_TEST_DATE, op.eq, pop.date),
+                                       (col.HIV_DIAGNOSED, op.eq, False)])
+
+        if len(general_pop) > 0:
+            # general diagnosis outcomes
+            diagnosed = pop.transform_group([col.PREP_TYPE, col.HIV_INFECTION_GE6M],
+                                            self.calc_general_diag_outcomes, sub_pop=general_pop)
+            # set outcomes
+            pop.set_present_variable(col.HIV_DIAGNOSED, diagnosed, general_pop)
+            pop.set_present_variable(col.HIV_DIAGNOSIS_DATE, pop.date,
+                                     sub_pop=pop.apply_bool_mask(diagnosed, general_pop))
+
+            # FIXME: should also include onart_tm1 and may need to be affected by date_most_recent_tb
+            # some people lost at diagnosis
+            lost = pop.transform_group([col.SEX_WORKER, col.NUM_PARTNERS, col.ADC, col.TB, col.NON_TB_WHO3],
+                                       self.calc_general_loss_at_diag,
+                                       sub_pop=pop.apply_bool_mask(diagnosed, general_pop))
+            pop.set_present_variable(col.UNDER_CARE, ~lost, sub_pop=pop.apply_bool_mask(diagnosed, general_pop))
+
+    def calc_prob_primary_diag(self, prep_type, prep_just_started):
+        """
+        Calculates the probability of an individual in primary infection getting
+        diagnosed with HIV based on test sensitivity and injectable PrEP usage.
+        """
+        eff_test_sens_primary = 0
+        prep_inj = prep_type == PrEPType.Cabotegravir or prep_type == PrEPType.Lenacapavir
+        # default Ab test type
+        if self.hiv_test_type == HIVTestType.Ab:
+            # injectable PrEP started before this time step
+            if prep_inj and not prep_just_started:
+                eff_test_sens_primary = self.test_sens_prep_inj_primary_ab
+            else:
+                eff_test_sens_primary = self.test_sens_primary_ab
+        # PCR test type
+        elif self.hiv_test_type == HIVTestType.PCR:
+            # injectable PrEP started before this time step
+            if prep_inj and not prep_just_started:
+                eff_test_sens_primary = self.test_sens_prep_inj_primary_pcr
+            else:
+                eff_test_sens_primary = 0.86
+        # Ag/Ab test type
+        elif self.hiv_test_type == HIVTestType.AgAb:
+            # injectable PrEP started before this time step
+            if prep_inj and not prep_just_started:
+                eff_test_sens_primary = 0
+            else:
+                eff_test_sens_primary = 0.75
+
+        return eff_test_sens_primary
+
+    def calc_primary_diag_outcomes(self, prep_type, prep_just_started, size):
+        """
+        Uses HIV test sensitivity and injectable PrEP usage to return
+        primary infection diagnosis outcomes.
+        """
+        prob_diag = self.calc_prob_primary_diag(prep_type, prep_just_started)
+        # outcomes
+        r = rng.uniform(size=size)
+        diagnosed = r < prob_diag
+
+        return diagnosed
+
+    def calc_prob_general_diag(self, prep_type, hiv_infection_ge6m):
+        """
+        Calculates the probability of an individual not in primary infection getting diagnosed
+        with HIV based on test sensitivity, injectable PrEP usage, and infection duration.
+        """
+        eff_test_sens_general = self.test_sens_general
+        # FIXME: does injectable use timing matter for general diagnosis?
+        # injectable PrEP in current use
+        if prep_type == PrEPType.Cabotegravir or prep_type == PrEPType.Lenacapavir:
+            if self.prep_inj_pcr:
+                # infected for 6 months or more
+                if hiv_infection_ge6m:
+                    eff_test_sens_general = self.test_sens_prep_inj_ge6m_pcr
+                else:
+                    eff_test_sens_general = self.test_sens_prep_inj_3m_pcr
+            else:
+                # infected for 6 months or more
+                if hiv_infection_ge6m:
+                    eff_test_sens_general = self.test_sens_prep_inj_ge6m_ab
+                else:
+                    eff_test_sens_general = self.test_sens_prep_inj_3m_ab
+
+        return eff_test_sens_general
+
+    def calc_general_diag_outcomes(self, prep_type, hiv_infection_ge6m, size):
+        """
+        Uses HIV test sensitivity, injectable PrEP usage, and infection duration
+        to return general diagnosis outcomes.
+        """
+        prob_diag = self.calc_prob_general_diag(prep_type, hiv_infection_ge6m)
+        # outcomes
+        r = rng.uniform(size=size)
+        diagnosed = r < prob_diag
+
+        return diagnosed
+
+    def calc_prob_loss_at_diag(self, sex_worker):
+        """
+        Calculates the generic probability of an individual diagnosed with HIV
+        exiting care after diagnosis based on sex worker status.
+        """
+        # FIXME: may need to be affected by lower future ART coverage and/or decr_prob_loss_at_diag_year_i
+        eff_prob_loss_at_diag = self.prob_loss_at_diag
+        if sex_worker:
+            # FIXME: use eff_sw_incr_prob_loss_at_diag after introducing ART and SW programs
+            eff_prob_loss_at_diag = min(1, eff_prob_loss_at_diag * self.sw_incr_prob_loss_at_diag)
+
+        return eff_prob_loss_at_diag
+
+    def calc_primary_loss_at_diag(self, sex_worker, size):
+        """
+        Uses sex worker status in individuals in primary infection after a
+        positive HIV diagnosis to return loss of care outcomes.
+        """
+        prob_loss = self.calc_prob_loss_at_diag(sex_worker)
+        # outcomes
+        r = rng.uniform(size=size)
+        lost = r < prob_loss
+
+        return lost
+
+    def calc_general_loss_at_diag(self, sex_worker, num_stp, adc, tb, non_tb_who3, size):
+        """
+        Uses sex worker, ADC, TB, and non-TB WHO3 status and number of short term partners in
+        individuals not in primary infection after a positive HIV diagnosis
+        to return loss of care outcomes.
+        """
+        # outcomes
+        r = rng.uniform(size=size)
+        # ADC, non-TB WHO3, and TB not present
+        if not adc and not non_tb_who3 and not tb:
+            generic_prob_loss = self.calc_prob_loss_at_diag(sex_worker)
+            # people with more partners less likely to be engaged with care
+            if self.higher_newp_less_engagement and num_stp > 1:
+                generic_prob_loss *= 1.5
+            lost = r < generic_prob_loss
+        # ADC or TB present
+        elif adc or tb:
+            lost = r < self.prob_loss_at_diag_adc_tb
+        # non-TB WHO3 present
+        elif non_tb_who3:
+            lost = r < self.prob_loss_at_diag_non_tb_who3
+
+        return lost

src/hivpy/population.py

@@ -225,7 +225,7 @@ def set_variable_by_group(self, target, groups, func, use_size=True, sub_pop=Non
             self.data.loc[sub_pop, target_col] = self.transform_group(groups, func,
                                                                       use_size, sub_pop)
 
-    def transform_group(self, param_list, func, use_size=True, sub_pop=None):
+    def transform_group(self, param_list, func, use_size=True, sub_pop=None, dropna=False):
         """
         Groups the data by a list of parameters and applies a function to each grouping.
 
@@ -238,7 +238,9 @@ def transform_group(self, param_list, func, use_size=True, sub_pop=None):
         of the group as an argument. \n
         `sub_pop` is `None` by default, in which case the transform acts upon the entire dataframe.
         If `sub_pop` is defined, then it acts only on the part of the dataframe defined
-        by `data.loc[sub_pop]`.
+        by `data.loc[sub_pop]`. \n
+        `dropna` is false by default to allow for the inclusion of missing values in groups, but
+        should be set to true if missing values should instead be dropped during groupby.
         """
         # Use Dummy column to in order to enable transform method and avoid any risks to data
         param_list = list(map(lambda x: self.get_correct_column(x), param_list))
@@ -257,7 +259,7 @@ def general_func(g):
             df = self.data.loc[sub_pop]
         else:
             df = self.data
-        return df.groupby(param_list)["Dummy"].transform(general_func)
+        return df.groupby(param_list, dropna=dropna)["Dummy"].transform(general_func)
 
     def evolve(self, time_step: timedelta):
         """

src/hivpy/prep.py

@@ -0,0 +1,16 @@
+from enum import IntEnum
+
+
+class PrEPType(IntEnum):
+    Oral = 0
+    Cabotegravir = 1  # injectable
+    Lenacapavir = 2  # injectable
+    VaginalRing = 3
+
+
+class PrEPModule:
+
+    def __init__(self, **kwargs):
+        # FIXME: move these to data file
+        self.rate_test_onprep_any = 1
+        self.prep_willing_threshold = 0.2