R/incubsim.R
In oswaldogressani/EpiLPS: A Fast and Flexible Bayesian Tool for Estimating Epidemiological Parameters
Defines functions
incubsim
Documented in
incubsim
#' Simulation of incubation times
#'
#' @description
#' This routine simulates symptom onset times, generation times and left and
#' right bounds of infecting exposure windows for source-recipient pairs based
#' on a given incubation distribution and a generation time distribution. The
#' generation time distribution is assumed to be a Weibull with shape 2.826 and
#' scale 5.665 (Ferretti et al. 2020). The incubation distribution can be chosen
#' by the user among:
#'
#' \itemize{
#' \item{A LogNormal distribution with a mean of 5.5 days and standard deviation of
#' 2.1 days (Ferretti et al. 2020).}
#' \item{A Weibull distribution (shape-scale parameterization) with a mean of 6.4 days and standard deviation of
#' 2.3 days (Backer et al. 2020).}
#' \item{An artificial bimodal distribution constructed from a mixture of two
#' Weibull distributions.}
#' \item{A Gamma distribution (shape-rate parameterization) with a mean of 3.8 days and standard deviation of
#' 2.9 days (Donnelly et al. 2003).}
#' }
#'
#' @usage incubsim(incubdist = c("LogNormal","Weibull","MixWeibull", "Gamma"), coarseness = 1,
#' n = 100, tmax = 20, tgridlen = 500, plotsim = FALSE)
#'
#' @param incubdist The distribution of the incubation period.
#' @param coarseness The average coarseness of the data. Default is 1 day.
#' @param n The sample size.
#' @param tmax The upper bound on which to evaluate the \code{incubdist} density.
#' @param plotsim Graphical visualization of the simulated data?
#' @param tgridlen The number of grid points on which to evaluate the density.
#'
#' @return A list including (among others) the left and right bounds of the
#' incubation period.
#'
#' @author Oswaldo Gressani \email{oswaldo_gressani@hotmail.fr}
#'
#' @references Ferretti, L., Wymant, C., Kendall, et al. (2020). Quantifying
#' SARS-CoV-2 transmission suggests epidemic control with digital contact
#' tracing. \emph{Science}, \strong{368}(6491), eabb6936.
#' @references Backer, J. A., Klinkenberg, D., & Wallinga, J. (2020). Incubation
#' period of 2019 novel coronavirus (2019-nCoV) infections among
#' travellers from Wuhan, China, 20–28 January 2020. \emph{Eurosurveillance},
#' \strong{25}(5), 2000062.
#' @references Donnelly, C. A., Ghani, A. C., Leung, G. M., et al. (2003).
#' Epidemiological determinants of spread of causal agent of severe acute
#' respiratory syndrome in Hong Kong. \emph{The Lancet}, \strong{361}(9371),
#' 1761-1766.
#'
#' @examples
#' simdat <- incubsim(n = 50) # Simulation of incubation times for 50 cases.
#' simdat$Dobsincub # Left and right bounds of incubation period.
#'
#' @export
incubsim <- function(incubdist = c("LogNormal","Weibull","MixWeibull", "Gamma"),
coarseness = 1, n = 100, tmax = 20, tgridlen = 500, plotsim = FALSE){
mixWeibull <- function(n, alpha, Ishape1 = 2.356, Iscale1 = 3.564,
Ishape2 = 9.351, Iscale2 = 12.563) {
#Probability weights in the mixture
w1 <- alpha
w2 <- 1 - alpha
#Sample scheme with replacement and weights on the elements being sampled
w <- sample(c(1, 2), n, replace = TRUE, prob = c(w1, w2))
#Sample of size n from the mixture of 2 Weibulls
sample <- (w == 1) * stats::rweibull(n = n, shape = Ishape1, scale = Iscale1) +
(w == 2) * stats::rweibull(n = n, shape = Ishape2, scale = Iscale2)
# set.seed(123)
# mixsample <- mixWeibull(n = 50000000, alpha = 0.5)
# summary(mixsample$sample)
# sd(mixsample$sample)
# quantile(mixsample$sample, probs = c(0.05,0.25,0.5,0.75,0.95))
outlist <- list(sample = sample, Ishape1 = Ishape1, Iscale1 = Iscale1,
Ishape2 = Ishape2, Iscale2 = Iscale2)
return(outlist)
}
# Domain of Generation time/Incubation time density
tt <- seq(0, tmax, length = tgridlen)
# Generation interval Weibull parameters #Ferretti et al. 2020
GIshape <- 2.826
# GIshape <- 3.4395
GIscale <- 5.665
if (match.arg(incubdist) == "LogNormal") { # Ferretti et al. 2020
Incubmean <- 1.644
Incubsd <- 0.363
if (coarseness == 1) {
dirforce <- 0.185
} else if (coarseness == 2) {
dirforce <- 0.375
}
Ipdf <- stats::dlnorm(tt, meanlog = Incubmean, sdlog = Incubsd)
meanI <- exp(Incubmean + 0.5 * (Incubsd ^ 2))
medianI <- exp(Incubmean)
sdI <- sqrt(exp(2 * Incubmean + Incubsd ^ 2) * (exp(Incubsd ^ 2) - 1))
Iparams <- c(Incubmean, Incubsd)
} else if (match.arg(incubdist) == "Weibull") { # Backer et al. 2020
Ishape <- 3.00060
Iscale <- 7.170345
if (coarseness == 1) {
dirforce <- 0.160
} else if (coarseness == 2) {
dirforce <- 0.320
}
Ipdf <- stats::dweibull(tt, shape = Ishape, scale = Iscale)
meanI <- Iscale * gamma(1 + 1 / Ishape)
medianI <- stats::qweibull(p = 0.5, shape = Ishape, scale = Iscale)
sdI <- sqrt((Iscale ^ 2) * ((gamma(1 + 2 / Ishape) - (gamma(1 + 1 / Ishape) ^ 2))))
Iparams <- c(Ishape, Iscale)
} else if (match.arg(incubdist) == "MixWeibull"){ # Artificial
if (coarseness == 1) {
dirforce <- 0.132
} else if (coarseness == 2) {
dirforce <- 0.264
}
mmix <- mixWeibull(n=1, alpha = 0.5)
Ishape1 <- mmix$Ishape1
Iscale1 <- mmix$Iscale1
Ishape2 <- mmix$Ishape2
Iscale2 <- mmix$Iscale2
Ipdf <- 0.5 * stats::dweibull(tt, shape = Ishape1, scale = Iscale1) +
0.5 * stats::dweibull(tt, shape = Ishape2, scale = Iscale2)
meanmix <- 0.5 * (Iscale1 * gamma(1 + 1 / Ishape1) +
Iscale2 * gamma(1 + 1 / Ishape2))
sdmix <- 4.622
qmix <- c(1.371, 1.885,2.301,2.680,3.050,3.434,3.856,4.361,5.075,7.191,
9.876,10.701,11.251,11.692,12.080,12.446,12.814,13.218,13.734)
} else if(match.arg(incubdist) == "Gamma") { # Donnelly 2003
Ishape <- 1.739646
Irate <- 0.4566
if (coarseness == 1) {
dirforce <- 0.262
} else if (coarseness == 2) {
dirforce <- 0.558
}
Ipdf <- stats::dgamma(tt, shape = Ishape, rate = Irate)
meanI <- Ishape / Irate
medianI <- stats::qgamma(p = 0.5, shape = Ishape, rate = Irate)
sdI <- sqrt(Ishape / (Irate ^ 2))
Iparams <- c(Ishape, Irate)
}
#-- Generation of symptom onset data and infection exposure intervals
GI <- c() # Vector to host generation times
SI <- c() # Vector to host serial interval times
TI1 <- c() # Vector to host infection times of source/infector
TI2 <- c() # Vector to host infection times of recipient/infectee
I1 <- c() # Vector to host incubation period for source/infector
I2 <- c() # Vector to host incubation period for recipient/infectee
SO1 <- c() # Vector to host symptom onset times of source/infector
SO2 <- c() # Vector to host symptom onset times of recipient/infectee
EL1 <- c() # Vector to host lower bound of exposure time of source/infector
ER1 <- c() # Vector to host upper bound of exposure time of source/infector
EL2 <- c() # Vector to host lower bound of exposure time of recipient/infectee
ER2 <- c() # Vector to host upper bound of exposure time of recipient/infectee
i <- 1
while(i <= n){
# Step 1: Draw generation time from its pdf (Weibull)
GI[i] <- stats::rweibull(n = 1, shape = GIshape, scale = GIscale)
# Step 2: Set time of infection of infector/source (0,365)
TI1[i] <- stats::runif(1, min = 0, max = 365)
# Step 3: Compute the time of infection of infectee/recipient
TI2[i] <- TI1[i] + GI[i]
# Step 4: Generate a pair of incubation periods from pdf (LogNormal)
if (match.arg(incubdist) == "LogNormal") {
I1[i] <- stats::rlnorm(1, meanlog = Incubmean, sdlog = Incubsd)
I2[i] <- stats::rlnorm(1, meanlog = Incubmean, sdlog = Incubsd)
} else if (match.arg(incubdist) == "Weibull") {
I1[i] <- stats::rweibull(1, shape = Ishape, scale = Iscale)
I2[i] <- stats::rweibull(1, shape = Ishape, scale = Iscale)
} else if (match.arg(incubdist) == "MixWeibull") {
I1[i] <- mixWeibull(n = 1, alpha = 0.5)$sample
I2[i] <- mixWeibull(n = 1, alpha = 0.5)$sample
} else if (match.arg(incubdist) == "Gamma") {
I1[i] <- stats::rgamma(1, shape = Ishape, rate = Irate)
I2[i] <- stats::rgamma(1, shape = Ishape, rate = Irate)
}
# Step 5: Compute times of symptom onset
SO1[i] <- TI1[i] + I1[i]
SO2[i] <- TI2[i] + I2[i]
SI[i] <- SO2[i] - SO1[i]
# Step 6: Compute lower/upper bounds of exposure times
minexpo1 <- TI1[i] - dirforce * I1[i]
minexpo2 <- TI2[i] - dirforce * I2[i]
exposure_draw <- function(tlow,tup){
norm_mean <- 0.5 * (tlow + tup)
sdnew <- 0
probtar <- 0.999
ptol <- stats::pnorm(q = tup, mean = norm_mean, sd = sdnew)
while(ptol >= probtar){
sdnew <- sdnew + 0.001
ptol <- stats::pnorm(q = tup, mean = norm_mean, sd = sdnew)
}
tdraw <- 1e+10
while(tdraw >= tup || tdraw <= tlow) {
tdraw <- stats::rnorm(n = 1, mean = norm_mean, sd = sdnew)
}
return(tdraw)
}
# Draw left bound of exposure for infector/infectee
EL1_temp <- exposure_draw(tlow = minexpo1, tup = TI1[i])
EL2_temp <- exposure_draw(tlow = minexpo2, tup = TI2[i])
while(EL2_temp < EL1_temp){# Ferretti et al. 2020 constraints
EL1_temp <- exposure_draw(tlow = minexpo1, tup = TI1[i])
EL2_temp <- exposure_draw(tlow = minexpo2, tup = TI2[i])
}
EL1[i] <- EL1_temp
EL2[i] <- EL2_temp
maxexpo1 <- TI1[i] + dirforce * I1[i]
maxexpo2 <- TI2[i] + dirforce * I2[i]
# Draw right bound of exposure for infector/infectee
ER1_temp <- exposure_draw(tlow = TI1[i], tup = maxexpo1)
ER2_temp <- exposure_draw(tlow = TI2[i], tup = maxexpo2)
while(ER1_temp > min(SO1[i],ER2_temp)){# Ferretti et al. 2020 constraints
ER1_temp <- exposure_draw(tlow = TI1[i], tup = maxexpo1)
ER2_temp <- exposure_draw(tlow = TI2[i], tup = maxexpo2)
}
ER1[i] <- ER1_temp
ER2[i] <- ER2_temp
ediff1 <- ER1_temp - EL1_temp
ediff2 <- ER2_temp - EL2_temp
if(ediff1 > 7 || ediff2 > 7){
next
}
i <- i + 1
}
# Gather data in matrix (observable)
Dobs <- matrix(0, nrow = n, ncol = 9)
colnames(Dobs) <- c("PairIndex", "SympOnset1", "SympOnset2", "SI",
"ExpoL1", "ExpoR1", "ExpoL2", "ExpoR2", "GI")
Dobs[,1] <- seq_len(n)
Dobs[,2] <- SO1
Dobs[,3] <- SO2
Dobs[,4] <- SO2-SO1
Dobs[,5] <- EL1
Dobs[,6] <- ER1
Dobs[,7] <- EL2
Dobs[,8] <- ER2
Dobs[,9] <- GI
Dobs <- data.frame(Dobs)
# Width of infection exposure windows
expoL <- c(Dobs$ExpoL1, Dobs$ExpoL2)
expoR <- c(Dobs$ExpoR1, Dobs$ExpoR2)
expo_diff <- expoR-expoL
exposure_mean <- mean(expo_diff)
Dobs_incub <- matrix(0, nrow = (2*n), ncol = 2)
colnames(Dobs_incub) <- c("tL","tR")
Dobs_incub[,1] <- c(Dobs$SympOnset1-Dobs$ExpoR1,
Dobs$SympOnset2-Dobs$ExpoR2)
Dobs_incub[,2] <- c(Dobs$SympOnset1-Dobs$ExpoL1,
Dobs$SympOnset2-Dobs$ExpoL2)
Dobs_incub <- Dobs_incub[sample(seq_len(2*n), replace = FALSE)[1:n],]
Dobs_incub <- as.data.frame(Dobs_incub)
# Check constraints on exposure intervals of Ferretti et al. (2020)
C1 <- all(ER1 <= apply(cbind(SO1,ER2), 1, min))
C2 <- all(EL2 >= EL1)
C3 <- all(SO2 >= ER2)
Ferretti_constr <- all(c(C1,C2,C3))
#-- Plot of generation interval and incubation time pdf
# Generation interval distribution
GTpdf <- stats::dweibull(tt, shape = GIshape, scale = GIscale)
meanGT <- GIscale * gamma(1 + 1 / GIshape)
sdGT <- sqrt((GIscale ^ 2) * (gamma(1 + 2 / GIshape) -
(gamma(1 + 1 / GIshape) ^ 2)))
# Plot
if(plotsim == TRUE){
graphics::par(mfrow = c(2,2))
plot(tt, GTpdf, type = "l", col = "black", lwd = 2,
xlab = "Time (days)", ylab = "Density", ylim = c(0,0.25))
graphics::rug(GI)
graphics::lines(tt, Ipdf, type = "l", col = "red", lwd = 2)
graphics::rug(I1, col = "red", side = 3)
graphics::legend("topright", c("Generation interval", "Incubation period"),
lwd = c(2,2), lty = c(1,1), col = c("black","red"), bty = "n")
graphics::hist(Dobs$SI, breaks = sqrt(n) + 5, freq = FALSE,
xlab = "Time (days)",
ylab = "Serial interval distribution", main = "")
graphics::lines(stats::density(Dobs$SI), type = "l",
col = "darkgreen", lwd = 2)
graphics::hist(Dobs$GI, breaks = sqrt(n) + 5, freq = FALSE,
xlab = "Time (days)",
ylab = "Generation interval distribution", main = "")
graphics::lines(stats::density(Dobs$GI), type = "l",
col = "purple", lwd = 2)
graphics::hist(expo_diff, xlab = "Time (days)", ylab = "",
main ="Duration of exposure", freq = TRUE, col = "cornflowerblue")
}
graphics::par(mfrow = c(1,1))
if (match.arg(incubdist) == "MixWeibull"){
outlist <- list(Dobs = Dobs,
Ferrconst = Ferretti_constr,
Constr = paste0("Ferretti constraints passed: ",
Ferretti_constr),
GIdist = paste0("Generation interval with mean: ",
round(meanGT,2), " days and standard deviation of: ",
round(sdGT,2), " days."),
Incubdist = paste0("Incubation distribution is a Weibull mixture"),
Expoduration = paste0("Mean duration of exposure window: ",
round(exposure_mean,2), " days."),
Expomean = exposure_mean,
Dobsincub = Dobs_incub,
coarseness = coarseness,
meanmix = meanmix,
sdmix = sdmix,
qmix = qmix,
Ipdf = Ipdf,
tdom = tt)
} else {
outlist <- list(Dobs = Dobs,
Ferrconst = Ferretti_constr,
Constr = paste0("Ferretti constraints passed: ",
Ferretti_constr),
GIdist = paste0("Generation interval with mean: ",
round(meanGT,2), " days and standard deviation of: ",
round(sdGT,2), " days."),
Incubdist = paste0("Incubation distribution is ",
match.arg(incubdist),"(",Iparams[1],",",Iparams[2],")",
" with mean ",
round(meanI,2), " days and standard deviation of: ",
round(sdI,2), " days."),
Iparams = Iparams,
Expoduration = paste0("Mean duration of exposure window: ",
round(exposure_mean,2), " days."),
Expomean = exposure_mean,
Dobsincub = Dobs_incub,
coarseness = coarseness,
tdom = tt)
}
return(outlist)
}
oswaldogressani/EpiLPS documentation built on Oct. 25, 2024, 8:15 p.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 incubsim
Documented in incubsim
#' Simulation of incubation times
#'
#' @description
#' This routine simulates symptom onset times, generation times and left and
#' right bounds of infecting exposure windows for source-recipient pairs based
#' on a given incubation distribution and a generation time distribution. The
#' generation time distribution is assumed to be a Weibull with shape 2.826 and
#' scale 5.665 (Ferretti et al. 2020). The incubation distribution can be chosen
#' by the user among:
#'
#' \itemize{
#' \item{A LogNormal distribution with a mean of 5.5 days and standard deviation of
#' 2.1 days (Ferretti et al. 2020).}
#' \item{A Weibull distribution (shape-scale parameterization) with a mean of 6.4 days and standard deviation of
#' 2.3 days (Backer et al. 2020).}
#' \item{An artificial bimodal distribution constructed from a mixture of two
#' Weibull distributions.}
#' \item{A Gamma distribution (shape-rate parameterization) with a mean of 3.8 days and standard deviation of
#' 2.9 days (Donnelly et al. 2003).}
#' }
#'
#' @usage incubsim(incubdist = c("LogNormal","Weibull","MixWeibull", "Gamma"), coarseness = 1,
#' n = 100, tmax = 20, tgridlen = 500, plotsim = FALSE)
#'
#' @param incubdist The distribution of the incubation period.
#' @param coarseness The average coarseness of the data. Default is 1 day.
#' @param n The sample size.
#' @param tmax The upper bound on which to evaluate the \code{incubdist} density.
#' @param plotsim Graphical visualization of the simulated data?
#' @param tgridlen The number of grid points on which to evaluate the density.
#'
#' @return A list including (among others) the left and right bounds of the
#' incubation period.
#'
#' @author Oswaldo Gressani \email{oswaldo_gressani@hotmail.fr}
#'
#' @references Ferretti, L., Wymant, C., Kendall, et al. (2020). Quantifying
#' SARS-CoV-2 transmission suggests epidemic control with digital contact
#' tracing. \emph{Science}, \strong{368}(6491), eabb6936.
#' @references Backer, J. A., Klinkenberg, D., & Wallinga, J. (2020). Incubation
#' period of 2019 novel coronavirus (2019-nCoV) infections among
#' travellers from Wuhan, China, 20–28 January 2020. \emph{Eurosurveillance},
#' \strong{25}(5), 2000062.
#' @references Donnelly, C. A., Ghani, A. C., Leung, G. M., et al. (2003).
#' Epidemiological determinants of spread of causal agent of severe acute
#' respiratory syndrome in Hong Kong. \emph{The Lancet}, \strong{361}(9371),
#' 1761-1766.
#'
#' @examples
#' simdat <- incubsim(n = 50) # Simulation of incubation times for 50 cases.
#' simdat$Dobsincub # Left and right bounds of incubation period.
#'
#' @export
incubsim <- function(incubdist = c("LogNormal","Weibull","MixWeibull", "Gamma"),
coarseness = 1, n = 100, tmax = 20, tgridlen = 500, plotsim = FALSE){
mixWeibull <- function(n, alpha, Ishape1 = 2.356, Iscale1 = 3.564,
Ishape2 = 9.351, Iscale2 = 12.563) {
#Probability weights in the mixture
w1 <- alpha
w2 <- 1 - alpha
#Sample scheme with replacement and weights on the elements being sampled
w <- sample(c(1, 2), n, replace = TRUE, prob = c(w1, w2))
#Sample of size n from the mixture of 2 Weibulls
sample <- (w == 1) * stats::rweibull(n = n, shape = Ishape1, scale = Iscale1) +
(w == 2) * stats::rweibull(n = n, shape = Ishape2, scale = Iscale2)
# set.seed(123)
# mixsample <- mixWeibull(n = 50000000, alpha = 0.5)
# summary(mixsample$sample)
# sd(mixsample$sample)
# quantile(mixsample$sample, probs = c(0.05,0.25,0.5,0.75,0.95))
outlist <- list(sample = sample, Ishape1 = Ishape1, Iscale1 = Iscale1,
Ishape2 = Ishape2, Iscale2 = Iscale2)
return(outlist)
}
# Domain of Generation time/Incubation time density
tt <- seq(0, tmax, length = tgridlen)
# Generation interval Weibull parameters #Ferretti et al. 2020
GIshape <- 2.826
# GIshape <- 3.4395
GIscale <- 5.665
if (match.arg(incubdist) == "LogNormal") { # Ferretti et al. 2020
Incubmean <- 1.644
Incubsd <- 0.363
if (coarseness == 1) {
dirforce <- 0.185
} else if (coarseness == 2) {
dirforce <- 0.375
}
Ipdf <- stats::dlnorm(tt, meanlog = Incubmean, sdlog = Incubsd)
meanI <- exp(Incubmean + 0.5 * (Incubsd ^ 2))
medianI <- exp(Incubmean)
sdI <- sqrt(exp(2 * Incubmean + Incubsd ^ 2) * (exp(Incubsd ^ 2) - 1))
Iparams <- c(Incubmean, Incubsd)
} else if (match.arg(incubdist) == "Weibull") { # Backer et al. 2020
Ishape <- 3.00060
Iscale <- 7.170345
if (coarseness == 1) {
dirforce <- 0.160
} else if (coarseness == 2) {
dirforce <- 0.320
}
Ipdf <- stats::dweibull(tt, shape = Ishape, scale = Iscale)
meanI <- Iscale * gamma(1 + 1 / Ishape)
medianI <- stats::qweibull(p = 0.5, shape = Ishape, scale = Iscale)
sdI <- sqrt((Iscale ^ 2) * ((gamma(1 + 2 / Ishape) - (gamma(1 + 1 / Ishape) ^ 2))))
Iparams <- c(Ishape, Iscale)
} else if (match.arg(incubdist) == "MixWeibull"){ # Artificial
if (coarseness == 1) {
dirforce <- 0.132
} else if (coarseness == 2) {
dirforce <- 0.264
}
mmix <- mixWeibull(n=1, alpha = 0.5)
Ishape1 <- mmix$Ishape1
Iscale1 <- mmix$Iscale1
Ishape2 <- mmix$Ishape2
Iscale2 <- mmix$Iscale2
Ipdf <- 0.5 * stats::dweibull(tt, shape = Ishape1, scale = Iscale1) +
0.5 * stats::dweibull(tt, shape = Ishape2, scale = Iscale2)
meanmix <- 0.5 * (Iscale1 * gamma(1 + 1 / Ishape1) +
Iscale2 * gamma(1 + 1 / Ishape2))
sdmix <- 4.622
qmix <- c(1.371, 1.885,2.301,2.680,3.050,3.434,3.856,4.361,5.075,7.191,
9.876,10.701,11.251,11.692,12.080,12.446,12.814,13.218,13.734)
} else if(match.arg(incubdist) == "Gamma") { # Donnelly 2003
Ishape <- 1.739646
Irate <- 0.4566
if (coarseness == 1) {
dirforce <- 0.262
} else if (coarseness == 2) {
dirforce <- 0.558
}
Ipdf <- stats::dgamma(tt, shape = Ishape, rate = Irate)
meanI <- Ishape / Irate
medianI <- stats::qgamma(p = 0.5, shape = Ishape, rate = Irate)
sdI <- sqrt(Ishape / (Irate ^ 2))
Iparams <- c(Ishape, Irate)
}
#-- Generation of symptom onset data and infection exposure intervals
GI <- c() # Vector to host generation times
SI <- c() # Vector to host serial interval times
TI1 <- c() # Vector to host infection times of source/infector
TI2 <- c() # Vector to host infection times of recipient/infectee
I1 <- c() # Vector to host incubation period for source/infector
I2 <- c() # Vector to host incubation period for recipient/infectee
SO1 <- c() # Vector to host symptom onset times of source/infector
SO2 <- c() # Vector to host symptom onset times of recipient/infectee
EL1 <- c() # Vector to host lower bound of exposure time of source/infector
ER1 <- c() # Vector to host upper bound of exposure time of source/infector
EL2 <- c() # Vector to host lower bound of exposure time of recipient/infectee
ER2 <- c() # Vector to host upper bound of exposure time of recipient/infectee
i <- 1
while(i <= n){
# Step 1: Draw generation time from its pdf (Weibull)
GI[i] <- stats::rweibull(n = 1, shape = GIshape, scale = GIscale)
# Step 2: Set time of infection of infector/source (0,365)
TI1[i] <- stats::runif(1, min = 0, max = 365)
# Step 3: Compute the time of infection of infectee/recipient
TI2[i] <- TI1[i] + GI[i]
# Step 4: Generate a pair of incubation periods from pdf (LogNormal)
if (match.arg(incubdist) == "LogNormal") {
I1[i] <- stats::rlnorm(1, meanlog = Incubmean, sdlog = Incubsd)
I2[i] <- stats::rlnorm(1, meanlog = Incubmean, sdlog = Incubsd)
} else if (match.arg(incubdist) == "Weibull") {
I1[i] <- stats::rweibull(1, shape = Ishape, scale = Iscale)
I2[i] <- stats::rweibull(1, shape = Ishape, scale = Iscale)
} else if (match.arg(incubdist) == "MixWeibull") {
I1[i] <- mixWeibull(n = 1, alpha = 0.5)$sample
I2[i] <- mixWeibull(n = 1, alpha = 0.5)$sample
} else if (match.arg(incubdist) == "Gamma") {
I1[i] <- stats::rgamma(1, shape = Ishape, rate = Irate)
I2[i] <- stats::rgamma(1, shape = Ishape, rate = Irate)
}
# Step 5: Compute times of symptom onset
SO1[i] <- TI1[i] + I1[i]
SO2[i] <- TI2[i] + I2[i]
SI[i] <- SO2[i] - SO1[i]
# Step 6: Compute lower/upper bounds of exposure times
minexpo1 <- TI1[i] - dirforce * I1[i]
minexpo2 <- TI2[i] - dirforce * I2[i]
exposure_draw <- function(tlow,tup){
norm_mean <- 0.5 * (tlow + tup)
sdnew <- 0
probtar <- 0.999
ptol <- stats::pnorm(q = tup, mean = norm_mean, sd = sdnew)
while(ptol >= probtar){
sdnew <- sdnew + 0.001
ptol <- stats::pnorm(q = tup, mean = norm_mean, sd = sdnew)
}
tdraw <- 1e+10
while(tdraw >= tup || tdraw <= tlow) {
tdraw <- stats::rnorm(n = 1, mean = norm_mean, sd = sdnew)
}
return(tdraw)
}
# Draw left bound of exposure for infector/infectee
EL1_temp <- exposure_draw(tlow = minexpo1, tup = TI1[i])
EL2_temp <- exposure_draw(tlow = minexpo2, tup = TI2[i])
while(EL2_temp < EL1_temp){# Ferretti et al. 2020 constraints
EL1_temp <- exposure_draw(tlow = minexpo1, tup = TI1[i])
EL2_temp <- exposure_draw(tlow = minexpo2, tup = TI2[i])
}
EL1[i] <- EL1_temp
EL2[i] <- EL2_temp
maxexpo1 <- TI1[i] + dirforce * I1[i]
maxexpo2 <- TI2[i] + dirforce * I2[i]
# Draw right bound of exposure for infector/infectee
ER1_temp <- exposure_draw(tlow = TI1[i], tup = maxexpo1)
ER2_temp <- exposure_draw(tlow = TI2[i], tup = maxexpo2)
while(ER1_temp > min(SO1[i],ER2_temp)){# Ferretti et al. 2020 constraints
ER1_temp <- exposure_draw(tlow = TI1[i], tup = maxexpo1)
ER2_temp <- exposure_draw(tlow = TI2[i], tup = maxexpo2)
}
ER1[i] <- ER1_temp
ER2[i] <- ER2_temp
ediff1 <- ER1_temp - EL1_temp
ediff2 <- ER2_temp - EL2_temp
if(ediff1 > 7 || ediff2 > 7){
next
}
i <- i + 1
}
# Gather data in matrix (observable)
Dobs <- matrix(0, nrow = n, ncol = 9)
colnames(Dobs) <- c("PairIndex", "SympOnset1", "SympOnset2", "SI",
"ExpoL1", "ExpoR1", "ExpoL2", "ExpoR2", "GI")
Dobs[,1] <- seq_len(n)
Dobs[,2] <- SO1
Dobs[,3] <- SO2
Dobs[,4] <- SO2-SO1
Dobs[,5] <- EL1
Dobs[,6] <- ER1
Dobs[,7] <- EL2
Dobs[,8] <- ER2
Dobs[,9] <- GI
Dobs <- data.frame(Dobs)
# Width of infection exposure windows
expoL <- c(Dobs$ExpoL1, Dobs$ExpoL2)
expoR <- c(Dobs$ExpoR1, Dobs$ExpoR2)
expo_diff <- expoR-expoL
exposure_mean <- mean(expo_diff)
Dobs_incub <- matrix(0, nrow = (2*n), ncol = 2)
colnames(Dobs_incub) <- c("tL","tR")
Dobs_incub[,1] <- c(Dobs$SympOnset1-Dobs$ExpoR1,
Dobs$SympOnset2-Dobs$ExpoR2)
Dobs_incub[,2] <- c(Dobs$SympOnset1-Dobs$ExpoL1,
Dobs$SympOnset2-Dobs$ExpoL2)
Dobs_incub <- Dobs_incub[sample(seq_len(2*n), replace = FALSE)[1:n],]
Dobs_incub <- as.data.frame(Dobs_incub)
# Check constraints on exposure intervals of Ferretti et al. (2020)
C1 <- all(ER1 <= apply(cbind(SO1,ER2), 1, min))
C2 <- all(EL2 >= EL1)
C3 <- all(SO2 >= ER2)
Ferretti_constr <- all(c(C1,C2,C3))
#-- Plot of generation interval and incubation time pdf
# Generation interval distribution
GTpdf <- stats::dweibull(tt, shape = GIshape, scale = GIscale)
meanGT <- GIscale * gamma(1 + 1 / GIshape)
sdGT <- sqrt((GIscale ^ 2) * (gamma(1 + 2 / GIshape) -
(gamma(1 + 1 / GIshape) ^ 2)))
# Plot
if(plotsim == TRUE){
graphics::par(mfrow = c(2,2))
plot(tt, GTpdf, type = "l", col = "black", lwd = 2,
xlab = "Time (days)", ylab = "Density", ylim = c(0,0.25))
graphics::rug(GI)
graphics::lines(tt, Ipdf, type = "l", col = "red", lwd = 2)
graphics::rug(I1, col = "red", side = 3)
graphics::legend("topright", c("Generation interval", "Incubation period"),
lwd = c(2,2), lty = c(1,1), col = c("black","red"), bty = "n")
graphics::hist(Dobs$SI, breaks = sqrt(n) + 5, freq = FALSE,
xlab = "Time (days)",
ylab = "Serial interval distribution", main = "")
graphics::lines(stats::density(Dobs$SI), type = "l",
col = "darkgreen", lwd = 2)
graphics::hist(Dobs$GI, breaks = sqrt(n) + 5, freq = FALSE,
xlab = "Time (days)",
ylab = "Generation interval distribution", main = "")
graphics::lines(stats::density(Dobs$GI), type = "l",
col = "purple", lwd = 2)
graphics::hist(expo_diff, xlab = "Time (days)", ylab = "",
main ="Duration of exposure", freq = TRUE, col = "cornflowerblue")
}
graphics::par(mfrow = c(1,1))
if (match.arg(incubdist) == "MixWeibull"){
outlist <- list(Dobs = Dobs,
Ferrconst = Ferretti_constr,
Constr = paste0("Ferretti constraints passed: ",
Ferretti_constr),
GIdist = paste0("Generation interval with mean: ",
round(meanGT,2), " days and standard deviation of: ",
round(sdGT,2), " days."),
Incubdist = paste0("Incubation distribution is a Weibull mixture"),
Expoduration = paste0("Mean duration of exposure window: ",
round(exposure_mean,2), " days."),
Expomean = exposure_mean,
Dobsincub = Dobs_incub,
coarseness = coarseness,
meanmix = meanmix,
sdmix = sdmix,
qmix = qmix,
Ipdf = Ipdf,
tdom = tt)
} else {
outlist <- list(Dobs = Dobs,
Ferrconst = Ferretti_constr,
Constr = paste0("Ferretti constraints passed: ",
Ferretti_constr),
GIdist = paste0("Generation interval with mean: ",
round(meanGT,2), " days and standard deviation of: ",
round(sdGT,2), " days."),
Incubdist = paste0("Incubation distribution is ",
match.arg(incubdist),"(",Iparams[1],",",Iparams[2],")",
" with mean ",
round(meanI,2), " days and standard deviation of: ",
round(sdI,2), " days."),
Iparams = Iparams,
Expoduration = paste0("Mean duration of exposure window: ",
round(exposure_mean,2), " days."),
Expomean = exposure_mean,
Dobsincub = Dobs_incub,
coarseness = coarseness,
tdom = tt)
}
return(outlist)
}
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.