R/plotCIN3Incidence.R

In mSimCC: Micro Simulation Model for Cervical Cancer Prevention

plotCIN3Incidence <- function(..., current=NULL, labels=NULL)
{
  if (dev.cur() != 1) dev.off()
  dots <- list(...)
  if (length(dots) < 1) stop("At least one scenario should be defined")
  
  CIN3Incidence1 <- function(i, x, aux){
    num <- aux[i]
    den <- sum(x[i, 1:4])
    num/den
  }
  newCin3Cases1     <- attr(dots[[1]], "newCIN3")
  newCin3CasesMean1 <- apply(newCin3Cases1, 2, mean)
  
  tmp <- lapply(as.list(unique(dots[[1]]$Ages)), function(i, dat){
    aux <- dots[[1]][dots[[1]]$Ages == i, ]
    apply(aux[,1:(ncol(aux)-2)], 2, mean)
  }, dat = dots[[1]])
  meanDat1 <- as.data.frame(t(as.matrix(as.data.frame(tmp))))
  names(newCin3CasesMean1) <- 1:length(newCin3CasesMean1)
  totCIN3Inc1 <- unlist(lapply(as.list(1:nrow(meanDat1)), CIN3Incidence1,
                               x = meanDat1, aux = newCin3CasesMean1))
  age <- seq(min(dots[[1]]$Ages), max(dots[[1]]$Ages), 0.5)
  par(xpd=TRUE)
  colorets <- "#000000"
  plot(age, totCIN3Inc1*attr(dots[[1]], "size"), pch = 20, xlab = "Age", ylab = paste0("CIN 3 incidence x ", attr(dots[[1]], "size")),
       main = "CIN 3 Incidence by age", axes = FALSE, type = "l", col=colorets)
  axis(2)
  axis(1)
  if (length(dots)>1)
  {
    colorets <- c(colorets, rainbow(length(dots)-1))
    for (j in 2:length(dots))
    {
      newCin3Cases1     <- attr(dots[[j]], "newCin3Cases")
      newCin3CasesMean1 <- apply(newCin3Cases1, 2, mean)
      
      tmp <- lapply(as.list(unique(dots[[j]]$Ages)), function(i, dat){
        aux <- dots[[j]][dots[[j]]$Ages == i, ]
        apply(aux[,1:(ncol(aux)-2)], 2, mean)
      }, dat = dots[[j]])
      meanDat1 <- as.data.frame(t(as.matrix(as.data.frame(tmp))))
      names(newCin3CasesMean1) <- 1:length(newCin3CasesMean1)
      totCIN3Inc1 <- unlist(lapply(as.list(1:nrow(meanDat1)), CIN3Incidence1,
                                   x = meanDat1, aux = newCin3CasesMean1))
      age <- min(dots[[j]]$age):max(dots[[j]]$age)
      lines(age, totCIN3Inc1*attr(dots[[1]], "size"), col=colorets[j])
    }
  }
  if (!is.null(current)) 
  {
    colorets  <- c(colorets, "#00ff00")
    h         <- as.numeric(substr(current[1, 1], 4, 5))[1]-as.numeric(substr(current[1, 1], 1, 2))[1]
    age       <- as.numeric(substr(current[1, 1], 1, 2)):as.numeric(substr(current[dim(current)[1], 1], 4, 5))
    CIN3incid <- rep(current[, 2], each=h+1)
    lines(age, CIN3incid, col=colorets[length(colorets)])
  }
  if (is.null(labels)) labels <- ""
  if (!is.null(current)) labels <- c(labels, "Current CIN3 inc.")
  legend("topright", labels, lty=1, col=colorets)
  return(list(data.frame(age=age, val=cbind(as.numeric(apply(newCin3Cases1, 2, mean)))), 
              data.frame(age=age, val=cbind(as.numeric(totCIN3Inc1*attr(dots[[1]], "size"))))))
}