mrc-ide · jeffeaton · Jan 2, 2025 · Jan 2, 2025 · Sep 25, 2024 · Jan 2, 2025
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: naomi
 Title: Naomi Model for Subnational HIV Estimates
-Version: 2.10.6
+Version: 2.10.9
 Authors@R:
     person(given = "Jeff",
            family = "Eaton",

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,19 @@
+# naomi 2.10.9
+
+* Remove `calendar_quarter_t5` two-year ahead projection.
+
+# naomi 2.10.8
+
+* Add national level aggregate to PEPFAR Target Setting Tool CSV. 
+  - This also adds a new column `country_or_psnu` that codes values `"Country"` or `"PSNU"`
+    accordingly. If the PSNU level is national, the value `"Country"` will be coded.
+* Remove two-year ahead projection (`*.T2`) indicators from PEPFAR Target Setting Tool output.
+
+# naomi 2.10.7
+
+* Add example datasets for 28 district, dropping the `District + Metro` level.
+  Datasets are saved in `extdata/demo-district28`.
+
 # naomi 2.10.6
 
 * Update `read_dp_art_dec31()` with new .DP file flags to ensure ART adjustment factor and ART patient reallocation counts are applied to number on ART extracted from Spectrum.

diff --git a/R/calibrate.R b/R/calibrate.R
@@ -503,11 +503,11 @@ calibrate_outputs <- function(output,
 
                   .expand(naomi_mf$calendar_quarter4, "plhiv"),
                   .expand(naomi_mf$calendar_quarter4, "plhiv_attend"),
-                  .expand(naomi_mf$calendar_quarter4, "infections"),
+                  .expand(naomi_mf$calendar_quarter4, "infections") ## ,
 
-                  .expand(naomi_mf$calendar_quarter5, "plhiv"),
-                  .expand(naomi_mf$calendar_quarter5, "plhiv_attend"),
-                  .expand(naomi_mf$calendar_quarter5, "infections")
+                  ## .expand(naomi_mf$calendar_quarter5, "plhiv"),
+                  ## .expand(naomi_mf$calendar_quarter5, "plhiv_attend"),
+                  ## .expand(naomi_mf$calendar_quarter5, "infections")
                 )
 
   byv <- c("indicator", "area_id", "sex", "age_group", "calendar_quarter")

diff --git a/R/model.R b/R/model.R
diff --git a/R/outputs.R b/R/outputs.R
@@ -114,10 +114,10 @@ extract_indicators <- function(naomi_fit, naomi_mf, na.rm = FALSE) {
                      "infections_t4_out" = "infections",
                      "lambda_t4_out" = "incidence")
 
-  indicators_t5 <- c("population_t5_out" = "population",
-                     "plhiv_t5_out" = "plhiv",
-                     "plhiv_attend_t5_out" = "plhiv_attend",
-                     "infections_t5_out" = "infections")
+  ## indicators_t5 <- c("population_t5_out" = "population",
+  ##                    "plhiv_t5_out" = "plhiv",
+  ##                    "plhiv_attend_t5_out" = "plhiv_attend",
+  ##                    "infections_t5_out" = "infections")
 
   if (naomi_mf$output_aware_plhiv) {
 
@@ -145,14 +145,14 @@ extract_indicators <- function(naomi_fit, naomi_mf, na.rm = FALSE) {
   indicator_est_t2 <- Map(get_est, names(indicators_t2), indicators_t2, naomi_mf$calendar_quarter2)
   indicator_est_t3 <- Map(get_est, names(indicators_t3), indicators_t3, naomi_mf$calendar_quarter3)
   indicator_est_t4 <- Map(get_est, names(indicators_t4), indicators_t4, naomi_mf$calendar_quarter4)
-  indicator_est_t5 <- Map(get_est, names(indicators_t5), indicators_t5, naomi_mf$calendar_quarter5)
+  ## indicator_est_t5 <- Map(get_est, names(indicators_t5), indicators_t5, naomi_mf$calendar_quarter5)
 
 
   indicator_est_t1 <- dplyr::bind_rows(indicator_est_t1)
   indicator_est_t2 <- dplyr::bind_rows(indicator_est_t2)
   indicator_est_t3 <- dplyr::bind_rows(indicator_est_t3)
   indicator_est_t4 <- dplyr::bind_rows(indicator_est_t4)
-  indicator_est_t5 <- dplyr::bind_rows(indicator_est_t5)
+  ## indicator_est_t5 <- dplyr::bind_rows(indicator_est_t5)
 
   indicators_anc_t1 <- c("anc_clients_t1_out" = "anc_clients",
                          "anc_plhiv_t1_out" = "anc_plhiv",
@@ -223,8 +223,8 @@ 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
+                  indicator_anc_est_t4 ## ,
+                  ## indicator_est_t5
                 )
 
   dplyr::select(out, names(naomi_mf$mf_out),
@@ -964,8 +964,8 @@ save_output <- function(filename, dir,
     c(naomi_output$fit$model_options$calendar_quarter_t1,
       naomi_output$fit$model_options$calendar_quarter_t2,
       naomi_output$fit$model_options$calendar_quarter_t3,
-      naomi_output$fit$model_options$calendar_quarter_t4,
-      naomi_output$fit$model_options$calendar_quarter_t5))
+      naomi_output$fit$model_options$calendar_quarter_t4)) ## ,
+      ## naomi_output$fit$model_options$calendar_quarter_t5))
   naomi_output$meta_period <- meta_period
 
   if (with_labels) {

diff --git a/R/pepfar-datapack.R b/R/pepfar-datapack.R
@@ -74,6 +74,10 @@ build_datapack_output <- function(naomi_output, psnu_level, dmppt2_output) {
     warning("PSNU level ", psnu_level, " not included in model outputs.")
   }
 
+  ## PEPFAR Target Setting Tool 2025: select both PSNU level and national aggregates
+  ## Assume that national aggregate is level 0
+  datapack_output_levels <- c(0L, psnu_level)
+
   datapack_indicator_map$calendar_quarter <- naomi_output$meta_period$calendar_quarter[datapack_indicator_map$time]
 
   datapack_indicator_map <- datapack_indicator_map %>%
@@ -120,9 +124,10 @@ build_datapack_output <- function(naomi_output, psnu_level, dmppt2_output) {
 
   dat <- indicators %>%
     dplyr::rename(sex_naomi = sex) %>%
-    dplyr::semi_join(
+    dplyr::inner_join(
       naomi_output$meta_area %>%
-        dplyr::filter(area_level == psnu_level),
+        dplyr::filter(area_level %in% datapack_output_levels) %>%
+        dplyr::select(area_id, area_level),
       by = "area_id"
     ) %>%
     dplyr::left_join(
@@ -142,6 +147,7 @@ build_datapack_output <- function(naomi_output, psnu_level, dmppt2_output) {
     ) %>%
     dplyr::transmute(
       area_id,
+      area_level,
       indicator,
       indicator_sort_order,
       sex_naomi,
@@ -156,7 +162,7 @@ build_datapack_output <- function(naomi_output, psnu_level, dmppt2_output) {
     dplyr::rename(age_sex_rse = rse) %>%
     dplyr::left_join(
       dplyr::filter(dat, age_group %in% c("Y000_999", "Y015_049")) %>%
-        dplyr::select(-indicator_sort_order, -age_group, -sex_naomi, -value) %>%
+        dplyr::select(-area_level, -indicator_sort_order, -age_group, -sex_naomi, -value) %>%
         dplyr::rename(district_rse = rse),
       by = c("area_id", "indicator", "calendar_quarter")
     ) %>%
@@ -165,7 +171,7 @@ build_datapack_output <- function(naomi_output, psnu_level, dmppt2_output) {
         dplyr::select(area_name, area_id),
       by = "area_id"
     ) %>%
-    dplyr::arrange(calendar_quarter, indicator_sort_order, area_id, sex_naomi, age_group)
+    dplyr::arrange(calendar_quarter, indicator_sort_order, area_level, area_id, sex_naomi, age_group)
 
 
   dat$district_rse[is.na(dat$district_rse) & dat$indicator %in% c("circ_new", "circ_ever")] <- 0.0
@@ -183,6 +189,9 @@ build_datapack_output <- function(naomi_output, psnu_level, dmppt2_output) {
   dat <- dplyr::left_join(dat, psnu_map, by = "area_id")
   dat$psnu <- ifelse(is.na(dat$map_name), dat$area_name, dat$map_name)
 
+  ## Recode area_level as "Country" or "PSNU"
+  dat$country_or_psnu <- ifelse(dat$area_level == 0, "Country", "PSNU")
+
   dat %>%
     dplyr::select(
       psnu,
@@ -197,18 +206,20 @@ build_datapack_output <- function(naomi_output, psnu_level, dmppt2_output) {
       calendar_quarter,
       value,
       age_sex_rse,
-      district_rse
+      district_rse,
+      country_or_psnu
     )
 }
 
 build_datapack_metadata <- function(naomi_output, ids) {
   cqs <- c(naomi_output$fit$model_options$calendar_quarter_t1,
            naomi_output$fit$model_options$calendar_quarter_t2,
            naomi_output$fit$model_options$calendar_quarter_t3,
-           naomi_output$fit$model_options$calendar_quarter_t4,
-           naomi_output$fit$model_options$calendar_quarter_t5)
+           naomi_output$fit$model_options$calendar_quarter_t4) ## ,
+           ## naomi_output$fit$model_options$calendar_quarter_t5)
   meta_period <- data.frame(
-    c("Time point", "t1", "t2", "t3", "t4", "t5"), c("Quarter", cqs)
+    ## c("Time point", "t1", "t2", "t3", "t4", "t5"), c("Quarter", cqs)
+    c("Time point", "t1", "t2", "t3", "t4"), c("Quarter", cqs)
   )
 
   info <- attr(naomi_output, "info")
@@ -249,7 +260,7 @@ datapack_aggregate_1to9 <- function(indicators) {
 
 
   indicators_keep <- c("plhiv", "plhiv_attend", "untreated_plhiv_attend", "infections",
-                       "population", "art_current", "art_current_residents", "aware_plhiv_num")
+                       "population", "art_current", "art_current_residents", "aware_plhiv_num", "aware_plhiv_attend")
 
   indicators1to9 <- indicators %>%
     dplyr::filter(

diff --git a/R/run-model.R b/R/run-model.R
@@ -351,13 +351,13 @@ naomi_prepare_data <- function(data, options) {
     calendar_quarter_t4 <- options$calendar_quarter_t4
   }
 
-  if(is.null(options$calendar_quarter_t5)) {
-    # If T5 is not specified, set T4 to 36-months after T3
-    t5 <- calendar_quarter_to_quarter_id(calendar_quarter_t4) + 3
-    calendar_quarter_t5 <- quarter_id_to_calendar_quarter(t5)
-  } else{
-    calendar_quarter_t5 <- options$calendar_quarter_t5
-  }
+  ## if(is.null(options$calendar_quarter_t5)) {
+  ##   # If T5 is not specified, set T4 to 36-months after T3
+  ##   t5 <- calendar_quarter_to_quarter_id(calendar_quarter_t4) + 3
+  ##   calendar_quarter_t5 <- quarter_id_to_calendar_quarter(t5)
+  ## } else{
+  ##   calendar_quarter_t5 <- options$calendar_quarter_t5
+  ## }
 
   prev_survey_ids  <- options$survey_prevalence
   recent_survey_ids <- options$survey_recently_infected
@@ -416,7 +416,7 @@ naomi_prepare_data <- function(data, options) {
     calendar_quarter2 = calendar_quarter_t2,
     calendar_quarter3 = calendar_quarter_t3,
     calendar_quarter4 = calendar_quarter_t4,
-    calendar_quarter5 = calendar_quarter_t5,
+    ## calendar_quarter5 = calendar_quarter_t5,
     spectrum_population_calibration = options$spectrum_population_calibration,
     output_aware_plhiv = as.logical(options$output_aware_plhiv),
     artattend = as.logical(options$artattend),

diff --git a/R/tmb-model-r-implementation.R b/R/tmb-model-r-implementation.R
@@ -821,55 +821,55 @@ naomi_objective_function_r <- function(d, p) {
                     anc_alpha_t4_out               = anc_alpha_t4_out)
 
 
-  ## Projection to time 5
+  ## ## Projection to time 5
 
-  mu_alpha_t5 <- mu_alpha_t4 + d$logit_alpha_t4t5_offset
-  alpha_t5 <- plogis(mu_alpha_t5)
+  ## mu_alpha_t5 <- mu_alpha_t4 + d$logit_alpha_t4t5_offset
+  ## alpha_t5 <- plogis(mu_alpha_t5)
 
-  infections_adult_t4t5 <- lambda_adult_t4 * (d$population_t4 - plhiv_t4)
-  plhiv_t5 <- as.vector(d$Lproj_hivpop_t4t5 %*% plhiv_t4 +
-                          d$Lproj_incid_t4t5 %*% infections_adult_t4t5 +
-                          d$Lproj_paed_t4t5 %*% plhiv_t4)
+  ## infections_adult_t4t5 <- lambda_adult_t4 * (d$population_t4 - plhiv_t4)
+  ## plhiv_t5 <- as.vector(d$Lproj_hivpop_t4t5 %*% plhiv_t4 +
+  ##                         d$Lproj_incid_t4t5 %*% infections_adult_t4t5 +
+  ##                         d$Lproj_paed_t4t5 %*% plhiv_t4)
 
-  rho_t5 <- plhiv_t5 / d$population_t5
-  prop_art_t5 <- rho_t5 * alpha_t5
-  artnum_t5 <- d$population_t5 * prop_art_t5
+  ## rho_t5 <- plhiv_t5 / d$population_t5
+  ## prop_art_t5 <- rho_t5 * alpha_t5
+  ## artnum_t5 <- d$population_t5 * prop_art_t5
 
-  plhiv_15to49_t5 <- as.vector(d$X_15to49 %*% plhiv_t5)
-  rho_15to49_t5 <- plhiv_15to49_t5 / as.vector(d$X_15to49 %*% d$population_t5)
-  alpha_15to49_t5 <- as.vector(d$X_15to49 %*% artnum_t5) / plhiv_15to49_t5
+  ## plhiv_15to49_t5 <- as.vector(d$X_15to49 %*% plhiv_t5)
+  ## rho_15to49_t5 <- plhiv_15to49_t5 / as.vector(d$X_15to49 %*% d$population_t5)
+  ## alpha_15to49_t5 <- as.vector(d$X_15to49 %*% artnum_t5) / plhiv_15to49_t5
 
-  mu_lambda_t5 <- d$X_lambda %*% p$beta_lambda +
-    d$log_lambda_t5_offset +
-    d$Z_x %*% (log(rho_15to49_t5) + log(1.0 - d$omega * alpha_15to49_t5)) +
-    d$Z_lambda_x %*% p$ui_lambda_x
+  ## mu_lambda_t5 <- d$X_lambda %*% p$beta_lambda +
+  ##   d$log_lambda_t5_offset +
+  ##   d$Z_x %*% (log(rho_15to49_t5) + log(1.0 - d$omega * alpha_15to49_t5)) +
+  ##   d$Z_lambda_x %*% p$ui_lambda_x
 
-  lambda_adult_t5 <- exp(mu_lambda_t5)
+  ## lambda_adult_t5 <- exp(mu_lambda_t5)
 
-  ## Add paediatric incidence
-  rho_15to49f_t5 <- d$X_15to49f %*% (plogis(mu_rho) * d$population_t5) / (d$X_15to49f %*% d$population_t5)
-  lambda_paed_t5 <- as.vector(d$X_paed_lambda_ratio_t5 %*% rho_15to49f_t5)
-  lambda_t5 <- lambda_adult_t5 + lambda_paed_t5
+  ## ## Add paediatric incidence
+  ## rho_15to49f_t5 <- d$X_15to49f %*% (plogis(mu_rho) * d$population_t5) / (d$X_15to49f %*% d$population_t5)
+  ## lambda_paed_t5 <- as.vector(d$X_paed_lambda_ratio_t5 %*% rho_15to49f_t5)
+  ## lambda_t5 <- lambda_adult_t5 + lambda_paed_t5
 
-  infections_t5 <- lambda_t5 * (d$population_t5 - plhiv_t5)
+  ## infections_t5 <- lambda_t5 * (d$population_t5 - plhiv_t5)
 
-  prop_art_ij_t5 <- as.vector(d$Xart_idx %*% prop_art_t5) * as.vector(d$Xart_gamma %*% gamma_art_t2)  ## Note: using same ART attendance as T2
-  population_ij_t5 <- as.vector(d$Xart_idx %*% d$population_t5)
-  artnum_ij_t5 <- population_ij_t5 * prop_art_ij_t5
+  ## prop_art_ij_t5 <- as.vector(d$Xart_idx %*% prop_art_t5) * as.vector(d$Xart_gamma %*% gamma_art_t2)  ## Note: using same ART attendance as T2
+  ## population_ij_t5 <- as.vector(d$Xart_idx %*% d$population_t5)
+  ## artnum_ij_t5 <- population_ij_t5 * prop_art_ij_t5
 
-  population_t5_out <- as.vector(d$A_out %*% d$population_t5)
-  plhiv_t5_out <- as.vector(d$A_out %*% plhiv_t5)
+  ## population_t5_out <- as.vector(d$A_out %*% d$population_t5)
+  ## plhiv_t5_out <- as.vector(d$A_out %*% plhiv_t5)
 
-  ## Calculate number of PLHIV who would attend facility in district i
-  plhiv_attend_ij_t5 <- as.vector(d$Xart_idx %*% plhiv_t5) * as.vector(d$Xart_gamma %*% gamma_art_t2)
-  plhiv_attend_t5_out <- as.vector(d$A_out %*% (d$A_artattend_mf %*% plhiv_attend_ij_t5))
+  ## ## Calculate number of PLHIV who would attend facility in district i
+  ## plhiv_attend_ij_t5 <- as.vector(d$Xart_idx %*% plhiv_t5) * as.vector(d$Xart_gamma %*% gamma_art_t2)
+  ## plhiv_attend_t5_out <- as.vector(d$A_out %*% (d$A_artattend_mf %*% plhiv_attend_ij_t5))
 
-  infections_t5_out <- as.vector(d$A_out %*% infections_t5)
+  ## infections_t5_out <- as.vector(d$A_out %*% infections_t5)
 
-  report_t5 <- list(population_t5_out              = population_t5_out,
-                    plhiv_t5_out                   = plhiv_t5_out,
-                    plhiv_attend_t5_out            = plhiv_attend_t5_out,
-                    infections_t5_out              = infections_t5_out)
+  ## report_t5 <- list(population_t5_out              = population_t5_out,
+  ##                   plhiv_t5_out                   = plhiv_t5_out,
+  ##                   plhiv_attend_t5_out            = plhiv_attend_t5_out,
+  ##                   infections_t5_out              = infections_t5_out)
 
   report_likelihood <- list(hhs_prev_ll           = hhs_prev_ll,
                             hhs_artcov_ll         = hhs_artcov_ll,
@@ -882,7 +882,7 @@ naomi_objective_function_r <- function(d, p) {
 
 
   v <- list(val = unname(val),
-            report = c(report_t1, report_t2, report_t3, report_t4, report_t5,
+            report = c(report_t1, report_t2, report_t3, report_t4, ## report_t5,
                        report_likelihood))
 }
 

diff --git a/R/tmb-model.R b/R/tmb-model.R
@@ -158,7 +158,7 @@ prepare_tmb_inputs <- function(naomi_data,
   X_paed_lambda_ratio_t2 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t2, df)
   X_paed_lambda_ratio_t3 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t3, df)
   X_paed_lambda_ratio_t4 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t4, df)
-  X_paed_lambda_ratio_t5 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t5, df)
+  ## X_paed_lambda_ratio_t5 <- sparse_model_matrix(~-1 + area_idf:paed_lambda_ratio_t5, df)
 
   f_rho_a <- if(all(is.na(df$rho_a_fct))) ~0 else ~0 + rho_a_fct
   f_alpha_a <- if(all(is.na(df$alpha_a_fct))) ~0 else ~0 + alpha_a_fct
@@ -345,7 +345,7 @@ prepare_tmb_inputs <- function(naomi_data,
     X_paed_lambda_ratio_t2 = X_paed_lambda_ratio_t2,
     X_paed_lambda_ratio_t3 = X_paed_lambda_ratio_t3,
     X_paed_lambda_ratio_t4 = X_paed_lambda_ratio_t4,
-    X_paed_lambda_ratio_t5 = X_paed_lambda_ratio_t5,
+    ## X_paed_lambda_ratio_t5 = X_paed_lambda_ratio_t5,
     ##
     ## Household survey input data
     x_prev = naomi_data$prev_dat$x_eff,
@@ -401,14 +401,14 @@ prepare_tmb_inputs <- function(naomi_data,
     Lproj_paed_t3t4 = naomi_data$Lproj_t3t4$Lproj_paed,
     logit_alpha_t3t4_offset = df$logit_alpha_t3t4_offset,
     log_lambda_t4_offset = df$log_lambda_t4_offset,
-    ##
-    ## Time 5 projection inputs
-    population_t5 = df$population_t5,
-    Lproj_hivpop_t4t5 = naomi_data$Lproj_t4t5$Lproj_hivpop,
-    Lproj_incid_t4t5 = naomi_data$Lproj_t4t5$Lproj_incid,
-    Lproj_paed_t4t5 = naomi_data$Lproj_t4t5$Lproj_paed,
-    logit_alpha_t4t5_offset = df$logit_alpha_t4t5_offset,
-    log_lambda_t5_offset = df$log_lambda_t5_offset,
+    ## ##
+    ## ## Time 5 projection inputs
+    ## population_t5 = df$population_t5,
+    ## Lproj_hivpop_t4t5 = naomi_data$Lproj_t4t5$Lproj_hivpop,
+    ## Lproj_incid_t4t5 = naomi_data$Lproj_t4t5$Lproj_incid,
+    ## Lproj_paed_t4t5 = naomi_data$Lproj_t4t5$Lproj_paed,
+    ## logit_alpha_t4t5_offset = df$logit_alpha_t4t5_offset,
+    ## log_lambda_t5_offset = df$log_lambda_t5_offset,
     ##
     A_out = naomi_data$A_out,
     A_anc_out = naomi_data$A_anc_out,

diff --git a/data-raw/demo-district28/README.md b/data-raw/demo-district28/README.md
@@ -0,0 +1,3 @@
+# District-level example datasets
+
+Create version of example datasets using 28 districts (area_level = 3) to avoid confusion about `District + Metro` level.
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,3 @@
		# District-level example datasets

		Create version of example datasets using 28 districts (area_level = 3) to avoid confusion about `District + Metro` level.