2009 A/H1N1 observed generation time distribution

Description

Observed generation time distribution for children in household for the 2009 A/H1N1 influenza pandemic.

Usage

data(GT.chld.hsld)

Format

A vector of unnamed numeric values.

References

S. Cauchemez, A. Bhattarai, T.L. Marchbanks, R.P. Fagan, S. Ostroff, N.M. Ferguson, et al., Role of Social Networks in Shaping Disease Transmission During a Community Outbreak of 2009 H1N1 Pandemic Influenza, Pnas. 108 (2011) 2825-2830.

Germany 1918 dataset

Description

This dataset is used as an example in the R0 package to estimate reproduction ratios using the available methods. It is the emporal distribution of Spanish flu in Prussia, Germany, from the 1918-19 influenza epidemic.

Usage

data(Germany.1918)

Format

A vector of numeric values named with date of record.

References

Peiper O: Die Grippe-Epidemie in Preussen im Jahre 1918/19. Veroeffentlichungen aus dem Gebiete der Medizinalverwaltung 1920, 10: 417-479 (in German).
Nishiura H. Time variations in the transmissibility of pandemic influenza in Prussia, Germany, from 1918-19. In: Theoretical Biology and Medical Modelling

H1N1 serial interval

Description

Serial interval (in its core definition, i.e. the time delay between symptom onset between pairs of infectors and infectees), taken from traced cases of H1N1 viruses.

Vector of values that represents the time lag between symptoms onset for pairs of infector/infectee, for a dataset of complete traced cases. Each value accounts for a pair of infector/infectee. This serial interval is often substitued for the generation time distribution, as it is easier to observe.

Usage

data(H1N1.serial.interval)

Format

A vector of unnamed numeric values.

References

S. Cauchemez, A. Bhattarai, T.L. Marchbanks, R.P. Fagan, S. Ostroff, N.M. Ferguson, et al., Role of Social Networks in Shaping Disease Transmission During a Community Outbreak of 2009 H1N1 Pandemic Influenza, Pnas. 108 (2011) 2825-2830.

Estimate reproduction number from exponential growth rate

Description

Calculates the basic reproduction number R_{0} from the exponential growth rate obtained by either Poisson or log-linear regression.

Usage

R.from.r(r, GT)

Arguments

Details

For internal use. Called by est.R0.EG().

R_{0} is calculated as the inverse of the discretized Laplace transform, making use of the discretized generation time distirbution (output from generation.time()).

Value

A numeric value for R.

Note

The formula for the discretized Laplace transform is taken from Wallinga & Lipsitch (2007).

Author(s)

Pierre-Yves Boelle, Thomas Obadia

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.

Deviation between theorietical incidence and observed data.

Description

When first records of incidence are unavailable, tries to impute censored cases to rebuild longer epidemic vector.

Usage

censored.deviation(optim.vect, epid, R0, GT)

Arguments

Details

For internal use. Called by impute.incid().

This function is not intended for stand-alone use. It computes the difference between theoretical incidence and observed epidemics, for a given vector of initial values. To find the find best-fitting incidence values, this same vector must be optimized, to minimize the value returned by SCE.

Value

The deviation between E(N_{t}) and N_{t} : sum((E(Nt) - Nt)^2).

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Check incidence vector in the input

Description

Checks incid in the input. For internal use only.

Usage

check.incid(incid, t = NULL, date.first.obs = NULL, time.step = 1)

Arguments

Details

For internal use. Called by estimation methods to format incidence input.

check.incid() handles everything related to incidence content integrity. It is designed to generate an output which comply with estimation functions requirements. Epidemic data can be provided as an epitools object (see below), or as vectors (incidence, dates, or both).

When dates are provided, they can be in a separate t vector, or computed with the first value and a time step. In the end, the function returns a list with componentsepid and t values. If you plan on using estimation functions on their own (and not throught estimate.R()), be aware that any incorrect input format will result in erratic behavior and/or crash.

Object incid is either a list or data.frame. Expect items/columns ⁠$dates⁠ and/or ⁠$stratum3⁠. This is expected to work with objects created by the epitools package (tested with v0.5-6).

epitools::epicurve.dates() returns (among other things) a list with $dates object. This list gives incidence per day. Other epicurve methods return $dates along with a ⁠$<time_period>⁠ object and a ⁠$stratum3⁠, which contains respectively daily incidence data agregated by the given time period, and the same data with colnames that comply with R standard time notation.

E.g.: epitools::epicurve.weeks() returns ⁠$dates⁠, ⁠$weeks⁠ and ⁠$stratum3.⁠ ⁠$stratum3⁠ object is a list of dates (correct syntax), where each date is repeated to reflect the incidence value at this time.

Incidence data should not contain negative or missing values. Incidence data and time vector should have the same length.

Value

A list with components incid and t.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Examples

#Loading package
library(R0)

## Data is taken from the paper by Nishiura for key transmission parameters of an institutional
## outbreak during 1918 influenza pandemic in Germany
data(Germany.1918)
Germany.1918

## check.incid will extract names from the vector and coerce them as dates
check.incid(Germany.1918)

## Had Germany.1918 not have names() set, output would have been with index dates
## To force such an output, we here impose t=1:126. 
## Erasing names(Germany.1918) would have produced the same
## If so, then the epid$t vector returned will be replacement values.
check.incid(Germany.1918, t=1:126)

## You can also choose not to provide a complete date vector, but to only
## indicated the first day of the observation, and the number of days between each
## observation. In this example we will assume a time step of 7 days.
check.incid(Germany.1918, date.first.obs="1918-01-01", time.step=7)

## Finally, if no names() are available for the dataset and date.first.obs is not provided,
## setting time.step to any integer value will generate a t vector starting 
## from 1 and incrementing by the time.step parameter.

Estimate generation time distribution

Description

Find the best-fitting generation time distribution from a series of serial interval.

Usage

est.GT(
  infector.onset.dates = NULL,
  infectee.onset.dates = NULL,
  serial.interval = NULL,
  request.plot = FALSE,
  ...
)

