diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 7760355..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 -sudo: false -cache: packages - -after_success: - - Rscript -e 'covr::codecov()' \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 3f4fa76..5b81bea 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: eppasm Title: Age-structured EPP Model for HIV Epidemic Estimates -Version: 0.7.4 +Version: 0.7.5 Authors@R: person("Jeff", "Eaton", email = "jeffrey.eaton@imperial.ac.uk", role = c("aut", "cre")) Description: What the package does (one paragraph). Depends: R (>= 3.1.0), @@ -11,6 +11,8 @@ Imports: epp (>= 0.4.1), fastmatch (>= 1.1), mvtnorm, + plyr, + readxl, reshape2, vroom, xml2 diff --git a/NEWS.md b/NEWS.md index b7bf6ce..643b334 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +## eppasm 0.7.5 + +* Qualify all package names and add all required packages into Imports section. + ## eppasm 0.7.4 * Implement recovery to next higher CD4 category following ART interruption for those on ART greater than one year. diff --git a/R/aggregate-fits.R b/R/aggregate-fits.R index 8d639a6..40a9d84 100644 --- a/R/aggregate-fits.R +++ b/R/aggregate-fits.R @@ -22,7 +22,7 @@ create_aggr_input <- function(inputlist){ ancrtcens <- do.call(rbind, lapply(eppdlist, "[[", "ancrtcens")) if(!is.null(ancrtcens) && nrow(ancrtcens)){ ancrtcens$x <- ancrtcens$prev * ancrtcens$n - ancrtcens <- aggregate(cbind(x,n) ~ year, ancrtcens, sum) + ancrtcens <- stats::aggregate(cbind(x,n) ~ year, ancrtcens, sum) ancrtcens$prev <- ancrtcens$x / ancrtcens$n ancrtcens <- ancrtcens[c("year", "prev", "n")] } diff --git a/R/cluster-functions.R b/R/cluster-functions.R index d3ef832..e21c69d 100644 --- a/R/cluster-functions.R +++ b/R/cluster-functions.R @@ -1,5 +1,5 @@ -pick_maxlpd <- function(set){set <- set[!sapply(set, is.null)]; set[[which.max(sapply(set, function(x) tail(x$stat,1)[1]))]]} +pick_maxlpd <- function(set){set <- set[!sapply(set, is.null)]; set[[which.max(sapply(set, function(x) utils::tail(x$stat,1)[1]))]]} get_imisiter <- function(x) nrow(x$stat) -get_logmargpost <- function(x) tail(x$stat, 1)[1] -get_nunique <- function(x) tail(x$stat, 1)[2] +get_logmargpost <- function(x) utils::tail(x$stat, 1)[1] +get_nunique <- function(x) utils::tail(x$stat, 1)[2] diff --git a/R/fit-model.R b/R/fit-model.R index 7582c53..84b8f36 100644 --- a/R/fit-model.R +++ b/R/fit-model.R @@ -48,7 +48,7 @@ prepare_spec_fit <- function(pjnz, proj.end=2016.5, popadjust = NULL, popupdate= ## output - val <- setNames(vector("list", length(eppd)), names(eppd)) + val <- stats::setNames(vector("list", length(eppd)), names(eppd)) set.list.attr <- function(obj, attrib, value.lst) mapply(function(set, value){ attributes(set)[[attrib]] <- value; set}, obj, value.lst) @@ -109,7 +109,7 @@ tidy_hhs_data <- function(eppd){ hhs } -create_subpop_specfp <- function(projp, demp, eppd, epp_t0=setNames(rep(1975, length(eppd)), names(eppd)), ..., popadjust=TRUE, popupdate=TRUE, perc_urban=NULL){ +create_subpop_specfp <- function(projp, demp, eppd, epp_t0=stats::setNames(rep(1975, length(eppd)), names(eppd)), ..., popadjust=TRUE, popupdate=TRUE, perc_urban=NULL){ country <- attr(eppd, "country") country_code <- attr(eppd, "country_code") @@ -195,7 +195,7 @@ create_subpop_specfp <- function(projp, demp, eppd, epp_t0=setNames(rep(1975, le subpop.dist <- prop.table(sapply(demp.subpop, get15to49pop, 2010)) if(nrow(subset(eppd[[1]]$hhs, used)) != 0){ # HH survey data available - hhsprev.means <- sapply(lapply(eppd, function(dat) na.omit(dat$hhs$prev[dat$hhs$used])), mean) + hhsprev.means <- sapply(lapply(eppd, function(dat) stats::na.omit(dat$hhs$prev[dat$hhs$used])), mean) art.dist <- prop.table(subpop.dist * hhsprev.means) } else { ## no HH survey data ## Apportion ART according to relative average ANC prevalence in each subpopulation @@ -244,7 +244,7 @@ prepare_national_fit <- function(pjnz, upd.path=NULL, proj.end=2013.5, hiv_steps epp.input <- epp::read_epp_input(pjnz) ## output - val <- setNames(vector("list", length(eppd)), names(eppd)) + val <- stats::setNames(vector("list", length(eppd)), names(eppd)) val <- list() ## aggregate census data across regions @@ -253,7 +253,7 @@ prepare_national_fit <- function(pjnz, upd.path=NULL, proj.end=2013.5, hiv_steps ancrtcens <- subset(ancrtcens, !is.na(prev) & !is.na(n)) if(nrow(ancrtcens)){ ancrtcens$x <- ancrtcens$prev * ancrtcens$n - ancrtcens <- aggregate(cbind(x,n) ~ year, ancrtcens, sum) + ancrtcens <- stats::aggregate(cbind(x,n) ~ year, ancrtcens, sum) ancrtcens$prev <- ancrtcens$x / ancrtcens$n ancrtcens <- ancrtcens[c("year", "prev", "n")] } @@ -284,9 +284,9 @@ fitmod <- function(obj, ..., epp=FALSE, B0 = 1e5, B = 1e4, B.re = 3000, number_k ## ... : updates to fixed parameters (fp) object to specify fitting options if(epp) - fp <- update(attr(obj, 'eppfp'), ...) + fp <- stats::update(attr(obj, 'eppfp'), ...) else - fp <- update(attr(obj, 'specfp'), ...) + fp <- stats::update(attr(obj, 'specfp'), ...) ## Prepare likelihood data @@ -347,13 +347,13 @@ fitmod <- function(obj, ..., epp=FALSE, B0 = 1e5, B = 1e4, B.re = 3000, number_k lpost0 <- likelihood(X0, fp, likdat, log=TRUE) + prior(X0, fp, log=TRUE) opt_init <- X0[which.max(lpost0)[1],] } - opt <- optim(opt_init, optfn, fp=fp, likdat=likdat, method=opt_method, control=list(fnscale=-1, trace=4, maxit=opt_maxit, ndeps=rep(opt_diffstep, length(opt_init)))) + opt <- stats::optim(opt_init, optfn, fp=fp, likdat=likdat, method=opt_method, control=list(fnscale=-1, trace=4, maxit=opt_maxit, ndeps=rep(opt_diffstep, length(opt_init)))) opt$fp <- fp opt$likdat <- likdat opt$param <- fnCreateParam(opt$par, fp) - opt$mod <- simmod(update(fp, list=opt$param)) + opt$mod <- simmod(stats::update(fp, list=opt$param)) if(opthess){ - opt$hessian <- optimHess(opt_init, optfn, fp=fp, likdat=likdat, + opt$hessian <- stats::optimHess(opt_init, optfn, fp=fp, likdat=likdat, control=list(fnscale=-1, trace=4, maxit=1e3, @@ -450,7 +450,7 @@ simfit.specfit <- function(fit, if(rwproj) fit <- rw_projection(fit) - fp_list <- lapply(fit$param, function(par) update(fit$fp, list=par)) + fp_list <- lapply(fit$param, function(par) stats::update(fit$fp, list=par)) mod.list <- lapply(fp_list, simmod) } else { @@ -598,7 +598,7 @@ simfit.eppfit <- function(fit, rwproj=fit$fp$eppmod == "rspline", pregprev=TRUE) fit$param <- lapply(fit$param, function(par){par$rvec <- sim_rvec_rwproj(par$rvec, firstidx, lastidx, fit$fp$dt); par}) } - fp.list <- lapply(fit$param, function(par) update(fit$fp, list=par)) + fp.list <- lapply(fit$param, function(par) stats::update(fit$fp, list=par)) mod.list <- lapply(fp.list, simmod) fit$rvec <- sapply(mod.list, attr, "rvec") @@ -636,7 +636,7 @@ sim_mod_list <- function(fit, rwproj=fit$fp$eppmod == "rspline"){ fit$param <- lapply(fit$param, function(par){par$rvec <- epp:::sim_rvec_rwproj(par$rvec, firstidx, lastidx, dt); par}) } - fp.list <- lapply(fit$param, function(par) update(fit$fp, list=par)) + fp.list <- lapply(fit$param, function(par) stats::update(fit$fp, list=par)) mod.list <- lapply(fp.list, simmod) ## strip unneeded attributes to preserve memory diff --git a/R/imis.R b/R/imis.R index 760d9ea..7bd8dd7 100644 --- a/R/imis.R +++ b/R/imis.R @@ -26,7 +26,7 @@ imis <- function(B0, B, B_re, number_k, opt_k=NULL, fp, likdat, ## Draw initial samples from prior distribution X_k <- sample_prior(B0, fp) # Draw initial samples from the prior distribution - cov_prior = cov(X_k) # estimate of the prior covariance + cov_prior = stats::cov(X_k) # estimate of the prior covariance ## Locations and covariance of mixture components center_all <- list() @@ -83,7 +83,7 @@ imis <- function(B0, B, B_re, number_k, opt_k=NULL, fp, likdat, max(weights), # the maximum weight 1/sum(weights^2), # the effictive sample size -sum(weights*log(weights), na.rm = TRUE) / log(length(weights)), # the entropy relative to uniform !!! NEEDS UDPATING FOR OMITTED SAMPLES !!! - var(weights/mean(weights)), # the variance of scaled weights !!! NEEDS UDPATING FOR OMITTED SAMPLES !!! + stats::var(weights/mean(weights)), # the variance of scaled weights !!! NEEDS UDPATING FOR OMITTED SAMPLES !!! as.numeric(iter_stop_time - iter_start_time)[3]) iter_start_time <- iter_stop_time @@ -108,25 +108,25 @@ imis <- function(B0, B, B_re, number_k, opt_k=NULL, fp, likdat, nlposterior <- function(theta){-prior(theta, fp, log=TRUE)-likelihood(theta, fp, likdat, log=TRUE)} ## opt <- optimization_step(theta_init, nlposterior, cov_prior) # Version by Bao uses prior covariance to parscale optimizer - opt <- optimization_step(theta_init, nlposterior, cov(X_all[1:n_all, ])) + opt <- optimization_step(theta_init, nlposterior, stats::cov(X_all[1:n_all, ])) center_all[[k]] <- opt$mu sigma_all[[k]] <- opt$sigma ## exclude the neighborhood of the local optima - distance_remain <- mahalanobis(X_all[seq_len(n_all),], center_all[[k]], diag(diag(sigma_all[[k]]))) + distance_remain <- stats::mahalanobis(X_all[seq_len(n_all),], center_all[[k]], diag(diag(sigma_all[[k]]))) idx_exclude <- union(idx_exclude, order(distance_remain)[seq_len(n_all/length(opt_k))]) } else { ## choose mixture component centered at input with current maximum weight center_all[[k]] <- X_all[which.max(weights),] - distance_all <- mahalanobis(X_all[1:n_all,], center_all[[k]], diag(diag(cov_prior))) # Raftery & Bao version - ## distance_all <- mahalanobis(X_all[1:n_all,], center_all[[k]], cov(X_all[1:n_all, ])) # Suggested by Fasiolo et al. + distance_all <- stats::mahalanobis(X_all[1:n_all,], center_all[[k]], diag(diag(cov_prior))) # Raftery & Bao version + ## distance_all <- stats::mahalanobis(X_all[1:n_all,], center_all[[k]], stats::cov(X_all[1:n_all, ])) # Suggested by Fasiolo et al. which_close <- order(distance_all)[seq_len(min(n_all, B))] # Choose B nearest inputs (use n_all if n_all < B) - sigma_all[[k]] <- cov.wt(X_all[which_close, , drop=FALSE], wt = weights[which_close]+1/n_all, center = center_all[[k]])$cov # Raftery & Bao version + sigma_all[[k]] <- stats::cov.wt(X_all[which_close, , drop=FALSE], wt = weights[which_close]+1/n_all, center = center_all[[k]])$cov # Raftery & Bao version if(any(is.na(sigma_all[[k]]))) sigma_all[[k]] <- cov_prior - ## sigma_all[[k]] <- cov.wt(X_all[which_close,], center = center_all[[k]])$cov # Suggested by Fasiolo et al. + ## sigma_all[[k]] <- stats::cov.wt(X_all[which_close,], center = center_all[[k]])$cov # Suggested by Fasiolo et al. } ## Update mixture weights according with new mixture component @@ -152,14 +152,14 @@ optimization_step <- function(theta, fn, cov){ ## The rough optimizer uses the Nelder-Mead algorithm. ptm.opt = proc.time() - optNM <- optim(theta, fn, method="Nelder-Mead", + optNM <- stats::optim(theta, fn, method="Nelder-Mead", control=list(maxit=5000, parscale=sqrt(diag(cov)))) ## The more efficient optimizer uses the BFGS algorithm optBFGS <- try(stop(""), TRUE) step_expon <- 0.2 while(inherits(optBFGS, "try-error") && step_expon <= 0.8){ - optBFGS <- try(optim(optNM$par, fn, method="BFGS", hessian=TRUE, + optBFGS <- try(stats::optim(optNM$par, fn, method="BFGS", hessian=TRUE, control=list(parscale=sqrt(diag(cov)), ndeps=rep(.Machine$double.eps^step_expon, length(optNM$par)), maxit=1000)), silent=TRUE) step_expon <- step_expon + 0.05 diff --git a/R/infections.R b/R/infections.R index 94c7816..0e714dd 100644 --- a/R/infections.R +++ b/R/infections.R @@ -10,17 +10,17 @@ calc_infections_eppspectrum <- function(fp, pop, hivpop, artpop, i, ii, r_ts){ hivn.ii <- sum(pop[p.age15to49.idx,,hivn.idx,i]) hivn.ii <- hivn.ii - sum(pop[p.age15to49.idx[1],,hivn.idx,i])*(1-DT*(ii-1)) - hivn.ii <- hivn.ii + sum(pop[tail(p.age15to49.idx,1)+1,,hivn.idx,i])*(1-DT*(ii-1)) + hivn.ii <- hivn.ii + sum(pop[utils::tail(p.age15to49.idx,1)+1,,hivn.idx,i])*(1-DT*(ii-1)) hivp.ii <- sum(pop[p.age15to49.idx,,hivp.idx,i]) hivp.ii <- hivp.ii - sum(pop[p.age15to49.idx[1],,hivp.idx,i])*(1-DT*(ii-1)) - hivp.ii <- hivp.ii + sum(pop[tail(p.age15to49.idx,1)+1,,hivp.idx,i])*(1-DT*(ii-1)) + hivp.ii <- hivp.ii + sum(pop[utils::tail(p.age15to49.idx,1)+1,,hivp.idx,i])*(1-DT*(ii-1)) art.ii <- sum(artpop[,,h.age15to49.idx,,i]) if(sum(hivpop[,h.age15to49.idx[1],,i]) + sum(artpop[,,h.age15to49.idx[1],,i]) > 0) art.ii <- art.ii - sum(pop[p.age15to49.idx[1],,hivp.idx,i] * colSums(artpop[,,h.age15to49.idx[1],,i],,2) / (colSums(hivpop[,h.age15to49.idx[1],,i],,1) + colSums(artpop[,,h.age15to49.idx[1],,i],,2))) * (1-DT*(ii-1)) - if(sum(hivpop[,tail(h.age15to49.idx, 1)+1,,i]) + sum(artpop[,,tail(h.age15to49.idx, 1)+1,,i]) > 0) - art.ii <- art.ii + sum(pop[tail(p.age15to49.idx,1)+1,,hivp.idx,i] * colSums(artpop[,,tail(h.age15to49.idx, 1)+1,,i],,2) / (colSums(hivpop[,tail(h.age15to49.idx, 1)+1,,i],,1) + colSums(artpop[,,tail(h.age15to49.idx, 1)+1,,i],,2))) * (1-DT*(ii-1)) + if(sum(hivpop[,utils::tail(h.age15to49.idx, 1)+1,,i]) + sum(artpop[,,utils::tail(h.age15to49.idx, 1)+1,,i]) > 0) + art.ii <- art.ii + sum(pop[utils::tail(p.age15to49.idx,1)+1,,hivp.idx,i] * colSums(artpop[,,utils::tail(h.age15to49.idx, 1)+1,,i],,2) / (colSums(hivpop[,utils::tail(h.age15to49.idx, 1)+1,,i],,1) + colSums(artpop[,,utils::tail(h.age15to49.idx, 1)+1,,i],,2))) * (1-DT*(ii-1)) transm_prev <- (hivp.ii - art.ii + fp$relinfectART*art.ii) / (hivn.ii+hivp.ii) @@ -79,4 +79,3 @@ create_beers <- function(n5yr){ ## Beers coefficient matrix beers_Amat <- create_beers(17)[16:81, 4:17] - diff --git a/R/likelihood.R b/R/likelihood.R index 97c3f26..c190961 100644 --- a/R/likelihood.R +++ b/R/likelihood.R @@ -75,27 +75,27 @@ ancrtsite.beta.pr.sd <- 1.0 #' pjnz <- system.file("extdata/testpjnz", "Botswana2017.PJNZ", package="eppasm") #' bw <- prepare_spec_fit(pjnz, proj.end=2021.5) #' -#' +#' #' bw_u_ancsite <- attr(bw$Urban, "eppd")$ancsitedat #' fp <- attr(bw$Urban, "specfp") #' #' ancsite_pred_df(bw_u_ancsite, fp) -#' +#' ancsite_pred_df <- function(ancsite_df, fp) { df <- ancsite_df anchor.year <- fp$ss$proj_start - + df$aidx <- df$age - fp$ss$AGE_START + 1L df$yidx <- df$year - anchor.year + 1 - + ## List of all unique agegroup / year combinations for which prevalence is needed datgrp <- unique(df[c("aidx", "yidx", "agspan")]) datgrp$qMidx <- seq_len(nrow(datgrp)) ## Indices for accessing prevalence offset from datgrp df <- merge(df, datgrp) - + list(df = df, datgrp = datgrp) } @@ -111,15 +111,15 @@ prepare_ancsite_likdat <- function(ancsitedat, fp){ df <- d$df[c("site", "year", "used", "type", "age", "agspan", "n", "prev", "aidx", "yidx", "qMidx")] - + ## Calculate probit transformed prevalence and variance approximation df$pstar <- (df$prev * df$n + 0.5) / (df$n + 1) - df$W <- qnorm(df$pstar) + df$W <- stats::qnorm(df$pstar) df$v <- 2 * pi * exp(df$W^2) * df$pstar * (1 - df$pstar) / df$n ## Design matrix for fixed effects portion df$type <- factor(df$type, c("ancss", "ancrt")) - Xancsite <- model.matrix(~type, df) + Xancsite <- stats::model.matrix(~type, df) ## Indices for observation df_idx.lst <- split(seq_len(nrow(df)), factor(df$site)) @@ -130,7 +130,7 @@ prepare_ancsite_likdat <- function(ancsitedat, fp){ df_idx.lst = df_idx.lst) } -#' @export +#' @export ll_ancsite <- function(mod, fp, coef=c(0, 0), vinfl=0, dat){ df <- dat$df @@ -144,24 +144,24 @@ ll_ancsite <- function(mod, fp, coef=c(0, 0), vinfl=0, dat){ ## approximate mid-year prevalence. ## NOTE: This will be inefficient if values for every year because it duplicates ## prevalence calculation for the same year. - + if (fp$projection_period == "calendar") { pregprevM_last <- agepregprev(mod, fp, dat$datgrp$aidx, dat$datgrp$yidx-1L, dat$datgrp$agspan) pregprevM <- 0.5 * (pregprevM + pregprevM_last) } - - qM <- suppressWarnings(qnorm(pregprevM)) + + qM <- suppressWarnings(stats::qnorm(pregprevM)) if(any(is.na(qM)) || any(qM == -Inf) || any(qM > 2)) ## prev < 0.977 return(-Inf) - + mu <- qM[df$qMidx] + dat$Xancsite %*% coef d <- df$W - mu v <- df$v + vinfl - + d.lst <- lapply(dat$df_idx.lst, function(idx) d[idx]) v.lst <- lapply(dat$df_idx.lst, function(idx) v[idx]) - + log(anclik::anc_resid_lik(d.lst, v.lst)) } @@ -182,9 +182,9 @@ ancrtcens.vinfl.pr.rate <- 1/0.015 prepare_ancrtcens_likdat <- function(dat, fp){ anchor.year <- fp$ss$proj_start - + x.ancrt <- (dat$prev*dat$n+0.5)/(dat$n+1) - dat$W.ancrt <- qnorm(x.ancrt) + dat$W.ancrt <- stats::qnorm(x.ancrt) dat$v.ancrt <- 2*pi*exp(dat$W.ancrt^2)*x.ancrt*(1-x.ancrt)/dat$n if(!exists("age", dat)) @@ -195,33 +195,33 @@ prepare_ancrtcens_likdat <- function(dat, fp){ dat$aidx <- dat$age - fp$ss$AGE_START + 1 dat$yidx <- dat$year - anchor.year + 1 - + return(dat) } -#' @export +#' @export ll_ancrtcens <- function(mod, dat, fp, pointwise = FALSE){ if(!nrow(dat)) return(0) qM.prev <- agepregprev(mod, fp, dat$aidx, dat$yidx, dat$agspan) - + ## If calendar year projection, average current and previous year prevalence to ## approximate mid-year prevalence. ## NOTE: This will be inefficient if values for every year because it duplicates ## prevalence calculation for the same year. - + if (fp$projection_period == "calendar") { qM.prev_last <- qM.prev <- agepregprev(mod, fp, dat$aidx, dat$yidx-1L, dat$agspan) qM.prev <- 0.5 * (qM.prev + qM.prev_last) } - - qM.prev <- suppressWarnings(qnorm(qM.prev)) + + qM.prev <- suppressWarnings(stats::qnorm(qM.prev)) if(any(is.na(qM.prev))) val <- rep(-Inf, nrow(dat)) else - val <- dnorm(dat$W.ancrt, qM.prev, sqrt(dat$v.ancrt + fp$ancrtcens.vinfl), log=TRUE) + val <- stats::dnorm(dat$W.ancrt, qM.prev, sqrt(dat$v.ancrt + fp$ancrtcens.vinfl), log=TRUE) if(pointwise) return(val) @@ -242,7 +242,7 @@ fnCreateParam <- function(theta, fp){ for(i in seq_along(fp$prior_args)) assign(names(fp$prior_args)[i], fp$prior_args[[i]]) } - + if(fp$eppmod %in% c("rspline", "logrw")){ epp_nparam <- fp$numKnots+1 @@ -262,7 +262,7 @@ fnCreateParam <- function(theta, fp){ param <- list(beta = beta, rvec = as.vector(fp$rvec.spldes %*% beta)) - + if(fp$eppmod %in% "logrw") param$rvec <- exp(param$rvec) @@ -290,7 +290,7 @@ fnCreateParam <- function(theta, fp){ param$rvec <- create_rvec(theta[1:fp$rt$n_param], fp$rt) param$iota <- transf_iota(theta[fp$rt$n_param+1], fp) } - + if(fp$ancsitedata){ param$ancbias <- theta[epp_nparam+1] if(!exists("v.infl", where=fp)){ @@ -301,14 +301,14 @@ fnCreateParam <- function(theta, fp){ } else anclik_nparam <- 0 - - + + paramcurr <- epp_nparam+anclik_nparam if(exists("ancrt", fp) && fp$ancrt %in% c("census", "both")){ param$log_frr_adjust <- theta[paramcurr+1] param$frr_cd4 <- fp$frr_cd4 * exp(param$log_frr_adjust) param$frr_art <- fp$frr_art * exp(param$log_frr_adjust) - + if(!exists("ancrtcens.vinfl", fp)){ param$ancrtcens.vinfl <- exp(theta[paramcurr+2]) paramcurr <- paramcurr+2 @@ -319,7 +319,7 @@ fnCreateParam <- function(theta, fp){ param$ancrtsite.beta <- theta[paramcurr+1] paramcurr <- paramcurr+1 } - + return(param) } @@ -334,7 +334,7 @@ fnCreateParam <- function(theta, fp){ prepare_hhsageprev_likdat <- function(hhsage, fp){ anchor.year <- floor(min(fp$proj.steps)) - hhsage$W.hhs <- qnorm(hhsage$prev) + hhsage$W.hhs <- stats::qnorm(hhsage$prev) hhsage$v.hhs <- 2*pi*exp(hhsage$W.hhs^2)*hhsage$se^2 hhsage$sd.W.hhs <- sqrt(hhsage$v.hhs) @@ -354,7 +354,7 @@ prepare_hhsageprev_likdat <- function(hhsage, fp){ startage <- as.integer(sub("([0-9]*)-([0-9]*)", "\\1", hhsage$agegr)) endage <- as.integer(sub("([0-9]*)-([0-9]*)", "\\2", hhsage$agegr)) - + hhsage$sidx <- match(hhsage$sex, c("both", "male", "female")) - 1L hhsage$aidx <- startage - fp$ss$AGE_START+1L hhsage$yidx <- as.integer(hhsage$year - (anchor.year - 1)) @@ -364,7 +364,7 @@ prepare_hhsageprev_likdat <- function(hhsage, fp){ } #' Log likelihood for age-specific household survey prevalence -#' @export +#' @export ll_hhsage <- function(mod, fp, dat, pointwise = FALSE){ prevM.age <- ageprev(mod, aidx = dat$aidx, sidx = dat$sidx, yidx = dat$yidx, agspan = dat$agspan) @@ -374,13 +374,13 @@ ll_hhsage <- function(mod, fp, dat, pointwise = FALSE){ prevM.age_last <- ageprev(mod, aidx = dat$aidx, sidx = dat$sidx, yidx = dat$yidx-1L, agspan = dat$agspan) prevM.age <- 0.5 * (prevM.age + prevM.age_last) } - - qM.age <- suppressWarnings(qnorm(prevM.age)) - + + qM.age <- suppressWarnings(stats::qnorm(prevM.age)) + if(any(is.na(qM.age))) val <- rep(-Inf, nrow(dat)) else - val <- dnorm(dat$W.hhs, qM.age, dat$sd.W.hhs, log=TRUE) + val <- stats::dnorm(dat$W.hhs, qM.age, dat$sd.W.hhs, log=TRUE) if(pointwise) return(val) @@ -389,7 +389,7 @@ ll_hhsage <- function(mod, fp, dat, pointwise = FALSE){ #' Log likelihood for age-specific household survey prevalence using binomial approximation -#' @export +#' @export ll_hhsage_binom <- function(mod, fp, dat, pointwise = FALSE){ prevM.age <- suppressWarnings(ageprev(mod, aidx = dat$aidx, sidx = dat$sidx, yidx = dat$yidx, agspan = dat$agspan)) @@ -417,10 +417,10 @@ ll_hhsage_binom <- function(mod, fp, dat, pointwise = FALSE){ #' ## Household survey ART coverage likelihood prepare_hhsartcov_likdat <- function(hhsartcov, fp){ - + anchor.year <- floor(min(fp$proj.steps)) - hhsartcov$W.hhs <- qnorm(hhsartcov$artcov) + hhsartcov$W.hhs <- stats::qnorm(hhsartcov$artcov) hhsartcov$v.hhs <- 2*pi*exp(hhsartcov$W.hhs^2)*hhsartcov$se^2 hhsartcov$sd.W.hhs <- sqrt(hhsartcov$v.hhs) @@ -432,7 +432,7 @@ prepare_hhsartcov_likdat <- function(hhsartcov, fp){ startage <- as.integer(sub("([0-9]*)-([0-9]*)", "\\1", hhsartcov$agegr)) endage <- as.integer(sub("([0-9]*)-([0-9]*)", "\\2", hhsartcov$agegr)) - + hhsartcov$sidx <- match(hhsartcov$sex, c("both", "male", "female")) - 1L hhsartcov$aidx <- startage - fp$ss$AGE_START+1L hhsartcov$yidx <- as.integer(hhsartcov$year - (anchor.year - 1)) @@ -453,13 +453,13 @@ ll_hhsartcov <- function(mod, fp, dat, pointwise = FALSE){ if (fp$projection_period == "calendar") { artcovM_obs <- 0.5 * (artcovM_obs + artcovM[dat$yidx-1L]) } - - qM <- qnorm(artcovM_obs) - + + qM <- stats::qnorm(artcovM_obs) + if(any(is.na(qM))) val <- rep(-Inf, nrow(dat)) else - val <- dnorm(dat$W.hhs, qM, dat$sd.W.hhs, log=TRUE) + val <- stats::dnorm(dat$W.hhs, qM, dat$sd.W.hhs, log=TRUE) if(pointwise) return(val) @@ -494,7 +494,7 @@ prepare_hhsincid_likdat <- function(hhsincid, fp){ #' @param hhsincid.dat prepared houshold survey incidence estimates (see perp ll_hhsincid <- function(mod, fp, hhsincid.dat){ logincid <- log(incid(mod, fp)) - ll.incid <- sum(dnorm(hhsincid.dat$log_incid, logincid[hhsincid.dat$idx], hhsincid.dat$log_incid.se, TRUE)) + ll.incid <- sum(stats::dnorm(hhsincid.dat$log_incid, logincid[hhsincid.dat$idx], hhsincid.dat$log_incid.se, TRUE)) return(ll.incid) } @@ -503,23 +503,23 @@ ll_hhsincid <- function(mod, fp, hhsincid.dat){ #### Likelihood function #### ############################### -#' @export +#' @export prepare_likdat <- function(eppd, fp){ likdat <- list() - + likdat$hhs.dat <- prepare_hhsageprev_likdat(eppd$hhs, fp) if(exists("ancsitedat", where=eppd)){ ancsitedat <- eppd$ancsitedat - + if(exists("ancrt", fp) && fp$ancrt %in% c("none", "census")) ancsitedat <- subset(ancsitedat, type == "ancss") - + likdat$ancsite.dat <- prepare_ancsite_likdat(ancsitedat, fp) } - + if(exists("ancrtcens", where=eppd)){ if(exists("ancrt", fp) && fp$ancrt %in% c("none", "site")) eppd$ancrtcens <- NULL @@ -538,7 +538,7 @@ prepare_likdat <- function(eppd, fp){ -#' @export +#' @export lprior <- function(theta, fp){ if(exists("prior_args", where = fp)){ @@ -557,34 +557,34 @@ lprior <- function(theta, fp){ lpr <- bayes_lmvt(theta[(1+fp$rtpenord):nk], tau2_prior_shape, tau2_prior_rate) if(exists("r0logiotaratio", fp) && fp$r0logiotaratio) - lpr <- lpr + dunif(theta[nk+1], r0logiotaratio.unif.prior[1], r0logiotaratio.unif.prior[2], log=TRUE) + lpr <- lpr + stats::dunif(theta[nk+1], r0logiotaratio.unif.prior[1], r0logiotaratio.unif.prior[2], log=TRUE) else lpr <- lpr + lprior_iota(theta[nk+1], fp) } else if(fp$eppmod == "rlogistic") { epp_nparam <- 5 - lpr <- sum(dnorm(theta[1:4], rlog_pr_mean, rlog_pr_sd, log=TRUE)) + lpr <- sum(stats::dnorm(theta[1:4], rlog_pr_mean, rlog_pr_sd, log=TRUE)) lpr <- lpr + lprior_iota(theta[5], fp) } else if(fp$eppmod == "rtrend"){ # rtrend - + epp_nparam <- 7 - - lpr <- dunif(round(theta[1]), t0.unif.prior[1], t0.unif.prior[2], log=TRUE) + - dnorm(round(theta[2]), t1.pr.mean, t1.pr.sd, log=TRUE) + - dnorm(theta[3], logr0.pr.mean, logr0.pr.sd, log=TRUE) + - sum(dnorm(theta[4:7], rtrend.beta.pr.mean, rtrend.beta.pr.sd, log=TRUE)) + + lpr <- stats::dunif(round(theta[1]), t0.unif.prior[1], t0.unif.prior[2], log=TRUE) + + stats::dnorm(round(theta[2]), t1.pr.mean, t1.pr.sd, log=TRUE) + + stats::dnorm(theta[3], logr0.pr.mean, logr0.pr.sd, log=TRUE) + + sum(stats::dnorm(theta[4:7], rtrend.beta.pr.mean, rtrend.beta.pr.sd, log=TRUE)) } else if(fp$eppmod == "rhybrid"){ epp_nparam <- fp$rt$n_param+1 - lpr <- sum(dnorm(theta[1:4], rlog_pr_mean, rlog_pr_sd, log=TRUE)) + - sum(dnorm(theta[4+1:fp$rt$n_rw], 0, rw_prior_sd, log=TRUE)) + lpr <- sum(stats::dnorm(theta[1:4], rlog_pr_mean, rlog_pr_sd, log=TRUE)) + + sum(stats::dnorm(theta[4+1:fp$rt$n_rw], 0, rw_prior_sd, log=TRUE)) lpr <- lpr + lprior_iota(theta[fp$rt$n_param+1], fp) } if(fp$ancsitedata){ - lpr <- lpr + dnorm(theta[epp_nparam+1], ancbias.pr.mean, ancbias.pr.sd, log=TRUE) + lpr <- lpr + stats::dnorm(theta[epp_nparam+1], ancbias.pr.mean, ancbias.pr.sd, log=TRUE) if(!exists("v.infl", where=fp)){ anclik_nparam <- 2 - lpr <- lpr + dexp(exp(theta[epp_nparam+2]), vinfl.prior.rate, TRUE) + theta[epp_nparam+2] # additional ANC variance + lpr <- lpr + stats::dexp(exp(theta[epp_nparam+2]), vinfl.prior.rate, TRUE) + theta[epp_nparam+2] # additional ANC variance } else anclik_nparam <- 1 } else @@ -592,15 +592,15 @@ lprior <- function(theta, fp){ paramcurr <- epp_nparam+anclik_nparam if(exists("ancrt", fp) && fp$ancrt %in% c("census", "both")){ - lpr <- lpr + dnorm(theta[paramcurr+1], log_frr_adjust.pr.mean, log_frr_adjust.pr.sd, log=TRUE) + lpr <- lpr + stats::dnorm(theta[paramcurr+1], log_frr_adjust.pr.mean, log_frr_adjust.pr.sd, log=TRUE) if(!exists("ancrtcens.vinfl", fp)){ - lpr <- lpr + dexp(exp(theta[paramcurr+2]), ancrtcens.vinfl.pr.rate, TRUE) + theta[paramcurr+2] + lpr <- lpr + stats::dexp(exp(theta[paramcurr+2]), ancrtcens.vinfl.pr.rate, TRUE) + theta[paramcurr+2] paramcurr <- paramcurr+2 } else paramcurr <- paramcurr+1 } if(exists("ancrt", fp) && fp$ancrt %in% c("site", "both")){ - lpr <- lpr + dnorm(theta[paramcurr+1], ancrtsite.beta.pr.mean, ancrtsite.beta.pr.sd, log=TRUE) + lpr <- lpr + stats::dnorm(theta[paramcurr+1], ancrtsite.beta.pr.mean, ancrtsite.beta.pr.sd, log=TRUE) paramcurr <- paramcurr+1 } @@ -608,13 +608,13 @@ lprior <- function(theta, fp){ } -#' @export +#' @export ll <- function(theta, fp, likdat){ theta.last <<- theta - fp <- update(fp, list=fnCreateParam(theta, fp)) + fp <- stats::update(fp, list=fnCreateParam(theta, fp)) if (fp$eppmod == "rspline") - if (any(is.na(fp$rvec)) || min(fp$rvec) < 0 || max(fp$rvec) > 20) + if (any(is.na(fp$rvec)) || min(fp$rvec) < 0 || max(fp$rvec) > 20) return(-Inf) mod <- simmod(fp) @@ -650,7 +650,7 @@ ll <- function(theta, fp, likdat){ else ll.hhsartcov <- 0 - + if(exists("equil.rprior", where=fp) && fp$equil.rprior){ if(fp$eppmod != "rspline") stop("error in ll(): equil.rprior is only for use with r-spline model") @@ -659,20 +659,20 @@ ll <- function(theta, fp, likdat){ likdat$hhs.dat$yidx, likdat$ancrtcens.dat$yidx, likdat$hhsincid.dat$idx) - - qM.all <- suppressWarnings(qnorm(prev(mod))) + + qM.all <- suppressWarnings(stats::qnorm(prev(mod))) if(any(is.na(qM.all[lastdata.idx - 9:0]))) { ll.rprior <- -Inf } else { rvec.ann <- fp$rvec[fp$proj.steps %% 1 == 0.5] - equil.rprior.mean <- epp:::muSS/(1-pnorm(qM.all[lastdata.idx])) + equil.rprior.mean <- epp:::muSS/(1-stats::pnorm(qM.all[lastdata.idx])) if(!is.null(fp$prior_args$equil.rprior.sd)) equil.rprior.sd <- fp$prior_args$equil.rprior.sd else - equil.rprior.sd <- sqrt(mean((epp:::muSS/(1-pnorm(qM.all[lastdata.idx - 9:0])) - rvec.ann[lastdata.idx - 9:0])^2)) # empirical sd based on 10 previous years - - ll.rprior <- sum(dnorm(rvec.ann[(lastdata.idx+1L):length(qM.all)], equil.rprior.mean, equil.rprior.sd, log=TRUE)) # prior starts year after last data + equil.rprior.sd <- sqrt(mean((epp:::muSS/(1-stats::pnorm(qM.all[lastdata.idx - 9:0])) - rvec.ann[lastdata.idx - 9:0])^2)) # empirical sd based on 10 previous years + + ll.rprior <- sum(stats::dnorm(rvec.ann[(lastdata.idx+1L):length(qM.all)], equil.rprior.mean, equil.rprior.sd, log=TRUE)) # prior starts year after last data } } else ll.rprior <- 0 @@ -728,7 +728,7 @@ sample.prior <- function(n, fp){ ancrt_nparam <- ancrt_nparam+1 nparam <- epp_nparam+anclik_nparam+ancrt_nparam - + ## Create matrix for storing samples mat <- matrix(NA, n, nparam) @@ -736,9 +736,9 @@ sample.prior <- function(n, fp){ epp_nparam <- fp$numKnots+1 if(fp$eppmod == "rspline") - mat[,1] <- rnorm(n, 1.5, 1) # u[1] + mat[,1] <- stats::rnorm(n, 1.5, 1) # u[1] else # logrw - mat[,1] <- rnorm(n, 0.2, 1) # u[1] + mat[,1] <- stats::rnorm(n, 0.2, 1) # u[1] if(fp$eppmod == "logrw"){ mat[,2:fp$rt$n_rw] <- bayes_rmvt(n, fp$rt$n_rw-1, rw_prior_shape, rw_prior_rate) # u[2:numKnots] } else { @@ -746,43 +746,43 @@ sample.prior <- function(n, fp){ } if(exists("r0logiotaratio", fp) && fp$r0logiotaratio) - mat[,fp$numKnots+1] <- runif(n, r0logiotaratio.unif.prior[1], r0logiotaratio.unif.prior[2]) # ratio r0 / log(iota) + mat[,fp$numKnots+1] <- stats::runif(n, r0logiotaratio.unif.prior[1], r0logiotaratio.unif.prior[2]) # ratio r0 / log(iota) else mat[,fp$numKnots+1] <- sample_iota(n, fp) } else if(fp$eppmod == "rlogistic"){ - mat[,1:4] <- t(matrix(rnorm(4*n, rlog_pr_mean, rlog_pr_sd), 4)) + mat[,1:4] <- t(matrix(stats::rnorm(4*n, rlog_pr_mean, rlog_pr_sd), 4)) mat[,5] <- sample_iota(n, fp) } else if(fp$eppmod == "rtrend"){ # r-trend - mat[,1] <- runif(n, t0.unif.prior[1], t0.unif.prior[2]) # t0 - mat[,2] <- rnorm(n, t1.pr.mean, t1.pr.sd) - mat[,3] <- rnorm(n, logr0.pr.mean, logr0.pr.sd) # r0 - mat[,4:7] <- t(matrix(rnorm(4*n, rtrend.beta.pr.mean, rtrend.beta.pr.sd), 4, n)) # beta + mat[,1] <- stats::runif(n, t0.unif.prior[1], t0.unif.prior[2]) # t0 + mat[,2] <- stats::rnorm(n, t1.pr.mean, t1.pr.sd) + mat[,3] <- stats::rnorm(n, logr0.pr.mean, logr0.pr.sd) # r0 + mat[,4:7] <- t(matrix(stats::rnorm(4*n, rtrend.beta.pr.mean, rtrend.beta.pr.sd), 4, n)) # beta } else if(fp$eppmod == "rhybrid") { - mat[,1:4] <- t(matrix(rnorm(4*n, rlog_pr_mean, rlog_pr_sd), 4)) - mat[,4+1:fp$rt$n_rw] <- rnorm(n*fp$rt$n_rw, 0, rw_prior_sd) # u[2:numKnots] + mat[,1:4] <- t(matrix(stats::rnorm(4*n, rlog_pr_mean, rlog_pr_sd), 4)) + mat[,4+1:fp$rt$n_rw] <- stats::rnorm(n*fp$rt$n_rw, 0, rw_prior_sd) # u[2:numKnots] mat[,fp$rt$n_param+1] <- sample_iota(n, fp) } ## sample ANC bias paramters if(fp$ancsitedata){ - mat[,epp_nparam+1] <- rnorm(n, ancbias.pr.mean, ancbias.pr.sd) # ancbias parameter + mat[,epp_nparam+1] <- stats::rnorm(n, ancbias.pr.mean, ancbias.pr.sd) # ancbias parameter if(!exists("v.infl", where=fp)) - mat[,epp_nparam+2] <- log(rexp(n, vinfl.prior.rate)) + mat[,epp_nparam+2] <- log(stats::rexp(n, vinfl.prior.rate)) } ## sample ANCRT parameters paramcurr <- epp_nparam+anclik_nparam if(exists("ancrt", where=fp) && fp$ancrt %in% c("census", "both")){ - mat[,paramcurr+1] <- rnorm(n, log_frr_adjust.pr.mean, log_frr_adjust.pr.sd) + mat[,paramcurr+1] <- stats::rnorm(n, log_frr_adjust.pr.mean, log_frr_adjust.pr.sd) if(!exists("ancrtcens.vinfl", fp)){ - mat[,paramcurr+2] <- log(rexp(n, ancrtcens.vinfl.pr.rate)) + mat[,paramcurr+2] <- log(stats::rexp(n, ancrtcens.vinfl.pr.rate)) paramcurr <- paramcurr+2 } else paramcurr <- paramcurr+1 } if(exists("ancrt", where=fp) && fp$ancrt %in% c("site", "both")){ - mat[,paramcurr+1] <- rnorm(n, ancrtsite.beta.pr.mean, ancrtsite.beta.pr.sd) + mat[,paramcurr+1] <- stats::rnorm(n, ancrtsite.beta.pr.mean, ancrtsite.beta.pr.sd) paramcurr <- paramcurr+1 } @@ -802,9 +802,9 @@ ldsamp <- function(theta, fp){ nk <- fp$numKnots if(fp$eppmod == "rspline") # u[1] - lpr <- dnorm(theta[1], 1.5, 1, log=TRUE) + lpr <- stats::dnorm(theta[1], 1.5, 1, log=TRUE) else # logrw - lpr <- dnorm(theta[1], 0.2, 1, log=TRUE) + lpr <- stats::dnorm(theta[1], 0.2, 1, log=TRUE) if(fp$eppmod == "logrw") bayes_lmvt(theta[2:fp$rt$n_rw], rw_prior_shape, rw_prior_rate) @@ -813,34 +813,34 @@ ldsamp <- function(theta, fp){ if(exists("r0logiotaratio", fp) && fp$r0logiotaratio) - lpr <- lpr + dunif(theta[nk+1], r0logiotaratio.unif.prior[1], r0logiotaratio.unif.prior[2], log=TRUE) + lpr <- lpr + stats::dunif(theta[nk+1], r0logiotaratio.unif.prior[1], r0logiotaratio.unif.prior[2], log=TRUE) else lpr <- lpr + ldsamp_iota(theta[nk+1], fp) } else if(fp$eppmod == "rlogistic") { epp_nparam <- 5 - lpr <- sum(dnorm(theta[1:4], rlog_pr_mean, rlog_pr_sd, log=TRUE)) + lpr <- sum(stats::dnorm(theta[1:4], rlog_pr_mean, rlog_pr_sd, log=TRUE)) lpr <- lpr + ldsamp_iota(theta[5], fp) } else if(fp$eppmod == "rtrend"){ # rtrend epp_nparam <- 7 - lpr <- dunif(round(theta[1]), t0.unif.prior[1], t0.unif.prior[2], log=TRUE) + - dnorm(round(theta[2]), t1.pr.mean, t1.pr.sd, log=TRUE) + - dnorm(theta[3], logr0.pr.mean, logr0.pr.sd, log=TRUE) + - sum(dnorm(theta[4:7], rtrend.beta.pr.mean, rtrend.beta.pr.sd, log=TRUE)) + lpr <- stats::dunif(round(theta[1]), t0.unif.prior[1], t0.unif.prior[2], log=TRUE) + + stats::dnorm(round(theta[2]), t1.pr.mean, t1.pr.sd, log=TRUE) + + stats::dnorm(theta[3], logr0.pr.mean, logr0.pr.sd, log=TRUE) + + sum(stats::dnorm(theta[4:7], rtrend.beta.pr.mean, rtrend.beta.pr.sd, log=TRUE)) } else if(fp$eppmod == "rhybrid"){ epp_nparam <- fp$rt$n_param+1 - lpr <- sum(dnorm(theta[1:4], rlog_pr_mean, rlog_pr_sd, log=TRUE)) + - sum(dnorm(theta[4+1:fp$rt$n_rw], 0, rw_prior_sd, log=TRUE)) + lpr <- sum(stats::dnorm(theta[1:4], rlog_pr_mean, rlog_pr_sd, log=TRUE)) + + sum(stats::dnorm(theta[4+1:fp$rt$n_rw], 0, rw_prior_sd, log=TRUE)) lpr <- lpr + ldsamp_iota(theta[fp$rt$n_param+1], fp) } if(fp$ancsitedata){ - lpr <- lpr + dnorm(theta[epp_nparam+1], ancbias.pr.mean, ancbias.pr.sd, log=TRUE) + lpr <- lpr + stats::dnorm(theta[epp_nparam+1], ancbias.pr.mean, ancbias.pr.sd, log=TRUE) if(!exists("v.infl", where=fp)){ anclik_nparam <- 2 - lpr <- lpr + dexp(exp(theta[epp_nparam+2]), vinfl.prior.rate, TRUE) + theta[epp_nparam+2] # additional ANC variance + lpr <- lpr + stats::dexp(exp(theta[epp_nparam+2]), vinfl.prior.rate, TRUE) + theta[epp_nparam+2] # additional ANC variance } else anclik_nparam <- 1 } else @@ -848,15 +848,15 @@ ldsamp <- function(theta, fp){ paramcurr <- epp_nparam+anclik_nparam if(exists("ancrt", fp) && fp$ancrt %in% c("census", "both")){ - lpr <- lpr + dnorm(theta[paramcurr+1], log_frr_adjust.pr.mean, log_frr_adjust.pr.sd, log=TRUE) + lpr <- lpr + stats::dnorm(theta[paramcurr+1], log_frr_adjust.pr.mean, log_frr_adjust.pr.sd, log=TRUE) if(!exists("ancrtcens.vinfl", fp)){ - lpr <- lpr + dexp(exp(theta[paramcurr+2]), ancrtcens.vinfl.pr.rate, TRUE) + theta[paramcurr+2] + lpr <- lpr + stats::dexp(exp(theta[paramcurr+2]), ancrtcens.vinfl.pr.rate, TRUE) + theta[paramcurr+2] paramcurr <- paramcurr+2 } else paramcurr <- paramcurr+1 } if(exists("ancrt", fp) && fp$ancrt %in% c("site", "both")){ - lpr <- lpr + dnorm(theta[paramcurr+1], ancrtsite.beta.pr.mean, ancrtsite.beta.pr.sd, log=TRUE) ## + + lpr <- lpr + stats::dnorm(theta[paramcurr+1], ancrtsite.beta.pr.mean, ancrtsite.beta.pr.sd, log=TRUE) ## + paramcurr <- paramcurr+1 } diff --git a/R/likelihood_ancsite.R b/R/likelihood_ancsite.R index c744e44..7dc9d8a 100644 --- a/R/likelihood_ancsite.R +++ b/R/likelihood_ancsite.R @@ -1,18 +1,18 @@ #' Sample from posterior ditribution for ANC site level random effects -#' +#' sample_b_site <- function(mod, fp, dat, resid=TRUE){ coef <- c(fp$ancbias, fp$ancrtsite.beta) vinfl <- fp$v.infl df <- dat$df - qM <- suppressWarnings(qnorm(agepregprev(mod, fp, dat$datgrp$aidx, dat$datgrp$yidx, dat$datgrp$agspan))) + qM <- suppressWarnings(stats::qnorm(agepregprev(mod, fp, dat$datgrp$aidx, dat$datgrp$yidx, dat$datgrp$agspan))) mu <- qM[df$qMidx] + dat$Xancsite %*% coef d <- df$W - mu v <- df$v + vinfl - + d.lst <- lapply(dat$df_idx.lst, function(idx) d[idx]) v.lst <- lapply(dat$df_idx.lst, function(idx) v[idx]) @@ -29,38 +29,38 @@ sample_b_site <- function(mod, fp, dat, resid=TRUE){ #' Sample posterior predictions for site-level ANC observations -#' +#' sample_ancsite_pred <- function(mod, fp, newdata, b_site){ coef <- c(fp$ancbias, fp$ancrtsite.beta) - + df <- newdata$df - qM <- suppressWarnings(qnorm(agepregprev(mod, fp, newdata$datgrp$aidx, newdata$datgrp$yidx, newdata$datgrp$agspan))) + qM <- suppressWarnings(stats::qnorm(agepregprev(mod, fp, newdata$datgrp$aidx, newdata$datgrp$yidx, newdata$datgrp$agspan))) ## Design matrix for fixed effects portion df$type <- factor(df$type, c("ancss", "ancrt")) - Xancsite <- model.matrix(~type, df) + Xancsite <- stats::model.matrix(~type, df) mu <- qM[df$qMidx] + Xancsite %*% coef + b_site[match(df$site, names(b_site))] - v <- 2 * pi * exp(mu^2) * pnorm(mu) * (1 - pnorm(mu))/df$n + fp$v.infl - rnorm(nrow(df), mu, sqrt(v)) -} + v <- 2 * pi * exp(mu^2) * stats::pnorm(mu) * (1 - stats::pnorm(mu))/df$n + fp$v.infl + stats::rnorm(nrow(df), mu, sqrt(v)) +} #' Pointwise likelihood for site-level ANC observations given site-level effects -#' +#' ll_ancsite_conditional <- function(mod, fp, newdata, b_site){ coef <- c(fp$ancbias, fp$ancrtsite.beta) - + df <- newdata$df - qM <- suppressWarnings(qnorm(agepregprev(mod, fp, newdata$datgrp$aidx, newdata$datgrp$yidx, newdata$datgrp$agspan))) + qM <- suppressWarnings(stats::qnorm(agepregprev(mod, fp, newdata$datgrp$aidx, newdata$datgrp$yidx, newdata$datgrp$agspan))) ## Design matrix for fixed effects portion df$type <- factor(df$type, c("ancss", "ancrt")) - Xancsite <- model.matrix(~type, df) + Xancsite <- stats::model.matrix(~type, df) mu <- qM[df$qMidx] + Xancsite %*% coef + b_site[match(df$site, names(b_site))] - v <- 2 * pi * exp(mu^2) * pnorm(mu) * (1 - pnorm(mu))/df$n + fp$v.infl - dnorm(df$W, mu, sqrt(v), log=TRUE) -} + v <- 2 * pi * exp(mu^2) * stats::pnorm(mu) * (1 - stats::pnorm(mu))/df$n + fp$v.infl + stats::dnorm(df$W, mu, sqrt(v), log=TRUE) +} diff --git a/R/model-outputs.R b/R/model-outputs.R index 08df696..3defdb4 100644 --- a/R/model-outputs.R +++ b/R/model-outputs.R @@ -2,8 +2,8 @@ estci2 <- function(x, na.rm=TRUE){ if(is.vector(x)) x <- matrix(x, 1) nd <- length(dim(x)) val <- apply(x, seq_len(nd-1), function(y) c(mean(y, na.rm=na.rm), - sd(y, na.rm=na.rm), - quantile(y, c(0.5, 0.025, 0.975), na.rm=na.rm))) + stats::sd(y, na.rm=na.rm), + stats::quantile(y, c(0.5, 0.025, 0.975), na.rm=na.rm))) val <- aperm(val, c(2:nd, 1)) dimnames(val)[[nd]] <- c("mean", "se", "median", "lower", "upper") val @@ -107,7 +107,7 @@ summary_outputs <- function(fit, modlist){ create_outputs <- function(fit){ paramlist <- lapply(seq_len(nrow(fit$resample)), function(ii) fnCreateParam(fit$resample[ii,], fit$fp)) - fplist <- lapply(paramlist, function(par) update(fit$fp, list=par)) + fplist <- lapply(paramlist, function(par) stats::update(fit$fp, list=par)) modlist <- lapply(fplist, simmod) out <- summary_outputs(fit, modlist) diff --git a/R/outpred.R b/R/outpred.R index 09158a7..0d0ac5f 100644 --- a/R/outpred.R +++ b/R/outpred.R @@ -12,25 +12,25 @@ outpred_prev <- function(fit, out, subsample=NULL){ ## Impute missing standard errors based on model prevalence na_i <- which(is.na(dat$sd.W.hhs)) na_p <- rowMeans(fit$ageprevdat)[na_i] - dat$sd.W.hhs[na_i] <- sqrt(na_p * (1 - na_p) / dat$n_eff[na_i]) / dnorm(qnorm(na_p)) + dat$sd.W.hhs[na_i] <- sqrt(na_p * (1 - na_p) / dat$n_eff[na_i]) / stats::dnorm(stats::qnorm(na_p)) M <- fit$ageprevdat - qM <- qnorm(M) - qpred <- array(rnorm(length(qM), qM, dat$sd.W.hhs), dim(qM)) + qM <- stats::qnorm(M) + qpred <- array(stats::rnorm(length(qM), qM, dat$sd.W.hhs), dim(qM)) elpd <- log(rowMeans(exp(ldbinom(dat$x_eff, dat$n_eff, M)))) resid <- rowMeans(M) - dat$prev - rse <- apply(M, 1, sd) / rowMeans(M) + rse <- apply(M, 1, stats::sd) / rowMeans(M) mae <- rowMeans(abs(M - dat$prev)) rmse <- sqrt(rowMeans((M - dat$prev)^2)) - elpd_q <- log(rowMeans(dnorm(dat$W.hhs, qM, dat$sd.W.hhs))) - rse_q <- apply(qM, 1, sd) / rowMeans(qM) + elpd_q <- log(rowMeans(stats::dnorm(dat$W.hhs, qM, dat$sd.W.hhs))) + rse_q <- apply(qM, 1, stats::sd) / rowMeans(qM) resid_q <- rowMeans(qM) - dat$W.hhs mae_q <- rowMeans(abs(qM - dat$W.hhs)) rmse_q<- sqrt(rowMeans((qM - dat$W.hhs)^2)) - qq <- mapply(function(f, x) f(x), apply(qpred, 1, ecdf), dat$W.hhs) + qq <- mapply(function(f, x) f(x), apply(qpred, 1, stats::ecdf), dat$W.hhs) vars <- intersect(c("country", "eppregion", "survyear", "year", "sex", "agegr", "prev", "se"), names(dat)) data.frame(dat[vars], @@ -50,7 +50,7 @@ outpred_hhs <- function(fit, newdata, subsample=NULL){ ## simulate model projections param_list <- lapply(seq_len(nrow(fit$resample)), function(ii) fnCreateParam(fit$resample[ii,], fit$fp)) - fp_list <- lapply(param_list, function(par) update(fit$fp, list=par)) + fp_list <- lapply(param_list, function(par) stats::update(fit$fp, list=par)) mod_list <- lapply(fp_list, simmod) ## new data for prediction @@ -63,16 +63,16 @@ outpred_hhs <- function(fit, newdata, subsample=NULL){ yidx = dat$yidx, agspan = dat$agspan) M <- matrix(M, nrow(dat)) - qM <- qnorm(M) - qpred <- array(rnorm(length(qM), qM, dat$sd.W.hhs), dim(qM)) + qM <- stats::qnorm(M) + qpred <- array(stats::rnorm(length(qM), qM, dat$sd.W.hhs), dim(qM)) - dat$elpd_q <- log(rowMeans(dnorm(dat$W.hhs, qM, dat$sd.W.hhs))) - dat$se_q <- apply(qM, 1, sd) + dat$elpd_q <- log(rowMeans(stats::dnorm(dat$W.hhs, qM, dat$sd.W.hhs))) + dat$se_q <- apply(qM, 1, stats::sd) dat$resid_q <- rowMeans(qM) - dat$W.hhs dat$pred <- rowMeans(M) - dat$se_p <- apply(M, 1, sd) + dat$se_p <- apply(M, 1, stats::sd) dat$resid <- rowMeans(M) - dat$prev - dat$qq <- mapply(function(f, x) f(x), apply(qpred, 1, ecdf), dat$W.hhs) + dat$qq <- mapply(function(f, x) f(x), apply(qpred, 1, stats::ecdf), dat$W.hhs) dat } @@ -88,7 +88,7 @@ outpred_ancsite <- function(fit, testdat){ ## simulate model projections param_list <- lapply(seq_len(nrow(fit$resample)), function(ii) fnCreateParam(fit$resample[ii,], fit$fp)) - fp_list <- lapply(param_list, function(par) update(fit$fp, list=par)) + fp_list <- lapply(param_list, function(par) stats::update(fit$fp, list=par)) mod_list <- lapply(fp_list, simmod) ## Site-level ANC data @@ -109,7 +109,7 @@ outpred_ancsite <- function(fit, testdat){ testdat$elpd <- log(rowMeans(exp(ancsite_ll))) testdat$resid <- rowMeans(ancsite_pred) - newdata$df$W - testdat$qq <- mapply(function(f, x) f(x), apply(ancsite_pred, 1, ecdf), newdata$df$W) + testdat$qq <- mapply(function(f, x) f(x), apply(ancsite_pred, 1, stats::ecdf), newdata$df$W) testdat } @@ -128,7 +128,7 @@ outpred_incid <- function(fit, newdata, subsample=NULL){ ## simulate model projections param_list <- lapply(seq_len(nrow(fit$resample)), function(ii) fnCreateParam(fit$resample[ii,], fit$fp)) - fp_list <- lapply(param_list, function(par) update(fit$fp, list=par)) + fp_list <- lapply(param_list, function(par) stats::update(fit$fp, list=par)) mod_list <- lapply(fp_list, simmod) ## new data for prediction @@ -140,15 +140,15 @@ outpred_incid <- function(fit, newdata, subsample=NULL){ incid_ll <- vapply(mod_list, ll_hhsincid, hhsincid.dat = dat, numeric(nrow(dat))) incid_ll <- matrix(incid_ll, nrow(dat)) - lpred <- array(rnorm(length(lM), lM, dat$log_incid.se), dim(lM)) + lpred <- array(stats::rnorm(length(lM), lM, dat$log_incid.se), dim(lM)) dat$elpd <- log(rowMeans(exp(incid_ll))) - dat$se_l <- apply(lM, 1, sd) + dat$se_l <- apply(lM, 1, stats::sd) dat$resid_l <- rowMeans(lM) - dat$log_incid dat$pred <- rowMeans(exp(lM)) - dat$se_p <- apply(exp(lM), 1, sd) + dat$se_p <- apply(exp(lM), 1, stats::sd) dat$resid <- rowMeans(exp(lM)) - dat$incid - dat$qq <- mapply(function(f, x) f(x), apply(lpred, 1, ecdf), dat$log_incid) + dat$qq <- mapply(function(f, x) f(x), apply(lpred, 1, stats::ecdf), dat$log_incid) dat } diff --git a/R/plot.R b/R/plot.R index a1fdf7b..a521578 100644 --- a/R/plot.R +++ b/R/plot.R @@ -1,13 +1,13 @@ cred.region <- function(x, y, ...){ if(nrow(y) != 2) y <- t(y) - polygon(c(x, rev(x)), c(y[1,], rev(y[2,])), border=NA, ...) + graphics::polygon(c(x, rev(x)), c(y[1,], rev(y[2,])), border=NA, ...) } -transp <- function(col, alpha=0.5) adjustcolor(col, alpha) +transp <- function(col, alpha=0.5) grDevices::adjustcolor(col, alpha) estci <- function(x){ if(is.vector(x)) x <- matrix(x, 1) - val <- cbind(rowMeans(x), t(apply(x, 1, quantile, c(0.5, 0.025, 0.975)))); colnames(val) <- c("mean", "median", "lower", "upper"); + val <- cbind(rowMeans(x), t(apply(x, 1, stats::quantile, c(0.5, 0.025, 0.975)))); colnames(val) <- c("mean", "median", "lower", "upper"); val } @@ -31,7 +31,7 @@ plot_compare_ageprev <- function(fit, fit2=NULL, fit3=NULL, specres=NULL, ylim=N survprev3 <- split(survprev3, factor(survprev3$survyear)) } ## - par(mfrow=c(4,2), mar=c(2, 3, 2, 1), tcl=-0.25, mgp=c(2, 0.5, 0), las=1, cex=1) + graphics::par(mfrow=c(4,2), mar=c(2, 3, 2, 1), tcl=-0.25, mgp=c(2, 0.5, 0), las=1, cex=1) for(isurv in names(survprev)) for(isex in c("male", "female")){ sp <- subset(survprev[[isurv]], sex==isex & as.integer(agegr) %in% 3:11) @@ -47,31 +47,31 @@ plot_compare_ageprev <- function(fit, fit2=NULL, fit3=NULL, specres=NULL, ylim=N paste0(sp$country[1], " ", survprev[[isurv]]$survyear[1], ", ", isex) plot(xx+0.5, sp$prev, type="n", xlim=c(4, 12), ylim=ylim, xaxt="n", main=main, xlab="", ylab="") - axis(1, xx+0.5, sp$agegr) + graphics::axis(1, xx+0.5, sp$agegr) ## - rect(xx+0.05, sp$lower, xx+0.95, sp$upper, + graphics::rect(xx+0.05, sp$lower, xx+0.95, sp$upper, col=transp(col[1]), border=NA) - segments(xx+0.05, sp$mean, xx+0.95, col=col[1], lwd=2) + graphics::segments(xx+0.05, sp$mean, xx+0.95, col=col[1], lwd=2) ## if(!is.null(fit2)){ - rect(xx+0.05, sp2$lower, xx+0.95, sp2$upper, + graphics::rect(xx+0.05, sp2$lower, xx+0.95, sp2$upper, col=transp(col[2]), border=NA) - segments(xx+0.05, sp2$mean, xx+0.95, col=col[2], lwd=2) + graphics::segments(xx+0.05, sp2$mean, xx+0.95, col=col[2], lwd=2) } if(!is.null(fit3)){ - rect(xx+0.05, sp3$lower, xx+0.95, sp3$upper, + graphics::rect(xx+0.05, sp3$lower, xx+0.95, sp3$upper, col=transp(col[3]), border=NA) - segments(xx+0.05, sp3$mean, xx+0.95, col=col[3], lwd=2) + graphics::segments(xx+0.05, sp3$mean, xx+0.95, col=col[3], lwd=2) } ## if(!is.null(specres)){ csex <- sub("(\\b[a-z]{1})", "\\U\\1" , isex, perl=TRUE) stryear <- as.character(survprev[[isurv]]$year[1]) specres.prev <- tapply(specres$hivpop[as.character(15:54), csex, stryear], rep(3:10, each=5), sum) / tapply(specres$totpop[as.character(15:54), csex, stryear], rep(3:10, each=5), sum) - segments(4:11+0.1, specres.prev, 4:11+0.9, lty=3, col="grey10", lwd=2) + graphics::segments(4:11+0.1, specres.prev, 4:11+0.9, lty=3, col="grey10", lwd=2) } - points(xx+0.5, sp$prev, pch=19) - segments(x0=xx+0.5, y0=sp$ci_l, y1=sp$ci_u) + graphics::points(xx+0.5, sp$prev, pch=19) + graphics::segments(x0=xx+0.5, y0=sp$ci_l, y1=sp$ci_u) } #### return(invisible()) @@ -81,52 +81,52 @@ plot_compare_ageprev <- function(fit, fit2=NULL, fit3=NULL, specres=NULL, ylim=N plot_prev <- function(fit, ..., ylim=NULL, xlim=c(1980, with(fit$fp$ss, proj_start+PROJ_YEARS-1)), col="blue", main="", ylab="prevalence", plotancdata=FALSE, plotprevdat=TRUE){ if(is.null(ylim)) - ylim <- c(0, 1.1*max(apply(fit$prev, 1, quantile, 0.975))) + ylim <- c(0, 1.1*max(apply(fit$prev, 1, stats::quantile, 0.975))) xx <- fit$fp$ss$proj_start-1+1:fit$fp$ss$PROJ_YEARS plot(xx, rowMeans(fit$prev), type="n", ylim=ylim, xlim=xlim, ylab=ylab, xlab="", yaxt="n", xaxt="n", main=main) - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) dots <- list(...) if(plotancdata){ - with(fit$likdat$anclik.dat, mapply(function(idx, W) points(idx+1970-1, pnorm(W), col=adjustcolor("grey", 0.5), pch=15), anc.idx.lst, W.lst)) - with(fit$likdat$anclik.dat, mapply(function(idx, W) lines(idx+1970-1, pnorm(W), col=adjustcolor("grey", 0.5), pch=15), anc.idx.lst, W.lst)) + with(fit$likdat$anclik.dat, mapply(function(idx, W) graphics::points(idx+1970-1, stats::pnorm(W), col=grDevices::adjustcolor("grey", 0.5), pch=15), anc.idx.lst, W.lst)) + with(fit$likdat$anclik.dat, mapply(function(idx, W) graphics::lines(idx+1970-1, stats::pnorm(W), col=grDevices::adjustcolor("grey", 0.5), pch=15), anc.idx.lst, W.lst)) } - + for(ii in seq_along(dots)) - cred.region(xx, apply(dots[[ii]]$prev, 1, quantile, c(0.025, 0.975)), col=transp(col[1+ii], 0.3)) - cred.region(xx, apply(fit$prev, 1, quantile, c(0.025, 0.975)), col=transp(col[1], 0.3)) + cred.region(xx, apply(dots[[ii]]$prev, 1, stats::quantile, c(0.025, 0.975)), col=transp(col[1+ii], 0.3)) + cred.region(xx, apply(fit$prev, 1, stats::quantile, c(0.025, 0.975)), col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, rowMeans(dots[[ii]]$prev), col=col[1+ii], lwd=1.5) - lines(xx, rowMeans(fit$prev), col=col[1], lwd=1.5) + graphics::lines(xx, rowMeans(dots[[ii]]$prev), col=col[1+ii], lwd=1.5) + graphics::lines(xx, rowMeans(fit$prev), col=col[1], lwd=1.5) ## if(plotprevdat){ - points(fit$likdat$hhslik.dat$year, fit$likdat$hhslik.dat$prev, pch=20) - segments(fit$likdat$hhslik.dat$year, - y0=pnorm(fit$likdat$hhslik.dat$W.hhs - qnorm(0.975)*fit$likdat$hhslik.dat$sd.W.hhs), - y1=pnorm(fit$likdat$hhslik.dat$W.hhs + qnorm(0.975)*fit$likdat$hhslik.dat$sd.W.hhs)) + graphics::points(fit$likdat$hhslik.dat$year, fit$likdat$hhslik.dat$prev, pch=20) + graphics::segments(fit$likdat$hhslik.dat$year, + y0=stats::pnorm(fit$likdat$hhslik.dat$W.hhs - stats::qnorm(0.975)*fit$likdat$hhslik.dat$sd.W.hhs), + y1=stats::pnorm(fit$likdat$hhslik.dat$W.hhs + stats::qnorm(0.975)*fit$likdat$hhslik.dat$sd.W.hhs)) } } plot_incid <- function(fit, ..., ylim=NULL, xlim=c(1980, with(fit$fp$ss, proj_start+PROJ_YEARS-1)), col="blue", main="", ylab="incidence rate"){ if(is.null(ylim)) - ylim <- c(0, 1.1*max(apply(fit$incid, 1, quantile, 0.975))) + ylim <- c(0, 1.1*max(apply(fit$incid, 1, stats::quantile, 0.975))) xx <- fit$fp$ss$proj_start-1+1:fit$fp$ss$PROJ_YEARS plot(xx, rowMeans(fit$incid), type="n", ylim=ylim, xlim=xlim, ylab=ylab, xlab="", yaxt="n", xaxt="n", main=main) - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) dots <- list(...) for(ii in seq_along(dots)) - cred.region(xx, apply(dots[[ii]]$incid, 1, quantile, c(0.025, 0.975)), col=transp(col[1+ii], 0.3)) - cred.region(xx, apply(fit$incid, 1, quantile, c(0.025, 0.975)), col=transp(col[1], 0.3)) + cred.region(xx, apply(dots[[ii]]$incid, 1, stats::quantile, c(0.025, 0.975)), col=transp(col[1+ii], 0.3)) + cred.region(xx, apply(fit$incid, 1, stats::quantile, c(0.025, 0.975)), col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, rowMeans(dots[[ii]]$incid), col=col[1+ii],lwd=1.5) - lines(xx, rowMeans(fit$incid), col=col, lwd=1.5) + graphics::lines(xx, rowMeans(dots[[ii]]$incid), col=col[1+ii],lwd=1.5) + graphics::lines(xx, rowMeans(fit$incid), col=col, lwd=1.5) ## if(exists("hhsincid.dat", where=fit$likdat)) - with(fit$likdat$hhsincid.dat,{ points(year, incid, pch=20); - segments(year, y0=exp(log_incid-qnorm(0.975)*log_incid.se), y1=exp(log_incid+qnorm(0.975)*log_incid.se))}) + with(fit$likdat$hhsincid.dat,{ graphics::points(year, incid, pch=20); + graphics::segments(year, y0=exp(log_incid-stats::qnorm(0.975)*log_incid.se), y1=exp(log_incid+stats::qnorm(0.975)*log_incid.se))}) } plot_rvec <- function(fit, ..., ylim=NULL, xlim=c(1980, with(fit$fp$ss, proj_start+PROJ_YEARS-1)), col="blue"){ @@ -137,21 +137,21 @@ plot_rvec <- function(fit, ..., ylim=NULL, xlim=c(1980, with(fit$fp$ss, proj_sta for(ii in seq_along(dots)) dots[[ii]]$rvec <- sapply(seq_len(ncol(dots[[ii]]$rvec)), function(i) replace(dots[[ii]]$rvec[,i], dots[[ii]]$fp$proj.steps < dots[[ii]]$param[[i]]$tsEpidemicStart, NA)) if(is.null(ylim)) - ylim <- c(0, min(quantile(apply(fit$rvec, 1, quantile, 0.975, na.rm=TRUE), 0.975, na.rm=TRUE), 5)) + ylim <- c(0, min(stats::quantile(apply(fit$rvec, 1, stats::quantile, 0.975, na.rm=TRUE), 0.975, na.rm=TRUE), 5)) plot(xx, rowMeans(fit$rvec, na.rm=TRUE)[idx], type="n", ylim=ylim, xlim=xlim, ylab="r(t)", yaxt="n", xlab="") - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) for(ii in seq_along(dots)) - cred.region(xx, apply(dots[[ii]]$rvec[idx,], 1, quantile, c(0.025, 0.975), na.rm=TRUE), col=transp(col[1+ii], 0.3)) - cred.region(xx, apply(fit$rvec[idx,], 1, quantile, c(0.025, 0.975), na.rm=TRUE), col=transp(col[1], 0.3)) + cred.region(xx, apply(dots[[ii]]$rvec[idx,], 1, stats::quantile, c(0.025, 0.975), na.rm=TRUE), col=transp(col[1+ii], 0.3)) + cred.region(xx, apply(fit$rvec[idx,], 1, stats::quantile, c(0.025, 0.975), na.rm=TRUE), col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, rowMeans(dots[[ii]]$rvec[idx,], na.rm=TRUE)[idx], col=col[1+ii]) - lines(xx, rowMeans(fit$rvec[idx,], na.rm=TRUE)[idx], col=col[1]) + graphics::lines(xx, rowMeans(dots[[ii]]$rvec[idx,], na.rm=TRUE)[idx], col=col[1+ii]) + graphics::lines(xx, rowMeans(fit$rvec[idx,], na.rm=TRUE)[idx], col=col[1]) return(invisible()) } plot_ancprev <- function(fit, ..., data=fit$likdat, ylim=NULL, xlim=c(1980, with(fit$fp$ss, proj_start+PROJ_YEARS-1)), col="blue", main="", ylab="ANC prevalence"){ - + ## Calculate parameters for bias term param <- apply(fit$resample, 1, fnCreateParam, fit$fp) ## ancrtcens.bias <- sapply(param, "[[", "ancrtcens.bias") @@ -165,27 +165,27 @@ plot_ancprev <- function(fit, ..., data=fit$likdat, ylim=NULL, xlim=c(1980, with dots_ancrt.prev <- lapply(dots, "[[", "pregprev") if(is.null(ylim)) - ylim <- c(0, 1.1*max(apply(ancrt.prev, 1, quantile, 0.975))) + ylim <- c(0, 1.1*max(apply(ancrt.prev, 1, stats::quantile, 0.975))) xx <- fit$fp$ss$proj_start-1+1:fit$fp$ss$PROJ_YEARS plot(xx, rowMeans(fit$pregprev), type="n", ylim=ylim, xlim=xlim, ylab=ylab, xlab="", yaxt="n", xaxt="n", main=main) - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) for(ii in seq_along(dots)) - cred.region(xx, apply(dots_ancrt.prev[[ii]], 1, quantile, c(0.025, 0.975)), col=transp(col[1+ii], 0.3)) - cred.region(xx, apply(ancrt.prev, 1, quantile, c(0.025, 0.975)), col=transp(col[1], 0.3)) + cred.region(xx, apply(dots_ancrt.prev[[ii]], 1, stats::quantile, c(0.025, 0.975)), col=transp(col[1+ii], 0.3)) + cred.region(xx, apply(ancrt.prev, 1, stats::quantile, c(0.025, 0.975)), col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, rowMeans(dots_ancrt.prev[[ii]]), col=col[1+ii], lwd=1.5) - lines(xx, rowMeans(ancrt.prev), col=col[1], lwd=1.5) + graphics::lines(xx, rowMeans(dots_ancrt.prev[[ii]]), col=col[1+ii], lwd=1.5) + graphics::lines(xx, rowMeans(ancrt.prev), col=col[1], lwd=1.5) ## if(!is.null(data$ancrtcens.dat)){ vinfl <- if(!is.null(fit$fp$ancrtcens.vinfl)) fit$fp$ancrtcens.vinfl else mean(ancrtcens.vinfl) - with(data$ancrtcens.dat, segments(x0=year, col="grey50", - y0=pnorm(qnorm(prev) - qnorm(0.975)*sqrt(v.ancrt+vinfl)), - y1=pnorm(qnorm(prev) + qnorm(0.975)*sqrt(v.ancrt+vinfl)))) - with(data$ancrtcens.dat, segments(year, - y0=pnorm(qnorm(prev) - qnorm(0.975)*sqrt(v.ancrt)), - y1=pnorm(qnorm(prev) + qnorm(0.975)*sqrt(v.ancrt)))) - with(data$ancrtcens.dat, points(year, prev, pch=20)) + with(data$ancrtcens.dat, graphics::segments(x0=year, col="grey50", + y0=stats::pnorm(stats::qnorm(prev) - stats::qnorm(0.975)*sqrt(v.ancrt+vinfl)), + y1=stats::pnorm(stats::qnorm(prev) + stats::qnorm(0.975)*sqrt(v.ancrt+vinfl)))) + with(data$ancrtcens.dat, graphics::segments(year, + y0=stats::pnorm(stats::qnorm(prev) - stats::qnorm(0.975)*sqrt(v.ancrt)), + y1=stats::pnorm(stats::qnorm(prev) + stats::qnorm(0.975)*sqrt(v.ancrt)))) + with(data$ancrtcens.dat, graphics::points(year, prev, pch=20)) } } @@ -198,16 +198,16 @@ plot_log_rvec <- function(fit, ..., ylim=NULL, xlim=c(1980, with(fit$fp$ss, proj for(ii in seq_along(dots)) dots[[ii]]$rvec <- sapply(seq_len(ncol(dots[[ii]]$rvec)), function(i) replace(dots[[ii]]$rvec[,i], dots[[ii]]$fp$proj.steps < dots[[ii]]$param[[i]]$tsEpidemicStart, NA)) if(is.null(ylim)) - ylim <- c(max(quantile(apply(log(fit$rvec), 1, quantile, 0.025, na.rm=TRUE), 0.025, na.rm=TRUE), -6), - min(quantile(apply(log(fit$rvec), 1, quantile, 0.975, na.rm=TRUE), 0.975, na.rm=TRUE), 5)) + ylim <- c(max(stats::quantile(apply(log(fit$rvec), 1, stats::quantile, 0.025, na.rm=TRUE), 0.025, na.rm=TRUE), -6), + min(stats::quantile(apply(log(fit$rvec), 1, stats::quantile, 0.975, na.rm=TRUE), 0.975, na.rm=TRUE), 5)) plot(xx, rowMeans(log(fit$rvec), na.rm=TRUE)[idx], type="n", ylim=ylim, xlim=xlim, ylab="log r(t)", yaxt="n", xlab="") - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) for(ii in seq_along(dots)) - cred.region(xx, apply(log(dots[[ii]]$rvec[idx,]), 1, quantile, c(0.025, 0.975), na.rm=TRUE), col=transp(col[1+ii], 0.3)) - cred.region(xx, apply(log(fit$rvec[idx,]), 1, quantile, c(0.025, 0.975), na.rm=TRUE), col=transp(col[1], 0.3)) + cred.region(xx, apply(log(dots[[ii]]$rvec[idx,]), 1, stats::quantile, c(0.025, 0.975), na.rm=TRUE), col=transp(col[1+ii], 0.3)) + cred.region(xx, apply(log(fit$rvec[idx,]), 1, stats::quantile, c(0.025, 0.975), na.rm=TRUE), col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, rowMeans(log(dots[[ii]]$rvec[idx,]), na.rm=TRUE)[idx], col=col[1+ii]) - lines(xx, rowMeans(log(fit$rvec[idx,]), na.rm=TRUE)[idx], col=col[1]) + graphics::lines(xx, rowMeans(log(dots[[ii]]$rvec[idx,]), na.rm=TRUE)[idx], col=col[1+ii]) + graphics::lines(xx, rowMeans(log(fit$rvec[idx,]), na.rm=TRUE)[idx], col=col[1]) return(invisible()) } diff --git a/R/plot2.R b/R/plot2.R index 0bb6e0a..e2fc311 100644 --- a/R/plot2.R +++ b/R/plot2.R @@ -4,17 +4,17 @@ plot_prev2 <- function(fit, ..., ylim=NULL, xlim=c(1980, max(as.integer(dimnames ylim <- c(0, 1.1*max(fit$prev[,"upper"])) xx <- as.integer(dimnames(fit$prev)[[1]]) plot(xx, fit$prev[,"mean"], type="n", ylim=ylim, xlim=xlim, ylab=ylab, xlab="", yaxt="n", xaxt="n", main=main) - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) dots <- list(...) dots <- dots[!sapply(dots, is.null)] - + for(ii in seq_along(dots)) cred.region(xx, dots[[ii]]$prev[,c("lower", "upper")], col=transp(col[1+ii], 0.3)) cred.region(xx, fit$prev[,c("lower", "upper")], col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, dots[[ii]]$prev[,"mean"], col=col[1+ii], lwd=1.5) - lines(xx, fit$prev[,"mean"], col=col[1], lwd=1.5) + graphics::lines(xx, dots[[ii]]$prev[,"mean"], col=col[1+ii], lwd=1.5) + graphics::lines(xx, fit$prev[,"mean"], col=col[1], lwd=1.5) } @@ -24,17 +24,17 @@ plot_incid2 <- function(fit, ..., ylim=NULL, xlim=c(1980, max(as.integer(dimname ylim <- c(0, 1.1*max(fit$incid[,"upper"])) xx <- as.integer(dimnames(fit$incid)[[1]]) plot(xx, fit$incid[,"mean"], type="n", ylim=ylim, xlim=xlim, ylab=ylab, xlab="", yaxt="n", xaxt="n", main=main) - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) dots <- list(...) dots <- dots[!sapply(dots, is.null)] - + for(ii in seq_along(dots)) cred.region(xx, dots[[ii]]$incid[,c("lower", "upper")], col=transp(col[1+ii], 0.3)) cred.region(xx, fit$incid[,c("lower", "upper")], col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, dots[[ii]]$incid[,"mean"], col=col[1+ii], lwd=1.5) - lines(xx, fit$incid[,"mean"], col=col[1], lwd=1.5) + graphics::lines(xx, dots[[ii]]$incid[,"mean"], col=col[1+ii], lwd=1.5) + graphics::lines(xx, fit$incid[,"mean"], col=col[1], lwd=1.5) } @@ -45,17 +45,17 @@ plot_log_transmrate <- function(fit, ..., ylim=NULL, xlim=c(1980, max(as.integer max(log(fit$transmrate[,"upper"])) + 0.2) xx <- as.integer(dimnames(fit$transmrate)[[1]]) plot(xx, fit$transmrate[,"mean"], type="n", ylim=ylim, xlim=xlim, ylab=ylab, xlab="", yaxt="n", xaxt="n", main=main) - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) dots <- list(...) dots <- dots[!sapply(dots, is.null)] - + for(ii in seq_along(dots)) cred.region(xx, log(dots[[ii]]$transmrate[,c("lower", "upper")]), col=transp(col[1+ii], 0.3)) cred.region(xx, log(fit$transmrate[,c("lower", "upper")]), col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, log(dots[[ii]]$transmrate[,"mean"]), col=col[1+ii], lwd=1.5) - lines(xx, log(fit$transmrate[,"mean"]), col=col[1], lwd=1.5) + graphics::lines(xx, log(dots[[ii]]$transmrate[,"mean"]), col=col[1+ii], lwd=1.5) + graphics::lines(xx, log(fit$transmrate[,"mean"]), col=col[1], lwd=1.5) } @@ -65,45 +65,45 @@ plot_incidsexratio <- function(fit, ..., ylim=NULL, xlim=c(1999, max(as.integer( ylim <- c(0, max(2.5, 1.1*max(fit$incidsexratio[,"upper"]))) xx <- as.integer(dimnames(fit$incidsexratio)[[1]]) plot(xx, fit$incidsexratio[,"mean"], type="n", ylim=ylim, xlim=xlim, ylab=ylab, xlab="", yaxt="n", xaxt="n", main=main) - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) dots <- list(...) dots <- dots[!sapply(dots, is.null)] - + for(ii in seq_along(dots)) cred.region(xx, dots[[ii]]$incidsexratio[,c("lower", "upper")], col=transp(col[1+ii], 0.3)) cred.region(xx, fit$incidsexratio[,c("lower", "upper")], col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, dots[[ii]]$incidsexratio[,"mean"], col=col[1+ii], lwd=1.5) - lines(xx, fit$incidsexratio[,"mean"], col=col[1], lwd=1.5) + graphics::lines(xx, dots[[ii]]$incidsexratio[,"mean"], col=col[1+ii], lwd=1.5) + graphics::lines(xx, fit$incidsexratio[,"mean"], col=col[1], lwd=1.5) } plot_pregprev <- function(fit, ..., likdat=NULL, ylim=NULL, xlim=c(1988, max(as.integer(dimnames(fit$pregprev)[[1]]))), col="blue", main="", ylab="Preg. prevalence"){ - + dots <- list(...) dots <- dots[!sapply(dots, is.null)] - + if(is.null(ylim)){ maxest <- max(fit$pregprev[,"upper"]) if(!is.null(likdat) && !is.null(likdat$anclik.dat) && length(likdat$anclik.dat$W.lst)) - maxdata <- max(pnorm(unlist(likdat$anclik.dat$W.lst))) + maxdata <- max(stats::pnorm(unlist(likdat$anclik.dat$W.lst))) else maxdata <- 0 ylim <- c(0, 1.1*max(maxest, min(maxdata, 2*maxest))) } xx <- as.integer(dimnames(fit$incidsexratio)[[1]]) - + plot(xx, fit$pregprev[,"mean"], type="n", ylim=ylim, xlim=xlim, ylab=ylab, xlab="", yaxt="n", xaxt="n", main=main) - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) ## if(!is.null(likdat)){ if(!is.null(likdat$anclik.dat) && length(likdat$anclik.dat$W.lst)){ - with(likdat$anclik.dat, mapply(function(idx, W) points(idx+1970-1, pnorm(W), col=adjustcolor("grey", 0.5), pch=15), anc.idx.lst, W.lst)) - with(likdat$anclik.dat, mapply(function(idx, W) lines(idx+1970-1, pnorm(W), col=adjustcolor("grey", 0.5)), anc.idx.lst, W.lst)) + with(likdat$anclik.dat, mapply(function(idx, W) graphics::points(idx+1970-1, stats::pnorm(W), col=grDevices::adjustcolor("grey", 0.5), pch=15), anc.idx.lst, W.lst)) + with(likdat$anclik.dat, mapply(function(idx, W) graphics::lines(idx+1970-1, stats::pnorm(W), col=grDevices::adjustcolor("grey", 0.5)), anc.idx.lst, W.lst)) } } ## @@ -111,14 +111,14 @@ plot_pregprev <- function(fit, ..., likdat=NULL, ylim=NULL, xlim=c(1988, max(as. cred.region(xx, dots[[ii]]$pregprev[,c("lower", "upper")], col=transp(col[1+ii], 0.3)) cred.region(xx, fit$pregprev[,c("lower", "upper")], col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, dots[[ii]]$pregprev[,"mean"], col=col[1+ii], lwd=1.5) - lines(xx, fit$pregprev[,"mean"], col=col[1], lwd=1.5) + graphics::lines(xx, dots[[ii]]$pregprev[,"mean"], col=col[1+ii], lwd=1.5) + graphics::lines(xx, fit$pregprev[,"mean"], col=col[1], lwd=1.5) ## if(!is.null(likdat$ancrtcens.dat) && length(likdat$ancrtcens.dat$W.ancrt)){ with(likdat$ancrtcens.dat, - segments(year, y0=pnorm(qnorm(prev) - qnorm(0.975)*sqrt(v.ancrt)), - y1=pnorm(qnorm(prev) + qnorm(0.975)*sqrt(v.ancrt)))) - with(likdat$ancrtcens.dat, points(year, prev, pch=15)) + graphics::segments(year, y0=stats::pnorm(stats::qnorm(prev) - stats::qnorm(0.975)*sqrt(v.ancrt)), + y1=stats::pnorm(stats::qnorm(prev) + stats::qnorm(0.975)*sqrt(v.ancrt)))) + with(likdat$ancrtcens.dat, graphics::points(year, prev, pch=15)) } } @@ -129,17 +129,17 @@ plot_artcov15plus <- function(fit, ..., ylim=NULL, xlim=c(2003, max(as.integer(d ylim <- c(0, 1) xx <- as.integer(dimnames(fit$artcov15plus)[[1]]) plot(xx, fit$artcov15plus[,"mean"], type="n", ylim=ylim, xlim=xlim, ylab=ylab, xlab="", yaxt="n", xaxt="n", main=main) - axis(1, labels=TRUE) - axis(2, labels=TRUE) + graphics::axis(1, labels=TRUE) + graphics::axis(2, labels=TRUE) dots <- list(...) dots <- dots[!sapply(dots, is.null)] - + for(ii in seq_along(dots)) cred.region(xx, dots[[ii]]$artcov15plus[,c("lower", "upper")], col=transp(col[1+ii], 0.3)) cred.region(xx, fit$artcov15plus[,c("lower", "upper")], col=transp(col[1], 0.3)) for(ii in seq_along(dots)) - lines(xx, dots[[ii]]$artcov15plus[,"mean"], col=col[1+ii], lwd=1.5) - lines(xx, fit$artcov15plus[,"mean"], col=col[1], lwd=1.5) + graphics::lines(xx, dots[[ii]]$artcov15plus[,"mean"], col=col[1+ii], lwd=1.5) + graphics::lines(xx, fit$artcov15plus[,"mean"], col=col[1], lwd=1.5) } @@ -160,14 +160,14 @@ plot_compare_ageprev2 <- function(fit, fit2=NULL, fit3=NULL, specres=NULL, likda survprev[c("ci_l", "ci_u")] <- with(survprev, binom::binom.exact(x_eff, n_eff))[c("lower", "upper")] } survprev$survyear <- with(survprev, factor(survyear, levels(survyear)[order(as.integer(substr(levels(survyear), 1, 4)))])) - survprev <- split(survprev, factor(survprev$survyear)) + survprev <- split(survprev, factor(survprev$survyear)) ## if(!is.null(fit2)) survprev2 <- split(fit2$ageprevdat, factor(fit2$ageprevdat$survyear)) if(!is.null(fit3)) survprev3 <- split(fit3$ageprevdat, factor(fit3$ageprevdat$survyear)) ## - par(mfrow=c(4,2), mar=c(2, 3, 2, 1), tcl=-0.25, mgp=c(2, 0.5, 0), las=1, cex=1) + graphics::par(mfrow=c(4,2), mar=c(2, 3, 2, 1), tcl=-0.25, mgp=c(2, 0.5, 0), las=1, cex=1) for(isurv in names(survprev)) for(isex in c("male", "female")){ sp <- subset(survprev[[isurv]], sex==isex & as.integer(agegr) %in% 3:11) @@ -183,32 +183,32 @@ plot_compare_ageprev2 <- function(fit, fit2=NULL, fit3=NULL, specres=NULL, likda paste0(sp$country[1], " ", survprev[[isurv]]$survyear[1], ", ", isex) plot(xx+0.5, sp$prev, type="n", xlim=c(4, 12), ylim=ylim, xaxt="n", main=main, xlab="", ylab="") - axis(1, xx+0.5, sp$agegr) + graphics::axis(1, xx+0.5, sp$agegr) ## - rect(xx+0.05, sp$lower, xx+0.95, sp$upper, + graphics::rect(xx+0.05, sp$lower, xx+0.95, sp$upper, col=transp(col[1]), border=NA) - segments(xx+0.05, sp$mean, xx+0.95, col=col[1], lwd=2) + graphics::segments(xx+0.05, sp$mean, xx+0.95, col=col[1], lwd=2) ## if(!is.null(fit2)){ - rect(xx+0.05, sp2$lower, xx+0.95, sp2$upper, + graphics::rect(xx+0.05, sp2$lower, xx+0.95, sp2$upper, col=transp(col[2]), border=NA) - segments(xx+0.05, sp2$mean, xx+0.95, col=col[2], lwd=2) + graphics::segments(xx+0.05, sp2$mean, xx+0.95, col=col[2], lwd=2) } if(!is.null(fit3)){ - rect(xx+0.05, sp3$lower, xx+0.95, sp3$upper, + graphics::rect(xx+0.05, sp3$lower, xx+0.95, sp3$upper, col=transp(col[3]), border=NA) - segments(xx+0.05, sp3$mean, xx+0.95, col=col[3], lwd=2) + graphics::segments(xx+0.05, sp3$mean, xx+0.95, col=col[3], lwd=2) } ## if(!is.null(specres)){ csex <- sub("(\\b[a-z]{1})", "\\U\\1" , isex, perl=TRUE) stryear <- as.character(survprev[[isurv]]$year[1]) specres.prev <- tapply(specres$hivpop[as.character(15:54), csex, stryear], rep(3:10, each=5), sum) / tapply(specres$totpop[as.character(15:54), csex, stryear], rep(3:10, each=5), sum) - segments(4:11+0.1, specres.prev, 4:11+0.9, lty=3, col="grey10", lwd=2) + graphics::segments(4:11+0.1, specres.prev, 4:11+0.9, lty=3, col="grey10", lwd=2) } if(exists("prev", sp)){ - points(xx+0.5, sp$prev, pch=19) - segments(x0=xx+0.5, y0=sp$ci_l, y1=sp$ci_u) + graphics::points(xx+0.5, sp$prev, pch=19) + graphics::segments(x0=xx+0.5, y0=sp$ci_l, y1=sp$ci_u) } } ## diff --git a/R/plot_agefit.R b/R/plot_agefit.R index c1f1325..b653118 100644 --- a/R/plot_agefit.R +++ b/R/plot_agefit.R @@ -1,7 +1,7 @@ plot_agefit <- function(icountry, eppmod, fitaggr, fitincrr, fit3=NULL, specres=NULL, cols=c("grey30", "darkred", "darkgreen"), pages=1:3, agegr3dat = subset(prev_agegr3sex_nat, country==icountry), age15to49dat=NULL){ if(1 %in% pages){ - par(oma=c(1, 1, 2.5, 1), mfrow=c(3, 2), mar=c(3.5, 3.5, 3, 1), tcl=-0.25, mgp=c(2.5, 0.5, 0), cex=1, las=1) + graphics::par(oma=c(1, 1, 2.5, 1), mfrow=c(3, 2), mar=c(3.5, 3.5, 3, 1), tcl=-0.25, mgp=c(2.5, 0.5, 0), cex=1, las=1) ## ## prevalence trend if(!is.null(fit3)) @@ -9,18 +9,18 @@ plot_agefit <- function(icountry, eppmod, fitaggr, fitincrr, fit3=NULL, specres= else plot_prev(fitaggr, fitincrr, col=cols, plotprevdat=is.null(age15to49dat)) if(!is.null(specres)) - lines(as.numeric(names(prev(specres))), prev(specres), lty=2, lwd=2, col="grey10") + graphics::lines(as.numeric(names(prev(specres))), prev(specres), lty=2, lwd=2, col="grey10") if(!is.null(age15to49dat)){ if(!is.null(age15to49dat$prev_spec17)){ - points(age15to49dat$year, age15to49dat$prev, pch=15, col="grey40") - points(age15to49dat$year, age15to49dat$prev_spec17, pch=19) - segments(age15to49dat$year, y0=age15to49dat$ci_l_spec17, y1=age15to49dat$ci_u_spec17) + graphics::points(age15to49dat$year, age15to49dat$prev, pch=15, col="grey40") + graphics::points(age15to49dat$year, age15to49dat$prev_spec17, pch=19) + graphics::segments(age15to49dat$year, y0=age15to49dat$ci_l_spec17, y1=age15to49dat$ci_u_spec17) } else { - points(age15to49dat$year, age15to49dat$prev, pch=19) - segments(age15to49dat$year, y0=age15to49dat$ci_l, y1=age15to49dat$ci_u) + graphics::points(age15to49dat$year, age15to49dat$prev, pch=19) + graphics::segments(age15to49dat$year, y0=age15to49dat$ci_l, y1=age15to49dat$ci_u) } } - mtext("HIV prevalence, age 15-49y", 3, 0.5, font=2, cex=1.2) + graphics::mtext("HIV prevalence, age 15-49y", 3, 0.5, font=2, cex=1.2) ## ## incidence trend if(!is.null(fit3)) @@ -28,8 +28,8 @@ plot_agefit <- function(icountry, eppmod, fitaggr, fitincrr, fit3=NULL, specres= else plot_incid(fitaggr, fitincrr, col=cols) if(!is.null(specres)) - lines(as.numeric(names(incid(specres))), incid(specres), lty=2, lwd=2, col="grey10") - mtext("HIV incidence rate, age 15-49y", 3, 0.5, font=2, cex=1.2) + graphics::lines(as.numeric(names(incid(specres))), incid(specres), lty=2, lwd=2, col="grey10") + graphics::mtext("HIV incidence rate, age 15-49y", 3, 0.5, font=2, cex=1.2) ## ## r-trend if(!is.null(fitaggr$rvec[[1]])){ @@ -37,7 +37,7 @@ plot_agefit <- function(icountry, eppmod, fitaggr, fitincrr, fit3=NULL, specres= plot_log_rvec(fitaggr, fitincrr, fit3, col=cols) else plot_log_rvec(fitaggr, fitincrr, col=cols) - mtext("log r(t)", 3, 0.5, font=2, cex=1.2) + graphics::mtext("log r(t)", 3, 0.5, font=2, cex=1.2) } ## ## vinfl @@ -45,31 +45,31 @@ plot_agefit <- function(icountry, eppmod, fitaggr, fitincrr, fit3=NULL, specres= ## sex incrr if(!is.null(fitincrr$param)){ yidx <- 2010-fitincrr$fp$ss$proj_start+1 - logincrrsex <- log(sapply(fitincrr$param, "[[", "incrr_sex")[yidx,]) # year 2010 - dens <- density(logincrrsex) - densCI <- which(dens$x >= quantile(logincrrsex, 0.025) & - dens$x <= quantile(logincrrsex, 0.975)) + logincrrsex <- log(sapply(fitincrr$param, "[[", "incrr_sex")[yidx,]) # year 2010 + dens <- stats::density(logincrrsex) + densCI <- which(dens$x >= stats::quantile(logincrrsex, 0.025) & + dens$x <= stats::quantile(logincrrsex, 0.975)) if(!is.null(fit3)){ logincrrsex3 <- log(sapply(fit3$param, "[[", "incrr_sex")[yidx,]) - dens3 <- density(logincrrsex3) - dens3CI <- which(dens3$x >= quantile(logincrrsex3, 0.025) & - dens3$x <= quantile(logincrrsex3, 0.975)) + dens3 <- stats::density(logincrrsex3) + dens3CI <- which(dens3$x >= stats::quantile(logincrrsex3, 0.025) & + dens3$x <= stats::quantile(logincrrsex3, 0.975)) } plot(dens, xlim=c(-0.1, 0.75), col=cols[2], main="F:M incidence RR (2010), posterior density", xlab="log F:M IRR") - polygon(dens$x[c(min(densCI), densCI, max(densCI))], c(0, dens$y[densCI], 0), + graphics::polygon(dens$x[c(min(densCI), densCI, max(densCI))], c(0, dens$y[densCI], 0), col=transp(cols[2]), border=NA) if(!is.null(fit3)){ - lines(dens3, col=cols[3]) - polygon(dens3$x[c(min(dens3CI), dens3CI, max(dens3CI))], c(0, dens3$y[dens3CI], 0), + graphics::lines(dens3, col=cols[3]) + graphics::polygon(dens3$x[c(min(dens3CI), dens3CI, max(dens3CI))], c(0, dens3$y[dens3CI], 0), col=transp(cols[3]), border=NA) } - segments(x0=mean(log(sapply(fitincrr$param, "[[", "incrr_sex")[yidx,])), y0=0, y1=6, col=cols[2], lwd=2) + graphics::segments(x0=mean(log(sapply(fitincrr$param, "[[", "incrr_sex")[yidx,])), y0=0, y1=6, col=cols[2], lwd=2) if(!is.null(fit3)) - segments(x0=mean(log(sapply(fit3$param, "[[", "incrr_sex")[yidx,])), y0=0, y1=6, col=cols[3], lwd=2) - segments(x0=log(tail(fitaggr$fp$incrr_sex, 1)), y0=0, y1=6, col=cols[1], lwd=2) - lines(seq(-0.1, 0.75, 0.01), dnorm(seq(-0.1, 0.75, 0.01), eppasm:::sexincrr.pr.mean, eppasm:::sexincrr.pr.sd), col="darkblue", lty=2) - segments(x0=eppasm:::sexincrr.pr.mean, y0=0, y1=2, col="darkblue", lwd=2, lty=2) + graphics::segments(x0=mean(log(sapply(fit3$param, "[[", "incrr_sex")[yidx,])), y0=0, y1=6, col=cols[3], lwd=2) + graphics::segments(x0=log(utils::tail(fitaggr$fp$incrr_sex, 1)), y0=0, y1=6, col=cols[1], lwd=2) + graphics::lines(seq(-0.1, 0.75, 0.01), stats::dnorm(seq(-0.1, 0.75, 0.01), eppasm:::sexincrr.pr.mean, eppasm:::sexincrr.pr.sd), col="darkblue", lty=2) + graphics::segments(x0=eppasm:::sexincrr.pr.mean, y0=0, y1=2, col="darkblue", lwd=2, lty=2) ## ## age incrr @@ -80,65 +80,65 @@ plot_agefit <- function(icountry, eppmod, fitaggr, fitincrr, fit3=NULL, specres= ## plot men plot(1:7+0.5, 1:7, type="n", xlim=c(1, 8), ylim=c(-2.5, 2), xlab="Age group", ylab="log IRR", main = "Male incidence RR (log), 2007", xaxt="n") - axis(1, 1:7+0.5, paste0(3:9*5, "-", 3:9*5+4)) - abline(h=0, col="grey") - points(3.5, 0, pch=4, lwd=2.5, col=1, cex=1.2) - rect(xx+0.1, - eppasm:::ageincrr.pr.mean[1:6]-qnorm(0.975)*eppasm:::ageincrr.pr.sd, + graphics::axis(1, 1:7+0.5, paste0(3:9*5, "-", 3:9*5+4)) + graphics::abline(h=0, col="grey") + graphics::points(3.5, 0, pch=4, lwd=2.5, col=1, cex=1.2) + graphics::rect(xx+0.1, + eppasm:::ageincrr.pr.mean[1:6]-stats::qnorm(0.975)*eppasm:::ageincrr.pr.sd, xx+0.9, - eppasm:::ageincrr.pr.mean[1:6]+qnorm(0.975)*eppasm:::ageincrr.pr.sd, + eppasm:::ageincrr.pr.mean[1:6]+stats::qnorm(0.975)*eppasm:::ageincrr.pr.sd, col=transp("darkblue", 0.1), border=transp("darkblue", 0.5), lty=2) ## - rect(xx+0.1, logincrrage[xx,3], xx+0.9, logincrrage[xx,4], col=transp(cols[2]), border=NA) + graphics::rect(xx+0.1, logincrrage[xx,3], xx+0.9, logincrrage[xx,4], col=transp(cols[2]), border=NA) if(!is.null(fit3)) - rect(xx+0.1, logincrrage3[xx,3], xx+0.9, logincrrage3[xx,4], col=transp(cols[3]), border=NA) + graphics::rect(xx+0.1, logincrrage3[xx,3], xx+0.9, logincrrage3[xx,4], col=transp(cols[3]), border=NA) defaultincrr <- log(fitaggr$fp$incrr_age[(xx-1)*5+1,1,yidx]) - segments(xx+0.1, defaultincrr, xx+0.9, col=cols[1], lwd=2) - segments(xx+0.1, eppasm:::ageincrr.pr.mean[1:6], xx+0.9, col=transp("darkblue", 0.5), lwd=1.5, lty=2) - segments(xx+0.1, logincrrage[xx,1], xx+0.9, col=cols[2], lwd=2) + graphics::segments(xx+0.1, defaultincrr, xx+0.9, col=cols[1], lwd=2) + graphics::segments(xx+0.1, eppasm:::ageincrr.pr.mean[1:6], xx+0.9, col=transp("darkblue", 0.5), lwd=1.5, lty=2) + graphics::segments(xx+0.1, logincrrage[xx,1], xx+0.9, col=cols[2], lwd=2) if(!is.null(fit3)) - segments(xx+0.1, logincrrage3[xx,1], xx+0.9, col=cols[3], lwd=2) + graphics::segments(xx+0.1, logincrrage3[xx,1], xx+0.9, col=cols[3], lwd=2) ## plot women plot(1:7+0.5, 1:7, type="n", xlim=c(1, 8), ylim=c(-2.5, 2), xlab="Age group", ylab="log IRR", main = "Female incidence RR (log), 2007", xaxt="n") - axis(1, 1:7+0.5, paste0(3:9*5, "-", 3:9*5+4)) - abline(h=0, col="grey") - points(3.5, 0, pch=4, lwd=2.5, col=1, cex=1.2) - rect(xx+0.1, - eppasm:::ageincrr.pr.mean[7:12]-qnorm(0.975)*eppasm:::ageincrr.pr.sd, + graphics::axis(1, 1:7+0.5, paste0(3:9*5, "-", 3:9*5+4)) + graphics::abline(h=0, col="grey") + graphics::points(3.5, 0, pch=4, lwd=2.5, col=1, cex=1.2) + graphics::rect(xx+0.1, + eppasm:::ageincrr.pr.mean[7:12]-stats::qnorm(0.975)*eppasm:::ageincrr.pr.sd, xx+0.9, - eppasm:::ageincrr.pr.mean[7:12]+qnorm(0.975)*eppasm:::ageincrr.pr.sd, + eppasm:::ageincrr.pr.mean[7:12]+stats::qnorm(0.975)*eppasm:::ageincrr.pr.sd, col=transp("darkblue", 0.1), border=transp("darkblue", 0.5), lty=2) ## offset <- nrow(logincrrage)/2 - rect(xx+0.1, logincrrage[xx+offset,3], xx+0.9, logincrrage[xx++offset,4], col=transp(cols[2]), border=NA) + graphics::rect(xx+0.1, logincrrage[xx+offset,3], xx+0.9, logincrrage[xx++offset,4], col=transp(cols[2]), border=NA) if(!is.null(fit3)){ offset3 <- nrow(logincrrage3)/2 - rect(xx+0.1, logincrrage3[xx+offset3,3], xx+0.9, logincrrage3[xx++offset3,4], col=transp(cols[3]), border=NA) + graphics::rect(xx+0.1, logincrrage3[xx+offset3,3], xx+0.9, logincrrage3[xx++offset3,4], col=transp(cols[3]), border=NA) } defaultincrr <- log(fitaggr$fp$incrr_age[(xx-1)*5+1,2,yidx]) - segments(xx+0.1, defaultincrr, xx+0.9, col=cols[1], lwd=2) - segments(xx+0.1, eppasm:::ageincrr.pr.mean[7:12], xx+0.9, col=transp("darkblue", 0.5), lwd=1.5, lty=2) - segments(xx+0.1, logincrrage[xx+offset,1], xx+0.9, col=cols[2], lwd=2) + graphics::segments(xx+0.1, defaultincrr, xx+0.9, col=cols[1], lwd=2) + graphics::segments(xx+0.1, eppasm:::ageincrr.pr.mean[7:12], xx+0.9, col=transp("darkblue", 0.5), lwd=1.5, lty=2) + graphics::segments(xx+0.1, logincrrage[xx+offset,1], xx+0.9, col=cols[2], lwd=2) if(!is.null(fit3)) - segments(xx+0.1, logincrrage3[xx+offset3,1], xx+0.9, col=cols[3], lwd=2) + graphics::segments(xx+0.1, logincrrage3[xx+offset3,1], xx+0.9, col=cols[3], lwd=2) ## } - mtext(paste0(icountry, ", ", eppmod, "; posterior distribution"), 3, 0.5, outer=TRUE, font=2, cex=1.3) + graphics::mtext(paste0(icountry, ", ", eppmod, "; posterior distribution"), 3, 0.5, outer=TRUE, font=2, cex=1.3) } ## ## Age specific prevalence compared to survey ## if(2 %in% pages){ plot_compare_ageprev(fitaggr, fitincrr, fit3, specres, col=cols) - mtext(paste0(icountry, ", ", eppmod, "; Age-specific prevalence"), 3, 0.5, outer=TRUE, font=2, cex=1.3) + graphics::mtext(paste0(icountry, ", ", eppmod, "; Age-specific prevalence"), 3, 0.5, outer=TRUE, font=2, cex=1.3) } ## ## ## Age prevalence trend by age group ## if(3 %in% pages){ - par(oma=c(1, 1, 2.5, 1), mfrow=c(3, 2), mar=c(3.5, 3.5, 3, 1), tcl=-0.25, mgp=c(2.5, 0.5, 0), cex=1, las=1) + graphics::par(oma=c(1, 1, 2.5, 1), mfrow=c(3, 2), mar=c(3.5, 3.5, 3, 1), tcl=-0.25, mgp=c(2.5, 0.5, 0), cex=1, las=1) ## for(iagegr in 1:3){ for(isex in 1:2){ @@ -157,28 +157,28 @@ plot_agefit <- function(icountry, eppmod, fitaggr, fitincrr, fit3=NULL, specres= cred.region(xx, t(fit.yprev[,3:4]), col=transp(cols[2], 0.4)) if(!is.null(fit3)) cred.region(xx, t(fit3.yprev[,3:4]), col=transp(cols[3], 0.4)) - lines(xx, aggr.yprev[,1], col=cols[1], lwd=2) - lines(xx, fit.yprev[,1], col=cols[2], lwd=2) + graphics::lines(xx, aggr.yprev[,1], col=cols[1], lwd=2) + graphics::lines(xx, fit.yprev[,1], col=cols[2], lwd=2) if(!is.null(fit3)) - lines(xx, fit3.yprev[,1], col=cols[3], lwd=2) + graphics::lines(xx, fit3.yprev[,1], col=cols[3], lwd=2) ## if(!is.null(specres)){ ages <- as.character(c(15, 25, 35)[iagegr] + 0:(c(10, 10, 15)[iagegr]-1)) specres.prev <- colSums(specres$hivpop[ages,isex,as.character(xx)]) / colSums(specres$totpop[ages,isex,as.character(xx)]) - lines(xx, specres.prev, lty=2, lwd=2, col="grey10") + graphics::lines(xx, specres.prev, lty=2, lwd=2, col="grey10") } ## if(!is.null(survdat$prev_spec17)){ - points(survdat$year, survdat$prev, pch=15, col="grey40") - points(survdat$year, survdat$prev_spec17, pch=19) - segments(survdat$year, y0=survdat$ci_l_spec17, y1=survdat$ci_u_spec17) + graphics::points(survdat$year, survdat$prev, pch=15, col="grey40") + graphics::points(survdat$year, survdat$prev_spec17, pch=19) + graphics::segments(survdat$year, y0=survdat$ci_l_spec17, y1=survdat$ci_u_spec17) } else { - points(survdat$year, survdat$prev, pch=19) - segments(survdat$year, y0=survdat$ci_l, y1=survdat$ci_u) + graphics::points(survdat$year, survdat$prev, pch=19) + graphics::segments(survdat$year, y0=survdat$ci_l, y1=survdat$ci_u) } } } - mtext(paste0(icountry, ", ", eppmod, "\nPrevalence trend by age group"), 3, 0, outer=TRUE, font=2, cex=1.3) + graphics::mtext(paste0(icountry, ", ", eppmod, "\nPrevalence trend by age group"), 3, 0, outer=TRUE, font=2, cex=1.3) } ## ## diff --git a/R/plot_output.R b/R/plot_output.R index 2b764be..75cb5ea 100644 --- a/R/plot_output.R +++ b/R/plot_output.R @@ -4,55 +4,55 @@ plot_output <- function(fit, fit2=NULL, fit3=NULL, specres=NULL, cols=c("navy", title="", agegr3dat=NULL, age15to49dat=NULL){ if(1 %in% pages){ - par(oma=c(1, 1, 2.5, 1), mfrow=c(3, 2), mar=c(3.5, 3.5, 3, 1), tcl=-0.25, mgp=c(2.5, 0.5, 0), cex=1, las=1) + graphics::par(oma=c(1, 1, 2.5, 1), mfrow=c(3, 2), mar=c(3.5, 3.5, 3, 1), tcl=-0.25, mgp=c(2.5, 0.5, 0), cex=1, las=1) ## ## prevalence trend plot_prev2(fit, fit2, fit3, col=cols) if(!is.null(specres)) - lines(as.numeric(names(prev(specres))), prev(specres), lty=2, lwd=2, col="grey10") + graphics::lines(as.numeric(names(prev(specres))), prev(specres), lty=2, lwd=2, col="grey10") if(!is.null(age15to49dat)){ if(!is.null(age15to49dat$prev_spec17)){ - points(age15to49dat$year, age15to49dat$prev, pch=15, col="grey40") - points(age15to49dat$year, age15to49dat$prev_spec17, pch=19) - segments(age15to49dat$year, y0=age15to49dat$ci_l_spec17, y1=age15to49dat$ci_u_spec17) + graphics::points(age15to49dat$year, age15to49dat$prev, pch=15, col="grey40") + graphics::points(age15to49dat$year, age15to49dat$prev_spec17, pch=19) + graphics::segments(age15to49dat$year, y0=age15to49dat$ci_l_spec17, y1=age15to49dat$ci_u_spec17) } else { - points(age15to49dat$year, age15to49dat$prev, pch=19) - segments(age15to49dat$year, y0=age15to49dat$ci_l, y1=age15to49dat$ci_u) + graphics::points(age15to49dat$year, age15to49dat$prev, pch=19) + graphics::segments(age15to49dat$year, y0=age15to49dat$ci_l, y1=age15to49dat$ci_u) } } - mtext("HIV prevalence, age 15-49y", 3, 0.5, font=2, cex=1.1) + graphics::mtext("HIV prevalence, age 15-49y", 3, 0.5, font=2, cex=1.1) ## ## incidence trend plot_incid2(fit, fit2, fit3, col=cols) if(!is.null(specres)) - lines(as.numeric(names(incid(specres))), incid(specres), lty=2, lwd=2, col="grey10") - mtext("HIV incidence rate, age 15-49y", 3, 0.5, font=2, cex=1.1) + graphics::lines(as.numeric(names(incid(specres))), incid(specres), lty=2, lwd=2, col="grey10") + graphics::mtext("HIV incidence rate, age 15-49y", 3, 0.5, font=2, cex=1.1) ## ## tranmission rate plot_log_transmrate(fit, fit2, fit3, col=cols) if(!is.null(specres)) - lines(as.numeric(names(incid(specres))), log(incid(specres) / prev(specres)), lty=2, lwd=2, col="grey10") - mtext("Log transmission rate", 3, 0.5, font=2, cex=1.1) + graphics::lines(as.numeric(names(incid(specres))), log(incid(specres) / prev(specres)), lty=2, lwd=2, col="grey10") + graphics::mtext("Log transmission rate", 3, 0.5, font=2, cex=1.1) ## ## Sex ratio of incidence plot_incidsexratio(fit, fit2, fit3, col=cols) if(!is.null(specres)) - lines(as.numeric(names(incid_sexratio(specres))), incid_sexratio(specres), lty=2, lwd=2, col="grey10") - mtext("Incidence sex ratio (F:M, 15-49y)", 3, 0.5, font=2, cex=1.1) + graphics::lines(as.numeric(names(incid_sexratio(specres))), incid_sexratio(specres), lty=2, lwd=2, col="grey10") + graphics::mtext("Incidence sex ratio (F:M, 15-49y)", 3, 0.5, font=2, cex=1.1) ## ## ART coverage age 15+ plot_artcov15plus(fit, fit2, fit3, col=cols) if(!is.null(specres) && exists("artnum.m", specres)) - lines(as.numeric(names(artcov15plus(specres))), artcov15plus(specres), lty=2, lwd=2, col="grey10") - mtext("ART coverage, age 15+", 3, 0.5, font=2, cex=1.1) - mtext(paste0(title, "; posterior distribution"), 3, 0.5, outer=TRUE, font=2, cex=1.3) + graphics::lines(as.numeric(names(artcov15plus(specres))), artcov15plus(specres), lty=2, lwd=2, col="grey10") + graphics::mtext("ART coverage, age 15+", 3, 0.5, font=2, cex=1.1) + graphics::mtext(paste0(title, "; posterior distribution"), 3, 0.5, outer=TRUE, font=2, cex=1.3) ## ## Pregnant women prevalence if(!is.null(fit$pregprev)){ plot_pregprev(fit, fit2, fit3, col=cols, likdat=likdat) if(!is.null(specres)) - lines(as.numeric(names(fnPregPrev(specres))), fnPregPrev(specres), lty=2, lwd=2, col="grey10") - mtext("Prevalence among pregnant women", 3, 0.5, font=2, cex=1.1) + graphics::lines(as.numeric(names(fnPregPrev(specres))), fnPregPrev(specres), lty=2, lwd=2, col="grey10") + graphics::mtext("Prevalence among pregnant women", 3, 0.5, font=2, cex=1.1) } } ## @@ -60,14 +60,14 @@ plot_output <- function(fit, fit2=NULL, fit3=NULL, specres=NULL, cols=c("navy", ## if(2 %in% pages){ plot_compare_ageprev2(fit, fit2, fit3, specres, col=cols, likdat=likdat) - mtext(paste0(title, "; Age-specific prevalence"), 3, 0.5, outer=TRUE, font=2, cex=1.3) + graphics::mtext(paste0(title, "; Age-specific prevalence"), 3, 0.5, outer=TRUE, font=2, cex=1.3) } ## ## ## Age prevalence trend by age group ## if(3 %in% pages){ - par(oma=c(1, 1, 2.5, 1), mfrow=c(3, 2), mar=c(3.5, 3.5, 3, 1), tcl=-0.25, mgp=c(2.5, 0.5, 0), cex=1, las=1) + graphics::par(oma=c(1, 1, 2.5, 1), mfrow=c(3, 2), mar=c(3.5, 3.5, 3, 1), tcl=-0.25, mgp=c(2.5, 0.5, 0), cex=1, las=1) ## for(iagegr in 1:3){ for(isex in 1:2){ @@ -97,31 +97,31 @@ plot_output <- function(fit, fit2=NULL, fit3=NULL, specres=NULL, cols=c("navy", cred.region(xx, t(fit2.yprev[,c("lower", "upper")]), col=transp(cols[2], 0.4)) if(!is.null(fit3)) cred.region(xx, t(fit3.yprev[,c("lower", "upper")]), col=transp(cols[3], 0.4)) - lines(xx, fit.yprev[,"mean"], col=cols[1], lwd=2) + graphics::lines(xx, fit.yprev[,"mean"], col=cols[1], lwd=2) if(!is.null(fit2)) - lines(xx, fit2.yprev[,"mean"], col=cols[2], lwd=2) + graphics::lines(xx, fit2.yprev[,"mean"], col=cols[2], lwd=2) if(!is.null(fit3)) - lines(xx, fit3.yprev[,"mean"], col=cols[3], lwd=2) + graphics::lines(xx, fit3.yprev[,"mean"], col=cols[3], lwd=2) ## if(!is.null(specres)){ ages <- as.character(c(15, 25, 35)[iagegr] + 0:(c(10, 10, 15)[iagegr]-1)) specres.prev <- colSums(specres$hivpop[ages,isex,as.character(xx)]) / colSums(specres$totpop[ages,isex,as.character(xx)]) - lines(xx, specres.prev, lty=2, lwd=2, col="grey10") + graphics::lines(xx, specres.prev, lty=2, lwd=2, col="grey10") } ## if(!is.null(survdat)){ if(!is.null(survdat$prev_spec17)){ - points(survdat$year, survdat$prev, pch=15, col="grey40") - points(survdat$year, survdat$prev_spec17, pch=19) - segments(survdat$year, y0=survdat$ci_l_spec17, y1=survdat$ci_u_spec17) + graphics::points(survdat$year, survdat$prev, pch=15, col="grey40") + graphics::points(survdat$year, survdat$prev_spec17, pch=19) + graphics::segments(survdat$year, y0=survdat$ci_l_spec17, y1=survdat$ci_u_spec17) } else { - points(survdat$year, survdat$prev, pch=19) - segments(survdat$year, y0=survdat$ci_l, y1=survdat$ci_u) + graphics::points(survdat$year, survdat$prev, pch=19) + graphics::segments(survdat$year, y0=survdat$ci_l, y1=survdat$ci_u) } } } } - mtext(paste0(title, "\nPrevalence trend by age group"), 3, 0, outer=TRUE, font=2, cex=1.3) + graphics::mtext(paste0(title, "\nPrevalence trend by age group"), 3, 0, outer=TRUE, font=2, cex=1.3) } ## ## diff --git a/R/read-spectrum-files.R b/R/read-spectrum-files.R index e81fbd8..e6677c4 100644 --- a/R/read-spectrum-files.R +++ b/R/read-spectrum-files.R @@ -31,9 +31,9 @@ get_dp_version <- function(dp){ read_dp <- function(pjnz, use_ep5 = FALSE){ if(use_ep5) { - dpfile <- grep(".ep5$", unzip(pjnz, list=TRUE)$Name, value=TRUE) + dpfile <- grep(".ep5$", utils::unzip(pjnz, list=TRUE)$Name, value=TRUE) } else { - dpfile <- grep(".DP$", unzip(pjnz, list=TRUE)$Name, value=TRUE) + dpfile <- grep(".DP$", utils::unzip(pjnz, list=TRUE)$Name, value=TRUE) } dp <- vroom::vroom(unz(pjnz, dpfile), delim = ",", @@ -46,7 +46,7 @@ read_dp <- function(pjnz, use_ep5 = FALSE){ #' @export read_pjn <- function(pjnz){ - pjnfile <- grep(".PJN$", unzip(pjnz, list=TRUE)$Name, value=TRUE) + pjnfile <- grep(".PJN$", utils::unzip(pjnz, list=TRUE)$Name, value=TRUE) pjn <- vroom::vroom(unz(pjnz, pjnfile), delim = ",", col_types = vroom::cols(.default = vroom::col_character()), .name_repair = "minimal", progress = FALSE) @@ -313,10 +313,10 @@ read_hivproj_output <- function(pjnz, single.age=TRUE){ specres$infections <- array(infections, lengths(dn1), dn1) } - specres$births <- setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) - specres$hivpregwomen <- setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) - specres$hivpregwomen <- setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) - specres$receivepmtct <- setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) + specres$births <- stats::setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) + specres$hivpregwomen <- stats::setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) + specres$hivpregwomen <- stats::setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) + specres$receivepmtct <- stats::setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) class(specres) <- "specres" attr(specres, "country") <- read_country(pjnz) @@ -453,15 +453,15 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){ ## sex/age-specific incidence ratios (time varying) incrr_age <- array(NA, c(AG, NG, length(proj.years)), list(0:(AG-1)*5, c("Male", "Female"), proj.years)) if(dp.vers == ""){ - incrr_sex <- setNames(as.numeric(dp[aids5.tidx+181,timedat.idx]), proj.years) # !!! Not sure what aids5.tidx+183 (Ratio of female to male prevalence) + incrr_sex <- stats::setNames(as.numeric(dp[aids5.tidx+181,timedat.idx]), proj.years) # !!! Not sure what aids5.tidx+183 (Ratio of female to male prevalence) incrr_age[,"Male",] <- sapply(dp[aids5.tidx+281:297,timedat.idx], as.numeric) incrr_age[,"Female",] <- sapply(dp[aids5.tidx+299:315,timedat.idx], as.numeric) } else if(dp.vers == ""){ - incrr_sex <- setNames(as.numeric(dp[hivsexrat.tidx+2, timedat.idx]), proj.years) + incrr_sex <- stats::setNames(as.numeric(dp[hivsexrat.tidx+2, timedat.idx]), proj.years) incrr_age[,"Male",] <- sapply(dp[hivagedist.tidx+3:19,timedat.idx], as.numeric) incrr_age[,"Female",] <- sapply(dp[hivagedist.tidx+21:37,timedat.idx], as.numeric) } else if(dp.vers == "Spectrum2016") { - incrr_sex <- setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) + incrr_sex <- stats::setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) incrr_age[,"Male",] <- sapply(dpsub("", 4:20, timedat.idx), as.numeric) incrr_age[,"Female",] <- sapply(dpsub("", 22:38, timedat.idx), as.numeric) } else if(dp.vers == "Spectrum2017") { @@ -482,7 +482,7 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){ if (length(data_row) != 1) { stop("DP tag not parsed correctly") } - incrr_sex <- setNames(as.numeric(dpsub("", data_row, timedat.idx)), proj.years) + incrr_sex <- stats::setNames(as.numeric(dpsub("", data_row, timedat.idx)), proj.years) } else if (exists_dptag("")) { sexincrr_idx <- as.integer(dpsub("", 2, 4)) # 0-based index @@ -567,13 +567,13 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){ if(dp.vers %in% c("", "")){ art15plus_numperc <- sapply(dp[adult.artnumperc.tidx+3:4, timedat.idx], as.numeric) art15plus_num <- sapply(dp[adult.art.tidx+3:4, timedat.idx], as.numeric) - art15plus_eligthresh <- setNames(as.numeric(dp[adult.arteligthresh.tidx+2, timedat.idx]), proj.years) - artelig_specpop <- setNames(dp[specpopelig.tidx+1:7,2:6], c("description", "pop", "elig", "percent", "year")) + art15plus_eligthresh <- stats::setNames(as.numeric(dp[adult.arteligthresh.tidx+2, timedat.idx]), proj.years) + artelig_specpop <- stats::setNames(dp[specpopelig.tidx+1:7,2:6], c("description", "pop", "elig", "percent", "year")) } else if(dp.vers %in% c("Spectrum2016", "Spectrum2017")) { art15plus_numperc <- sapply(dpsub("", 4:5, timedat.idx), as.numeric) art15plus_num <- sapply(dpsub("", 4:5, timedat.idx), as.numeric) - art15plus_eligthresh <- setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) - artelig_specpop <- setNames(dpsub("", 3:9, 2:6), c("description", "pop", "elig", "percent", "year")) + art15plus_eligthresh <- stats::setNames(as.numeric(dpsub("", 2, timedat.idx)), proj.years) + artelig_specpop <- stats::setNames(dpsub("", 3:9, 2:6), c("description", "pop", "elig", "percent", "year")) } if(exists_dptag("")) @@ -633,7 +633,7 @@ read_hivproj_param <- function(pjnz, use_ep5=FALSE){ ## vertical transmission and paediatric survival - verttrans <- setNames(sapply(dpsub("", 2, timedat.idx), as.numeric)/100, proj.years) + verttrans <- stats::setNames(sapply(dpsub("", 2, timedat.idx), as.numeric)/100, proj.years) if(exists_dptag("")) hivpop <- array(sapply(dpsub("", c(3:83, 85:165), timedat.idx), as.numeric), @@ -772,7 +772,7 @@ read_demog_param <- function(upd.file, age.intervals = 1){ stop("Invalid age interval interval") age.groups <- c(rep(seq(0, 79, age.intervals), each=age.intervals), 80) age.intervals <- c(rep(age.intervals, 80 / age.intervals), 1) - ## } else if(length(age.groups) < 81 && age.groups[1] == 0 && tail(age.groups, 1) == 80){ + ## } else if(length(age.groups) < 81 && age.groups[1] == 0 && utils::tail(age.groups, 1) == 80){ ## stop("not defined yet") ## define thie one } else if(length(age.intervals) < 81){ if(sum(age.intervals != 81)) @@ -781,7 +781,7 @@ read_demog_param <- function(upd.file, age.intervals = 1){ } ## Read and parse udp file - upd <- read.csv(upd.file, header=FALSE, as.is=TRUE) + upd <- utils::read.csv(upd.file, header=FALSE, as.is=TRUE) wpp.version <- ifelse(any(upd[,1] == "RevisionYear=2015"), "wpp2015", "wpp2012") @@ -803,12 +803,12 @@ read_demog_param <- function(upd.file, age.intervals = 1){ srb.tidx <- which(upd[,1] == "") srb.eidx <- which(upd[,1] == "") - bp <- setNames(data.frame(upd[(bp.tidx+2):(bp.eidx-1),1:4]), upd[bp.tidx+1,1:4]) - lt <- setNames(data.frame(upd[(lt.tidx+2):(lt.eidx-1),]), upd[lt.tidx+1,]) - pasfrs <- setNames(data.frame(upd[(pasfrs.tidx+2):(pasfrs.eidx-1),1:3]), upd[pasfrs.tidx+1,1:3]) - migration <- setNames(data.frame(upd[(migration.tidx+2):(migration.eidx-1),1:4]), upd[migration.tidx+1,1:4]) - tfr <- setNames(as.numeric(upd[(tfr.tidx+2):(tfr.eidx-1),2]), upd[(tfr.tidx+2):(tfr.eidx-1),1]) - srb <- setNames(as.numeric(upd[(srb.tidx+2):(srb.eidx-1),2]), upd[(srb.tidx+2):(srb.eidx-1),1]) + bp <- stats::setNames(data.frame(upd[(bp.tidx+2):(bp.eidx-1),1:4]), upd[bp.tidx+1,1:4]) + lt <- stats::setNames(data.frame(upd[(lt.tidx+2):(lt.eidx-1),]), upd[lt.tidx+1,]) + pasfrs <- stats::setNames(data.frame(upd[(pasfrs.tidx+2):(pasfrs.eidx-1),1:3]), upd[pasfrs.tidx+1,1:3]) + migration <- stats::setNames(data.frame(upd[(migration.tidx+2):(migration.eidx-1),1:4]), upd[migration.tidx+1,1:4]) + tfr <- stats::setNames(as.numeric(upd[(tfr.tidx+2):(tfr.eidx-1),2]), upd[(tfr.tidx+2):(tfr.eidx-1),1]) + srb <- stats::setNames(as.numeric(upd[(srb.tidx+2):(srb.eidx-1),2]), upd[(srb.tidx+2):(srb.eidx-1),1]) ## Aggregate into specified age groups @@ -917,7 +917,7 @@ read_specdp_demog_param <- function(pjnz, use_ep5=FALSE){ tfr.tidx <- which(dp[,1] == "") asfd.tidx <- which(dp[,1] == "") - tfr <- setNames(as.numeric(dp[tfr.tidx + 2, timedat.idx]), proj.years) + tfr <- stats::setNames(as.numeric(dp[tfr.tidx + 2, timedat.idx]), proj.years) asfd <- sapply(dp[asfd.tidx + 3:9, timedat.idx], as.numeric)/100 asfd <- apply(asfd / 5, 2, rep, each=5) @@ -930,11 +930,11 @@ read_specdp_demog_param <- function(pjnz, use_ep5=FALSE){ asfr <- sweep(asfd, 2, tfr, "*") births.tidx <- which(dp[,1] == "") - births <- setNames(as.numeric(dp[births.tidx + 2, timedat.idx]), proj.years) + births <- stats::setNames(as.numeric(dp[births.tidx + 2, timedat.idx]), proj.years) ## srb srb.tidx <- which(dp[,1] == "") - srb <- setNames(as.numeric(dp[srb.tidx + 2, timedat.idx]), proj.years) + srb <- stats::setNames(as.numeric(dp[srb.tidx + 2, timedat.idx]), proj.years) ## migration if(dp.vers == "Spectrum2016"){ @@ -1030,7 +1030,7 @@ read_specdp_demog_param <- function(pjnz, use_ep5=FALSE){ #' @export read_epp_perc_urban <- function(pjnz){ - xmlfile <- grep(".xml", unzip(pjnz, list=TRUE)$Name, value=TRUE) + xmlfile <- grep(".xml", utils::unzip(pjnz, list=TRUE)$Name, value=TRUE) con <- unz(pjnz, xmlfile) epp.xml <- xml2::read_xml(con) @@ -1045,7 +1045,7 @@ read_epp_perc_urban <- function(pjnz){ yr_start <- xml2::xml_integer(r[["worksetStartYear"]]) yr_end <- xml2::xml_integer(r[["worksetEndYear"]]) - return(setNames(perc_urban, yr_start:yr_end)) + return(stats::setNames(perc_urban, yr_start:yr_end)) } ## Read epidemic start year from EPP XML file @@ -1056,7 +1056,7 @@ read_epp_perc_urban <- function(pjnz){ #' @export read_epp_t0 <- function(pjnz){ - xmlfile <- grep(".xml", unzip(pjnz, list=TRUE)$Name, value=TRUE) + xmlfile <- grep(".xml", utils::unzip(pjnz, list=TRUE)$Name, value=TRUE) con <- unz(pjnz, xmlfile) epp.xml <- xml2::read_xml(con) @@ -1096,12 +1096,12 @@ read_subp_file <- function(filepath){ datidx <- seq(startidx+1, endidx-1, AG+1) - header <- setNames(read.csv(text=dat[datidx], header=FALSE), + header <- stats::setNames(utils::read.csv(text=dat[datidx], header=FALSE), c("country_code", "country", "region", "sex")) header$sex <- factor(header$sex, c("Male", "Female")) - data <- lapply(datidx, function(idx) read.csv(text=dat[idx+1:AG], header=FALSE)) - data <- lapply(tapply(setNames(data, header$sex), header$region, abind::abind, along=0), + data <- lapply(datidx, function(idx) utils::read.csv(text=dat[idx+1:AG], header=FALSE)) + data <- lapply(tapply(stats::setNames(data, header$sex), header$region, abind::abind, along=0), aperm, c(2,1,3)) data <- lapply(data, function(x) x[,c("Male", "Female"),]) data <- lapply(data, "dimnames<-", list(Age=0:(AG-1), Sex=c("Male", "Female"), Year=startyear+1:years-1L)) @@ -1175,7 +1175,7 @@ read_incid_input <- function(pjnz){ val <- NULL } - val <- setNames(val, proj_years) + val <- stats::setNames(val, proj_years) attr(val, "incidpopage") <- as.integer(dpsub("", 2, 4)) # Adults 15-49 = 0; Adults 15+ = 1 val <- val / 100 diff --git a/R/read-spectrum-pop1.R b/R/read-spectrum-pop1.R index 955c5bb..8dc4b6e 100644 --- a/R/read-spectrum-pop1.R +++ b/R/read-spectrum-pop1.R @@ -33,7 +33,7 @@ read_pop1 <- function(pop1file, country, years = 2000:2021){ pop1 <- as.data.frame(readxl::read_excel(pop1file, sheet=as.character(yr))) pop1 <- pop1[-1,] names(pop1)[1:3] <- c("sex", "cd4", "artdur") - pop1 <- reshape(pop1, idvar=c("sex", "cd4", "artdur"), varying=c(0:79, "80+"), times=c(0:79, "80+"), timevar="age", v.names="pop", direction="long") + pop1 <- stats::reshape(pop1, idvar=c("sex", "cd4", "artdur"), varying=c(0:79, "80+"), times=c(0:79, "80+"), timevar="age", v.names="pop", direction="long") pop1$year <- as.integer(yr) pop1$age <- as.integer(gsub("\\+", "", pop1$age)) pop1$cd4 <- as.integer(gsub("CD4 Count: ", "", pop1$cd4)) diff --git a/R/rt-models.R b/R/rt-models.R index d1cae59..95eb5d8 100644 --- a/R/rt-models.R +++ b/R/rt-models.R @@ -88,7 +88,7 @@ prepare_rhybrid <- function(fp, ## Linearly interpolate between 0 and 1 over the period (rw_start, rw_start + rw_trans) ## Add a small value to avoid R error in approx() if rw_trans = 0 - rt$rw_transition <- approx(c(rw_start, rw_start + rw_trans + 0.001), c(0, 1), rt$rw_steps[-1], rule = 2)$y + rt$rw_transition <- stats::approx(c(rw_start, rw_start + rw_trans + 0.001), c(0, 1), rt$rw_steps[-1], rule = 2)$y rt$dt <- 1 / fp$ss$hiv_steps_per_year @@ -129,7 +129,7 @@ create_rvec <- function(theta, rt){ sample_invgamma_post <- function(x, prior_shape, prior_rate){ ## x: n_samples, n_knots if(is.vector(x)) x <- matrix(x, 1) - 1/rgamma(nrow(x), shape=prior_shape + ncol(x)/2, + 1/stats::rgamma(nrow(x), shape=prior_shape + ncol(x)/2, rate=prior_rate + 0.5*rowSums(x^2)) } @@ -167,7 +167,7 @@ extend_projection <- function(fit, proj_years){ if(nsteps > 0){ thetanew <- matrix(nrow=nrow(theta), ncol=fpnew$rt$n_rw) thetanew[,1:ncol(theta)] <- theta - thetanew[,ncol(theta)+1:nsteps] <- rnorm(nrow(theta)*nsteps, sd=fit$rw_sigma) + thetanew[,ncol(theta)+1:nsteps] <- stats::rnorm(nrow(theta)*nsteps, sd=fit$rw_sigma) if(idx1 > 1) fit$resample <- cbind(fit$resample[,1:(idx1-1), drop=FALSE], thetanew, fit$resample[,(idx2+1):ncol(fit$resample), drop=FALSE]) @@ -192,11 +192,11 @@ calc_rtrend_rt <- function(t, fp, rveclast, prevlast, pop, i, ii){ hivn.ii <- sum(pop[p.age15to49.idx,,hivn.idx,i]) hivn.ii <- hivn.ii - sum(pop[p.age15to49.idx[1],,hivn.idx,i])*(1-DT*(ii-1)) - hivn.ii <- hivn.ii + sum(pop[tail(p.age15to49.idx,1)+1,,hivn.idx,i])*(1-DT*(ii-1)) + hivn.ii <- hivn.ii + sum(pop[utils::tail(p.age15to49.idx,1)+1,,hivn.idx,i])*(1-DT*(ii-1)) hivp.ii <- sum(pop[p.age15to49.idx,,hivp.idx,i]) hivp.ii <- hivp.ii - sum(pop[p.age15to49.idx[1],,hivp.idx,i])*(1-DT*(ii-1)) - hivp.ii <- hivp.ii + sum(pop[tail(p.age15to49.idx,1)+1,,hivp.idx,i])*(1-DT*(ii-1)) + hivp.ii <- hivp.ii + sum(pop[utils::tail(p.age15to49.idx,1)+1,,hivp.idx,i])*(1-DT*(ii-1)) prevcurr <- hivp.ii / (hivn.ii + hivp.ii) @@ -247,7 +247,7 @@ lprior_iota <- function(par, fp){ if(exists("logitiota", fp) && fp$logitiota) ldinvlogit(par) # Note: parameter is defined on range logiota.unif.prior, so no need to check bound else - dunif(par, logiota.unif.prior[1], logiota.unif.prior[2], log=TRUE) + stats::dunif(par, logiota.unif.prior[1], logiota.unif.prior[2], log=TRUE) } sample_iota <- function(n, fp){ @@ -256,9 +256,9 @@ sample_iota <- function(n, fp){ assign(names(fp$prior_args)[i], fp$prior_args[[i]]) } if(exists("logitiota", fp) && fp$logitiota) - return(logit(runif(n))) + return(logit(stats::runif(n))) else - runif(n, logiota.unif.prior[1], logiota.unif.prior[2]) + stats::runif(n, logiota.unif.prior[1], logiota.unif.prior[2]) } ldsamp_iota <- lprior_iota diff --git a/R/spectrum.R b/R/spectrum.R index 1edbc0f..c37e3f6 100644 --- a/R/spectrum.R +++ b/R/spectrum.R @@ -238,7 +238,7 @@ create_spectrum_fixpar <- function(projp, demp, hiv_steps_per_year = 10L, proj_s ## Vertical transmission and survival to AGE_START for lagged births - fp$verttrans_lag <- setNames(c(rep(0, AGE_START), projp$verttrans[1:(PROJ_YEARS-AGE_START)]), proj_start:proj_end) + fp$verttrans_lag <- stats::setNames(c(rep(0, AGE_START), projp$verttrans[1:(PROJ_YEARS-AGE_START)]), proj_start:proj_end) ## calculate probability of HIV death in each year hivqx <- apply(projp$hivdeaths[1:AGE_START,,], c(1,3), sum) / apply(projp$hivpop[1:AGE_START,,], c(1,3), sum) @@ -247,7 +247,7 @@ create_spectrum_fixpar <- function(projp, demp, hiv_steps_per_year = 10L, proj_s ## probability of surviving to AGE_START for each cohort (product along diagonal) cumhivsurv <- sapply(1:(PROJ_YEARS - AGE_START), function(i) prod(1-hivqx[cbind(1:15, i-1+1:15)])) - fp$paedsurv_lag <- setNames(c(rep(1, AGE_START), cumhivsurv), proj_start:proj_end) + fp$paedsurv_lag <- stats::setNames(c(rep(1, AGE_START), cumhivsurv), proj_start:proj_end) ## ## EQUIVALENT CODE, easier to read ## fp$paedsurv_lag <- rep(1.0, PROJ_YEARS) diff --git a/R/tidy-outputs.R b/R/tidy-outputs.R index 88ac325..9d5688e 100644 --- a/R/tidy-outputs.R +++ b/R/tidy-outputs.R @@ -26,7 +26,7 @@ tidy_output <- function(fit, modlab, country=NA, eppregion=NA, ancsite = TRUE){ if(fit$fp$eppmod == "rspline") fit <- rw_projection(fit) - fp_list <- lapply(param_list, function(par) update(fit$fp, list=par)) + fp_list <- lapply(param_list, function(par) stats::update(fit$fp, list=par)) mod_list <- lapply(fp_list, simmod) @@ -214,7 +214,7 @@ plot_output <- function(out, data_color = "grey20", th = theme_light()){ out$pregprev$indicator <- "Prevalence among pregnant women" df <- rbind(subset(out$core, year %in% 1990:2018), subset(out$pregprev, agegr == "15-49", - c(country, eppregion, year, indicator, model, mean, se, median, lower, upper))) + c(country, eppregion, year, indicator, model, mean, se, stats::median, lower, upper))) obs <- data.frame(indicator = "Prevalence 15-49y", subset(prev_15to49_eppregion, @@ -257,7 +257,7 @@ plot_output <- function(out, data_color = "grey20", th = theme_light()){ ## ANC prevalence df <- out$pregprev - df$agegr <- relevel(factor(df$agegr), "15-49") + df$agegr <- stats::relevel(factor(df$agegr), "15-49") obs <- subset(ancrtcens, country == out$country & eppregion == out$eppregion) if(nrow(obs)) @@ -265,10 +265,10 @@ plot_output <- function(out, data_color = "grey20", th = theme_light()){ if(nrow(obs)) obs$agegr <- factor(obs$agegr, levels(df$agegr)) x.ancrt <- (obs$prev*obs$n+0.5)/(obs$n+1) - obs$W.ancrt <- qnorm(x.ancrt) + obs$W.ancrt <- stats::qnorm(x.ancrt) obs$v.ancrt <- 2*pi*exp(obs$W.ancrt^2)*x.ancrt*(1-x.ancrt)/obs$n - obs$ci_l <- with(obs, pnorm(W.ancrt - qnorm(0.975) * sqrt(v.ancrt))) - obs$ci_u <- with(obs, pnorm(W.ancrt + qnorm(0.975) * sqrt(v.ancrt))) + obs$ci_l <- with(obs, stats::pnorm(W.ancrt - stats::qnorm(0.975) * sqrt(v.ancrt))) + obs$ci_u <- with(obs, stats::pnorm(W.ancrt + stats::qnorm(0.975) * sqrt(v.ancrt))) ancsiteobs <- subset(ancsitedata, country == out$country & eppregion == out$eppregion & @@ -303,7 +303,7 @@ get_pointwise_ll <- function(fit, newdata = fit$likdat){ if(fit$fp$eppmod == "rhybrid") fit <- extend_projection(fit, proj_years = fit$fp$ss$PROJ_YEARS) - fp_list <- lapply(param_list, function(par) update(fit$fp, list=par)) + fp_list <- lapply(param_list, function(par) stats::update(fit$fp, list=par)) mod_list <- lapply(fp_list, simmod) diff --git a/data-raw/spectrum-country-list.R b/data-raw/spectrum-country-list.R index 39f9df7..d53231d 100644 --- a/data-raw/spectrum-country-list.R +++ b/data-raw/spectrum-country-list.R @@ -4,10 +4,10 @@ load("../R/sysdata.rda") -spectrum5_countrylist <- read.csv("CountryListMaster.csv", as.is=TRUE, encoding = "UTF-8") +spectrum5_countrylist <- utils::read.csv("CountryListMaster.csv", as.is=TRUE, encoding = "UTF-8") spectrum5_countrylist$Country[spectrum5_countrylist$Code == 384] <- "Côte d'Ivoire" -iso <- read.csv("iso.csv", na.strings = "", as.is=TRUE) +iso <- utils::read.csv("iso.csv", na.strings = "", as.is=TRUE) names(iso) <- sub("alpha.2", "iso2", names(iso)) names(iso) <- sub("alpha.3", "iso3", names(iso)) diff --git a/data-raw/subp-gfr-estimates.R b/data-raw/subp-gfr-estimates.R index b37f21f..f853064 100644 --- a/data-raw/subp-gfr-estimates.R +++ b/data-raw/subp-gfr-estimates.R @@ -119,8 +119,8 @@ bh <- merge(bh, wm[c("HH1", "HH2", "LN", "start", "end")], by=c("HH1", "HH2", "L bh <- subset(bh, BH4C > start & BH4C <= end) bh$births_wt <- bh$wmweight -gfr <- merge(aggregate(births_wt ~ HH6, bh, sum), - aggregate(pys_wt ~ HH6, wm, sum)) +gfr <- merge(stats::aggregate(births_wt ~ HH6, bh, sum), + stats::aggregate(pys_wt ~ HH6, wm, sum)) gfr$gfr <- with(gfr, births_wt / pys_wt) subset(subp_gfr, country == "South Sudan") @@ -153,8 +153,8 @@ bh <- merge(bh, wm[c("HH1", "HH2", "LN", "start", "end")], by=c("HH1", "HH2", "L bh <- subset(bh, BH4C > start & BH4C <= end) bh$births_wt <- bh$wmweight -gfr <- merge(aggregate(births_wt ~ HH6, bh, sum), - aggregate(pys_wt ~ HH6, wm, sum)) +gfr <- merge(stats::aggregate(births_wt ~ HH6, bh, sum), + stats::aggregate(pys_wt ~ HH6, wm, sum)) gfr$gfr <- with(gfr, births_wt / pys_wt) subset(subp_gfr, country == "Guinea-Bissau") @@ -182,8 +182,8 @@ bh <- merge(bh, wm[c("HH1", "HH2", "LN", "start", "end")], by=c("HH1", "HH2", "L bh <- subset(bh, BH4C > start & BH4C <= end) bh$births_wt <- bh$WMWEIGHT -gfr <- merge(aggregate(births_wt ~ HH6, bh, sum), - aggregate(pys_wt ~ HH6, wm, sum)) +gfr <- merge(stats::aggregate(births_wt ~ HH6, bh, sum), + stats::aggregate(pys_wt ~ HH6, wm, sum)) gfr$gfr <- with(gfr, births_wt / pys_wt) sum(gfr$births_wt) / sum(gfr$pys_wt) @@ -212,8 +212,8 @@ bh <- merge(bh, wm[c("HH1", "HH2", "LN", "start", "end")], by=c("HH1", "HH2", "L bh <- subset(bh, BH4C > start & BH4C <= end) bh$births_wt <- bh$wmweight -gfr <- merge(aggregate(births_wt ~ HH6, bh, sum), - aggregate(pys_wt ~ HH6, wm, sum)) +gfr <- merge(stats::aggregate(births_wt ~ HH6, bh, sum), + stats::aggregate(pys_wt ~ HH6, wm, sum)) gfr$gfr <- with(gfr, births_wt / pys_wt) gfr sum(gfr$births_wt) / sum(gfr$pys_wt) @@ -248,8 +248,8 @@ bh$births_wt <- bh$wmweight bh$HH6[bh$HH6=="Nomadic"] <- "Rural" -gfr <- merge(aggregate(births_wt ~ HH7, bh, sum), - aggregate(pys_wt ~ HH7, wm, sum)) +gfr <- merge(stats::aggregate(births_wt ~ HH7, bh, sum), + stats::aggregate(pys_wt ~ HH7, wm, sum)) gfr$gfr <- with(gfr, births_wt / pys_wt) gfr sum(gfr$births_wt) / sum(gfr$pys_wt) @@ -280,8 +280,8 @@ bh$births_wt <- bh$wmweight bh$HH6[bh$HH6=="Nomadic"] <- "Rural" -gfr <- merge(aggregate(births_wt ~ HH7, bh, sum), - aggregate(pys_wt ~ HH7, wm, sum)) +gfr <- merge(stats::aggregate(births_wt ~ HH7, bh, sum), + stats::aggregate(pys_wt ~ HH7, wm, sum)) gfr$gfr <- with(gfr, births_wt / pys_wt) gfr sum(gfr$births_wt) / sum(gfr$pys_wt) @@ -317,8 +317,8 @@ bh$births_wt <- bh$wmweight bh$residence[bh$residence == "City/Town"] <- "Urban Village" wm$residence[wm$residence == "City/Town"] <- "Urban Village" -gfr <- merge(aggregate(births_wt ~ residence, bh, sum), - aggregate(pys_wt ~ residence, wm, sum)) +gfr <- merge(stats::aggregate(births_wt ~ residence, bh, sum), + stats::aggregate(pys_wt ~ residence, wm, sum)) gfr$gfr <- with(gfr, births_wt / pys_wt) gfr sum(gfr$births_wt) / sum(gfr$pys_wt) diff --git a/data-raw/test-pjnz.R b/data-raw/test-pjnz.R index 6f9358e..f02a389 100644 --- a/data-raw/test-pjnz.R +++ b/data-raw/test-pjnz.R @@ -17,12 +17,12 @@ file.copy("~/Documents/Data/Spectrum files/2017 final/2017 Final Spectrum files/ "~/Documents/Code/R/eppasm/inst/extdata/testpjnz/Botswana2017.PJNZ") fp <- "~/Documents/Code/R/eppasm/inst/extdata/testpjnz/Botswana2017.PJNZ" -unzip(fp, list=TRUE) +utils::unzip(fp, list=TRUE) ## Remove files to reduce size -zip(fp, grep("bm2$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") -zip(fp, grep("SPU$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") -zip(fp, grep("UAC$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("bm2$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("SPU$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("UAC$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") #' ## Dominican Republic 2017 file @@ -31,12 +31,12 @@ file <- "~/Documents/Code/R/eppasm/inst/extdata/testpjnz/DominicanRepublic2017.P file.copy("~/Documents/Data/Spectrum files/2017 final/CAR/Dominican Republic_2017_final.PJNZ", file) -unzip(fp, list=TRUE) +utils::unzip(fp, list=TRUE) ## Remove files to reduce size -zip(file, grep("bm2$", unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") -zip(file, grep("SPU$", unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") -zip(file, grep("UAC$", unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") +zip(file, grep("bm2$", utils::unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") +zip(file, grep("SPU$", utils::unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") +zip(file, grep("UAC$", utils::unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") #' ## Dominican Republic 2017 file @@ -45,12 +45,12 @@ file <- "~/Documents/Code/R/eppasm/inst/extdata/testpjnz/Netherlands2017.PJNZ" file.copy("~/Documents/Data/Spectrum files/2017 final/WCENA/Netherlands_2017_final.PJNZ", file) -unzip(fp, list=TRUE) +utils::unzip(fp, list=TRUE) ## Remove files to reduce size -zip(file, grep("bm2$", unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") -zip(file, grep("SPU$", unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") -zip(file, grep("UAC$", unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") +zip(file, grep("bm2$", utils::unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") +zip(file, grep("SPU$", utils::unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") +zip(file, grep("UAC$", utils::unzip(file, list=TRUE)$Name, value=TRUE), flags="-d") @@ -59,13 +59,13 @@ file.copy("~/Documents/Data/Spectrum files/2018 final/SSA/Botswana_ 2018 updated "~/Documents/Code/R/eppasm/inst/extdata/testpjnz/Botswana2018.PJNZ") fp <- "~/Documents/Code/R/eppasm/inst/extdata/testpjnz/Botswana2018.PJNZ" -unzip(fp, list=TRUE) +utils::unzip(fp, list=TRUE) ## Remove files to reduce size -zip(fp, grep("bm2$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") -zip(fp, grep("SPU$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") -zip(fp, grep("UAC$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") -zip(fp, grep("DPUAD_AUA$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("bm2$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("SPU$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("UAC$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("DPUAD_AUA$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") ## Mozambique Maputo Cidade 2018 file @@ -73,11 +73,11 @@ file.copy("~/Documents/Data/Spectrum files/2018 final/SSA/11_MZ_Maputo Cidade_v5 here::here("inst/extdata/testpjnz/Mozambique_Maputo_Cidade2018.PJNZ")) fp <- here::here("inst/extdata/testpjnz/Mozambique_Maputo_Cidade2018.PJNZ") -unzip(fp, list=TRUE) +utils::unzip(fp, list=TRUE) ## Remove files to reduce size -zip(fp, grep("bm2$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") -zip(fp, grep("SPU$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") -zip(fp, grep("UAC$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") -zip(fp, grep("DPUAD_AUA$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") -zip(fp, grep("ep5.bak$", unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("bm2$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("SPU$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("UAC$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("DPUAD_AUA$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") +zip(fp, grep("ep5.bak$", utils::unzip(fp, list=TRUE)$Name, value=TRUE), flags="-d") diff --git a/dev/dev-ageincid.R b/dev/dev-ageincid.R index 9316aa3..8b07f8c 100644 --- a/dev/dev-ageincid.R +++ b/dev/dev-ageincid.R @@ -13,7 +13,7 @@ fp <- create_spectrum_fixpar(projp, demp, proj_end = 2016, time_epi_start = proj theta.rspline <- c(1.31820011, -0.09884313, -0.40054248, 0.06277183, 0.16923859, 0.41277390, -0.17640756, -14.13863910, 0.09765759, -3.73232668, -5.12046650) param <- fnCreateParam(theta.rspline, fp) -fp <- update(fp, list=param) +fp <- stats::update(fp, list=param) fp$iota <- fp$iota*15 modR <- simmod(fp, VERSION="R") diff --git a/man/create_spectrum_fixpar.Rd b/man/create_spectrum_fixpar.Rd index 8f5b4ea..6a48b36 100644 --- a/man/create_spectrum_fixpar.Rd +++ b/man/create_spectrum_fixpar.Rd @@ -19,7 +19,8 @@ create_spectrum_fixpar( who34percelig = 0, frr_art6mos = projp$frr_art6mos, frr_art1yr = projp$frr_art6mos, - projection_period = NULL + projection_period = NULL, + art_dropout_recover_cd4 = NULL ) } \description{ diff --git a/tests/testthat/test_ll.R b/tests/testthat/test_ll.R index fb0a5b9..86c4644 100644 --- a/tests/testthat/test_ll.R +++ b/tests/testthat/test_ll.R @@ -24,7 +24,7 @@ theta_rspline <- c(0.753890895839314, 0.700227807059524, -1.92682607738255, 1.12 1.97574135543477, 0.24367501614327, -3.80194029103096, -0.28992227570734, -4.03773908970723) -fp_rspline_eq <- update(fp_rspline, equil.rprior = TRUE) +fp_rspline_eq <- stats::update(fp_rspline, equil.rprior = TRUE) fp_rtrend <- prepare_rtrend_model(fp) theta_rtrend <- c(1977.68776203971, 16.2122044236609, 0.196472777334891, 0.455548771721825, diff --git a/tests/testthat/test_simmod.R b/tests/testthat/test_simmod.R index 2d0b1fe..86aab77 100644 --- a/tests/testthat/test_simmod.R +++ b/tests/testthat/test_simmod.R @@ -17,7 +17,7 @@ bw_theta <- c(-0.407503322169364, -2.76794181367538, -1.26018073624346, 1995.964 -6.10051517060137) param <- fnCreateParam(bw_theta, bw_fp) -bw_fp <- update(bw_fp, list=param) +bw_fp <- stats::update(bw_fp, list=param) bw_prev_mod <- c(0.00045, 0.00080, 0.00140, 0.00245, 0.00424, 0.00725, 0.01214, 0.01985, 0.03147, 0.04804, 0.07013, 0.09750, 0.12857, 0.16083, @@ -49,7 +49,7 @@ theta <- c(-0.407503322169364, -2.76794181367538, -1.26018073624346, 1995.964477 -6.10051517060137) param <- fnCreateParam(theta, mp_fp) -mp_fp <- update(mp_fp, list=param) +mp_fp <- stats::update(mp_fp, list=param) mp_prev_mod <- c(0.00049, 0.00087, 0.00154, 0.00271, 0.00468, 0.00792, 0.01299, 0.0207, 0.03197, 0.04795, 0.0696, 0.0971, 0.12906, 0.16258, 0.19452, diff --git a/vignettes_src/eppasm-basic-example.Rmd b/vignettes_src/eppasm-basic-example.Rmd index d041676..f98bf3d 100644 --- a/vignettes_src/eppasm-basic-example.Rmd +++ b/vignettes_src/eppasm-basic-example.Rmd @@ -63,7 +63,7 @@ fp$rw_start <- 2005 fp$incidmod <- "eppspectrum" param <- fnCreateParam(theta_ur, fp) -fp_par <- update(fp, list = param) +fp_par <- stats::update(fp, list = param) ``` Simulate the model once. diff --git a/vignettes_src/eppasm-spectrum-workflow.Rmd b/vignettes_src/eppasm-spectrum-workflow.Rmd index 7120545..4ef0624 100644 --- a/vignettes_src/eppasm-spectrum-workflow.Rmd +++ b/vignettes_src/eppasm-spectrum-workflow.Rmd @@ -120,7 +120,7 @@ outputs_age$outcome <- factor(outputs_age$outcome, c("totpop", "hivpop", "infect dcast(outputs_age[sampleid == 1], year + sex + age ~ outcome) %>% write.csv("bwaggr_ageoutputs_single-resample_v2.csv", row.names=FALSE) -dcast(outputs_age[ , .(value = median(value)), .(year, sex, age, outcome)], +dcast(outputs_age[ , .(value = stats::median(value)), .(year, sex, age, outcome)], year + sex + age ~ outcome) %>% write.csv("bwaggr_ageoutputs_median_v2.csv", row.names=FALSE) @@ -161,7 +161,7 @@ outputs_aggr$outcome <- factor(outputs_aggr$outcome, dcast(outputs_aggr[sampleid == 1], year ~ outcome) %>% write.csv("bwaggr_aggroutputs_single-resample_v2.csv", row.names=FALSE) -dcast(outputs_aggr[ , .(value = median(value)), .(year, outcome)], +dcast(outputs_aggr[ , .(value = stats::median(value)), .(year, outcome)], year ~ outcome) %>% write.csv("bwaggr_aggroutputs_median_v2.csv", row.names=FALSE) @@ -200,7 +200,7 @@ bwfit <- fitmod(bw, bwfit <- extend_projection(bwfit, 52) bwfit$param <- lapply(seq_len(nrow(bwfit$resample)), function(ii) fnCreateParam(bwfit$resample[ii,], bwfit$fp)) -fplist <- lapply(bwfit$param, function(par) update(bwfit$fp, list=par)) +fplist <- lapply(bwfit$param, function(par) stats::update(bwfit$fp, list=par)) bwnat <- lapply(fplist, simmod) bwnat <- sim_mod_list(bwfit) @@ -245,7 +245,7 @@ outputs_age$outcome <- factor(outputs_age$outcome, c("totpop", "hivpop", "infect dcast(outputs_age[sampleid == 1], year + sex + age ~ outcome) %>% write.csv("bwnat_ageoutputs_single-resample_v2.csv", row.names=FALSE) -dcast(outputs_age[ , .(value = median(value)), .(year, sex, age, outcome)], +dcast(outputs_age[ , .(value = stats::median(value)), .(year, sex, age, outcome)], year + sex + age ~ outcome) %>% write.csv("bwnat_ageoutputs_median_v2.csv", row.names=FALSE) @@ -278,7 +278,7 @@ outputs_aggr$outcome <- factor(outputs_aggr$outcome, dcast(outputs_aggr[sampleid == 1], year ~ outcome) %>% write.csv("bwnat_aggroutputs_single-resample_v2.csv", row.names=FALSE) -dcast(outputs_aggr[ , .(value = median(value)), .(year, outcome)], +dcast(outputs_aggr[ , .(value = stats::median(value)), .(year, outcome)], year ~ outcome) %>% write.csv("bwnat_aggroutputs_median_v2.csv", row.names=FALSE)