diff --git a/src/hivpy/column_names.py b/src/hivpy/column_names.py
index 72aa66d..b44097f 100644
--- a/src/hivpy/column_names.py
+++ b/src/hivpy/column_names.py
@@ -59,8 +59,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
@@ -81,6 +83,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
 
diff --git a/src/hivpy/data/hiv_diagnosis.yaml b/src/hivpy/data/hiv_diagnosis.yaml
new file mode 100644
index 0000000..b82b849
--- /dev/null
+++ b/src/hivpy/data/hiv_diagnosis.yaml
@@ -0,0 +1,22 @@
+# test sensitivities
+test_sens_general: 0.98
+# antibody only test
+test_sens_primary_ab:
+  Value: [0.5, 0.75]
+test_sens_prep_inj_primary_ab:
+  Value: [0, 0.1]
+test_sens_prep_inj_3m_ab:
+  Value: [0, 0.2]
+test_sens_prep_inj_ge6m_ab:
+  Value: [0.10, 0.25, 0.50]
+# nucleic acid based test
+tests_sens_prep_inj:
+  Value: [0, 1, 2, 3]
+test_sens_prep_inj_primary_na: [0.7, 0.5, 0.3, 0.2]
+test_sens_prep_inj_3m_na: [0.85, 0.7, 0.5, 0.3]
+test_sens_prep_inj_ge6m_na: [0.95, 0.8, 0.7, 0.5]
+
+# loss of care at diagnosis
+prob_loss_at_diag:
+  Value: [0.02, 0.05, 0.15, 0.35, 0.50]
+  Probability: [0.60, 0.30, 0.05, 0.04, 0.01]
diff --git a/src/hivpy/hiv_diagnosis.py b/src/hivpy/hiv_diagnosis.py
new file mode 100644
index 0000000..5695955
--- /dev/null
+++ b/src/hivpy/hiv_diagnosis.py
@@ -0,0 +1,235 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from .population import Population
+
+import importlib.resources
+import operator as op
+from enum import IntEnum
+
+import hivpy.column_names as col
+
+from .common import rng, timedelta
+from .hiv_diagnosis_data import HIVDiagnosisData
+from .prep import PrEPType
+
+
+# Ab (default), NA (RNA VL / PCR), Ag/Ab
+class HIVTestType(IntEnum):
+    Ab = 0  # antibody only assay
+    NA = 1  # nucleic acid based
+    AgAb = 2  # antigen/antibody combined assay
+
+
+class HIVDiagnosisModule:
+
+    def __init__(self, **kwargs):
+
+        # init hiv diagnosis data
+        with importlib.resources.path("hivpy.data", "hiv_diagnosis.yaml") as data_path:
+            self.hd_data = HIVDiagnosisData(data_path)
+
+        # FIXME: consider how to move more variables to data file
+        self.hiv_test_type = HIVTestType.Ab
+        self.init_prep_inj_na = rng.choice([True, False])
+        self.prep_inj_na = rng.choice([True, False]) if self.init_prep_inj_na else 0
+
+        self.test_sens_general = self.hd_data.test_sens_general
+        self.test_sens_primary_ab = self.hd_data.test_sens_primary_ab.sample()
+        self.test_sens_prep_inj_primary_ab = self.hd_data.test_sens_prep_inj_primary_ab.sample()
+        self.test_sens_prep_inj_3m_ab = self.hd_data.test_sens_prep_inj_3m_ab.sample()
+        self.test_sens_prep_inj_ge6m_ab = self.hd_data.test_sens_prep_inj_ge6m_ab.sample()
+        self.tests_sens_prep_inj = self.hd_data.tests_sens_prep_inj.sample()
+        # test_sens_prep_inj used as index for further sensitivity selection
+        self.test_sens_prep_inj_primary_na = self.hd_data.test_sens_prep_inj_primary_na[self.tests_sens_prep_inj]
+        self.test_sens_prep_inj_3m_na = self.hd_data.test_sens_prep_inj_3m_na[self.tests_sens_prep_inj]
+        self.test_sens_prep_inj_ge6m_na = self.hd_data.test_sens_prep_inj_ge6m_na[self.tests_sens_prep_inj]
+
+        self.prob_loss_at_diag = self.hd_data.prob_loss_at_diag.sample()
+        # 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 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 NA (RNA VL / PCR) 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.HIV_STATUS, op.eq, True),
+                                       (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:
+            # set infection timeframe bool
+            general_infection_dates = pop.get_variable(col.DATE_HIV_INFECTION, general_pop)
+            pop.set_present_variable(col.HIV_INFECTION_GE6M,
+                                     [pop.date]*len(general_infection_dates) - general_infection_dates >=
+                                     timedelta(months=6), general_pop)
+            # 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
+        # NA test type
+        elif self.hiv_test_type == HIVTestType.NA:
+            # 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_na
+            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_na:
+                # infected for 6 months or more
+                if hiv_infection_ge6m:
+                    eff_test_sens_general = self.test_sens_prep_inj_ge6m_na
+                else:
+                    eff_test_sens_general = self.test_sens_prep_inj_3m_na
+            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
diff --git a/src/hivpy/hiv_diagnosis_data.py b/src/hivpy/hiv_diagnosis_data.py
new file mode 100644
index 0000000..f84d5a2
--- /dev/null
+++ b/src/hivpy/hiv_diagnosis_data.py
@@ -0,0 +1,28 @@
+from hivpy.exceptions import DataLoadException
+
+from .data_reader import DataReader
+
+
+class HIVDiagnosisData(DataReader):
+    """
+    Class to hold and interpret HIV diagnosis data loaded from a yaml file.
+    """
+
+    def __init__(self, filename):
+        super().__init__(filename)
+
+        try:
+            self.test_sens_general = self.data["test_sens_general"]
+            self.test_sens_primary_ab = self._get_discrete_dist("test_sens_primary_ab")
+            self.test_sens_prep_inj_primary_ab = self._get_discrete_dist("test_sens_prep_inj_primary_ab")
+            self.test_sens_prep_inj_3m_ab = self._get_discrete_dist("test_sens_prep_inj_3m_ab")
+            self.test_sens_prep_inj_ge6m_ab = self._get_discrete_dist("test_sens_prep_inj_ge6m_ab")
+            self.tests_sens_prep_inj = self._get_discrete_dist("tests_sens_prep_inj")
+            self.test_sens_prep_inj_primary_na = self.data["test_sens_prep_inj_primary_na"]
+            self.test_sens_prep_inj_3m_na = self.data["test_sens_prep_inj_3m_na"]
+            self.test_sens_prep_inj_ge6m_na = self.data["test_sens_prep_inj_ge6m_na"]
+            self.prob_loss_at_diag = self._get_discrete_dist("prob_loss_at_diag")
+
+        except KeyError as ke:
+            print(ke.args)
+            raise DataLoadException
diff --git a/src/hivpy/hiv_status.py b/src/hivpy/hiv_status.py
index cee99c7..92bf430 100644
--- a/src/hivpy/hiv_status.py
+++ b/src/hivpy/hiv_status.py
@@ -86,10 +86,12 @@ 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)
diff --git a/src/hivpy/population.py b/src/hivpy/population.py
index fa78f4d..0146461 100644
--- a/src/hivpy/population.py
+++ b/src/hivpy/population.py
@@ -8,9 +8,11 @@
 from .circumcision import CircumcisionModule
 from .common import LogicExpr, date, timedelta
 from .demographics import DemographicsModule