Arguments

Details

Generation Time (GT) distribution can be estimated by two inputs methods. User can either provide two vectors of dates or a unique vector of reported serial intervals. If two vectors are provided, both ⁠*.onset.dates⁠ vectors should have the same length. The i-th element is the onset date for individual i. This implies that infector k (symptoms on day infector.onset.dates[k]) infected infectee k (symptoms on day infectee.onset.dates[k]). If only serial.interval is provided, each value is assumed to be the time elapsed between each pair of infector and infectee (i.e. the difference between the two ⁠*.onset.dates⁠ vectors).

When request.plot is set to TRUE, a graphical output provides standardized histogram of observed data along with the best-fitting adjusted model.

Value

A object of class R0.GT that complies with generation.time() distribution requirements of the R0 package.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Examples

#Loading package
library(R0)

# Data taken from traced cases of H1N1 viruses.
data(H1N1.serial.interval)
est.GT(serial.interval=H1N1.serial.interval)

## Best fitting GT distribution is a gamma distribution with mean = 3.039437 and sd = 1.676551 .
## Discretized Generation Time distribution
## mean: 3.070303 , sd: 1.676531 
## [1] 0.0000000000 0.1621208802 0.2704857362 0.2358751176 0.1561845680 0.0888997193 0.0459909903 
## 0.0222778094 0.0102848887 0.0045773285 0.0019791984 0.0008360608 0.0003464431 0.0001412594


# The same result can be achieved with two vectors of dates of onset.
# Here we use the same data, but trick the function into thinking onset dates are all "0".
data(H1N1.serial.interval)
est.GT(infector.onset.dates=rep(0,length(H1N1.serial.interval)), 
       infectee.onset.dates=H1N1.serial.interval)

Estimate R0 from attack rate of an epidemic

Description

Estimate basic reproduction number R0 from Attack Rate (AR) of an epidemic.

Usage

est.R0.AR(
  AR = NULL,
  incid = NULL,
  pop.size = NULL,
  S0 = 1,
  checked = FALSE,
  ...
)

Arguments

Details

For internal use. Called by estimate.R().

In the simple SIR model, the relation between R0 and the Attack Rate is in the form R0 = -ln((1-AR)/S0) / (AR - (1-S0)). If the population size is provided, the variance of R0 is estimated using the delta method. The hypothesis are that of homogeneous mixing, no more transmission (epidemic ended), no change in transmission or interventions during the epidemic. This estimate may be correct in closed populations, and may be less valid in other cases.

The correction for incomplete susceptibility is based on the SIR model equations.

A 95% confidence interval is computed for the attack rate considering the total population size (CI(AR) = AR +/- 1.96*sqrt(AR*(1-AR)/n)), and so the confidence interval for the reproduction number is computed with these extreme values.

Value

A list with components:

Note

This is the implementation of the formula by Dietz (1993).

Author(s)

Pierre-Yves Boelle, Thomas Obadia

References

Dietz, K. "The Estimation of the Basic Reproduction Number for Infectious Diseases." Statistical Methods in Medical Research 2, no. 1 (March 1, 1993): 23-41.

Examples

#Loading package
library(R0)

## Woodall reported an attack rate of 0.31 in a population of 1732 during
## the 1957 H2N2 influenza pandemic ('Age and Asian Influenza, 1957', BMJ, 1958)

est.R0.AR(pop.size=1732, AR=0.31)
# Reproduction number estimate using Attack Rate method
# R :  1.19698[ 1.179606 , 1.215077 ]

est.R0.AR(AR=0.31)
# Reproduction number estimate using  Attack Rate  method.
# R :  1.19698

est.R0.AR(pop.size=1732, incid=31)
# Reproduction number estimate using Attack Rate method
# R :  1.009057[ 1.005873 , 1.012269 ]

est.R0.AR(pop.size=1732, incid=c(2,3,4,7,4,2,4,5))
# Reproduction number estimate using Attack Rate method
# R :  1.009057[ 1.005873 , 1.012269 ]

est.R0.AR(pop.size=1732, incid=c(2,3,0,7,4,2,0,5))
# Reproduction number estimate using Attack Rate method
# R :  1.006699[ 1.003965 , 1.009453 ]

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.

Usage

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

Arguments

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.

Value

A list with components:

Note

This is the implementation of the method provided by Wallinga & Lipsitch (2007).

Author(s)

Pierre-Yves Boelle, Thomas Obadia

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.

Examples

#Loading package
library(R0)

## Data is taken from the paper by Nishiura for key transmission parameters of an institutional
## outbreak during 1918 influenza pandemic in Germany)

data(Germany.1918)
mGT<-generation.time("gamma", c(3, 1.5))

est.R0.EG(Germany.1918, mGT, begin=1, end=27)
## Reproduction number estimate using  Exponential Growth 
## R :  1.525895[ 1.494984 , 1.557779 ]

Estimate R0 by Maximum Likelihood

Description

Estimate R0 by maximum likelihood, assuming Poisson distribution of offsprings generated from infected individuals at each time step.

Usage

est.R0.ML(
  epid,
  GT,
  import = NULL,
  t = NULL,
  begin = NULL,
  end = NULL,
  date.first.obs = NULL,
  time.step = 1,
  range = c(0.01, 50),
  unknown.GT = FALSE,
  impute.incid = FALSE,
  checked = FALSE,
  ...
)

Arguments

Details

For internal use. Called by estimate.R().

White & Pagano (2009) detail two maximum likelihood methods for estimatig the reproduction ratio. The principle of the methods described by White & all is to compute the expected number of cases in the future, and optimise to get R using a Poisson distribution.

The first (and used by default in this package) assumes that the serial interval distirbution is known, and subsequently the likelihood is only maximised depending on the value of R.

The second method can be used if the serial interval distribution is unknown: in that case, the generation time is set to follow a Gamma distribution with two parameters (size, shape), and the optimization routine finds the values of R, size and shape that maximize the likelihood. However, the epidemic curve must be long enough to account for a whole generation. The authors showed that this is achieved when the cumulated amount of incident cases reaches approximately 150. When using this method, the flag unknown.GT must be set to TRUE. GT must still be provided with a R0.GT-class object, though its mean and sd will be recycled as starting value for the optimization routine.

