R/est.R0.EG.R

In tobadia/R0: Estimation of R0 and Real-Time Reproduction Number from Epidemics

# Name   : est.R0.EG
# Desc   : Estimation of basic Reproduction Number using Exponential Growth method,
#          as presented by Wallinga & Lipsitch
# Date   : 2011/11/09
# Update : 2023/03/01
# Author : Boelle, Obadia
###############################################################################


#' @title
#' Estimate R0 from exponential growth rate of an epidemic
#' 
#' @description
#' Estimate R0 from the initial phase of an epidemic when incident cases are 
#' presumed to follow an exponential distribution.
#' 
#' @note
#' This is the implementation of the method provided by Wallinga & Lipsitch (2007).
#' 
#' @references
#' Wallinga, J., and M. Lipsitch. "How Generation Intervals Shape the Relationship Between Growth Rates and Reproductive Numbers." Proceedings of the Royal Society B: Biological Sciences 274, no. 1609 (2007): 599.
#' 
#' @details
#' For internal use. Called by [estimate.R()].
#' 
#' Method "poisson" uses Poisson regression of incidence to estimate the exponential 
#' growth rate. Method "linear" uses linear regression of log(incidence) against time.
#' 
#' The 95% confidence interval is computed from the 1/M(-r) formula, using bounds 
#' on r from the Poisson or linear regression.
#' 
#' @param epid Object containing epidemic curve data.
#' @param GT Generation time distribution from [generation.time()].
#' @param t Vector of dates at which incidence was observed.
#' @param begin At what time estimation begins.
#' @param end Time at which to end computation.
#' @param date.first.obs Optional date of first observation, if `t` not specified
#' @param time.step Optional. If date of first observation is specified, number of day between each incidence observation
#' @param reg.met Regression method used. Default is "poisson" (for GLM), but can be forced to "linear".
#' @param checked Internal flag used to check whether integrity checks were ran or not.
#' @param ... Parameters passed to inner functions.
#' 
#' @return
#' A list with components:
#' \item{R}{The estimate of the reproduction ratio.}
#' \item{conf.int}{The 95% confidence interval for the R estimate.}
#' \item{r}{Exponential growth rate of the epidemic.}
#' \item{conf.int.r}{Confidence interval of the exponential growth rate of the epidemic.}
#' \item{Rsquared}{The deviance R-squared measure for the considered dates and model.}
#' \item{epid}{Original epidemic data.}
#' \item{GT}{Generation time distribution used in the computation.}
#' \item{data.name}{Name of the data used in the fit.}
#' \item{begin}{Starting date for the fit.}
#' \item{begin.nb}{The number of the first day used in the fit.}
#' \item{end}{The end date for the fit.}
#' \item{end.nb}{The number of the las day used for the fit.}
#' \item{fit}{Method used for fitting.}
#' \item{pred}{Prediction on the period used for the fit.}
#' \item{method}{Method for estimation.}
#' \item{method.code}{Internal code used to designate method.}
#' 
#' @importFrom stats poisson predict lm glm coefficients
#' 
#' @export
#'
#' @keywords internal
#' 
#' @example tests/est.R0.EG.R
#' 
#' @author Pierre-Yves Boelle, Thomas Obadia



# Function declaration

est.R0.EG <-function( 
    epid, 
    GT, 
    t              = NULL, 
    begin          = NULL, 
    end            = NULL, 
    date.first.obs = NULL, 
    time.step      = 1, 
    reg.met        = "poisson", 
    checked        = FALSE,  
    ...
)


# Code

{
  DNAME <-  deparse(substitute(epid))
  #CALL = sapply(match.call()[-1], deparse)
  CALL <- match.call()
  
  # Various class and integrity check
  if (checked == FALSE) {
    parameters <- integrity.checks(epid=epid, GT=GT, t=t, begin=begin, end=end, date.first.obs=date.first.obs, time.step=time.step, AR=NULL, S0=NULL, methods="EG")
    begin <- parameters$begin
    end <- parameters$end
  }
  epid <- check.incid(epid, t, date.first.obs, time.step)
  begin.nb <- which(epid$t == begin)
  end.nb <- which(epid$t == end)
  
  
  
  #Backup original epidemic data
  epid.orig <- epid
  
  #Epidemic data used for fit is truncted to keep only values within [begin,end]
  #epid <- list(incid=epid$incid[begin.nb:end.nb], t=epid$t[begin.nb:end.nb], t.glm=seq(from=begin.nb, to=end.nb, by=1))
  #epid <- list(incid=epid$incid[begin.nb:end.nb], t=epid$t[begin.nb:end.nb], t.glm=epid$t[begin.nb:end.nb])
  epid <- list(incid=epid$incid[begin.nb:end.nb], t=epid$t[begin.nb:end.nb])
  
  # Different methods to estimate epidemic growth rate (r) from data
  # Method 1: Linear regression
  
  reg.met <- match.arg(reg.met)
  if (reg.met == "linear") {
    #tmp <-lm((log(incid)) ~ t.glm, data=epid)
    tmp <-lm((log(incid)) ~ t, data=epid)
    Rsquared = summary(tmp)$r.squared
    r <- coefficients(tmp)[2]
    conf.int <- confint(tmp)[2,]
    pred <- exp(predict(tmp,type="response"))
  } 
  
  # Method 2: Poisson regression
  else if (reg.met == "poisson") {
    #tmp <- glm(incid ~ t.glm, family=poisson(), data=epid)
    tmp <- glm(incid ~ t, family=poisson(), data=epid)
    Rsquared = (tmp$null.deviance-tmp$deviance)/(tmp$null.deviance)
    r <- coefficients(tmp)[2]
    confint <- confint(tmp)[2,]
    pred <- predict(tmp,type="response")
  }
  
  # GT represents Generation Time distribution -> #WE DEFINE A GI type SO THAT NO CHECK REQUIRED BUT type/class
  # Apply method of Wallinga / Lipsitch for discretized Laplace transform
  
  R <- as.numeric(R.from.r(r,GT))
  R.inf <- as.numeric(R.from.r(confint[1],GT))
  R.sup <- as.numeric(R.from.r(confint[2],GT))
  
  #return(structure(list(R=R, conf.int = c(R.inf,R.sup), r=r, conf.int.r=confint,epid=epid.orig, GT=GT, data.name =DNAME, begin=begin,end=end,method="Exponential Growth",pred=pred,fit=met),class="R"))
  return(structure(list(R=R, conf.int = c(R.inf,R.sup), r=r, conf.int.r=confint, Rsquared=Rsquared, epid=epid.orig, GT=GT, data.name=DNAME, call=CALL, begin=begin, begin.nb=begin.nb, end=end, end.nb=end.nb, method="Exponential Growth",pred=pred,fit=reg.met, method.code="EG"),class="R0.R"))
  
}