In mrc-ide/eppasm: Age-structured EPP Model for HIV Epidemic Estimates

knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  comment = "#>"
)

This vignette develops the workflow for fitting the EPP-ASM model and generating outputs to provide as Spectrum inputs.

## eppasm version >= 0.5.2
## devtool::install_github("mrc-ide/eppasm@dev")
## library(eppasm)
devtools::load_all()

Preparing EPP-ASM inputs

pjnz <- system.file("extdata/testpjnz", "Botswana2018.PJNZ", package="eppasm")
bw <- prepare_spec_fit(pjnz, proj.end=2021.5)

Fitting the EPP-ASM model

bwfit <- list()

bwfit$Urban <- fitmod(bw$Urban,
                      eppmod = "rlogistic_rw", rw_start = 2005, fitincrr = FALSE,
                      B0=1e4, B=1e3, opt_iter = 1:2*5, number_k=50)
bwfit$Rural <- fitmod(bw$Rural,
                      eppmod = "rlogistic_rw", rw_start = 2005, fitincrr = FALSE,
                      B0=1e4, B=1e3, opt_iter = 1:2*5, number_k=50)

When fitting, the random-walk based models only simulate through the end of the data period. The extend_projection() function extends the random walk for $r(t)$ through the end of the projection period.

bwfit <- lapply(bwfit[1], extend_projection, proj_years = 52)

Pooling EPP subpopulation results

The function aggr_specfit() This involves simulating the model for all resamples in each subregion and summing the following pop, hivpop, and artpop arrays for each of the 3000 resamples to generate 3000 national outputs.

bwaggr <- aggr_specfit(bwfit)

Generating outputs

Now generate outputs for the following by age, sex, and year:

• Total population
• HIV positive population
• New infections
• Non-HIV deaths
• HIV deaths

This is all a bit manual right now. In due course, will write better functions to assist with this.

library(magrittr)
library(data.table)

dimnm <- list(age = 15:80, sex = c("male", "female"), year = 1970:2021)

totpop <- lapply(bwaggr, apply, c(1, 2, 4), sum) %>%
  abind::abind(., rev.along=0, use.dnns=TRUE,
               new.names = c(dimnm, list(sampleid = seq_along(.)))) %>%
  melt %>%
  data.table(., outcome = "totpop")

hivpop <- lapply(bwaggr, function(x) x[ , , 2, ]) %>%
  abind::abind(., rev.along=0, use.dnns=TRUE,
               new.names = c(dimnm, list(sampleid = seq_along(.)))) %>%
  melt %>%
  data.table(., outcome = "hivpop")

infections <- lapply(bwaggr, attr, "infections") %>%
  abind::abind(., rev.along=0, use.dnns=TRUE,
               new.names = c(dimnm, list(sampleid = seq_along(.)))) %>%
  melt %>%
  data.table(., outcome = "infections")

natdeaths <- lapply(bwaggr, attr, "natdeaths") %>%
  abind::abind(., rev.along=0, use.dnns=TRUE,
               new.names = c(dimnm, list(sampleid = seq_along(.)))) %>%
  melt %>%
  data.table(., outcome = "natdeaths")

hivdeaths <- lapply(bwaggr, attr, "hivdeaths") %>%
  abind::abind(., rev.along=0, use.dnns=TRUE,
               new.names = c(dimnm, list(sampleid = seq_along(.)))) %>%
  melt %>%
  data.table(., outcome = "hivdeaths")

outputs_age <- rbind(totpop, hivpop, infections, natdeaths, hivdeaths)

outputs_age$outcome <- factor(outputs_age$outcome, c("totpop", "hivpop", "infections", "natdeaths", "hivdeaths"))

dcast(outputs_age[sampleid == 1], year + sex + age ~ outcome) %>%
  write.csv("bwaggr_ageoutputs_single-resample_v2.csv", row.names=FALSE)

dcast(outputs_age[ , .(value = stats::median(value)), .(year, sex, age, outcome)],
      year + sex + age ~ outcome) %>%
  write.csv("bwaggr_ageoutputs_median_v2.csv", row.names=FALSE)

dcast(outputs_age[ , .(value = mean(value)), .(year, sex, age, outcome)],
      year + sex + age ~ outcome) %>%
  write.csv("bwaggr_ageoutputs_mean_v2.csv", row.names=FALSE)

For comparisons, the following annual aggregate outputs are generated:

• HIV prevalence among age 15-49
• HIV incidence rate among age 15-49
• ART coverage among age 15+

hivprev <- sapply(bwaggr, prev) %>%
  "dimnames<-"(list(year = 1970:2021, sampleid = seq_len(ncol(.)))) %>%
  melt %>%
  data.table(., outcome = "prev15to49")

hivincid <- sapply(bwaggr, incid) %>%
  "dimnames<-"(list(year = 1970:2021, sampleid = seq_len(ncol(.)))) %>%
  melt %>%
  data.table(., outcome = "incid15to49")

artcov <- sapply(bwaggr, artcov15plus) %>%
  "dimnames<-"(list(year = 1970:2021, sampleid = seq_len(ncol(.)))) %>%
  melt %>%
  data.table(., outcome = "artcov15plus")

outputs_aggr <- rbind(hivprev, hivincid, artcov)

outputs_aggr$outcome <- factor(outputs_aggr$outcome,
                               c("prev15to49", "incid15to49", "artcov15plus"))

dcast(outputs_aggr[sampleid == 1], year ~ outcome) %>%
  write.csv("bwaggr_aggroutputs_single-resample_v2.csv", row.names=FALSE)