The 95% confidence interval is achieved by profiling the likelihood.

In addition to the estimation method, we implemented a framework to impute unobserved incidence data, when the epidemic curve is found to not be available from the start of the outbreak. The details of this method are availble in the Supplementary Material S1 from doi:10.1186/1472-6947-12-147Obadia et al., 2012.

Value

A list with components:

Note

This is the implementation of the method provided by White & Pagano (2009).

Author(s)

Pierre-Yves Boelle, Thomas Obadia

References

White, L.F., J. Wallinga, L. Finelli, C. Reed, S. Riley, M. Lipsitch, and M. Pagano. "Estimation of the Reproductive Number and the Serial Interval in Early Phase of the 2009 Influenza A/H1N1 Pandemic in the USA." Influenza and Other Respiratory Viruses 3, no. 6 (2009): 267-276.

Examples

#Loading package
library(R0)

## Data is taken from paper by Nishiura for key transmission parameters of an institutional
## outbreak during the 1918 influenza pandemic in Germany)

data(Germany.1918)
mGT <- generation.time("gamma", c(2.45, 1.38))
est.R0.ML(Germany.1918, mGT, begin=1, end=27, range=c(0.01,50))
# Reproduction number estimate using  Maximum Likelihood  method.
# R :  1.307222[ 1.236913 , 1.380156 ]

res <- est.R0.ML(Germany.1918, mGT, begin=1, end=27, range=c(0.01,50))
plot(res)

## no change in R with varying range 
## (dates here are the same index as before. Just to illustrate different use)
est.R0.ML(Germany.1918, mGT, begin="1918-09-29", end="1918-10-25", range=c(0.01,100))
# Reproduction number estimate using  Maximum Likelihood  method.
# R :  1.307249[ 1.236913 , 1.380185 ]

Estimate the time dependent R using a Bayesian method

Description

Estimate R by a sequential Bayesian method, using known data up to a point in time as a Bayesian prior for the next iteration (see details).

Usage

est.R0.SB(
  epid,
  GT,
  t = NULL,
  begin = NULL,
  end = NULL,
  date.first.obs = NULL,
  time.step = 1,
  force.prior = FALSE,
  checked = FALSE,
  ...
)

Arguments

Details

For internal use. Called by estimate.R().

Initial prior is an unbiased uniform distribution for R, between 0 and the maximum of incid(t+1) - incid(t). For each subsequent iteration, a new distribution is computed for R, using the previous output as new prior.

The 95% confidence intervan is achieved by a cumulative sum of the posterior distirbution of R and corresponds to the 2.5-th and 97.5-th percentiles.

Value

A list with components:

Note

This is the implementation of the method provided by Bettencourt & Ribeiro (2008).

Author(s)

Pierre-Yves Boelle, Thomas Obadia

References

Bettencourt, L.M.A., and R.M. Ribeiro. "Real Time Bayesian Estimation of the Epidemic Potential of Emerging Infectious Diseases." PLoS One 3, no. 5 (2008): e2185.

Examples

#Loading package
library(R0)

## Data is taken from the paper by Nishiura for key transmission parameters of an institutional
## outbreak during 1918 influenza pandemic in Germany)

data(Germany.1918)
mGT <- generation.time("gamma", c(3,1.5))
SB <- est.R0.SB(Germany.1918, mGT)

## Results will include "most likely R(t)" (ie. the R(t) value for which the computed probability 
## is the highest), along with 95% CI, in a data.frame object
SB
# Reproduction number estimate using  Real Time Bayesian  method.
# 0 0 2.02 0.71 1.17 1.7 1.36 1.53 1.28 1.43 ...

SB$Rt.quant
# Date R.t. CI.lower. CI.upper.
# 1  1918-09-29 0.00      0.01      1.44
# 2  1918-09-30 0.00      0.01      1.42
# 3  1918-10-01 2.02      0.97      2.88
# 4  1918-10-02 0.71      0.07      1.51
# 5  1918-10-03 1.17      0.40      1.84
# 6  1918-10-04 1.70      1.09      2.24
# 7  1918-10-05 1.36      0.84      1.83
# 8  1918-10-06 1.53      1.08      1.94
# 9  1918-10-07 1.28      0.88      1.66
# 10 1918-10-08 1.43      1.08      1.77
# ...

## "Plot" will provide the most-likely R value at each time unit, along with 95CI
plot(SB)
## "Plotfit" will show the complete distribution of R for 9 time unit throughout the outbreak
plotfit(SB)

Estimate the Time-Dependent reproduction number

Description

Estimate the Time-Dependent reproduction number, R(t), as defined by Wallinga & Teunis, by exploring all possible pairs of infectors/infectees across likely transmission trees.

Usage

est.R0.TD(
  epid,
  GT,
  import = NULL,
  n.t0 = NULL,
  t = NULL,
  begin = NULL,
  end = NULL,
  date.first.obs = NULL,
  time.step = 1,
  q = c(0.025, 0.975),
  correct = TRUE,
  nsim = 10000,
  checked = FALSE,
  ...
)

Arguments

Details

For internal use. Called by estimate.R().

The confidence interval is computed by multinomial simulations at each time step, using the expected value of R.

Value

A list with components:

Note

This is the implementation of the method provided by Wallinga & Teunis (2004). Correction for estimation in real time is implemented as in Cauchemez et al, AJE (2006).

If imported cases are provided, they are counted in addition to autonomous cases. The final plot will show overall incidence.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

References

Wallinga, J., and Teunis P. "Different Epidemic Curves for Severe Acute Respiratory Syndrome Reveal Similar Impacts of Control Measures." American Journal of Epidemiology 160, no. 6 (2004): 509.
Cauchemez S., and Valleron AJ. "Estimating in Real Time the Efficacy of Measures to Control Emerging Communicable Diseases" American Journal of Epidemiology 164, no. 6 (2006): 591.

