R/memmodel.R
In lozalojo/mem: The Moving Epidemic Method
Defines functions
plot.mem
summary.mem
print.mem
memmodel
Documented in
memmodel
#' @title Methods for influenza modelization
#'
#' @description
#' Function \code{memmodel} is used to calculate the threshold for influenza epidemic using historical
#' records (surveillance rates).
#'
#' The method to calculate the threshold is described in the Moving Epidemics Method (MEM) used to
#' monitor influenza activity in a weekly surveillance system.
#'
#' @name memmodel
#'
#' @aliases summary.flu,plot.flu,print.flu
#'
#' @param i.data Data frame of input data.
#' @param i.seasons Maximum number of seasons to use.
#' @param i.type.threshold Type of confidence interval to calculate the threshold.
#' @param i.level.threshold Level of confidence interval to calculate the threshold.
#' @param i.tails.threshold Tails for the confidence interval to calculate the threshold.
#' @param i.type.intensity Type of confidence interval to calculate the intensity thresholds.
#' @param i.level.intensity Levels of confidence interval to calculate the intensity thresholds.
#' @param i.tails.intensity Tails for the confidence interval to calculate the threshold.
#' @param i.type.curve Type of confidence interval to calculate the modelled curve.
#' @param i.level.curve Level of confidence interval to calculate the modelled curve.
#' @param i.type.other Type of confidence interval to calculate length, start and percentages.
#' @param i.level.other Level of confidence interval to calculate length, start and percentages.
#' @param i.method Method to calculate the optimal timing of the epidemic.
#' @param i.param Parameter to calculate the optimal timing of the epidemic.
#' @param i.centering Number of weeks to center the moving seasons.
#' @param i.n.max Number of pre-epidemic values used to calculate the threshold.
#' @param i.type.boot Type of bootstrap technique.
#' @param i.iter.boot Number of bootstrap iterations.
#' @param i.mem.info include information about the package in the graph.
#'
#' @return
#' \code{memmodel} returns an object of class \code{mem}.
#' An object of class \code{mem} is a list containing at least the following components:
#' \itemize{
#' \item{i.data} {input data}
#' \item{pre.post.intervals} {Pre/post confidence intervals (Threhold is the upper limit
#' of the confidence interval).}
#' \item{ci.length} {Mean epidemic length confidence interval.}
#' \item{ci.percent} {Mean covered percentage confidence interval.}
#' \item{mean.length} {Mean length.}
#' \item{moving.epidemics} {Moving epidemic rates.}
#' \item{mean.start} {Mean epidemic start.}
#' \item{epi.intervals} {Epidemic levels of intensity.}
#' \item{typ.curve} {Typical epidemic curve.}
#' \item{n.max} {Effective number of pre epidemic values.}
#' }
#'
#' @details
#' Input data is a data frame containing rates that represent historical influenza surveillance
#' data. It can start and end at any given week (tipically at week 40th), and rates can be
#' expressed as per 100,000 inhabitants (or per consultations, if population is not
#' available) or any other scale.
#'
#' Parameters \code{i.type}, \code{i.type.threshold} and \code{i.type.curve} defines how to
#' calculate confidence intervals along the process.
#'
#' \code{i.type.curve} is used for calculating the typical influenza curve,
#' \code{i.type.threshold} is used to calculate the pre and post epidemic threshold and
#' \code{i.type} is used for any other confidende interval used in the method.
#'
#' All three parameters must be a number between \code{1} and \code{6}:
#'
#' \itemize{
#' \item{1} {Arithmetic mean and mean confidence interval.}
#' \item{2} {Geometric mean and mean confidence interval.}
#' \item{3} {Median and the Nyblom/normal aproximation confidence interval.}
#' \item{4} {Median and bootstrap confidence interval.}
#' \item{5} {Arithmetic mean and point confidence interval (standard deviations).}
#' \item{6} {Geometric mean and point confidence interval (standard deviations).}
#' }
#'
#' Option \code{3} uses the Hettmansperger and Sheather (1986) and Nyblom (1992) method,
#' when there is enough sample size. If sample size is small, then the normal aproximation
#' will be used as described in Conover, 1980, p. 112. Refer to EnvStats package for
#' more information.
#'
#' Option \code{4} uses two more parameters: \code{i.type.boot} indicates which bootstrap
#' method to use. The values are the same of those of the \code{\link{boot.ci}} function.
#' Parameter \code{i.iter.boot} indicates the number of bootstrap samples to use. See
#' \code{\link{boot}} for more information about this topic.
#'
#' Parameters \code{i.level}, \code{i.level.threshold} and \code{i.level.curve} indicates,
#' respectively, the level of the confidence intervals described above.
#'
#' The \code{i.n.max} parameter indicates how many pre epidemic values to use to calculate
#' the threshold. A value of -1 indicates the program to use an appropiate number of points
#' depending on the number of seasons provided as input. \code{i.tails} tells the program
#' to use {1} or {2} tailed confidence intervals when calculating the threshold (1 is
#' recommended).
#'
#' Parameters \code{i.method} and \code{i.param} indicates how to find the optimal timing
#' of the epidemics. See \code{\link{memtiming}} for details on the values this parameters
#' can have.
#'
#' It is important to know how to arrange information in order to use with memapp. The key points are:
#'
#' \itemize{
#' \item One single epidemic wave each season.
#' \item Never delete a rate inside an epidemic.
#' \item Accommodate week 53.
#' \item Do not inflate missing values with zeroes.
#' }
#'
#' Data must contain information from the historical series. Surveillance period can start and end at
#' any given week (typically start at week 40th and ends at week 20th), and data can have any units and
#' can be expressed in any scale (typically rates per 100,000 inhabitants or consultations).
#'
#' The table must have one row per epidemiological week and one column per surveillance season. A season
#' is a full surveillance period from the beginning to the end, where occurs at some point one single
#' epidemic wave on it. No epidemic wave can be spared in two consecutive seasons. If so, you have to
#' redefine the start and end of the season defined in your dataset. If a season have two waves, it must
#' be split in two periods and must be named accordingly with the seasons name conventions described
#' below. Each cell contains the value for a given week in a given season.
#'
#' The first column should contain the names of the weeks. When the season contains two different calendar
#' years, the week will go from 40th of the first year to 52nd, and then from 1st to 20th. When the season
#' contains one year, the weeks will go from 1st to 52nd.
#'
#' Note: If there is no column with week names, the application will name the weeks numbering from 1 to
#' the number of rows.
#'
#' \itemize{
#' \item In the northern hemisphere countries, the surveillance period usually goes from
#' week 40 to 20 of the following year (notation: season 2016/2017).
#' \item In the southern hemisphere countries, the surveillance period usually goes from
#' week 18 to 39 same year (notation: season 2017).
#' }
#'
#' The first row must contain the names of the seasons. This application understand the naming of a
#' season when it contains one or two four digits year separated by / and one one-digit number
#' between parenthesis to identify the wave number. The wave number part in a name of a season is
#' used when a single surveillance period has two epidemic waves that have to be separated in order
#' to have reliable results. In this case, each wave is placed in different columns and named ending
#' with (1) for the first period, (2) for the second, and so on.
#'
#' @examples
#' # Castilla y Leon Influenza Rates data
#' data(flucyl)
#' # Finds the timing of the first season: 2001/2002
#' epi <- memmodel(flucyl)
#' print(epi)
#' summary(epi)
#' plot(epi)
#' @author Jose E. Lozano \email{lozalojo@@gmail.com}
#'
#' @references
#' Vega T, Lozano JE, Ortiz de Lejarazu R, Gutierrez Perez M. Modelling influenza epidemic - can we
#' detect the beginning and predict the intensity and duration? Int Congr Ser. 2004 Jun;1263:281-3.
#'
#' Vega T, Lozano JE, Meerhoff T, Snacken R, Mott J, Ortiz de Lejarazu R, et al. Influenza surveillance
#' in Europe: establishing epidemic thresholds by the moving epidemic method. Influenza Other Respir
#' Viruses. 2013 Jul;7(4):546-58. DOI:10.1111/j.1750-2659.2012.00422.x.
#'
#' Vega T, Lozano JE, Meerhoff T, Snacken R, Beaute J, Jorgensen P, et al. Influenza surveillance in
#' Europe: comparing intensity levels calculated using the moving epidemic method. Influenza Other
#' Respir Viruses. 2015 Sep;9(5):234-46. DOI:10.1111/irv.12330.
#'
#' Lozano JE. lozalojo/mem: Second release of the MEM R library. Zenodo [Internet]. [cited 2017 Feb 1];
#' Available from: \url{https://zenodo.org/record/165983}. DOI:10.5281/zenodo.165983
#'
#' Hettmansperger, T. P., and S. J Sheather. 1986. Confidence Intervals Based on Interpolated Order
#' Statistics. Statistics and Probability Letters 4: 75-79. doi:10.1016/0167-7152(86)90021-0.
#'
#' Nyblom, J. 1992. Note on Interpolated Order Statistics. Statistics and Probability Letters 14:
#' 129-31. doi:10.1016/0167-7152(92)90076-H.
#'
#' Conover, W.J. (1980). Practical Nonparametric Statistics. Second Edition. John Wiley and Sons, New York.
#'
#' @keywords influenza
#'
#' @export
#' @importFrom stats loess qt qnorm quantile runif var
#' @importFrom grDevices colorRampPalette colorRamp dev.off rgb tiff
#' @importFrom RColorBrewer brewer.pal
#' @importFrom graphics abline axis legend matplot mtext par points text lines
memmodel <- function(i.data,
i.seasons = 10,
i.type.threshold = 5,
i.level.threshold = 0.95,
i.tails.threshold = 1,
i.type.intensity = 6,
i.level.intensity = c(0.40, 0.90, 0.975),
i.tails.intensity = 1,
i.type.curve = 2,
i.level.curve = 0.95,
i.type.other = 3,
i.level.other = 0.95,
i.method = 2,
i.param = 2.8,
i.centering = -1,
i.n.max = -1,
i.type.boot = "norm",
i.iter.boot = 10000,
i.mem.info = T) {
if (is.null(dim(i.data))) {
memmodel.output <- NULL
cat("Incorrect number of dimensions, input must be a data.frame\n")
} else {
# Test if first column is week.name
if (all(i.data[,1] %in% 1:53)){
datos <- i.data[-1]
rownames(datos)<-i.data[,1]
}else{
datos <- i.data
if (!all(rownames(i.data) %in% 1:53)) rownames(datos)<-1:NROW(i.data)
}
datos <- datos[apply(datos, 2, function(x) sum(x, na.rm = T) > 0)]
if (!(ncol(datos) > 1)) {
memmodel.output <- NULL
cat("Incorrect number of dimensions, at least two seasons of valid data (non-zero values) required\n")
} else {
if (is.matrix(datos)) datos <- as.data.frame(datos)
if (is.null(i.seasons)) i.seasons <- NA
if (!is.na(i.seasons)) if (i.seasons > 0) datos <- datos[(max((dim(datos)[2]) - i.seasons + 1, 1)):(dim(datos)[2])]
datos <- as.data.frame(apply(datos, 2, fill.missing))
semanas <- dim(datos)[1]
anios <- NCOL(datos)
# Calcular el optimo
if (is.na(i.n.max)) {
n.max <- max(1, round(30 / anios, 0))
} else {
if (i.n.max == -1) {
n.max <- max(1, round(30 / anios, 0))
# if (anios>=10) n.max=3 else n.max=5
} else if (i.n.max == 0) {
n.max <- semanas
} else {
n.max <- i.n.max
}
}
if (is.null(i.method)) i.method <- 2
if (is.na(i.method)) i.method <- 2
if (is.null(i.param)) i.param <- 2.8
if (is.na(i.param)) i.param <- 2.8
optimo <- apply(datos, 2, memtiming, i.n.values = n.max, i.method = i.method, i.param = i.param)
datos.duracion.real <- extraer.datos.optimo.map(optimo)
# por defecto la aritmetica en ambos
ic.duracion <- iconfianza(as.numeric(datos.duracion.real[1, ]), i.level.other, i.type.other, T, i.type.boot, i.iter.boot, 2)
ic.duracion <- rbind(ic.duracion, c(floor(ic.duracion[1]), round(ic.duracion[2]), ceiling(ic.duracion[3])))
ic.porcentaje <- iconfianza(as.numeric(datos.duracion.real[2, ]), i.level.other, i.type.other, T, i.type.boot, i.iter.boot, 2)
duracion.media <- ic.duracion[2, 2]
if (is.na(i.centering) | i.centering == -1) duracion.centrado <- duracion.media else duracion.centrado <- i.centering
########################################################################################
# Calcula todos los parametros asociados a la temporada gripal en base a la estimacion #
# del numero de semanas de duracion total. #
########################################################################################
# Datos de la gripe en la temporada REAL y en la temporada MEDIA, esta ultima sirve para unir las temporadas.
semana.inicio <- as.integer(rownames(datos)[1])
gripe <- array(dim = c(7, anios, 3), dimnames = list(
"indicators" = c("epidemic start", "epidemic end", "sum of epidemic values", "sum of values", "covered percentage", "sum of pre-epidemic values", "sum of post-epidemic values"),
"seasons" = names(datos),
"type" = c("real epidemic", "mean length epidemic", "centered length epidemic")
))
datos.duracion.media <- extraer.datos.curva.map(optimo, duracion.media)
datos.duracion.centrado <- extraer.datos.curva.map(optimo, duracion.centrado)
## Estimacion del numero de semana en el comienza la gripe
ic.inicio <- iconfianza(as.numeric(datos.duracion.real[4, ]), i.level.other, i.type.other, T, i.type.boot, i.iter.boot, 2)
ic.inicio <- rbind(ic.inicio, c(floor(ic.inicio[1]), round(ic.inicio[2]), ceiling(ic.inicio[3])))
ic.inicio <- rbind(ic.inicio, c(semana.absoluta(ic.inicio[2, 1], semana.inicio), semana.absoluta(ic.inicio[2, 2], semana.inicio), semana.absoluta(ic.inicio[2, 3], semana.inicio)))
ic.inicio <- ic.inicio[c(2, 3, 1), ]
inicio.medio <- ic.inicio[1, 2]
ic.inicio.centrado <- iconfianza(as.numeric(datos.duracion.centrado[4, ]), i.level.other, i.type.other, T, i.type.boot, i.iter.boot, 2)
ic.inicio.centrado <- rbind(ic.inicio.centrado, c(floor(ic.inicio.centrado[1]), round(ic.inicio.centrado[2]), ceiling(ic.inicio.centrado[3])))
ic.inicio.centrado <- rbind(ic.inicio.centrado, c(semana.absoluta(ic.inicio.centrado[2, 1], semana.inicio), semana.absoluta(ic.inicio.centrado[2, 2], semana.inicio), semana.absoluta(ic.inicio.centrado[2, 3], semana.inicio)))
ic.inicio.centrado <- ic.inicio.centrado[c(2, 3, 1), ]
inicio.centrado <- ic.inicio.centrado[1, 2]
for (j in 1:anios) {
## Semana en la q comienza la temporada.real
gripe[1, j, 1] <- datos.duracion.real[4, j]
## Semana en la q acaba
gripe[2, j, 1] <- datos.duracion.real[5, j]
## numero de casos en la temporada.real de gripe
gripe[3, j, 1] <- datos.duracion.real[3, j]
## Numero de casos totales de la temporada global
gripe[4, j, 1] <- sum(datos[, j], na.rm = TRUE)
## Porcentaje cubierto de casos de la temporada gripal
gripe[5, j, 1] <- datos.duracion.real[2, j]
## Casos Acumulados justo antes de comenzar la temporada gripal
if (gripe[1, j, 1] > 1) gripe[6, j, 1] <- sum(datos[1:(gripe[1, j, 1] - 1), j], na.rm = TRUE) else gripe[6, j, 1] <- 0
## Casos Despues de la temporada gripal
if (gripe[2, j, 1] < semanas) gripe[7, j, 1] <- sum(datos[(gripe[2, j, 1] + 1):semanas, j], na.rm = TRUE) else gripe[7, j, 1] <- 0
gripe[1, j, 2] <- datos.duracion.media[4, j]
gripe[2, j, 2] <- datos.duracion.media[5, j]
gripe[3, j, 2] <- datos.duracion.media[3, j]
gripe[4, j, 2] <- sum(datos[, j], na.rm = TRUE)
gripe[5, j, 2] <- datos.duracion.media[2, j]
if (gripe[1, j, 2] > 1) gripe[6, j, 2] <- sum(datos[1:(gripe[1, j, 2] - 1), j], na.rm = TRUE) else gripe[6, j, 2] <- 0
if (gripe[2, j, 2] < semanas) gripe[7, j, 2] <- sum(datos[(gripe[2, j, 2] + 1):semanas, j], na.rm = TRUE) else gripe[7, j, 2] <- 0
gripe[1, j, 3] <- datos.duracion.centrado[4, j]
gripe[2, j, 3] <- datos.duracion.centrado[5, j]
gripe[3, j, 3] <- datos.duracion.centrado[3, j]
gripe[4, j, 3] <- sum(datos[, j], na.rm = TRUE)
gripe[5, j, 3] <- datos.duracion.centrado[2, j]
if (gripe[1, j, 3] > 1) gripe[6, j, 3] <- sum(datos[1:(gripe[1, j, 3] - 1), j], na.rm = TRUE) else gripe[6, j, 3] <- 0
if (gripe[2, j, 3] < semanas) gripe[7, j, 3] <- sum(datos[(gripe[2, j, 3] + 1):semanas, j], na.rm = TRUE) else gripe[7, j, 3] <- 0
}
# Indice de estado temporal. Si es un 1 me dice que estamos en epoca pretemporada,
# si pone un 2 indica que estamos en temporada gripal. Si pone un 3 indica que
# estamos en posttemporada.
indices.temporada <- array(dim = c(semanas, anios, 3), dimnames = list(
"weeks" = rownames(datos),
"seasons" = names(datos),
"type" = c("real epidemic", "mean length epidemic", "centered length epidemic")
))
for (j in 1:anios) {
indices.temporada[-((datos.duracion.real[4, j]):semanas), j, 1] <- 1
indices.temporada[(datos.duracion.real[4, j]):(datos.duracion.real[5, j]), j, 1] <- 2
indices.temporada[-(1:(datos.duracion.real[5, j])), j, 1] <- 3
indices.temporada[-((datos.duracion.media[4, j]):semanas), j, 2] <- 1
indices.temporada[(datos.duracion.media[4, j]):(datos.duracion.media[5, j]), j, 2] <- 2
indices.temporada[-(1:(datos.duracion.media[5, j])), j, 2] <- 3
indices.temporada[-((datos.duracion.centrado[4, j]):semanas), j, 3] <- 1
indices.temporada[(datos.duracion.centrado[4, j]):(datos.duracion.centrado[5, j]), j, 3] <- 2
indices.temporada[-(1:(datos.duracion.centrado[5, j])), j, 3] <- 3
}
# Grafico del esquema de las temporadas gripales para su union
longitud.esquema <- semanas + maxFixNA(datos.duracion.centrado[4, ]) - minFixNA(datos.duracion.centrado[4, ])
inicio.epidemia.esquema <- maxFixNA(datos.duracion.centrado[4, ])
fin.epidemia.esquema <- maxFixNA(datos.duracion.centrado[5, ])
esquema.temporadas <- array(dim = c(longitud.esquema, anios + 2, 3), dimnames = list(
"dummy weeks" = as.character(1:longitud.esquema),
"seasons" = c(names(datos), "mean season week", "mean season index"),
"type" = c("absolute week", "relative week", "values")
))
for (j in 1:anios) {
diferencia.anual <- inicio.epidemia.esquema - datos.duracion.centrado[4, j]
for (i in 1:semanas) {
esquema.temporadas[i + diferencia.anual, j, 1] <- i
esquema.temporadas[i + diferencia.anual, j, 2] <- semana.absoluta(i, semana.inicio)
esquema.temporadas[i + diferencia.anual, j, 3] <- datos[i, j]
}
}
diferencia.global <- inicio.epidemia.esquema - inicio.centrado
for (i in 1:semanas) {
if ((i + diferencia.global) >= 1 & (i + diferencia.global) <= longitud.esquema) {
esquema.temporadas[i + diferencia.global, anios + 1, 1] <- i
esquema.temporadas[i + diferencia.global, anios + 1, 2] <- semana.absoluta(i, semana.inicio)
esquema.temporadas[i + diferencia.global, anios + 1, 3] <- i
}
}
esquema.temporadas[, anios + 2, ] <- 3
esquema.temporadas[1:(diferencia.global + inicio.medio - 1), anios + 2, ] <- 1
esquema.temporadas[(diferencia.global + inicio.medio):(diferencia.global + inicio.medio + duracion.media - 1), anios + 2, ] <- 2
# esquema.temporadas[1:(inicio.epidemia.esquema-1),anios+2,2]<-1
# esquema.temporadas[(inicio.epidemia.esquema):(inicio.epidemia.esquema+duracion.media-1),anios+2,2]<-2
## Temporadas moviles
temporadas.moviles <- esquema.temporadas[, c(1:anios), 3]
## Limites de la temporada
# limites.temporada<-c(inicio.medio,inicio.medio+duracion.media-1)
# limites.temporada<-rbind(limites.temporada,c(semana.absoluta(inicio.medio,semana.inicio),semana.absoluta(inicio.medio+duracion.media-1,semana.inicio)))
# Limites relativos
# limites.esquema<-c(inicio.epidemia.esquema,inicio.epidemia.esquema+duracion.media-1)
## En ttotal reordenamos la gripe para q todas las temporadas comiencen la semana "estinicio" (es decir, quitamos informacion
## por delante y por detras de esquema.temporada[,,3] para que queden 35 semanas (las que hay).
temporadas.moviles.recortada <- array(dim = c(semanas, anios), dimnames = list(
"dummy weeks" = as.character(1:semanas),
"seasons" = names(datos)
))
## Temporada gripal comienza "estinicio" y termina en "estinicio+temporada-1"
for (j in 1:anios) {
for (i in 1:semanas) {
if ((i + diferencia.global) >= 1 & (i + diferencia.global) <= longitud.esquema) {
temporadas.moviles.recortada[i, j] <- temporadas.moviles[i + diferencia.global, j]
} else {
temporadas.moviles.recortada[i, j] <- NA
}
}
}
## Dibujamos la curva epidemica tipo
# Importante: ?Quitamos los ceros?
# iy<-(is.na(temporadas.moviles.recortada) | temporadas.moviles.recortada==0)
iy <- is.na(temporadas.moviles.recortada)
temporadas.moviles.recortada.no.ceros <- temporadas.moviles.recortada
temporadas.moviles.recortada.no.ceros[iy] <- NA
curva.tipo <- t(apply(temporadas.moviles.recortada.no.ceros, 1, iconfianza, nivel = i.level.curve, tipo = i.type.curve, ic = T, tipo.boot = i.type.boot, iteraciones.boot = i.iter.boot, colas = 2))
## Seleccionamos los periodos pre, epidemia, y post
## PRE y POST-TEMPORADA GRIPAL
## Como no se registra mas q de la semana 40 a la 20, tenemos q muchas de las semanas tienen tasa 0, eliminamos esas semanas,
## ya q podrian llevar a una infraestimacion de la tasa base fuera de temporada
pre.datos <- as.vector(as.matrix(extraer.datos.pre.epi(optimo)))
post.datos <- as.vector(as.matrix(extraer.datos.post.epi(optimo)))
epi.datos <- as.vector(as.matrix(extraer.datos.epi(optimo)))
epi.datos.2 <- as.vector(as.matrix(apply(datos, 2, maxnvalores, n.max = n.max)))
pre.post.datos <- rbind(pre.datos, post.datos)
# IC de la linea basica de pre y post temporada
# por defecto estaba la geometrica
# pre.d<-pre.datos[!(is.na(pre.datos) | pre.datos==0)]
# post.d<-post.datos[!(is.na(post.datos) | post.datos==0)]
# epi.d<-epi.datos[!(is.na(epi.datos) | epi.datos==0)]
# epi.d.2<-epi.datos.2[!(is.na(epi.datos.2) | epi.datos.2==0)]
pre.d <- pre.datos[!is.na(pre.datos)]
post.d <- post.datos[!is.na(post.datos)]
epi.d <- epi.datos[!is.na(epi.datos)]
epi.d.2 <- epi.datos.2[!is.na(epi.datos.2)]
pre.i <- iconfianza(pre.d, nivel = i.level.threshold, tipo = i.type.threshold, ic = T, tipo.boot = i.type.boot, iteraciones.boot = i.iter.boot, colas = i.tails.threshold)
post.i <- iconfianza(post.d, nivel = i.level.threshold, tipo = i.type.threshold, ic = T, tipo.boot = i.type.boot, iteraciones.boot = i.iter.boot, colas = i.tails.threshold)
epi.intervalos <- numeric()
for (niv in i.level.intensity) epi.intervalos <- rbind(epi.intervalos, c(niv, iconfianza(epi.d, nivel = niv, tipo = i.type.intensity, ic = T, colas = i.tails.intensity)))
epi.intervalos.2 <- numeric()
for (niv in i.level.intensity) epi.intervalos.2 <- rbind(epi.intervalos.2, c(niv, iconfianza(epi.d.2, nivel = niv, tipo = i.type.intensity, ic = T, colas = i.tails.intensity)))
pre.post.intervalos <- rbind(pre.i, post.i)
memmodel.output <- list(
pre.post.intervals = pre.post.intervalos,
ci.length = ic.duracion,
ci.percent = ic.porcentaje,
ci.start = ic.inicio,
mean.length = duracion.media,
mean.start = inicio.medio,
centered.length = duracion.centrado,
centered.start = inicio.centrado,
moving.epidemics = temporadas.moviles.recortada,
epi.intervals = epi.intervalos,
epi.intervals.full = epi.intervalos.2,
typ.curve = curva.tipo,
n.max = n.max,
n.seasons = anios,
n.weeks = semanas,
data = datos,
start.week = semana.inicio,
epi.data = epi.datos,
epi.data.nona = epi.d,
epi.data.full = epi.datos.2,
epi.data.nona.full = epi.d.2,
optimum = optimo,
real.length.data = datos.duracion.real,
seasons.data = gripe,
season.indexes = indices.temporada,
season.scheme = esquema.temporadas,
seasons.moving = temporadas.moviles,
pre.post.data = pre.post.datos,
pre.data.nona = pre.d,
post.data.nona = post.d,
epidemic.thresholds = as.numeric(pre.post.intervalos[, 3]),
intensity.thresholds = as.numeric(epi.intervalos[, 4]),
param.data = i.data,
param.seasons = i.seasons,
param.type.threshold = i.type.threshold,
param.level.threshold = i.level.threshold,
param.tails.threshold = i.tails.threshold,
param.type.intensity = i.type.intensity,
param.level.intensity = i.level.intensity,
param.tails.intensity = i.tails.intensity,
param.type.curve = i.type.curve,
param.level.curve = i.level.curve,
param.type.other = i.type.other,
param.level.other = i.level.other,
param.method = i.method,
param.param = i.param,
param.centering = i.centering,
param.n.max = i.n.max,
param.type.boot = i.type.boot,
param.iter.boot = i.iter.boot,
param.mem.info = i.mem.info
)
memmodel.output$call <- match.call()
class(memmodel.output) <- "mem"
}
}
return(memmodel.output)
}
#' @export
print.mem <- function(x, ...) {
threshold <- as.data.frame(round(t(x$pre.post.intervals[1:2, 3]), 2))
colnames(threshold) <- c("Pre", "Post")
rownames(threshold) <- "Threshold"
values <- as.data.frame(round(x$epi.intervals[, 4], 2))
colnames(values) <- "Threshold"
rownames(values) <- paste(c("Medium", "High", "Very high"), " (", round(x$epi.intervals[, 1] * 100, 1), "%)", sep = "")
cat("Call:\n")
print(x$call)
cat("\nEpidemic threshold:\n")
print(threshold)
cat("\nIntensity thresholds:\n")
print(values)
}
#' @export
summary.mem <- function(object, ...) {
threshold <- as.data.frame(round(t(object$pre.post.intervals[1:2, 3]), 2))
colnames(threshold) <- c("Pre", "Post")
rownames(threshold) <- "Threshold"
values <- as.data.frame(round(object$epi.intervals[, 4], 2))
colnames(values) <- "Threshold"
rownames(values) <- paste(c("Medium", "High", "Very high"), " (", round(object$epi.intervals[, 1] * 100, 1), "%)", sep = "")
cat("Call:\n", sep = "")
print(object$call)
cat("\nParameters:\n", sep = "")
cat("\t- General:\n", sep = "")
cat("\t\t+ Number of seasons restriction: ", if (object$param.seasons == -1 | is.na(object$param.seasons)) "Unrestricted (maximum)" else paste("Restricted to ", object$param.seasons, sep = ""), "\n", sep = "")
cat("\t\t+ Number of seasons used: ", object$n.seasons, "\n", sep = "")
cat("\t\t+ Seasons used: ", paste(names(object$data), collapse = ", "), "\n", sep = "")
cat("\t\t+ Number of weeks: ", object$n.weeks, "\n", sep = "")
cat("\t- Confidence intervals:\n", sep = "")
cat("\t\t+ Epidemic threshold: ", output.ci(object$param.type.threshold, object$param.level.threshold, object$param.tails.threshold), "\n", sep = "")
cat("\t\t+ Intensity: ", output.ci(object$param.type.intensity, object$param.level.intensity, object$param.tails.intensity), "\n", sep = "")
cat("\t\t+ Curve: ", output.ci(object$param.type.curve, object$param.level.curve, 2), "\n", sep = "")
cat("\t\t+ Others: ", output.ci(object$param.type.other, object$param.level.other, 2), "\n", sep = "")
cat("\t- Epidemic timing calculation:\n", sep = "")
cat("\t\t+ Method: ", object$param.method, "\n", sep = "")
cat("\t\t+ Parameter: ", object$param.param, "\n", sep = "")
cat("\t- Epidemic threshold calculation:\n", sep = "")
cat("\t\t+ Pre-epidemic values: ", if (is.na(object$param.n.max)) {
paste("Optimized: ", object$n.max, sep = "")
} else {
if (object$param.n.max == -1) paste("Optimized: ", object$n.max, sep = "") else object$n.max
}, "\n", sep = "")
cat("\t\t+ Tails of CI: ", object$param.tails.threshold, "\n", sep = "")
cat("\t- Intensity thresholds calculation:\n", sep = "")
cat("\t\t+ Number of values: ", if (is.na(object$param.n.max)) {
paste("Optimized: ", object$n.max, sep = "")
} else {
if (object$param.n.max == -1) paste("Optimized: ", object$n.max, sep = "") else object$n.max
}, "\n", sep = "")
cat("\t\t+ Tails of CI: ", object$param.tails.intensity, "\n", sep = "")
cat("\t\t+ Levels of CI: ", paste0(paste0(c("Medium", "High", "Very high"), ": ", round(object$epi.intervals[, 1] * 100, 1), "%"), collapse = ", "), "\n")
cat("\t- Bootstrap (if used):\n", sep = "")
cat("\t\t+ Technique: ", if (is.na(object$param.type.boot)) "-" else object$param.type.boot, "\n", sep = "")
cat("\t\t+ Bootstrap samples: ", if (is.na(object$param.iter.boot)) "-" else object$param.iter.boot, "\n", sep = "")
cat("\nEpidemic description:\n", sep = "")
cat("\t- Typical influenza epidemics starts at week ", round(object$ci.start[2, 2], 2), ". ", 100 * object$param.level.other, "% CI: [", round(object$ci.start[2, 1], 2), ",", round(object$ci.start[2, 3], 2), "]\n", sep = "")
cat("\t- Typical influenza season lasts ", round(object$ci.length[1, 2], 2), " weeks. ", 100 * object$param.level.other, "% CI: [", round(object$ci.length[1, 1], 2), ",", round(object$ci.length[1, 3], 2), "]\n", sep = "")
cat("\t- This optimal ", object$mean.length, " weeks influenza season includes the ", round(object$ci.percent[2], 2), "% of the total sum of rates\n\n", sep = "")
cat("\nEpidemic threshold:\n", sep = "")
print(threshold)
cat("\n\nIntensity thresholds:\n", sep = "")
print(values)
}
#' @export
plot.mem <- function(x, ...) {
# opar<-par(mfrow=c(1,2))
opar <- par(mfrow = c(1, 2), mar = c(4, 3, 1, 2) + 0.1, mgp = c(3, 0.5, 0), xpd = T)
# Graph 1
semanas <- dim(x$data)[1]
anios <- dim(x$data)[2]
names.seasons <- sub("^.*\\(([^\\)]*)\\)$", "\\1", names(x$data), perl = T)
datos.graf <- x$moving.epidemics
colnames(datos.graf) <- names(x$data)
# lab.graf<-(1:semanas)+x$mean.start[2]-x$mean.start[1]
lab.graf <- (1:semanas) + x$ci.start[2, 2] - x$ci.start[1, 2]
lab.graf[lab.graf > 52] <- (lab.graf - 52)[lab.graf > 52]
lab.graf[lab.graf < 1] <- (lab.graf + 52)[lab.graf < 1]
tipos <- rep(1, anios)
anchos <- rep(2, anios)
pal <- colorRampPalette(brewer.pal(4, "Spectral"))
colores <- pal(anios)
# colores<-c(rgb(runif(anios-1),runif(anios-1),runif(anios-1)),"#FF0000")
# limite.superior<-c(1.05*maxFixNA(datos.graf))
otick <- optimal.tickmarks(0, maxFixNA(datos.graf), 10)
range.y <- c(otick$range[1], otick$range[2] + otick$by / 2)
matplot(1:semanas,
datos.graf,
type = "l",
sub = paste("Weekly incidence rates of the seasons matching their relative position in the model."),
lty = tipos,
lwd = anchos,
col = colores,
xlim = c(1, semanas),
# xlab="Week",
# ylab="Rate",
xlab = "",
ylab = "",
axes = F,
# font.axis=1,
# font.lab=1,
font.main = 2,
font.sub = 1,
ylim = range.y
)
# Ejes
axis(1, at = seq(1, semanas, 1), labels = F, cex.axis = 0.7, col.axis = "#404040", col = "#C0C0C0")
axis(1,
at = seq(1, semanas, 2), tick = F,
labels = as.character(lab.graf)[seq(1, semanas, 2)], cex.axis = 0.7, col.axis = "#404040", col = "#C0C0C0"
)
axis(1,
at = seq(2, semanas, 2), tick = F,
labels = as.character(lab.graf)[seq(2, semanas, 2)], cex.axis = 0.7, line = 0.60, col.axis = "#404040", col = "#C0C0C0"
)
mtext(1, text = "Week", line = 2, cex = 0.8, col = "#000040")
axis(2, at = otick$tickmarks, lwd = 1, cex.axis = 0.6, col.axis = "#404040", col = "#C0C0C0")
mtext(2, text = "Weekly rate", line = 1.3, cex = 0.8, col = "#000040")
if (x$param.mem.info) mtext(4, text = paste("mem R library - Jose E. Lozano - https://github.com/lozalojo/mem", sep = ""), line = 0.75, cex = 0.6, col = "#404040")
i.temporada <- x$ci.start[1, 2]
f.temporada <- x$ci.start[1, 2] + x$mean.length - 1
lines(x = rep(i.temporada - 0.5, 2), y = range.y, col = "#FF0000", lty = 2, lwd = 2)
lines(x = rep(f.temporada + 0.5, 2), y = range.y, col = "#40FF40", lty = 2, lwd = 2)
lines(x = c(1, semanas), y = rep(x$epidemic.thresholds[1], 2), col = "#8c6bb1", lty = 2, lwd = 2)
# abline(v=c(i.temporada-0.5,f.temporada+0.5),col=c("#00C000","#FFB401"),lty=c(2,2))
# ya<-range.y[2]*0.975
ya <- otick$range[2]
if ((i.temporada - 1) <= (semanas - f.temporada)) xa <- f.temporada + 1 else xa <- 1
legend(
x = xa, y = ya, inset = c(0, 0), xjust = 0, seg.len = 1,
legend = names.seasons,
bty = "n",
lty = tipos,
lwd = anchos,
col = colores,
cex = 0.75,
x.intersp = 0.5,
y.intersp = 0.6,
text.col = "#000000",
ncol = 1
)
temp1 <- memsurveillance(
i.current = data.frame(rates = x$typ.curve[, 2], row.names = rownames(x$data)),
i.epidemic.thresholds = x$epidemic.thresholds,
i.intensity.thresholds = x$intensity.thresholds,
i.mean.length = x$mean.length,
i.force.length = T,
i.graph.file = F,
i.week.report = NA,
i.equal = F,
i.pos.epidemic = T,
i.no.epidemic = F,
i.no.intensity = F,
i.epidemic.start = x$ci.start[2, 2],
i.range.x = as.numeric(c(rownames(x$data)[1], rownames(x$data)[semanas])),
i.range.x.53 = F,
i.range.y = NA,
i.no.labels = F,
i.start.end.marks = T
)
par(opar)
}
lozalojo/mem documentation built on July 7, 2023, 10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plot.mem summary.mem print.mem memmodel
Documented in memmodel
#' @title Methods for influenza modelization
#'
#' @description
#' Function \code{memmodel} is used to calculate the threshold for influenza epidemic using historical
#' records (surveillance rates).
#'
#' The method to calculate the threshold is described in the Moving Epidemics Method (MEM) used to
#' monitor influenza activity in a weekly surveillance system.
#'
#' @name memmodel
#'
#' @aliases summary.flu,plot.flu,print.flu
#'
#' @param i.data Data frame of input data.
#' @param i.seasons Maximum number of seasons to use.
#' @param i.type.threshold Type of confidence interval to calculate the threshold.
#' @param i.level.threshold Level of confidence interval to calculate the threshold.
#' @param i.tails.threshold Tails for the confidence interval to calculate the threshold.
#' @param i.type.intensity Type of confidence interval to calculate the intensity thresholds.
#' @param i.level.intensity Levels of confidence interval to calculate the intensity thresholds.
#' @param i.tails.intensity Tails for the confidence interval to calculate the threshold.
#' @param i.type.curve Type of confidence interval to calculate the modelled curve.
#' @param i.level.curve Level of confidence interval to calculate the modelled curve.
#' @param i.type.other Type of confidence interval to calculate length, start and percentages.
#' @param i.level.other Level of confidence interval to calculate length, start and percentages.
#' @param i.method Method to calculate the optimal timing of the epidemic.
#' @param i.param Parameter to calculate the optimal timing of the epidemic.
#' @param i.centering Number of weeks to center the moving seasons.
#' @param i.n.max Number of pre-epidemic values used to calculate the threshold.
#' @param i.type.boot Type of bootstrap technique.
#' @param i.iter.boot Number of bootstrap iterations.
#' @param i.mem.info include information about the package in the graph.
#'
#' @return
#' \code{memmodel} returns an object of class \code{mem}.
#' An object of class \code{mem} is a list containing at least the following components:
#' \itemize{
#' \item{i.data} {input data}
#' \item{pre.post.intervals} {Pre/post confidence intervals (Threhold is the upper limit
#' of the confidence interval).}
#' \item{ci.length} {Mean epidemic length confidence interval.}
#' \item{ci.percent} {Mean covered percentage confidence interval.}
#' \item{mean.length} {Mean length.}
#' \item{moving.epidemics} {Moving epidemic rates.}
#' \item{mean.start} {Mean epidemic start.}
#' \item{epi.intervals} {Epidemic levels of intensity.}
#' \item{typ.curve} {Typical epidemic curve.}
#' \item{n.max} {Effective number of pre epidemic values.}
#' }
#'
#' @details
#' Input data is a data frame containing rates that represent historical influenza surveillance
#' data. It can start and end at any given week (tipically at week 40th), and rates can be
#' expressed as per 100,000 inhabitants (or per consultations, if population is not
#' available) or any other scale.
#'
#' Parameters \code{i.type}, \code{i.type.threshold} and \code{i.type.curve} defines how to
#' calculate confidence intervals along the process.
#'
#' \code{i.type.curve} is used for calculating the typical influenza curve,
#' \code{i.type.threshold} is used to calculate the pre and post epidemic threshold and
#' \code{i.type} is used for any other confidende interval used in the method.
#'
#' All three parameters must be a number between \code{1} and \code{6}:
#'
#' \itemize{
#' \item{1} {Arithmetic mean and mean confidence interval.}
#' \item{2} {Geometric mean and mean confidence interval.}
#' \item{3} {Median and the Nyblom/normal aproximation confidence interval.}
#' \item{4} {Median and bootstrap confidence interval.}
#' \item{5} {Arithmetic mean and point confidence interval (standard deviations).}
#' \item{6} {Geometric mean and point confidence interval (standard deviations).}
#' }
#'
#' Option \code{3} uses the Hettmansperger and Sheather (1986) and Nyblom (1992) method,
#' when there is enough sample size. If sample size is small, then the normal aproximation
#' will be used as described in Conover, 1980, p. 112. Refer to EnvStats package for
#' more information.
#'
#' Option \code{4} uses two more parameters: \code{i.type.boot} indicates which bootstrap
#' method to use. The values are the same of those of the \code{\link{boot.ci}} function.
#' Parameter \code{i.iter.boot} indicates the number of bootstrap samples to use. See
#' \code{\link{boot}} for more information about this topic.
#'
#' Parameters \code{i.level}, \code{i.level.threshold} and \code{i.level.curve} indicates,
#' respectively, the level of the confidence intervals described above.
#'
#' The \code{i.n.max} parameter indicates how many pre epidemic values to use to calculate
#' the threshold. A value of -1 indicates the program to use an appropiate number of points
#' depending on the number of seasons provided as input. \code{i.tails} tells the program
#' to use {1} or {2} tailed confidence intervals when calculating the threshold (1 is
#' recommended).
#'
#' Parameters \code{i.method} and \code{i.param} indicates how to find the optimal timing
#' of the epidemics. See \code{\link{memtiming}} for details on the values this parameters
#' can have.
#'
#' It is important to know how to arrange information in order to use with memapp. The key points are:
#'
#' \itemize{
#' \item One single epidemic wave each season.
#' \item Never delete a rate inside an epidemic.
#' \item Accommodate week 53.
#' \item Do not inflate missing values with zeroes.
#' }
#'
#' Data must contain information from the historical series. Surveillance period can start and end at
#' any given week (typically start at week 40th and ends at week 20th), and data can have any units and
#' can be expressed in any scale (typically rates per 100,000 inhabitants or consultations).
#'
#' The table must have one row per epidemiological week and one column per surveillance season. A season
#' is a full surveillance period from the beginning to the end, where occurs at some point one single
#' epidemic wave on it. No epidemic wave can be spared in two consecutive seasons. If so, you have to
#' redefine the start and end of the season defined in your dataset. If a season have two waves, it must
#' be split in two periods and must be named accordingly with the seasons name conventions described
#' below. Each cell contains the value for a given week in a given season.
#'
#' The first column should contain the names of the weeks. When the season contains two different calendar
#' years, the week will go from 40th of the first year to 52nd, and then from 1st to 20th. When the season
#' contains one year, the weeks will go from 1st to 52nd.
#'
#' Note: If there is no column with week names, the application will name the weeks numbering from 1 to
#' the number of rows.
#'
#' \itemize{
#' \item In the northern hemisphere countries, the surveillance period usually goes from
#' week 40 to 20 of the following year (notation: season 2016/2017).
#' \item In the southern hemisphere countries, the surveillance period usually goes from
#' week 18 to 39 same year (notation: season 2017).
#' }
#'
#' The first row must contain the names of the seasons. This application understand the naming of a
#' season when it contains one or two four digits year separated by / and one one-digit number
#' between parenthesis to identify the wave number. The wave number part in a name of a season is
#' used when a single surveillance period has two epidemic waves that have to be separated in order
#' to have reliable results. In this case, each wave is placed in different columns and named ending
#' with (1) for the first period, (2) for the second, and so on.
#'
#' @examples
#' # Castilla y Leon Influenza Rates data
#' data(flucyl)
#' # Finds the timing of the first season: 2001/2002
#' epi <- memmodel(flucyl)
#' print(epi)
#' summary(epi)
#' plot(epi)
#' @author Jose E. Lozano \email{lozalojo@@gmail.com}
#'
#' @references
#' Vega T, Lozano JE, Ortiz de Lejarazu R, Gutierrez Perez M. Modelling influenza epidemic - can we
#' detect the beginning and predict the intensity and duration? Int Congr Ser. 2004 Jun;1263:281-3.
#'
#' Vega T, Lozano JE, Meerhoff T, Snacken R, Mott J, Ortiz de Lejarazu R, et al. Influenza surveillance
#' in Europe: establishing epidemic thresholds by the moving epidemic method. Influenza Other Respir
#' Viruses. 2013 Jul;7(4):546-58. DOI:10.1111/j.1750-2659.2012.00422.x.
#'
#' Vega T, Lozano JE, Meerhoff T, Snacken R, Beaute J, Jorgensen P, et al. Influenza surveillance in
#' Europe: comparing intensity levels calculated using the moving epidemic method. Influenza Other
#' Respir Viruses. 2015 Sep;9(5):234-46. DOI:10.1111/irv.12330.
#'
#' Lozano JE. lozalojo/mem: Second release of the MEM R library. Zenodo [Internet]. [cited 2017 Feb 1];
#' Available from: \url{https://zenodo.org/record/165983}. DOI:10.5281/zenodo.165983
#'
#' Hettmansperger, T. P., and S. J Sheather. 1986. Confidence Intervals Based on Interpolated Order
#' Statistics. Statistics and Probability Letters 4: 75-79. doi:10.1016/0167-7152(86)90021-0.
#'
#' Nyblom, J. 1992. Note on Interpolated Order Statistics. Statistics and Probability Letters 14:
#' 129-31. doi:10.1016/0167-7152(92)90076-H.
#'
#' Conover, W.J. (1980). Practical Nonparametric Statistics. Second Edition. John Wiley and Sons, New York.
#'
#' @keywords influenza
#'
#' @export
#' @importFrom stats loess qt qnorm quantile runif var
#' @importFrom grDevices colorRampPalette colorRamp dev.off rgb tiff
#' @importFrom RColorBrewer brewer.pal
#' @importFrom graphics abline axis legend matplot mtext par points text lines
memmodel <- function(i.data,
i.seasons = 10,
i.type.threshold = 5,
i.level.threshold = 0.95,
i.tails.threshold = 1,
i.type.intensity = 6,
i.level.intensity = c(0.40, 0.90, 0.975),
i.tails.intensity = 1,
i.type.curve = 2,
i.level.curve = 0.95,
i.type.other = 3,
i.level.other = 0.95,
i.method = 2,
i.param = 2.8,
i.centering = -1,
i.n.max = -1,
i.type.boot = "norm",
i.iter.boot = 10000,
i.mem.info = T) {
if (is.null(dim(i.data))) {
memmodel.output <- NULL
cat("Incorrect number of dimensions, input must be a data.frame\n")
} else {
# Test if first column is week.name
if (all(i.data[,1] %in% 1:53)){
datos <- i.data[-1]
rownames(datos)<-i.data[,1]
}else{
datos <- i.data
if (!all(rownames(i.data) %in% 1:53)) rownames(datos)<-1:NROW(i.data)
}
datos <- datos[apply(datos, 2, function(x) sum(x, na.rm = T) > 0)]
if (!(ncol(datos) > 1)) {
memmodel.output <- NULL
cat("Incorrect number of dimensions, at least two seasons of valid data (non-zero values) required\n")
} else {
if (is.matrix(datos)) datos <- as.data.frame(datos)
if (is.null(i.seasons)) i.seasons <- NA
if (!is.na(i.seasons)) if (i.seasons > 0) datos <- datos[(max((dim(datos)[2]) - i.seasons + 1, 1)):(dim(datos)[2])]
datos <- as.data.frame(apply(datos, 2, fill.missing))
semanas <- dim(datos)[1]
anios <- NCOL(datos)
# Calcular el optimo
if (is.na(i.n.max)) {
n.max <- max(1, round(30 / anios, 0))
} else {
if (i.n.max == -1) {
n.max <- max(1, round(30 / anios, 0))
# if (anios>=10) n.max=3 else n.max=5
} else if (i.n.max == 0) {
n.max <- semanas
} else {
n.max <- i.n.max
}
}
if (is.null(i.method)) i.method <- 2
if (is.na(i.method)) i.method <- 2
if (is.null(i.param)) i.param <- 2.8
if (is.na(i.param)) i.param <- 2.8
optimo <- apply(datos, 2, memtiming, i.n.values = n.max, i.method = i.method, i.param = i.param)
datos.duracion.real <- extraer.datos.optimo.map(optimo)
# por defecto la aritmetica en ambos
ic.duracion <- iconfianza(as.numeric(datos.duracion.real[1, ]), i.level.other, i.type.other, T, i.type.boot, i.iter.boot, 2)
ic.duracion <- rbind(ic.duracion, c(floor(ic.duracion[1]), round(ic.duracion[2]), ceiling(ic.duracion[3])))
ic.porcentaje <- iconfianza(as.numeric(datos.duracion.real[2, ]), i.level.other, i.type.other, T, i.type.boot, i.iter.boot, 2)
duracion.media <- ic.duracion[2, 2]
if (is.na(i.centering) | i.centering == -1) duracion.centrado <- duracion.media else duracion.centrado <- i.centering
########################################################################################
# Calcula todos los parametros asociados a la temporada gripal en base a la estimacion #
# del numero de semanas de duracion total. #
########################################################################################
# Datos de la gripe en la temporada REAL y en la temporada MEDIA, esta ultima sirve para unir las temporadas.
semana.inicio <- as.integer(rownames(datos)[1])
gripe <- array(dim = c(7, anios, 3), dimnames = list(
"indicators" = c("epidemic start", "epidemic end", "sum of epidemic values", "sum of values", "covered percentage", "sum of pre-epidemic values", "sum of post-epidemic values"),
"seasons" = names(datos),
"type" = c("real epidemic", "mean length epidemic", "centered length epidemic")
))
datos.duracion.media <- extraer.datos.curva.map(optimo, duracion.media)
datos.duracion.centrado <- extraer.datos.curva.map(optimo, duracion.centrado)
## Estimacion del numero de semana en el comienza la gripe
ic.inicio <- iconfianza(as.numeric(datos.duracion.real[4, ]), i.level.other, i.type.other, T, i.type.boot, i.iter.boot, 2)
ic.inicio <- rbind(ic.inicio, c(floor(ic.inicio[1]), round(ic.inicio[2]), ceiling(ic.inicio[3])))
ic.inicio <- rbind(ic.inicio, c(semana.absoluta(ic.inicio[2, 1], semana.inicio), semana.absoluta(ic.inicio[2, 2], semana.inicio), semana.absoluta(ic.inicio[2, 3], semana.inicio)))
ic.inicio <- ic.inicio[c(2, 3, 1), ]
inicio.medio <- ic.inicio[1, 2]
ic.inicio.centrado <- iconfianza(as.numeric(datos.duracion.centrado[4, ]), i.level.other, i.type.other, T, i.type.boot, i.iter.boot, 2)
ic.inicio.centrado <- rbind(ic.inicio.centrado, c(floor(ic.inicio.centrado[1]), round(ic.inicio.centrado[2]), ceiling(ic.inicio.centrado[3])))
ic.inicio.centrado <- rbind(ic.inicio.centrado, c(semana.absoluta(ic.inicio.centrado[2, 1], semana.inicio), semana.absoluta(ic.inicio.centrado[2, 2], semana.inicio), semana.absoluta(ic.inicio.centrado[2, 3], semana.inicio)))
ic.inicio.centrado <- ic.inicio.centrado[c(2, 3, 1), ]
inicio.centrado <- ic.inicio.centrado[1, 2]
for (j in 1:anios) {
## Semana en la q comienza la temporada.real
gripe[1, j, 1] <- datos.duracion.real[4, j]
## Semana en la q acaba
gripe[2, j, 1] <- datos.duracion.real[5, j]
## numero de casos en la temporada.real de gripe
gripe[3, j, 1] <- datos.duracion.real[3, j]
## Numero de casos totales de la temporada global
gripe[4, j, 1] <- sum(datos[, j], na.rm = TRUE)
## Porcentaje cubierto de casos de la temporada gripal
gripe[5, j, 1] <- datos.duracion.real[2, j]
## Casos Acumulados justo antes de comenzar la temporada gripal
if (gripe[1, j, 1] > 1) gripe[6, j, 1] <- sum(datos[1:(gripe[1, j, 1] - 1), j], na.rm = TRUE) else gripe[6, j, 1] <- 0
## Casos Despues de la temporada gripal
if (gripe[2, j, 1] < semanas) gripe[7, j, 1] <- sum(datos[(gripe[2, j, 1] + 1):semanas, j], na.rm = TRUE) else gripe[7, j, 1] <- 0
gripe[1, j, 2] <- datos.duracion.media[4, j]
gripe[2, j, 2] <- datos.duracion.media[5, j]
gripe[3, j, 2] <- datos.duracion.media[3, j]
gripe[4, j, 2] <- sum(datos[, j], na.rm = TRUE)
gripe[5, j, 2] <- datos.duracion.media[2, j]
if (gripe[1, j, 2] > 1) gripe[6, j, 2] <- sum(datos[1:(gripe[1, j, 2] - 1), j], na.rm = TRUE) else gripe[6, j, 2] <- 0
if (gripe[2, j, 2] < semanas) gripe[7, j, 2] <- sum(datos[(gripe[2, j, 2] + 1):semanas, j], na.rm = TRUE) else gripe[7, j, 2] <- 0
gripe[1, j, 3] <- datos.duracion.centrado[4, j]
gripe[2, j, 3] <- datos.duracion.centrado[5, j]
gripe[3, j, 3] <- datos.duracion.centrado[3, j]
gripe[4, j, 3] <- sum(datos[, j], na.rm = TRUE)
gripe[5, j, 3] <- datos.duracion.centrado[2, j]
if (gripe[1, j, 3] > 1) gripe[6, j, 3] <- sum(datos[1:(gripe[1, j, 3] - 1), j], na.rm = TRUE) else gripe[6, j, 3] <- 0
if (gripe[2, j, 3] < semanas) gripe[7, j, 3] <- sum(datos[(gripe[2, j, 3] + 1):semanas, j], na.rm = TRUE) else gripe[7, j, 3] <- 0
}
# Indice de estado temporal. Si es un 1 me dice que estamos en epoca pretemporada,
# si pone un 2 indica que estamos en temporada gripal. Si pone un 3 indica que
# estamos en posttemporada.
indices.temporada <- array(dim = c(semanas, anios, 3), dimnames = list(
"weeks" = rownames(datos),
"seasons" = names(datos),
"type" = c("real epidemic", "mean length epidemic", "centered length epidemic")
))
for (j in 1:anios) {
indices.temporada[-((datos.duracion.real[4, j]):semanas), j, 1] <- 1
indices.temporada[(datos.duracion.real[4, j]):(datos.duracion.real[5, j]), j, 1] <- 2
indices.temporada[-(1:(datos.duracion.real[5, j])), j, 1] <- 3
indices.temporada[-((datos.duracion.media[4, j]):semanas), j, 2] <- 1
indices.temporada[(datos.duracion.media[4, j]):(datos.duracion.media[5, j]), j, 2] <- 2
indices.temporada[-(1:(datos.duracion.media[5, j])), j, 2] <- 3
indices.temporada[-((datos.duracion.centrado[4, j]):semanas), j, 3] <- 1
indices.temporada[(datos.duracion.centrado[4, j]):(datos.duracion.centrado[5, j]), j, 3] <- 2
indices.temporada[-(1:(datos.duracion.centrado[5, j])), j, 3] <- 3
}
# Grafico del esquema de las temporadas gripales para su union
longitud.esquema <- semanas + maxFixNA(datos.duracion.centrado[4, ]) - minFixNA(datos.duracion.centrado[4, ])
inicio.epidemia.esquema <- maxFixNA(datos.duracion.centrado[4, ])
fin.epidemia.esquema <- maxFixNA(datos.duracion.centrado[5, ])
esquema.temporadas <- array(dim = c(longitud.esquema, anios + 2, 3), dimnames = list(
"dummy weeks" = as.character(1:longitud.esquema),
"seasons" = c(names(datos), "mean season week", "mean season index"),
"type" = c("absolute week", "relative week", "values")
))
for (j in 1:anios) {
diferencia.anual <- inicio.epidemia.esquema - datos.duracion.centrado[4, j]
for (i in 1:semanas) {
esquema.temporadas[i + diferencia.anual, j, 1] <- i
esquema.temporadas[i + diferencia.anual, j, 2] <- semana.absoluta(i, semana.inicio)
esquema.temporadas[i + diferencia.anual, j, 3] <- datos[i, j]
}
}
diferencia.global <- inicio.epidemia.esquema - inicio.centrado
for (i in 1:semanas) {
if ((i + diferencia.global) >= 1 & (i + diferencia.global) <= longitud.esquema) {
esquema.temporadas[i + diferencia.global, anios + 1, 1] <- i
esquema.temporadas[i + diferencia.global, anios + 1, 2] <- semana.absoluta(i, semana.inicio)
esquema.temporadas[i + diferencia.global, anios + 1, 3] <- i
}
}
esquema.temporadas[, anios + 2, ] <- 3
esquema.temporadas[1:(diferencia.global + inicio.medio - 1), anios + 2, ] <- 1
esquema.temporadas[(diferencia.global + inicio.medio):(diferencia.global + inicio.medio + duracion.media - 1), anios + 2, ] <- 2
# esquema.temporadas[1:(inicio.epidemia.esquema-1),anios+2,2]<-1
# esquema.temporadas[(inicio.epidemia.esquema):(inicio.epidemia.esquema+duracion.media-1),anios+2,2]<-2
## Temporadas moviles
temporadas.moviles <- esquema.temporadas[, c(1:anios), 3]
## Limites de la temporada
# limites.temporada<-c(inicio.medio,inicio.medio+duracion.media-1)
# limites.temporada<-rbind(limites.temporada,c(semana.absoluta(inicio.medio,semana.inicio),semana.absoluta(inicio.medio+duracion.media-1,semana.inicio)))
# Limites relativos
# limites.esquema<-c(inicio.epidemia.esquema,inicio.epidemia.esquema+duracion.media-1)
## En ttotal reordenamos la gripe para q todas las temporadas comiencen la semana "estinicio" (es decir, quitamos informacion
## por delante y por detras de esquema.temporada[,,3] para que queden 35 semanas (las que hay).
temporadas.moviles.recortada <- array(dim = c(semanas, anios), dimnames = list(
"dummy weeks" = as.character(1:semanas),
"seasons" = names(datos)
))
## Temporada gripal comienza "estinicio" y termina en "estinicio+temporada-1"
for (j in 1:anios) {
for (i in 1:semanas) {
if ((i + diferencia.global) >= 1 & (i + diferencia.global) <= longitud.esquema) {
temporadas.moviles.recortada[i, j] <- temporadas.moviles[i + diferencia.global, j]
} else {
temporadas.moviles.recortada[i, j] <- NA
}
}
}
## Dibujamos la curva epidemica tipo
# Importante: ?Quitamos los ceros?
# iy<-(is.na(temporadas.moviles.recortada) | temporadas.moviles.recortada==0)
iy <- is.na(temporadas.moviles.recortada)
temporadas.moviles.recortada.no.ceros <- temporadas.moviles.recortada
temporadas.moviles.recortada.no.ceros[iy] <- NA
curva.tipo <- t(apply(temporadas.moviles.recortada.no.ceros, 1, iconfianza, nivel = i.level.curve, tipo = i.type.curve, ic = T, tipo.boot = i.type.boot, iteraciones.boot = i.iter.boot, colas = 2))
## Seleccionamos los periodos pre, epidemia, y post
## PRE y POST-TEMPORADA GRIPAL
## Como no se registra mas q de la semana 40 a la 20, tenemos q muchas de las semanas tienen tasa 0, eliminamos esas semanas,
## ya q podrian llevar a una infraestimacion de la tasa base fuera de temporada
pre.datos <- as.vector(as.matrix(extraer.datos.pre.epi(optimo)))
post.datos <- as.vector(as.matrix(extraer.datos.post.epi(optimo)))
epi.datos <- as.vector(as.matrix(extraer.datos.epi(optimo)))
epi.datos.2 <- as.vector(as.matrix(apply(datos, 2, maxnvalores, n.max = n.max)))
pre.post.datos <- rbind(pre.datos, post.datos)
# IC de la linea basica de pre y post temporada
# por defecto estaba la geometrica
# pre.d<-pre.datos[!(is.na(pre.datos) | pre.datos==0)]
# post.d<-post.datos[!(is.na(post.datos) | post.datos==0)]
# epi.d<-epi.datos[!(is.na(epi.datos) | epi.datos==0)]
# epi.d.2<-epi.datos.2[!(is.na(epi.datos.2) | epi.datos.2==0)]
pre.d <- pre.datos[!is.na(pre.datos)]
post.d <- post.datos[!is.na(post.datos)]
epi.d <- epi.datos[!is.na(epi.datos)]
epi.d.2 <- epi.datos.2[!is.na(epi.datos.2)]
pre.i <- iconfianza(pre.d, nivel = i.level.threshold, tipo = i.type.threshold, ic = T, tipo.boot = i.type.boot, iteraciones.boot = i.iter.boot, colas = i.tails.threshold)
post.i <- iconfianza(post.d, nivel = i.level.threshold, tipo = i.type.threshold, ic = T, tipo.boot = i.type.boot, iteraciones.boot = i.iter.boot, colas = i.tails.threshold)
epi.intervalos <- numeric()
for (niv in i.level.intensity) epi.intervalos <- rbind(epi.intervalos, c(niv, iconfianza(epi.d, nivel = niv, tipo = i.type.intensity, ic = T, colas = i.tails.intensity)))
epi.intervalos.2 <- numeric()
for (niv in i.level.intensity) epi.intervalos.2 <- rbind(epi.intervalos.2, c(niv, iconfianza(epi.d.2, nivel = niv, tipo = i.type.intensity, ic = T, colas = i.tails.intensity)))
pre.post.intervalos <- rbind(pre.i, post.i)
memmodel.output <- list(
pre.post.intervals = pre.post.intervalos,
ci.length = ic.duracion,
ci.percent = ic.porcentaje,
ci.start = ic.inicio,
mean.length = duracion.media,
mean.start = inicio.medio,
centered.length = duracion.centrado,
centered.start = inicio.centrado,
moving.epidemics = temporadas.moviles.recortada,
epi.intervals = epi.intervalos,
epi.intervals.full = epi.intervalos.2,
typ.curve = curva.tipo,
n.max = n.max,
n.seasons = anios,
n.weeks = semanas,
data = datos,
start.week = semana.inicio,
epi.data = epi.datos,
epi.data.nona = epi.d,
epi.data.full = epi.datos.2,
epi.data.nona.full = epi.d.2,
optimum = optimo,
real.length.data = datos.duracion.real,
seasons.data = gripe,
season.indexes = indices.temporada,
season.scheme = esquema.temporadas,
seasons.moving = temporadas.moviles,
pre.post.data = pre.post.datos,
pre.data.nona = pre.d,
post.data.nona = post.d,
epidemic.thresholds = as.numeric(pre.post.intervalos[, 3]),
intensity.thresholds = as.numeric(epi.intervalos[, 4]),
param.data = i.data,
param.seasons = i.seasons,
param.type.threshold = i.type.threshold,
param.level.threshold = i.level.threshold,
param.tails.threshold = i.tails.threshold,
param.type.intensity = i.type.intensity,
param.level.intensity = i.level.intensity,
param.tails.intensity = i.tails.intensity,
param.type.curve = i.type.curve,
param.level.curve = i.level.curve,
param.type.other = i.type.other,
param.level.other = i.level.other,
param.method = i.method,
param.param = i.param,
param.centering = i.centering,
param.n.max = i.n.max,
param.type.boot = i.type.boot,
param.iter.boot = i.iter.boot,
param.mem.info = i.mem.info
)
memmodel.output$call <- match.call()
class(memmodel.output) <- "mem"
}
}
return(memmodel.output)
}
#' @export
print.mem <- function(x, ...) {
threshold <- as.data.frame(round(t(x$pre.post.intervals[1:2, 3]), 2))
colnames(threshold) <- c("Pre", "Post")
rownames(threshold) <- "Threshold"
values <- as.data.frame(round(x$epi.intervals[, 4], 2))
colnames(values) <- "Threshold"
rownames(values) <- paste(c("Medium", "High", "Very high"), " (", round(x$epi.intervals[, 1] * 100, 1), "%)", sep = "")
cat("Call:\n")
print(x$call)
cat("\nEpidemic threshold:\n")
print(threshold)
cat("\nIntensity thresholds:\n")
print(values)
}
#' @export
summary.mem <- function(object, ...) {
threshold <- as.data.frame(round(t(object$pre.post.intervals[1:2, 3]), 2))
colnames(threshold) <- c("Pre", "Post")
rownames(threshold) <- "Threshold"
values <- as.data.frame(round(object$epi.intervals[, 4], 2))
colnames(values) <- "Threshold"
rownames(values) <- paste(c("Medium", "High", "Very high"), " (", round(object$epi.intervals[, 1] * 100, 1), "%)", sep = "")
cat("Call:\n", sep = "")
print(object$call)
cat("\nParameters:\n", sep = "")
cat("\t- General:\n", sep = "")
cat("\t\t+ Number of seasons restriction: ", if (object$param.seasons == -1 | is.na(object$param.seasons)) "Unrestricted (maximum)" else paste("Restricted to ", object$param.seasons, sep = ""), "\n", sep = "")
cat("\t\t+ Number of seasons used: ", object$n.seasons, "\n", sep = "")
cat("\t\t+ Seasons used: ", paste(names(object$data), collapse = ", "), "\n", sep = "")
cat("\t\t+ Number of weeks: ", object$n.weeks, "\n", sep = "")
cat("\t- Confidence intervals:\n", sep = "")
cat("\t\t+ Epidemic threshold: ", output.ci(object$param.type.threshold, object$param.level.threshold, object$param.tails.threshold), "\n", sep = "")
cat("\t\t+ Intensity: ", output.ci(object$param.type.intensity, object$param.level.intensity, object$param.tails.intensity), "\n", sep = "")
cat("\t\t+ Curve: ", output.ci(object$param.type.curve, object$param.level.curve, 2), "\n", sep = "")
cat("\t\t+ Others: ", output.ci(object$param.type.other, object$param.level.other, 2), "\n", sep = "")
cat("\t- Epidemic timing calculation:\n", sep = "")
cat("\t\t+ Method: ", object$param.method, "\n", sep = "")
cat("\t\t+ Parameter: ", object$param.param, "\n", sep = "")
cat("\t- Epidemic threshold calculation:\n", sep = "")
cat("\t\t+ Pre-epidemic values: ", if (is.na(object$param.n.max)) {
paste("Optimized: ", object$n.max, sep = "")
} else {
if (object$param.n.max == -1) paste("Optimized: ", object$n.max, sep = "") else object$n.max
}, "\n", sep = "")
cat("\t\t+ Tails of CI: ", object$param.tails.threshold, "\n", sep = "")
cat("\t- Intensity thresholds calculation:\n", sep = "")
cat("\t\t+ Number of values: ", if (is.na(object$param.n.max)) {
paste("Optimized: ", object$n.max, sep = "")
} else {
if (object$param.n.max == -1) paste("Optimized: ", object$n.max, sep = "") else object$n.max
}, "\n", sep = "")
cat("\t\t+ Tails of CI: ", object$param.tails.intensity, "\n", sep = "")
cat("\t\t+ Levels of CI: ", paste0(paste0(c("Medium", "High", "Very high"), ": ", round(object$epi.intervals[, 1] * 100, 1), "%"), collapse = ", "), "\n")
cat("\t- Bootstrap (if used):\n", sep = "")
cat("\t\t+ Technique: ", if (is.na(object$param.type.boot)) "-" else object$param.type.boot, "\n", sep = "")
cat("\t\t+ Bootstrap samples: ", if (is.na(object$param.iter.boot)) "-" else object$param.iter.boot, "\n", sep = "")
cat("\nEpidemic description:\n", sep = "")
cat("\t- Typical influenza epidemics starts at week ", round(object$ci.start[2, 2], 2), ". ", 100 * object$param.level.other, "% CI: [", round(object$ci.start[2, 1], 2), ",", round(object$ci.start[2, 3], 2), "]\n", sep = "")
cat("\t- Typical influenza season lasts ", round(object$ci.length[1, 2], 2), " weeks. ", 100 * object$param.level.other, "% CI: [", round(object$ci.length[1, 1], 2), ",", round(object$ci.length[1, 3], 2), "]\n", sep = "")
cat("\t- This optimal ", object$mean.length, " weeks influenza season includes the ", round(object$ci.percent[2], 2), "% of the total sum of rates\n\n", sep = "")
cat("\nEpidemic threshold:\n", sep = "")
print(threshold)
cat("\n\nIntensity thresholds:\n", sep = "")
print(values)
}
#' @export
plot.mem <- function(x, ...) {
# opar<-par(mfrow=c(1,2))
opar <- par(mfrow = c(1, 2), mar = c(4, 3, 1, 2) + 0.1, mgp = c(3, 0.5, 0), xpd = T)
# Graph 1
semanas <- dim(x$data)[1]
anios <- dim(x$data)[2]
names.seasons <- sub("^.*\\(([^\\)]*)\\)$", "\\1", names(x$data), perl = T)
datos.graf <- x$moving.epidemics
colnames(datos.graf) <- names(x$data)
# lab.graf<-(1:semanas)+x$mean.start[2]-x$mean.start[1]
lab.graf <- (1:semanas) + x$ci.start[2, 2] - x$ci.start[1, 2]
lab.graf[lab.graf > 52] <- (lab.graf - 52)[lab.graf > 52]
lab.graf[lab.graf < 1] <- (lab.graf + 52)[lab.graf < 1]
tipos <- rep(1, anios)
anchos <- rep(2, anios)
pal <- colorRampPalette(brewer.pal(4, "Spectral"))
colores <- pal(anios)
# colores<-c(rgb(runif(anios-1),runif(anios-1),runif(anios-1)),"#FF0000")
# limite.superior<-c(1.05*maxFixNA(datos.graf))
otick <- optimal.tickmarks(0, maxFixNA(datos.graf), 10)
range.y <- c(otick$range[1], otick$range[2] + otick$by / 2)
matplot(1:semanas,
datos.graf,
type = "l",
sub = paste("Weekly incidence rates of the seasons matching their relative position in the model."),
lty = tipos,
lwd = anchos,
col = colores,
xlim = c(1, semanas),
# xlab="Week",
# ylab="Rate",
xlab = "",
ylab = "",
axes = F,
# font.axis=1,
# font.lab=1,
font.main = 2,
font.sub = 1,
ylim = range.y
)
# Ejes
axis(1, at = seq(1, semanas, 1), labels = F, cex.axis = 0.7, col.axis = "#404040", col = "#C0C0C0")
axis(1,
at = seq(1, semanas, 2), tick = F,
labels = as.character(lab.graf)[seq(1, semanas, 2)], cex.axis = 0.7, col.axis = "#404040", col = "#C0C0C0"
)
axis(1,
at = seq(2, semanas, 2), tick = F,
labels = as.character(lab.graf)[seq(2, semanas, 2)], cex.axis = 0.7, line = 0.60, col.axis = "#404040", col = "#C0C0C0"
)
mtext(1, text = "Week", line = 2, cex = 0.8, col = "#000040")
axis(2, at = otick$tickmarks, lwd = 1, cex.axis = 0.6, col.axis = "#404040", col = "#C0C0C0")
mtext(2, text = "Weekly rate", line = 1.3, cex = 0.8, col = "#000040")
if (x$param.mem.info) mtext(4, text = paste("mem R library - Jose E. Lozano - https://github.com/lozalojo/mem", sep = ""), line = 0.75, cex = 0.6, col = "#404040")
i.temporada <- x$ci.start[1, 2]
f.temporada <- x$ci.start[1, 2] + x$mean.length - 1
lines(x = rep(i.temporada - 0.5, 2), y = range.y, col = "#FF0000", lty = 2, lwd = 2)
lines(x = rep(f.temporada + 0.5, 2), y = range.y, col = "#40FF40", lty = 2, lwd = 2)
lines(x = c(1, semanas), y = rep(x$epidemic.thresholds[1], 2), col = "#8c6bb1", lty = 2, lwd = 2)
# abline(v=c(i.temporada-0.5,f.temporada+0.5),col=c("#00C000","#FFB401"),lty=c(2,2))
# ya<-range.y[2]*0.975
ya <- otick$range[2]
if ((i.temporada - 1) <= (semanas - f.temporada)) xa <- f.temporada + 1 else xa <- 1
legend(
x = xa, y = ya, inset = c(0, 0), xjust = 0, seg.len = 1,
legend = names.seasons,
bty = "n",
lty = tipos,
lwd = anchos,
col = colores,
cex = 0.75,
x.intersp = 0.5,
y.intersp = 0.6,
text.col = "#000000",
ncol = 1
)
temp1 <- memsurveillance(
i.current = data.frame(rates = x$typ.curve[, 2], row.names = rownames(x$data)),
i.epidemic.thresholds = x$epidemic.thresholds,
i.intensity.thresholds = x$intensity.thresholds,
i.mean.length = x$mean.length,
i.force.length = T,
i.graph.file = F,
i.week.report = NA,
i.equal = F,
i.pos.epidemic = T,
i.no.epidemic = F,
i.no.intensity = F,
i.epidemic.start = x$ci.start[2, 2],
i.range.x = as.numeric(c(rownames(x$data)[1], rownames(x$data)[semanas])),
i.range.x.53 = F,
i.range.y = NA,
i.no.labels = F,
i.start.end.marks = T
)
par(opar)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.