+from .hiv_diagnosis import HIVDiagnosisModule
 from .hiv_status import HIVStatusModule
 from .hiv_testing import HIVTestingModule
 from .pregnancy import PregnancyModule
+from .prep import PrEPModule
 from .sexual_behaviour import SexualBehaviourModule
 
 HIV_APPEARANCE = date(1989, 1, 1)
@@ -42,6 +44,8 @@ def __init__(self, size, start_date):
         self.pregnancy = PregnancyModule()
         self.hiv_status = HIVStatusModule()
         self.hiv_testing = HIVTestingModule()
+        self.hiv_diagnosis = HIVDiagnosisModule()
+        self.prep = PrEPModule()
         self.HIV_introduced = False
         self._sample_parameters()
         self._create_population_data()
@@ -74,6 +78,7 @@ def _create_population_data(self):
         self.init_variable(col.AGE_GROUP, 0)
 
         self.hiv_status.init_HIV_variables(self)
+        self.prep.init_prep_variables(self)
         self.init_variable(col.TEST_MARK, False)
         self.init_variable(col.EVER_TESTED, False)
         self.init_variable(col.LAST_TEST_DATE, None)
@@ -226,7 +231,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.
 
@@ -239,7 +244,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.
         """
         def general_func(g):
             if len(param_list) == 1:
@@ -256,7 +263,7 @@ def general_func(g):
         else:
             df = self.data
         # Use Dummy column to in order to enable transform method and avoid any risks to 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):
         """
