diff --git a/.Rbuildignore b/.Rbuildignore index 6b2aa7c..c5d67ae 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,5 +1,11 @@ ^README\.Rmd$ +^README\.md$ +^\.idea$ ^data-raw$ ^LICENSE\.md$ ^dev$ -^.travis.yml$ \ No newline at end of file +^.*\.Rproj$ +^\.Rproj\.user$ +^.github$ +^.ProgramDataGuidance\.md$ +^.SurveyDataGuidance\.md$ diff --git a/.github/.gitignore b/.github/.gitignore new file mode 100644 index 0000000..2d19fc7 --- /dev/null +++ b/.github/.gitignore @@ -0,0 +1 @@ +*.html diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml new file mode 100644 index 0000000..0c42ded --- /dev/null +++ b/.github/workflows/R-CMD-check.yaml @@ -0,0 +1,58 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master] + pull_request: + branches: [main, master] + +name: R-CMD-check + +jobs: + R-CMD-check: + runs-on: ${{ matrix.config.os }} + + name: ${{ matrix.config.os }} (${{ matrix.config.r }}) + + strategy: + fail-fast: false + matrix: + config: + - {os: macos-latest, r: 'release'} + - {os: windows-latest, r: 'release'} + - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} + - {os: ubuntu-latest, r: 'release'} + - {os: ubuntu-latest, r: 'oldrel-1'} + + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + R_KEEP_PKG_SOURCE: yes + + steps: + - uses: actions/checkout@v3 + + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 + with: + r-version: ${{ matrix.config.r }} + http-user-agent: ${{ matrix.config.http-user-agent }} + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::rcmdcheck + needs: check + + - uses: actions/checkout@master + with: + repository: mrc-ide/shiny90_sample_files + ## ${{ github.token }} is scoped to the current repository, so we + ## need to provide our own PAT + token: ${{ secrets.SAMPLE_FILES_GH_PAT }} + path: tests/testthat/sample_files + + - uses: r-lib/actions/check-r-package@v2 + with: + upload-snapshots: true + error-on: '"error"' diff --git a/.gitignore b/.gitignore index 4666775..ce0729a 100644 --- a/.gitignore +++ b/.gitignore @@ -7,7 +7,8 @@ .DS_Store .idea tests/testthat/sample_files +tests/testthat/Rplots.pdf *.tar.gz *.Rcheck -README.html \ No newline at end of file +README.html diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 359af69..0000000 --- a/.travis.yml +++ /dev/null @@ -1,8 +0,0 @@ -# R for travis: see documentation at https://docs.travis-ci.com/user/languages/r - -language: R -warnings_are_errors: false -before_install: - - git clone https://$github_pat@github.com/mrc-ide/shiny90_sample_files tests/testthat/sample_files -cache: - packages: true diff --git a/DESCRIPTION b/DESCRIPTION index 05c7ffc..18a3906 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: first90 -Version: 1.6.9 +Version: 1.6.10 Title: The first90 model Description: Implements the Shiny90 model for estimating progress towards the UNAIDS "first 90" target for HIV awareness of status in sub-Saharan Africa. Authors@R: @@ -25,7 +25,7 @@ Encoding: UTF-8 LazyData: true ByteCompile: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 Depends: R (>= 2.10) Imports: diff --git a/NAMESPACE b/NAMESPACE index 7b1af33..f905a74 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,5 @@ # Generated by roxygen2: do not edit by hand -export() export(add_ss_indices) export(artcov15to49) export(combine_inputs) diff --git a/NEWS.md b/NEWS.md index 4f6cd3c..39dc33b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# first90 1.6.10 + +* Qualify package names, fix some R CMD warnings and notes + # first90 1.6.9 * Bug fix: account for end-year net migration in the ART population in the first year of ART start (implemented in v1.6.0). diff --git a/R/create_hts_param.R b/R/create_hts_param.R index 5ea4a1a..6f26b90 100644 --- a/R/create_hts_param.R +++ b/R/create_hts_param.R @@ -60,16 +60,16 @@ create_hts_param <- function(theta, fp) { n_k1 <- length(knots) - 1 # we remove 1995 n_k2 <- n_k1 + length(knots_rr_dxunt) rate_f <- exp(theta[1:n_k1]) - rr_dxunt <- plogis(theta[(n_k1 + 1):n_k2]) * 8 # Range is 0-8 - rr_m <- plogis(theta[(n_k2 + 1):(n_k2 + 2)]) * 1.1 - rr_test <- 0.95 + plogis(theta[(n_k2 + 3):(n_k2 + 4)]) * 7.05 # Range is 1-8 - rr_plhiv <- 0.05 + plogis(theta[n_k2 + 5]) * (1.95 - 0.05) - rr_dxart <- plogis(theta[n_k2 + 6]) - rr_25m <- 0.1 + plogis(theta[n_k2 + 7]) * (6 - 0.1) - rr_35m <- 0.1 + plogis(theta[n_k2 + 8]) * (6 - 0.1) - rr_25f <- 0.1 + plogis(theta[n_k2 + 9]) * (6 - 0.1) - rr_35f <- 0.1 + plogis(theta[n_k2 + 10]) * (6 - 0.1) - pr_oidx <- 0.25 + (plogis(theta[n_k2 + 11]) * (1.75 - 0.25)) + rr_dxunt <- stats::plogis(theta[(n_k1 + 1):n_k2]) * 8 # Range is 0-8 + rr_m <- stats::plogis(theta[(n_k2 + 1):(n_k2 + 2)]) * 1.1 + rr_test <- 0.95 + stats::plogis(theta[(n_k2 + 3):(n_k2 + 4)]) * 7.05 # Range is 1-8 + rr_plhiv <- 0.05 + stats::plogis(theta[n_k2 + 5]) * (1.95 - 0.05) + rr_dxart <- stats::plogis(theta[n_k2 + 6]) + rr_25m <- 0.1 + stats::plogis(theta[n_k2 + 7]) * (6 - 0.1) + rr_35m <- 0.1 + stats::plogis(theta[n_k2 + 8]) * (6 - 0.1) + rr_25f <- 0.1 + stats::plogis(theta[n_k2 + 9]) * (6 - 0.1) + rr_35f <- 0.1 + stats::plogis(theta[n_k2 + 10]) * (6 - 0.1) + pr_oidx <- 0.25 + (stats::plogis(theta[n_k2 + 11]) * (1.75 - 0.25)) # Age function for males and females agefn_m <- c(rep(1, 3), rep(rr_25m, 2), rep(rr_35m, 3), rr_35m * 0.8112) # last age group is 50+ @@ -77,15 +77,15 @@ create_hts_param <- function(theta, fp) { # Time trends in testing rates fp$t_hts_start <- as.integer(1995 - fp$ss$proj_start + 1L) - base_rate_f <- approx(knots, c(0, rate_f), seq_len(fp$ss$PROJ_YEARS), rule = 2)$y + base_rate_f <- stats::approx(knots, c(0, rate_f), seq_len(fp$ss$PROJ_YEARS), rule = 2)$y knots_rr_m <- c(2005, 2012) - fp$ss$proj_start + 1L - rr_m <- approx(knots_rr_m, rr_m, seq_len(fp$ss$PROJ_YEARS), rule = 2)$y + rr_m <- stats::approx(knots_rr_m, rr_m, seq_len(fp$ss$PROJ_YEARS), rule = 2)$y base_rate_m <- base_rate_f * rr_m # Retest with time knots_rr_test <- c(2005, 2010, 2015) - fp$ss$proj_start + 1L - rate_rr_test <- exp(approx(knots_rr_test, log(c(1, rr_test)), seq_len(fp$ss$PROJ_YEARS), rule = 2)$y) + rate_rr_test <- exp(stats::approx(knots_rr_test, log(c(1, rr_test)), seq_len(fp$ss$PROJ_YEARS), rule = 2)$y) ## Testing rate for HIV negative population: ## 1: never tested @@ -147,7 +147,7 @@ create_hts_param <- function(theta, fp) { ## Relative testing among aware # Re-testing rate among PLHIV aware (not on ART) - rate_dxunt <- approx(knots_rr_dxunt, rr_dxunt, seq_len(fp$ss$PROJ_YEARS), rule = 2)$y + rate_dxunt <- stats::approx(knots_rr_dxunt, rr_dxunt, seq_len(fp$ss$PROJ_YEARS), rule = 2)$y diagn_rate_dxunt <- sweep(diagn_rate[,,,3,, drop = FALSE], 5, rate_dxunt, "*") diagn_rate[,,,3,] <- diagn_rate_dxunt @@ -162,4 +162,3 @@ create_hts_param <- function(theta, fp) { fp } - diff --git a/R/extract-pjnz.R b/R/extract-pjnz.R index b92c4d9..1f11849 100644 --- a/R/extract-pjnz.R +++ b/R/extract-pjnz.R @@ -134,9 +134,9 @@ extract_pjnz <- function(pjnz = NULL, dp_file= NULL, pjn_file = NULL){ v$art15plus_num <- v$art15plus_num * adult_artadj_factor } - v$art15plus_eligthresh <- setNames(as.numeric(dpsub(dp, "", 2, col_idx)), proj_years) + v$art15plus_eligthresh <- stats::setNames(as.numeric(dpsub(dp, "", 2, col_idx)), proj_years) - artelig_specpop <- setNames(dpsub(dp, "", 3:9, 2:6), + artelig_specpop <- stats::setNames(dpsub(dp, "", 3:9, 2:6), c("description", "pop", "elig", "percent", "year")) artelig_specpop$pop <- c("PW", "TBHIV", "DC", "FSW", "MSM", "IDU", "OTHER") artelig_specpop$elig <- as.logical(as.integer(artelig_specpop$elig)) @@ -146,8 +146,8 @@ extract_pjnz <- function(pjnz = NULL, dp_file= NULL, pjn_file = NULL){ rownames(artelig_specpop) <- artelig_specpop$pop v$artelig_specpop <- artelig_specpop - v$median_cd4init <- setNames(as.numeric(dpsub(dp, "", 2, col_idx)), proj_years) - v$art_dropout <- setNames(as.numeric(dpsub(dp, "", 2, col_idx)), proj_years) + v$median_cd4init <- stats::setNames(as.numeric(dpsub(dp, "", 2, col_idx)), proj_years) + v$art_dropout <- stats::setNames(as.numeric(dpsub(dp, "", 2, col_idx)), proj_years) v$art_alloc_method <- get_dp_art_alloc_method(dp) v$art_prop_alloc <- get_dp_art_prop_alloc(dp) @@ -763,7 +763,7 @@ get_pjn_iso3 <- function(pjn){ #' @export file_in_zip <- function(zfile, ext){ - ff <- grep(ext, unzip(zfile, list = TRUE)$Name, value = TRUE) + ff <- grep(ext, utils::unzip(zfile, list = TRUE)$Name, value = TRUE) unz(zfile, ff) } diff --git a/R/first90-package.R b/R/first90-package.R index b7c13e5..bc134b9 100644 --- a/R/first90-package.R +++ b/R/first90-package.R @@ -5,3 +5,13 @@ NULL #' @importFrom Rcpp sourceCpp ## usethis namespace: end NULL + +## Suppress R CMD check NOTEs about "no visible binding for +## global variable". These NOTEs aren't that important as we're using +## but the thousands of lines of NOTE they print does obscure real issues. +utils::globalVariables(c( + "DT", "PROJ_YEARS", "ag.idx", "agegr", "aglast.idx", "country", "f.idx", + "h.ag.span", "h.age15plus.idx", "h.fert.idx", "hAG", "hDS", "hTS", + "hiv_steps_per_year", "hivn.idx", "hivp.idx", "hivstatus", "logit", "m.idx", + "nearPD", "outcome", "p.age15plus.idx", "p.age15to49.idx", "p.fert.idx", + "pAG", "pDS", "quantile", "sex", "tot", "totpos", "year")) diff --git a/R/inputs.R b/R/inputs.R index 0ffea3b..323b9dc 100644 --- a/R/inputs.R +++ b/R/inputs.R @@ -1,5 +1,5 @@ #' Process programmatic data on number of tests -#' +#' #' @export select_prgmdata <- function(prgm_dat, cnt, age_group) { # any country with HTS data? if so, we select them @@ -24,19 +24,19 @@ select_prgmdata <- function(prgm_dat, cnt, age_group) { anc = prg.i$anc[prg.i$sex == 'female'], ancpos = prg.i$ancpos[prg.i$sex == 'female']) value_verif <- prg_dat_both$totpos[prg_dat_both$year == yr_sex[i]] - if (length(value_verif) > 0 & is.na(tot.i$totpos)) { + if (length(value_verif) > 0 & is.na(tot.i$totpos)) { tot.i$totpos <- prg_dat_both$totpos[prg_dat_both$year == yr_sex[i]] } prg_dat_sex_both <- rbind(prg_dat_sex_both, tot.i) } prg_dat <- rbind(prg_dat_both[!(prg_dat_both$year %in% yr_sex), ], prg_dat_sex_both, prg_dat_sex) - } } - + } } + if (!any(prgm_dat$country == cnt)) { ## -- UPDATE HERE -- ## * year vector needs to be extended to output results to current year - - prg_dat <- data.frame(country = cnt, + + prg_dat <- data.frame(country = cnt, year = 2010:2022, agegr = '15-99', sex = 'both', tot = NA, totpos = NA, @@ -48,65 +48,65 @@ select_prgmdata <- function(prgm_dat, cnt, age_group) { } #' Process survey data on hiv testing behaviors -#' +#' #' @export select_hts <- function(survey_hts, cnt, age_group) { require(first90) # We first select stratified by HIV status for females - dat_f <- survey_hts[survey_hts$country == cnt & (survey_hts$hivstatus != "all") & - (survey_hts$agegr %in% age_group) & survey_hts$sex == "female" & + dat_f <- survey_hts[survey_hts$country == cnt & (survey_hts$hivstatus != "all") & + (survey_hts$agegr %in% age_group) & survey_hts$sex == "female" & survey_hts$outcome == "evertest", ] # We first select stratified by HIV status for males - dat_m <- survey_hts[survey_hts$country == cnt & (survey_hts$hivstatus != "all") & - (survey_hts$agegr %in% age_group) & survey_hts$sex == "male" & + dat_m <- survey_hts[survey_hts$country == cnt & (survey_hts$hivstatus != "all") & + (survey_hts$agegr %in% age_group) & survey_hts$sex == "male" & survey_hts$outcome == "evertest", ] if (cnt == 'Equatorial Guinea') { # special case of GE b/c data for HIV entered manually (not age disaggregated) - dat_f <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "all" & + dat_f <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "all" & (survey_hts$agegr %in% age_group) & survey_hts$sex == "female" & survey_hts$outcome == "evertest", ] - - dat_m <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "all" & + + dat_m <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "all" & (survey_hts$agegr %in% age_group) & survey_hts$sex == "male" & survey_hts$outcome == "evertest", ] - - dat_f_hiv <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "positive" & + + dat_f_hiv <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "positive" & survey_hts$agegr == '15-49' & survey_hts$sex == "female" & survey_hts$outcome == "evertest", ] - - dat_m_hiv <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "positive" & + + dat_m_hiv <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "positive" & survey_hts$agegr == '15-49' & survey_hts$sex == "male" & survey_hts$outcome == "evertest", ] - + dat_f <- rbind(dat_f, dat_f_hiv) dat_m <- rbind(dat_m, dat_m_hiv) } - + # We find surveys that collected information on HIV status but for different age groups all_svy_hiv <- unique(survey_hts$surveyid[survey_hts$country == cnt & (survey_hts$hivstatus != "all") & survey_hts$outcome == "evertest"]) other_svy_hiv <- all_svy_hiv[!(all_svy_hiv %in% unique(rbind(dat_f, dat_m)$surveyid))] - + dat_both <- NULL if (length(other_svy_hiv) > 0) { - dat_f2 <- survey_hts[(survey_hts$surveyid %in% other_svy_hiv) & (survey_hts$hivstatus != "all") & - survey_hts$agegr == '15-49' & survey_hts$sex == "female" & + dat_f2 <- survey_hts[(survey_hts$surveyid %in% other_svy_hiv) & (survey_hts$hivstatus != "all") & + survey_hts$agegr == '15-49' & survey_hts$sex == "female" & survey_hts$outcome == "evertest", ] - dat_m2 <- survey_hts[(survey_hts$surveyid %in% other_svy_hiv) & (survey_hts$hivstatus != "all") & - survey_hts$agegr == '15-49' & survey_hts$sex == "male" & + dat_m2 <- survey_hts[(survey_hts$surveyid %in% other_svy_hiv) & (survey_hts$hivstatus != "all") & + survey_hts$agegr == '15-49' & survey_hts$sex == "male" & survey_hts$outcome == "evertest", ] if (dim(dat_f2)[1] > 0 | dim(dat_m2)[1] > 0) { dat_f <- rbind(dat_f, dat_f2) - dat_m <- rbind(dat_m, dat_m2) + dat_m <- rbind(dat_m, dat_m2) } else { - dat_both <- survey_hts[(survey_hts$surveyid %in% other_svy_hiv) & (survey_hts$hivstatus != "all") & - survey_hts$agegr == '15-49' & survey_hts$sex == "both" & + dat_both <- survey_hts[(survey_hts$surveyid %in% other_svy_hiv) & (survey_hts$hivstatus != "all") & + survey_hts$agegr == '15-49' & survey_hts$sex == "both" & survey_hts$outcome == "evertest", ] } } - + # We find surveys that did not collect information on HIV status all_survey <- unique(survey_hts$surveyid[survey_hts$country == cnt & survey_hts$outcome == "evertest"]) other_survey <- all_survey[!(all_survey %in% unique(rbind(dat_f, dat_m, dat_both)$surveyid))] - + dat_fall <- NULL dat_mall <- NULL @@ -115,7 +115,7 @@ select_hts <- function(survey_hts, cnt, age_group) { age_strat_f <- unique(survey_hts$agegr[survey_hts$country == cnt & survey_hts$hivstatus == "all" & survey_hts$surveyid == other_survey[i] & survey_hts$sex == "female" & survey_hts$outcome == "evertest"]) age_strat_m <- unique(survey_hts$agegr[survey_hts$country == cnt & survey_hts$hivstatus == "all" & survey_hts$surveyid == other_survey[i] & - survey_hts$sex == "male" & survey_hts$outcome == "evertest"]) + survey_hts$sex == "male" & survey_hts$outcome == "evertest"]) ## Female - If we have the correct age disagregation, we stick with it if (all(age_group %in% age_strat_f)) { dat_f_i <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "all" & survey_hts$surveyid == other_survey[i] & @@ -129,7 +129,7 @@ select_hts <- function(survey_hts, cnt, age_group) { } else { dat_f_i <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "all" & survey_hts$surveyid == other_survey[i] & (survey_hts$agegr %in% c('15-99')) & survey_hts$sex == "female" & survey_hts$outcome == "evertest",] - } + } ## Female - If we have the correct age disagregation, we stick with it if (all(age_group %in% age_strat_m)) { dat_m_i <- survey_hts[survey_hts$country == cnt & survey_hts$hivstatus == "all" & survey_hts$surveyid == other_survey[i] & @@ -154,7 +154,7 @@ select_hts <- function(survey_hts, cnt, age_group) { } #' Function to present the survey data on shinny interface -#' +#' #' @export svy_hts_interface <- function(survey_hts, cnt, age_group = c('15-24','25-34','35-49')) { # Data used in the fitting @@ -172,13 +172,13 @@ svy_hts_interface <- function(survey_hts, cnt, age_group = c('15-24','25-34','35 sex == 'both' & (outcome %in% c('aware', 'art_both', 'art_self') | hivstatus == 'positive' & outcome == 'evertest')) - + svy_plot <- rbind(svy_plot1549, svy_plotage, svy_plot90s) svy_plot <- svy_plot[order(svy_plot$year, svy_plot$sex, svy_plot$hivstatus, svy_plot$agegr), ] svy_plot$color <- "Plotting" - + svy <- rbind(svy_dat, svy_plot) - + return(svy) } @@ -201,9 +201,9 @@ add_ss_indices <- function(dat, ss) { ## Convert 15+ to 15-99 for parsing below agegr_tmp <- sub("15\\+", "15-99", df$agegr) - - df$agestart <- type.convert(sub("([0-9]+)-([0-9]+)", "\\1", agegr_tmp), as.is = TRUE) - ageend <- type.convert(sub("([0-9]+)-([0-9]+)", "\\2", agegr_tmp), as.is = TRUE)+1L + + df$agestart <- utils::type.convert(sub("([0-9]+)-([0-9]+)", "\\1", agegr_tmp), as.is = TRUE) + ageend <- utils::type.convert(sub("([0-9]+)-([0-9]+)", "\\2", agegr_tmp), as.is = TRUE)+1L df$aidx <- df$agestart - ss$AGE_START + 1L df$agspan <- ageend - df$agestart @@ -215,14 +215,14 @@ add_ss_indices <- function(dat, ss) { ## if max age is 65 or greater, assume age 50+ end_haidx[ageend >= 65] <- length(hag_breaks) + 1L - + df$hagspan <- end_haidx - df$haidx - + ## testing history is homogenous among oldest age group df$haidx[df$agestart >= max(hag_breaks) & ageend >= max(hag_breaks)] <- length(hag_breaks) df$hagspan[df$agestart >= max(hag_breaks) & ageend >= max(hag_breaks)] <- 1L - + df$sidx <- match(df$sex, c("both", "male", "female")) - 1L df$hvidx <- match(df$hivstatus, c("all", "negative", "positive")) - 1L df$yidx <- df$year - ss$proj_start + 1L @@ -242,32 +242,32 @@ add_ss_indices <- function(dat, ss) { prepare_inputs_from_extracts <- function(pjnz_in){ pjnz_aggr <- combine_inputs(pjnz_in) - + demp <- list(basepop = pjnz_aggr$totpop, Sx = pjnz_aggr$Sx, asfr = pjnz_aggr$asfr, srb = pjnz_aggr$srb, netmigr = pjnz_aggr$netmigr) - + projp <- pjnz_aggr projp$age14totpop <- projp$totpop["14",,] dimnames(projp$fert_rat)[[1]] <- gsub("(.*)-.*", "\\1", dimnames(projp$fert_rat)[[1]]) - + fp <- create_fp(projp, demp) ## Set pop adjust fp$popadjust <- pjnz_aggr$popadjust - + ## Calculate incidence input fp$infections <- pjnz_aggr$infections[fp$ss$AGE_START + fp$ss$p.age15plus.idx, , ] fp$eppmod <- "directinfections_hts" - + ## initialise to no testing fp$hts_rate <- array(0.0, c(fp$ss$hAG, fp$ss$NG, fp$ss$pDS, fp$ss$PROJ_YEARS)) fp$diagn_rate <- array(0.0, c(fp$ss$hDS, fp$ss$hAG, fp$ss$NG, 4, fp$ss$PROJ_YEARS)) - + fp$t_hts_start <- fp$ss$PROJ_YEARS+1L - + fp } @@ -285,9 +285,10 @@ prepare_inputs_from_extracts <- function(pjnz_in){ #' #' @examples #' -#' pjnzlist <- list.files("~/Documents/Data/Spectrum files/2018 final/SSA/", "CotedIvoire.*PJNZ$", full.names=TRUE, ignore.case=TRUE) +#' pjnzlist <- list.files("~/Documents/Data/Spectrum files/2018 final/SSA/", +#' "CotedIvoire.*PJNZ$", full.names=TRUE, ignore.case=TRUE) #' pjnzlist <- "~/Documents/Data/Spectrum files/2018 final/SSA/Malawi_2018_version_8.PJNZ" -#' +#' #' @export prepare_inputs <- function(pjnzlist){ @@ -300,13 +301,13 @@ add_fp <- function(v, fp){ fp[names(v)] <- v fp } - + create_paedsurv_inputs <- function(projp, ss, artcd4elig_idx) { proj_start <- ss$proj_start proj_end <- proj_start + ss$PROJ_YEARS - 1L - + val <- list() ## HIV prevalence and ART coverage among age 15 entrants @@ -337,7 +338,7 @@ create_paedsurv_inputs <- function(projp, ss, artcd4elig_idx) { ## sum to 1 in each sex and year. for(g in 1:ss$NG) for(i in 2:ss$PROJ_YEARS){ - + if((hiv14[g,i-1] - art14[g,i-1]) > 0) val$paedsurv_cd4dist[,g,i] <- hiv_noart14[,g,i-1] %*% cd4convert / (hiv14[g,i-1] - art14[g,i-1]) if(art14[g,i-1]){ @@ -359,7 +360,7 @@ create_paedsurv_inputs <- function(projp, ss, artcd4elig_idx) { .fp_targetpop <- function(totpop, ss){ v <- list() - + proj_start <- ss$proj_start proj_end <- proj_start + ss$PROJ_YEARS - 1L diff --git a/R/likelihood.R b/R/likelihood.R index 35809a1..735f9dd 100644 --- a/R/likelihood.R +++ b/R/likelihood.R @@ -1,9 +1,9 @@ #' @export ll_hts <- function(theta, fp, likdat) { - + fp <- create_hts_param(theta, fp) mod <- simmod(fp) - + val1 <- ll_evertest(mod, fp, likdat$evertest) if (!is.null(likdat$hts)) { val2 <- ll_prgdat(mod, fp, likdat$hts) @@ -11,11 +11,11 @@ ll_hts <- function(theta, fp, likdat) { if (!is.null(likdat$hts_pos)) { val3 <- ll_prgdat(mod, fp, likdat$hts_pos) } else { val3 <- 0 } - + val_prior <- lprior_hts(theta, mod, fp) val <- val1 + val2 + val3 + val_prior - + return(val) } @@ -34,11 +34,11 @@ ll_evertest <- function(mod, fp, dat) { } if (any(is.na(mu)) || any(mu < 0) || any(mu > 1)) { - llk <- log(0) + llk <- log(0) } else { - llk <- sum(dbinom(x = dat$nsuccess, size = dat$neff, prob = mu, log = TRUE)) + llk <- sum(stats::dbinom(x = dat$nsuccess, size = dat$neff, prob = mu, log = TRUE)) } - + return(llk) } @@ -47,16 +47,17 @@ ll_evertest <- function(mod, fp, dat) { ll_prgdat <- function(mod, fp, dat) { if (as.character(dat$hivstatus[1]) == 'all') { mu <- total_tests(mod, df = add_ss_indices(dat, fp$ss)) - llk <- sum(dnorm(x = dat$tot, mean = mu, sd = dat$l_est_se, log = TRUE)) - } + llk <- sum(stats::dnorm(x = dat$tot, mean = mu, sd = dat$l_est_se, log = TRUE)) + } if (as.character(dat$hivstatus[1]) == 'positive') { mu <- total_tests(mod, df = add_ss_indices(dat, fp$ss)) - llk <- sum(dnorm(x = dat$tot, mean = mu, sd = dat$l_est_se, log = TRUE)) + llk <- sum(stats::dnorm(x = dat$tot, mean = mu, sd = dat$l_est_se, log = TRUE)) } if (!(as.character(dat$hivstatus[1]) %in% c('all', 'positive'))) { - print('Error - HIV status is incorrect'); break + print('Error - HIV status is incorrect') + return() } - + return(llk) } @@ -68,56 +69,56 @@ art_constraint_penalty <- function(mod, fp, max_year = 2022) { tot_late <- apply(attr(mod, "late_diagnoses")[,,, ind_year], 4, sum) tot_untreated_pop <- apply(attr(mod, "hivpop")[,,, ind_year], 4, sum) ratio_late_per_1000_untreated <- tot_late / tot_untreated_pop * 1000 - penalty <- sum(dnorm(x = 0, mean = ratio_late_per_1000_untreated, + penalty <- sum(stats::dnorm(x = 0, mean = ratio_late_per_1000_untreated, sd = 20, log = TRUE)) return(penalty) } # Include this in ll_hts if you want to incorporate the likelihood constraint on ART. # val_art_penalty <- art_constraint_penalty(mod, fp, max_year = 2022) # val <- val1 + val2 + val3 + val_prior + val_art_penalty - + # Function to prepare the data for input in the likelihood function. #' @export prepare_hts_likdat <- function(dat_evertest, dat_prg, fp) { - + # Testing behavior data dat_evertest$est <- ifelse(is.na(dat_evertest$est), NA, dat_evertest$est) - dat_evertest <- dat_evertest[complete.cases(dat_evertest$est),] - dat_evertest$l_est <- qlogis(dat_evertest$est) + dat_evertest <- dat_evertest[stats::complete.cases(dat_evertest$est),] + dat_evertest$l_est <- stats::qlogis(dat_evertest$est) dat_evertest$l_est_se <- dat_evertest$se / (dat_evertest$est * (1 - dat_evertest$est)) - + # if using binomial, we calculate Neff from SE - dat_evertest$neff <- (dat_evertest$est * (1 - dat_evertest$est)) / dat_evertest$se^2 + dat_evertest$neff <- (dat_evertest$est * (1 - dat_evertest$est)) / dat_evertest$se^2 dat_evertest$neff <- ifelse(dat_evertest$est < 1e-5, dat_evertest$counts, dat_evertest$neff) dat_evertest$neff <- ifelse(dat_evertest$est > 0.999, dat_evertest$counts, dat_evertest$neff) dat_evertest$neff <- ifelse(dat_evertest$neff == Inf, dat_evertest$counts, dat_evertest$neff) dat_evertest$nsuccess <- round(dat_evertest$est * dat_evertest$neff, 0) dat_evertest$neff <- ifelse(round(dat_evertest$neff, 0) < 1, 1, round(dat_evertest$neff, 0)) dat_evertest <- add_ss_indices(dat_evertest, fp$ss) - + # Verifying if input data is OK if(any(dat_prg$agegr != '15-99')) { stop(print('Error, HTS program data should be for the 15+ age group only /n other age-grouping not supported at the moment'))} - + # We remove years with NA from programmatic data dat_prg_pos <- dat_prg - dat_prg_pos <- dat_prg_pos[complete.cases(dat_prg_pos$totpos),] + dat_prg_pos <- dat_prg_pos[stats::complete.cases(dat_prg_pos$totpos),] dat_prg_pos_sex <- subset(dat_prg_pos, sex != 'both') yr_sex <- unique(dat_prg_pos_sex$year) dat_prg_pos1 <- subset(dat_prg_pos, sex == 'both' & !(year %in% yr_sex)) dat_prg_pos2 <- subset(dat_prg_pos, sex != 'both' & (year %in% yr_sex)) dat_prg_pos <- rbind(dat_prg_pos1, dat_prg_pos2) - + dat_prg_pos$tot <- dat_prg_pos$totpos - + dat_prg_sex <- subset(dat_prg, sex != 'both') yr_sex <- unique(dat_prg_sex$year) dat_prg1 <- subset(dat_prg, sex == 'both' & !(year %in% yr_sex)) dat_prg2 <- subset(dat_prg, sex != 'both' & (year %in% yr_sex)) dat_prg <- rbind(dat_prg1, dat_prg2) - dat_prg <- dat_prg[complete.cases(dat_prg$tot),] - - # Programmatic data - nb of test (total: VCT + PMTCT + other; ideally among 15+) + dat_prg <- dat_prg[stats::complete.cases(dat_prg$tot),] + + # Programmatic data - nb of test (total: VCT + PMTCT + other; ideally among 15+) if (dim(dat_prg)[1] > 0) { dat_prg$l_est_se <- ifelse(dat_prg$sex == 'both', dat_prg$tot * 0.05, dat_prg$tot * 0.05 * sqrt(1 + 1 + 2*1)) # the sqrt() is to maintain same SE as if fitting on both sex - assuming perfect correlation dat_prg$hivstatus <- 'all' @@ -130,7 +131,7 @@ prepare_hts_likdat <- function(dat_evertest, dat_prg, fp) { dat_prg_pos$hivstatus <- 'positive' dat_prg_pos <- add_ss_indices(dat_prg_pos, fp$ss) } else { dat_prg_pos <- NULL } - + return(list(evertest = dat_evertest, hts = dat_prg, hts_pos = dat_prg_pos)) } @@ -143,47 +144,46 @@ lprior_hts <- function(theta, mod, fp) { ## * Extend knots by 1 year to current year knots <- 2000:2023 - fp$ss$proj_start + 1L ## -- UPDATE ABOVE -- - + n_k1 <- length(knots) - n_k2 <- n_k1 * 2 - 10 + n_k2 <- n_k1 * 2 - 10 penalty_f <- log(fp$hts_rate[1,2,1, knots[-1]]) - log(fp$hts_rate[1,2,1, knots[-n_k1]]) penalty_rr_dxunt <- theta[c((n_k1 + 2):n_k2)] - theta[c((n_k1 + 1):(n_k2 - 1))] - penalty_rr_m <- log(plogis(theta[n_k2 + 2]) * 1.1) - log(plogis(theta[n_k2 + 1]) * 1.1) + penalty_rr_m <- log(stats::plogis(theta[n_k2 + 2]) * 1.1) - log(stats::plogis(theta[n_k2 + 1]) * 1.1) penalty_rr_test <- theta[n_k2 + 4] - theta[n_k2 + 3] - lprior <- + lprior <- ## Prior for first baseline rate for females # exp(log(0.001) + 1.96*0.25) - dnorm(x = theta[1], mean = log(0.005), sd = 0.25, log = TRUE) + - ## Relative testing among PLHIV diagnosed, untreated. 1.50 (95%CI: 0.14-6.00) # plogis(qlogis(1.5/8) - 1.96*1.31)*8 - dnorm(x = theta[n_k2], mean = qlogis(1.5/8), sd = 1.31, log = TRUE) + - ## Prior for male RR 0.6 (95%CI: 0.07-1.05)# plogis(qlogis(0.6/1.1) + 1.96*1.46) * 1.1 - sum(dnorm(x = theta[n_k2 + 1], mean = qlogis(0.6/1.1), sd = 1.46, log = TRUE)) + - ## Relative increase among previously tested. 1.93 (95%CI: 1.08-5.00) # 0.95 + plogis(qlogis((1.93 - 0.95)/7.05) - qnorm(0.975)*1.084)*7.05 - dnorm(x = theta[n_k2 + 3], mean = qlogis((1.93 - 0.95)/7.05), sd = 1.084, log = TRUE) + - ## Relative factor for re-testing among PLHIV. 1.00 (95%CI: 0.10-1.90) # 0.05 + plogis(qlogis(0.95 / 1.90) - 1.96*1.85) * (1.95 - 0.05) - dnorm(x = theta[n_k2 + 5], mean = qlogis(0.95 / 1.90), sd = 1.85, log = TRUE) + - ## Relative testing among PLHIV already on ART (95%CI: 0.01-0.90) # plogis(qlogis(0.25) - 1.96*1.68) - dnorm(x = theta[n_k2 + 6], mean = qlogis(0.25), sd = 1.68, log = TRUE) + + stats::dnorm(x = theta[1], mean = log(0.005), sd = 0.25, log = TRUE) + + ## Relative testing among PLHIV diagnosed, untreated. 1.50 (95%CI: 0.14-6.00) # stats::plogis(stats::qlogis(1.5/8) - 1.96*1.31)*8 + stats::dnorm(x = theta[n_k2], mean = stats::qlogis(1.5/8), sd = 1.31, log = TRUE) + + ## Prior for male RR 0.6 (95%CI: 0.07-1.05)# stats::plogis(stats::qlogis(0.6/1.1) + 1.96*1.46) * 1.1 + sum(stats::dnorm(x = theta[n_k2 + 1], mean = stats::qlogis(0.6/1.1), sd = 1.46, log = TRUE)) + + ## Relative increase among previously tested. 1.93 (95%CI: 1.08-5.00) # 0.95 + stats::plogis(stats::qlogis((1.93 - 0.95)/7.05) - stats::qnorm(0.975)*1.084)*7.05 + stats::dnorm(x = theta[n_k2 + 3], mean = stats::qlogis((1.93 - 0.95)/7.05), sd = 1.084, log = TRUE) + + ## Relative factor for re-testing among PLHIV. 1.00 (95%CI: 0.10-1.90) # 0.05 + stats::plogis(stats::qlogis(0.95 / 1.90) - 1.96*1.85) * (1.95 - 0.05) + stats::dnorm(x = theta[n_k2 + 5], mean = stats::qlogis(0.95 / 1.90), sd = 1.85, log = TRUE) + + ## Relative testing among PLHIV already on ART (95%CI: 0.01-0.90) # stats::plogis(stats::qlogis(0.25) - 1.96*1.68) + stats::dnorm(x = theta[n_k2 + 6], mean = stats::qlogis(0.25), sd = 1.68, log = TRUE) + ## Prior for age (95% CI is 0.14-5.0)# 0.1 + invlogit(logit(0.9/5.9) - 1.96*1.5) * (6 - 0.1) - sum(dnorm(x = theta[c((n_k2 + 7):(n_k2 + 10))], mean = qlogis(0.9/5.9), sd = 1.685, log = TRUE)) + - ## RR OI diangosed for HIV relative to ART coverage 1.0 (0.3-1.7) # 0.25 + plogis(qlogis(0.5) - 1.96*1.75) * (1.75 - 0.25) - dnorm(x = theta[n_k2 + 11], mean = qlogis(0.5), sd = 1.75, log = TRUE) - + sum(stats::dnorm(x = theta[c((n_k2 + 7):(n_k2 + 10))], mean = stats::qlogis(0.9/5.9), sd = 1.685, log = TRUE)) + + ## RR OI diangosed for HIV relative to ART coverage 1.0 (0.3-1.7) # 0.25 + stats::plogis(stats::qlogis(0.5) - 1.96*1.75) * (1.75 - 0.25) + stats::dnorm(x = theta[n_k2 + 11], mean = stats::qlogis(0.5), sd = 1.75, log = TRUE) + prior_sd_f <- 0.205 # exp(log(0.5) - 1.96*0.205); previous of 0.35 when double-counting prior_sd_rr_m <- 0.26 # exp(log(0.6) + 1.96*0.26) prior_sd_rr_dxunt <- 0.25 prior_sd_rr_test <- 0.25 - - hyperprior <- + + hyperprior <- ## Penalty for 1st degree difference (female baseline rate) - sum(dnorm(x = penalty_f, mean = 0, sd = prior_sd_f, log = TRUE)) + - ## Penalty for 1st degree difference (male RR) - sum(dnorm(x = penalty_rr_m, mean = 0, sd = prior_sd_rr_m, log = TRUE)) + - ## Penalty for 1st degree difference (RR_dxunt) - sum(dnorm(x = penalty_rr_dxunt, mean = 0, sd = prior_sd_rr_dxunt, log = TRUE)) + - ## Penalty for 1st degree difference (RR_test) - sum(dnorm(x = penalty_rr_test, mean = 0, sd = prior_sd_rr_test, log = TRUE)) + sum(stats::dnorm(x = penalty_f, mean = 0, sd = prior_sd_f, log = TRUE)) + + ## Penalty for 1st degree difference (male RR) + sum(stats::dnorm(x = penalty_rr_m, mean = 0, sd = prior_sd_rr_m, log = TRUE)) + + ## Penalty for 1st degree difference (RR_dxunt) + sum(stats::dnorm(x = penalty_rr_dxunt, mean = 0, sd = prior_sd_rr_dxunt, log = TRUE)) + + ## Penalty for 1st degree difference (RR_test) + sum(stats::dnorm(x = penalty_rr_test, mean = 0, sd = prior_sd_rr_test, log = TRUE)) return(lprior + hyperprior) } - diff --git a/R/plot_functions.R b/R/plot_functions.R index e21656e..015ea38 100644 --- a/R/plot_functions.R +++ b/R/plot_functions.R @@ -20,61 +20,61 @@ get_pjnz_summary_data <- function(fp) { ## -- UPDATE HERE -- ## * update yr_pred to current year plot_pjnz_prv <- function(pjnz_summary, yr_pred = 2022) { - pjnz_summary <- na.omit(data.frame(year = pjnz_summary[["year"]], prv = pjnz_summary[["prevalence"]]*100)) + pjnz_summary <- stats::na.omit(data.frame(year = pjnz_summary[["year"]], prv = pjnz_summary[["prevalence"]]*100)) pjnz_summary$year <- pjnz_summary$year + 0.5 - plot(pjnz_summary$prv ~ pjnz_summary$year, + plot(pjnz_summary$prv ~ pjnz_summary$year, type = 'l', lwd = 2, xlim = c(2000, yr_pred), ylim = c(0, max(pjnz_summary$prv) * 1.2), main = 'HIV Prevalence (15-49 years)', ylab = 'HIV Prevalence (%)', xlab = 'Year', col='grey40') - polygon(x = c(pjnz_summary$year, rev(pjnz_summary$year)), + graphics::polygon(x = c(pjnz_summary$year, rev(pjnz_summary$year)), y = c(rep(0, length(pjnz_summary$prv)), rev(pjnz_summary$prv)), - col = rgb(155, 0, 50, 150, max = 255), border = NA) + col = grDevices::rgb(155, 0, 50, 150, max = 255), border = NA) } #' @export ## -- UPDATE HERE -- ## * update yr_pred to current year plot_pjnz_inc <- function(pjnz_summary, yr_pred = 2022) { - pjnz_summary <- na.omit(data.frame(year = pjnz_summary[["year"]], inc = pjnz_summary[["incidence"]]*1000)) + pjnz_summary <- stats::na.omit(data.frame(year = pjnz_summary[["year"]], inc = pjnz_summary[["incidence"]]*1000)) pjnz_summary$year <- pjnz_summary$year + 0.5 pjnz_summary <- subset(pjnz_summary, year >= 2000) - plot(pjnz_summary$inc ~ pjnz_summary$year, + plot(pjnz_summary$inc ~ pjnz_summary$year, type = 'l', lwd = 2, xlim = c(2000, yr_pred), ylim = c(0, max(pjnz_summary$inc) * 1.1), main='HIV Incidence (15-49 years)', ylab = 'HIV Incidence (per 1000)', xlab='Year', col='grey40') - polygon(x = c(pjnz_summary$year, rev(pjnz_summary$year)), + graphics::polygon(x = c(pjnz_summary$year, rev(pjnz_summary$year)), y = c(rep(0, length(pjnz_summary$inc)), rev(pjnz_summary$inc)), - col = rgb(255, 155, 0, 200, max = 255), border = NA) + col = grDevices::rgb(255, 155, 0, 200, max = 255), border = NA) } #' @export ## -- UPDATE HERE -- ## * update yr_pred to current year plot_pjnz_pop <- function(pjnz_summary, yr_pred = 2022) { - pjnz_summary <- na.omit(data.frame(year = pjnz_summary[["year"]], pop = pjnz_summary[["pop"]]/1000)) + pjnz_summary <- stats::na.omit(data.frame(year = pjnz_summary[["year"]], pop = pjnz_summary[["pop"]]/1000)) pjnz_summary$year <- pjnz_summary$year + 0.5 - plot(pjnz_summary$pop ~ pjnz_summary$year, + plot(pjnz_summary$pop ~ pjnz_summary$year, type = 'l', lwd = 2, xlim = c(2000, yr_pred), ylim = c(0, max(pjnz_summary$pop) * 1.2), main = 'Total Population (15-49 years)', ylab = 'Population (in 1,000s)', xlab='Year', col='grey40') - polygon(x = c(pjnz_summary$year, rev(pjnz_summary$year)), + graphics::polygon(x = c(pjnz_summary$year, rev(pjnz_summary$year)), y = c(rep(0, length(pjnz_summary$pop)), rev(pjnz_summary$pop)), - col = rgb(90, 170, 240, 150, max = 255), border = NA) + col = grDevices::rgb(90, 170, 240, 150, max = 255), border = NA) } #' @export ## -- UPDATE HERE -- ## * update yr_pred to current year plot_pjnz_plhiv <- function(pjnz_summary, yr_pred = 2022) { - pjnz_summary <- na.omit(data.frame(year = pjnz_summary[["year"]], plhiv = pjnz_summary[["plhiv"]]/1000)) + pjnz_summary <- stats::na.omit(data.frame(year = pjnz_summary[["year"]], plhiv = pjnz_summary[["plhiv"]]/1000)) pjnz_summary$year <- pjnz_summary$year + 0.5 - plot(pjnz_summary$plhiv ~ pjnz_summary$year, + plot(pjnz_summary$plhiv ~ pjnz_summary$year, type = 'l', lwd = 2, xlim = c(2000, yr_pred), ylim = c(0, max(pjnz_summary$plhiv) * 1.2), main = 'Total Number of PLHIV (15-49 years)', ylab = 'Population of PLHIV (in 1,000s)', xlab='Year', col='grey40') - polygon(x = c(pjnz_summary$year, rev(pjnz_summary$year)), + graphics::polygon(x = c(pjnz_summary$year, rev(pjnz_summary$year)), y = c(rep(0, length(pjnz_summary$plhiv)), rev(pjnz_summary$plhiv)), - col = rgb(250, 50, 10, 150, max = 255), border = NA) + col = grDevices::rgb(250, 50, 10, 150, max = 255), border = NA) } @@ -85,7 +85,7 @@ plot_pjnz_plhiv <- function(pjnz_summary, yr_pred = 2022) { ## * update yr_pred to current year plot_pjnz <- function(fp, yr_pred = 2022) { summary <- get_pjnz_summary_data(fp) - par(mfrow=c(2,2)) + graphics::par(mfrow=c(2,2)) plot_pjnz_prv(summary, yr_pred) plot_pjnz_inc(summary, yr_pred) plot_pjnz_pop(summary, yr_pred) @@ -142,15 +142,15 @@ plot_input_tot <- function(prgdat, fp, yr_pred = 2022) { xlim = c(2005, yr_pred), lwd = 2, col = 'steelblue3', xlab = 'Year', ylab = 'Total number of tests (in 1,000s)', pch = 16, main = expression(bold(paste("Total ", N^o, " of Tests Performed (VCT & ANC)")))) - lines(I(prg_t$tot/1000) ~ prg_t$year, lwd = 2, col = 'steelblue4') + graphics::lines(I(prg_t$tot/1000) ~ prg_t$year, lwd = 2, col = 'steelblue4') p_test_pop <- round(tot_prg$tot/pop[(tot_prg$year - start)]*100, 0) p_test_pop_r <- c(min(p_test_pop, na.rm = TRUE), max(p_test_pop, na.rm = TRUE)) - text(paste('No Test/Pop = ', p_test_pop_r[1], '-', p_test_pop_r[2], '%'), + graphics::text(paste('No Test/Pop = ', p_test_pop_r[1], '-', p_test_pop_r[2], '%'), x = 2005, y = max(tot_prg$tot/1000, na.rm = TRUE)/20, pos = 4) - lines(I(prg_f$tot/1000) ~ prg_f$year, lty = 2, lwd = 2, col = 'steelblue4') - points(I(prg_f$tot/1000) ~ prg_f$year, pch = 'f') - lines(I(prg_m$tot/1000) ~ prg_m$year, lty = 2, lwd = 2, col = 'steelblue4') - points(I(prg_m$tot/1000) ~ prg_m$year, pch = 'm') + graphics::lines(I(prg_f$tot/1000) ~ prg_f$year, lty = 2, lwd = 2, col = 'steelblue4') + graphics::points(I(prg_f$tot/1000) ~ prg_f$year, pch = 'f') + graphics::lines(I(prg_m$tot/1000) ~ prg_m$year, lty = 2, lwd = 2, col = 'steelblue4') + graphics::points(I(prg_m$tot/1000) ~ prg_m$year, pch = 'm') } } @@ -179,21 +179,21 @@ plot_input_totpos <- function(prgdat, fp, yr_pred = 2022) { ylab = 'Total of HIV+ tests (in 1,000s)', main = expression(bold(paste("Total ", N^o, " of HIV+ Tests (VCT & ANC)")))) - lines(I(prg_t$totpos/1000) ~ prg_t$year, lwd = 2, col = 'violetred4') + graphics::lines(I(prg_t$totpos/1000) ~ prg_t$year, lwd = 2, col = 'violetred4') yield <- round(tot_prg$totpos/tot_prg$tot * 100, 1) y_range <- c(min(yield, na.rm = TRUE), max(yield, na.rm = TRUE)) p_pos_pop <- round(tot_prg$totpos/plhiv[(tot_prg$year - start)]*100, 0) p_pos_pop_r <- c(min(p_pos_pop, na.rm = TRUE), max(p_pos_pop, na.rm = TRUE)) - text(paste('Positivity = ', y_range[1], '-', y_range[2], '%', sep = ''), + graphics::text(paste('Positivity = ', y_range[1], '-', y_range[2], '%', sep = ''), x = 2005, y = max(tot_prg$totpos/1000, na.rm = TRUE)/6, pos = 4) - text(paste('No Positive/PLHIV = ', p_pos_pop_r[1], '-', + graphics::text(paste('No Positive/PLHIV = ', p_pos_pop_r[1], '-', p_pos_pop_r[2], '%', sep=''), x = 2005, y = max(tot_prg$totpos/1000, na.rm = TRUE)/20, pos = 4) - lines(I(prg_f$totpos/1000) ~ prg_f$year, lty = 2, lwd = 2, col = 'violetred4') - points(I(prg_f$totpos/1000) ~ prg_f$year, pch = 'f') - lines(I(prg_m$totpos/1000) ~ prg_m$year, lty = 2, lwd = 2, col = 'violetred4') - points(I(prg_m$totpos/1000) ~ prg_m$year, pch = 'm') + graphics::lines(I(prg_f$totpos/1000) ~ prg_f$year, lty = 2, lwd = 2, col = 'violetred4') + graphics::points(I(prg_f$totpos/1000) ~ prg_f$year, pch = 'f') + graphics::lines(I(prg_m$totpos/1000) ~ prg_m$year, lty = 2, lwd = 2, col = 'violetred4') + graphics::points(I(prg_m$totpos/1000) ~ prg_m$year, pch = 'm') } } } @@ -210,7 +210,7 @@ plot_input_anctot <- function(prgdat, fp, yr_pred = 2022) { lwd = 2, col = 'darkorange2', xlab = 'Year', ylab = 'Total number of ANC tests (in 1,000s)', pch=16, main = expression(bold(paste("Total ", N^o, " of HIV Tests Performed at ANC")))) - lines(I(prgdat$anc/1000) ~ prgdat$year, lwd = 2, col = 'darkorange3') + graphics::lines(I(prgdat$anc/1000) ~ prgdat$year, lwd = 2, col = 'darkorange3') } } @@ -221,12 +221,12 @@ plot_input_ancpos <- function(prgdat, fp, yr_pred = 2022) { if (sum(prgdat$ancpos, na.rm = TRUE) > 0) { prgdat <- subset(prgdat, sex != 'male') - plot(I(prgdat$ancpos/1000) ~ prgdat$year, + plot(I(prgdat$ancpos/1000) ~ prgdat$year, ylim = c(0, max(prgdat$ancpos/1000, na.rm = TRUE)*1.1), xlim = c(2005, yr_pred ), - lwd = 2, col = 'deeppink2', xlab='Year', + lwd = 2, col = 'deeppink2', xlab='Year', ylab = 'Total number of positive ANC tests (in 1,000s)', pch = 17, main = expression(bold(paste("Total ", N^o, " of HIV+ Tests at ANC")))) - lines(I(prgdat$ancpos/1000) ~ prgdat$year, lwd = 2, col = 'deeppink2') + graphics::lines(I(prgdat$ancpos/1000) ~ prgdat$year, lwd = 2, col = 'deeppink2') } } @@ -236,42 +236,42 @@ plot_input_ancpos <- function(prgdat, fp, yr_pred = 2022) { ## -- UPDATE HERE -- ## * update yr_pred to current year plot_inputdata <- function(prgm_dat, fp, yr_pred = 2022) { - par(mfrow = c(2,2)) + graphics::par(mfrow = c(2,2)) plot_input_tot(prgm_dat, fp, yr_pred) plot_input_totpos(prgm_dat, fp, yr_pred) plot_input_anctot(prgm_dat, fp, yr_pred) - plot_input_ancpos(prgm_dat, fp, yr_pred) + plot_input_ancpos(prgm_dat, fp, yr_pred) } plot_prior <- function(plotprior = TRUE, sd = c(1.625, 1.85, 1.31, 1.68)){ if (plotprior == TRUE){ x <- logit(seq(0.00001, 1, 0.001)) - par(mfrow = c(1,4)) - rr_test <- dnorm(x, mean = logit(0.7/7.2), sd = sd[1]) - rr_plhiv <- dnorm(x, mean = logit(0.5), sd = sd[2]) - rr_dxunt <- dnorm(x, mean = logit(1.5/8), sd = sd[3]) - rr_dxart <- dnorm(x, mean = logit(0.25), sd = sd[4]) - plot(rr_test ~ I(0.8 + plogis(x) * 7.2), type = 'l', xlab = 'RR re-testing', ylab = 'Density', xlim = c(0,8)) - abline(v = 0.8 + 0.7/7.2 * 7.2, lty = 2, col = 'grey50') - lci <- round(0.8 + plogis(logit(0.7/7.2) - qnorm(0.975) * sd[1]) * 7.2, 1) - uci <- round(0.8 + plogis(logit(0.7/7.2) + qnorm(0.975) * sd[1]) * 7.2, 1) - text(x = 8, y = max(rr_test), paste('Prior: ', 1.5, '(95%: ', lci, '-', uci, ')'), pos = 2) - plot(rr_plhiv ~ I(0.05 + plogis(x) * 1.9), type = 'l', xlab = 'RR PLHIV', ylab = 'Density', xlim = c(0,2)) - abline(v = 0.05 + 0.5 * 1.9, lty = 2, col = 'grey50') - lci <- round(0.05 + plogis(logit(0.5) - qnorm(0.975) * sd[2]) * 1.9, 1) - uci <- round(0.05 + plogis(logit(0.5) + qnorm(0.975) * sd[2]) * 1.9, 1) - text(x = 2, y = max(rr_plhiv), paste('Prior: ', 1.0, '(95%: ', lci, '-', uci, ')'), pos = 2) - plot(rr_dxunt ~ I(plogis(x) * 8), type = 'l', xlab = 'RR re-testing for untreated diagnosed', ylab = 'Density', xlim = c(0,8)) - abline(v = 1.5/8 * 8, lty = 2, col = 'grey50') - lci <- round(plogis(logit(1.5/8) - qnorm(0.975) * sd[3]) * 8, 1) - uci <- round(plogis(logit(1.5/8) + qnorm(0.975) * sd[3]) * 8, 1) - text(x = 8, y = max(rr_test), paste('Prior: ', 1, '(95%: ', lci, '-', uci, ')'), pos = 2) - plot(rr_dxart ~ I(plogis(x)), type = 'l', xlab = 'RR re-testing on ART', ylab = 'Density', xlim = c(0,1)) - abline(v = 0.25, lty = 2, col = 'grey50') - lci <- round(plogis(logit(0.25) - qnorm(0.975) * sd[4]), 1) - uci <- round(plogis(logit(0.25) + qnorm(0.975) * sd[4]), 1) - text(x = 8, y = max(rr_test), paste('Prior: ', 0.25, '(95%: ', lci, '-', uci, ')'), pos = 2) + graphics::par(mfrow = c(1,4)) + rr_test <- stats::dnorm(x, mean = logit(0.7/7.2), sd = sd[1]) + rr_plhiv <- stats::dnorm(x, mean = logit(0.5), sd = sd[2]) + rr_dxunt <- stats::dnorm(x, mean = logit(1.5/8), sd = sd[3]) + rr_dxart <- stats::dnorm(x, mean = logit(0.25), sd = sd[4]) + plot(rr_test ~ I(0.8 + stats::plogis(x) * 7.2), type = 'l', xlab = 'RR re-testing', ylab = 'Density', xlim = c(0,8)) + graphics::abline(v = 0.8 + 0.7/7.2 * 7.2, lty = 2, col = 'grey50') + lci <- round(0.8 + stats::plogis(logit(0.7/7.2) - stats::qnorm(0.975) * sd[1]) * 7.2, 1) + uci <- round(0.8 + stats::plogis(logit(0.7/7.2) + stats::qnorm(0.975) * sd[1]) * 7.2, 1) + graphics::text(x = 8, y = max(rr_test), paste('Prior: ', 1.5, '(95%: ', lci, '-', uci, ')'), pos = 2) + plot(rr_plhiv ~ I(0.05 + stats::plogis(x) * 1.9), type = 'l', xlab = 'RR PLHIV', ylab = 'Density', xlim = c(0,2)) + graphics::abline(v = 0.05 + 0.5 * 1.9, lty = 2, col = 'grey50') + lci <- round(0.05 + stats::plogis(logit(0.5) - stats::qnorm(0.975) * sd[2]) * 1.9, 1) + uci <- round(0.05 + stats::plogis(logit(0.5) + stats::qnorm(0.975) * sd[2]) * 1.9, 1) + graphics::text(x = 2, y = max(rr_plhiv), paste('Prior: ', 1.0, '(95%: ', lci, '-', uci, ')'), pos = 2) + plot(rr_dxunt ~ I(stats::plogis(x) * 8), type = 'l', xlab = 'RR re-testing for untreated diagnosed', ylab = 'Density', xlim = c(0,8)) + graphics::abline(v = 1.5/8 * 8, lty = 2, col = 'grey50') + lci <- round(stats::plogis(logit(1.5/8) - stats::qnorm(0.975) * sd[3]) * 8, 1) + uci <- round(stats::plogis(logit(1.5/8) + stats::qnorm(0.975) * sd[3]) * 8, 1) + graphics::text(x = 8, y = max(rr_test), paste('Prior: ', 1, '(95%: ', lci, '-', uci, ')'), pos = 2) + plot(rr_dxart ~ I(stats::plogis(x)), type = 'l', xlab = 'RR re-testing on ART', ylab = 'Density', xlim = c(0,1)) + graphics::abline(v = 0.25, lty = 2, col = 'grey50') + lci <- round(stats::plogis(logit(0.25) - stats::qnorm(0.975) * sd[4]), 1) + uci <- round(stats::plogis(logit(0.25) + stats::qnorm(0.975) * sd[4]), 1) + graphics::text(x = 8, y = max(rr_test), paste('Prior: ', 0.25, '(95%: ', lci, '-', uci, ')'), pos = 2) } } @@ -279,19 +279,19 @@ plot_oi <- function(plotprior = TRUE, par = c(logit(0.5), 0.56)) { if (plotprior == TRUE){ r <- 0.5 tdx <- seq(1995, 2017, by = 0.1) - date_oidx <- 2005 + plogis(rnorm(2500, mean = par[1], sd = par[2])) * (2015 - 2005) - plot(I(rnorm(2)) ~ 1, pch = '', xlab = 'Year', ylab = '% OI Diagnosed', + date_oidx <- 2005 + stats::plogis(stats::rnorm(2500, mean = par[1], sd = par[2])) * (2015 - 2005) + plot(I(stats::rnorm(2)) ~ 1, pch = '', xlab = 'Year', ylab = '% OI Diagnosed', ylim = c(0, 1), xlim = c(1995, 2017)) for (i in 1:length(date_oidx)) { pr_oi <- 0.95 / (1 + exp(-r * (tdx - date_oidx[i]))) - lines(pr_oi ~ tdx, col = rgb(150, 100, 200, 25, max = 255), lwd = 1) + graphics::lines(pr_oi ~ tdx, col = grDevices::rgb(150, 100, 200, 25, max = 255), lwd = 1) } pr_oimin <- 0.95 / (1 + exp(-r * (tdx - 2005))) pr_oimax <- 0.95 / (1 + exp(-r * (tdx - 2015))) - pr_median <- 0.95 / (1 + exp(-r * (tdx - median(date_oidx)))) - lines(pr_oimin ~ tdx, col = rgb(100, 100, 100, 150, max = 255), lwd = 1, lty = 2) - lines(pr_oimax ~ tdx, col = rgb(100, 100, 100, 150, max = 255), lwd = 1, lty = 2) - lines(pr_median ~ tdx, col = 'steelblue4', lwd = 2, lty = 1) + pr_median <- 0.95 / (1 + exp(-r * (tdx - stats::median(date_oidx)))) + graphics::lines(pr_oimin ~ tdx, col = grDevices::rgb(100, 100, 100, 150, max = 255), lwd = 1, lty = 2) + graphics::lines(pr_oimax ~ tdx, col = grDevices::rgb(100, 100, 100, 150, max = 255), lwd = 1, lty = 2) + graphics::lines(pr_median ~ tdx, col = 'steelblue4', lwd = 2, lty = 1) } } @@ -299,16 +299,15 @@ plot_oi <- function(plotprior = TRUE, par = c(logit(0.5), 0.56)) { optimized_par <- function(opt, param = NULL) { n_k2 <- length(opt$par) - 11 n_k1 <- n_k2 - 10 - require(Matrix) - if (is.null(opt$hessian)) { - print('You forgot to specify Hessian = TRUE; 95%CrI not calculated.') - se <- rep(NA, length(opt$par)) + if (is.null(opt$hessian)) { + print('You forgot to specify Hessian = TRUE; 95%CrI not calculated.') + se <- rep(NA, length(opt$par)) } else { # We first test is system is singular (TEMPORARY FIX) if (any(diag(opt$hessian) == 0)) { index <- which(diag(opt$hessian) == 0) opt$hessian[index, index] <- opt$hessian[index, index] - 1e-3 - } + } # From the hessian, simulate the model vcova <- solve(-opt$hessian) # Test if positive semi-definite @@ -317,95 +316,90 @@ optimized_par <- function(opt, param = NULL) { if (!all(ev >= -1e-06 * abs(ev[1L]))) { vcova <- nearPD(vcova, corr = FALSE)$mat } se <- sqrt(diag(as.matrix(vcova))) } - + param_names <- c("RRm_05", "RRm_12", "RR_Test10","RR_Test15", - "RR_PLHIV", "RR_DxUnt10", "RR_DxUnt17", + "RR_PLHIV", "RR_DxUnt10", "RR_DxUnt17", "RR_DxART", "RR_25p_m", "RR_35p_m", "RR_25p_f", "RR_35p_f", "RR OI Dx") - - description <- c('RR testing: men in 2005', + + description <- c('RR testing: men in 2005', 'RR testing: men in 2012', - 'RR re-testing 2010', - 'RR re-testing 2015', + 'RR re-testing 2010', + 'RR re-testing 2015', 'RR testing: PLHIV unaware', 'RR re-testing: PLHIV aware (not ART) 2010', ## -- UPDATE HERE -- ## * update label to current year - 'RR re-testing: PLHIV aware (not ART) 2022', - 'RR re-testing: PLHIV on ART (*RR not ART)', + 'RR re-testing: PLHIV aware (not ART) 2022', + 'RR re-testing: PLHIV on ART (*RR not ART)', 'RR among 25-34 men', 'RR among 35+ men', 'RR among 25-34 women', - 'RR among 35+ women', + 'RR among 35+ women', 'RR OI Dx (ART Cov)') - - pt <- c(round(plogis(opt$par[n_k2 + 1]) * 1.1, 2), - round(plogis(opt$par[n_k2 + 2]) * 1.1, 2), - round(0.95 + plogis(opt$par[n_k2 + 3]) * 7.05, 2), - round(0.95 + plogis(opt$par[n_k2 + 4]) * 7.05, 2), - round(0.05 + plogis(opt$par[n_k2 + 5]) * (1.95 - 0.05), 2), - round(plogis(opt$par[n_k1 + 1]) * 8, 2), - round(plogis(opt$par[n_k2 - 1]) * 8, 2), - round(plogis(opt$par[n_k2 + 6]), 2), - round(0.1 + plogis(opt$par[n_k2 + 7]) * (6 - 0.1), 2), - round(0.1 + plogis(opt$par[n_k2 + 8]) * (6 - 0.1), 2), - round(0.1 + plogis(opt$par[n_k2 + 9]) * (6 - 0.1), 2), - round(0.1 + plogis(opt$par[n_k2 + 10]) * (6 - 0.1), 2), - round(0.25 + (plogis(opt$par[n_k2 + 11])) * (1.75 - 0.25), 1)) - + + pt <- c(round(stats::plogis(opt$par[n_k2 + 1]) * 1.1, 2), + round(stats::plogis(opt$par[n_k2 + 2]) * 1.1, 2), + round(0.95 + stats::plogis(opt$par[n_k2 + 3]) * 7.05, 2), + round(0.95 + stats::plogis(opt$par[n_k2 + 4]) * 7.05, 2), + round(0.05 + stats::plogis(opt$par[n_k2 + 5]) * (1.95 - 0.05), 2), + round(stats::plogis(opt$par[n_k1 + 1]) * 8, 2), + round(stats::plogis(opt$par[n_k2 - 1]) * 8, 2), + round(stats::plogis(opt$par[n_k2 + 6]), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 7]) * (6 - 0.1), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 8]) * (6 - 0.1), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 9]) * (6 - 0.1), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 10]) * (6 - 0.1), 2), + round(0.25 + (stats::plogis(opt$par[n_k2 + 11])) * (1.75 - 0.25), 1)) + if (is.null(param)) { - lci <- c(round(plogis(opt$par[n_k2 + 1] - qnorm(0.975) * se[n_k2 + 1]) * 1.1, 2), - round(plogis(opt$par[n_k2 + 2] - qnorm(0.975) * se[n_k2 + 2]) * 1.1, 2), - round(0.95 + plogis(opt$par[n_k2 + 3] - qnorm(0.975) * se[n_k2 + 3]) * 7.05, 2), - round(0.95 + plogis(opt$par[n_k2 + 4] - qnorm(0.975) * se[n_k2 + 4]) * 7.05, 2), - round(0.05 + plogis(opt$par[n_k2 + 5] - qnorm(0.975) * se[n_k2 + 5]) * (1.95 - 0.05), 2), - round(plogis(opt$par[n_k1 + 1] - qnorm(0.975) * se[n_k1 + 1]) * 8, 2), - round(plogis(opt$par[n_k2 - 1] - qnorm(0.975) * se[n_k2 - 1]) * 8, 2), - round(plogis(opt$par[n_k2 + 6] - qnorm(0.975) * se[n_k2 + 6]), 2), - round(0.1 + plogis(opt$par[n_k2 + 7] - qnorm(0.975) * se[n_k2 + 7]) * (6 - 0.1), 2), - round(0.1 + plogis(opt$par[n_k2 + 8] - qnorm(0.975) * se[n_k2 + 8]) * (6 - 0.1), 2), - round(0.1 + plogis(opt$par[n_k2 + 9] - qnorm(0.975) * se[n_k2 + 9]) * (6 - 0.1), 2), - round(0.1 + plogis(opt$par[n_k2 + 10] - qnorm(0.975) * se[n_k2 + 10]) * (6 - 0.1), 2), - round(0.25 + (plogis(opt$par[n_k2 + 11] - qnorm(0.975) * se[n_k2 + 11])) * (1.75 - 0.25), 1)) - - uci <- c(round(plogis(opt$par[n_k2 + 1] + qnorm(0.975) * se[n_k2 + 1]) * 1.1, 2), - round(plogis(opt$par[n_k2 + 2] + qnorm(0.975) * se[n_k2 + 2]) * 1.1, 2), - round(0.95 + plogis(opt$par[n_k2 + 3] + qnorm(0.975) * se[n_k2 + 3]) * 7.05, 2), - round(0.95 + plogis(opt$par[n_k2 + 4] + qnorm(0.975) * se[n_k2 + 4]) * 7.05, 2), - round(0.05 + plogis(opt$par[n_k2 + 5] + qnorm(0.975) * se[n_k2 + 5]) * (1.95 - 0.05), 2), - round(plogis(opt$par[n_k1 + 1] + qnorm(0.975) * se[n_k1 + 1]) * 8, 2), - round(plogis(opt$par[n_k2 - 1] + qnorm(0.975) * se[n_k2 - 1]) * 8, 2), - round(plogis(opt$par[n_k2 + 6] + qnorm(0.975) * se[n_k2 + 6]), 2), - round(0.1 + plogis(opt$par[n_k2 + 7] + qnorm(0.975) * se[n_k2 + 7]) * (6 - 0.1), 2), - round(0.1 + plogis(opt$par[n_k2 + 8] + qnorm(0.975) * se[n_k2 + 8]) * (6 - 0.1), 2), - round(0.1 + plogis(opt$par[n_k2 + 9] + qnorm(0.975) * se[n_k2 + 9]) * (6 - 0.1), 2), - round(0.1 + plogis(opt$par[n_k2 + 10] + qnorm(0.975) * se[n_k2 + 10]) * (6 - 0.1), 2), - round(0.25 + (plogis(opt$par[n_k2 + 11] + qnorm(0.975) * se[n_k2 + 11])) * (1.75 - 0.25), 1)) + lci <- c(round(stats::plogis(opt$par[n_k2 + 1] - stats::qnorm(0.975) * se[n_k2 + 1]) * 1.1, 2), + round(stats::plogis(opt$par[n_k2 + 2] - stats::qnorm(0.975) * se[n_k2 + 2]) * 1.1, 2), + round(0.95 + stats::plogis(opt$par[n_k2 + 3] - stats::qnorm(0.975) * se[n_k2 + 3]) * 7.05, 2), + round(0.95 + stats::plogis(opt$par[n_k2 + 4] - stats::qnorm(0.975) * se[n_k2 + 4]) * 7.05, 2), + round(0.05 + stats::plogis(opt$par[n_k2 + 5] - stats::qnorm(0.975) * se[n_k2 + 5]) * (1.95 - 0.05), 2), + round(stats::plogis(opt$par[n_k1 + 1] - stats::qnorm(0.975) * se[n_k1 + 1]) * 8, 2), + round(stats::plogis(opt$par[n_k2 - 1] - stats::qnorm(0.975) * se[n_k2 - 1]) * 8, 2), + round(stats::plogis(opt$par[n_k2 + 6] - stats::qnorm(0.975) * se[n_k2 + 6]), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 7] - stats::qnorm(0.975) * se[n_k2 + 7]) * (6 - 0.1), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 8] - stats::qnorm(0.975) * se[n_k2 + 8]) * (6 - 0.1), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 9] - stats::qnorm(0.975) * se[n_k2 + 9]) * (6 - 0.1), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 10] - stats::qnorm(0.975) * se[n_k2 + 10]) * (6 - 0.1), 2), + round(0.25 + (stats::plogis(opt$par[n_k2 + 11] - stats::qnorm(0.975) * se[n_k2 + 11])) * (1.75 - 0.25), 1)) + + uci <- c(round(stats::plogis(opt$par[n_k2 + 1] + stats::qnorm(0.975) * se[n_k2 + 1]) * 1.1, 2), + round(stats::plogis(opt$par[n_k2 + 2] + stats::qnorm(0.975) * se[n_k2 + 2]) * 1.1, 2), + round(0.95 + stats::plogis(opt$par[n_k2 + 3] + stats::qnorm(0.975) * se[n_k2 + 3]) * 7.05, 2), + round(0.95 + stats::plogis(opt$par[n_k2 + 4] + stats::qnorm(0.975) * se[n_k2 + 4]) * 7.05, 2), + round(0.05 + stats::plogis(opt$par[n_k2 + 5] + stats::qnorm(0.975) * se[n_k2 + 5]) * (1.95 - 0.05), 2), + round(stats::plogis(opt$par[n_k1 + 1] + stats::qnorm(0.975) * se[n_k1 + 1]) * 8, 2), + round(stats::plogis(opt$par[n_k2 - 1] + stats::qnorm(0.975) * se[n_k2 - 1]) * 8, 2), + round(stats::plogis(opt$par[n_k2 + 6] + stats::qnorm(0.975) * se[n_k2 + 6]), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 7] + stats::qnorm(0.975) * se[n_k2 + 7]) * (6 - 0.1), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 8] + stats::qnorm(0.975) * se[n_k2 + 8]) * (6 - 0.1), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 9] + stats::qnorm(0.975) * se[n_k2 + 9]) * (6 - 0.1), 2), + round(0.1 + stats::plogis(opt$par[n_k2 + 10] + stats::qnorm(0.975) * se[n_k2 + 10]) * (6 - 0.1), 2), + round(0.25 + (stats::plogis(opt$par[n_k2 + 11] + stats::qnorm(0.975) * se[n_k2 + 11])) * (1.75 - 0.25), 1)) } if (!is.null(param)){ - est <- rbind(plogis(param[, n_k2 + 1]) * 1.1, - plogis(param[, n_k2 + 2]) * 1.1, - 0.95 + plogis(param[, n_k2 + 3]) * 7.05, - 0.95 + plogis(param[, n_k2 + 4]) * 7.05, - 0.05 + plogis(param[, n_k2 + 5]) * (1.95 - 0.05), - plogis(param[, n_k1 + 1]) * 8, - plogis(param[, n_k2 - 1]) * 8, - plogis(param[, n_k2 + 6]), - 0.1 + plogis(param[, n_k2 + 7]) * (6 - 0.1), - 0.1 + plogis(param[, n_k2 + 8]) * (6 - 0.1), - 0.1 + plogis(param[, n_k2 + 9]) * (6 - 0.1), - 0.1 + plogis(param[, n_k2 + 10]) * (6 - 0.1), - 0.25 + (plogis(param[, n_k2 + 11])) * (1.75 - 0.25)) - lci <- round(apply(est, 1, quantile, probs = 0.025), 2) - uci <- round(apply(est, 1, quantile, probs = 0.975), 2) - + est <- rbind(stats::plogis(param[, n_k2 + 1]) * 1.1, + stats::plogis(param[, n_k2 + 2]) * 1.1, + 0.95 + stats::plogis(param[, n_k2 + 3]) * 7.05, + 0.95 + stats::plogis(param[, n_k2 + 4]) * 7.05, + 0.05 + stats::plogis(param[, n_k2 + 5]) * (1.95 - 0.05), + stats::plogis(param[, n_k1 + 1]) * 8, + stats::plogis(param[, n_k2 - 1]) * 8, + stats::plogis(param[, n_k2 + 6]), + 0.1 + stats::plogis(param[, n_k2 + 7]) * (6 - 0.1), + 0.1 + stats::plogis(param[, n_k2 + 8]) * (6 - 0.1), + 0.1 + stats::plogis(param[, n_k2 + 9]) * (6 - 0.1), + 0.1 + stats::plogis(param[, n_k2 + 10]) * (6 - 0.1), + 0.25 + (stats::plogis(param[, n_k2 + 11])) * (1.75 - 0.25)) + lci <- round(apply(est, 1, stats::quantile, probs = 0.025), 2) + uci <- round(apply(est, 1, stats::quantile, probs = 0.975), 2) + } RR_opt <- data.frame(Parameter_Name = description, Estimate = pt, LCI = lci, UCI = uci) - + return(RR_opt) } - - - - - diff --git a/R/plot_outputs.R b/R/plot_outputs.R index 43c1627..bc3331b 100644 --- a/R/plot_outputs.R +++ b/R/plot_outputs.R @@ -3,7 +3,7 @@ get_out_evertest <- function(mod, fp, agegr = c("15-24", '25-34','35-49', '15-49', '15-64', '15+'), sex = c("both", "female", "male"), hivstatus = c("all", "negative", "positive") ) { - + end_date <- fp$ss$proj_start + fp$ss$PROJ_YEARS - 1L out_evertest <- expand.grid(year = 2000:end_date, outcome = "evertest", @@ -12,13 +12,13 @@ get_out_evertest <- function(mod, fp, hivstatus = hivstatus) out_evertest$value <- evertest(mod, fp, add_ss_indices(out_evertest, fp$ss)) - + out_evertest } #' @export get_out_nbtest <- function(mod, fp) { - + end_date <- fp$ss$proj_start + fp$ss$PROJ_YEARS - 1L out_nbtest <- expand.grid(year = 2000:end_date, outcome = "numbertests", @@ -32,7 +32,7 @@ get_out_nbtest <- function(mod, fp) { #' @export get_out_nbtest_sex <- function(mod, fp) { - + end_date <- fp$ss$proj_start + fp$ss$PROJ_YEARS - 1L out_nbtest <- expand.grid(year = 2000:end_date, outcome = "numbertests", @@ -54,7 +54,7 @@ get_out_nbtest_pos <- function(mod, fp) { hivstatus = 'positive') out_nbtest_pos$hivstatus <- as.character(out_nbtest_pos$hivstatus) out_nbtest_pos$value <- total_tests(mod, add_ss_indices(out_nbtest_pos, fp$ss)) - + out_nbtest_pos } @@ -68,7 +68,7 @@ get_out_nbtest_pos_sex <- function(mod, fp) { hivstatus = 'positive') out_nbtest_pos$hivstatus <- as.character(out_nbtest_pos$hivstatus) out_nbtest_pos$value <- total_tests(mod, add_ss_indices(out_nbtest_pos, fp$ss)) - + out_nbtest_pos } @@ -82,7 +82,7 @@ get_out_aware <- function(mod, fp, agegr = c('15-49', '15-64', "15+"), sex = sex, hivstatus = 'positive') out_aware$value <- diagnosed(mod, fp, add_ss_indices(out_aware, fp$ss)) - + out_aware } @@ -96,7 +96,7 @@ get_out_nbaware <- function(mod, fp, agegr = '15-49', sex = sex, hivstatus = 'positive') out_nbaware$value <- nb_aware(mod, fp, add_ss_indices(out_nbaware, fp$ss)) - + out_nbaware } @@ -107,7 +107,7 @@ get_out_art <- function(mod, fp, out_artcov <- data.frame(year = 1970:end_date, outcome = "artcov", agegr = "15-49", sex = gender, hivstatus = "positive", value = artcov15to49(mod, sex = gender)) out_artcov <- subset(out_artcov) - + out_artcov } @@ -122,18 +122,18 @@ get_out_pregprev <- function(mod, fp) { undiagn_w <- colSums((attr(mod, "hivpop") - attr(mod, "diagnpop"))[, fp$ss$h.fert.idx, fp$ss$f.idx, ]) diagn_w <- colSums(attr(mod, "diagnpop")[, fp$ss$h.fert.idx, fp$ss$f.idx, ]) artpop_w <- colSums(attr(mod, "artpop")[, , fp$ss$h.fert.idx, fp$ss$f.idx, ],,2) - + tot_w <- (hivneg_w + undiagn_w + diagn_w + artpop_w) - + ## Number of births by HIV, awareness, and ART status hivneg_births <- colSums(births_ha * hivneg_w / tot_w) undiagn_births <- colSums(births_ha * undiagn_w / tot_w) diagn_births <- colSums(births_ha * diagn_w / tot_w) artpop_births <- colSums(births_ha * artpop_w / tot_w) - + births <- colSums(births_ha) hivp_births <- births - hivneg_births - + data.frame(year = fp$ss$proj_start - 1 + seq_len(fp$ss$PROJ_YEARS), prev = hivp_births / births, aware = (diagn_births + artpop_births) / hivp_births, @@ -143,7 +143,7 @@ get_out_pregprev <- function(mod, fp) { # ---- Individuals functions to plot model outputs ---- #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year plot_out_nbtest <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = 2022, @@ -151,30 +151,30 @@ plot_out_nbtest <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = 2022, # if fitting with HTS program data stratified by sex, we add both sex back ld <- likdat$hts if (!is.null(ld)) { - likdat$hts <- aggregate(ld$tot, by = list(ld$year), FUN = sum) - colnames(likdat$hts) <- c('year','tot') + likdat$hts <- stats::aggregate(ld$tot, by = list(ld$year), FUN = sum) + colnames(likdat$hts) <- c('year','tot') } # redact <- c('Namibia','Uganda','Zambia','Zimbabwe') redact <- c('XXX') - + cnt_to_plot <- ifelse(cnt == "Cote d'Ivoire", "Côte d'Ivoire", - ifelse(cnt == "Swaziland", "eSwatini", cnt)) + ifelse(cnt == "Swaziland", "eSwatini", cnt)) start <- fp$ss$proj_start mod <- simmod(fp) - plhiv <- apply(attr(mod, "hivpop")[,1:8,,], 4, FUN = sum) + + plhiv <- apply(attr(mod, "hivpop")[,1:8,,], 4, FUN = sum) + apply(attr(mod, "artpop")[,,1:8,,], 5, FUN = sum) pop <- apply(mod[1:35,,,], 4, FUN = sum) out_nbtest <- subset(get_out_nbtest(mod, fp), year <= yr_pred) out_nbtest$year <- out_nbtest$year + 0.5 - + if (!is.null(ld) & any(ld$sex != 'both')) { ld_sex <- subset(ld, sex != 'both') ld_sex_pred <- total_tests(mod, df = add_ss_indices(ld_sex, fp$ss)) } - + # Decide if we plot CI or not if (!is.null(simul)) { ci <- getCI(simul$number.test) @@ -186,76 +186,76 @@ plot_out_nbtest <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = 2022, } else { plot.ci <- FALSE } - + if (plot_title == TRUE) { main_title <- expression(bold(paste("Total ", N^o, " of Tests"))) } else { main_title <- "" } - - col_ci <- rgb(150, 150, 150, 125, max = 255) - + + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + # Total tested if (plot.ci == T){ plot(I(out_nbtest$value / 1000) ~ out_nbtest$year, pch = '', ylim = c(0, max(out_nbtest$value, na.rm = TRUE) / 1000 * 1.25), - cex = 1, ylab = expression(paste(N^o, " of Tests (in 1,000)")), - xlab = 'Year', xlim = c(2000, yr_pred), + cex = 1, ylab = expression(paste(N^o, " of Tests (in 1,000)")), + xlab = 'Year', xlim = c(2000, yr_pred), main = main_title) - polygon(x = c(ci$year, rev(ci$year)), + graphics::polygon(x = c(ci$year, rev(ci$year)), y = c(ci$upper / 1000, rev(ci$lower / 1000)), col = col_ci, border = NA) - lines(I(out_nbtest$value / 1000) ~ out_nbtest$year, lwd = 1, col = 'seagreen3') - text(x = 2000, y = max(out_nbtest$value, na.rm = TRUE) / 1000 * 1.15, cnt_to_plot, cex = 1.25, pos = 4) - + graphics::lines(I(out_nbtest$value / 1000) ~ out_nbtest$year, lwd = 1, col = 'seagreen3') + graphics::text(x = 2000, y = max(out_nbtest$value, na.rm = TRUE) / 1000 * 1.15, cnt_to_plot, cex = 1.25, pos = 4) + if(length(likdat$hts$tot) > 0) { if (cnt %in% redact) { pchpt <- '' } else { pchpt <- 16 } - points(I(likdat$hts$tot/1000) ~ I(likdat$hts$year + 0.5), pch = pchpt) + graphics::points(I(likdat$hts$tot/1000) ~ I(likdat$hts$year + 0.5), pch = pchpt) p_test_pop <- round(likdat$hts$tot / pop[likdat$hts$year - start] * 100, 0) p_test_pop_r <- c(min(p_test_pop, na.rm = TRUE), max(p_test_pop, na.rm = TRUE)) - text(paste('No Test/Pop = ', p_test_pop_r[1], '-', p_test_pop_r[2], '%'), + graphics::text(paste('No Test/Pop = ', p_test_pop_r[1], '-', p_test_pop_r[2], '%'), x = yr_pred, y = max(likdat$hts$tot / 1000, na.rm = TRUE)/25, pos = 2) } } else { plot(I(out_nbtest$value/1000) ~ out_nbtest$year, pch = '', ylim = c(0, max(out_nbtest$value, na.rm = TRUE) / 1000 + max(out_nbtest$value, na.rm = TRUE)/1000*0.25), - cex = 1, ylab = expression(paste(N^o, " of Tests (in 1,000)")), - xlab = 'Year', xlim = c(2000, yr_pred), + cex = 1, ylab = expression(paste(N^o, " of Tests (in 1,000)")), + xlab = 'Year', xlim = c(2000, yr_pred), main = main_title) - lines(I(out_nbtest$value/1000) ~ out_nbtest$year, lwd = 1, col = 'seagreen3') - text(x = 2000, y = max(out_nbtest$value, na.rm = TRUE) / 1000, cnt_to_plot, cex = 1.5, pos = 2) + graphics::lines(I(out_nbtest$value/1000) ~ out_nbtest$year, lwd = 1, col = 'seagreen3') + graphics::text(x = 2000, y = max(out_nbtest$value, na.rm = TRUE) / 1000, cnt_to_plot, cex = 1.5, pos = 2) if(length(likdat$hts$tot) > 0) { if (cnt %in% redact) { pchpt <- '' } else { pchpt <- 16 } - points(I(likdat$hts$tot / 1000) ~ I(likdat$hts$year + 0.5), pch = pchpt) + graphics::points(I(likdat$hts$tot / 1000) ~ I(likdat$hts$year + 0.5), pch = pchpt) p_test_pop <- round(likdat$hts$tot/pop[(likdat$hts$year - start)] * 100, 0) p_test_pop_r <- c(min(p_test_pop, na.rm = TRUE), max(p_test_pop, na.rm = TRUE)) - text(paste('No Test/Pop = ', p_test_pop_r[1], '-', p_test_pop_r[2], '%'), + graphics::text(paste('No Test/Pop = ', p_test_pop_r[1], '-', p_test_pop_r[2], '%'), x = yr_pred, y = max(likdat$hts$tot / 1000, na.rm = TRUE) / 20, pos = 2) } } } -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year plot_out_nbtest_sex <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = 2022) { - + # redact <- c('Namibia','Uganda','Zambia','Zimbabwe') redact <- c('XXX') - + cnt_to_plot <- ifelse(cnt == "Cote d'Ivoire", "Côte d'Ivoire", - ifelse(cnt == "Swaziland", "eSwatini", cnt)) + ifelse(cnt == "Swaziland", "eSwatini", cnt)) start <- fp$ss$proj_start mod <- simmod(fp) - plhiv <- apply(attr(mod, "hivpop")[,1:8,,], 4, FUN=sum) + + plhiv <- apply(attr(mod, "hivpop")[,1:8,,], 4, FUN=sum) + apply(attr(mod, "artpop")[,,1:8,,], 5, FUN=sum) pop <- apply(mod[1:35,,,], 4, FUN=sum) - + out_nbtest <- subset(get_out_nbtest_sex(mod, fp), year <= yr_pred) out_nbtest$year <- out_nbtest$year + 0.5 out_nbtest_f <- subset(out_nbtest, sex == 'female') out_nbtest_m <- subset(out_nbtest, sex == 'male') max_ylim <- max(c(out_nbtest_f$value, out_nbtest_m$value), na.rm = TRUE) - + lik_sex <- likdat$hts if (!is.null(lik_sex)){ lik_f <- subset(lik_sex, sex == 'female') @@ -274,51 +274,51 @@ plot_out_nbtest_sex <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = 20 } else { plot.ci <- FALSE } - - col_ci <- rgb(150, 150, 150, 125, max = 255) - + + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + # Total tested if (plot.ci == T){ plot(I(out_nbtest_f$value / 1000) ~ out_nbtest_f$year, pch = '', ylim = c(0, max_ylim / 1000 * 1.25), - cex = 1, ylab = expression(paste(N^o, " of Tests (in 1,000)")), - xlab = 'Year', xlim = c(2000, yr_pred), + cex = 1, ylab = expression(paste(N^o, " of Tests (in 1,000)")), + xlab = 'Year', xlim = c(2000, yr_pred), main = expression(bold(paste("Total ", N^o, " of Tests")))) - polygon(x = c(ci_f$year, rev(ci_f$year)), + graphics::polygon(x = c(ci_f$year, rev(ci_f$year)), y = c(ci_f$upper / 1000, rev(ci_f$lower / 1000)), col = col_ci, border = NA) - polygon(x = c(ci_m$year, rev(ci_m$year)), + graphics::polygon(x = c(ci_m$year, rev(ci_m$year)), y = c(ci_m$upper / 1000, rev(ci_m$lower / 1000)), col = col_ci, border = NA) - lines(I(out_nbtest_f$value / 1000) ~ out_nbtest_f$year, lwd = 1, col = 'seagreen3') - lines(I(out_nbtest_m$value / 1000) ~ out_nbtest_m$year, lwd = 1, col = 'seagreen3') - text(x = 2000, y = max(out_nbtest$value, na.rm = TRUE) / 1000 * 1.15, cnt_to_plot, cex = 1.25, pos = 4) - + graphics::lines(I(out_nbtest_f$value / 1000) ~ out_nbtest_f$year, lwd = 1, col = 'seagreen3') + graphics::lines(I(out_nbtest_m$value / 1000) ~ out_nbtest_m$year, lwd = 1, col = 'seagreen3') + graphics::text(x = 2000, y = max(out_nbtest$value, na.rm = TRUE) / 1000 * 1.15, cnt_to_plot, cex = 1.25, pos = 4) + if(length(lik_f$tot) > 0 | length(lik_m$tot) > 0) { if (cnt %in% redact) { pchpt <- '' } else { pchpt <- 16 } - points(I(lik_f$tot / 1000) ~ I(lik_f$year + 0.5), pch = pchpt) - points(I(lik_m$tot / 1000) ~ I(lik_m$year + 0.5), pch = pchpt + 1) + graphics::points(I(lik_f$tot / 1000) ~ I(lik_f$year + 0.5), pch = pchpt) + graphics::points(I(lik_m$tot / 1000) ~ I(lik_m$year + 0.5), pch = pchpt + 1) } } else { plot(I(out_nbtest_f$value / 1000) ~ out_nbtest_f$year, pch = '', ylim = c(0, max_ylim / 1000 * 1.25), - cex = 1, ylab = expression(paste(N^o, " of Tests (in 1,000)")), - xlab = 'Year', xlim = c(2000, yr_pred), + cex = 1, ylab = expression(paste(N^o, " of Tests (in 1,000)")), + xlab = 'Year', xlim = c(2000, yr_pred), main = expression(bold(paste("Total ", N^o, " of Tests")))) - lines(I(out_nbtest_f$value / 1000) ~ out_nbtest_f$year, lwd = 1, col = 'seagreen3') - lines(I(out_nbtest_m$value / 1000) ~ out_nbtest_m$year, lwd = 1, col = 'seagreen3') - text(x = 2000, y = max(out_nbtest$value, na.rm = TRUE) / 1000, cnt_to_plot, cex = 1.5, pos = 2) - + graphics::lines(I(out_nbtest_f$value / 1000) ~ out_nbtest_f$year, lwd = 1, col = 'seagreen3') + graphics::lines(I(out_nbtest_m$value / 1000) ~ out_nbtest_m$year, lwd = 1, col = 'seagreen3') + graphics::text(x = 2000, y = max(out_nbtest$value, na.rm = TRUE) / 1000, cnt_to_plot, cex = 1.5, pos = 2) + if(length(likdat$hts$tot) > 0) { if (cnt %in% redact) { pchpt <- '' } else { pchpt <- 16 } - points(I(lik_f$tot / 1000) ~ I(lik_f$year + 0.5), pch = pchpt) - points(I(lik_m$tot / 1000) ~ I(lik_m$year + 0.5), pch = pchpt + 1) } + graphics::points(I(lik_f$tot / 1000) ~ I(lik_f$year + 0.5), pch = pchpt) + graphics::points(I(lik_m$tot / 1000) ~ I(lik_m$year + 0.5), pch = pchpt + 1) } } } #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year plot_out_nbpostest <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = 2022, @@ -326,29 +326,29 @@ plot_out_nbpostest <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = 202 # if fitting with HTS program data stratified by sex, we add both sex back ld <- likdat$hts if (!is.null(ld)) { - likdat$hts <- aggregate(ld$tot, by = list(ld$year), FUN = sum) + likdat$hts <- stats::aggregate(ld$tot, by = list(ld$year), FUN = sum) colnames(likdat$hts) <- c('year','tot') - likdat$hts$totpos <- aggregate(ld$totpos, by = list(ld$year), FUN = sum)$x + likdat$hts$totpos <- stats::aggregate(ld$totpos, by = list(ld$year), FUN = sum)$x } ldp <- likdat$hts_pos if (!is.null(ldp)) { - likdat$hts_pos <- aggregate(ldp$tot, by = list(ldp$year), FUN = sum) + likdat$hts_pos <- stats::aggregate(ldp$tot, by = list(ldp$year), FUN = sum) colnames(likdat$hts_pos) <- c('year','tot') } - + # redact <- c('Namibia','Uganda','Zambia','Zimbabwe') redact <- c('XXX') - + start <- fp$ss$proj_start mod <- simmod(fp) - plhiv <- apply(attr(mod, "hivpop")[, 1:8,,], 4, FUN = sum) + + plhiv <- apply(attr(mod, "hivpop")[, 1:8,,], 4, FUN = sum) + apply(attr(mod, "artpop")[,, 1:8,,], 5, FUN = sum) pop <- apply(mod[1:35,,,], 4, FUN = sum) out_nbtest_pos <- get_out_nbtest_pos(mod, fp) out_nbtest_pos <- subset(out_nbtest_pos, year <= yr_pred) - out_nbtest_pos$year <- out_nbtest_pos$year + 0.5 - + out_nbtest_pos$year <- out_nbtest_pos$year + 0.5 + # Decide if we plot CI or not if (!is.null(simul)){ ci <- getCI(simul$number.test.pos) @@ -366,77 +366,77 @@ plot_out_nbpostest <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = 202 } else { main_title <- "" } - - col_ci <- rgb(150, 150, 150, 125, max = 255) - + + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + # Postive tests if (plot.ci == T){ plot(I(out_nbtest_pos$value/1000) ~ out_nbtest_pos$year, pch = '', ylim = c(0, max(out_nbtest_pos$value, na.rm = TRUE) / 1000 * 1.25), - cex = 1, ylab = expression(paste(N^o, " of Positive Tests (in 1,000)")), - xlab='Year', xlim = c(2000, yr_pred), + cex = 1, ylab = expression(paste(N^o, " of Positive Tests (in 1,000)")), + xlab='Year', xlim = c(2000, yr_pred), main = main_title) - polygon(x = c(ci$year, rev(ci$year)), + graphics::polygon(x = c(ci$year, rev(ci$year)), y = c(ci$upper / 1000, rev(ci$lower / 1000)), col = col_ci, border = NA) - lines(I(out_nbtest_pos$value / 1000) ~ out_nbtest_pos$year, lwd = 1, col = 'orangered2') + graphics::lines(I(out_nbtest_pos$value / 1000) ~ out_nbtest_pos$year, lwd = 1, col = 'orangered2') if (!is.null(likdat$hts_pos$tot)) { if (cnt %in% redact) { pchpt <- '' } else { pchpt <- 16 } - points(I(likdat$hts_pos$tot / 1000) ~ I(likdat$hts_pos$year + 0.5), pch = pchpt) + graphics::points(I(likdat$hts_pos$tot / 1000) ~ I(likdat$hts_pos$year + 0.5), pch = pchpt) yield <- round(likdat$hts$totpos / likdat$hts$tot * 100, 1) y_range <- c(min(yield, na.rm = TRUE), max(yield, na.rm = TRUE)) p_pos_pop <- round(likdat$hts_pos$tot / plhiv[(likdat$hts_pos$year - start)] * 100, 0) p_pos_pop_r <- c(min(p_pos_pop, na.rm = TRUE), max(p_pos_pop, na.rm = TRUE)) - text(paste('Positivity = ', y_range[1], '-', y_range[2], '%', sep=''), + graphics::text(paste('Positivity = ', y_range[1], '-', y_range[2], '%', sep=''), x = yr_pred, y = max(likdat$hts_pos$tot / 1000, na.rm = TRUE) / 6, pos = 2) - text(paste('No Positive/PLHIV = ', p_pos_pop_r[1], '-', - p_pos_pop_r[2], '%', sep=''), x = yr_pred, + graphics::text(paste('No Positive/PLHIV = ', p_pos_pop_r[1], '-', + p_pos_pop_r[2], '%', sep=''), x = yr_pred, y = max(likdat$hts_pos$tot / 1000, na.rm = TRUE)/20, pos = 2) } } else { plot(I(out_nbtest_pos$value / 1000) ~ out_nbtest_pos$year, pch = '', ylim = c(0, max(out_nbtest_pos$value, na.rm = TRUE) / 1000 * 1.25), - cex = 1, ylab = expression(paste(N^o, " of Positive Tests (in 1,000)")), - xlab = 'Year', xlim = c(2000, yr_pred), + cex = 1, ylab = expression(paste(N^o, " of Positive Tests (in 1,000)")), + xlab = 'Year', xlim = c(2000, yr_pred), main = main_title) - lines(I(out_nbtest_pos$value / 1000) ~ out_nbtest_pos$year, lwd = 1, col = 'orangered2') + graphics::lines(I(out_nbtest_pos$value / 1000) ~ out_nbtest_pos$year, lwd = 1, col = 'orangered2') if (!is.null(likdat$hts_pos$tot)) { if (cnt %in% redact) { pchpt <- '' } else { pchpt <- 16 } - points(I(likdat$hts_pos$tot / 1000) ~ I(likdat$hts_pos$year + 0.5), pch = pchpt) + graphics::points(I(likdat$hts_pos$tot / 1000) ~ I(likdat$hts_pos$year + 0.5), pch = pchpt) yield <- round(likdat$hts$totpos / likdat$hts$tot * 100, 1) y_range <- c(min(yield, na.rm = TRUE), max(yield, na.rm = TRUE)) p_pos_pop <- round(likdat$hts_pos$tot / plhiv[(likdat$hts_pos$year - start)] * 100, 0) p_pos_pop_r <- c(min(p_pos_pop, na.rm = TRUE), max(p_pos_pop, na.rm = TRUE)) - text(paste('Positivity = ', y_range[1], '-', y_range[2], '%', sep=''), + graphics::text(paste('Positivity = ', y_range[1], '-', y_range[2], '%', sep=''), x = yr_pred, y = max(likdat$hts_pos$tot / 1000, na.rm = TRUE) / 6, pos = 2) - text(paste('No Positive/PLHIV = ', p_pos_pop_r[1], '-', - p_pos_pop_r[2], '%', sep = ''), x = yr_pred, + graphics::text(paste('No Positive/PLHIV = ', p_pos_pop_r[1], '-', + p_pos_pop_r[2], '%', sep = ''), x = yr_pred, y = max(likdat$hts_pos$tot / 1000, na.rm = TRUE) / 20, pos = 2) } } } #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year plot_out_nbpostest_sex <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = 2022) { - + # redact <- c('Namibia','Uganda','Zambia','Zimbabwe') redact <- c('XXX') - + start <- fp$ss$proj_start mod <- simmod(fp) - plhiv <- apply(attr(mod, "hivpop")[,1:8,,], 4, FUN = sum) + + plhiv <- apply(attr(mod, "hivpop")[,1:8,,], 4, FUN = sum) + apply(attr(mod, "artpop")[,, 1:8,,], 5, FUN = sum) pop <- apply(mod[1:35,,,], 4, FUN = sum) - + out_nbtest_pos <- subset(get_out_nbtest_pos_sex(mod, fp), year <= yr_pred) out_nbtest_pos$year <- out_nbtest_pos$year + 0.5 out_nbtest_pos_f <- subset(out_nbtest_pos, sex == 'female') out_nbtest_pos_m <- subset(out_nbtest_pos, sex == 'male') max_ylim <- max(c(out_nbtest_pos_f$value, out_nbtest_pos_m$value), na.rm = TRUE) - + lik_sex <- likdat$hts_pos if (!is.null(lik_sex)){ lik_f <- subset(lik_sex, sex == 'female') @@ -455,57 +455,57 @@ plot_out_nbpostest_sex <- function(mod, fp, likdat, cnt, simul = NULL, yr_pred = } else { plot.ci <- FALSE } - - col_ci <- rgb(150, 150, 150, 125, max = 255) - + + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + # Postive tests if (plot.ci == T){ plot(I(out_nbtest_pos_f$value / 1000) ~ out_nbtest_pos_f$year, pch = '', ylim = c(0, max_ylim / 1000 * 1.25), - cex = 1, ylab = expression(paste(N^o, " of Positive Tests (in 1,000)")), - xlab = 'Year', xlim = c(2000, yr_pred), + cex = 1, ylab = expression(paste(N^o, " of Positive Tests (in 1,000)")), + xlab = 'Year', xlim = c(2000, yr_pred), main = expression(bold(paste("Total ", N^o, " of Positive Tests")))) - polygon(x = c(ci_f$year, rev(ci_f$year)), + graphics::polygon(x = c(ci_f$year, rev(ci_f$year)), y = c(ci_f$upper / 1000, rev(ci_f$lower / 1000)), col = col_ci, border = NA) - polygon(x = c(ci_m$year, rev(ci_m$year)), + graphics::polygon(x = c(ci_m$year, rev(ci_m$year)), y = c(ci_m$upper / 1000, rev(ci_m$lower / 1000)), col = col_ci, border = NA) - lines(I(out_nbtest_pos_f$value / 1000) ~ out_nbtest_pos_f$year, lwd = 1, col = 'orangered2') - lines(I(out_nbtest_pos_m$value / 1000) ~ out_nbtest_pos_m$year, lwd = 1, col = 'orangered2') - + graphics::lines(I(out_nbtest_pos_f$value / 1000) ~ out_nbtest_pos_f$year, lwd = 1, col = 'orangered2') + graphics::lines(I(out_nbtest_pos_m$value / 1000) ~ out_nbtest_pos_m$year, lwd = 1, col = 'orangered2') + if (length(lik_f) > 0 || length(lik_m) > 0) { if (cnt %in% redact) { pchpt <- '' } else { pchpt <- 16 } - points(I(lik_f$totpos / 1000) ~ I(lik_f$year + 0.5), pch = pchpt) - points(I(lik_m$totpos / 1000) ~ I(lik_m$year + 0.5), pch = pchpt + 1) + graphics::points(I(lik_f$totpos / 1000) ~ I(lik_f$year + 0.5), pch = pchpt) + graphics::points(I(lik_m$totpos / 1000) ~ I(lik_m$year + 0.5), pch = pchpt + 1) } } else { plot(I(out_nbtest_pos_f$value / 1000) ~ out_nbtest_pos_f$year, pch = '', ylim = c(0, max_ylim / 1000 * 1.25), - cex = 1, ylab = expression(paste(N^o, " of Positive Tests (in 1,000)")), - xlab = 'Year', xlim = c(2000, yr_pred), + cex = 1, ylab = expression(paste(N^o, " of Positive Tests (in 1,000)")), + xlab = 'Year', xlim = c(2000, yr_pred), main = expression(bold(paste("Total ", N^o, " of Positive Tests")))) - lines(I(out_nbtest_pos_f$value / 1000) ~ out_nbtest_pos_f$year, lwd = 1, col = 'orangered2') - lines(I(out_nbtest_pos_m$value / 1000) ~ out_nbtest_pos_m$year, lwd = 1, col = 'orangered2') + graphics::lines(I(out_nbtest_pos_f$value / 1000) ~ out_nbtest_pos_f$year, lwd = 1, col = 'orangered2') + graphics::lines(I(out_nbtest_pos_m$value / 1000) ~ out_nbtest_pos_m$year, lwd = 1, col = 'orangered2') if (!is.null(likdat$hts_pos$tot)) { if (cnt %in% redact) { pchpt <- '' } else { pchpt <- 16 } - points(I(lik_f$totpos / 1000) ~ I(lik_f$year + 0.5), pch = pchpt) - points(I(lik_m$totpos / 1000) ~ I(lik_m$year + 0.5), pch = pchpt + 1) + graphics::points(I(lik_f$totpos / 1000) ~ I(lik_f$year + 0.5), pch = pchpt) + graphics::points(I(lik_m$totpos / 1000) ~ I(lik_m$year + 0.5), pch = pchpt + 1) } } } #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year -plot_out_evertestneg <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, +plot_out_evertestneg <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, simul = NULL, plot_legend = TRUE, yr_pred = 2022) { - + out_evertest <- subset(out_evertest, year <= yr_pred) out_evertest$year <- out_evertest$year + 0.5 survey_hts$year <- survey_hts$year + 0.5 - + # Decide if we plot CI or not if (!is.null(simul)){ ci <- getCI(simul$ever.test) @@ -517,14 +517,14 @@ plot_out_evertestneg <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, } if (is.null(simul)) plot.ci <- FALSE - col_ci <- rgb(150, 150, 150, 125, max = 255) - + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + # Plots of Ever Test by Sex (among negatives) out_f <- subset(out_evertest, agegr == '15-49' & outcome == 'evertest' & sex =='female' & hivstatus == 'negative') out_m <- subset(out_evertest, agegr == '15-49' & outcome == 'evertest' & sex =='male' & hivstatus == 'negative') dat_f <- subset(survey_hts, country == cnt & outcome == 'evertest' & agegr == '15-49' & sex == 'female' & hivstatus == 'negative') dat_m <- subset(survey_hts, country == cnt & outcome == 'evertest' & agegr == '15-49' & sex == 'male' & hivstatus == 'negative') - + if (plot.ci == T){ ci_f <- subset(ci, agegr == '15-49' & outcome == 'evertest' & sex =='female' & hivstatus == 'negative') ci_m <- subset(ci, agegr == '15-49' & outcome == 'evertest' & sex =='male' & hivstatus == 'negative') @@ -532,51 +532,51 @@ plot_out_evertestneg <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, plot(I(out_f$value * 100) ~ out_f$year, pch = "", ylim = c(0, 100), col = "maroon", main = "Negative Ever Tested", xlim = c(2000, yr_pred), ylab = "Proportion Ever Tested Among Susceptibles", xlab = "Year", lwd = 1) - polygon( + graphics::polygon( x = c(ci_f$year, rev(ci_f$year)), y = c(I(ci_f$upper * 100), rev(I(ci_f$lower * 100))), col = col_ci, border = NA) - polygon( + graphics::polygon( x = c(ci_m$year, rev(ci_m$year)), y = c(I(ci_m$upper * 100), rev(I(ci_m$lower * 100))), col = col_ci, border = NA) - lines(I(out_f$value * 100) ~ out_f$year, col = "maroon", lwd = 1) - lines(I(out_m$value * 100) ~ out_m$year, col = "navy", lwd = 1) - points(I(dat_f$est * 100) ~ dat_f$year, pch = 15, col = "maroon") - points(I(dat_m$est * 100) ~ dat_m$year, pch = 16, col = "navy") - segments( + graphics::lines(I(out_f$value * 100) ~ out_f$year, col = "maroon", lwd = 1) + graphics::lines(I(out_m$value * 100) ~ out_m$year, col = "navy", lwd = 1) + graphics::points(I(dat_f$est * 100) ~ dat_f$year, pch = 15, col = "maroon") + graphics::points(I(dat_m$est * 100) ~ dat_m$year, pch = 16, col = "navy") + graphics::segments( x0 = dat_f$year, y0 = I(dat_f$ci_l * 100), x1 = dat_f$year, y1 = I(dat_f$ci_u * 100), lwd = 1, col = "maroon") - segments( + graphics::segments( x0 = dat_m$year, y0 = I(dat_m$ci_l * 100), x1 = dat_m$year, y1 = I(dat_m$ci_u * 100), lwd = 1, col = "navy") - if (plot_legend) { - legend(x = 2000, y = 100, - legend = c('Women (15-49 years)','Men (15-49 years)'), + if (plot_legend) { + graphics::legend(x = 2000, y = 100, + legend = c('Women (15-49 years)','Men (15-49 years)'), col = c('maroon','navy'), bty = 'n', lwd = 1, pch = c(15,16)) } } else { plot(I(out_f$value * 100) ~ out_f$year, type = "l", ylim = c(0, 100), col = "maroon", main = "Negative Ever Tested", xlab = "Year", lwd = 1, xlim = c(2000, yr_pred), ylab = "Proportion Ever Tested Among Susceptibles") - lines(I(out_m$value * 100) ~ out_m$year, col = "navy", lwd = 1) - points(I(dat_f$est * 100) ~ dat_f$year, pch = 15, col = "maroon") - points(I(dat_m$est * 100) ~ dat_m$year, pch = 16, col = "navy") - segments( + graphics::lines(I(out_m$value * 100) ~ out_m$year, col = "navy", lwd = 1) + graphics::points(I(dat_f$est * 100) ~ dat_f$year, pch = 15, col = "maroon") + graphics::points(I(dat_m$est * 100) ~ dat_m$year, pch = 16, col = "navy") + graphics::segments( x0 = dat_f$year, y0 = I(dat_f$ci_l * 100), x1 = dat_f$year, y1 = I(dat_f$ci_u * 100), lwd = 1, col = "maroon") - segments( + graphics::segments( x0 = dat_m$year, y0 = I(dat_m$ci_l * 100), x1 = dat_m$year, y1 = I(dat_m$ci_u * 100), lwd = 1, col = "navy") - if (plot_legend) { - legend(x = 2000, y = 100, - legend = c('Women (15-49 years)','Men (15-49 years)'), + if (plot_legend) { + graphics::legend(x = 2000, y = 100, + legend = c('Women (15-49 years)','Men (15-49 years)'), col = c('maroon','navy'), bty = 'n', lwd = 2, pch = c(15, 16)) } } } #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year plot_out_evertestpos <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, @@ -586,7 +586,7 @@ plot_out_evertestpos <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, out_evertest <- subset(out_evertest, year <= yr_pred) out_evertest$year <- out_evertest$year + 0.5 survey_hts$year <- survey_hts$year + 0.5 - + # Decide if we plot CI or not if (!is.null(simul)){ ci <- getCI(simul$ever.test) @@ -597,71 +597,71 @@ plot_out_evertestpos <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, plot.ci <- TRUE } if (is.null(simul)) plot.ci <- FALSE - + if (plot_title == TRUE) { main_title <- "PLHIV Ever Tested" } else { main_title <- "" } - col_ci <- rgb(150, 150, 150, 125, max = 255) - + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + # Plots of Ever Test by Sex (among PLHIV) out_f <- subset(out_evertest, agegr == '15-49' & outcome == 'evertest' & sex =='female' & hivstatus == 'positive') out_m <- subset(out_evertest, agegr == '15-49' & outcome == 'evertest' & sex =='male' & hivstatus == 'positive') dat_f <- subset(survey_hts, country == cnt & agegr == '15-49' & sex == 'female' & hivstatus == 'positive' & outcome == 'evertest') dat_m <- subset(survey_hts, country == cnt & agegr == '15-49' & sex == 'male' & hivstatus == 'positive' & outcome == 'evertest') - + if (plot.ci == T){ ci_f <- subset(ci, agegr == '15-49' & outcome == 'evertest' & sex =='female' & hivstatus == 'positive') ci_m <- subset(ci, agegr == '15-49' & outcome == 'evertest' & sex =='male' & hivstatus == 'positive') plot(I(out_f$value * 100) ~ out_f$year, pch = '', ylim = c(0, 100), col = 'maroon', main = main_title, xlim=c(2000, yr_pred), ylab = 'Proportion PLHIV Ever Tested', xlab = 'Year', lwd = 1) - polygon(x = c(ci_f$year, rev(ci_f$year)), + graphics::polygon(x = c(ci_f$year, rev(ci_f$year)), y = c(I(ci_f$upper * 100), rev(I(ci_f$lower*100))), col = col_ci, border = NA) - polygon(x = c(ci_m$year, rev(ci_m$year)), + graphics::polygon(x = c(ci_m$year, rev(ci_m$year)), y = c(I(ci_m$upper * 100), rev(I(ci_m$lower * 100))), col = col_ci, border = NA) - lines(I(out_f$value * 100) ~ out_f$year, col = 'maroon', lwd = 1) - lines(I(out_m$value * 100) ~ out_m$year, col = 'navy', lwd = 1) - points(I(dat_f$est * 100) ~ I(dat_f$year - 0.1), pch = 15, col = 'maroon') - points(I(dat_m$est * 100) ~ I(dat_m$year + 0.1), pch = 16, col = 'navy') - segments(x0 = dat_f$year - 0.1, y0 = I(dat_f$ci_l * 100), + graphics::lines(I(out_f$value * 100) ~ out_f$year, col = 'maroon', lwd = 1) + graphics::lines(I(out_m$value * 100) ~ out_m$year, col = 'navy', lwd = 1) + graphics::points(I(dat_f$est * 100) ~ I(dat_f$year - 0.1), pch = 15, col = 'maroon') + graphics::points(I(dat_m$est * 100) ~ I(dat_m$year + 0.1), pch = 16, col = 'navy') + graphics::segments(x0 = dat_f$year - 0.1, y0 = I(dat_f$ci_l * 100), x1 = dat_f$year - 0.1, y1 = I(dat_f$ci_u * 100), lwd = 1, col = 'maroon') - segments(x0 = dat_m$year + 0.1, y0 = I(dat_m$ci_l * 100), + graphics::segments(x0 = dat_m$year + 0.1, y0 = I(dat_m$ci_l * 100), x1 = dat_m$year + 0.1, y1 = I(dat_m$ci_u * 100), lwd = 1, col = 'navy') - if (plot_legend) { - legend(x = 2000, y = 100, - legend = c('Women (15-49 years)','Men (15-49 years)'), + if (plot_legend) { + graphics::legend(x = 2000, y = 100, + legend = c('Women (15-49 years)','Men (15-49 years)'), col = c('maroon','navy'), bty = 'n', lwd = 1, pch = c(15, 16)) } } else { plot(I(out_f$value * 100) ~ out_f$year, type='l', ylim = c(0,100), col = 'maroon', main = main_title, xlim = c(2000, yr_pred), ylab = 'Proportion PLHIV Ever Tested', xlab = 'Year', lwd = 1) - lines(I(out_m$value * 100) ~ out_m$year, col = 'navy', lwd = 1) - points(I(dat_f$est * 100) ~ I(dat_f$year - 0.1), pch = 15, col = 'maroon') - points(I(dat_m$est * 100) ~ I(dat_m$year + 0.1), pch = 16, col = 'navy') - segments(x0 = dat_f$year - 0.1, y0 = I(dat_f$ci_l * 100), + graphics::lines(I(out_m$value * 100) ~ out_m$year, col = 'navy', lwd = 1) + graphics::points(I(dat_f$est * 100) ~ I(dat_f$year - 0.1), pch = 15, col = 'maroon') + graphics::points(I(dat_m$est * 100) ~ I(dat_m$year + 0.1), pch = 16, col = 'navy') + graphics::segments(x0 = dat_f$year - 0.1, y0 = I(dat_f$ci_l * 100), x1 = dat_f$year - 0.1, y1 = I(dat_f$ci_u * 100), lwd = 1, col = 'maroon') - segments(x0 = dat_m$year + 0.1, y0 = I(dat_m$ci_l * 100), + graphics::segments(x0 = dat_m$year + 0.1, y0 = I(dat_m$ci_l * 100), x1 = dat_m$year + 0.1, y1 = I(dat_m$ci_u * 100), lwd = 1, col = 'navy') - if (plot_legend) { - legend(x = 2000, y = 100, - legend = c('Women (15-49 years)', 'Men (15-49 years)'), + if (plot_legend) { + graphics::legend(x = 2000, y = 100, + legend = c('Women (15-49 years)', 'Men (15-49 years)'), col = c('maroon', 'navy'), bty = 'n', lwd = 2, pch = c(15, 16)) } } } #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year -plot_out_evertest <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, +plot_out_evertest <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, simul = NULL, plot_legend = TRUE, yr_pred = 2022, plot_title = TRUE) { out_evertest <- subset(out_evertest, year <= yr_pred) out_evertest$year <- out_evertest$year + 0.5 survey_hts$year <- survey_hts$year + 0.5 - + # Decide if we plot CI or not if (!is.null(simul)){ ci <- getCI(simul$ever.test) @@ -672,7 +672,7 @@ plot_out_evertest <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, plot.ci <- TRUE } if (is.null(simul)) plot.ci <- FALSE - + if (plot_title == TRUE) { main_title <- "Population Ever Tested" } else { main_title <- "" } @@ -681,57 +681,57 @@ plot_out_evertest <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, out_m <- subset(out_evertest, agegr == '15-49' & outcome == 'evertest' & sex =='male' & hivstatus == 'all') dat_f <- subset(survey_hts, country == cnt & agegr == '15-49' & sex == 'female' & hivstatus == 'all' & outcome == 'evertest') dat_m <- subset(survey_hts, country == cnt & agegr == '15-49' & sex == 'male' & hivstatus == 'all' & outcome == 'evertest') - - col_ci <- rgb(150, 150, 150, 125, max = 255) - + + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + if (plot.ci == T){ ci_f <- subset(ci, agegr == '15-49' & outcome == 'evertest' & sex =='female' & hivstatus == 'all') ci_m <- subset(ci, agegr == '15-49' & outcome == 'evertest' & sex =='male' & hivstatus == 'all') plot(I(out_f$value * 100) ~ out_f$year, pch='', ylim = c(0,100), col = 'maroon', main = main_title, xlim=c(2000, yr_pred), ylab = 'Proportion Ever Tested', xlab = 'Year', lwd = 1) - polygon(x = c(ci_f$year, rev(ci_f$year)), + graphics::polygon(x = c(ci_f$year, rev(ci_f$year)), y = c(I(ci_f$upper * 100), rev(I(ci_f$lower*100))), col = col_ci, border = NA) - lines(I(out_f$value * 100) ~ out_m$year, col='maroon', lwd=1) - polygon(x = c(ci_m$year, rev(ci_m$year)), + graphics::lines(I(out_f$value * 100) ~ out_m$year, col='maroon', lwd=1) + graphics::polygon(x = c(ci_m$year, rev(ci_m$year)), y = c(I(ci_m$upper * 100), rev(I(ci_m$lower*100))), col = col_ci, border = NA) - lines(I(out_m$value * 100) ~ out_m$year, col = 'navy', lwd = 1) - points(I(dat_f$est * 100) ~ I(dat_f$year + 0.1), pch = 15, col='maroon') - points(I(dat_m$est * 100) ~ I(dat_m$year - 0.1), pch = 16, col = 'navy') - segments(x0 = dat_f$year + 0.1, y0 = I(dat_f$ci_l*100), + graphics::lines(I(out_m$value * 100) ~ out_m$year, col = 'navy', lwd = 1) + graphics::points(I(dat_f$est * 100) ~ I(dat_f$year + 0.1), pch = 15, col='maroon') + graphics::points(I(dat_m$est * 100) ~ I(dat_m$year - 0.1), pch = 16, col = 'navy') + graphics::segments(x0 = dat_f$year + 0.1, y0 = I(dat_f$ci_l*100), x1 = dat_f$year + 0.1, y1 = I(dat_f$ci_u*100), lwd = 1, col = 'maroon') - segments(x0 = dat_m$year - 0.1, y0 = I(dat_m$ci_l*100), + graphics::segments(x0 = dat_m$year - 0.1, y0 = I(dat_m$ci_l*100), x1 = dat_m$year - 0.1, y1 = I(dat_m$ci_u*100), lwd = 1, col = 'navy') - if (plot_legend) { - legend(x = 2000, y = 100, - legend = c('Women (15-49 years)','Men (15-49 years)'), + if (plot_legend) { + graphics::legend(x = 2000, y = 100, + legend = c('Women (15-49 years)','Men (15-49 years)'), col = c('maroon','navy'), bty = 'n', lwd = 1, pch = c(15,16)) } } else { - plot(I(out_f$value * 100) ~ out_f$year, type = 'l', ylim = c(0, 100), + plot(I(out_f$value * 100) ~ out_f$year, type = 'l', ylim = c(0, 100), col = 'maroon', main = main_title, - xlim = c(2000, yr_pred), ylab = 'Proportion Ever Tested', + xlim = c(2000, yr_pred), ylab = 'Proportion Ever Tested', xlab = 'Year', lwd = 1) - lines(I(out_m$value*100) ~ out_m$year, col = 'navy', lwd = 1) - points(I(dat_f$est*100) ~ I(dat_f$year + 0.1), pch = 15, col = 'maroon') - points(I(dat_m$est*100) ~ I(dat_m$year - 0.1), pch = 16, col = 'navy') - segments(x0 = dat_f$year + 0.1, y0 = I(dat_f$ci_l*100), + graphics::lines(I(out_m$value*100) ~ out_m$year, col = 'navy', lwd = 1) + graphics::points(I(dat_f$est*100) ~ I(dat_f$year + 0.1), pch = 15, col = 'maroon') + graphics::points(I(dat_m$est*100) ~ I(dat_m$year - 0.1), pch = 16, col = 'navy') + graphics::segments(x0 = dat_f$year + 0.1, y0 = I(dat_f$ci_l*100), x1 = dat_f$year + 0.1, y1 = I(dat_f$ci_u*100), lwd = 1, col = 'maroon') - segments(x0 = dat_m$year - 0.1, y0 = I(dat_m$ci_l*100), + graphics::segments(x0 = dat_m$year - 0.1, y0 = I(dat_m$ci_l*100), x1 = dat_m$year - 0.1, y1 = I(dat_m$ci_u*100), lwd = 1, col = 'navy') - if (plot_legend) { - legend(x = 2000, y = 100, - legend = c('Women (15-49 years)','Men (15-49 years)'), + if (plot_legend) { + graphics::legend(x = 2000, y = 100, + legend = c('Women (15-49 years)','Men (15-49 years)'), col = c('maroon','navy'), bty = 'n', lwd = 2, pch = c(15,16)) } } } #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year -plot_out_90s <- function(mod, fp, likdat, cnt, out_evertest, survey_hts, +plot_out_90s <- function(mod, fp, likdat, cnt, out_evertest, survey_hts, simul = NULL, plot_legend = TRUE, yr_pred = 2022) { phia_aware <- subset(survey_hts, country == cnt & agegr == '15-49' & @@ -740,16 +740,16 @@ plot_out_90s <- function(mod, fp, likdat, cnt, out_evertest, survey_hts, sex == 'both' & outcome == 'art_both') phia_aware$year <- phia_aware$year + 0.5 phia_art$year <- phia_art$year + 0.5 - + if (all(is.na(phia_art$est))) { phia_art <- subset(survey_hts, country == cnt & agegr == '15-49' & - sex == 'both' & outcome == 'art_self') + sex == 'both' & outcome == 'art_self') phia_art$year <- phia_art$year + 0.5 } svy_evertest <- subset(survey_hts, country == cnt & agegr == '15-49' & hivstatus == 'positive' & sex == 'both' & outcome == 'evertest') svy_evertest$year <- svy_evertest$year + 0.5 - + out_evertest <- subset(out_evertest, year <= yr_pred) out_aware <- subset(get_out_aware(mod, fp), year <= yr_pred) out_art <- subset(get_out_art(mod, fp), year >= 2000 & year <= yr_pred) @@ -771,15 +771,15 @@ plot_out_90s <- function(mod, fp, likdat, cnt, out_evertest, survey_hts, } if (is.null(simul)) plot.ci <- FALSE - col_ci <- rgb(150, 150, 150, 125, max = 255) - + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + # First 90 out_f <- subset(out_aware, agegr == '15-49' & outcome == 'aware' & sex =='female') out_m <- subset(out_aware, agegr == '15-49' & outcome == 'aware' & sex =='male') out_a <- subset(out_aware, agegr == '15-49' & outcome == 'aware' & sex =='both') out_ever_all <- subset(out_evertest, agegr == '15-49' & outcome == 'evertest' & sex == 'both' & hivstatus == 'positive') - + if (plot.ci == TRUE){ ci_f <- subset(ci$diagnoses, agegr == '15-49' & outcome == 'aware' & sex =='female') ci_m <- subset(ci$diagnoses, agegr == '15-49' & outcome == 'aware' & sex =='male') @@ -787,78 +787,78 @@ plot_out_90s <- function(mod, fp, likdat, cnt, out_evertest, survey_hts, ci_ever_all <- subset(ci$ever.test, agegr == '15-49' & outcome == 'evertest' & sex == 'both' & hivstatus == 'positive') - plot(I(out_a$value*100) ~ out_a$year, pch = '', ylim = c(0,100), + plot(I(out_a$value*100) ~ out_a$year, pch = '', ylim = c(0,100), col = 'darkslategrey', main = 'PLHIV Ever Tested, 1st 90, and ART Coverage', xlim = c(2000, yr_pred), ylab = 'Proportion', xlab='Year', lwd = 1) - polygon(x = c(ci_a$year, rev(ci_a$year)), + graphics::polygon(x = c(ci_a$year, rev(ci_a$year)), y = c(I(ci_a$upper*100), rev(I(ci_a$lower*100))), col = col_ci, border = NA) - lines(I(out_a$value*100) ~ out_a$year, col = 'darkslategrey', lwd = 1) - lines(I(out_art$value*100) ~ out_art$year, col = 'deepskyblue2', lwd = 1, lty = 1) - polygon(x = c(ci_ever_all$year, rev(ci_ever_all$year)), + graphics::lines(I(out_a$value*100) ~ out_a$year, col = 'darkslategrey', lwd = 1) + graphics::lines(I(out_art$value*100) ~ out_art$year, col = 'deepskyblue2', lwd = 1, lty = 1) + graphics::polygon(x = c(ci_ever_all$year, rev(ci_ever_all$year)), y = c(I(ci_ever_all$upper*100), rev(I(ci_ever_all$lower*100))), col = col_ci,border = NA) - lines(I(out_ever_all$value*100) ~ out_ever_all$year, col = 'orange3', lwd = 1, lty = 1) - legend(x = 2000, y = 100, legend = c('PLHIV Ever Tested','PLHIV Aware','ART Coverage'), + graphics::lines(I(out_ever_all$value*100) ~ out_ever_all$year, col = 'orange3', lwd = 1, lty = 1) + graphics::legend(x = 2000, y = 100, legend = c('PLHIV Ever Tested','PLHIV Aware','ART Coverage'), col = c('orange3','darkslategrey','deepskyblue2'), pch = c(25, 1, 2), lwd = 1, bty = 'n') # We plot ever tested if (dim(svy_evertest)[1] > 0) { - points(I(svy_evertest$est*100) ~ svy_evertest$year, pch = 25, bg = 'grey15') - segments(x0 = svy_evertest$year, y0 = I(svy_evertest$ci_l*100), + graphics::points(I(svy_evertest$est*100) ~ svy_evertest$year, pch = 25, bg = 'grey15') + graphics::segments(x0 = svy_evertest$year, y0 = I(svy_evertest$ci_l*100), x1 = svy_evertest$year, y1 = I(svy_evertest$ci_u*100), lwd = 1, col = 'bisque4') } # We plot PHIA if any if (dim(phia_aware)[1] > 0) { - points(I(phia_aware$est*100) ~ phia_aware$year, pch = 1) - segments(x0 = phia_aware$year, y0 = I(phia_aware$ci_l*100), + graphics::points(I(phia_aware$est*100) ~ phia_aware$year, pch = 1) + graphics::segments(x0 = phia_aware$year, y0 = I(phia_aware$ci_l*100), x1 = phia_aware$year, y1 = I(phia_aware$ci_u*100), lwd = 1, col = 'bisque4') - text('Empty symbols are \n direct survey measures \n (not included in fitting)', x = 2015, y = 12.5, cex = 0.8) + graphics::text('Empty symbols are \n direct survey measures \n (not included in fitting)', x = 2015, y = 12.5, cex = 0.8) if (any(cnt == c('Swaziland','Zambia','Zimbabwe'))) { - text('(not ARV-corrected)', x = 2015, y = 1, cex = 0.8) } + graphics::text('(not ARV-corrected)', x = 2015, y = 1, cex = 0.8) } } if (dim(phia_art)[1] > 0) { - points(I(phia_art$est*100) ~ phia_art$year, pch = 2) - segments(x0 = phia_art$year, y0 = I(phia_art$ci_l*100), + graphics::points(I(phia_art$est*100) ~ phia_art$year, pch = 2) + graphics::segments(x0 = phia_art$year, y0 = I(phia_art$ci_l*100), x1 = phia_art$year, y1 = I(phia_art$ci_u*100), lwd = 1, col = 'bisque4') } } else { - plot(I(out_a$value*100) ~ out_a$year, type = 'l', ylim = c(0,100), + plot(I(out_a$value*100) ~ out_a$year, type = 'l', ylim = c(0,100), col = 'darkslategrey', main='PLHIV Ever Tested, 1st 90, and ART Coverage', xlim = c(2000, yr_pred), ylab = 'Proportion', xlab = 'Year', lwd = 1) - lines(I(out_art$value*100) ~ out_art$year, col = 'deepskyblue2', lwd = 1, lty = 1) - lines(I(out_ever_all$value*100) ~ out_ever_all$year, col='orange3', lwd = 1, lty = 1) - legend(x = 2000, y = 100, legend = c('PLHIV Ever Tested','PLHIV Aware','ART Coverage'), + graphics::lines(I(out_art$value*100) ~ out_art$year, col = 'deepskyblue2', lwd = 1, lty = 1) + graphics::lines(I(out_ever_all$value*100) ~ out_ever_all$year, col='orange3', lwd = 1, lty = 1) + graphics::legend(x = 2000, y = 100, legend = c('PLHIV Ever Tested','PLHIV Aware','ART Coverage'), col = c('orange3','darkslategrey','deepskyblue2'), pch = c(25, 1, 2), lwd = 1, bty = 'n') # We plot ever tested if (dim(svy_evertest)[1] > 0) { - points(I(svy_evertest$est*100) ~ svy_evertest$year, pch = 25, bg = 'grey15') - segments(x0 = svy_evertest$year, y0 = I(svy_evertest$ci_l*100), + graphics::points(I(svy_evertest$est*100) ~ svy_evertest$year, pch = 25, bg = 'grey15') + graphics::segments(x0 = svy_evertest$year, y0 = I(svy_evertest$ci_l*100), x1 = svy_evertest$year, y1 = I(svy_evertest$ci_u*100), lwd = 1, col = 'bisque4') } # We plot PHIA if any if (dim(phia_aware)[1]>0) { - points(I(phia_aware$est*100) ~ phia_aware$year, pch = 1) - segments(x0 = phia_aware$year, y0 = I(phia_aware$ci_l*100), + graphics::points(I(phia_aware$est*100) ~ phia_aware$year, pch = 1) + graphics::segments(x0 = phia_aware$year, y0 = I(phia_aware$ci_l*100), x1 = phia_aware$year, y1 = I(phia_aware$ci_u*100), lwd = 1, col = 'bisque4') - text('X-Validation', x = 2015, y = 12.5, cex = 0.8) + graphics::text('X-Validation', x = 2015, y = 12.5, cex = 0.8) if (any(cnt == c('Swaziland','Zambia','Zimbabwe'))) { - text('(not ARV-corrected)', x = 2015, y = 1, cex = 0.8) } + graphics::text('(not ARV-corrected)', x = 2015, y = 1, cex = 0.8) } } if (dim(phia_art)[1]>0) { - points(I(phia_art$est*100) ~ phia_art$year, pch = 2) - segments(x0 = phia_art$year, y0 = I(phia_art$ci_l*100), + graphics::points(I(phia_art$est*100) ~ phia_art$year, pch = 2) + graphics::segments(x0 = phia_art$year, y0 = I(phia_art$ci_l*100), x1 = phia_art$year, y1 = I(phia_art$ci_u*100), lwd = 1, col = 'bisque4') } } } #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year -plot_out_evertest_fbyage <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, +plot_out_evertest_fbyage <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, simul = NULL, plot_legend = TRUE, yr_pred = 2022) { out_evertest <- subset(out_evertest, year <= yr_pred) out_evertest$year <- out_evertest$year + 0.5 survey_hts$year <- survey_hts$year + 0.5 - + # Decide if we plot CI or not if (!is.null(simul)){ ci <- getCI(simul$ever.test) @@ -880,79 +880,79 @@ plot_out_evertest_fbyage <- function(mod, fp, likdat, cnt, survey_hts, out_evert agegr == "25-34" & sex == "female" & outcome == "evertest") datage35f <- subset(survey_hts, country == cnt & hivstatus == "all" & agegr == "35-49" & sex == "female" & outcome == "evertest") - - col_ci <- rgb(150, 150, 150, 125, max = 255) - + + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + if (plot.ci == T){ ci_1524f <- subset(ci, agegr == '15-24' & outcome == 'evertest' & hivstatus=='all' & sex=='female') ci_2534f <- subset(ci, agegr == '25-34' & outcome == 'evertest' & hivstatus=='all' & sex=='female') ci_3549f <- subset(ci, agegr == '35-49' & outcome == 'evertest' & hivstatus=='all' & sex=='female') - plot(I(out_1524f$value*100) ~ out_1524f$year, pch = '', ylim = c(0,100), + plot(I(out_1524f$value*100) ~ out_1524f$year, pch = '', ylim = c(0,100), xlim=c(2000, yr_pred), - ylab = 'Proportion ever tested', xlab = 'Year', + ylab = 'Proportion ever tested', xlab = 'Year', main = 'Women Ever Tested, by Age') - polygon(x = c(ci_1524f$year, rev(ci_1524f$year)), + graphics::polygon(x = c(ci_1524f$year, rev(ci_1524f$year)), y = c(I(ci_1524f$upper*100), rev(I(ci_1524f$lower*100))), col = col_ci, border = NA) - polygon(x = c(ci_2534f$year, rev(ci_2534f$year)), + graphics::polygon(x = c(ci_2534f$year, rev(ci_2534f$year)), y = c(I(ci_2534f$upper*100), rev(I(ci_2534f$lower*100))), col = col_ci, border = NA) - polygon(x = c(ci_3549f$year, rev(ci_3549f$year)), + graphics::polygon(x = c(ci_3549f$year, rev(ci_3549f$year)), y = c(I(ci_3549f$upper*100), rev(I(ci_3549f$lower*100))), col = col_ci, border = NA) - lines(I(out_1524f$value*100) ~ out_1524f$year, col = 'maroon', lwd = 1, lty = 3) - lines(I(out_2534f$value*100) ~ out_2534f$year, col = 'maroon', lwd = 1, lty = 2) - lines(I(out_3549f$value*100) ~ out_3549f$year, col = 'maroon', lwd = 1, lty = 1) - if (plot_legend) { - legend(x = 2000, y = 100, - legend = c('15-24 years', '25-34 years','35-49 years'), + graphics::lines(I(out_1524f$value*100) ~ out_1524f$year, col = 'maroon', lwd = 1, lty = 3) + graphics::lines(I(out_2534f$value*100) ~ out_2534f$year, col = 'maroon', lwd = 1, lty = 2) + graphics::lines(I(out_3549f$value*100) ~ out_3549f$year, col = 'maroon', lwd = 1, lty = 1) + if (plot_legend) { + graphics::legend(x = 2000, y = 100, + legend = c('15-24 years', '25-34 years','35-49 years'), col = c('maroon','maroon'), lwd = 2, lty = c(3,2,1), pch = c(16,15,17), pt.cex = 0.7, bty = 'n') } - points(I(datage15f$est*100) ~ I(datage15f$year-0.1), pch = 16, col = 'goldenrod3', cex = 0.8) - points(I(datage25f$est*100) ~ datage25f$year, pch=15, col = 'goldenrod3', cex = 0.8) - points(I(datage35f$est*100) ~ I(datage35f$year+0.1), pch = 17, col = 'goldenrod3', cex = 0.8) - segments(x0 = datage15f$year-0.1, y0 = I(datage15f$ci_l*100), + graphics::points(I(datage15f$est*100) ~ I(datage15f$year-0.1), pch = 16, col = 'goldenrod3', cex = 0.8) + graphics::points(I(datage25f$est*100) ~ datage25f$year, pch=15, col = 'goldenrod3', cex = 0.8) + graphics::points(I(datage35f$est*100) ~ I(datage35f$year+0.1), pch = 17, col = 'goldenrod3', cex = 0.8) + graphics::segments(x0 = datage15f$year-0.1, y0 = I(datage15f$ci_l*100), x1 = datage15f$year-0.1, y1 = I(datage15f$ci_u*100), col = 'goldenrod3', lwd=1) - segments(x0 = datage25f$year, y0 = I(datage25f$ci_l*100), + graphics::segments(x0 = datage25f$year, y0 = I(datage25f$ci_l*100), x1 = datage25f$year, y1 = I(datage25f$ci_u*100), col = 'goldenrod3', lwd = 1) - segments(x0 = datage35f$year+0.1, y0 = I(datage35f$ci_l*100), + graphics::segments(x0 = datage35f$year+0.1, y0 = I(datage35f$ci_l*100), x1 = datage35f$year+0.1, y1 = I(datage35f$ci_u*100), col = 'goldenrod3', lwd = 1) } else { plot(I(out_1524f$value*100) ~ out_1524f$year, type = 'l', col = 'maroon', ylim = c(0,100), xlim = c(2000, yr_pred), lwd = 1, lty = 3, - ylab = 'Proportion ever tested', xlab = 'Year', + ylab = 'Proportion ever tested', xlab = 'Year', main = 'Women Ever Tested, by Age') - lines(I(out_2534f$value*100) ~ out_2534f$year, col = 'maroon', lwd = 1, lty = 2) - lines(I(out_3549f$value*100) ~ out_3549f$year, col = 'maroon', lwd = 1, lty = 1) - if (plot_legend) { - legend(x = 2000, y = 100, - legend = c('15-24 years', '25-34 years','35-49 years'), + graphics::lines(I(out_2534f$value*100) ~ out_2534f$year, col = 'maroon', lwd = 1, lty = 2) + graphics::lines(I(out_3549f$value*100) ~ out_3549f$year, col = 'maroon', lwd = 1, lty = 1) + if (plot_legend) { + graphics::legend(x = 2000, y = 100, + legend = c('15-24 years', '25-34 years','35-49 years'), col = c('maroon','maroon'), lwd = 2, lty = c(3,2,1), pch = c(16,15,17), pt.cex = 0.7, bty = 'n') } - points(I(datage15f$est*100) ~ I(datage15f$year-0.1), pch = 16, col = 'goldenrod3', cex = 0.8) - points(I(datage25f$est*100) ~ datage25f$year, pch = 15, col = 'goldenrod3', cex = 0.8) - points(I(datage35f$est*100) ~ I(datage35f$year+0.1), pch = 17, col = 'goldenrod3', cex = 0.8) - segments(x0 = datage15f$year-0.1, y0 = I(datage15f$ci_l*100), + graphics::points(I(datage15f$est*100) ~ I(datage15f$year-0.1), pch = 16, col = 'goldenrod3', cex = 0.8) + graphics::points(I(datage25f$est*100) ~ datage25f$year, pch = 15, col = 'goldenrod3', cex = 0.8) + graphics::points(I(datage35f$est*100) ~ I(datage35f$year+0.1), pch = 17, col = 'goldenrod3', cex = 0.8) + graphics::segments(x0 = datage15f$year-0.1, y0 = I(datage15f$ci_l*100), x1 = datage15f$yea-0.1, y1 = I(datage15f$ci_u*100), col = 'goldenrod3', lwd = 1) - segments(x0 = datage25f$year, y0 = I(datage25f$ci_l*100), + graphics::segments(x0 = datage25f$year, y0 = I(datage25f$ci_l*100), x1 = datage25f$year, y1 = I(datage25f$ci_u*100), col = 'goldenrod3', lwd = 1) - segments(x0 = datage35f$year+0.1, y0 = I(datage35f$ci_l*100), + graphics::segments(x0 = datage35f$year+0.1, y0 = I(datage35f$ci_l*100), x1 = datage35f$year+0.1, y1 = I(datage35f$ci_u*100), col = 'goldenrod3', lwd = 1) } } #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year -plot_out_evertest_mbyage <- function(mod, fp, likdat, cnt, survey_hts, - out_evertest, simul = NULL, +plot_out_evertest_mbyage <- function(mod, fp, likdat, cnt, survey_hts, + out_evertest, simul = NULL, plot_legend = TRUE, yr_pred = 2022) { out_evertest <- subset(out_evertest, year <= yr_pred) out_evertest$year <- out_evertest$year + 0.5 survey_hts$year <- survey_hts$year + 0.5 - + # Decide if we plot CI or not if (!is.null(simul)){ ci <- getCI(simul$ever.test) @@ -974,61 +974,61 @@ plot_out_evertest_mbyage <- function(mod, fp, likdat, cnt, survey_hts, datage35m <- subset(survey_hts, country == cnt & hivstatus == "all" & agegr == "35-49" & sex == "male" & outcome == "evertest") - col_ci <- rgb(150, 150, 150, 125, max = 255) - + col_ci <- grDevices::rgb(150, 150, 150, 125, max = 255) + if (plot.ci == T){ ci_1524m <- subset(ci, agegr == '15-24' & outcome == 'evertest' & hivstatus=='all' & sex=='male') ci_2534m <- subset(ci, agegr == '25-34' & outcome == 'evertest' & hivstatus=='all' & sex=='male') ci_3549m <- subset(ci, agegr == '35-49' & outcome == 'evertest' & hivstatus=='all' & sex=='male') - plot(I(out_1524m$value*100) ~ out_1524m$year, pch = '', + plot(I(out_1524m$value*100) ~ out_1524m$year, pch = '', ylim = c(0,100), xlim = c(2000, yr_pred), ylab = 'Proportion ever tested', xlab = 'Year', main = 'Men Ever Tested, by Age') - polygon(x = c(ci_1524m$year, rev(ci_1524m$year)), + graphics::polygon(x = c(ci_1524m$year, rev(ci_1524m$year)), y = c(I(ci_1524m$upper*100), rev(I(ci_1524m$lower*100))), col = col_ci, border = NA) - polygon(x = c(ci_2534m$year, rev(ci_2534m$year)), + graphics::polygon(x = c(ci_2534m$year, rev(ci_2534m$year)), y = c(I(ci_2534m$upper*100), rev(I(ci_2534m$lower*100))), col = col_ci, border = NA) - polygon(x = c(ci_3549m$year, rev(ci_3549m$year)), + graphics::polygon(x = c(ci_3549m$year, rev(ci_3549m$year)), y = c(I(ci_3549m$upper*100), rev(I(ci_3549m$lower*100))), col = col_ci, border = NA) - lines(I(out_1524m$value*100) ~ I(out_1524m$year-0.1), col = 'navy', lwd = 1, lty = 3) - lines(I(out_2534m$value*100) ~ out_2534m$year, col = 'navy', lwd = 1, lty = 2) - lines(I(out_3549m$value*100) ~ I(out_3549m$year+0.1), col = 'navy', lwd = 1, lty = 1) - if (plot_legend) { - legend(x = 2000, y = 100, + graphics::lines(I(out_1524m$value*100) ~ I(out_1524m$year-0.1), col = 'navy', lwd = 1, lty = 3) + graphics::lines(I(out_2534m$value*100) ~ out_2534m$year, col = 'navy', lwd = 1, lty = 2) + graphics::lines(I(out_3549m$value*100) ~ I(out_3549m$year+0.1), col = 'navy', lwd = 1, lty = 1) + if (plot_legend) { + graphics::legend(x = 2000, y = 100, legend = c('15-24', '25-34', '35-49'), col = c('navy','navy'), lwd = 2, lty = c(3,2,1), pch = c(16,15,17), pt.cex = 0.7, bty = 'n') } - points(I(datage15m$est*100) ~ I(datage15m$year-0.1), pch = 16, col = 'goldenrod3', cex = 0.8) - points(I(datage25m$est*100) ~ datage25m$year, pch = 15, col = 'goldenrod3', cex = 0.8) - points(I(datage35m$est*100) ~ I(datage25m$year+0.1), pch = 17, col = 'goldenrod3', cex = 0.8) - segments(x0 = datage15m$year-0.1, y0 = I(datage15m$ci_l*100), + graphics::points(I(datage15m$est*100) ~ I(datage15m$year-0.1), pch = 16, col = 'goldenrod3', cex = 0.8) + graphics::points(I(datage25m$est*100) ~ datage25m$year, pch = 15, col = 'goldenrod3', cex = 0.8) + graphics::points(I(datage35m$est*100) ~ I(datage25m$year+0.1), pch = 17, col = 'goldenrod3', cex = 0.8) + graphics::segments(x0 = datage15m$year-0.1, y0 = I(datage15m$ci_l*100), x1 = datage15m$year-0.1, y1 = I(datage15m$ci_u*100), col = 'goldenrod3', lwd = 1) - segments(x0 = datage25m$year, y0 = I(datage25m$ci_l*100), + graphics::segments(x0 = datage25m$year, y0 = I(datage25m$ci_l*100), x1=datage25m$year, y1 = I(datage25m$ci_u*100), col = 'goldenrod3', lwd = 1) - segments(x0 = datage35m$year+0.1, y0 = I(datage35m$ci_l*100), + graphics::segments(x0 = datage35m$year+0.1, y0 = I(datage35m$ci_l*100), x1 = datage35m$year+0.1, y1 = I(datage35m$ci_u*100), col = 'goldenrod3', lwd = 1) } else { - plot(I(out_1524m$value*100) ~ out_1524m$year, type = 'l', col = 'navy', + plot(I(out_1524m$value*100) ~ out_1524m$year, type = 'l', col = 'navy', ylim = c(0,100), xlim = c(2000, yr_pred), lwd = 1, lty = 3, - ylab = 'Proportion ever tested', xlab = 'Year', + ylab = 'Proportion ever tested', xlab = 'Year', main = 'Men Ever Tested, by Age') - lines(I(out_2534m$value*100) ~ out_2534m$year, col = 'navy', lwd = 1, lty = 2) - lines(I(out_3549m$value*100) ~ out_3549m$year, col = 'navy', lwd = 1, lty = 1) - if (plot_legend) { - legend(x = 2000, y = 100, + graphics::lines(I(out_2534m$value*100) ~ out_2534m$year, col = 'navy', lwd = 1, lty = 2) + graphics::lines(I(out_3549m$value*100) ~ out_3549m$year, col = 'navy', lwd = 1, lty = 1) + if (plot_legend) { + graphics::legend(x = 2000, y = 100, legend = c('15-24', '25-34', '35-49'), col = c('navy','navy'), lwd = 2, lty = c(3,2,1), pch = c(16,15,17), pt.cex = 0.7, bty = 'n') } - points(I(datage15m$est*100) ~ I(datage15m$year-0.1), pch = 16, col = 'goldenrod3', cex = 0.8) - points(I(datage25m$est*100) ~ datage25m$year, pch = 15, col = 'goldenrod3', cex = 0.8) - points(I(datage35m$est*100) ~ I(datage25m$year+0.1), pch = 17, col = 'goldenrod3', cex = 0.8) - segments(x0=datage15m$year-0.1, y0=I(datage15m$ci_l*100), + graphics::points(I(datage15m$est*100) ~ I(datage15m$year-0.1), pch = 16, col = 'goldenrod3', cex = 0.8) + graphics::points(I(datage25m$est*100) ~ datage25m$year, pch = 15, col = 'goldenrod3', cex = 0.8) + graphics::points(I(datage35m$est*100) ~ I(datage25m$year+0.1), pch = 17, col = 'goldenrod3', cex = 0.8) + graphics::segments(x0=datage15m$year-0.1, y0=I(datage15m$ci_l*100), x1=datage15m$year-0.1, y1=I(datage15m$ci_u*100), col = 'goldenrod3', lwd = 1) - segments(x0=datage25m$year, y0=I(datage25m$ci_l*100), + graphics::segments(x0=datage25m$year, y0=I(datage25m$ci_l*100), x1=datage25m$year, y1=I(datage25m$ci_u*100), col = 'goldenrod3', lwd = 1) - segments(x0=datage35m$year+0.1, y0=I(datage35m$ci_l*100), + graphics::segments(x0=datage35m$year+0.1, y0=I(datage35m$ci_l*100), x1=datage35m$year+0.1, y1=I(datage25m$ci_u*100), col = 'goldenrod3', lwd = 1) } } @@ -1039,49 +1039,49 @@ plot_out_evertest_mbyage <- function(mod, fp, likdat, cnt, survey_hts, #' @export ## -- UPDATE HERE -- ## * update yr_pred to current year -plot_out <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, simul = NULL, +plot_out <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, simul = NULL, plot_legend = TRUE, yr_pred = 2022) { -par(mfrow = c(3,2), mar = c(4,4,2,2)) +graphics::par(mfrow = c(3,2), mar = c(4,4,2,2)) plot_out_nbtest(mod, fp, likdat, cnt, simul, yr_pred) - plot_out_nbpostest(mod, fp, likdat, cnt, simul, yr_pred) + plot_out_nbpostest(mod, fp, likdat, cnt, simul, yr_pred) plot_out_evertestneg(mod, fp, likdat, cnt, survey_hts, out_evertest, simul, plot_legend, yr_pred) plot_out_evertestpos(mod, fp, likdat, cnt, survey_hts, out_evertest, simul, plot_legend, yr_pred) plot_out_evertest(mod, fp, likdat, cnt, survey_hts, out_evertest, simul, plot_legend, yr_pred) - plot_out_90s(mod, fp, likdat, cnt, out_evertest, survey_hts, simul, plot_legend, yr_pred) + plot_out_90s(mod, fp, likdat, cnt, out_evertest, survey_hts, simul, plot_legend, yr_pred) } #' @export -## +## ## -- UPDATE HERE -- ## * update yr_pred to current year -plot_out_strat <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, simul = NULL, +plot_out_strat <- function(mod, fp, likdat, cnt, survey_hts, out_evertest, simul = NULL, plot_legend = TRUE, yr_pred = 2022) { - par(mfrow = c(1,2), mar = c(4,4,2,2)) - plot_out_evertest_mbyage(mod, fp, likdat, cnt, survey_hts, out_evertest, simul, plot_legend, yr_pred) - plot_out_evertest_fbyage(mod, fp, likdat, cnt, survey_hts, out_evertest, simul, plot_legend, yr_pred) + graphics::par(mfrow = c(1,2), mar = c(4,4,2,2)) + plot_out_evertest_mbyage(mod, fp, likdat, cnt, survey_hts, out_evertest, simul, plot_legend, yr_pred) + plot_out_evertest_fbyage(mod, fp, likdat, cnt, survey_hts, out_evertest, simul, plot_legend, yr_pred) } # ---- Output Tabular Format ---- #' @export end_of_year <- function(year, value){ - if (length(unique(year)) != length(year)) { print('non unique years'); break } + if (length(unique(year)) != length(year)) { print('non unique years'); return() } new_x <- year + 0.5 - new_value <- approx(year, value, new_x, method = 'linear', rule = 2)$y + new_value <- stats::approx(year, value, new_x, method = 'linear', rule = 2)$y return(new_value) -} +} #' @export -## +## ## -- UPDATE HERE -- ## * update year_range to include current year -tab_out_evertest <- function(mod, fp, age_grp = '15-49', gender = 'both', - hiv = 'all', year_range = c(2010, 2022), +tab_out_evertest <- function(mod, fp, age_grp = '15-49', gender = 'both', + hiv = 'all', year_range = c(2010, 2022), simul = NULL, end_year = TRUE) { interpolate_output <- end_year && fp$projection_period == "midyear" || !end_year && fp$projection_period == "calendar" - + if (length(year_range) == 1) { year_range <- c(year_range, year_range) } if (is.null(simul)) { out <- get_out_evertest(mod, fp, age_grp, gender, hiv) @@ -1104,7 +1104,7 @@ tab_out_evertest <- function(mod, fp, age_grp = '15-49', gender = 'both', } outci$lower <- round(outci$lower * 100, 1) outci$upper <- round(outci$upper * 100, 1) - outci <- subset(outci, year >= year_range[1] & year <= year_range[2]) + outci <- subset(outci, year >= year_range[1] & year <= year_range[2]) tab_evertest <- merge(out, outci) } row.names(tab_evertest) <- NULL @@ -1112,11 +1112,11 @@ tab_out_evertest <- function(mod, fp, age_grp = '15-49', gender = 'both', } #' @export -## +## ## -- UPDATE HERE -- ## * update year_range to include current year -tab_out_aware <- function(mod, fp, age_grp = '15-49', gender = 'both', - year_range = c(2010, 2022), simul = NULL, +tab_out_aware <- function(mod, fp, age_grp = '15-49', gender = 'both', + year_range = c(2010, 2022), simul = NULL, end_year = TRUE) { interpolate_output <- end_year && fp$projection_period == "midyear" || @@ -1125,7 +1125,7 @@ tab_out_aware <- function(mod, fp, age_grp = '15-49', gender = 'both', if (age_grp == "15-99") { age_grp <- "15+" } - + if (length(year_range) == 1) { year_range <- c(year_range, year_range) } @@ -1144,30 +1144,30 @@ tab_out_aware <- function(mod, fp, age_grp = '15-49', gender = 'both', } out$value <- round(out$value * 100, 1) outci <- getCI(simul$diagnoses) - outci <- subset(outci, agegr == age_grp & sex == gender) + outci <- subset(outci, agegr == age_grp & sex == gender) - if (interpolate_output) { + if (interpolate_output) { outci$lower <- end_of_year(outci$year, outci$lower) outci$upper <- end_of_year(outci$year, outci$upper) } outci$lower <- round(outci$lower * 100, 1) outci$upper <- round(outci$upper * 100, 1) - outci <- subset(outci, year >= year_range[1] & year <= year_range[2]) + outci <- subset(outci, year >= year_range[1] & year <= year_range[2]) tab_aware <- merge(out, outci) } - + row.names(tab_aware) <- NULL tab_aware } #' @export -## +## ## -- UPDATE HERE -- ## * update year_range to include current year -tab_out_nbaware <- function(mod, fp, age_grp = '15-49', - gender = 'both', year_range = c(2010, 2022), +tab_out_nbaware <- function(mod, fp, age_grp = '15-49', + gender = 'both', year_range = c(2010, 2022), end_year = TRUE) { interpolate_output <- end_year && fp$projection_period == "midyear" || @@ -1191,13 +1191,13 @@ tab_out_nbaware <- function(mod, fp, age_grp = '15-49', #' @export -## +## ## -- UPDATE HERE -- ## * update year_range to include current year -tab_out_artcov <- function(mod, fp, gender = 'both', +tab_out_artcov <- function(mod, fp, gender = 'both', year_range = c(2010, 2022)) { ## ART coverage is already end-of-year, no need to adjust - + if (length(year_range) == 1) { year_range <- c(year_range, year_range) } @@ -1232,10 +1232,10 @@ tab_out_artcov <- function(mod, fp, gender = 'both', } #' @export -## +## ## -- UPDATE HERE -- ## * update year_range to include current year -tab_out_pregprev <- function(mod, fp, year_range = c(2010, 2022), +tab_out_pregprev <- function(mod, fp, year_range = c(2010, 2022), end_year = TRUE) { if (length(year_range) == 1) { diff --git a/R/retest.R b/R/retest.R index e3b8476..ac3ce16 100644 --- a/R/retest.R +++ b/R/retest.R @@ -120,8 +120,8 @@ plot_retest_test_neg <- function(mod, fp, likdat, cnt, relative = F, add_ss_indices(out_retest, fp$ss))$tests out_retest$first_test <- out_retest$tests - out_retest$retests - col1 <- rgb(230, 200, 0, 200, max=255) - col2 <- rgb(50, 100, 155, 200, max=255) + col1 <- grDevices::rgb(230, 200, 0, 200, max=255) + col2 <- grDevices::rgb(50, 100, 155, 200, max=255) if(relative == F) { if (plot_title == TRUE) { main_title <- expression(bold(paste("Total ", N^o, " of HIV- Tests (in 1,000)"))) @@ -131,13 +131,13 @@ plot_retest_test_neg <- function(mod, fp, likdat, cnt, relative = F, main = main_title, xlim = c(2000, yr_pred), ylim = c(0, max(out_retest$tests/1000)*1.2), ylab=expression(paste(N^o, " of HIV- Tests (in 1,000)")), xlab = 'Year', lty = 1, lwd = 1) - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c(rep(0, length = length(out_retest$tests)), rev(out_retest$first_test/1000)), col = col1, border = NA) - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c(out_retest$tests/1000, rev(out_retest$first_test/1000)), col = col2, border = NA) - legend('topleft', inset = 0.02, legend = c('Repeat Testers','First-Time Testers'), + graphics::legend('topleft', inset = 0.02, legend = c('Repeat Testers','First-Time Testers'), col = c(col2, col1), bty = 'n', pch = c(15,15), pt.cex = 2, cex = 0.9) } else if(relative == T) { if (plot_title == TRUE) { main_title <- "Proportion of HIV- Tests \n Among Never and Ever Tested" @@ -146,13 +146,13 @@ plot_retest_test_neg <- function(mod, fp, likdat, cnt, relative = F, plot((out_retest$tests/out_retest$tests)*100 ~ out_retest$year, type='l', ylim=c(0,100), col=col1, main = main_title, xlim = c(2000, yr_pred), ylab='Proportion of Total HIV- Tests (%)', xlab='Year', lty=1, lwd=1) - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c((out_retest$retest/out_retest$tests)*100, rep(100, length = length(out_retest$tests))), col = col1, border = NA) - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c(rep(0, length = length(out_retest$tests)), rev((out_retest$retest/out_retest$tests)*100)), col = col2, border = NA) - legend('topleft', inset = 0.02, legend = c('% Repeat Testers','% First-Time Testers'), + graphics::legend('topleft', inset = 0.02, legend = c('% Repeat Testers','% First-Time Testers'), col = c(col2, col1), bg = 'grey90', box.col = 'grey90', pch = c(15,15), pt.cex = 2, cex = 0.8) } else { print('True for relative scale, False for absolute scale') @@ -184,9 +184,9 @@ plot_retest_test_pos <- function(mod, fp, likdat, cnt, relative = F, ylim <- c(0, max(out_retest$tests/1000, na.rm = TRUE) * 1.3) - col1 <- rgb(230, 180, 205, 200, max=255) - col2 <- rgb(150, 0, 50, 200, max=255) - col3 <- rgb(100, 0, 100, 200, max=255) + col1 <- grDevices::rgb(230, 180, 205, 200, max=255) + col2 <- grDevices::rgb(150, 0, 50, 200, max=255) + col3 <- grDevices::rgb(100, 0, 100, 200, max=255) if (plot_title == TRUE) { main_title <- "Distribution of HIV+ Tests by \n Awareness and ART Status" @@ -199,32 +199,32 @@ plot_retest_test_pos <- function(mod, fp, likdat, cnt, relative = F, main = main_title, xlim = c(2000, yr_pred), ylab = expression(paste(N^o, " of HIV+ Tests (in 1,000)")), xlab = 'Year') - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c(out_retest$newuna/1000, rev(out_retest$tests/1000)), col = col1, border = NA) - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c(out_retest$newdiag/1000, rev(out_retest$newuna/1000)), col = col2, border = NA) - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c(rep(0, length = length(out_retest$tests)), rev(out_retest$newdiag/1000)), col = col3, border = NA) if (plot_legend) { - legend('topleft', inset = 0.02, legend = c('Retests - PLHIV on ART','Retests - Aware Not on ART','New Diagnoses'), + graphics::legend('topleft', inset = 0.02, legend = c('Retests - PLHIV on ART','Retests - Aware Not on ART','New Diagnoses'), col=c(col1, col2, col3), bty = 'n', pch = c(15,15, 15), pt.cex = 2, cex = 0.9) } } else if(relative == T) { plot((out_retest$tests/out_retest$tests)*100 ~ out_retest$year, type='n', ylim = c(0,100), main = "Distribution of HIV+ Tests by \n Awareness and ART Status", xlim = c(2000, yr_pred), ylab = "Proportion of Total HIV+ Tests (%)", xlab = 'Year', lty=1, lwd=1) - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c((out_retest$newuna/out_retest$tests)*100, rev((out_retest$tests/out_retest$tests)*100)), col = col1, border = NA) - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c((out_retest$newdiag/out_retest$tests)*100, rev((out_retest$newuna/out_retest$tests)*100)), col = col2, border = NA) - polygon(x = c(out_retest$year, rev(out_retest$year)), + graphics::polygon(x = c(out_retest$year, rev(out_retest$year)), y = c(rep(0, length = length(out_retest$tests)), rev((out_retest$newdiag/out_retest$tests)*100)), col = col3, border = NA) - legend('bottomleft', inset = 0.02, legend = c('% Retests - PLHIV on ART','% Retests - Aware Not on ART','% New Diagnoses'), + graphics::legend('bottomleft', inset = 0.02, legend = c('% Retests - PLHIV on ART','% Retests - Aware Not on ART','% New Diagnoses'), col = c(col1, col2, col3), bg = 'grey90', box.col = 'grey90', pch = c(15,15, 15), pt.cex = 2, cex = 0.8) } else { print('TRUE for relative scale; FALSE for absolute scale') @@ -277,9 +277,9 @@ plot_prv_pos_yld <- function(mod, fp, likdat, cnt, yr_pred = 2022, ylim <- c(0, ifelse(max(c(prv, yld, ndx) * 1.3, na.rm = TRUE) > 1, 1, max(c(prv, yld, ndx) * 1.3, na.rm = TRUE))) * 100 - col1 <- rgb(255, 155, 50, 250, max = 255) - col2 <- rgb(255, 0, 130, 250, max = 255) - col3 <- rgb(100, 0, 100, 250, max = 255) + col1 <- grDevices::rgb(255, 155, 50, 250, max = 255) + col2 <- grDevices::rgb(255, 0, 130, 250, max = 255) + col3 <- grDevices::rgb(100, 0, 100, 250, max = 255) if (plot_title == TRUE) { main_title <- "HIV Prevalence (15+), Positivity, and \n Yield of New HIV Diagnoses" @@ -291,12 +291,10 @@ plot_prv_pos_yld <- function(mod, fp, likdat, cnt, yr_pred = 2022, main = main_title, xlim = c(2000, yr_pred), col = col1, lty = 2, lwd = 1.5, ylab = "Proportion (%)", xlab = 'Year') - lines(I(yld * 100) ~ out_test$year, col = col2, lwd = 1.5) - lines(I(ndx * 100) ~ out_postest$year, col = col3, lwd = 1.5) + graphics::lines(I(yld * 100) ~ out_test$year, col = col2, lwd = 1.5) + graphics::lines(I(ndx * 100) ~ out_postest$year, col = col3, lwd = 1.5) if (plot_legend) { - legend('topright', inset = 0.01, legend = c('HIV Prevalence (15+)','Positivity','Yield of New Diagnoses'), + graphics::legend('topright', inset = 0.01, legend = c('HIV Prevalence (15+)','Positivity','Yield of New Diagnoses'), col=c(col1, col2, col3), bty = 'n', lty = c(2, 1, 1), lwd = 1.5, cex = 0.7) } } - - diff --git a/R/simul.R b/R/simul.R index ae50735..1563111 100644 --- a/R/simul.R +++ b/R/simul.R @@ -19,8 +19,8 @@ simul.sample <- function(hessian, par, fp, sim = 3000, likdat = NULL, SIR = FALS ev <- eS$values if (!all(ev >= -1e-06 * abs(ev[1L]))) { vcova <- Matrix::nearPD(vcova, corr = FALSE)$mat - vcova <- matrix(vcova@x, - nrow = vcova@Dim[1], + vcova <- matrix(vcova@x, + nrow = vcova@Dim[1], ncol = vcova@Dim[2] ) } @@ -33,7 +33,10 @@ simul.sample <- function(hessian, par, fp, sim = 3000, likdat = NULL, SIR = FALS # Sampling Importance Resampling # http://www.sumsar.net/blog/2013/12/shaping_up_laplace_approximation/ if (SIR == TRUE) { - if (is.null(likdat)) { print('SIR needs to include '); break } + if (is.null(likdat)) { + print('SIR needs to include ') + return() + } # Sample parameters from multivariate t-distribution to get thicker tails par_sir <- mvtnorm::rmvt(n = nsir, delta = par, sigma = as.matrix(vcova), df = 2) # We had the mode to the resampled values (to be sure that it is included in the CI) @@ -99,7 +102,7 @@ simul.run <- function(samp, fp, progress = NULL){ nbtestpos_ss <- add_ss_indices(out_nbtest_pos, fp$ss) nbtest_ss$hivstatus <- as.character(nbtest_ss$hivstatus) nbtestpos_ss$hivstatus <- as.character(nbtestpos_ss$hivstatus) - + # Create parameters (proper scale, etc.), and simulate model for (i in 1:nrow(samp)){ fp <- create_hts_param(samp[i,], fp) @@ -123,7 +126,7 @@ simul.run <- function(samp, fp, progress = NULL){ # Simulation of the dataset #' @export simul.test <- function(opt, fp, sim = 3000, likdat = NULL, - SIR = FALSE, nsir = 50000, with_replacement = TRUE, + SIR = FALSE, nsir = 50000, with_replacement = TRUE, sir_progress = NULL, run_progress = NULL){ samp <- simul.sample(opt$hessian, opt$par, fp, sim = sim, likdat = likdat, @@ -131,13 +134,13 @@ simul.test <- function(opt, fp, sim = 3000, likdat = NULL, simul.run(samp, fp, progress = run_progress) } -# Get the confidence interval for the data, first for parameters, second for +# Get the confidence interval for the data, first for parameters, second for # the data itself getCI <- function(df) { n_col <- ncol(df) # Function for the confidence interval of the estimates CI <- apply(X = df[,(6:n_col)], MARGIN = 1, - FUN = quantile, probs = c(0.025, 0.975), na.rm = TRUE) + FUN = stats::quantile, probs = c(0.025, 0.975), na.rm = TRUE) lower <- CI[1,] upper <- CI[2,] df1 <- data.frame(df[c(1:5)], lower, upper) diff --git a/R/time_functions.R b/R/time_functions.R index fc27a96..27a2977 100644 --- a/R/time_functions.R +++ b/R/time_functions.R @@ -54,9 +54,9 @@ prb_dx_one_yr <- function(fp, year = c(2000:2022), age = "15-24", sex = "male", ind_yr <- year[n] - fp$ss$proj_start + 1 # what is the average age of incident HIV infection - avg_age15 <- 14 + weighted.mean(1:10, w = fp$infections[1:10, ind_sex, ind_yr]) - avg_age25 <- 14 + weighted.mean(11:20, w = fp$infections[11:20, ind_sex, ind_yr]) - avg_age35 <- 14 + weighted.mean(21:35, w = fp$infections[21:35, ind_sex, ind_yr]) + avg_age15 <- 14 + stats::weighted.mean(1:10, w = fp$infections[1:10, ind_sex, ind_yr]) + avg_age25 <- 14 + stats::weighted.mean(11:20, w = fp$infections[11:20, ind_sex, ind_yr]) + avg_age35 <- 14 + stats::weighted.mean(21:35, w = fp$infections[21:35, ind_sex, ind_yr]) # Testing rates by CD4 category test_rate15 <- fp$diagn_rate[, 1, ind_sex, ind_th, ind_yr] @@ -186,7 +186,7 @@ pool_prb_dx_one_yr <- function(mod, fp, year = c(2000:2022), denom$w_e[denom$agegr == time_mat$age[i] & denom$sex == time_mat$sex[i]] } } time_mat_sel <- time_mat[time_mat$age %in% age & time_mat$sex %in% sex, ] - val$prb1yr[val$year == year[n]] <- weighted.mean(time_mat_sel$prb1, w = time_mat_sel$w) + val$prb1yr[val$year == year[n]] <- stats::weighted.mean(time_mat_sel$prb1, w = time_mat_sel$w) } return(val) } @@ -217,7 +217,7 @@ simul_pool_prb_dx_one_yr <- function(samp, mod, fp, year = c(2010:2022), prb1yr = NA, prb1yr_lci = NA, prb1yr_uci = NA) for (n in 1:n_year) { - uncertainty <- quantile(val_df$prb1yr[val_df$year == year[n]], probs = c(0.5, 0.025, 0.975)) + uncertainty <- stats::quantile(val_df$prb1yr[val_df$year == year[n]], probs = c(0.5, 0.025, 0.975)) val$prb1yr[val$year == year[n]] <- uncertainty[1] val$prb1yr_lci[val$year == year[n]] <- uncertainty[2] val$prb1yr_uci[val$year == year[n]] <- uncertainty[3] diff --git a/README.Rmd b/README.Rmd index 8a0e0ba..70a77d6 100644 --- a/README.Rmd +++ b/README.Rmd @@ -18,7 +18,7 @@ library(data.table) ``` # first90 -[![Travis-CI Build Status](https://travis-ci.com/mrc-ide/first90release.svg?branch=master)](https://travis-ci.com/mrc-ide/first90release) +![R CMD check](https://github.com/mrc-ide/first90release/actions/workflows/R-CMD-check.yml/badge.svg) UNAIDS put forward the ambitious 90-90-90 target to end the AIDS epidemic by 2030. This target aims for 90% of people living with HIV (PLHIV) to be aware of their HIV-positive status, 90% of those diagnosed to receive antiretroviral therapy, and 90% of those on treatment to have a suppressed viral load by 2020 (each reaching 95% by 2030). HIV testing remains an important bottleneck in this cascade, however, and obtaining reliable epidemiological data on the proportion of PLHIV aware of their status is difficult. Such information is nevertheless crucial to effectively monitor HIV prevention efforts. Tracking progress towards achievement of this “first 90” target could be improved by combining population-based surveys and programmatic data on the number of HIV tests performed (and yield) in a coherent deterministic/statistical model. This type of integrative systems modelling is especially useful to fully consider HIV incidence, mortality, testing behaviours, as well as to coherently combine different sources of data. @@ -146,18 +146,18 @@ first90::plot_out_evertest_mbyage(mod, fp, likdat, cnt, survey_hts, out_evertest We can compare HIV tests' positivity through time, the estimated *true* yield of new HIV diagnoses, and compare those to population-level HIV prevalence. ```{r, fig.height = 5, fig.width = 5, warning = FALSE} -par(mfrow = c(1,1)) +graphics::par(mfrow = c(1,1)) first90::plot_prv_pos_yld(mod, fp, likdat, cnt, yr_pred = 2018) ``` We can also examine some ouptuts related to the distribution of HIV tests performed in those susceptibles to HIV and PLHIV by different awareness and treatment status. (First sets of graphs is on the absolute scale, second one on the relative scale.) ```{r, fig.height = 5, fig.width = 8.5, warning = FALSE} -par(mfrow = c(1,2)) +graphics::par(mfrow = c(1,2)) first90::plot_retest_test_neg(mod, fp, likdat, cnt) first90::plot_retest_test_pos(mod, fp, likdat, cnt) -par(mfrow = c(1,2)) +graphics::par(mfrow = c(1,2)) first90::plot_retest_test_neg(mod, fp, likdat, cnt, relative = TRUE) first90::plot_retest_test_pos(mod, fp, likdat, cnt, relative = TRUE) ``` diff --git a/README.md b/README.md index 50a4f3d..92eb9f5 100644 --- a/README.md +++ b/README.md @@ -4,8 +4,7 @@ Please edit that file and run `knitr::knit("README.Rmd") from the root of this d # first90 -[![Travis-CI Build -Status](https://travis-ci.com/mrc-ide/first90release.svg?branch=master)](https://travis-ci.com/mrc-ide/first90release) +![R CMD check](https://github.com/mrc-ide/first90release/actions/workflows/R-CMD-check.yml/badge.svg) UNAIDS put forward the ambitious 90-90-90 target to end the AIDS epidemic by 2030. This target aims for 90% of people living with HIV diff --git a/data-raw/initial-theta0-values.R b/data-raw/initial-theta0-values.R index f7fe1c3..6490a11 100644 --- a/data-raw/initial-theta0-values.R +++ b/data-raw/initial-theta0-values.R @@ -28,13 +28,13 @@ theta0 <- c(par_med_f15to24rate, rep(par_med_f15to24rate[n_k_2021], n_k - n_k_20 # old theta - keeping for comparisons theta1 <- c(log(seq(0.01, 0.2, length.out = n_k)), # Females 15-24 baseline testing rate 2000 to 2020 - qlogis(rep(1.5 / 8, n_k - 10)), # Factor testing among diagnosed, untreated from 2010 to 2020 - qlogis(rep(0.6 / 1.1, 2)), # RR for males in 2005, 2012 - qlogis(rep((1.93 - 0.95) / 7.05, 2)), # Factor increase among previously tested - qlogis(0.5), # RR re-test among PLHIV - qlogis(0.25), # Factor testing among diagnosed, on ART - rep(qlogis(0.9 / 5.9), 4), # Rate ratio for 25-34 - qlogis(0.5)) # Date of inflection point for logistic growth curve OI + stats::qlogis(rep(1.5 / 8, n_k - 10)), # Factor testing among diagnosed, untreated from 2010 to 2020 + stats::qlogis(rep(0.6 / 1.1, 2)), # RR for males in 2005, 2012 + stats::qlogis(rep((1.93 - 0.95) / 7.05, 2)), # Factor increase among previously tested + stats::qlogis(0.5), # RR re-test among PLHIV + stats::qlogis(0.25), # Factor testing among diagnosed, on ART + rep(stats::qlogis(0.9 / 5.9), 4), # Rate ratio for 25-34 + stats::qlogis(0.5)) # Date of inflection point for logistic growth curve OI stopifnot(length(theta0) == n_k + n_k-10 + 11) stopifnot(length(theta1) == n_k + n_k-10 + 11) diff --git a/first90release.Rproj b/first90release.Rproj new file mode 100644 index 0000000..497f8bf --- /dev/null +++ b/first90release.Rproj @@ -0,0 +1,20 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source diff --git a/man/get_pjn_country.Rd b/man/get_pjn_country.Rd index abcbd9f..48924a3 100644 --- a/man/get_pjn_country.Rd +++ b/man/get_pjn_country.Rd @@ -13,5 +13,5 @@ get_pjn_country(pjn) Get country name from parsed PJN } \details{ -\code{pjn} should be via \code{read.csv(pjn_file, as.is = TRUE)} +\code{pjn} should be via \code{first90_read_csv_character(pjn_file)} } diff --git a/man/get_pjn_region.Rd b/man/get_pjn_region.Rd index 28bd522..314048d 100644 --- a/man/get_pjn_region.Rd +++ b/man/get_pjn_region.Rd @@ -13,5 +13,5 @@ get_pjn_region(pjn) Get subnational region from parsed PJN } \details{ -\code{pjn} should be via \code{read.csv(pjn_file, as.is = TRUE)} +\code{pjn} should be via \code{first90_read_csv_character(pjn_file)} } diff --git a/man/prepare_inputs.Rd b/man/prepare_inputs.Rd index 62f984d..11fed31 100644 --- a/man/prepare_inputs.Rd +++ b/man/prepare_inputs.Rd @@ -21,7 +21,8 @@ The aggregation makes a number of assumptions: } \examples{ -pjnzlist <- list.files("~/Documents/Data/Spectrum files/2018 final/SSA/", "CotedIvoire.*PJNZ$", full.names=TRUE, ignore.case=TRUE) +pjnzlist <- list.files("~/Documents/Data/Spectrum files/2018 final/SSA/", + "CotedIvoire.*PJNZ$", full.names=TRUE, ignore.case=TRUE) pjnzlist <- "~/Documents/Data/Spectrum files/2018 final/SSA/Malawi_2018_version_8.PJNZ" }