fix age15plus indexing

DESCRIPTION

@@ -1,10 +1,10 @@
-Package: eppspectrum
+Package: eppasm
 Title: Age-structured EPP Model for HIV Epidemic Estimates
-Version: 0.1.1
+Version: 0.4.0
 Authors@R: person("Jeff", "Eaton", email = "[email protected]", role = c("aut", "cre"))
 Description: What the package does (one paragraph).
-Depends: R (>= 3.1.0), epp (>= 0.3), 
-Imports: fastmatch (>= 1.1), abind (>= 1.4)
+Depends: R (>= 3.1.0),
+Imports: epp (>= 0.3.3), fastmatch (>= 1.1), abind (>= 1.4), mgcv (>= 1.8), mvtnorm
 LinkingTo: BH
 License: GPL-3
 Encoding: UTF-8

NAMESPACE

@@ -1,3 +1,3 @@
 # Generated by roxygen2: fake comment so roxygen2 overwrites silently.
 exportPattern("^[^\\.]")
-useDynLib(eppspectrum, checkBoostAsserts, spectrumC)
+useDynLib(eppasm, checkBoostAsserts, spectrumC)

NEWS.md

@@ -2,7 +2,70 @@
 
 - Update to interface with `epp v0.3` which implements ANC routine testing (ANC-RT) likelihood in EPP 2017.
 
+
 # eppspectrum 0.1.1
 
 - Add option to use .ep5 file inputs instead of .DP file when preparing fit from PJNZ file. (use_ep5=FALSE by default.)
-- update prepare_spec_fit() to not use popadjust if EPP projection has a single population.
+- update prepare_spec_fit() to not use popadjust if EPP projection has a single population.
+
+# eppspectrum 0.1.3
+
+- EPP model setup occurs in fitmod()
+- Allow order of spline model penalty to vary (fp$rtpenord)
+- Allow number of knots to vary (fp$numKnots)
+- Return spline coefficients in fnCreateParam()
+
+# eppspectrum 0.1.5
+
+- Implemented O-spline model (penalized B-spline) using mgcv to construct
+- Added `r0logiotaratio` option to estimate ratio of r(t0) and log of initial seed, reduces correation. Not set as default.
+
+# eppspectrum 0.1.6
+
+- Bug fix: add check in ll_hhage_binom() that ldbinom() does not return NaN.
+- Bug fix: in create_spectrum_fixpar(), handle case where ART initiation starts in different years for men and women.
+
+
+# eppasm 0.3.0
+
+- Rename package `eppasm`.
+- Add IMIS function to eppasm package. In this version of IMIS, at each iteration a new mixture component is constructed either centered on input with greatest weight or based on optimizer. Allows optimizers to be run at arbitrary iterations, but doesn't implment multiple optimizers which might still be usful. Covariance for optimizer mixture component when Hessian is degenerate still needs further work.
+
+- Change from dependency to imports epp package.
+- Add ANC-RT data to aggregate national fits (prepare_national_fit() and create_aggr_input()).
+
+- Implement dsamp() to calculate density for initial IMIS samples. This allows different initial sampling density from prior distribution.
+
+- Specify knot locations for o-spline model
+- Estimate HIV+:HIV- FRR adjustment.
+- Specify knot locations for o-spline model
+
+- Add functions to read country and subnational region from Spectrum PJN file.
+
+# eppasm 0.3.2
+
+- Integrate random effects variance parameter out of likelihood.
+- Implement hybrid r-spline and RW model for r(t).
+- Implement logistic model for log r(t) (rlogistic), and version with random walk (rlogistic_rw).
+- Add optimization option to fitmod().
+
+- Add option to fit without site-level ANC data.
+
+# eppasm 0.3.3
+
+- Modularize code for calculating new infections.
+
+- Implement age 15 entrant prevalence and ART coverage inputs.
+- Implement logit transformation for sampling from uniform logiota prior distribution.
+- Set age IRR penalty variance at a fixed constant value (9003)
+
+- Implement piecewise-linear time-varying M:F incidence/transmission ratio
+  and 15-24 incidence rate ratio. Knots fixed at 2002, 2007, and 2012.
+- Modularize code for incidence rate ratio parameterization, now in infections.R.
+
+- Updated default IRR priors to informative priors based on fits countries with several surveys.
+- Add function to generate summary model outputs
+
+
+# eppasm 0.4.0
+- Add annual direct incidence input option.

R/aggregate-fits.R

