-
Notifications
You must be signed in to change notification settings - Fork 2
/
simulate_chains.Rd
188 lines (169 loc) · 7.29 KB
/
simulate_chains.Rd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/simulate.R
\name{simulate_chains}
\alias{simulate_chains}
\title{Simulate transmission chains}
\usage{
simulate_chains(
n_chains,
statistic = c("size", "length"),
offspring_dist,
...,
stat_threshold = Inf,
pop = Inf,
percent_immune = 0,
generation_time = NULL,
t0 = 0,
tf = Inf
)
}
\arguments{
\item{n_chains}{Number of chains to simulate.}
\item{statistic}{The chain statistic to track as the
stopping criteria for each chain being simulated when \code{stat_threshold} is not
\code{Inf}; A \verb{<string>}. It can be one of:
\itemize{
\item "size": the total number of cases produced by a chain before it goes
extinct.
\item "length": the total number of generations reached by a chain before
it goes extinct.
}}
\item{offspring_dist}{Offspring distribution: a \verb{<function>} like the ones
provided by R to generate random numbers from given distributions (e.g.,
\code{\link{rpois}} for Poisson). More specifically, the function needs to
accept at least one argument, \code{n}, which is the number of random
numbers to generate. It can accept further arguments, which will be passed
on to the random number generating functions. Examples that can be provided
here are \code{rpois} for Poisson distributed offspring, \code{rnbinom} for negative
binomial offspring, or custom functions.}
\item{...}{Parameters of the offspring distribution as required by R.}
\item{stat_threshold}{A stopping criterion for individual chain simulations;
a positive number coercible to integer. When any chain's cumulative statistic
reaches or surpasses \code{stat_threshold}, that chain ends. Defaults to \code{Inf}.
For example, if \code{statistic = "size"} and \code{stat_threshold = 10}, then any
chain that produces 10 or more cases will stop. Note that setting
\code{stat_threshold} does not guarantee that all chains will stop at the same
value.}
\item{pop}{Population size; An \verb{<Integer>}. Used alongside \code{percent_immune}
to define the susceptible population. Defaults to \code{Inf}.}
\item{percent_immune}{Percent of the population immune to
infection at the start of the simulation; A \verb{<numeric>} between 0 and 1.
Used alongside \code{pop} to initialise the susceptible population. Defaults to
0.}
\item{generation_time}{The generation time function; the name
of a user-defined named or anonymous function with only one argument \code{n},
representing the number of generation times to sample.}
\item{t0}{Start time (if generation time is given); either a single value
or a vector of same length as \code{n_chains} (number of initial cases) with
corresponding initial times. Defaults to 0, meaning all cases started at
time 0.}
\item{tf}{A number for the cut-off for the infection times (if generation
time is given); Defaults to \code{Inf}.}
}
\value{
An \verb{<epichains>} object, which is basically a \verb{<data.frame>}
with columns:
\itemize{
\item \code{chain} - an ID for active/ongoing chains,
\item \code{infectee} - a unique ID for each infectee.
\item \code{infector} - an ID for the infector of each infectee.
\item \code{generation} - a discrete time point during which infection occurs, and
optionally,
\item \code{time} - the time of infection.
}
}
\description{
It generates independent transmission chains starting with a single case
per chain, using a simple branching process model (See details for
definition of "chains" and assumptions). Offspring for each
chain are generated with an offspring distribution, and an optional
generation time distribution function.
The individual chain simulations are controlled by customisable stopping
criteria, including a threshold chain size or length, and a generation time
cut off. The function also optionally accepts population related inputs
such as the population size (defaults to Inf) and percentage of the
population initially immune (defaults to 0).
}
\section{Definition of a transmission chain}{
A transmission chain as used here is an independent case and all
the secondary cases linked to it through transmission. The chain starts with
a single case, and each case in the chain generates secondary cases
according to the offspring distribution. The chain ends when no more
secondary cases are generated.
}
\section{Calculating chain sizes and lengths}{
The function simulates the chain size for chain \eqn{i} at time
\eqn{t}, \eqn{I_{t, i}}, as:
\deqn{I_{t, i} = \sum_{i}^{I_{t-1}}X_{t, i},}
and the chain length/duration for chain \eqn{i} at time \eqn{t},
\eqn{L_{t, i}}, as:
\deqn{L_{t, i} = {\sf min}(1, X_{t, i}), }
where \eqn{X_{t, i}} is the secondary cases generated by chain \eqn{i}
at time \eqn{t}, and \eqn{I_{0, i} = L_{0, i} = 1}.
The distribution of secondary cases, \eqn{X_{t, i}} is modelled by the
offspring distribution (\code{offspring_dist}).
}
\section{Specifying \code{generation_time}}{
The argument \code{generation_time} must be specified as a function with
one argument, \code{n}.
For example, assuming we want to specify the generation time as a random
log-normally distributed variable with \code{meanlog = 0.58} and \code{sdlog = 1.58},
we could define a named function, let's call it "generation_time_fn",
with only one argument representing the number of generation times to
sample: \code{generation_time_fn <- function(n){rlnorm(n, 0.58, 1.38)}},
and assign the name of the function to \code{generation_time} in
the simulation function, i.e.
\code{`simulate_*`(..., generation_time = generation_time_fn)},
where \code{...} are the other arguments to \verb{simulate_*()} and * is a placeholder
for the rest of simulation function's name.
Alternatively, we could assign an anonymous function to \code{generation_time}
in the \verb{simulate_*()} call, i.e.
\code{simulate_*(..., generation_time = function(n){rlnorm(n, 0.58, 1.38)})}
OR \code{simulate_*(..., generation_time = \(n){rlnorm(n, 0.58, 1.38)})},
where \code{...} are the other arguments to \verb{simulate_*()}.
}
\examples{
# Using a Poisson offspring distribution and simulating from an infinite
# population up to chain size 10.
set.seed(32)
chains_pois_offspring <- simulate_chains(
n_chains = 10,
statistic = "size",
offspring_dist = rpois,
stat_threshold = 10,
generation_time = function(n) rep(3, n),
lambda = 2
)
chains_pois_offspring
# Using a Negative binomial offspring distribution and simulating from a
# finite population up to chain size 10.
set.seed(32)
chains_nbinom_offspring <- simulate_chains(
n_chains = 10,
pop = 100,
percent_immune = 0,
statistic = "size",
offspring_dist = rnbinom,
stat_threshold = 10,
generation_time = function(n) rep(3, n),
mu = 2,
size = 0.2
)
chains_nbinom_offspring
}
\references{
Jacob C. (2010). Branching processes: their role in epidemiology.
International journal of environmental research and public health, 7(3),
1186–1204. \doi{https://doi.org/10.3390/ijerph7031204}
Blumberg, S., and J. O. Lloyd-Smith. 2013. "Comparing Methods for
Estimating R0 from the Size Distribution of Subcritical Transmission
Chains." Epidemics 5 (3): 131–45.
\doi{https://doi.org/10.1016/j.epidem.2013.05.002}.
Farrington, C. P., M. N. Kanaan, and N. J. Gay. 2003.
"Branching Process Models for Surveillance of Infectious Diseases
Controlled by Mass Vaccination.” Biostatistics (Oxford, England)
4 (2): 279–95. \doi{https://doi.org/10.1093/biostatistics/4.2.279}.
}
\author{
James M. Azam, Sebastian Funk
}