Examples

#Loading package
library(R0)

## Data is taken from the paper by Nishiura for key transmission parameters of an institutional
## outbreak during 1918 influenza pandemic in Germany)

data(Germany.1918)
mGT <- generation.time("gamma", c(3, 1.5))
TD <- est.R0.TD(Germany.1918, mGT, begin=1, end=126, nsim=100)
# Warning messages:
# 1: In est.R0.TD(Germany.1918, mGT) : Simulations may take several minutes.
# 2: In est.R0.TD(Germany.1918, mGT) : Using initial incidence as initial number of cases.
TD
# Reproduction number estimate using  Time-Dependent  method.
# 2.322239 2.272013 1.998474 1.843703 2.019297 1.867488 1.644993 1.553265 1.553317 1.601317 ...

## An interesting way to look at these results is to agregate initial data by longest time unit,
## such as weekly incidence. This gives a global overview of the epidemic.
TD.weekly <- smooth.Rt(TD, 7)
TD.weekly
# Reproduction number estimate using  Time-Dependant  method.
# 1.878424 1.580976 1.356918 1.131633 0.9615463 0.8118902 0.8045254 0.8395747 0.8542518 0.8258094..
plot(TD.weekly)

Estimate reproduction number (R0 or Rt) for one incidence dataset using available methods

Description

Estimate R_{0} or R(t) for an incidence dataset using the methods implemented in the R0 package.

Usage

estimate.R(
  epid = NULL,
  GT = NULL,
  t = NULL,
  begin = NULL,
  end = NULL,
  date.first.obs = NULL,
  time.step = 1,
  AR = NULL,
  pop.size = NULL,
  S0 = 1,
  methods = NULL,
  checked = TRUE,
  ...
)

Arguments

Details

Currently, supported methods are Exponential Growth (EG), Maximum Likelihood (ML), Attack Rate (AR), Time-Dependant (TD), and Sequential Bayesian (SB). The corresponding references from the literature are available below.

This function acts as a front-end and will prepare relevant inputs to pass them to internal estimation routines. In particular, all inputs will undergo validation through integrity.checks() and the checked flag (defaulting as TRUE here) will be passed to internal estimation routines. Any warning raised by integrity.checks() should warrant careful thinking and investigation.

Value

A list with components:

Author(s)

Pierre-Yves Boelle, Thomas Obadia

References

est.R0.AR() : Dietz, K. "The Estimation of the Basic Reproduction Number for Infectious Diseases." Statistical Methods in Medical Research 2, no. 1 (March 1, 1993): 23-41.
est.R0.EG() : 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.
est.R0.ML() : White, L.F., J. Wallinga, L. Finelli, C. Reed, S. Riley, M. Lipsitch, and M. Pagano. "Estimation of the Reproductive Number and the Serial Interval in Early Phase of the 2009 Influenza A/H1N1 Pandemic in the USA." Influenza and Other Respiratory Viruses 3, no. 6 (2009): 267-276.
est.R0.SB() : Bettencourt, L.M.A., and R.M. Ribeiro. "Real Time Bayesian Estimation of the Epidemic Potential of Emerging Infectious Diseases." PLoS One 3, no. 5 (2008): e2185.
est.R0.TD() : Wallinga, J., and P. Teunis. "Different Epidemic Curves for Severe Acute Respiratory Syndrome Reveal Similar Impacts of Control Measures." American Journal of Epidemiology 160, no. 6 (2004): 509; Cauchemez S., and Valleron AJ. "Estimating in Real Time the Efficacy of Measures to Control Emerging Communicable Diseases" American Journal of Epidemiology 164, no. 6 (2006): 591.

Examples

#Loading package
library(R0)

## Outbreak during 1918 influenza pandemic in Germany)
data(Germany.1918)
mGT <- generation.time("gamma", c(3, 1.5))
estR0 <- estimate.R(Germany.1918, mGT, begin=1, end=27, methods=c("EG", "ML", "TD", "AR", "SB"), 
                  pop.size=100000, nsim=100)

attributes(estR0)
## $names
## [1] "epid"      "GT"        "begin"     "end"       "estimates"
## 
## $class
## [1] "R0.sR"

## Estimates results are stored in the $estimates object
estR0
## Reproduction number estimate using  Exponential Growth  method.
## R :  1.525895[ 1.494984 , 1.557779 ]
## 
## Reproduction number estimate using  Maximum Likelihood  method.
## R :  1.383996[ 1.309545 , 1.461203 ]
## 
## Reproduction number estimate using  Attack Rate  method.
## R :  1.047392[ 1.046394 , 1.048393 ]
## 
## Reproduction number estimate using  Time-Dependent  method.
## 2.322239 2.272013 1.998474 1.843703 2.019297 1.867488 1.644993 1.553265 1.553317 1.601317 ...
## 
## Reproduction number estimate using  Sequential Bayesian  method.
## 0 0 2.22 0.66 1.2 1.84 1.43 1.63 1.34 1.52 ...


## If no date vector nor date of first observation are provided, results are the same
## except time values in $t are replaced by index

Poisson log-likelihood for an observed epidemic

Description

Computes the Poisson log-likelihood of an observed epidemic, compared to the expected incidence values for a given R and generation time distirbution.

Usage

fit.epid(log.R, epid, GT, import, pred = FALSE, offset = 0)

Arguments

Details

For internal use. Called from est.R0.ML().

Computes and returns the Poisson log-likelihood of an observed epidemic, conditional on a given value of R and a generation time distribution.

Value

The Poisson log-likelihood of the observed epidemic passed as argument epid, or the Poisson-predicted incidence series given R and GT values.

Note

This function is also used in fit.epid.optim(), which is a wrapper for the joint optimization routine to determine the best-fitting R and GT values when the est.R0.ML() method is called and asked to jointly estimate reproduction ratios and generation time distirbution.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Joint estimation of R and generation time distribution for the ML method

Description

