R/plot_agefit.R

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

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){
    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))
      plot_prev(fitaggr, fitincrr, fit3, col=cols, plotprevdat=is.null(age15to49dat))
    else
      plot_prev(fitaggr, fitincrr, col=cols, plotprevdat=is.null(age15to49dat))
    if(!is.null(specres))
      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)){
          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 {
          graphics::points(age15to49dat$year, age15to49dat$prev, pch=19)
          graphics::segments(age15to49dat$year, y0=age15to49dat$ci_l, y1=age15to49dat$ci_u)
        }
    }
    graphics::mtext("HIV prevalence, age 15-49y", 3, 0.5, font=2, cex=1.2)
    ##
    ## incidence trend
    if(!is.null(fit3))
      plot_incid(fitaggr, fitincrr, fit3, col=cols)
    else
      plot_incid(fitaggr, fitincrr, col=cols)
    if(!is.null(specres))
      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]])){
      if(!is.null(fit3))
        plot_log_rvec(fitaggr, fitincrr, fit3, col=cols)
      else
        plot_log_rvec(fitaggr, fitincrr, col=cols)
      graphics::mtext("log r(t)", 3, 0.5, font=2, cex=1.2)
    }
    ##
    ## vinfl
    ##
    ## 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 <- 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 <- 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")
      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)){
        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)
      }
      graphics::segments(x0=mean(log(sapply(fitincrr$param, "[[", "incrr_sex")[yidx,])), y0=0, y1=6, col=cols[2], lwd=2)
      if(!is.null(fit3))
        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
      xx <- c(1:2, 4:7)
      logincrrage <- estci(sapply(fitincrr$param, "[[", "logincrr_age"))
      if(!is.null(fit3))
        logincrrage3 <- estci(sapply(fit3$param, "[[", "logincrr_age"))
      ## 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")
      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]+stats::qnorm(0.975)*eppasm:::ageincrr.pr.sd,
           col=transp("darkblue", 0.1), border=transp("darkblue", 0.5), lty=2)
      ##
      graphics::rect(xx+0.1,  logincrrage[xx,3], xx+0.9, logincrrage[xx,4], col=transp(cols[2]), border=NA)
      if(!is.null(fit3))
        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])
      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))
        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")
      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]+stats::qnorm(0.975)*eppasm:::ageincrr.pr.sd,
           col=transp("darkblue", 0.1), border=transp("darkblue", 0.5), lty=2)
      ##
      offset <- nrow(logincrrage)/2
      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
        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])
      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))
        graphics::segments(xx+0.1, logincrrage3[xx+offset3,1], xx+0.9, col=cols[3], lwd=2)
      ##
    }
    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)
    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){
    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){
        strsex <- c("male", "female")[isex]
        stragegr <- c("15-24", "25-34", "35-49")[iagegr]
        aggr.yprev <- estci(fitaggr$agegr3prev[iagegr,isex,,])
        fit.yprev <- estci(fitincrr$agegr3prev[iagegr,isex,,])
        if(!is.null(fit3))
          fit3.yprev <- estci(fit3$agegr3prev[iagegr,isex,,])
        survdat <- subset(agegr3dat, sex == strsex & agegr3==stragegr)
        ##
        xx <- 1998+seq_len(nrow(fit.yprev))
        plot(xx, fit.yprev[,1], type="n", ylim=c(0, 0.05*ceiling(max(survdat$ci_u, aggr.yprev[,1], fit.yprev[,1])/0.05)),
             ylab="", xlab="", main=paste(strsex, stragegr))
        cred.region(xx, t(aggr.yprev[,3:4]), col=transp(cols[1], 0.4))
        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))
        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))
          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)])
          graphics::lines(xx, specres.prev, lty=2, lwd=2, col="grey10")
        }
        ##
        if(!is.null(survdat$prev_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 {
          graphics::points(survdat$year, survdat$prev, pch=19)
          graphics::segments(survdat$year, y0=survdat$ci_l, y1=survdat$ci_u)
        }
      }
    }
    graphics::mtext(paste0(icountry, ", ", eppmod, "\nPrevalence trend by age group"), 3, 0, outer=TRUE, font=2, cex=1.3)
  }
  ##
  ##
  return(invisible())
}