diff --git a/src/hivpy/prep.py b/src/hivpy/prep.py
new file mode 100644
index 0000000..574164b
--- /dev/null
+++ b/src/hivpy/prep.py
@@ -0,0 +1,29 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from .population import Population
+
+from enum import IntEnum
+
+import hivpy.column_names as col
+
+
+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
+
+    def init_prep_variables(self, pop: Population):
+        pop.init_variable(col.PREP_TYPE, None)
+        pop.init_variable(col.PREP_JUST_STARTED, False)
diff --git a/src/tests/test_hiv_diagnosis.py b/src/tests/test_hiv_diagnosis.py
new file mode 100644
index 0000000..0ea3099
--- /dev/null
+++ b/src/tests/test_hiv_diagnosis.py
@@ -0,0 +1,311 @@
+import operator as op
+from math import sqrt
+
+import hivpy.column_names as col
+from hivpy.common import date
+from hivpy.hiv_diagnosis import HIVTestType
+from hivpy.population import Population
+from hivpy.prep import PrEPType
+
+
+def test_primary_infection_diagnosis():
+    N = 1000
+    pop = Population(size=N, start_date=date(1989, 1, 1))
+    pop.data[col.IN_PRIMARY_INFECTION] = True
+    pop.data[col.DATE_HIV_INFECTION] = pop.date
+    pop.data[col.LAST_TEST_DATE] = pop.date
+    pop.data[col.HIV_DIAGNOSED] = False
+    pop.data[col.PREP_TYPE] = None
+    # test sensitivities
+    pop.hiv_diagnosis.test_sens_primary_ab = 0.50
+    test_sens_primary_na = 0.86
+    test_sens_primary_agab = 0.75
+
+    # Ab primary infection outcomes
+    pop.hiv_diagnosis.hiv_test_type = HIVTestType.Ab
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * pop.hiv_diagnosis.test_sens_primary_ab
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.test_sens_primary_ab))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # NA primary infection outcomes
+    pop.hiv_diagnosis.hiv_test_type = HIVTestType.NA
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * test_sens_primary_na
+    stdev = sqrt(mean * (1 - test_sens_primary_na))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # AgAb primary infection outcomes
+    pop.hiv_diagnosis.hiv_test_type = HIVTestType.AgAb
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * test_sens_primary_agab
+    stdev = sqrt(mean * (1 - test_sens_primary_agab))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+
+def test_primary_infection_prep_diagnosis():
+    N = 1000
+    pop = Population(size=N, start_date=date(1989, 1, 1))
+    pop.data[col.IN_PRIMARY_INFECTION] = True
+    pop.data[col.DATE_HIV_INFECTION] = pop.date
+    pop.data[col.LAST_TEST_DATE] = pop.date
+    pop.data[col.HIV_DIAGNOSED] = False
+    pop.data[col.PREP_TYPE] = PrEPType.Cabotegravir
+    pop.data[col.PREP_JUST_STARTED] = False
+    # test sensitivities
+    pop.hiv_diagnosis.test_sens_prep_inj_primary_ab = 0.1
+    pop.hiv_diagnosis.test_sens_prep_inj_primary_na = 0.3
+    test_sens_prep_inj_primary_agab = 0
+
+    # Ab + PrEP primary infection outcomes
+    pop.hiv_diagnosis.hiv_test_type = HIVTestType.Ab
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * pop.hiv_diagnosis.test_sens_prep_inj_primary_ab
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.test_sens_prep_inj_primary_ab))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # NA + PrEP primary infection outcomes
+    pop.hiv_diagnosis.hiv_test_type = HIVTestType.NA
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * pop.hiv_diagnosis.test_sens_prep_inj_primary_na
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.test_sens_prep_inj_primary_na))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # AgAb + PrEP primary infection outcomes
+    pop.hiv_diagnosis.hiv_test_type = HIVTestType.AgAb
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * test_sens_prep_inj_primary_agab
+    stdev = sqrt(mean * (1 - test_sens_prep_inj_primary_agab))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+
+def test_general_population_diagnosis():
+    N = 1000
+    pop = Population(size=N, start_date=date(1989, 1, 1))
+    pop.data[col.HIV_STATUS] = True
+    pop.data[col.IN_PRIMARY_INFECTION] = False
+    pop.data[col.DATE_HIV_INFECTION] = date(1988, 1, 1)
+    pop.data[col.LAST_TEST_DATE] = pop.date
+    pop.data[col.HIV_DIAGNOSED] = False
+    pop.data[col.PREP_TYPE] = None
+
+    # general outcomes
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * pop.hiv_diagnosis.test_sens_general
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.test_sens_general))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+
+def test_general_population_prep_diagnosis():
+    N = 1000
+    pop = Population(size=N, start_date=date(1989, 1, 1))
+    pop.data[col.HIV_STATUS] = True
+    pop.data[col.IN_PRIMARY_INFECTION] = False
+    pop.data[col.LAST_TEST_DATE] = pop.date
+    pop.data[col.HIV_DIAGNOSED] = False
+    pop.data[col.PREP_TYPE] = PrEPType.Lenacapavir
+    # test sensitivities
+    pop.hiv_diagnosis.test_sens_prep_inj_3m_ab = 0.2
+    pop.hiv_diagnosis.test_sens_prep_inj_ge6m_ab = 0.5
+    pop.hiv_diagnosis.test_sens_prep_inj_3m_na = 0.7
+    pop.hiv_diagnosis.test_sens_prep_inj_ge6m_na = 0.8
+
+    # Ab + PrEP general outcomes (recent infection)
+    pop.hiv_diagnosis.prep_inj_na = False
+    pop.data[col.DATE_HIV_INFECTION] = date(1988, 9, 1)
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * pop.hiv_diagnosis.test_sens_prep_inj_3m_ab
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.test_sens_prep_inj_3m_ab))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # Ab + PrEP general outcomes (older infection)
+    pop.data[col.DATE_HIV_INFECTION] = date(1988, 6, 1)
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * pop.hiv_diagnosis.test_sens_prep_inj_ge6m_ab
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.test_sens_prep_inj_ge6m_ab))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # NA + PrEP general outcomes (recent infection)
+    pop.hiv_diagnosis.prep_inj_na = True
+    pop.data[col.DATE_HIV_INFECTION] = date(1988, 9, 1)
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * pop.hiv_diagnosis.test_sens_prep_inj_3m_na
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.test_sens_prep_inj_3m_na))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # NA + PrEP general outcomes (older infection)
+    pop.data[col.DATE_HIV_INFECTION] = date(1988, 6, 1)
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    diag_pop = len(pop.get_sub_pop([(col.HIV_DIAGNOSED, op.eq, True)]))
+    mean = N * pop.hiv_diagnosis.test_sens_prep_inj_ge6m_na
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.test_sens_prep_inj_ge6m_na))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= diag_pop <= mean + 3 * stdev
+
+
+def test_primary_loss_at_diagnosis():
+    N = 1000
+    pop = Population(size=N, start_date=date(1989, 1, 1))
+    pop.data[col.IN_PRIMARY_INFECTION] = True
+    pop.data[col.DATE_HIV_INFECTION] = pop.date
+    pop.data[col.LAST_TEST_DATE] = pop.date
+    pop.data[col.HIV_DIAGNOSED] = False
+    pop.data[col.PREP_TYPE] = None
+    pop.data[col.SEX_WORKER] = False
+    # adjust probabilities
+    pop.hiv_diagnosis.test_sens_primary_ab = 1  # diagnose everyone
+    pop.hiv_diagnosis.prob_loss_at_diag = 0.50
+    pop.hiv_diagnosis.sw_incr_prob_loss_at_diag = 1.4
+    sw_prob_loss_at_diag = pop.hiv_diagnosis.prob_loss_at_diag * pop.hiv_diagnosis.sw_incr_prob_loss_at_diag
+
+    # Ab primary infection loss of care
+    pop.hiv_diagnosis.hiv_test_type = HIVTestType.Ab
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    lost = len(pop.get_sub_pop([(col.UNDER_CARE, op.eq, False)]))
+    mean = N * pop.hiv_diagnosis.prob_loss_at_diag
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.prob_loss_at_diag))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= lost <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # Ab primary infection loss of care (sex workers)
+    pop.data[col.SEX_WORKER] = True
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    lost = len(pop.get_sub_pop([(col.UNDER_CARE, op.eq, False)]))
+    mean = N * sw_prob_loss_at_diag
+    stdev = sqrt(mean * (1 - sw_prob_loss_at_diag))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= lost <= mean + 3 * stdev
+
+
+def test_general_loss_at_diagnosis():
+    N = 1000
+    pop = Population(size=N, start_date=date(1989, 1, 1))
+    pop.data[col.HIV_STATUS] = True
+    pop.data[col.IN_PRIMARY_INFECTION] = False
+    pop.data[col.DATE_HIV_INFECTION] = date(1988, 1, 1)
+    pop.data[col.LAST_TEST_DATE] = pop.date
+    pop.data[col.HIV_DIAGNOSED] = False
+    pop.data[col.PREP_TYPE] = None
+    pop.data[col.SEX_WORKER] = False
+    pop.data[col.ADC] = False
+    pop.data[col.TB] = False
+    pop.data[col.NON_TB_WHO3] = False
+    pop.data[col.NUM_PARTNERS] = 2
+    # adjust probabilities
+    pop.hiv_diagnosis.test_sens_general = 1  # diagnose everyone
+    pop.hiv_diagnosis.prob_loss_at_diag = 0.50
+    pop.hiv_diagnosis.prob_loss_at_diag_adc_tb = 0.30
+    pop.hiv_diagnosis.prob_loss_at_diag_non_tb_who3 = 0.45
+
+    # general outcomes (less engaged)
+    pop.hiv_diagnosis.higher_newp_less_engagement = True
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    lost = len(pop.get_sub_pop([(col.UNDER_CARE, op.eq, False)]))
+    mean = N * pop.hiv_diagnosis.prob_loss_at_diag * 1.5
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.prob_loss_at_diag * 1.5))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= lost <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # general outcomes (less engaged but no stp)
+    pop.data[col.NUM_PARTNERS] = 0
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    lost = len(pop.get_sub_pop([(col.UNDER_CARE, op.eq, False)]))
+    mean = N * pop.hiv_diagnosis.prob_loss_at_diag
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.prob_loss_at_diag))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= lost <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # general outcomes (adc or tb)
+    pop.data[col.ADC] = True
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    lost = len(pop.get_sub_pop([(col.UNDER_CARE, op.eq, False)]))
+    mean = N * pop.hiv_diagnosis.prob_loss_at_diag_adc_tb
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.prob_loss_at_diag_adc_tb))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= lost <= mean + 3 * stdev
+
+    # reset diagnosis
+    pop.data[col.HIV_DIAGNOSED] = False
+    # general outcomes (non-tb who3)
+    pop.data[col.ADC] = False
+    pop.data[col.NON_TB_WHO3] = True
+    pop.hiv_diagnosis.update_HIV_diagnosis(pop)
+
+    # get stats
+    lost = len(pop.get_sub_pop([(col.UNDER_CARE, op.eq, False)]))
+    mean = N * pop.hiv_diagnosis.prob_loss_at_diag_non_tb_who3
+    stdev = sqrt(mean * (1 - pop.hiv_diagnosis.prob_loss_at_diag_non_tb_who3))
+    # check tested value is within 3 standard deviations
+    assert mean - 3 * stdev <= lost <= mean + 3 * stdev
diff --git a/src/tests/test_hiv_status.py b/src/tests/test_hiv_status.py
index 99f378f..a91e1ca 100644
--- a/src/tests/test_hiv_status.py
+++ b/src/tests/test_hiv_status.py
@@ -1,4 +1,3 @@
-import operator
 import operator as op
 
 import numpy as np