An objective function, wrapped around fit.epid(), that will be passed to stats::optim() to jointly estimate the best-fitting R and GT distributions in the context of the Maximum-Likelihood method from est.R0.ML().

Usage

fit.epid.optim(par = c(1, 1, 1), ...)

Arguments

Details

For internal use. Called from est.R0.ML().

This is a wrapper function used to pass proper arguments to fit.epid() when the ML method is used to estimate simultaneously R and GT. This function is used as objective to maximize by stats::optim().

Value

A Poisson log-likelihood from fit.epid().

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Generation Time distribution

Description

Create an object of class GT representing a discretized Generation Time distribution, that can be subsequently passed to estimation routines.

Usage

generation.time(
  type = c("empirical", "gamma", "weibull", "lognormal"),
  val = NULL,
  truncate = NULL,
  step = 1,
  first.half = TRUE,
  p0 = TRUE
)

Arguments

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.

Value

A list with components:

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Examples

#Loading package
library(R0)

# GT for children at house(from Cauchemez PNAS 2011)

GT.chld.hsld1 <- generation.time("empirical", c(0,0.25,0.2,0.15,0.1,0.09,0.05,0.01))
plot(GT.chld.hsld1, col="green")
GT.chld.hsld1
# Discretized Generation Time distribution
# mean: 2.729412 , sd: 1.611636 
# [1] 0.00000000 0.29411765 0.23529412 0.17647059 0.11764706 0.10588235 0.05882353
# [8] 0.01176471

GT.chld.hsld2 <- generation.time("gamma", c(2.45, 1.38))
GT.chld.hsld2
# Discretized Generation Time distribution
# mean: 2.504038 , sd: 1.372760
# [1] 0.0000000000 0.2553188589 0.3247178420 0.2199060781 0.1144367560
# [6] 0.0515687896 0.0212246257 0.0082077973 0.0030329325 0.0010825594
#[11] 0.0003760069 0.0001277537


# GT for school & community
GTs1 <- generation.time("empirical", c(0,0.95,0.05))
plot(GTs1, col='blue')


plot(GT.chld.hsld1, ylim=c(0,0.5), col="red")
par(new=TRUE)
plot(GT.chld.hsld2, xlim=c(0,7), ylim=c(0,0.5), col="black")

Scaling of x-axis

Description

Internal scaling function to display proper x-axis labels.

Usage

get.scale(scale)

Arguments

Details

Builds the x-axis labels corresponding to a human-friendly level (day, week...).

Value

An integer corresponding to the number of days between each x-axis tickmark.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Impute censored cases to rebuild longer epidemic vector

Description

When first records of incidence are unavailable, tries to impute censored cases to rebuild a longer epidemic vector.

Usage

impute.incid(CD.optim.vect, CD.epid, CD.R0, CD.GT)

Arguments

Details

This function is not intended for stand-alone use. It optimizes the values of vect, based upon minimization of deviation between actual epidemics data and observed generation time. The optimized function is censored.deviation(), which returns the deviation used for minimization. Stand-alone use can be conducted, however this assumes data are all of the correct format.

Value

A vector with both imputed incidence and source available data.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Audit input data for common issues

Description

Inspect input data and look for common mistakes. The function does not return any object but yields warnings wupon detecting possible inconsistencies, along with suggestions as to how to clean inputs before running estimation routines.

Usage

inspect.data(incid, GT = NULL, t = NULL)

Arguments

Details

inspect.data() looks for common issues that could affect estimation routines. Such issues include too low incidence counts, leading/trailing zeros, non-integer values...

Before any checks are conducted, the data are passed to check.incid() to try and guess the format of the data.

A not-so-uncommon issue is to provide non-integer counts for incidence, for example when working with aggregated data that represent averaged number of cases across different communities. This however does not agree well with parametric likelihood that assume exponential growth over the early stage of an epidemic or Poisson distribution of cases, where non-integer values will cause calculations to fail.

Missing values may cause issues if not handled properly. By default, check.incid() will recast missing values to zero. Leading and trailing NA's should be omitted entirely from the input. Gaps found between available data may also cause issues if they span over a period that's longer than the total generation time. A warning is raised to inform on these possible issues.

Likewise, leading and tailing zeros would cause similar issues. Begin will default to the first value and end to the peak one. Just in case, these will be inspected here too. Sequence of 0s exceeding the length of the generation time will also yield a warning.

Scarce data may also cause errors when optimizing likelihood functions. A time-series of incidence spanning for a duration shorter than that of the generation time distribution is likely to correspond to an index case that hasn't yet infected all its offsprings. This would biais estimates downwards and should be taken into account when interpreting results.

Value

No object is returned. Instead, warnings are thrown upon detecting inconsistences.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Integrity checks for input parameters

Description

Before any requested estimation routine is ran, integrity.checks() is called to ensure the data passed as arguments meet the proper format and can be properly interpreted by subsequent functions.

Usage

integrity.checks(
  epid,
  GT,
  t,
  begin,
  end,
  date.first.obs,
  time.step,
  AR,
  S0,
  methods
)

Arguments

Details

For internal use. Called by all implemented estimation methods. All integrity/class checks are handled by this core function. GT must be an object of class R0.GT, and epidemic curve along with time values are handled here. If you plan on calling manually any other estimation function, make sure data are provided with correct format.

The epidemic curve epid may be provided as a vector. In that case, a vector t may be provided with the dates of observation. If t is not numeric, an attempt is made to convert to dates with as.Date(). If t is not provided, dates are obtained from the names of incid, and if not available, index values are used. Finally, one can provide an epidemic curve object generated by the epitools package (see check.incid() for more details).

A quick note on t, begin and end : when a date vector is provided (t), it will be used instead of index values to establish a date-related incidence. If no date vector is provided, then begin and end can still be forced to numeric values. It then links to the corresponding index values for incidence data. If a date vector is provided, begin and end can either be numeric values or dates. If numeric, they will link to the correspondig index values for incidence, and be afterward interpreted as the associated date. If date, they will be directly associated to incidence data.

Basicly, if specified, begin and end must always have the same class.

Value

