X-Git-Url: https://git.nmode.ca/Rnaught/blobdiff_plain/ff76a13dbabf3a0e8a0be8c82293597762df55cc..12ec5983556b900dfe0abbd8b55b4fb2bc84a050:/reference/seq_bayes.html diff --git a/reference/seq_bayes.html b/reference/seq_bayes.html new file mode 100644 index 0000000..f96bc04 --- /dev/null +++ b/reference/seq_bayes.html @@ -0,0 +1,283 @@ + +Sequential Bayes (seqB) — seq_bayes • Rnaught + Skip to contents + + +
+
+
+ +
+

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.

+
+ +
+

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
+
+
+
+ + +
+ + + + + + +