dcast(outputs_aggr[ , .(value = stats::median(value)), .(year, outcome)],
      year ~ outcome) %>%
  write.csv("bwaggr_aggroutputs_median_v2.csv", row.names=FALSE)

dcast(outputs_aggr[ , .(value = mean(value)), .(year, outcome)],
      year ~ outcome) %>%
  write.csv("bwaggr_aggroutputs_mean_v2.csv", row.names=FALSE)

Parameter vectors for single resample results

Storing the paremter vectors for the single resample outputs so that we can reproduce and step through the simulations if necessary.

theta_urban1 <- dput(c(bwfit$Urban$resample[1,]))

theta_rural1 <- dput(c(bwfit$Rural$resample[1,]))

National EPP fit outputs

Reproduce the above with a national EPP fit rather than Urban / Rural fits

```r

devtools::load_all("~/Documents/Code/R/eppasmdata")

bw <- prepare_national_fit(pjnz, proj.end=2021.5) attr(bw, "eppd") <- create_eppd("Botswana")

bwfit <- fitmod(bw, eppmod = "rlogistic_rw", rw_start = 2005, fitincrr = "linincrr", B0=1e4, B=1e3, opt_iter = 1:2*5, number_k=50) bwfit <- extend_projection(bwfit, 52)

bwfit$param <- lapply(seq_len(nrow(bwfit$resample)), function(ii) fnCreateParam(bwfit$resample[ii,], bwfit$fp)) fplist <- lapply(bwfit$param, function(par) stats::update(bwfit$fp, list=par)) bwnat <- lapply(fplist, simmod)

bwnat <- sim_mod_list(bwfit)

dimnm <- list(age = 15:80, sex = c("male", "female"), year = 1970:2021)

totpop <- lapply(bwnat, apply, c(1, 2, 4), sum) %>% abind::abind(., rev.along=0, use.dnns=TRUE, new.names = c(dimnm, list(sampleid = seq_along(.)))) %>% melt %>% data.table(., outcome = "totpop")

hivpop <- lapply(bwnat, function(x) x[ , , 2, ]) %>% abind::abind(., rev.along=0, use.dnns=TRUE, new.names = c(dimnm, list(sampleid = seq_along(.)))) %>% melt %>% data.table(., outcome = "hivpop")

infections <- lapply(bwnat, attr, "infections") %>% abind::abind(., rev.along=0, use.dnns=TRUE, new.names = c(dimnm, list(sampleid = seq_along(.)))) %>% melt %>% data.table(., outcome = "infections")

natdeaths <- lapply(bwnat, attr, "natdeaths") %>% abind::abind(., rev.along=0, use.dnns=TRUE, new.names = c(dimnm, list(sampleid = seq_along(.)))) %>% melt %>% data.table(., outcome = "natdeaths")

hivdeaths <- lapply(bwnat, attr, "hivdeaths") %>% abind::abind(., rev.along=0, use.dnns=TRUE, new.names = c(dimnm, list(sampleid = seq_along(.)))) %>% melt %>% data.table(., outcome = "hivdeaths")

outputs_age <- rbind(totpop, hivpop, infections, natdeaths, hivdeaths)

outputs_age$outcome <- factor(outputs_age$outcome, c("totpop", "hivpop", "infections", "natdeaths", "hivdeaths"))

dcast(outputs_age[sampleid == 1], year + sex + age ~ outcome) %>% write.csv("bwnat_ageoutputs_single-resample_v2.csv", row.names=FALSE)

dcast(outputs_age[ , .(value = stats::median(value)), .(year, sex, age, outcome)], year + sex + age ~ outcome) %>% write.csv("bwnat_ageoutputs_median_v2.csv", row.names=FALSE)

dcast(outputs_age[ , .(value = mean(value)), .(year, sex, age, outcome)], year + sex + age ~ outcome) %>% write.csv("bwnat_ageoutputs_mean_v2.csv", row.names=FALSE)

hivprev <- sapply(bwnat, prev) %>% "dimnames<-"(list(year = 1970:2021, sampleid = seq_len(ncol(.)))) %>% melt %>% data.table(., outcome = "prev15to49")

hivincid <- sapply(bwnat, incid) %>% "dimnames<-"(list(year = 1970:2021, sampleid = seq_len(ncol(.)))) %>% melt %>% data.table(., outcome = "incid15to49")

artcov <- sapply(bwnat, artcov15plus) %>% "dimnames<-"(list(year = 1970:2021, sampleid = seq_len(ncol(.)))) %>% melt %>% data.table(., outcome = "artcov15plus")

outputs_aggr <- rbind(hivprev, hivincid, artcov)

outputs_aggr$outcome <- factor(outputs_aggr$outcome, c("prev15to49", "incid15to49", "artcov15plus"))

dcast(outputs_aggr[sampleid == 1], year ~ outcome) %>% write.csv("bwnat_aggroutputs_single-resample_v2.csv", row.names=FALSE)

dcast(outputs_aggr[ , .(value = stats::median(value)), .(year, outcome)], year ~ outcome) %>% write.csv("bwnat_aggroutputs_median_v2.csv", row.names=FALSE)

dcast(outputs_aggr[ , .(value = mean(value)), .(year, outcome)], year ~ outcome) %>% write.csv("bwnat_aggroutputs_mean_v2.csv", row.names=FALSE)

mrc-ide/eppasm documentation built on Dec. 10, 2024, 8:19 a.m.