Nothing
# Name : GT
# Desc : Generate a discrete Generation Time distribution given a distribution
# type and specifications (mean, sd).
# Date : 2011/11/09
# Update : 2023/03/02
# Author : Boelle, Obadia
###############################################################################
#' @title
#' Generation Time distribution
#'
#' @description
#' Create an object of class `GT` representing a discretized Generation Time
#' distribution, that can be subsequently passed to estimation routines.
#'
#' @details
#' How the GT is discretized may have some impact on the shape of the distribution.
#' For example, the distribution may be discretized in intervals of 1 time step
#' starting at time 0, i.e. [0,1), [1,2), and so on. Or it may be discretized as
#' [0,0.5), [0.5, 1.5), ... (the default).
#'
#' If the GT is discretized from a given continuous distribution, the expected
#' duration of the Generation Time will be less than the nominal, it will be in
#' better agreement using the second discretization (default behavior).
#'
#' If p0 is set to `TRUE` (default), the generation time distribution is set to 0 for
#' day 0, meaning that the infectees generated by an infected individual will not
#' become incident on the same day.
#'
#' If no truncation is provided, the distribution will be truncated at 99.99%
#' probability.
#'
#' @param type Type of distribution (can be any of `"empirical"`, `"gamma"`, `"weibull"`, or `"lognormal"`)
#' @param val Vector of values used for the empirical distribution, or `c(mean, sd)` if parametric.
#' @param truncate Maximum extent of the GT distribution.
#' @param step Time step used in discretization.
#' @param first.half Boolean. When set to `TRUE` (default), the first probability is computed on a half period.
#' @param p0 Boolen. When set to `TRUE` the probability on day 0 is forced to 0.
#'
#' @return
#' A list with components:
#' \item{GT}{The probabilities for each time unit, starting at time 0.}
#' \item{time}{The time at which probabilities are calculated.}
#' \item{mean}{The mean of the discretized GT.}
#' \item{sd}{The standard deviation of the discretized GT.}
#'
#' @importFrom stats pweibull pgamma plnorm
#'
#' @export
#'
#' @example tests/generation.time.R
#'
#' @author Pierre-Yves Boelle, Thomas Obadia
# Function declaration
generation.time <- function(
type = c("empirical", "gamma", "weibull", "lognormal"),
val = NULL,
truncate = NULL,
step = 1,
first.half = TRUE,
p0 = TRUE
)
# Code
{
# all tests required
type=match.arg(type)
#if empirical, check values
if (type=="empirical") {
GT <- val
if (any(GT <0))
stop("Values in 'val' must be positive")
if (sum(GT) >1)
warning("Values will be standardized to sum to 1")
if (!is.null(truncate)) {
if (truncate < length(val)) {
warning(paste("Empirical distribution truncated at length ",truncate))
GT <- GT[1:truncate]
}
}
# if parametric
} else {
if (length(val)<2 )
stop("val= c(mean,sd) must be provided for parametric GT")
mean <- val[1]
sd <- val[2]
if (any(c(mean,sd) <= 0))
stop("'mean' and 'sd' must be positive")
if (is.null(truncate)) {
tmax <- ceiling(mean+10*sd); # sufficiently large
} else {
tmax <- truncate
}
if (first.half) {
t.scale <- c(0,0.5+c(0:tmax))
} else {
t.scale <- c(0:tmax)
}
if (type == "gamma") {
a <- mean*mean/(sd*sd)
s <- sd*sd / mean
GT <- diff(pgamma(t.scale,shape=a,scale=s))
#other cases; weibull, ...
} else if (type == "lognormal") {
meanlog <- log(mean^2/sqrt(mean^2+sd^2))
sdlog <- sqrt(2)*sqrt(log(sqrt(mean^2+sd^2)/mean))
GT <- diff(plnorm(t.scale,meanlog=meanlog,sdlog=sdlog))
} else if (type == "weibull") {
cv <- sd/(mean)
if (cv < 1e-06) {
nu <- cv/(sqrt(trigamma(1)) - cv * digamma(1))
shape <- 1/nu
scale <- (mean)/(1 + nu * digamma(1))
} else {
aa <- log(cv^2 + 1)
nu <- 2 * cv/(1 + cv)
repeat {
gb <- (lgamma(1 + 2 * nu) - 2 * lgamma(1 + nu) - aa)/(2 * (digamma(1 + 2 * nu) - digamma(1 + nu)))
nu <- nu - gb
if (abs(gb) < 1e-12) break
}
shape <- 1/nu
scale <- exp(log(mean) - lgamma(1 + nu))
}
GT <- diff(pweibull(t.scale,shape=shape,scale=scale))
}
if (is.null(truncate)) {
# truncate when GI distribution >0.9999
GT.cum <- cumsum(GT)
if(length(GT.cum[GT.cum > 0.9999]) != 0){
truncate <- (GT.cum > 0.9999)*(1:length(GT.cum))
truncate <- min(truncate[truncate > 0])
if (truncate == 0) warning(paste('provide truncate larger than ',mean + 10 * sd))
GT <- GT[1:truncate]
}
}
}
if (p0 == TRUE) GT[1]=0
time <- 0:(length(GT)-1)
GT <- GT/sum(GT)
mu <- sum(GT * time)
sigma <- sqrt(sum(GT*time^2) - mu^2)
return(structure(list(GT=GT,time=time,mean=mu,sd=sigma),class="R0.GT"))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.