A list with two components, begin and end.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot a generation time distribution

Description

Plots the distribution of a generation time object from generation.time().

Usage

## S3 method for class 'R0.GT'
plot(x, ...)

Arguments

Details

For internal use. Called by base::print().

Value

Called for side effect. Shows a plot of the generation time distribution.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot the R0/Rt value along with confidence interval

Description

Generates the graphical output for an estimated R or R(t). The graph depends on the estimation method. In general, the estimated value is presented along with its confidence interval.

Usage

## S3 method for class 'R0.R'
plot(x, xscale = "w", TD.split = FALSE, ...)

Arguments

Details

For internal use. This function is called when using plot.R0.sR().

This is a wrapper function that will call the correct plotting routine depending on which method was used to estimate a reproduction number.

Value

This function does not return any data. Called for side effect. Draws all R_{0} or /eqnR(t) values from one estimation method.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot of sensitivity analyses.

Description

Generates the graphical output for an object generated through sensitivity.analysis().

Usage

## S3 method for class 'R0.S'
plot(x, what = "heatmap", time.step = 1, skip = 5, ...)

Arguments

Details

For internal use. Called by base::plot() when applied to R0.S objects. A plot will be shown and the best model fit will be returned.

Value

A list with best R0 measure for each possible time period, along with corresponding begin/end dates.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot the R0/Rt value along with confidence interval

Description

Generates the graphical output for estimated R or R(t) fits.

Usage

## S3 method for class 'R0.sR'
plot(x, xscale = "w", TD.split = FALSE, ...)

Arguments

Details

For internal use. Called by [base::plot())].

Tweaked [base::plot()] function that draws the reproduction number values for each method contained in the object constructed by [estimate.R()].

[base::plot())]: R:base::plot()) [base::plot()]: R:base::plot() [estimate.R()]: R:estimate.R()

Value

This function does not return any data. Called for side effect. Draws all R_{0} or /eqnR(t) values from an output of estimate.R(), built with one or more methods.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot R0 for Attack Rate method

Description

Internal plot method for AR estimates.

Usage

plotRAR(x, ...)

Arguments

Details

This plots the R0 reproduction number along with confidence interval.

Value

This function does not return any data. Called for side effect.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot R0 for Exponential Growth method

Description

Internal plot method for EG estimates.

Usage

plotREG(x, ...)

Arguments

Details

This plots the exponential growth parameter r, along with R and corresponding confidence interval.

Value

This function does not return any data. Called for side effect.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot R0 for Maximum Likelihood method

Description

Internal plot method for ML estimates.

Usage

plotRML(x, ...)

Arguments

Details

This plots the R0 reproduction number along with confidence interval.

Value

This function does not return any data. Called for side effect.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot R0 for Sequential Bayesian method

Description

Internal plot method for SB estimates.

Usage

plotRSB(x, xscale, ...)

Arguments

Details

This plots the bayesian estimation of R_{0} along with its credible interval.

Value

This function does not return any data. Called for side effect.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot R0 for Time-Dependent method

Description

Internal plot method for TD estimates.

Usage

plotRTD(x, xscale, TD.split, ...)

Arguments

Details

This plots the time-dependent R(t) reproduction number along with confidence band.

Value

This function does not return any data. Called for side effect.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

S3 method for objects of class R0.R or R0.sR

Description

This method ensures compliance with CRAN checks etc. It is used to properly call plotfit.R0.R() and plotfit.R0.sR() and build corresponding plots.

Usage

plotfit(x, xscale = "w", SB.dist = TRUE, ...)

Arguments

Details

For internal use.

This method is designed to either call plotfit.R0.R() or plotfit.R0.sR(), and complies with S3 requirements at the time of writing.

It allows plotting the goodness-of-fit of an estimated model to the original epidemic curve provided the user.

Depending on the method of estimation, the graphical output will vary :

• EG, ML and TD methods will show the original epidemic curve, along with the best-fitting prediction model,

• AR will only show the epidemic curve, since no actual model is computed,

• SB will display 9 density curves for the R distribution throughout the epidemic.

Value

This function does not return any data.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot a model fit for R0.R objects

Description

Plot the fit of a single model output to epidemic data.

Usage

## S3 method for class 'R0.R'
plotfit(x, xscale = "w", SB.dist = TRUE, ...)

Arguments

Details

For internal use. This function is called by the plotfit() S3 method. Depending on the estimation method, either plotfitRxx(), plotfitRAR() or plotfitRSB() will be called.

Value

This function does not return any data. Called for side effect. Draws the fit of one estimation method to the data.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Plot a model fit for R0.sR objects

Description

Plot the fit of model outputs to epidemic data.

Usage

## S3 method for class 'R0.sR'
plotfit(x, xscale = "w", SB.dist = TRUE, ...)

Arguments

Details

For internal use. This function is called by the plotfit() S3 method. It will sequentialy call the necesary plotfit() methods for each estimation contained in the x argument.

Value

This function does not return any data. Called for side effect. Draws the fit of one or more estimation method to the data.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Internal plotfit method for AR estimates

Description

Plot the fit of a single model output to epidemic data when the method is AR.

Usage

plotfitRAR(x, xscale, ...)

Arguments

Details

For internal use. This function is called by the plotfit.R0.R().

Value

This function does not return any data. Called for side effect. Draws the fit of one estimation method to the data.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Internal plotfit method for AR estimates

Description

Plot the fit of a single model output to epidemic data when the method is SB.

Usage

plotfitRSB(x, xscale, SB.dist, ...)

Arguments

Details

For internal use. This function is called by the plotfit.R0.R().

Value

This function does not return any data. Called for side effect. Draws the fit of one estimation method to the data.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Internal plotfit method for EG, ML and TD estimates

Description

Plot the fit of a single model output to epidemic data when the method is EG, ML or TD.

Usage

plotfitRxx(x, xscale, ...)

Arguments

Details

For internal use. This function is called by the plotfit.R0.R().

Value

This function does not return any data. Called for side effect. Draws the fit of one estimation method to the data.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Print method for objects of class R0.GT

