R/calcular.indicadores.R

In lozalojo/mem: The Moving Epidemic Method

#' Calculates specificity and sensitivity
#'
#' @keywords internal
#' @importFrom graphics rect
#' @importFrom grDevices png
#' @importFrom RcppRoll roll_sum
calcular.indicadores <- function(i.current,
                                 i.umbral.pre,
                                 i.umbral.pos = NA,
                                 i.intensidades,
                                 i.duracion.intensidad = 10,
                                 i.equal = F,
                                 i.metodo.calculo = "default",
                                 i.semanas.por.encima = 1,
                                 i.valores.parametro.deteccion = seq(0.1, 5.0, 0.1),
                                 i.output = ".",
                                 i.graph = F,
                                 i.graph.name = "",
                                 i.mem.info = T) {
  if (is.na(i.umbral.pre)) i.umbral.pre <- Inf

  semanas <- dim(i.current)[1]
  nombre.semana <- rownames(i.current)
  numero.semana <- 1:semanas

  # Preparacion de datos necesarios
  umbral.pre <- i.umbral.pre
  if (i.equal) umbral.pos <- i.umbral.pre else umbral.pos <- i.umbral.pos
  if (is.na(umbral.pos)) umbral.pos <- i.umbral.pre
  duracion.media <- i.duracion.intensidad
  punto.maximo <- maxFixNA(i.current)
  semana.punto.maximo <- min(numero.semana[i.current == punto.maximo], na.rm = T)
  semanas.encima <- i.current[, 1] > umbral.pre
  suma.semanas.encima <- c(rep(NA, i.semanas.por.encima - 1), roll_sum(semanas.encima, i.semanas.por.encima))
  if (any(suma.semanas.encima == i.semanas.por.encima, na.rm = T)) {
    semana.inicio.umbral <- min((1:semanas)[suma.semanas.encima == i.semanas.por.encima], na.rm = T)
  } else {
    semana.inicio.umbral <- NA
  }
  if (!is.na(semana.inicio.umbral)) {
    semanas.debajo <- i.current[, 1] < umbral.pos & semana.inicio.umbral < (1:semanas)
    if (any(semanas.debajo, na.rm = T)) {
      semana.fin.umbral <- min((1:semanas)[semanas.debajo], na.rm = T)
    } else {
      semana.fin.umbral <- NA
    }
  } else {
    semana.fin.umbral <- NA
  }
  # M?todo 0: Mediante optimo mem, datos definitivos
  curva.map <- calcular.map(as.vector(as.matrix(i.current)))
  resultado.1 <- numeric()
  n.parametros <- length(i.valores.parametro.deteccion)
  optimos <- numeric()
  for (i in 1:n.parametros) {
    optimo.map <- calcular.optimo(curva.map, 2, i.valores.parametro.deteccion[i])$resultados
    semana.i.optimo <- optimo.map[4]
    if (optimo.map[5] < semanas) semana.f.optimo <- optimo.map[5] + 1 else semana.f.optimo <- NA
    resultado.1.1 <- rep(NA, semanas)
    if (!is.na(semana.i.optimo)) {
      if (semana.i.optimo > 1) resultado.1.1[1:(semana.i.optimo - 1)] <- 1
      if (!is.na(semana.f.optimo)) {
        if (semana.i.optimo < semana.f.optimo) resultado.1.1[semana.i.optimo:(semana.f.optimo - 1)] <- 2
        resultado.1.1[semana.f.optimo:semanas] <- 3
      } else {
        resultado.1.1[semana.i.optimo:semanas] <- 2
      }
    } else {
      resultado.1.1[1:semanas] <- 1
    }
    resultado.1.1[is.na(i.current[, 1])] <- 0
    resultado.1 <- rbind(resultado.1, resultado.1.1)
    optimos <- rbind(optimos, c(semana.i.optimo, semana.f.optimo))
  }

  if (!(i.metodo.calculo == "alternative")) {
    # M?todo 1: Tasas por encima y por debajo del umbral.
    resultado.2 <- rep(NA, semanas)
    if (semana.punto.maximo < semanas) {
      resultado.2.1 <- i.current[(1:semana.punto.maximo), ] > umbral.pre
      resultado.2[1:semana.punto.maximo][resultado.2.1] <- 2
      resultado.2[1:semana.punto.maximo][!resultado.2.1] <- 1
      resultado.2.2 <- i.current[((semana.punto.maximo + 1):semanas), ] >= umbral.pos
      resultado.2[(semana.punto.maximo + 1):semanas][resultado.2.2] <- 2
      resultado.2[(semana.punto.maximo + 1):semanas][!resultado.2.2] <- 3
    } else {
      resultado.2.1 <- i.current[(1:semana.punto.maximo), ] > umbral.pre
      resultado.2[1:semana.punto.maximo][resultado.2.1] <- 2
      resultado.2[1:semana.punto.maximo][!resultado.2.1] <- 1
    }
    resultado.2[is.na(i.current[, 1])] <- 0
  } else {
    # M?todo 2: Hasta iniciar temporada (x tasas por encima), pre, hasta finalizar, epi, despues pos, independientemente de que las tasas est?n por encima o por debajo
    resultado.2 <- rep(NA, semanas)
    if (!is.na(semana.inicio.umbral)) {
      if (semana.inicio.umbral > 1) resultado.2[1:(semana.inicio.umbral - 1)] <- 1
      if (!is.na(semana.fin.umbral)) {
        if (semana.inicio.umbral < semana.fin.umbral) resultado.2[semana.inicio.umbral:(semana.fin.umbral - 1)] <- 2
        resultado.2[semana.fin.umbral:semanas] <- 3
      } else {
        resultado.2[semana.inicio.umbral:semanas] <- 2
      }
    } else {
      resultado.2[1:semanas] <- 1
    }
    resultado.2[is.na(i.current[, 1])] <- 0
  }

  resultado.2 <- matrix(rep(resultado.2, n.parametros), nrow = n.parametros, byrow = T)

  resultado.3 <- numeric()
  for (i in 1:n.parametros) {
    resultado.3.1 <- apply(rbind(resultado.1[i, ], resultado.2[i, ]), 2, comparar.metodos)
    resultado.3 <- rbind(resultado.3, resultado.3.1)
  }

  true.pos.t <- sum(resultado.3 == "TP", na.rm = T)
  true.pos <- apply(resultado.3 == "TP", 1, sum, na.rm = T)
  false.neg.t <- sum(resultado.3 == "FN", na.rm = T)
  false.neg <- apply(resultado.3 == "FN", 1, sum, na.rm = T)
  false.pos.t <- sum(resultado.3 == "FP", na.rm = T)
  false.pos <- apply(resultado.3 == "FP", 1, sum, na.rm = T)
  true.neg.t <- sum(resultado.3 == "TN", na.rm = T)
  true.neg <- apply(resultado.3 == "TN", 1, sum, na.rm = T)

  sensibilidad <- true.pos / (true.pos + false.neg)
  sensibilidad[is.nan(sensibilidad)] <- NA
  especificidad <- true.neg / (true.neg + false.pos)
  especificidad[is.nan(especificidad)] <- NA
  ppv <- true.pos / (true.pos + false.pos)
  ppv[is.nan(ppv)] <- NA
  npv <- true.neg / (true.neg + false.neg)
  npv[is.nan(npv)] <- NA
  pos.likehood.ratio <- sensibilidad / (1 - especificidad)
  pos.likehood.ratio[is.nan(pos.likehood.ratio)] <- NA
  neg.likehood.ratio <- (1 - sensibilidad) / especificidad
  neg.likehood.ratio[is.nan(neg.likehood.ratio)] <- NA
  percent.agreement <- (true.pos + true.neg) / (true.pos + true.neg + false.pos + false.neg)
  percent.agreement[is.nan(percent.agreement)] <- NA
  # Matthews correlation coefficient
  mcc <- (true.pos * true.neg - false.pos * false.neg) / sqrt((true.pos + false.pos) * (true.pos + false.neg) * (true.neg + false.pos) * (true.neg + false.neg))
  mcc[is.nan(mcc)] <- NA
  # Youden's Index
  youden <- sensibilidad + especificidad - 1

  if (true.pos.t + false.neg.t > 0) sensibilidad.t <- true.pos.t / (true.pos.t + false.neg.t) else sensibilidad.t <- NA
  if (true.neg.t + false.pos.t > 0) especificidad.t <- true.neg.t / (true.neg.t + false.pos.t) else especificidad.t <- NA
  if (true.pos.t + false.pos.t > 0) ppv.t <- true.pos.t / (true.pos.t + false.pos.t) else ppv.t <- NA
  if (true.neg.t + false.neg.t > 0) npv.t <- true.neg.t / (true.neg.t + false.neg.t) else npv.t <- NA
  pos.likehood.ratio.t <- NA
  if (!is.na(especificidad.t)) if (1 - especificidad.t > 0) pos.likehood.ratio.t <- sensibilidad.t / (1 - especificidad.t) else pos.likehood.ratio.t <- NA
  neg.likehood.ratio.t <- NA
  if (!is.na(especificidad.t)) if (especificidad.t > 0) neg.likehood.ratio.t <- (1 - sensibilidad.t) / especificidad.t else neg.likehood.ratio.t <- NA
  if (true.pos.t + true.neg.t + false.pos.t + false.neg.t > 0) percent.agreement.t <- (true.pos.t + true.neg.t) / (true.pos.t + true.neg.t + false.pos.t + false.neg.t) else percent.agreement.t <- NA
  if ((true.pos.t + false.pos.t) > 0 & (true.pos.t + false.neg.t) > 0 & (true.neg.t + false.pos.t) > 0 & (true.neg.t + false.neg.t) > 0) mcc.t <- (true.pos.t * true.neg.t - false.pos.t * false.neg.t) / (sqrt(true.pos.t + false.pos.t) * sqrt(true.pos.t + false.neg.t) * sqrt(true.neg.t + false.pos.t) * sqrt(true.neg.t + false.neg.t)) else mcc.t <- NA
  youden.t <- sensibilidad.t + especificidad.t - 1

  semanas.not.na <- sum(!is.na(i.current))

  indicadores.t <- as.matrix(c(
    semanas, semanas.not.na, true.pos.t, false.pos.t, true.neg.t,
    false.neg.t, sensibilidad.t, especificidad.t, ppv.t, npv.t,
    pos.likehood.ratio.t, neg.likehood.ratio.t, percent.agreement.t,
    mcc.t, youden.t
  ))
  rownames(indicadores.t) <- c(
    "Weeks", "Non-missing weeks", "True positives", "False positives",
    "True negatives", "False negatives", "Sensitivity", "Specificity",
    "Positive predictive value", "Negative predictive value",
    "Positive likehood ratio", "Negative likehood ratio",
    "Percent agreement", "Matthews correlation coefficient", "Youdens Index"
  )
  colnames(indicadores.t) <- "values"

  indicadores <- data.frame(
    parametro = i.valores.parametro.deteccion, semanas = semanas, semanas.not.na = semanas.not.na,
    true.pos = true.pos, false.pos = false.pos, true.neg = true.neg, false.neg = false.neg,
    sensibilidad = sensibilidad, especificidad = especificidad, ppv = ppv, npv = npv,
    pos.likehood.ratio = pos.likehood.ratio, neg.likehood.ratio = neg.likehood.ratio,
    percent.agreement = percent.agreement, mcc = mcc, youden = youden
  )

  if (i.graph) {

    ######## Gr?fico

    limites.niveles <- as.vector(i.intensidades)
    # nombres.niveles<-as.character(i.flu$epi.intervalos[,1])
    limites.niveles[limites.niveles < 0] <- 0

    # Datos para el grafico
    if (is.na(semana.inicio.umbral)) {
      # No iniciada
      umbrales.1 <- rep(umbral.pre, semanas)
      umbrales.2 <- rep(NA, semanas)
      intensidades.1 <- array(dim = c(semanas, 3))
      intensidades.2 <- array(dim = c(semanas, 3))
    } else {
      if (is.na(semana.fin.umbral)) {
        # Iniciada y no finalizada
        umbrales.1 <- rep(umbral.pre, semana.inicio.umbral - 1)
        umbrales.2 <- rep(NA, max(duracion.media, semanas - semana.inicio.umbral + 2))
        intensidades.1 <- array(dim = c(max(semana.inicio.umbral - 2, 0), 3))
        intensidades.2 <- matrix(rep(limites.niveles, max(duracion.media, semanas - semana.inicio.umbral + 2)), ncol = 3, byrow = T)
      } else {
        # Iniciada y finalizada
        umbrales.1 <- rep(umbral.pre, semana.inicio.umbral - 1)
        umbrales.2 <- rep(NA, semana.fin.umbral - semana.inicio.umbral)
        intensidades.1 <- array(dim = c(max(semana.inicio.umbral - 2, 0), 3))
        intensidades.2 <- matrix(rep(limites.niveles, semana.fin.umbral - semana.inicio.umbral + 2), ncol = 3, byrow = T)
      }
    }
    umbrales.3 <- rep(umbral.pos, semanas)

    umbrales <- c(umbrales.1, umbrales.2, umbrales.3)[1:semanas]
    intensidades.3 <- array(dim = c(semanas, 3))
    intensidades <- rbind(intensidades.1, intensidades.2, intensidades.3)[1:semanas, ]
    dgraf <- as.data.frame(cbind(i.current, umbrales, intensidades))
    names(dgraf) <- c("Rate", "Epidemic threshold", paste("Intensidad", 1:3))

    etiquetas <- c("Weekly rates", "Epidemic thr", "Medium thr", "High thr", "Very high thr")
    tipos <- c(1, 2, 2, 2, 2)
    anchos <- c(3, 2, 2, 2, 2)
    colores <- c("#808080", "#8c6bb1", "#88419d", "#810f7c", "#4d004b")
    colores.epi <- c("#00C000", "#980043", "#FFB401")

    range.x <- 1:semanas

    # calculo el rango y para que tenga 10 marcas o este cerca

    maximo.y <- maxFixNA(dgraf)
    posicion.ticks <- optimal.tickmarks(0, maximo.y, 10)$by
    range.y <- c(-1.5 * posicion.ticks, ceiling(maximo.y / posicion.ticks) * posicion.ticks)
    range.y.seq <- seq(0, ceiling(maximo.y / posicion.ticks) * posicion.ticks, posicion.ticks)

    for (i in 1:n.parametros) {
      if (i.graph.name == "") graph.name <- paste("surveillance graph (", format(round(i.valores.parametro.deteccion[i], 1), digits = 3, nsmall = 1), ")", sep = "") else graph.name <- paste(i.graph.name, " (", format(round(i.valores.parametro.deteccion[i], 1), digits = 3, nsmall = 1), ")", sep = "")

      png(
        filename = file.path(i.output, paste0(graph.name, ".png")),
        width = 8, height = 6, units = "in", pointsize = "12",
        bg = "white", res = 300, antialias = "none"
      )
      opar <- par(mar = c(4, 4, 1, 8) + 0.1, xpd = TRUE)
      # ,mgp=c(3,0.5,0),xpd=T)
      # Grafico principal
      matplot(range.x,
        dgraf,
        type = "l",
        lty = tipos, lwd = anchos, col = colores,
        xlab = "", ylab = "", axes = F,
        ylim = range.y
      )
      # Puntos de la serie de tasas
      # points(1:semanas,dgraf[,1],pch=19,type="p",col="#000000",cex=1)
      # pre
      puntos.1 <- i.current[, 1]
      puntos.1[!(resultado.1[i, ] == 1)] <- NA
      # if (!is.na(optimos[i,1])) puntos.1[optimos[i,1]:semanas]<-NA
      points(puntos.1, pch = 19, type = "p", col = colores.epi[1], cex = 1.5)
      # epi
      puntos.2 <- i.current[, 1]
      puntos.2[!(resultado.1[i, ] == 2)] <- NA
      # if (is.na(optimos[i,1])){
      #   puntos.2[1:semanas]<-NA
      # }else{
      #   if (optimos[i,1]>1) puntos.2[1:(optimos[i,1]-1)]<-NA
      #   if (!is.na(optimos[i,2])) puntos.2[optimos[i,2]:semanas]<-NA
      # }
      points(puntos.2, pch = 19, type = "p", col = colores.epi[2], cex = 1.5)
      # post
      puntos.3 <- i.current[, 1]
      puntos.3[!(resultado.1[i, ] == 3)] <- NA
      # if (is.na(optimos[i,1])){
      #   puntos.3[1:semanas]<-NA
      # }else{
      #   if (is.na(optimos[i,2])) puntos.3[1:semanas]<-NA else puntos.3[1:(optimos[i,2]-1)]<-NA
      # }
      points(puntos.3, pch = 19, type = "p", col = colores.epi[3], cex = 1.5)
      # Ejes
      axis(2, at = range.y.seq, lwd = 1, cex.axis = 0.6, col.axis = "#404040", col = "#C0C0C0", mgp = c(3, 0.5, 0))
      mtext(2, text = "Weekly rate", line = 2, cex = 0.8, col = "#000040")
      axis(1, pos = 0, at = seq(1, semanas, 1), labels = F, cex.axis = 0.7, col.axis = "#404040", col = "#C0C0C0")
      axis(1, at = seq(0.5, semanas + 0.5, 1), labels = F, cex.axis = 0.7, col.axis = "#404040", col = "#C0C0C0")
      axis(1,
        at = seq(1, semanas, 2), tick = F, mgp = c(3, 0.5, 0),
        labels = nombre.semana[seq(1, semanas, 2)], cex.axis = 0.6, col.axis = "#404040", col = "#C0C0C0"
      )
      axis(1,
        at = seq(2, semanas, 2), tick = F, mgp = c(3, 0.5, 0),
        labels = nombre.semana[seq(2, semanas, 2)], cex.axis = 0.6, line = 0.60, col.axis = "#404040", col = "#C0C0C0"
      )
      mtext(1, text = "Week", line = 2.5, cex = 0.8, col = "#000040")
      if (i.mem.info) {
        mtext(4,
          text = paste("mem R library - Jose E. Lozano - https://github.com/lozalojo/mem", sep = ""),
          line = 7, cex = 0.6, col = "#404040"
        )
      }
      # Etiquetas de los 4 umbrales

      if (is.na(semana.inicio.umbral)) {
        # No iniciada
        text.x <- semanas - semanas / 25
        text.y <- umbral.pre
        text.l <- round(umbral.pre, 2)
        text.p <- 3
        text.s <- 1
        text.c <- colores[2]
      } else {
        if (is.na(semana.fin.umbral)) {
          # Iniciada y no finalizada
          lugar.intensidad <- min(semanas, semana.inicio.umbral + max(duracion.media, semanas - semana.inicio.umbral + 1))
        } else {
          # Iniciada y finalizada
          lugar.intensidad <- semana.fin.umbral
        }
        text.x <- c(semana.inicio.umbral, rep(lugar.intensidad, 3), semanas) - semanas / 25
        text.y <- c(umbral.pre, limites.niveles, umbral.pos)
        text.l <- round(c(umbral.pre, limites.niveles, umbral.pos), 2)
        text.p <- rep(3, 5)
        text.s <- rep(1, 5)
        text.c <- colores[c(2:5, 2)]
      }
      text(text.x, text.y, text.l, pos = text.p, col = text.c, cex = text.s)

      # Etiquetas de la leyenda
      etiquetas.leyenda <- c("Pre", "Epidemic", "Post", etiquetas)
      tipos.leyenda <- c(1, 1, 1, tipos)
      anchos.leyenda <- c(1, 1, 1, anchos)
      colores.leyenda <- c("#C0C0C0", "#C0C0C0", "#C0C0C0", colores)
      puntos.leyenda <- c(21, 21, 21, rep(NA, 5))
      bg.leyenda <- c("#FFFFFF", "#FFFFFF", "#FFFFFF", rep(NA, 5))
      pt.bg.leyenda <- c("#00C000", "#980043", "#FFB401", rep(NA, 5))

      legend("topright",
        inset = c(-0.25, 0),
        legend = rev(etiquetas.leyenda),
        bty = "n",
        lty = rev(tipos.leyenda),
        lwd = rev(anchos.leyenda),
        col = rev(colores.leyenda),
        pch = rev(puntos.leyenda),
        bg = rev(bg.leyenda),
        pt.bg = rev(pt.bg.leyenda),
        cex = 0.9,
        text.col = "#000000",
        ncol = 1
      )

      colores.cuadros <- c("#C0C0C0", colores.epi)
      rect(seq(0.5, semanas - 0.5, 1),
        -3 * posicion.ticks / 2,
        seq(1.5, (semanas + 0.5), 1),
        -2 * posicion.ticks / 2,
        density = NULL,
        col = colores.cuadros[1 + resultado.2[i, ]],
        border = "white"
      )
      rect(seq(0.5, semanas - 0.5, 1),
        -2 * posicion.ticks / 2,
        seq(1.5, (semanas + 0.5), 1),
        -1 * posicion.ticks / 2,
        density = NULL,
        col = colores.cuadros[1 + resultado.1[i, ]],
        border = "white"
      )

      text.xy <- expand.grid(x = seq(1, semanas, 1), y = seq(-2.5, -1.5, 1) * posicion.ticks / 2)
      text.l <- c(resultado.3, rep(NA, semanas))
      text.p <- rep(1, semanas * 2)
      text.s <- rep(0.25, semanas * 2)
      text.c <- rep("#FFFFFF", semanas * 2)
      text(text.xy$x, text.xy$y, text.l, col = text.c, cex = text.s)

      if (!(i.metodo.calculo == "alternative")) {
        text.method <- paste0("* Method used: ", i.metodo.calculo, " (weeks above/below the threshold)")
      } else {
        text.method <- paste0("* Method used: ", i.metodo.calculo, " (", i.semanas.por.encima, " week above the threshold)")
      }

      text(semanas, -3.5 * posicion.ticks / 2, pos = 2, label = paste("Sensitivity: ", format(round(sensibilidad[i], 2), nsmall = 2, align = "right"), ", Specificity: ", format(round(especificidad[i], 2), nsmall = 2, align = "right"), sep = ""), cex = 0.5)
      text(1, -3.5 * posicion.ticks / 2, pos = 4, label = text.method, cex = 0.5)

      axis(2,
        at = seq(-2.5, -1.5, 1) * posicion.ticks / 2,
        labels = c("threshold *", paste("algorithm (", format(round(i.valores.parametro.deteccion[i], 1), digits = 3, nsmall = 1), ")", sep = "")),
        tick = F,
        las = 1,
        lwd = 1,
        cex.axis = 0.6, col.axis = "#404040", col = "#C0C0C0", mgp = c(3, 0, 0)
      )

      # Etiquetas de la leyenda
      etiquetas.leyenda <- c("Pre", "Epidemic", "Post")
      tipos.leyenda <- c(0, 0, 0)
      anchos.leyenda <- c(0, 0, 0)

      colores.leyenda <- c("#C0C0C0", "#C0C0C0", "#C0C0C0")
      puntos.leyenda <- c(21, 21, 21)
      bg.leyenda <- c("#FFFFFF", "#FFFFFF", "#FFFFFF")
      pt.bg.leyenda <- c("#00C000", "#980043", "#FFB401")
      legend("bottomright",
        inset = c(-0.18, 0.025),
        x.intersp = -1,
        # legend=rev(etiquetas.leyenda),
        legend = etiquetas.leyenda,
        bty = "n",
        lty = rev(tipos.leyenda),
        lwd = rev(anchos.leyenda),
        fill = pt.bg.leyenda,
        cex = 0.9,
        text.col = "#000000",
        ncol = 1
      )

      par(opar)
      dev.off()
    }
  }

  return(list(resultado.1 = resultado.1, resultado.2 = resultado.2, resultado.3 = resultado.3, indicadores = indicadores, indicadores.t = indicadores.t))
}