@@ -98,11 +97,11 @@ def test_hiv_initial_ages(pop_with_initial_hiv: Population):
     """
     Check that HIV is not introduced to anyone <= 15 or > 65.
     """
-    under_15s = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, operator.eq, True),
-                                                  (col.AGE, operator.le, 15)])
-    over_65s = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, operator.eq, True),
-                                                 (col.AGE, operator.gt, 65)])
-    HIV_pos = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, operator.eq, True)])
+    under_15s = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, op.eq, True),
+                                                  (col.AGE, op.le, 15)])
+    over_65s = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, op.eq, True),
+                                                 (col.AGE, op.gt, 65)])
+    HIV_pos = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, op.eq, True)])
     assert not any(under_15s)
     assert not any(over_65s)
     assert any(HIV_pos)
@@ -111,7 +110,7 @@ def test_hiv_initial_ages(pop_with_initial_hiv: Population):
 def test_hiv_update(pop_with_initial_hiv: Population):
     pd.set_option('display.max_columns', None)
     prev_status = pop_with_initial_hiv.get_variable(col.HIV_STATUS).copy()
-    initial_infected = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, operator.eq, True)])
+    initial_infected = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, op.eq, True)])
     for i in range(10):
         pop_with_initial_hiv.date += timedelta(days=30)
         pop_with_initial_hiv.hiv_status.set_viral_load_groups(pop_with_initial_hiv)