@@ -0,0 +1,59 @@
+create_aggr_input <- function(inputlist){
+
+  val <- inputlist[[1]]
+
+  eppdlist <- lapply(inputlist, attr, "eppd")
+
+  ## Only keep intersecting years
+  anc.prev.list <- lapply(eppdlist, "[[", "anc.prev")
+  anc.prev.list <- lapply(anc.prev.list, function(x) x[,Reduce(intersect, lapply(anc.prev.list, colnames))])
+
+  anc.n.list <- lapply(eppdlist, "[[", "anc.n")
+  anc.n.list <- lapply(anc.n.list, function(x) x[,Reduce(intersect, lapply(anc.n.list, colnames))])
+
+  ancrtsite.prev.list <- lapply(eppdlist, "[[", "ancrtsite.prev")
+  ancrtsite.prev.list <- lapply(ancrtsite.prev.list, function(x) x[,Reduce(intersect, lapply(ancrtsite.prev.list, colnames))])
+
+  ancrtsite.n.list <- lapply(eppdlist, "[[", "ancrtsite.n")
+  ancrtsite.n.list <- lapply(ancrtsite.n.list, function(x) x[,Reduce(intersect, lapply(ancrtsite.n.list, colnames))])
+
+
+  ## aggregate census data across regions
+  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$prev <- ancrtcens$x / ancrtcens$n
+    ancrtcens <- ancrtcens[c("year", "prev", "n")]
+  }
+
+  attr(val, "eppd") <- list(anc.used = do.call(c, lapply(eppdlist, "[[", "anc.used")),
+                            anc.prev = do.call(rbind, anc.prev.list),
+                            anc.n = do.call(rbind, anc.n.list),
+                            ancrtsite.prev = do.call(rbind, ancrtsite.prev.list),
+                            ancrtsite.n = do.call(rbind, ancrtsite.n.list),
+                            ancrtcens = ancrtcens)
+
+  artnumperc <- !attr(inputlist[[1]], "specfp")$art15plus_isperc
+  artnumlist <- lapply(lapply(inputlist, attr, "specfp"), "[[", "art15plus_num")
+
+  art15plus_num <- artnumlist[[1]]
+  art15plus_num[artnumperc] <- Reduce("+", lapply(artnumlist, "[", artnumperc))
+
+  attr(val, "specfp")$art15plus_num <- art15plus_num
+
+  return(val)
+}
+
+
+fnAddHHSLikDat <- function(obj){
+  objcountry <- attr(obj, "country")
+  fp <- attr(obj, "specfp")
+  anchor.year <- as.integer(floor(min(fp$proj.steps)))
+
+  attr(obj, "eppd")$hhs <- subset(prev.15to49.nat, country==objcountry)
+  attr(obj, "eppd")$hhsage <- subset(prev.agesex.nat, country==objcountry)
+  ## attr(obj, "eppd")$sibmx <- subset(sib.mx.tips, country==objcountry)
+
+  return(obj)
+}

R/cluster-functions.R

@@ -0,0 +1,5 @@
+
+pick_maxlpd <- function(set){set <- set[!sapply(set, is.null)]; set[[which.max(sapply(set, function(x) 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]

R/directincid.R

@@ -0,0 +1,48 @@
+calc_entrantprev <- function(specres){
+  ageprev1 <- function(specres, str){colSums(specres$hivpop[str,,,drop=FALSE],,2) / colSums(specres$totpop[str,,,drop=FALSE],,2)}
+  prev14 <- ageprev1(specres, "14")
+  prev13 <- ageprev1(specres, "13")
+  Sx <- c(0, prev14[-1] / prev13[-length(prev13)])  # assume that survival from age 14 to 15 is the same as 13 to 14
+  Sx <- ifelse(is.na(Sx), 0, Sx)
+  return(Sx*prev14)
+}
+
+calc_entrantartcov <- function(specres){
+  ## Assume uniform ART coverage among age 10-14
+  artcov <- with(specres, ((artnum.f+artnum.m)/(hivnum.f+hivnum.m))["10-15",])
+  artcov <- ifelse(is.nan(artcov), 0, artcov)
+  artcov
+}
+
+create_specfp <- function(pjnz, upd.path=NULL, hiv_steps_per_year = 10L, use_ep5=FALSE){
+
+  ## demographic inputs
+  if(!is.null(upd.path))
+    demp <- read_demog_param(upd.path) # from UPD file
+  else
+    demp <- read_specdp_demog_param(pjnz, use_ep5=use_ep5)  # from DP file or ep5 file
+
+  ## HIV projection parameters (most AIM adult parameters)
+  projp <- read_hivproj_param(pjnz, use_ep5=use_ep5)
+
+  ## HIV projection results
+  specres <- read_hivproj_output(pjnz)
+
+  specfp <- create_spectrum_fixpar(projp, demp, hiv_steps_per_year= hiv_steps_per_year)
+  specfp$entrantprev <- calc_entrantprev(specres)
+  specfp$entrantartcov <- calc_entrantartcov(specres)
+
+  attr(specfp, "country") <- read_country(pjnz)
+  attr(specfp, "country") <- read_country(pjnz)
+  attr(specfp, "region") <- read_region(pjnz)
+
+  return(specfp)
+}
+
+prepare_directincid <- function(pjnz){
+  specfp <- create_specfp(pjnz)
+  specfp$eppmod <- "directincid"
+  specfp$incidinput <- read_incid_input(pjnz)
+  specfp$incidpopage <- attr(specfp$incidinput, "incidpopage")
+  return(specfp)
+}