Nothing
R/RNoughtIndices.R
In CSIndicators: Climate Services' Indicators Based on Sub-Seasonal to Decadal Predictions
Defines functions
.empirical
.wesolowski
.mordecai
.liuhelmerssohn
.caminade
RNoughtIndices
CST_RNoughtIndices
Documented in
CST_RNoughtIndices
RNoughtIndices
#'R0 Index computation on s2dv_cube objects
#'
#'@author Javier Corvillo, \email{javier.corvillo@bsc.es}
#'@description This function computes the environmental contribution to Aedes-borne
#' disease transmissibility. Values higher than 1 imply that the environmental conditions allow for
#' the disease to be transmitted to more than 1 person (the infection can spread) while R0 values lower
#' than 1 imply that the environmental conditions are not suitable for the disease to spread.
#'
#' This function utilizes four possible ento-epidemiological models, each of them stated
#' in Caminade et al., 2015; Liu-Helmerssohn et al., 2014; and Mordecai et
#' al., 2017 and Wesolowski et al., 2015 respectively. Additionally, an adjustment to
#' real life data, using recorded DENV data recorded between 2014-2017 in the Caribbean, can also be performed.
#'
#'@references Caminade, Cyril et al. (Jan. 2017). ‘Global risk model for vector-borne
#' transmission of Zika virus reveals the role of El Niño 2015’.
#' In: Proceedings of the National Academy of Sciences 114.1.
#' Publisher: Proceedings of the National Academy of
#' Sciences, pp. 119–124. doi: 10.1073/pnas.1614303114.
#' url: https://www.pnas.org/doi/full/10.1073/pnas.1614303114.
#'
#'@references Liu-Helmerssohn, Jing et al. (Mar. 2014). ‘Vectorial Capacity of Aedes
#' aegypti: Effects of Temperature and Implications for Global Dengue Epidemic Potential’. en.
#' In: PLOS ONE 9.3.
#' Publisher: Public Library of Science, e89783.
#' issn: 1932-6203. doi: 10.1371/journal.pone.0089783.
#' url: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0089783.
#'
#' @references Mordecai, Erin A. et al. (Apr. 2017). ‘Detecting the impact of temperature on
#' transmission of Zika, dengue, and chikungunya using mechanistic models’. en.
#' In: PLOS Neglected Tropical Diseases 11.4.
#' Publisher: Public Library of Science, e0005568.
#' issn: 1935-2735. doi: 10.1371/journal.pntd.0005568.
#' url: https://journals.plos.org/plosntds/article?id=10.1371/journal.pntd.0005568.
#'
#' @references Wesolowski, Amy et al. (Sept. 2015). ‘Impact of human mobility on the emergence
#' of dengue epidemics in Pakistan’.
#' In: Proceedings of the National Academy of Sciences 112.38.
#' Publisher: Proceedings of the National Academy of Sciences, pp. 11887–11892.
#' doi: 10.1073/pnas.1504964112.
#' url: https://www.pnas.org/doi/full/10.1073/pnas.1504964112.
#'
#'@param temp An s2dv_cube object with temperature values, expressed in degrees Celsius. If
#' "caminade" is selected in "method", named spatial dimensions are required in the s2dv_cube object,
#' ordered as ("latitude", "longitude"). Valid dimension names are "lon", "longitude", "lat", "latitude".
#'
#' @param method a string indicating the R0 ento-epidemiological
#' model used to obtain the R0 outputs. Methods include "caminade",
#' "wesolowski", "liuhelmerssohn" and "mordecai". Alternatively, an
#' extrapolation to real life data can be made by setting "method" to "empirical".
#'
#' @param mm A Kramer probability matrix of dimensions (2, longitude, latitude) to be passed
#' if "caminade" is selected as the computation method. The spatial dimensions must
#' match those of temp and be equal in length. Valid dimension names are "lon", "longitude",
#' "lat", "latitude".
#'
#' @param lon_dim a character vector indicating the longitude dimension name in
#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
#'
#' @param lat_dim a character vector indicating the latitude dimension name in
#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
#'
#'@param ncores An integer indicating the number of cores to use in parallel
#' computation for temporal subsetting.
#'@return An s2dv_cube object containing the R0 Indices (unitless).
#'
#'@examples
#'dims <- c(time = 100, lat = 5, lon = 4)
#'temp_cube <- list(
#'data = array(
#' runif(prod(dims), min = 15, max = 35),
#' dim = dims
#' ),
#' attrs = list(
#' Variable = list(varName = '2m Temperature')
#' )
#' )
#' class(temp_cube) <- 's2dv_cube'
#'
#' R0 <- CST_RNoughtIndices(temp_cube, method = "mordecai", mm = NULL,
#' lon_dim = NULL, lat_dim = NULL, ncores = NULL)
#'
#' # Caminade method (requires Kramer matrix)
#' mm <- array(runif(2 * 5 * 4, 0.1, 1.0),
#' dim = c(2, 4, 5))
#' names(dim(mm)) <- c("", "lon", "lat")
#' R0 <- CST_RNoughtIndices(temp_cube, method = "caminade", mm = mm,
#' lon_dim = "lon", lat_dim = "lat", ncores = NULL)
#'
#'@import multiApply
#'@export
CST_RNoughtIndices <- function(temp, method, mm, lon_dim, lat_dim, ncores = NULL) {
# Check "s2dv_cube"
if (!inherits(temp, "s2dv_cube")) {
stop("Parameter 'temp' must be of the class 's2dv_cube'.")
}
RNought <- RNoughtIndices(temp = temp$data, method = method, mm = mm,
lon_dim = lon_dim, lat_dim = lat_dim,
ncores = ncores)
temp$data <- RNought
if ("Variable" %in% names(temp$attrs)) {
if ("varName" %in% names(temp$attrs$Variable)) {
temp$attrs$Variable$varName <- "RNoughtIndex"
}
}
return(temp)
}
#'R0 Index computation on s2dv_cube objects
#'
#'@author Javier Corvillo, \email{javier.corvillo@bsc.es}
#'@description This function computes the environmental contribution to Aedes-borne
#' disease transmissibility. Values higher than 1 imply that the environmental conditions allow for
#' the disease to be transmitted to more than 1 person (the infection can spread) while R0 values lower
#' than 1 imply that the environmental conditions are not suitable for the disease to spread.
#'
#' This function utilizes four possible ento-epidemiological models, each of them stated
#' in Caminade et al., 2015; Liu-Helmerssohn et al., 2014; and Mordecai et
#' al., 2017 and Wesolowski et al., 2015 respectively. Additionally, an adjustment to
#' real life data, using recorded DENV data recorded between 2014-2017 in the Caribbean, can also be performed.
#'
#'@references Caminade, Cyril et al. (Jan. 2017). ‘Global risk model for vector-borne
#' transmission of Zika virus reveals the role of El Niño 2015’.
#' In: Proceedings of the National Academy of Sciences 114.1.
#' Publisher: Proceedings of the National Academy of
#' Sciences, pp. 119–124. doi: 10.1073/pnas.1614303114.
#' url: https://www.pnas.org/doi/full/10.1073/pnas.1614303114.
#'
#'@references Liu-Helmerssohn, Jing et al. (Mar. 2014). ‘Vectorial Capacity of Aedes
#' aegypti: Effects of Temperature and Implications for Global Dengue Epidemic Potential’. en.
#' In: PLOS ONE 9.3.
#' Publisher: Public Library of Science, e89783.
#' issn: 1932-6203. doi: 10.1371/journal.pone.0089783.
#' url: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0089783.
#'
#' @references Mordecai, Erin A. et al. (Apr. 2017). ‘Detecting the impact of temperature on
#' transmission of Zika, dengue, and chikungunya using mechanistic models’. en.
#' In: PLOS Neglected Tropical Diseases 11.4.
#' Publisher: Public Library of Science, e0005568.
#' issn: 1935-2735. doi: 10.1371/journal.pntd.0005568.
#' url: https://journals.plos.org/plosntds/article?id=10.1371/journal.pntd.0005568.
#'
#' @references Wesolowski, Amy et al. (Sept. 2015). ‘Impact of human mobility on the emergence
#' of dengue epidemics in Pakistan’.
#' In: Proceedings of the National Academy of Sciences 112.38.
#' Publisher: Proceedings of the National Academy of Sciences, pp. 11887–11892.
#' doi: 10.1073/pnas.1504964112.
#' url: https://www.pnas.org/doi/full/10.1073/pnas.1504964112.
#'
#'@param temp An array with temperature values, expressed in degrees Celsius. If
#' "caminade" is selected in "method", named spatial dimensions are required in the s2dv_cube object,
#' ordered as ("latitude", "longitude"). Valid dimension names are "lon", "longitude", "lat", "latitude".
#'
#' @param method a string indicating the R0 ento-epidemiological
#' model used to obtain the R0 outputs. Methods include "caminade",
#' "wesolowski", "liuhelmerssohn" and "mordecai". Alternatively, an
#' extrapolation to real life data can be made by setting "method" to "empirical".
#'
#' @param mm A Kramer probability matrix of dimensions (2, longitude, latitude) to be passed
#' if "caminade" is selected as the computation method. The spatial dimensions must
#' match those of temp and be equal in length. Valid dimension names are "lon", "longitude",
#' "lat", "latitude".
#'
#' @param lon_dim a character vector indicating the longitude dimension name in
#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
#'
#' @param lat_dim a character vector indicating the latitude dimension name in
#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
#'
#'@param ncores An integer indicating the number of cores to use in parallel
#' computation for temporal subsetting.
#'@return An s2dv_cube object containing the R0 Indices (unitless).
#'
#'@examples
#' dims <- c(time = 100, lat = 5, lon = 4)
#' temp <- array(
#' runif(prod(dims), min = 15, max = 35),
#' dim = dims
#' )
#'
#' R0 <- RNoughtIndices(temp, method = "mordecai", mm = NULL, lon_dim = NULL,
#' lat_dim = NULL, ncores = NULL)
#' # Caminade method (requires Kramer matrix)
#' mm <- array(runif(2 * 5 * 4, 0.1, 1.0),
#' dim = c(2, 4, 5))
#' names(dim(mm)) <- c("", "lon", "lat")
#' R0 <- RNoughtIndices(temp, method = "caminade", mm = mm, lon_dim = "lon",
#' lat_dim = "lat", ncores = NULL)
#'
#'@importFrom utils read.table
#'@import multiApply
#'@export
RNoughtIndices <- function(temp, method, mm, lon_dim,
lat_dim, ncores = NULL) {
if (!(method %in% c("caminade", "wesolowski", "liuhelmerssohn", "mordecai",
"empirical"))) {
stop("Invalid value for method. Accepted methods include 'caminade',
'wesolowski', 'liuhelmerssohn' and 'mordecai'. Alternatively, an empirical
calculation can be performed with 'empirical'")
} else if (method == "caminade") {
if (is.null(mm)) {
stop("The method 'caminade' has been selected but no Kramer probability matrix",
" in 'mm' has been provided.")
}
# Ensure the input probability matrix has dimensions (2, lon_dim, lat_dim)
if (length(dim(mm)) != 3 || dim(mm)[1] != 2) {
stop("The probability matrix 'mm' must be a 3D array with dimensions (2, lon_dim, lat_dim).")
}
# Ensure the input temperature array is at least 2D with lon_dim and lat_dim
if (!is.array(temp) || length(dim(temp)) < 2) {
stop("The temperature array must be at least 2D with spatial dimensions (lat_dim, lon_dim).")
}
if (is.null(names(dim(mm))) || is.null(names(dim(temp))) ||
!all(c(lon_dim, lat_dim) %in% names(dim(mm))) ||
!all(c(lon_dim, lat_dim) %in% names(dim(temp)))) {
stop("Both 'mm' and 'temperature' must have the spatial dimensions set in 'lon_dim'",
" and 'lat_dim'.")
}
if (dim(mm)[which(names(dim(mm)) == lat_dim)] != dim(temp)[which(names(dim(temp)) == lat_dim)] ||
dim(mm)[which(names(dim(mm)) == lon_dim)] != dim(temp)[which(names(dim(temp)) == lon_dim)]) {
stop("The arrays 'mm' and 'temperature' don't have matching spatial dimension lengths")
}
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c(lat_dim, lon_dim),
fun = .caminade,
mm = mm,
nlat = dim(mm)[which(names(dim(mm)) == lat_dim)],
nlon = dim(mm)[which(names(dim(mm)) == lon_dim)],
ncores = ncores
)$output1
env_index <- drop(aperm(env_index, c(3:length(dim(env_index)), 1, 2)))
} else if (method == "liuhelmerssohn") {
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c("index"),
fun = .liuhelmerssohn,
ncores = ncores
)$output1
env_index <- drop(array(env_index, dim = dim(env_index)[-1]))
} else if (method == "mordecai") {
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c("index"),
fun = .mordecai,
ncores = ncores
)$output1
env_index <- drop(array(env_index, dim = dim(env_index)[-1]))
} else if (method == "wesolowski") {
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c("index"),
fun = .wesolowski,
ncores = ncores
)$output1
env_index <- drop(array(env_index, dim = dim(env_index)[-1]))
} else if (method == "empirical") {
empirical <- read.table(
system.file(
"index_empirical", "empirical.txt", package = "CSIndicators", mustWork = TRUE
)
)
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c("index"),
fun = .empirical,
empirical = empirical,
ncores = ncores
)$output1
}
return(env_index)
}
.caminade <- function(temp, mm, nlat, nlon) {
# Setting parameter values for Caminade et al., 2015:
a <- array(NA, dim = c(2, nlat, nlon))
a[1, , ] <- 0.0043 * temp + 0.0943 #biting rates per day
a[2, , ] <- 0.5 * (0.0043 * temp + 0.0943)
# Vector preference (0-1); typical: 0.88-1.0
phi <- NULL
phi[1] <- 1
#typical: 0.24-1.0
phi[2] <- 0.5
# Transmission probability (vector to host); typical: 0.1-0.75
b <- NULL
b[1] <- 0.5
b[2] <- 0.5
# Transmission probability (host to vector); typical: 0.1-0.75
beta <- NULL
beta[1] <- 0.1
beta[2] <- 0.033
mu <- array(NA, dim = c(2, nlat, nlon))
for (ix in 1:nlon) {
for (iy in 1:nlat) {
# Including option for NA
if (is.na(temp[iy, ix]) == TRUE) {
mu[1, iy, ix] <- NA
mu[2, iy, ix] <- NA
} else {
if (temp[iy, ix] < 22) {
mu[1, iy, ix] <- 1 / (1.22 + exp(-3.05 + 0.72 * temp[iy, ix])) + 0.196
} else {
mu[1, iy, ix] <- 1 / (1.24 + exp(35.14 - 0.905 * temp[iy, ix])) + 0.195
}
if (temp[iy, ix] < 15) {
mu[2, iy, ix] <- 1 / (1.1 + exp(-4.04 + 0.576 * temp[iy, ix])) + 0.12
} else if (15 <= temp[iy, ix] && temp[iy, ix] < 26.3) {
mu[2, iy, ix] <- 0.000339 * temp[iy, ix]^2 - 0.0189 * temp[iy, ix] + 0.336
} else if (temp[iy, ix] >= 26.3) {
mu[2, iy, ix] <- 1 / (1.065 + exp(32.2 - 0.92 * temp[iy, ix])) + 0.0747
}
}
}
}
# Inverse of the EIP
eip <- array(NA, dim = c(2, nlat, nlon))
eip[1, , ] <- 4 + exp(5.15 - 0.123 * temp)
eip[2, , ] <- 1.03 * (4 + exp(5.15 - 0.123 * temp))
#Recovery rate
r <- 1 / 7
# Calculating the R0 parameters:
env_index <- array(NA, dim = c(2, nlat, nlon))
for (i in 1:2) {
for (ix in 1:nlon) {
for (iy in 1:nlat) {
if (is.na(temp[iy, ix]) == TRUE) {
env_index[i, iy, ix] <- NA
} else {
env_index[i, iy, ix] <- (b[i] * beta[i] * a[i, iy, ix]^2 / mu[i, iy, ix]) *
(1 / eip[i, iy, ix] / (1 / eip[i, iy, ix] + mu[i, iy, ix])) * (phi[i]^2 *
mm[i, ix, iy] / r)
}
}
}
}
env_index <- sqrt(env_index[1, , ] + env_index[2, , ])
dim(env_index) <- dim(temp)
return(env_index)
}
.liuhelmerssohn <- function(temp) {
if (is.na(temp) == TRUE) {
env_index <- NA
} else {
if (temp < 12.4 || temp > 32) {
a <- 0
} else {
a <- 0.0043 * temp + 0.0943
}
if (temp < 12.4 || temp > 32.5) {
b <- 0
} else if (temp >= 12.4 && temp <= 26.1) {
b <- 0.0729 * temp - 0.9037
} else {
b <- 1
}
if (temp < 12.3 || temp > 32.5) {
c <- 0
} else {
c <- 0.001044 * temp * (temp - 12.286) * (32.461 - temp)^(1 / 2)
if (is.nan(c)) {
c <- 0
}
}
bc <- b * c
if (temp < 14 || temp > 32) {
mu <- 0.04168548
} else {
mu <- 0.8692 - 0.1590 * temp + 0.01116 * temp^2 - 3.408e-04 * temp^3 + 3.809e-06 * temp^4
}
if (temp < 12 || temp > 36) {
eip <- 0
} else {
eip <- 4 + exp(-0.123 * temp + 5.15)
}
env_index <- (a^2 * bc * exp(-mu * eip) / mu)
}
dim(env_index) <- dim(temp)
return(env_index)
}
.mordecai <- function(temp) {
ec <- 0.00001
rr <- 8.3144598
if (is.na(temp) == TRUE) {
env_index <- NA
} else {
a <- -5.4 + 1.8 * temp - 0.2124 * temp^2 + 0.01015 * temp^3 - 1.515e-04 * temp^4
if (a < 0) {
a <- 0
}
b <- 0.0729 * temp - 0.97
if (b < 0) {
b <- 0
}
c <- 1.044e-03 * temp * (temp - 12.286) * (32.461 - temp)^(1 / 2)
if (is.nan(c)) {
c <- 0
}
bc <- b * c
x1 <- (0.8692 - 0.1599 * temp + 0.01116 * temp^2 - 3.408e-04 * temp^3 + 3.809e-06 * temp^4)
if (x1 < 0.01) {
x1 <- 0.01
}
lf <- 1 / x1
mu <- 1 / (lf + ec)
e2a <- exp(-(2.13 - 0.3787 * temp + 0.02457 * temp^2 - 6.778e-04 * temp^3 + 6.794e-06 * temp^4))
mdr <- 0.131 - 0.05723 * temp + 0.01164 * temp^2 - 0.001341 * temp^3 +
0.00008723 * temp^4 - 3.017e-06 * temp^5 + 5.153e-08 * temp^6 - 3.42e-10 * temp^7
pdr <- (1 + exp((6.203e21) / rr * (1 / (-2.176e30) - 1 / (temp + 273.15)))) /
((3.3589e-03 * (temp + 273.15)) / 298 * exp((1.500 / rr) * (1 / 298 - 1 / (temp + 273.15))))
if (temp < 8.02 || temp > 35.65) {
efd <- 0
} else {
efd <- 4.88E-02 * (temp - 8.02) * (35.65 - temp)^0.5
}
env_index <- sqrt((a^2 * bc * (efd * e2a * mdr / (mu)^2) *
exp((-mu / (pdr + ec)))) / (mu)) / 1000
}
dim(env_index) <- dim(temp)
return(env_index)
}
.wesolowski <- function(temp) {
# Setting parameter values for Wesolowski et al., 2015:
ec <- 0.00001
rr <- 8.3144598
if (is.na(temp) == TRUE) {
env_index <- NA
} else {
a <- -5.4 + 1.8 * temp - 0.2124 * temp^2 + 0.01015 * temp^3 - 1.515e-04 * temp^4
if (a < 0) {
a <- 0
}
b <- 0.0729 * temp - 0.97
if (b < 0) {
b <- 0
}
c <- 1.044e-03 * temp * (temp - 12.286) * (32.461 - temp)^(1 / 2)
if (is.nan(c)) {
c <- 0
}
bc <- b * c
e2a <- exp(-(2.13 - 0.3787 * temp + 0.02457 * temp^2 - 6.778e-04 * temp^3 + 6.794e-06 * temp^4))
x1 <- (0.8692 - 0.1599 * temp + 0.01116 * temp^2 - 3.408e-04 * temp^3 + 3.809e-06 * temp^4)
if (x1 < 0.01) {
x1 <- 0.01
}
lf <- 1 / x1
mdr <- 0.131 - 0.05723 * temp + 0.01164 * temp^2 - 0.001341 * temp^3 +
0.00008723 * temp^4 - 3.017e-06 * temp^5 + 5.153e-08 * temp^6 - 3.42e-10 * temp^7
pdr <- (1 + exp((6.203e21) / rr * (1 / (-2.176e30) - 1 / (temp + 273.15)))) /
((3.3589e-03 * (temp + 273.15)) / 298 * exp((1.500 / rr) * (1 / 298 - 1 / (temp + 273.15))))
mu <- 1 / (lf + ec)
env_index <- sqrt((a^2 * bc * (e2a * mdr / (mu)^2) * exp((-mu / (pdr + ec)))) / (mu)) / 1000
}
dim(env_index) <- dim(temp)
return(env_index)
}
.empirical <- function(temp, empirical) {
env_index <- ifelse(
is.na(temp) == TRUE || temp < min(empirical$temp) || temp > max(empirical$temp),
NA,
sqrt(empirical$aegy[which.min(abs(empirical$temp - temp)) / 10] +
empirical$albo[which.min(abs(empirical$temp - temp))]
)
)
return(env_index)
}
Try the CSIndicators package in your browser
Any scripts or data that you put into this service are public.
CSIndicators documentation built on Nov. 25, 2025, 1:07 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 .empirical .wesolowski .mordecai .liuhelmerssohn .caminade RNoughtIndices CST_RNoughtIndices
Documented in CST_RNoughtIndices RNoughtIndices
#'R0 Index computation on s2dv_cube objects
#'
#'@author Javier Corvillo, \email{javier.corvillo@bsc.es}
#'@description This function computes the environmental contribution to Aedes-borne
#' disease transmissibility. Values higher than 1 imply that the environmental conditions allow for
#' the disease to be transmitted to more than 1 person (the infection can spread) while R0 values lower
#' than 1 imply that the environmental conditions are not suitable for the disease to spread.
#'
#' This function utilizes four possible ento-epidemiological models, each of them stated
#' in Caminade et al., 2015; Liu-Helmerssohn et al., 2014; and Mordecai et
#' al., 2017 and Wesolowski et al., 2015 respectively. Additionally, an adjustment to
#' real life data, using recorded DENV data recorded between 2014-2017 in the Caribbean, can also be performed.
#'
#'@references Caminade, Cyril et al. (Jan. 2017). ‘Global risk model for vector-borne
#' transmission of Zika virus reveals the role of El Niño 2015’.
#' In: Proceedings of the National Academy of Sciences 114.1.
#' Publisher: Proceedings of the National Academy of
#' Sciences, pp. 119–124. doi: 10.1073/pnas.1614303114.
#' url: https://www.pnas.org/doi/full/10.1073/pnas.1614303114.
#'
#'@references Liu-Helmerssohn, Jing et al. (Mar. 2014). ‘Vectorial Capacity of Aedes
#' aegypti: Effects of Temperature and Implications for Global Dengue Epidemic Potential’. en.
#' In: PLOS ONE 9.3.
#' Publisher: Public Library of Science, e89783.
#' issn: 1932-6203. doi: 10.1371/journal.pone.0089783.
#' url: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0089783.
#'
#' @references Mordecai, Erin A. et al. (Apr. 2017). ‘Detecting the impact of temperature on
#' transmission of Zika, dengue, and chikungunya using mechanistic models’. en.
#' In: PLOS Neglected Tropical Diseases 11.4.
#' Publisher: Public Library of Science, e0005568.
#' issn: 1935-2735. doi: 10.1371/journal.pntd.0005568.
#' url: https://journals.plos.org/plosntds/article?id=10.1371/journal.pntd.0005568.
#'
#' @references Wesolowski, Amy et al. (Sept. 2015). ‘Impact of human mobility on the emergence
#' of dengue epidemics in Pakistan’.
#' In: Proceedings of the National Academy of Sciences 112.38.
#' Publisher: Proceedings of the National Academy of Sciences, pp. 11887–11892.
#' doi: 10.1073/pnas.1504964112.
#' url: https://www.pnas.org/doi/full/10.1073/pnas.1504964112.
#'
#'@param temp An s2dv_cube object with temperature values, expressed in degrees Celsius. If
#' "caminade" is selected in "method", named spatial dimensions are required in the s2dv_cube object,
#' ordered as ("latitude", "longitude"). Valid dimension names are "lon", "longitude", "lat", "latitude".
#'
#' @param method a string indicating the R0 ento-epidemiological
#' model used to obtain the R0 outputs. Methods include "caminade",
#' "wesolowski", "liuhelmerssohn" and "mordecai". Alternatively, an
#' extrapolation to real life data can be made by setting "method" to "empirical".
#'
#' @param mm A Kramer probability matrix of dimensions (2, longitude, latitude) to be passed
#' if "caminade" is selected as the computation method. The spatial dimensions must
#' match those of temp and be equal in length. Valid dimension names are "lon", "longitude",
#' "lat", "latitude".
#'
#' @param lon_dim a character vector indicating the longitude dimension name in
#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
#'
#' @param lat_dim a character vector indicating the latitude dimension name in
#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
#'
#'@param ncores An integer indicating the number of cores to use in parallel
#' computation for temporal subsetting.
#'@return An s2dv_cube object containing the R0 Indices (unitless).
#'
#'@examples
#'dims <- c(time = 100, lat = 5, lon = 4)
#'temp_cube <- list(
#'data = array(
#' runif(prod(dims), min = 15, max = 35),
#' dim = dims
#' ),
#' attrs = list(
#' Variable = list(varName = '2m Temperature')
#' )
#' )
#' class(temp_cube) <- 's2dv_cube'
#'
#' R0 <- CST_RNoughtIndices(temp_cube, method = "mordecai", mm = NULL,
#' lon_dim = NULL, lat_dim = NULL, ncores = NULL)
#'
#' # Caminade method (requires Kramer matrix)
#' mm <- array(runif(2 * 5 * 4, 0.1, 1.0),
#' dim = c(2, 4, 5))
#' names(dim(mm)) <- c("", "lon", "lat")
#' R0 <- CST_RNoughtIndices(temp_cube, method = "caminade", mm = mm,
#' lon_dim = "lon", lat_dim = "lat", ncores = NULL)
#'
#'@import multiApply
#'@export
CST_RNoughtIndices <- function(temp, method, mm, lon_dim, lat_dim, ncores = NULL) {
# Check "s2dv_cube"
if (!inherits(temp, "s2dv_cube")) {
stop("Parameter 'temp' must be of the class 's2dv_cube'.")
}
RNought <- RNoughtIndices(temp = temp$data, method = method, mm = mm,
lon_dim = lon_dim, lat_dim = lat_dim,
ncores = ncores)
temp$data <- RNought
if ("Variable" %in% names(temp$attrs)) {
if ("varName" %in% names(temp$attrs$Variable)) {
temp$attrs$Variable$varName <- "RNoughtIndex"
}
}
return(temp)
}
#'R0 Index computation on s2dv_cube objects
#'
#'@author Javier Corvillo, \email{javier.corvillo@bsc.es}
#'@description This function computes the environmental contribution to Aedes-borne
#' disease transmissibility. Values higher than 1 imply that the environmental conditions allow for
#' the disease to be transmitted to more than 1 person (the infection can spread) while R0 values lower
#' than 1 imply that the environmental conditions are not suitable for the disease to spread.
#'
#' This function utilizes four possible ento-epidemiological models, each of them stated
#' in Caminade et al., 2015; Liu-Helmerssohn et al., 2014; and Mordecai et
#' al., 2017 and Wesolowski et al., 2015 respectively. Additionally, an adjustment to
#' real life data, using recorded DENV data recorded between 2014-2017 in the Caribbean, can also be performed.
#'
#'@references Caminade, Cyril et al. (Jan. 2017). ‘Global risk model for vector-borne
#' transmission of Zika virus reveals the role of El Niño 2015’.
#' In: Proceedings of the National Academy of Sciences 114.1.
#' Publisher: Proceedings of the National Academy of
#' Sciences, pp. 119–124. doi: 10.1073/pnas.1614303114.
#' url: https://www.pnas.org/doi/full/10.1073/pnas.1614303114.
#'
#'@references Liu-Helmerssohn, Jing et al. (Mar. 2014). ‘Vectorial Capacity of Aedes
#' aegypti: Effects of Temperature and Implications for Global Dengue Epidemic Potential’. en.
#' In: PLOS ONE 9.3.
#' Publisher: Public Library of Science, e89783.
#' issn: 1932-6203. doi: 10.1371/journal.pone.0089783.
#' url: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0089783.
#'
#' @references Mordecai, Erin A. et al. (Apr. 2017). ‘Detecting the impact of temperature on
#' transmission of Zika, dengue, and chikungunya using mechanistic models’. en.
#' In: PLOS Neglected Tropical Diseases 11.4.
#' Publisher: Public Library of Science, e0005568.
#' issn: 1935-2735. doi: 10.1371/journal.pntd.0005568.
#' url: https://journals.plos.org/plosntds/article?id=10.1371/journal.pntd.0005568.
#'
#' @references Wesolowski, Amy et al. (Sept. 2015). ‘Impact of human mobility on the emergence
#' of dengue epidemics in Pakistan’.
#' In: Proceedings of the National Academy of Sciences 112.38.
#' Publisher: Proceedings of the National Academy of Sciences, pp. 11887–11892.
#' doi: 10.1073/pnas.1504964112.
#' url: https://www.pnas.org/doi/full/10.1073/pnas.1504964112.
#'
#'@param temp An array with temperature values, expressed in degrees Celsius. If
#' "caminade" is selected in "method", named spatial dimensions are required in the s2dv_cube object,
#' ordered as ("latitude", "longitude"). Valid dimension names are "lon", "longitude", "lat", "latitude".
#'
#' @param method a string indicating the R0 ento-epidemiological
#' model used to obtain the R0 outputs. Methods include "caminade",
#' "wesolowski", "liuhelmerssohn" and "mordecai". Alternatively, an
#' extrapolation to real life data can be made by setting "method" to "empirical".
#'
#' @param mm A Kramer probability matrix of dimensions (2, longitude, latitude) to be passed
#' if "caminade" is selected as the computation method. The spatial dimensions must
#' match those of temp and be equal in length. Valid dimension names are "lon", "longitude",
#' "lat", "latitude".
#'
#' @param lon_dim a character vector indicating the longitude dimension name in
#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
#'
#' @param lat_dim a character vector indicating the latitude dimension name in
#' the element 'temp' if "caminade" is chosen in "method". Otherwise, can be set to NULL
#'
#'@param ncores An integer indicating the number of cores to use in parallel
#' computation for temporal subsetting.
#'@return An s2dv_cube object containing the R0 Indices (unitless).
#'
#'@examples
#' dims <- c(time = 100, lat = 5, lon = 4)
#' temp <- array(
#' runif(prod(dims), min = 15, max = 35),
#' dim = dims
#' )
#'
#' R0 <- RNoughtIndices(temp, method = "mordecai", mm = NULL, lon_dim = NULL,
#' lat_dim = NULL, ncores = NULL)
#' # Caminade method (requires Kramer matrix)
#' mm <- array(runif(2 * 5 * 4, 0.1, 1.0),
#' dim = c(2, 4, 5))
#' names(dim(mm)) <- c("", "lon", "lat")
#' R0 <- RNoughtIndices(temp, method = "caminade", mm = mm, lon_dim = "lon",
#' lat_dim = "lat", ncores = NULL)
#'
#'@importFrom utils read.table
#'@import multiApply
#'@export
RNoughtIndices <- function(temp, method, mm, lon_dim,
lat_dim, ncores = NULL) {
if (!(method %in% c("caminade", "wesolowski", "liuhelmerssohn", "mordecai",
"empirical"))) {
stop("Invalid value for method. Accepted methods include 'caminade',
'wesolowski', 'liuhelmerssohn' and 'mordecai'. Alternatively, an empirical
calculation can be performed with 'empirical'")
} else if (method == "caminade") {
if (is.null(mm)) {
stop("The method 'caminade' has been selected but no Kramer probability matrix",
" in 'mm' has been provided.")
}
# Ensure the input probability matrix has dimensions (2, lon_dim, lat_dim)
if (length(dim(mm)) != 3 || dim(mm)[1] != 2) {
stop("The probability matrix 'mm' must be a 3D array with dimensions (2, lon_dim, lat_dim).")
}
# Ensure the input temperature array is at least 2D with lon_dim and lat_dim
if (!is.array(temp) || length(dim(temp)) < 2) {
stop("The temperature array must be at least 2D with spatial dimensions (lat_dim, lon_dim).")
}
if (is.null(names(dim(mm))) || is.null(names(dim(temp))) ||
!all(c(lon_dim, lat_dim) %in% names(dim(mm))) ||
!all(c(lon_dim, lat_dim) %in% names(dim(temp)))) {
stop("Both 'mm' and 'temperature' must have the spatial dimensions set in 'lon_dim'",
" and 'lat_dim'.")
}
if (dim(mm)[which(names(dim(mm)) == lat_dim)] != dim(temp)[which(names(dim(temp)) == lat_dim)] ||
dim(mm)[which(names(dim(mm)) == lon_dim)] != dim(temp)[which(names(dim(temp)) == lon_dim)]) {
stop("The arrays 'mm' and 'temperature' don't have matching spatial dimension lengths")
}
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c(lat_dim, lon_dim),
fun = .caminade,
mm = mm,
nlat = dim(mm)[which(names(dim(mm)) == lat_dim)],
nlon = dim(mm)[which(names(dim(mm)) == lon_dim)],
ncores = ncores
)$output1
env_index <- drop(aperm(env_index, c(3:length(dim(env_index)), 1, 2)))
} else if (method == "liuhelmerssohn") {
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c("index"),
fun = .liuhelmerssohn,
ncores = ncores
)$output1
env_index <- drop(array(env_index, dim = dim(env_index)[-1]))
} else if (method == "mordecai") {
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c("index"),
fun = .mordecai,
ncores = ncores
)$output1
env_index <- drop(array(env_index, dim = dim(env_index)[-1]))
} else if (method == "wesolowski") {
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c("index"),
fun = .wesolowski,
ncores = ncores
)$output1
env_index <- drop(array(env_index, dim = dim(env_index)[-1]))
} else if (method == "empirical") {
empirical <- read.table(
system.file(
"index_empirical", "empirical.txt", package = "CSIndicators", mustWork = TRUE
)
)
env_index <- Apply(
data = array(temp, c(index = 1, dim(temp))),
target_dims = c("index"),
fun = .empirical,
empirical = empirical,
ncores = ncores
)$output1
}
return(env_index)
}
.caminade <- function(temp, mm, nlat, nlon) {
# Setting parameter values for Caminade et al., 2015:
a <- array(NA, dim = c(2, nlat, nlon))
a[1, , ] <- 0.0043 * temp + 0.0943 #biting rates per day
a[2, , ] <- 0.5 * (0.0043 * temp + 0.0943)
# Vector preference (0-1); typical: 0.88-1.0
phi <- NULL
phi[1] <- 1
#typical: 0.24-1.0
phi[2] <- 0.5
# Transmission probability (vector to host); typical: 0.1-0.75
b <- NULL
b[1] <- 0.5
b[2] <- 0.5
# Transmission probability (host to vector); typical: 0.1-0.75
beta <- NULL
beta[1] <- 0.1
beta[2] <- 0.033
mu <- array(NA, dim = c(2, nlat, nlon))
for (ix in 1:nlon) {
for (iy in 1:nlat) {
# Including option for NA
if (is.na(temp[iy, ix]) == TRUE) {
mu[1, iy, ix] <- NA
mu[2, iy, ix] <- NA
} else {
if (temp[iy, ix] < 22) {
mu[1, iy, ix] <- 1 / (1.22 + exp(-3.05 + 0.72 * temp[iy, ix])) + 0.196
} else {
mu[1, iy, ix] <- 1 / (1.24 + exp(35.14 - 0.905 * temp[iy, ix])) + 0.195
}
if (temp[iy, ix] < 15) {
mu[2, iy, ix] <- 1 / (1.1 + exp(-4.04 + 0.576 * temp[iy, ix])) + 0.12
} else if (15 <= temp[iy, ix] && temp[iy, ix] < 26.3) {
mu[2, iy, ix] <- 0.000339 * temp[iy, ix]^2 - 0.0189 * temp[iy, ix] + 0.336
} else if (temp[iy, ix] >= 26.3) {
mu[2, iy, ix] <- 1 / (1.065 + exp(32.2 - 0.92 * temp[iy, ix])) + 0.0747
}
}
}
}
# Inverse of the EIP
eip <- array(NA, dim = c(2, nlat, nlon))
eip[1, , ] <- 4 + exp(5.15 - 0.123 * temp)
eip[2, , ] <- 1.03 * (4 + exp(5.15 - 0.123 * temp))
#Recovery rate
r <- 1 / 7
# Calculating the R0 parameters:
env_index <- array(NA, dim = c(2, nlat, nlon))
for (i in 1:2) {
for (ix in 1:nlon) {
for (iy in 1:nlat) {
if (is.na(temp[iy, ix]) == TRUE) {
env_index[i, iy, ix] <- NA
} else {
env_index[i, iy, ix] <- (b[i] * beta[i] * a[i, iy, ix]^2 / mu[i, iy, ix]) *
(1 / eip[i, iy, ix] / (1 / eip[i, iy, ix] + mu[i, iy, ix])) * (phi[i]^2 *
mm[i, ix, iy] / r)
}
}
}
}
env_index <- sqrt(env_index[1, , ] + env_index[2, , ])
dim(env_index) <- dim(temp)
return(env_index)
}
.liuhelmerssohn <- function(temp) {
if (is.na(temp) == TRUE) {
env_index <- NA
} else {
if (temp < 12.4 || temp > 32) {
a <- 0
} else {
a <- 0.0043 * temp + 0.0943
}
if (temp < 12.4 || temp > 32.5) {
b <- 0
} else if (temp >= 12.4 && temp <= 26.1) {
b <- 0.0729 * temp - 0.9037
} else {
b <- 1
}
if (temp < 12.3 || temp > 32.5) {
c <- 0
} else {
c <- 0.001044 * temp * (temp - 12.286) * (32.461 - temp)^(1 / 2)
if (is.nan(c)) {
c <- 0
}
}
bc <- b * c
if (temp < 14 || temp > 32) {
mu <- 0.04168548
} else {
mu <- 0.8692 - 0.1590 * temp + 0.01116 * temp^2 - 3.408e-04 * temp^3 + 3.809e-06 * temp^4
}
if (temp < 12 || temp > 36) {
eip <- 0
} else {
eip <- 4 + exp(-0.123 * temp + 5.15)
}
env_index <- (a^2 * bc * exp(-mu * eip) / mu)
}
dim(env_index) <- dim(temp)
return(env_index)
}
.mordecai <- function(temp) {
ec <- 0.00001
rr <- 8.3144598
if (is.na(temp) == TRUE) {
env_index <- NA
} else {
a <- -5.4 + 1.8 * temp - 0.2124 * temp^2 + 0.01015 * temp^3 - 1.515e-04 * temp^4
if (a < 0) {
a <- 0
}
b <- 0.0729 * temp - 0.97
if (b < 0) {
b <- 0
}
c <- 1.044e-03 * temp * (temp - 12.286) * (32.461 - temp)^(1 / 2)
if (is.nan(c)) {
c <- 0
}
bc <- b * c
x1 <- (0.8692 - 0.1599 * temp + 0.01116 * temp^2 - 3.408e-04 * temp^3 + 3.809e-06 * temp^4)
if (x1 < 0.01) {
x1 <- 0.01
}
lf <- 1 / x1
mu <- 1 / (lf + ec)
e2a <- exp(-(2.13 - 0.3787 * temp + 0.02457 * temp^2 - 6.778e-04 * temp^3 + 6.794e-06 * temp^4))
mdr <- 0.131 - 0.05723 * temp + 0.01164 * temp^2 - 0.001341 * temp^3 +
0.00008723 * temp^4 - 3.017e-06 * temp^5 + 5.153e-08 * temp^6 - 3.42e-10 * temp^7
pdr <- (1 + exp((6.203e21) / rr * (1 / (-2.176e30) - 1 / (temp + 273.15)))) /
((3.3589e-03 * (temp + 273.15)) / 298 * exp((1.500 / rr) * (1 / 298 - 1 / (temp + 273.15))))
if (temp < 8.02 || temp > 35.65) {
efd <- 0
} else {
efd <- 4.88E-02 * (temp - 8.02) * (35.65 - temp)^0.5
}
env_index <- sqrt((a^2 * bc * (efd * e2a * mdr / (mu)^2) *
exp((-mu / (pdr + ec)))) / (mu)) / 1000
}
dim(env_index) <- dim(temp)
return(env_index)
}
.wesolowski <- function(temp) {
# Setting parameter values for Wesolowski et al., 2015:
ec <- 0.00001
rr <- 8.3144598
if (is.na(temp) == TRUE) {
env_index <- NA
} else {
a <- -5.4 + 1.8 * temp - 0.2124 * temp^2 + 0.01015 * temp^3 - 1.515e-04 * temp^4
if (a < 0) {
a <- 0
}
b <- 0.0729 * temp - 0.97
if (b < 0) {
b <- 0
}
c <- 1.044e-03 * temp * (temp - 12.286) * (32.461 - temp)^(1 / 2)
if (is.nan(c)) {
c <- 0
}
bc <- b * c
e2a <- exp(-(2.13 - 0.3787 * temp + 0.02457 * temp^2 - 6.778e-04 * temp^3 + 6.794e-06 * temp^4))
x1 <- (0.8692 - 0.1599 * temp + 0.01116 * temp^2 - 3.408e-04 * temp^3 + 3.809e-06 * temp^4)
if (x1 < 0.01) {
x1 <- 0.01
}
lf <- 1 / x1
mdr <- 0.131 - 0.05723 * temp + 0.01164 * temp^2 - 0.001341 * temp^3 +
0.00008723 * temp^4 - 3.017e-06 * temp^5 + 5.153e-08 * temp^6 - 3.42e-10 * temp^7
pdr <- (1 + exp((6.203e21) / rr * (1 / (-2.176e30) - 1 / (temp + 273.15)))) /
((3.3589e-03 * (temp + 273.15)) / 298 * exp((1.500 / rr) * (1 / 298 - 1 / (temp + 273.15))))
mu <- 1 / (lf + ec)
env_index <- sqrt((a^2 * bc * (e2a * mdr / (mu)^2) * exp((-mu / (pdr + ec)))) / (mu)) / 1000
}
dim(env_index) <- dim(temp)
return(env_index)
}
.empirical <- function(temp, empirical) {
env_index <- ifelse(
is.na(temp) == TRUE || temp < min(empirical$temp) || temp > max(empirical$temp),
NA,
sqrt(empirical$aegy[which.min(abs(empirical$temp - temp)) / 10] +
empirical$albo[which.min(abs(empirical$temp - temp))]
)
)
return(env_index)
}
Try the CSIndicators package in your browser
Any scripts or data that you put into this service are public.
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.