@@ -124,8 +123,8 @@ def test_hiv_update(pop_with_initial_hiv: Population):
                    == pop_with_initial_hiv.date)
     print("Num new HIV+ = ", sum(new_cases))
     miracles = (~current_status) & (prev_status)
-    under_15s_idx = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, operator.eq, True),
-                                                      (col.AGE, operator.le, 15)])
+    under_15s_idx = pop_with_initial_hiv.get_sub_pop([(col.HIV_STATUS, op.eq, True),
+                                                      (col.AGE, op.le, 15)])
     assert not any(miracles)
     assert any(new_cases)
     assert not any(under_15s_idx)
@@ -337,7 +336,7 @@ def check_init_cd4_by_sex_age(pop, hivpos_subpop, hivneg_subpop, hiv_module, sex
     neg_cd4_counts = pop.get_variable(col.CD4, hivneg_subpop)
     assert np.allclose(neg_cd4_counts, 0.0)
     hiv_pop_by_sex = pop.get_sub_pop_intersection(
-        pop.get_sub_pop([(col.SEX, operator.eq, sex)]),
+        pop.get_sub_pop([(col.SEX, op.eq, sex)]),
         hivpos_subpop
     )
     cd4_counts = pop.get_variable(col.CD4, hiv_pop_by_sex)
diff --git a/tutorials/general.md b/tutorials/general.md
index b3b1bdd..244c701 100644
--- a/tutorials/general.md
+++ b/tutorials/general.md
@@ -12,6 +12,7 @@ Welcome to the HIVpy package. Below is an overview of the repository contents:
 
 - [Circumcision Tutorial](circumcision.md)
 - [HIV Testing Tutorial](hiv_testing.md)
+- [HIV Diagnosis Tutorial](hiv_diagnosis.md)
 - [Pregnancy Tutorial](pregnancy.md)
 - [Sexual Behaviour Tutorial](sexual_behaviour.md)
 
diff --git a/tutorials/hiv_diagnosis.md b/tutorials/hiv_diagnosis.md
new file mode 100644
index 0000000..36aab66
--- /dev/null
+++ b/tutorials/hiv_diagnosis.md
@@ -0,0 +1,10 @@
+## HIV Diagnosis Tutorial
+
+The HIV diagnosis module tracks the current HIV test type, test sensitivities, as well as loss of care at diagnosis. The most relevant files are listed below:
+
+- `src/hivpy/hiv_diagnosis.py` - The HIV diagnosis module.
+- `src/tests/test_hiv_diagnosis.py` - Tests for the HIV diagnosis module.
+- `src/hivpy/data/hiv_diagnosis.yaml` - HIV diagnosis data and variables.
+- `src/hivpy/hiv_diagnosis_data.py` - A class for storing data loaded from `hiv_diagnosis.yaml`.
+
+If there are any testing-related variables you would like to change before running your simulation, please change them in `hiv_diagnosis.yaml`.