Skip to content

Commit

Permalink
Merge pull request #42 from mrc-ide/import-unzip
Browse files Browse the repository at this point in the history
Qualify package names
  • Loading branch information
r-ash authored Mar 18, 2024
2 parents 7f62e30 + 4cc870a commit fd37b5e
Show file tree
Hide file tree
Showing 30 changed files with 502 additions and 504 deletions.
8 changes: 0 additions & 8 deletions .travis.yml

This file was deleted.

4 changes: 3 additions & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -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 = "[email protected]", role = c("aut", "cre"))
Description: What the package does (one paragraph).
Depends: R (>= 3.1.0),
Expand All @@ -11,6 +11,8 @@ Imports:
epp (>= 0.4.1),
fastmatch (>= 1.1),
mvtnorm,
plyr,
readxl,
reshape2,
vroom,
xml2
Expand Down
4 changes: 4 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
@@ -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.
Expand Down
2 changes: 1 addition & 1 deletion R/aggregate-fits.R
Original file line number Diff line number Diff line change
Expand Up @@ -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")]
}
Expand Down
6 changes: 3 additions & 3 deletions R/cluster-functions.R
Original file line number Diff line number Diff line change
@@ -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]
26 changes: 13 additions & 13 deletions R/fit-model.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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")]
}
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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 {
Expand Down Expand Up @@ -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")
Expand Down Expand Up @@ -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
Expand Down
20 changes: 10 additions & 10 deletions R/imis.R
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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
Expand Down
9 changes: 4 additions & 5 deletions R/infections.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -79,4 +79,3 @@ create_beers <- function(n5yr){

## Beers coefficient matrix
beers_Amat <- create_beers(17)[16:81, 4:17]

Loading

0 comments on commit fd37b5e

Please sign in to comment.