Sequential Bayes (seqB) — seq_bayes • Rnaught

This function implements a sequential Bayesian estimation method of R0 due to Bettencourt and Riberio (PloS One, 2008). See details for important implementation notes.

Usage

seq_bayes(cases, mu, kappa = 20, post = FALSE)

Arguments

cases: Vector of case counts. The vector must only contain non-negative integers, and have at least two positive integers.
mu: Mean of the serial distribution. This must be a positive number. The value should match the case counts in time units. For example, if case counts are weekly and the serial distribution has a mean of seven days, then mu should be set to 1. If case counts are daily and the serial distribution has a mean of seven days, then mu should be set to 7.
kappa: Largest possible value of the uniform prior (defaults to 20). This must be a number greater than or equal to 1. It describes the prior belief on the ranges of R0, and should be set to a higher value if R0 is believed to be larger.
post: Whether to return the posterior distribution of R0 instead of the estimate of R0 (defaults to FALSE). This must be a value identical to TRUE or FALSE.

Value

If post is identical to TRUE, a list containing the following components is returned:

supp - the support of the posterior distribution of R0
pmf - the probability mass function of the posterior distribution of R0

Otherwise, if post is identical to FALSE, only the estimate of R0 is returned. Note that the estimate is equal to sum(supp * pmf) (i.e., the posterior mean).

Details

The method sets a uniform prior distribution on R0 with possible values between 0 and kappa, discretized to a fine grid. The distribution of R0 is then updated sequentially, with one update for each new case count observation. The final estimate of R0 is the mean of the (last) posterior distribution. The prior distribution is the initial belief of the distribution of R0, which is the uninformative uniform distribution with values between 0 and kappa. Users can change the value of kappa only (i.e., the prior distribution cannot be changed from the uniform). As more case counts are observed, the influence of the prior distribution should lessen on the final estimate.

This method is based on an approximation of the SIR model, which is most valid at the beginning of an epidemic. The method assumes that the mean of the serial distribution (sometimes called the serial interval) is known. The final estimate can be quite sensitive to this value, so sensitivity testing is strongly recommended. Users should be careful about units of time (e.g., are counts observed daily or weekly?) when implementing.

Our code has been modified to provide an estimate even if case counts equal to zero are present in some time intervals. This is done by grouping the counts over such periods of time. Without grouping, and in the presence of zero counts, no estimate can be provided.

References

Bettencourt and Riberio (PloS One, 2008)

See also

vignette("seq_bayes_post", package = "Rnaught") for examples of using the posterior distribution.

Examples

# Weekly data.
cases <- c(1, 4, 10, 5, 3, 4, 19, 3, 3, 14, 4)

# Obtain R0 when the serial distribution has a mean of five days.
seq_bayes(cases, mu = 5 / 7)
#> [1] 1.026563

# Obtain R0 when the serial distribution has a mean of three days.
seq_bayes(cases, mu = 3 / 7)
#> [1] 1.015938

# Obtain R0 when the serial distribution has a mean of seven days, and R0 is
# believed to be at most 4.
estimate <- seq_bayes(cases, mu = 1, kappa = 4)

# Same as above, but return the posterior distribution of R0 instead of the
# estimate.
posterior <- seq_bayes(cases, mu = 1, kappa = 4, post = TRUE)

