Merge pull request #449 from mrc-ide/datapack-2025

PEPFAR Datapack 2025
mrc-ide · Nov 21, 2024 · 31ba516 · 31ba516
2 parents c89c12d + 1938b30
commit 31ba516
Show file tree

Hide file tree

Showing 10 changed files with 125 additions and 6 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: naomi
 Title: Naomi Model for Subnational HIV Estimates
-Version: 2.10.1
+Version: 2.10.2
 Authors@R:
     person(given = "Jeff",
            family = "Eaton",

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,10 @@
+# naomi 2.10.2
+
+* Add ANC testing outputs to T4 projection period for including in PEPFAR datapack output.
+* Rename Datapack input CSV in the output ZIP folder for 2025 to `"pepfar_datapack_indicators_2025.csv"`.
+* Add ANC testing indicators to T4 projection reprsenting the end of one year COP planning
+  period. New indicators are `PMTCT_STAT*.T` and `PMTCT_ART*.T`
+
 # naomi 2.10.1
 
 * Patch error in reading `anc_already_art` from Spectrum PJNZ file (was errantly

diff --git a/R/model.R b/R/model.R
@@ -567,6 +567,10 @@ naomi_model_frame <- function(area_merged,
                       spec_prev_t4 = prevalence,
                       spec_incid_t4 = incidence,
                       spec_artcov_t4 = art_coverage,
+                      spec_unaware_untreated_prop_t4 = unaware_untreated_prop,
+                      asfr_t4 = asfr,
+                      frr_plhiv_t4 = frr_plhiv,
+                      frr_already_art_t4 = frr_already_art                      
                     ),
              by = c("spectrum_region_code", "sex", "age_group")
            ) %>%

diff --git a/R/outputs.R b/R/outputs.R
@@ -184,18 +184,31 @@ extract_indicators <- function(naomi_fit, naomi_mf, na.rm = FALSE) {
                          "anc_rho_t3_out" = "anc_prevalence",
                          "anc_alpha_t3_out" = "anc_art_coverage")
 
+  indicators_anc_t4 <- c("anc_clients_t4_out" = "anc_clients",
+                         "anc_plhiv_t4_out" = "anc_plhiv",
+                         "anc_already_art_t4_out" = "anc_already_art",
+                         "anc_art_new_t4_out" = "anc_art_new",
+                         "anc_known_pos_t4_out" = "anc_known_pos",
+                         "anc_tested_pos_t4_out" = "anc_tested_pos",
+                         "anc_tested_neg_t4_out" = "anc_tested_neg",
+                         "anc_rho_t4_out" = "anc_prevalence",
+                         "anc_alpha_t4_out" = "anc_art_coverage")
+
 
   indicator_anc_est_t1 <- Map(get_est, names(indicators_anc_t1), indicators_anc_t1,
                               naomi_mf$calendar_quarter1, list(naomi_mf$mf_anc_out))
   indicator_anc_est_t2 <- Map(get_est, names(indicators_anc_t2), indicators_anc_t2,
                               naomi_mf$calendar_quarter2, list(naomi_mf$mf_anc_out))
   indicator_anc_est_t3 <- Map(get_est, names(indicators_anc_t3), indicators_anc_t3,
                               naomi_mf$calendar_quarter3, list(naomi_mf$mf_anc_out))
+  indicator_anc_est_t4 <- Map(get_est, names(indicators_anc_t4), indicators_anc_t4,
+                              naomi_mf$calendar_quarter4, list(naomi_mf$mf_anc_out))
 
 
   indicator_anc_est_t1 <- dplyr::bind_rows(indicator_anc_est_t1)
   indicator_anc_est_t2 <- dplyr::bind_rows(indicator_anc_est_t2)
   indicator_anc_est_t3 <- dplyr::bind_rows(indicator_anc_est_t3)
+  indicator_anc_est_t4 <- dplyr::bind_rows(indicator_anc_est_t4)
 
   mf_anc_out <- naomi_mf$mf_areas %>%
     dplyr::transmute(area_id,
@@ -210,6 +223,7 @@ extract_indicators <- function(naomi_fit, naomi_mf, na.rm = FALSE) {
                   indicator_est_t3,
                   indicator_anc_est_t3,
                   indicator_est_t4,
+                  indicator_anc_est_t4,
                   indicator_est_t5
                 )
 

diff --git a/R/pepfar-datapack.R b/R/pepfar-datapack.R
@@ -1,4 +1,4 @@
-PEPFAR_DATAPACK_FILENAME <- "pepfar_datapack_indicators_2024.csv"
+PEPFAR_DATAPACK_FILENAME <- "pepfar_datapack_indicators_2025.csv"
 
 #' Export naomi outputs to PEPFAR Data Pack format
 #'

diff --git a/R/tmb-model-r-implementation.R b/R/tmb-model-r-implementation.R
@@ -758,6 +758,26 @@ naomi_objective_function_r <- function(d, p) {
 
   infections_t4 <- lambda_t4 * (d$population_t4 - plhiv_t4)
 
+  ## Note: currently assuming same district effects parameters from t2 for t4
+  mu_anc_rho_t4 <- qlogis(rho_t4) +
+    d$logit_anc_rho_t2_offset +
+    d$X_ancrho %*% (p$beta_anc_rho + p$beta_anc_rho_t2) +
+    d$Z_ancrho_x %*% (p$ui_anc_rho_x + p$ui_anc_rho_xt)
+  mu_anc_rho_t4 <- as.vector(mu_anc_rho_t4)
+  anc_rho_t4 <- stats::plogis(mu_anc_rho_t4)
+
+  mu_anc_alpha_t4 <- mu_alpha_t4 +
+    d$logit_anc_alpha_t4_offset +
+    d$X_ancalpha %*% (p$beta_anc_alpha + p$beta_anc_alpha_t2) +
+    d$Z_ancalpha_x %*% (p$ui_anc_alpha_x + p$ui_anc_alpha_xt)
+  mu_anc_alpha_t4 <- as.vector(mu_anc_alpha_t4)
+  anc_alpha_t4 <- plogis(mu_anc_alpha_t4)
+
+  anc_clients_t4 <- d$population_t4 * exp(d$log_asfr_t4_offset + mu_asfr)
+  anc_plhiv_t4 <- anc_clients_t4 * anc_rho_t4
+  anc_already_art_t4 <- anc_plhiv_t4 * anc_alpha_t4
+
+
   prop_art_ij_t4 <- as.vector(d$Xart_idx %*% prop_art_t4) * as.vector(d$Xart_gamma %*% gamma_art_t2)  ## Note: using same ART attendance as T2
   population_ij_t4 <- as.vector(d$Xart_idx %*% d$population_t4)
   artnum_ij_t4 <- population_ij_t4 * prop_art_ij_t4
@@ -774,11 +794,31 @@ naomi_objective_function_r <- function(d, p) {
   infections_t4_out <- as.vector(d$A_out %*% infections_t4)
   lambda_t4_out <- infections_t4_out / (population_t4_out - plhiv_t4_out)
 
+  anc_clients_t4_out <- as.vector(d$A_anc_out %*% anc_clients_t4)
+  anc_plhiv_t4_out <- as.vector(d$A_anc_out %*% anc_plhiv_t4)
+  anc_already_art_t4_out <- as.vector(d$A_anc_out %*% anc_already_art_t4)
+  anc_art_new_t4_out <- anc_plhiv_t4_out - anc_already_art_t4_out
+  anc_known_pos_t4_out <- anc_already_art_t4_out
+  anc_tested_pos_t4_out <- anc_plhiv_t4_out - anc_known_pos_t4_out
+  anc_tested_neg_t4_out <- anc_clients_t4_out - anc_plhiv_t4_out
+
+  anc_rho_t4_out <- anc_plhiv_t4_out / anc_clients_t4_out
+  anc_alpha_t4_out <- anc_already_art_t4_out / anc_plhiv_t4_out
+
   report_t4 <- list(population_t4_out              = population_t4_out,
                     plhiv_t4_out                   = plhiv_t4_out,
                     plhiv_attend_t4_out            = plhiv_attend_t4_out,
                     lambda_t4_out                  = lambda_t4_out,
-                    infections_t4_out              = infections_t4_out)
+                    infections_t4_out              = infections_t4_out,
+                    anc_clients_t4_out             = anc_clients_t4_out,
+                    anc_plhiv_t4_out               = anc_plhiv_t4_out,
+                    anc_already_art_t4_out         = anc_already_art_t4_out,
+                    anc_art_new_t4_out             = anc_art_new_t4_out,
+                    anc_known_pos_t4_out           = anc_known_pos_t4_out,
+                    anc_tested_pos_t4_out          = anc_tested_pos_t4_out,
+                    anc_tested_neg_t4_out          = anc_tested_neg_t4_out,
+                    anc_rho_t4_out                 = anc_rho_t4_out,
+                    anc_alpha_t4_out               = anc_alpha_t4_out)
 
 
   ## Projection to time 5

diff --git a/R/tmb-model.R b/R/tmb-model.R
@@ -294,12 +294,15 @@ prepare_tmb_inputs <- function(naomi_data,
     log_asfr_t1_offset = log(df$asfr_t1),
     log_asfr_t2_offset = log(df$asfr_t2),
     log_asfr_t3_offset = log(df$asfr_t3),
+    log_asfr_t4_offset = log(df$asfr_t4),
     logit_anc_rho_t1_offset = log(df$frr_plhiv_t1),
     logit_anc_rho_t2_offset = log(df$frr_plhiv_t2),
     logit_anc_rho_t3_offset = log(df$frr_plhiv_t3),
+    logit_anc_rho_t4_offset = log(df$frr_plhiv_t4),
     logit_anc_alpha_t1_offset = log(df$frr_already_art_t1),
     logit_anc_alpha_t2_offset = log(df$frr_already_art_t2),
     logit_anc_alpha_t3_offset = log(df$frr_already_art_t3),
+    logit_anc_alpha_t4_offset = log(df$frr_already_art_t4),
     ##
     logit_rho_offset = naomi_data$mf_model$logit_rho_offset * naomi_data$mf_model$bin_rho_model,
     logit_alpha_offset = naomi_data$mf_model$logit_alpha_offset,

diff --git a/inst/datapack/datapack_indicator_mapping.csv b/inst/datapack/datapack_indicator_mapping.csv
@@ -11,6 +11,13 @@ anc_clients,RM8gRoxtsNw,PMTCT_STAT_SUBNAT.D.T_1,# Pregnant Women who present at
 anc_known_pos,tAE7ZD7p9zu,PMTCT_STAT_SUBNAT.N.Known.Pos.T_1,# Pregnant Women who are known HIV-positive at ANC1,TRUE,3
 anc_tested_neg,tAE7ZD7p9zu,PMTCT_STAT_SUBNAT.N.New.Neg.T_1,# Pregnant Women who are newly tested and found HIV-negative at ANC1,TRUE,3
 anc_tested_pos,tAE7ZD7p9zu,PMTCT_STAT_SUBNAT.N.New.Pos.T_1,# Pregnant Women who are newly tested and found HIV-positive at ANC1,TRUE,3
+anc_plhiv,PMTCT_ART_SUBNAT.D.T,PMTCT_ART_SUBNAT.D.T,# HIV-positIve Pregnant Women identified at ANC1,TRUE,4
+anc_already_art,PMTCT_ART_SUBNAT.N.Already.T,PMTCT_ART_SUBNAT.N.Already.T,# HIV-positive Pregnant Women already on ART at ANC1,TRUE,4
+anc_art_new,PMTCT_ART_SUBNAT.N.New.T,PMTCT_ART_SUBNAT.N.New.T,# HIV-positive Pregnant Women newly diagnosed at ANC1 and placed on ART,TRUE,4
+anc_clients,PMTCT_STAT_SUBNAT.D.T,PMTCT_STAT_SUBNAT.D.T,# Pregnant Women who present at ANC1,TRUE,4
+anc_known_pos,PMTCT_STAT_SUBNAT.N.Known.Pos.T,PMTCT_STAT_SUBNAT.N.Known.Pos.T,# Pregnant Women who are known HIV-positive at ANC1,TRUE,4
+anc_tested_neg,PMTCT_STAT_SUBNAT.N.New.Neg.T,PMTCT_STAT_SUBNAT.N.New.Neg.T,# Pregnant Women who are newly tested and found HIV-negative at ANC1,TRUE,4
+anc_tested_pos,PMTCT_STAT_SUBNAT.N.New.Pos.T,PMTCT_STAT_SUBNAT.N.New.Pos.T,# Pregnant Women who are newly tested and found HIV-positive at ANC1,TRUE,4
 art_current,xghQXueYJxu,TX_CURR_SUBNAT.T_1,Total # of PLHIV on ART,TRUE,3
 art_coverage,PopART.Rt.T_1,PopART.Rt.T_1,ART coverage,FALSE,3
 aware_plhiv_num,nF19GOjcnoD,DIAGNOSED_SUBNAT.T_1,Number of PLHIV aware of HIV positive status,TRUE,3

diff --git a/src/tmb.cpp b/src/tmb.cpp
@@ -104,14 +104,17 @@ Type objective_function<Type>::operator() ()
   DATA_VECTOR(log_asfr_t1_offset);
   DATA_VECTOR(log_asfr_t2_offset);
   DATA_VECTOR(log_asfr_t3_offset);
+  DATA_VECTOR(log_asfr_t4_offset);
 
   DATA_VECTOR(logit_anc_rho_t1_offset);
   DATA_VECTOR(logit_anc_rho_t2_offset);
   DATA_VECTOR(logit_anc_rho_t3_offset);
+  DATA_VECTOR(logit_anc_rho_t4_offset);
 
   DATA_VECTOR(logit_anc_alpha_t1_offset);
   DATA_VECTOR(logit_anc_alpha_t2_offset);
   DATA_VECTOR(logit_anc_alpha_t3_offset);
+  DATA_VECTOR(logit_anc_alpha_t4_offset);
 
   DATA_SPARSE_MATRIX(Z_asfr_x);
   DATA_SPARSE_MATRIX(Z_ancrho_x);
@@ -1019,6 +1022,24 @@ Type objective_function<Type>::operator() ()
 
     vector<Type> infections_t4(lambda_t4 * (population_t4 - plhiv_t4));
 
+    // Note: currently assuming same district effects parameters from t2 for t4
+    vector<Type> mu_anc_rho_t4(logit(rho_t4) +
+			       logit_anc_rho_t2_offset +
+			       X_ancrho * vector<Type>(beta_anc_rho + beta_anc_rho_t2) +
+			       Z_ancrho_x * vector<Type>(ui_anc_rho_x + ui_anc_rho_xt));
+    vector<Type> anc_rho_t4(invlogit(mu_anc_rho_t4));
+
+    vector<Type> mu_anc_alpha_t4(mu_alpha_t4 +
+				 logit_anc_alpha_t4_offset +
+				 X_ancalpha * vector<Type>(beta_anc_alpha + beta_anc_alpha_t2) +
+				 Z_ancalpha_x * vector<Type>(ui_anc_alpha_x + ui_anc_alpha_xt));
+    vector<Type> anc_alpha_t4(invlogit(mu_anc_alpha_t4));
+
+    vector<Type> anc_clients_t4(population_t4 * exp(log_asfr_t4_offset + mu_asfr));
+    vector<Type> anc_plhiv_t4(anc_clients_t4 * anc_rho_t4);
+    vector<Type> anc_already_art_t4(anc_plhiv_t4 * anc_alpha_t4);
+
+
     vector<Type> prop_art_ij_t4((Xart_idx * prop_art_t4) * (Xart_gamma * gamma_art_t2));  // Note: using same ART attendance as T2
     vector<Type> population_ij_t4(Xart_idx * population_t4);
     vector<Type> artnum_ij_t4(population_ij_t4 * prop_art_ij_t4);
@@ -1040,6 +1061,18 @@ Type objective_function<Type>::operator() ()
     vector<Type> infections_t4_out(A_out * infections_t4);
     vector<Type> lambda_t4_out(infections_t4_out / (population_t4_out - plhiv_t4_out));
 
+    vector<Type> anc_clients_t4_out(A_anc_out * anc_clients_t4);
+    vector<Type> anc_plhiv_t4_out(A_anc_out * anc_plhiv_t4);
+    vector<Type> anc_already_art_t4_out(A_anc_out * anc_already_art_t4);
+    vector<Type> anc_art_new_t4_out(anc_plhiv_t4_out - anc_already_art_t4_out);
+    vector<Type> anc_known_pos_t4_out(anc_already_art_t4_out);
+    vector<Type> anc_tested_pos_t4_out(anc_plhiv_t4_out - anc_known_pos_t4_out);
+    vector<Type> anc_tested_neg_t4_out(anc_clients_t4_out - anc_plhiv_t4_out);
+
+    vector<Type> anc_rho_t4_out(anc_plhiv_t4_out / anc_clients_t4_out);
+    vector<Type> anc_alpha_t4_out(anc_already_art_t4_out / anc_plhiv_t4_out);
+
+
     REPORT(population_t4_out);
     // REPORT(rho_t4_out);
     REPORT(plhiv_t4_out);
@@ -1053,6 +1086,17 @@ Type objective_function<Type>::operator() ()
     REPORT(lambda_t4_out);
     REPORT(infections_t4_out);
 
+    REPORT(anc_clients_t4_out);
+    REPORT(anc_plhiv_t4_out);
+    REPORT(anc_already_art_t4_out);
+    REPORT(anc_art_new_t4_out);
+    REPORT(anc_known_pos_t4_out);
+    REPORT(anc_tested_pos_t4_out);
+    REPORT(anc_tested_neg_t4_out);
+    REPORT(anc_rho_t4_out);
+    REPORT(anc_alpha_t4_out);
+
+
     // Projection to time 5
     // Only PLHIV, ART and infections calculated. No ANC or awareness indicators
 

diff --git a/tests/testthat/test-01-run-model.R b/tests/testthat/test-01-run-model.R
@@ -311,12 +311,12 @@ test_that("model run can be calibrated", {
   ## * 1 output time - 4 indicators   [NOT INCLUDED IN PLOT OUTPUTS]
   ##
   ## ANC indicators outputs
-  ## 3 = number or output times
+  ## 4 = number or output times
   ## 9 = number of ANC indicators
   ## 9 = number of areas
   ## 12 = number of ANC age groups
   expect_equal(nrow(calibrated_output_obj$output_package$indicators),
-               33 * 3 * 9 * (3 * 16 + 1 * 5 + 1 * 4) + 3 * 9 * 9 * 12)
+               33 * 3 * 9 * (3 * 16 + 1 * 5 + 1 * 4) + 4 * 9 * 9 * 12)
 
   ## Plot data output: T3 and T4 indicators not included -> fewer rows
   plot_data_output <- read_hintr_output(calibrated_output$plot_data_path)
@@ -505,7 +505,7 @@ test_that("Model can be run without .shiny90 file", {
   ## Check there is some data
   ## 11 indicators (5 fewer because missing awareness of status indicators)
   expect_equal(nrow(indicators_output$output_package$indicators),
-               33 * 3 * 9 * (3 * 11 + 1 * 5 + 1 * 4) + 3 * 9 * 9 * 12)
+               33 * 3 * 9 * (3 * 11 + 1 * 5 + 1 * 4) + 4 * 9 * 9 * 12)
 })
 
 test_that("hintr_run_model can skip validation", {