Description

Prints summary statistics for a generation time distribution from generation.time().

Usage

## S3 method for class 'R0.GT'
print(x, ...)

Arguments

Details

For internal use.

Displays the mean and standard deviation for generation time distributions created with generation.time(), along with the probabilities for each time interval.

Value

This function does not return any data. Called for side effects.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Print method for objects of class R0.R

Description

Prints summary statistics for an estimated reproduction ratio.

Usage

## S3 method for class 'R0.R'
print(x, ...)

Arguments

Details

For internal use.

Displays the estimated reproduction ratio along with its confidence interval. For the TD method, the time-series of R(t) is printed.

Value

This function does not return any data. Called for side effects.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Print method for objects of class R0.sR

Description

Prints summary statistics from estimate.R() outputs.

Usage

## S3 method for class 'R0.sR'
print(x, ...)

Arguments

Details

For internal use.

Displays the estimated reproduction ratio of one or more estimation results generated by estimate.R().

Value

This function does not return any data. Called for side effects.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Sensitivity of R0 to varying generation time distributions

Description

Sensitivity analysis to estimate the variation of reproduction numbers according to the disitrbution of generation time.

Usage

sa.GT(
  incid,
  GT.type,
  GT.mean.range,
  GT.sd.range,
  begin = NULL,
  end = NULL,
  est.method,
  t = NULL,
  date.first.obs = NULL,
  time.step = 1,
  ...
)

Arguments

Details

By using different Generation Time (GT) distribution, different estimates of the reproduction ratio can be analyzed.

Value

A data.frame s.a with following components :

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Sensitivity of R0 to time estimation windows

Description

Sensitivity analysis to estimate the variation of reproduction numbers according to period over which the incidence is analyzed.

Usage

sa.time(
  incid,
  GT,
  begin = NULL,
  end = NULL,
  est.method,
  t = NULL,
  date.first.obs = NULL,
  time.step = 1,
  res = NULL,
  ...
)

Arguments

Details

By varying different pairs of begin and end dates,different estimates of reproduction ratio can be analyzed.

'begin' and 'end' vector must have the same length for the sensitivity analysis to run. They can be provided either as dates or numeric values, depending on the other parameters (see check.incid()). If some begin/end dates overlap, they are ignored, and corresponding uncomputed data are set to NA.

Also, note that unreliable Rsquared values are achieved for very small time periods (begin ~ end). These values are not representative of the epidemic outbreak behaviour.

Value

A list with components :

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Sensitivity analysis of basic reproduction ratio to begin/end dates

Description

Sensitivity analysis of reproduction ratio using supported estimation methods.

Usage

sensitivity.analysis(
  incid,
  GT = NULL,
  begin = NULL,
  end = NULL,
  est.method = NULL,
  sa.type,
  res = NULL,
  GT.type = NULL,
  GT.mean.range = NULL,
  GT.sd.range = NULL,
  t = NULL,
  date.first.obs = NULL,
  time.step = 1,
  ...
)

Arguments

Details

This is a generic call function to use either sa.time or sa.GT. Argument must be chosen accordingly to sa.type. Please refer to sa.time() and sa.GT() for further details about arguments.