# Display the support and probability mass function of the posterior.
posterior$supp
#>   [1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14
#>  [16] 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29
#>  [31] 0.30 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41 0.42 0.43 0.44
#>  [46] 0.45 0.46 0.47 0.48 0.49 0.50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59
#>  [61] 0.60 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.70 0.71 0.72 0.73 0.74
#>  [76] 0.75 0.76 0.77 0.78 0.79 0.80 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89
#>  [91] 0.90 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04
#> [106] 1.05 1.06 1.07 1.08 1.09 1.10 1.11 1.12 1.13 1.14 1.15 1.16 1.17 1.18 1.19
#> [121] 1.20 1.21 1.22 1.23 1.24 1.25 1.26 1.27 1.28 1.29 1.30 1.31 1.32 1.33 1.34
#> [136] 1.35 1.36 1.37 1.38 1.39 1.40 1.41 1.42 1.43 1.44 1.45 1.46 1.47 1.48 1.49
#> [151] 1.50 1.51 1.52 1.53 1.54 1.55 1.56 1.57 1.58 1.59 1.60 1.61 1.62 1.63 1.64
#> [166] 1.65 1.66 1.67 1.68 1.69 1.70 1.71 1.72 1.73 1.74 1.75 1.76 1.77 1.78 1.79
#> [181] 1.80 1.81 1.82 1.83 1.84 1.85 1.86 1.87 1.88 1.89 1.90 1.91 1.92 1.93 1.94
#> [196] 1.95 1.96 1.97 1.98 1.99 2.00 2.01 2.02 2.03 2.04 2.05 2.06 2.07 2.08 2.09
#> [211] 2.10 2.11 2.12 2.13 2.14 2.15 2.16 2.17 2.18 2.19 2.20 2.21 2.22 2.23 2.24
#> [226] 2.25 2.26 2.27 2.28 2.29 2.30 2.31 2.32 2.33 2.34 2.35 2.36 2.37 2.38 2.39
#> [241] 2.40 2.41 2.42 2.43 2.44 2.45 2.46 2.47 2.48 2.49 2.50 2.51 2.52 2.53 2.54
#> [256] 2.55 2.56 2.57 2.58 2.59 2.60 2.61 2.62 2.63 2.64 2.65 2.66 2.67 2.68 2.69
#> [271] 2.70 2.71 2.72 2.73 2.74 2.75 2.76 2.77 2.78 2.79 2.80 2.81 2.82 2.83 2.84
#> [286] 2.85 2.86 2.87 2.88 2.89 2.90 2.91 2.92 2.93 2.94 2.95 2.96 2.97 2.98 2.99
#> [301] 3.00 3.01 3.02 3.03 3.04 3.05 3.06 3.07 3.08 3.09 3.10 3.11 3.12 3.13 3.14
#> [316] 3.15 3.16 3.17 3.18 3.19 3.20 3.21 3.22 3.23 3.24 3.25 3.26 3.27 3.28 3.29
#> [331] 3.30 3.31 3.32 3.33 3.34 3.35 3.36 3.37 3.38 3.39 3.40 3.41 3.42 3.43 3.44
#> [346] 3.45 3.46 3.47 3.48 3.49 3.50 3.51 3.52 3.53 3.54 3.55 3.56 3.57 3.58 3.59
#> [361] 3.60 3.61 3.62 3.63 3.64 3.65 3.66 3.67 3.68 3.69 3.70 3.71 3.72 3.73 3.74
#> [376] 3.75 3.76 3.77 3.78 3.79 3.80 3.81 3.82 3.83 3.84 3.85 3.86 3.87 3.88 3.89
#> [391] 3.90 3.91 3.92 3.93 3.94 3.95 3.96 3.97 3.98 3.99 4.00
posterior$pmf
#>   [1]  4.396081e-14  6.866777e-14  1.069979e-13  1.663113e-13  2.578585e-13
#>   [6]  3.987896e-13  6.151736e-13  9.465244e-13  1.452563e-12  2.223289e-12
#>  [11]  3.393931e-12  5.167074e-12  7.845296e-12  1.187914e-11  1.793742e-11
#>  [16]  2.700983e-11  4.055633e-11  6.072366e-11  9.065827e-11  1.349567e-10
#>  [21]  2.003117e-10  2.964355e-10  4.373744e-10  6.433725e-10  9.435054e-10
#>  [26]  1.379386e-09  2.010358e-09  2.920742e-09  4.229916e-09  6.106257e-09
#>  [31]  8.786364e-09  1.260143e-08  1.801328e-08  2.566338e-08  3.643912e-08
#>  [36]  5.156328e-08  7.271379e-08  1.021837e-07  1.430935e-07  1.996715e-07
#>  [41]  2.776227e-07  3.846102e-07  5.308815e-07  7.300785e-07  1.000278e-06
#>  [46]  1.365319e-06  1.856496e-06  2.514686e-06  3.393017e-06  4.560205e-06
#>  [51]  6.104659e-06  8.139541e-06  1.080892e-05  1.429523e-05  1.882819e-05
#>  [56]  2.469542e-05  3.225496e-05  4.194988e-05  5.432506e-05  7.004650e-05
#>  [61]  8.992291e-05  1.149299e-04  1.462361e-04  1.852322e-04  2.335598e-04
#>  [66]  2.931434e-04  3.662201e-04  4.553696e-04  5.635406e-04  6.940727e-04
#>  [71]  8.507116e-04  1.037615e-03  1.259347e-03  1.520857e-03  1.827440e-03
#>  [76]  2.184679e-03  2.598363e-03  3.074378e-03  3.618573e-03  4.236597e-03
#>  [81]  4.933708e-03  5.714562e-03  6.582973e-03  7.541662e-03  8.591990e-03
#>  [86]  9.733699e-03  1.096466e-02  1.228062e-02  1.367503e-02  1.513890e-02
#>  [91]  1.666067e-02  1.822624e-02  1.981901e-02  2.142007e-02  2.300844e-02
#>  [96]  2.456146e-02  2.605520e-02  2.746508e-02  2.876640e-02  2.993509e-02
#> [101]  3.094837e-02  3.178548e-02  3.242832e-02  3.286214e-02  3.307605e-02
#> [106]  3.306344e-02  3.282231e-02  3.235539e-02  3.167013e-02  3.077854e-02
#> [111]  2.969682e-02  2.844487e-02  2.704571e-02  2.552471e-02  2.390888e-02
#> [116]  2.222595e-02  2.050365e-02  1.876888e-02  1.704700e-02  1.536123e-02
#> [121]  1.373213e-02  1.217724e-02  1.071084e-02  9.343864e-03  8.083917e-03
#> [126]  6.935429e-03  5.899896e-03  4.976201e-03  4.160989e-03  3.449076e-03
#> [131]  2.833858e-03  2.307727e-03  1.862441e-03  1.489475e-03  1.180313e-03
#> [136]  9.266886e-04  7.207803e-04  5.553459e-04  4.238132e-04  3.203273e-04
#> [141]  2.397617e-04  1.777009e-04  1.304008e-04  9.473459e-05  6.812875e-05
#> [146]  4.849552e-05  3.416469e-05  2.381841e-05  1.643092e-05  1.121447e-05
#> [151]  7.572111e-06  5.057422e-06  3.340936e-06  2.182658e-06  1.410047e-06
#> [156]  9.006636e-07  5.687529e-07  3.550312e-07  2.190488e-07  1.335660e-07
#> [161]  8.047851e-08  4.791165e-08  2.817909e-08  1.637135e-08  9.394205e-09
#> [166]  5.323528e-09  2.978849e-09  1.645703e-09  8.975385e-10  4.831664e-10
#> [171]  2.566999e-10  1.345806e-10  6.961600e-11  3.552605e-11  1.788285e-11
#> [176]  8.878079e-12  4.346434e-12  2.098061e-12  9.984180e-13  4.683320e-13
#> [181]  2.165108e-13  9.863385e-14  4.427200e-14  1.957601e-14  8.526007e-15
#> [186]  3.657023e-15  1.544555e-15  6.422519e-16  2.628849e-16  1.059047e-16
#> [191]  4.198412e-17  1.637588e-17  6.283525e-18  2.371425e-18  8.801378e-19
#> [196]  3.211841e-19  1.152248e-19  4.063045e-20  1.407979e-20  4.794057e-21
#> [201]  1.603601e-21  5.268629e-22  1.699923e-22  5.385323e-23  1.674809e-23
#> [206]  5.112219e-24  1.531308e-24  4.500307e-25  1.297374e-25  3.668155e-26
#> [211]  1.016962e-26  2.764082e-27  7.363758e-28  1.922485e-28  4.917601e-29
#> [216]  1.232200e-29  3.023823e-30  7.265876e-31  1.709166e-31  3.935069e-32
#> [221]  8.865403e-33  1.954015e-33  4.212549e-34  8.880821e-35  1.830437e-35
#> [226]  3.687665e-36  7.260114e-37  1.396466e-37  2.623677e-38  4.813725e-39
#> [231]  8.622619e-40  1.507577e-40  2.572153e-41  4.281381e-42  6.950773e-43
#> [236]  1.100362e-43  1.698170e-44  2.554224e-45  3.743317e-46  5.343929e-47
#> [241]  7.429450e-48  1.005609e-48  1.324832e-49  1.698377e-50  2.118019e-51
#> [246]  2.568782e-52  3.029042e-53  3.471688e-54  3.866425e-55  4.182989e-56
#> [251]  4.394865e-57  4.482891e-58  4.438084e-59  4.263117e-60  3.972112e-61
#> [256]  3.588768e-62  3.143145e-63  2.667736e-64  2.193531e-65  1.746742e-66
#> [261]  1.346659e-67  1.004822e-68  7.254056e-70  5.065108e-71  3.419533e-72
#> [266]  2.231354e-73  1.406840e-74  8.567343e-76  5.037567e-77  2.859011e-78
#> [271]  1.565585e-79  8.268889e-81  4.210849e-82  2.066735e-83  9.773085e-85
#> [276]  4.450901e-86  1.951504e-87  8.234368e-89  3.342436e-90  1.304664e-91
#> [281]  4.895156e-93  1.764793e-94  6.110914e-96  2.031551e-97  6.481579e-99
#> [286] 1.983744e-100 5.821848e-102 1.637658e-103 4.413549e-105 1.139116e-106
#> [291] 2.814324e-108 6.652979e-110 1.504185e-111 3.251137e-113 6.714617e-115
#> [296] 1.324527e-116 2.494329e-118 4.482263e-120 7.682227e-122 1.255211e-123
#> [301] 1.954236e-125 2.897728e-127 4.090211e-129 5.493210e-131 7.015875e-133
#> [306] 8.517133e-135 9.822891e-137 1.075711e-138 1.117987e-140 1.102133e-142
#> [311] 1.030047e-144 9.121635e-147 7.649726e-149 6.072101e-151 4.559445e-153
#> [316] 3.236849e-155 2.171329e-157 1.375543e-159 8.224662e-162 4.638796e-164
#> [321] 2.466499e-166 1.235627e-168 5.828612e-171 2.587326e-173 1.080141e-175
#> [326] 4.238233e-178 1.562047e-180 5.404238e-183 1.754003e-185 5.337064e-188
#> [331] 1.521491e-190 4.061129e-193 1.014255e-195 2.368537e-198 5.168362e-201
#> [336] 1.053099e-203 2.002304e-206 3.550040e-209 5.865079e-212 9.022847e-215
#> [341] 1.291610e-217 1.719185e-220 2.126178e-223 2.441418e-226 2.600913e-229
#> [346] 2.568760e-232 2.350200e-235 1.990379e-238 1.559119e-241 1.128742e-244
#> [351] 7.546371e-248 4.655454e-251 2.647985e-254 1.387535e-257 6.692519e-261
#> [356] 2.968866e-264 1.210268e-267 4.529956e-271 1.555446e-274 4.895396e-278
#> [361] 1.410954e-281 3.720881e-285 8.970133e-289 1.975060e-292 3.968208e-296
#> [366] 7.268446e-300 1.212601e-303 1.840840e-307 2.540527e-311 3.184384e-315
#> [371] 3.621946e-319  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
#> [376]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
#> [381]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
#> [386]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
#> [391]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
#> [396]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
#> [401]  0.000000e+00

# Note that the following always holds:
estimate == sum(posterior$supp * posterior$pmf)
#> [1] TRUE