| Title: | Estimation of R0 and Real-Time Reproduction Number from Epidemics |
|---|---|
| Description: | Estimation of reproduction numbers for disease outbreak, based on incidence data. The R0 package implements several documented methods. It is therefore possible to compare estimations according to the methods used. Depending on the methods requested by user, basic reproduction number (commonly denoted as R0) or real-time reproduction number (referred to as R(t)) is computed, along with a 95% Confidence Interval. Plotting outputs will give different graphs depending on the methods requested : basic reproductive number estimations will only show the epidemic curve (collected data) and an adjusted model, whereas real-time methods will also show the R(t) variations throughout the outbreak time period. Sensitivity analysis tools are also provided, and allow for investigating effects of varying Generation Time distribution or time window on estimates. |
| Authors: | Pierre-Yves Boelle, Thomas Obadia |
| Maintainer: | Thomas Obadia <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.3-0 |
| Built: | 2024-10-10 04:13:43 UTC |
| Source: | https://github.com/tobadia/R0 |
Find the best-fitting generation time distribution from a series of serial interval.
est.GT( infector.onset.dates = NULL, infectee.onset.dates = NULL, serial.interval = NULL, request.plot = FALSE, ... )est.GT( infector.onset.dates = NULL, infectee.onset.dates = NULL, serial.interval = NULL, request.plot = FALSE, ... )
infector.onset.dates |
Vector of dates for infector symptoms onset. |
infectee.onset.dates |
Vector of dates for infectee symptoms onset. |
serial.interval |
Vector of reported serial interval. |
request.plot |
Should data adjustment be displayed at the end? |
... |
Parameters passed to other functions (useful for hidden parameters of |
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.
A object of class R0.GT that complies with generation.time() distribution
requirements of the R0 package.
Pierre-Yves Boelle, Thomas Obadia
#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)#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 or for an incidence dataset using
the methods implemented in the R0 package.
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, ... )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, ... )
epid |
Object containing epidemic curve data. |
GT |
Generation time distribution from |
t |
Vector of dates at which incidence was observed. |
begin |
Begin date for estimation. Can be an integer or a date (YYYY-mm-dd or YYYY/mm/dd). |
end |
End date for estimation. Can be an integer or a date (YYYY-mm-dd or YYYY/mm/dd). |
date.first.obs |
Optional date of first observation, if |
time.step |
Optional. If date of first observation is specified, number of day between each incidence observation. |
AR |
Attack rate as a percentage from total population. |
pop.size |
Population size in which the incident cases were observed. See more details in |
S0 |
Initial proportion of the population considered susceptible. |
methods |
Vector of methods to be used for R/R0/Rt estimation. Must be provided as |
checked |
Internal flag used to check whether integrity checks were ran or not. |
... |
Parameters passed to inner functions. |
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.
A list with components:
estimates |
List containing all results from called methods. |
epid |
Epidemic curve. |
GT |
Generation Time distribution function. |
t |
Date vector. |
begin |
Begin date for estimation. |
end |
End date for estimation. |
Pierre-Yves Boelle, Thomas Obadia
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.
#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#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
Create an object of class GT representing a discretized Generation Time
distribution, that can be subsequently passed to estimation routines.
generation.time( type = c("empirical", "gamma", "weibull", "lognormal"), val = NULL, truncate = NULL, step = 1, first.half = TRUE, p0 = TRUE )generation.time( type = c("empirical", "gamma", "weibull", "lognormal"), val = NULL, truncate = NULL, step = 1, first.half = TRUE, p0 = TRUE )
type |
Type of distribution (can be any of |
val |
Vector of values used for the empirical distribution, or |
truncate |
Maximum extent of the GT distribution. |
step |
Time step used in discretization. |
first.half |
Boolean. When set to |
p0 |
Boolen. When set to |
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.
A list with components:
GT |
The probabilities for each time unit, starting at time 0. |
time |
The time at which probabilities are calculated. |
mean |
The mean of the discretized GT. |
sd |
The standard deviation of the discretized GT. |
Pierre-Yves Boelle, Thomas Obadia
#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")#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")
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.
data(Germany.1918)data(Germany.1918)
A vector of numeric values named with date of record.
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
Observed generation time distribution for children in household for the 2009 A/H1N1 influenza pandemic.
data(GT.chld.hsld)data(GT.chld.hsld)
A vector of unnamed numeric values.
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.
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.
data(H1N1.serial.interval)data(H1N1.serial.interval)
A vector of unnamed numeric values.
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.
When first records of incidence are unavailable, tries to impute censored cases to rebuild a longer epidemic vector.
impute.incid(CD.optim.vect, CD.epid, CD.R0, CD.GT)impute.incid(CD.optim.vect, CD.epid, CD.R0, CD.GT)
CD.optim.vect |
Vector of two elements ( |
CD.epid |
Original epidemic vector, output of |
CD.R0 |
Assumed R0 value for the original epidemic vector. |
CD.GT |
Generation time distribution to be used for computations. |
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.
A vector with both imputed incidence and source available data.
Pierre-Yves Boelle, Thomas Obadia
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.
inspect.data(incid, GT = NULL, t = NULL)inspect.data(incid, GT = NULL, t = NULL)
incid |
An object (vector, data.frame, list) storing incidence. |
GT |
Generation time distribution from |
t |
Vector of dates at which incidence was observed (optional). |
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.
No object is returned. Instead, warnings are thrown upon detecting inconsistences.
Pierre-Yves Boelle, Thomas Obadia
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.
integrity.checks( epid, GT, t, begin, end, date.first.obs, time.step, AR, S0, methods )integrity.checks( epid, GT, t, begin, end, date.first.obs, time.step, AR, S0, methods )
epid |
Epidemic dataset, expecting incidence counts in a varity of possible formats (see |
GT |
Generation time distribution from |
t |
Vector of dates at which incidence was observed. |
begin |
Begin date for estimation. Can be an integer or a date (YYYY-mm-dd or YYYY/mm/dd). |
end |
End date for estimation. Can be an integer or a date (YYYY-mm-dd or YYYY/mm/dd). |
date.first.obs |
Optional date of first observation, if |
time.step |
Optional. If date of first observation is specified, number of day between each incidence observation. |
AR |
Attack rate as a percentage from total population. |
S0 |
Initial proportion of the population considered susceptible. |
methods |
Vector of methods to be used for R/R0/Rt estimation. Must be provided as |
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.
A list with two components, begin and end.
Pierre-Yves Boelle, Thomas Obadia
Sensitivity analysis to estimate the variation of reproduction numbers according to the disitrbution of generation time.
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, ... )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, ... )
incid |
A vector of incident cases. |
GT.type |
Type of distribution for GT (see GT.R for details). |
GT.mean.range |
Range of mean values used for all GT distributions throughout the simulation. Must be provided as a vector. |
GT.sd.range |
Range of standard deviation values used for GT distributions. Must be provided as a vector. |
begin |
Vector of begin dates for the estimation of epidemic. |
end |
Vector of end dates for estimation of the epidemic. |
est.method |
Estimation method used for sensitivity analysis. |
t |
Dates vector to be passed to estimation function. |
date.first.obs |
Optional date of first observation, if t not specified. |
time.step |
Optional. If date of first observation is specified, number of day between each incidence observation. |
... |
Parameters passed to inner functions |
By using different Generation Time (GT) distribution, different estimates of the reproduction ratio can be analyzed.
A data.frame s.a with following components :
GT.type |
Type of distribution for GT. |
GT.mean |
Range of means used for tested GTs. |
GT.sd |
Range of standard deviations used for tested GTs. |
R |
Computed value for Reproduction Number given GT.type, GT.mean and GT.sd. |
CI.lowe |
The lower limit of 95% CI for R. |
CI.upper |
The upper limit of 95% CI for R. |
Pierre-Yves Boelle, Thomas Obadia
Sensitivity analysis to estimate the variation of reproduction numbers according to period over which the incidence is analyzed.
sa.time( incid, GT, begin = NULL, end = NULL, est.method, t = NULL, date.first.obs = NULL, time.step = 1, res = NULL, ... )sa.time( incid, GT, begin = NULL, end = NULL, est.method, t = NULL, date.first.obs = NULL, time.step = 1, res = NULL, ... )
incid |
A vector of incident cases. |
GT |
Generation time distribution from |
begin |
Vector of begin dates for the estimation of epidemic. |
end |
Vector of end dates for estimation of the epidemic. |
est.method |
Estimation method used for sensitivity analysis. |
t |
Dates vector to be passed to estimation function. |
date.first.obs |
Optional date of first observation, if t not specified. |
time.step |
Optional. If date of first observation is specified, number of day between each incidence observation. |
res |
If specified, will extract most of data from a |
... |
Parameters passed to inner functions |
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.
A list with components :
df |
data.frame object with all results from sensitivity analysis. |
df.clean |
the same object, with NA rows removed. Used only for easy export of results. |
mat.sen |
Matrix with values of R0 given begin (rows) and end (columns) dates. |
begin |
A range of begin dates in epidemic. |
end |
A range of end dates in epidemic. |
Pierre-Yves Boelle, Thomas Obadia
Sensitivity analysis of reproduction ratio using supported estimation methods.
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, ... )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, ... )
incid |
A vector of incident cases. |
GT |
Generation time distribution from |
begin |
Vector of begin dates for the estimation of epidemic. |
end |
Vector of end dates for estimation of the epidemic. |
est.method |
Estimation method used for sensitivity analysis. |
sa.type |
String argument to choose between |
res |
If specified, will extract most of data from a |
GT.type |
Type of distribution for GT (see GT.R for details). |
GT.mean.range |
Range of mean values used for all GT distributions throughout the simulation. Must be provided as a vector. |
GT.sd.range |
Range of standard deviation values used for GT distributions. Must be provided as a vector. |
t |
Dates vector to be passed to estimation function. |
date.first.obs |
Optional date of first observation, if t not specified. |
time.step |
Optional. If date of first observation is specified, number of day between each incidence observation. |
... |
Parameters passed to inner functions |
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.
A sensitivity analysis object of class R0.S with components depending on sensitivity analysis sa.type.
Pierre-Yves Boelle, Thomas Obadia
#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#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
Generates several epidemic curves with specified distribution and reproduction number.
sim.epid( epid.nb, GT, R0, epid.length, family, negbin.size = NULL, peak.value = 50 )sim.epid( epid.nb, GT, R0, epid.length, family, negbin.size = NULL, peak.value = 50 )
epid.nb |
Number of epidemics to be simulated (defaults to 1) |
GT |
Generation time distribution from |
R0 |
Basic reproduction number, in its core definition. |
epid.length |
Maximum length of the epidemic (cases infected after this length will be truncated). |
family |
Distribution of offspring. Can be either |
negbin.size |
If family is set to "negbin", sets the size parameter of the negative binomial distribution. |
peak.value |
Threashold value for incidence before epidemics begins decreasing |
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, implies that
.
A matrix with epidemics stored as columns (incidence count).
Pierre-Yves Boelle, Thomas Obadia
#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.#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.
Generates several epidemic curves with an individual-based model.
sim.epid.indiv(beta, Tmax, n = 1, family = "poisson", negbin.size = NULL)sim.epid.indiv(beta, Tmax, n = 1, family = "poisson", negbin.size = NULL)
beta |
Contact rate in the SEIR model. |
Tmax |
Maximum length of the epidemic (cases infected after this length will be truncated). |
n |
Number of epidemics to be simulated (defaults to 1) |
family |
Distribution of offspring. Can be either |
negbin.size |
If family is set to "negbin", sets the size parameter of the negative binomial distribution. |
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 and variance
, to allow for overdispersion.
A matrix with epidemics stored as columns (incidence count).
This is the exact function as used in the manuscript (Obadia et al., 2012).
Pierre-Yves Boelle, Thomas Obadia
Aggregate real-time estimates of (from the TD or even the SB
methods) over larger time windows, while still accounting for the generation
time distribution.
smooth.Rt(res, time.period)smooth.Rt(res, time.period)
res |
An object of class |
time.period |
Time period to be used for aggregation. |
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.
A list with components :
R |
The estimate of the reproduction ratio. |
conf.int |
The 95% confidence interval for the R estimate. |
GT |
Generation time distribution uised in the computation. |
epid |
Original or augmented epidemic data, depending whether impute.values is set to FALSE or TRUE. |
begin |
Starting date for the fit. |
begin.nb |
The number of the first day used in the fit. |
end |
The end date for the fit. |
end.nb |
The number of the las day used for the fit. |
data.name |
The name of the dataset used. |
call |
Call used for the function. |
method |
Method used for fitting. |
method.code |
Internal code used to designate method. |
Pierre-Yves Boelle, Thomas Obadia
#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#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