begin and end vectors must have the same length for the sensitivity analysis to run. They can be provided either as dates or numeric values, depending on the other parameters (see check.incid(). If some begin or end dates overlap, they are ignored and corresponding uncomputed data are set to NA. Also, note that unreliable Rsquared values are achieved for very small time periods (begin ~ end). These values are not representative of the epidemic outbreak behaviour.

Value

A sensitivity analysis object of class R0.S with components depending on sensitivity analysis sa.type.

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Examples

#Loading package
library(R0)

## Data is taken from the paper by Nishiura for key transmission parameters of an institutional
## outbreak during 1918 influenza pandemic in Germany)
data(Germany.1918)

## For this exemple, we use the exact same call as for the internal sensitivity analysis function

## sa.type = "GT"

## Here we will test GT with means of 1 to 5, each time with SD constant (1)
## GT and SD can be either fixed value or vectors of values
## Actual value in simulations may differ, as they are adapted according to the distribution type
tmp <- sensitivity.analysis(
     sa.type     = "GT", 
     incid       = Germany.1918, 
     GT.type     = "gamma", 
     GT.mean     = seq(1,5,1), 
     GT.sd.range = 1, 
     begin       = 1, 
     end         = 27, 
     est.method  = "EG")

## Results are stored in a matrix, each line dedicated to a (mean,sd) couple
plot(
     x    = tmp[,"GT.Mean"], 
     xlab = "mean GT (days)", 
     y    = tmp[,"R"], 
     ylim = c(1.2, 2.1), 
     ylab = "R0 (95% CI)", 
     type = "p", 
     pch  = 19, 
     col  = "black", 
     main = "Sensitivity of R0 to mean GT"
)
arrows(
     x0 = as.numeric(tmp[,"GT.Mean"]), 
     y0 = as.numeric(tmp[,"CI.lower"]), 
     y1 = as.numeric(tmp[,"CI.upper"]), 
     angle = 90, 
     code = 3, 
     col = "black", 
     length = 0.05
)
## One could tweak this example to change sorting of values (per mean, or per standard deviation)
## eg: 'x=tmp[,c('GT.Mean')]' could become 'x=tmp[,c('GT.SD')]'


## sa.type="time"

mGT <- generation.time("gamma", c(2.6,1))
sen <- sensitivity.analysis(
     sa.type    = "time", 
     incid      = Germany.1918, 
     GT         = mGT, 
     begin      = 1:10, 
     end        = 16:25, 
     est.method = "EG"
)
# ...
# Warning message:
# If 'begin' and 'end' overlap, cases where begin >= end are skipped.
# These cases often return Rsquared = 1 and are thus ignored.
## A list with different estimates of reproduction ratio, exponential growth rate and 95%CI 
## wtih different pairs of begin and end dates in form of data frame is returned.
## If method is "EG", results will include growth rate and deviance R-squared measure
## Else, if "ML" method is used, growth rate and R-squared will be set as NA

## Interesting results include the variation of R0 given specific begin/end dates.
## Such results can be plot as a colored matrix and display Rsquared=f(time period)
plot(sen, what=c("criterion","heatmap"))
## Returns complete data.frame of best R0 value for each time period 
## (allows for quick visualization)
## The "best.fit" is the time period over which the estimate is the more robust

# $best.fit
#    Time.period Begin.dates  End.dates       R Growth.rate  Rsquared CI.lower. CI.upper.
# 92          15  1970-01-08 1970-01-23 1.64098   0.1478316 0.9752564  1.574953  1.710209

Epidemic outbreak simulation

Description

Generates several epidemic curves with specified distribution and reproduction number.

Usage

sim.epid(
  epid.nb,
  GT,
  R0,
  epid.length,
  family,
  negbin.size = NULL,
  peak.value = 50
)

Arguments

Details

This function is only used for simulation purposes. The output is a matrix of n columns (number of outbreaks) by m rows (maximum length of an outbreak).

When using rnbinom with mean and size moments, the variance is given by mean + mean^2/size (see stats::rnbinom()). One should determine the size accordingly to the R0 value to increase the dispersion. From the previous formula for the variance, Var(X) = k*R_{0} implies that size = R0/(k-1).

Value

A matrix with epidemics stored as columns (incidence count).

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Examples

#Loading package
library(R0)

## In this example we simulate n=100 epidemic curves, with peak value at 150 incident cases, 
## and maximum epidemic length of 30 time units.
## Only the outbreak phase is computed. When the peak value is reached, the process is stopped 
## and another epidemic is generated.
sim.epid(epid.nb=100, GT=generation.time("gamma",c(3,1.5)), R0=1.5, 
         epid.length=30, family="poisson", peak.value=150)

# Here, a 30*100 matrix is returned. Each column is a single epidemic.

Influenza-like illness simulation (individual-based model)

Description

Generates several epidemic curves with an individual-based model.

Usage

sim.epid.indiv(beta, Tmax, n = 1, family = "poisson", negbin.size = NULL)

Arguments

Details

The epidemic is simulated using a branching process, with infinite number of susceptibles to allow for exponential growth. The model used follows the Crump-Mode-Jagers description, with S/E/I/R description of the natural history.

Latent and infectious period follow parametrized Gamma distributions typical of influenza. An index case is first introduced, and offspring is sampled from a negative binomial distribution, with mean beta*I and variance negbin.size*beta*I, to allow for overdispersion.

Value

A matrix with epidemics stored as columns (incidence count).

Note

This is the exact function as used in the manuscript (Obadia et al., 2012).

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Smooth real-time reproduction number over larger time periods

Description

Aggregate real-time estimates of R(t) (from the TD or even the SB methods) over larger time windows, while still accounting for the generation time distribution.

Usage

smooth.Rt(res, time.period)

Arguments

Details

Regrouping Time-Dependant R(t) values, or even Real Time Bayesian most-likely R values (according to R distributions) should take into account the generation time. Results can be plotted exactly the same was as input estimations, except they won't show any goodness-of-fit curve.

Value

A list with components :

Author(s)

Pierre-Yves Boelle, Thomas Obadia

Examples

#Loading package
library(R0)

## This script allows for generating a new estimation for RTB and TD methods.
## Estimations used as input are agregated by a time period provided by user.
## Results can be plotted exactly the same was as input estimations,
## except they won't show any goodness of fit curve.
data(Germany.1918)
mGT <- generation.time("gamma", c(3,1.5))
TD <- estimate.R(Germany.1918, mGT, begin=1, end=126, methods="TD", nsim=100)
TD
# Reproduction number estimate using  Time-Dependant  method.
# 2.322239 2.272013 1.998474 1.843703 2.019297 1.867488 1.644993 1.553265 1.553317 1.601317 ...
TD$estimates$TD$Rt.quant
#     Date      R.t. CI.lower.  CI.upper.
# 1      1 2.3222391 1.2000000  2.4000000
# 2      2 2.2720131 2.7500000  6.2500000
# 3      3 1.9984738 2.7500000  6.5000000
# 4      4 1.8437031 0.7368421  1.5789474
# 5      5 2.0192967 3.1666667  6.1666667
# 6      6 1.8674878 1.6923077  3.2307692
# 7      7 1.6449928 0.8928571  1.6428571
# 8      8 1.5532654 1.3043478  2.2608696
# 9      9 1.5533172 1.0571429  1.7428571
# 10    10 1.6013169 1.6666667  2.6666667
# ...

TD.weekly <- smooth.Rt(TD$estimates$TD, 7)
TD.weekly
# Reproduction number estimate using  Time-Dependant  method.
# 1.878424 1.580976 1.356918 1.131633 0.9615463 0.8118902 0.8045254 0.8395747 0.8542518 0.8258094..

TD.weekly$Rt.quant
#    Date      R.t. CI.lower. CI.upper.
# 1     1 1.8784240 1.3571429 2.7380952
# 2     8 1.5809756 1.3311037 2.0100334
# 3    15 1.3569175 1.1700628 1.5308219
# 4    22 1.1316335 0.9961229 1.2445302
# 5    29 0.9615463 0.8365561 1.0453074
# 6    36 0.8118902 0.7132668 0.9365193
# 7    43 0.8045254 0.6596685 0.9325967
# 8    50 0.8395747 0.6776557 1.0402930
# 9    57 0.8542518 0.6490251 1.1086351
# 10   64 0.8258094 0.5836735 1.1142857
# 11   71 0.8543877 0.5224719 1.1460674
# 12   78 0.9776385 0.6228070 1.4912281
# 13   85 0.9517133 0.5304348 1.3652174
# 14   92 0.9272833 0.5045045 1.3423423
# 15   99 0.9635479 0.4875000 1.5125000
# 16  106 0.9508951 0.5000000 1.6670455
# 17  113 0.9827432 0.5281989 1.8122157
# 18  120 0.5843895 0.1103040 0.9490928