diff --git a/compare-interventions.md b/compare-interventions.md index 1a428cec..6f47272c 100644 --- a/compare-interventions.md +++ b/compare-interventions.md @@ -86,7 +86,7 @@ The diagram below describes the flow of individuals through the different compar To run the model with no vaccination in place we can *either* create two vaccination objects (one for each dose) using `vaccination()` with the time start, time end and vaccination rate all set to 0, or we can use the `no_vaccination()` function to create a vaccination object for two doses with all values set to 0. -```r +``` r no_vaccination <- no_vaccination(population = uk_population, doses = 2) ``` :::::::::::::::::::::: @@ -99,7 +99,7 @@ We can run the Vacamole model with [default parameter values](https://epiverse-t -```r +``` r output <- model_vacamole_cpp( population = uk_population, vaccination = no_vaccination, @@ -118,7 +118,7 @@ output <- model_vacamole_cpp( 1. Run the model -```r +``` r polymod <- socialmixr::polymod contact_data <- socialmixr::contact_matrix( survey = polymod, @@ -128,15 +128,11 @@ contact_data <- socialmixr::contact_matrix( ) ``` -```{.output} -Using POLYMOD social contact data. To cite this in a publication, use the 'get_citation()' function -``` - -```{.output} +``` output Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option ``` -```r +``` r # prepare contact matrix contact_matrix <- t(contact_data$matrix) @@ -171,19 +167,19 @@ uk_population <- population( ) ``` -```{.error} +``` error Error in population(name = "UK", contact_matrix = contact_matrix, demography_vector = demography_vector, : could not find function "population" ``` -```r +``` r no_vaccination <- no_vaccination(population = uk_population, doses = 2) ``` -```{.error} +``` error Error in no_vaccination(population = uk_population, doses = 2): could not find function "no_vaccination" ``` -```r +``` r # run model output <- model_vacamole_cpp( population = uk_population, @@ -192,14 +188,14 @@ output <- model_vacamole_cpp( ) ``` -```{.error} +``` error Error in model_vacamole_cpp(population = uk_population, vaccination = no_vaccination, : could not find function "model_vacamole_cpp" ``` 2. Plot the number of deaths through time -```r +``` r ggplot(output[output$compartment == "dead", ]) + geom_line( aes(time, value, colour = demography_group), @@ -225,7 +221,7 @@ ggplot(output[output$compartment == "dead", ]) + ) ``` -```{.error} +``` error Error in ggplot(output[output$compartment == "dead", ]): could not find function "ggplot" ``` diff --git a/create-forecast.md b/create-forecast.md index 1354ca1c..7e062ba8 100644 --- a/create-forecast.md +++ b/create-forecast.md @@ -47,10 +47,18 @@ Given case data, we can create estimates of the current and future number of cas The function `epinow()` described in the previous tutorial is a wrapper for the function `estimate_infections()` used to estimate cases by date of infection. Using the same code in the previous tutorial we can extract the short-term forecast using: +``` warning +Warning: `dist_spec()` was deprecated in EpiNow2 1.5.0. +ℹ Please use distribution functions such as `Gamma()` or `Lognormal()` instead. +ℹ The function will become internal only in the next version. +This warning is displayed once every 8 hours. +Call `lifecycle::last_lifecycle_warnings()` to see where this warning was +generated. +``` -```r +``` r reported_cases <- cases[1:90, ] estimates <- epinow( reported_cases = reported_cases, @@ -60,12 +68,30 @@ estimates <- epinow( ) ``` -```{.output} -WARN [2024-07-02 02:08:31] epinow: There were 7 divergent transitions after warmup. See +``` warning +Warning: The `reported_cases` argument of `epinow()` is deprecated as of EpiNow2 1.5.0. +ℹ Please use the `data` argument instead. +ℹ The argument will be removed completely in the next version. +This warning is displayed once every 8 hours. +Call `lifecycle::last_lifecycle_warnings()` to see where this warning was +generated. +``` + +``` output +WARN [2024-07-02 02:25:21] epinow: There were 2000 divergent transitions after warmup. See https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup to find out why this is a problem and how to eliminate them. - -WARN [2024-07-02 02:08:31] epinow: Examine the pairs() plot to diagnose sampling problems +WARN [2024-07-02 02:25:21] epinow: Examine the pairs() plot to diagnose sampling problems - +WARN [2024-07-02 02:25:22] epinow: The largest R-hat is NA, indicating chains have not mixed. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#r-hat - +WARN [2024-07-02 02:25:23] epinow: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#bulk-ess - +WARN [2024-07-02 02:25:26] epinow: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#tail-ess - ``` @@ -79,7 +105,7 @@ We can visualise the estimates of the effective reproduction number and the esti -```r +``` r plot(estimates) ``` @@ -94,7 +120,7 @@ By default, the short-term forecasts are created using the latest estimate of th The reproduction number that is projected into the future can be changed to a less recent estimate based on more data using `rt_opts()`: -```r +``` r rt_opts(future = "estimate") ``` @@ -114,14 +140,14 @@ We will pass another input into `epinow()` called `obs` to define an observation Let's say we believe the COVID-19 outbreak data from the previous tutorial do not include all reported cases. We believe that we only observe 40% of cases. To specify this in the observation model, we must pass a scaling factor with a mean and standard deviation. If we assume that 40% of cases are in the case data (with standard deviation 1%), then we specify the `scale` input to `obs_opts()` as follows: -```r +``` r obs_scale <- list(mean = 0.4, sd = 0.01) ``` To run the inference framework with this observation process, we add `obs = obs_opts(scale = obs_scale)` to the input arguments of `epinow()`: -```r +``` r obs_scale <- list(mean = 0.4, sd = 0.01) reported_cases <- cases[1:90, ] estimates <- epinow( @@ -133,26 +159,35 @@ estimates <- epinow( ) ``` -```{.output} -WARN [2024-07-02 02:15:47] epinow: There were 2 divergent transitions after warmup. See +``` output +WARN [2024-07-02 02:30:27] epinow: There were 2000 divergent transitions after warmup. See https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup to find out why this is a problem and how to eliminate them. - -WARN [2024-07-02 02:15:47] epinow: Examine the pairs() plot to diagnose sampling problems +WARN [2024-07-02 02:30:27] epinow: Examine the pairs() plot to diagnose sampling problems - +WARN [2024-07-02 02:30:28] epinow: The largest R-hat is NA, indicating chains have not mixed. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#r-hat - +WARN [2024-07-02 02:30:29] epinow: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#bulk-ess - +WARN [2024-07-02 02:30:31] epinow: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#tail-ess - ``` -```r +``` r summary(estimates) ``` -```{.output} - measure estimate - -1: New confirmed cases by infection date 17817 (9909 -- 31327) -2: Expected change in daily cases Likely decreasing -3: Effective reproduction no. 0.88 (0.56 -- 1.3) -4: Rate of growth -0.016 (-0.067 -- 0.038) -5: Doubling/halving time (days) -42 (18 -- -10) +``` output + measure estimate + +1: New infections per day 2487 (257 -- 26341) +2: Expected change in daily reports Likely decreasing +3: Effective reproduction no. 0.81 (0.56 -- 1.1) +4: Rate of growth -0.21 (-0.58 -- 0.064) +5: Doubling/halving time (days) -3.3 (11 -- -1.2) ``` @@ -172,7 +207,7 @@ First, we must format our data to have the following columns: + `secondary` : number of secondary observations date, in this example **deaths**. -```r +``` r reported_cases_deaths <- aggregate( cbind(cases_new, deaths_new) ~ date, data = @@ -198,7 +233,7 @@ There are further function inputs to `estimate_secondary()` which can be specifi To find the model fit between cases and deaths : -```r +``` r estimate_cases_to_deaths <- estimate_secondary( reports = reported_cases_deaths[1:30, ], secondary = secondary_opts(type = "incidence"), @@ -209,6 +244,24 @@ estimate_cases_to_deaths <- estimate_secondary( ) ``` +``` warning +Warning: The `reports` argument of `estimate_secondary()` is deprecated as of EpiNow2 +1.5.0. +ℹ Please use the `data` argument instead. +ℹ The argument will be removed completely in the next version. +This warning is displayed once every 8 hours. +Call `lifecycle::last_lifecycle_warnings()` to see where this warning was +generated. +``` + +``` output +WARN [2024-07-02 02:30:42] estimate_secondary (chain: 1): There were 2 divergent transitions after warmup. See +https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup +to find out why this is a problem and how to eliminate them. - +WARN [2024-07-02 02:30:42] estimate_secondary (chain: 1): Examine the pairs() plot to diagnose sampling problems + - +``` + ::::::::::::::::::::::::::::::::::::: callout ### Be cautious of timescale @@ -220,7 +273,7 @@ In the early stages of an outbreak there can be substantial changes in testing a We plot the model fit (shaded ribbons) with the secondary observations (bar plot) and primary observations (dotted line) as follows: -```r +``` r plot(estimate_cases_to_deaths, primary = TRUE) ``` @@ -230,7 +283,7 @@ To use this model fit to forecast deaths, we pass a data frame consisting of the *Note : in this tutorial we are using data where we know the deaths and cases, so we create a data frame by extracting the cases. But in practice, this would be a different data set consisting of cases only.* -```r +``` r cases_to_forecast <- reported_cases_deaths[31:60, c("date", "primary")] colnames(cases_to_forecast) <- c("date", "value") ``` @@ -238,7 +291,7 @@ colnames(cases_to_forecast) <- c("date", "value") To forecast, we use the model fit `estimate_cases_to_deaths`: -```r +``` r deaths_forecast <- forecast_secondary( estimate = estimate_cases_to_deaths, primary = cases_to_forecast @@ -246,8 +299,13 @@ deaths_forecast <- forecast_secondary( plot(deaths_forecast) ``` -```{.warning} -Warning: Removed 30 rows containing missing values (`position_stack()`). +``` warning +Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning +-Inf +Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning +-Inf +Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning +-Inf ``` @@ -298,7 +356,7 @@ As the data consists of date of symptom onset, we only need to specify a delay d We specify the distributions with some uncertainty around the mean and standard deviation of the log normal distribution for the incubation period and the Gamma distribution for the generation time. -```r +``` r ebola_incubation_period <- dist_spec( mean = 2.487, sd = 0.330, mean_sd = 0.5, sd_sd = 0.5, @@ -312,12 +370,20 @@ ebola_generation_time <- dist_spec( ) ``` +``` warning +Warning in new_dist_spec(params, "gamma"): Uncertain gamma distribution +specified in terms of parameters that are not the "natural" parameters of the +distribution (shape, rate). Converting using a crude and very approximate +method that is likely to produce biased results. If possible, it is preferable +to specify the distribution directly in terms of the natural parameters. +``` + As we want to also create a two week forecast, we specify `horizon = 14` to forecast 14 days instead of the default 7 days. -```r +``` r # read data # e.g.: if path to file is data/raw-data/ebola_cases.csv then: ebola_cases <- @@ -325,7 +391,7 @@ ebola_cases <- ``` -```r +``` r # format date column ebola_cases$date <- as.Date(ebola_cases$date) @@ -338,34 +404,37 @@ ebola_estimates <- epinow( ) ``` -```{.output} -WARN [2024-07-02 02:17:51] epinow: There were 24 divergent transitions after warmup. See +``` output +WARN [2024-07-02 02:32:43] epinow: There were 6 divergent transitions after warmup. See https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup to find out why this is a problem and how to eliminate them. - -WARN [2024-07-02 02:17:51] epinow: Examine the pairs() plot to diagnose sampling problems +WARN [2024-07-02 02:32:43] epinow: Examine the pairs() plot to diagnose sampling problems - +WARN [2024-07-02 02:32:46] epinow: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#tail-ess - ``` -```r +``` r summary(ebola_estimates) ``` -```{.output} - measure estimate - -1: New confirmed cases by infection date 100 (46 -- 263) -2: Expected change in daily cases Increasing -3: Effective reproduction no. 1.7 (1.1 -- 2.9) -4: Rate of growth 0.042 (0.0039 -- 0.092) -5: Doubling/halving time (days) 17 (7.6 -- 180) +``` output + measure estimate + +1: New infections per day 107 (46 -- 271) +2: Expected change in daily reports Increasing +3: Effective reproduction no. 1.8 (1 -- 3.2) +4: Rate of growth 0.049 (-0.016 -- 0.12) +5: Doubling/halving time (days) 14 (5.7 -- -42) ``` -The effective reproduction number $R_t$ estimate (on the last date of the data) is 1.7 (1.1 -- 2.9). The exponential growth rate of case numbers is 0.042 (0.0039 -- 0.092). +The effective reproduction number $R_t$ estimate (on the last date of the data) is 1.8 (1 -- 3.2). The exponential growth rate of case numbers is 0.049 (-0.016 -- 0.12). Visualize the estimates: -```r +``` r plot(ebola_estimates) ``` diff --git a/fig/create-forecast-rendered-unnamed-chunk-11-1.png b/fig/create-forecast-rendered-unnamed-chunk-11-1.png index 83e208f3..1959ee54 100644 Binary files a/fig/create-forecast-rendered-unnamed-chunk-11-1.png and b/fig/create-forecast-rendered-unnamed-chunk-11-1.png differ diff --git a/fig/create-forecast-rendered-unnamed-chunk-3-1.png b/fig/create-forecast-rendered-unnamed-chunk-3-1.png index 696f0a16..31265049 100644 Binary files a/fig/create-forecast-rendered-unnamed-chunk-3-1.png and b/fig/create-forecast-rendered-unnamed-chunk-3-1.png differ diff --git a/fig/create-forecast-rendered-unnamed-chunk-9-1.png b/fig/create-forecast-rendered-unnamed-chunk-9-1.png index e59d6b33..a9d24a8c 100644 Binary files a/fig/create-forecast-rendered-unnamed-chunk-9-1.png and b/fig/create-forecast-rendered-unnamed-chunk-9-1.png differ diff --git a/fig/introduction-rendered-unnamed-chunk-5-1.png b/fig/introduction-rendered-unnamed-chunk-5-1.png index 0f6e9c4d..f8e9aaf0 100644 Binary files a/fig/introduction-rendered-unnamed-chunk-5-1.png and b/fig/introduction-rendered-unnamed-chunk-5-1.png differ diff --git a/fig/model-choices-rendered-unnamed-chunk-3-1.png b/fig/model-choices-rendered-unnamed-chunk-3-1.png deleted file mode 100644 index c76d4111..00000000 Binary files a/fig/model-choices-rendered-unnamed-chunk-3-1.png and /dev/null differ diff --git a/fig/model-choices-rendered-unnamed-chunk-4-1.png b/fig/model-choices-rendered-unnamed-chunk-4-1.png deleted file mode 100644 index ce04c708..00000000 Binary files a/fig/model-choices-rendered-unnamed-chunk-4-1.png and /dev/null differ diff --git a/fig/modelling-interventions-rendered-baseline-1.png b/fig/modelling-interventions-rendered-baseline-1.png deleted file mode 100644 index 8bde486a..00000000 Binary files a/fig/modelling-interventions-rendered-baseline-1.png and /dev/null differ diff --git a/fig/modelling-interventions-rendered-plot_masks-1.png b/fig/modelling-interventions-rendered-plot_masks-1.png deleted file mode 100644 index 6669bb0d..00000000 Binary files a/fig/modelling-interventions-rendered-plot_masks-1.png and /dev/null differ diff --git a/fig/modelling-interventions-rendered-plot_vaccinate-1.png b/fig/modelling-interventions-rendered-plot_vaccinate-1.png deleted file mode 100644 index f598e879..00000000 Binary files a/fig/modelling-interventions-rendered-plot_vaccinate-1.png and /dev/null differ diff --git a/fig/quantify-transmissibility-rendered-unnamed-chunk-15-1.png b/fig/quantify-transmissibility-rendered-unnamed-chunk-15-1.png index a3aa907f..8756258d 100644 Binary files a/fig/quantify-transmissibility-rendered-unnamed-chunk-15-1.png and b/fig/quantify-transmissibility-rendered-unnamed-chunk-15-1.png differ diff --git a/fig/quantify-transmissibility-rendered-unnamed-chunk-16-1.png b/fig/quantify-transmissibility-rendered-unnamed-chunk-16-1.png index a358f186..60f2937c 100644 Binary files a/fig/quantify-transmissibility-rendered-unnamed-chunk-16-1.png and b/fig/quantify-transmissibility-rendered-unnamed-chunk-16-1.png differ diff --git a/fig/quantify-transmissibility-rendered-unnamed-chunk-19-1.png b/fig/quantify-transmissibility-rendered-unnamed-chunk-19-1.png index 09186acd..d5f25e1e 100644 Binary files a/fig/quantify-transmissibility-rendered-unnamed-chunk-19-1.png and b/fig/quantify-transmissibility-rendered-unnamed-chunk-19-1.png differ diff --git a/fig/read-delays-rendered-unnamed-chunk-20-1.png b/fig/read-delays-rendered-unnamed-chunk-20-1.png index c0464797..522dda04 100644 Binary files a/fig/read-delays-rendered-unnamed-chunk-20-1.png and b/fig/read-delays-rendered-unnamed-chunk-20-1.png differ diff --git a/fig/simulating-transmission-rendered-plot-1.png b/fig/simulating-transmission-rendered-plot-1.png deleted file mode 100644 index 81e76c24..00000000 Binary files a/fig/simulating-transmission-rendered-plot-1.png and /dev/null differ diff --git a/fig/simulating-transmission-rendered-traj-1.png b/fig/simulating-transmission-rendered-traj-1.png deleted file mode 100644 index 6de308ad..00000000 Binary files a/fig/simulating-transmission-rendered-traj-1.png and /dev/null differ diff --git a/fig/simulating-transmission-rendered-visualise-1.png b/fig/simulating-transmission-rendered-visualise-1.png deleted file mode 100644 index 6de308ad..00000000 Binary files a/fig/simulating-transmission-rendered-visualise-1.png and /dev/null differ diff --git a/introduction.md b/introduction.md index fab36112..f53e4bba 100644 --- a/introduction.md +++ b/introduction.md @@ -68,20 +68,31 @@ Reflect on your experiences. The `{EpiNow2}` package provides a three-step solution to _quantify the transmissibility_. Let's see how to do this with a minimal example. First, load the package: -```r +``` r library(EpiNow2) ``` +``` output + +Attaching package: 'EpiNow2' +``` + +``` output +The following object is masked from 'package:stats': + + Gamma +``` + ### First, get your case data Case incidence data must be stored in a data frame with the observed number of cases per day. We can read an example from the package: -```r +``` r example_confirmed ``` -```{.output} +``` output date confirm 1: 2020-02-22 14 @@ -102,7 +113,7 @@ example_confirmed Not all primary cases have the same probability of generating a secondary case. The onset and cessation of [infectiousness](../learners/reference.md#infectiousness) may occur gradually. For `{EpiNow2}`, we can specify it as a probability `distribution` with `mean`, standard deviation `sd`, and maximum value `max`: -```r +``` r generation_time <- dist_spec( mean = 3.6, sd = 3.1, @@ -111,6 +122,15 @@ generation_time <- dist_spec( ) ``` +``` warning +Warning: `dist_spec()` was deprecated in EpiNow2 1.5.0. +ℹ Please use distribution functions such as `Gamma()` or `Lognormal()` instead. +ℹ The function will become internal only in the next version. +This warning is displayed once every 8 hours. +Call `lifecycle::last_lifecycle_warnings()` to see where this warning was +generated. +``` + ### Let's calculate the reproduction number! In the `epinow()` function we can add: @@ -120,7 +140,7 @@ In the `epinow()` function we can add: - the computation `stan` parameters for this calculation: -```r +``` r epinow_estimates <- epinow( # cases reported_cases = example_confirmed[1:60], @@ -134,11 +154,19 @@ epinow_estimates <- epinow( ) ``` -```{.output} -WARN [2024-07-02 01:12:06] epinow: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. +``` output +WARN [2024-07-02 01:51:26] epinow: There were 940 divergent transitions after warmup. See +https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup +to find out why this is a problem and how to eliminate them. - +WARN [2024-07-02 01:51:26] epinow: Examine the pairs() plot to diagnose sampling problems + - +WARN [2024-07-02 01:51:27] epinow: The largest R-hat is NA, indicating chains have not mixed. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#r-hat - +WARN [2024-07-02 01:51:27] epinow: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#bulk-ess - -WARN [2024-07-02 01:12:06] epinow: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. +WARN [2024-07-02 01:51:28] epinow: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#tail-ess - ``` @@ -146,7 +174,7 @@ https://mc-stan.org/misc/warnings.html#tail-ess - As an output, we get the time-varying (or [effective](../learners/reference.md#effectiverepro)) reproduction number, as well as the cases by date of report and date of infection: -```r +``` r base::plot(epinow_estimates) ``` diff --git a/md5sum.txt b/md5sum.txt index 2671ef78..19efb661 100644 --- a/md5sum.txt +++ b/md5sum.txt @@ -16,4 +16,4 @@ "learners/reference.md" "3c3b0a87a98f91231d47f5a2e24ad27d" "site/built/reference.md" "2024-07-02" "learners/setup.md" "9030a5d0d8e636cb89ee8d48e864469e" "site/built/setup.md" "2024-07-02" "profiles/learner-profiles.md" "31b503c4b5bd1f0960ada730eca4a25e" "site/built/learner-profiles.md" "2024-07-02" -"renv/profiles/lesson-requirements/renv.lock" "56f73fb7f1e0898f2a946b441fe4ccdb" "site/built/renv.lock" "2024-07-02" +"renv/profiles/lesson-requirements/renv.lock" "3ea8ae92d2aac126735753a34373f1c2" "site/built/renv.lock" "2024-07-02" diff --git a/model-choices.md b/model-choices.md index 52f315f6..ff9c6180 100644 --- a/model-choices.md +++ b/model-choices.md @@ -214,7 +214,7 @@ You have been tasked to generate initial trajectories of an Ebola outbreak in Gu ### Code for initial conditions -```r +``` r # set population size population_size <- 14e6 @@ -256,7 +256,7 @@ Adapt the code from the [accounting for uncertainty](../episodes/simulating-tran -```r +``` r output <- model_ebola_r( population = guinea_population, transmissibility = 1.1 / 12, @@ -267,7 +267,13 @@ output <- model_ebola_r( funeral_risk = 0.5, time_end = 100 ) +``` + +``` error +Error in model_ebola_r(population = guinea_population, transmissibility = 1.1/12, : could not find function "model_ebola_r" +``` +``` r ggplot(output[output$compartment == "infectious", ]) + geom_line( aes(time, value), @@ -285,14 +291,16 @@ ggplot(output[output$compartment == "infectious", ]) + ) ``` - +``` error +Error in eval(expr, envir, enclos): object 'output' not found +``` 2. Run model 100 times and plot the mean, upper and lower 95% quantiles of the number of infectious individuals through time We run the model 100 times with the *same* parameter values. -```r +``` r output_samples <- Map( x = seq(100), f = function(x) { @@ -311,9 +319,21 @@ output_samples <- Map( output } ) +``` +``` error +Error in model_ebola_r(population = guinea_population, transmissibility = 1.1/12, : could not find function "model_ebola_r" +``` + +``` r output_samples <- bind_rows(output_samples) # requires the dplyr package +``` +``` error +Error in eval(expr, envir, enclos): object 'output_samples' not found +``` + +``` r ggplot( output_samples[output_samples$compartment == "infectious", ], aes(time, value) @@ -338,7 +358,9 @@ ggplot( ) ``` - +``` error +Error in eval(expr, envir, enclos): object 'output_samples' not found +``` ::::::::::::::::::::::::::: diff --git a/modelling-interventions.md b/modelling-interventions.md index 7d32922b..45387549 100644 --- a/modelling-interventions.md +++ b/modelling-interventions.md @@ -50,7 +50,7 @@ In this tutorial different types of intervention and how they can be modelled ar We will investigate the effect of interventions on a COVID-19 outbreak using an SEIR model (`model_default_cpp()` in the R package `{epidemics}`). We will set $R_0 = 2.7$, latent period or preinfectious period $= 4$ and the infectious_period $= 5.5$ (parameters adapted from [Davies et al. (2020)](https://doi.org/10.1016/S2468-2667(20)30133-X)). We load a contact matrix with age bins 0-18, 18-65, 65 years and older using `{socialmixr}` and assume that one in every 1 million in each age group is infectious at the start of the epidemic. -```r +``` r polymod <- socialmixr::polymod contact_data <- socialmixr::contact_matrix( polymod, @@ -95,18 +95,18 @@ The first NPI we will consider is the effect of school closures on reducing the To include an intervention in our model we must create an `intervention` object. The inputs are the name of the intervention (`name`), the type of intervention (`contacts` or `rate`), the start time (`time_begin`), the end time (`time_end`) and the reduction (`reduction`). The values of the reduction matrix are specified in the same order as the age groups in the contact matrix. -```r +``` r rownames(contact_matrix) ``` -```{.output} +``` output [1] "[0,15)" "[15,65)" "65+" ``` Therefore, we specify ` reduction = matrix(c(0.5, 0.01, 0.01))`. We assume that the school closures start on day 50 and are in place for a further 100 days. Therefore our intervention object is : -```r +``` r close_schools <- intervention( name = "School closure", type = "contacts", @@ -122,7 +122,7 @@ close_schools <- intervention( In `epidemics`, the contact matrix is scaled down by proportions for the period in which the intervention is in place. To understand how the reduction is calculated within the model functions, consider a contact matrix for two age groups with equal number of contacts: -```{.output} +``` output [,1] [,2] [1,] 1 1 [2,] 1 1 @@ -131,7 +131,7 @@ In `epidemics`, the contact matrix is scaled down by proportions for the period If the reduction is 50% in group 1 and 10% in group 2, the contact matrix during the intervention will be: -```{.output} +``` output [,1] [,2] [1,] 0.25 0.45 [2,] 0.45 0.81 @@ -144,7 +144,7 @@ The contacts within group 1 are reduced by 50% twice to accommodate for a 50% re Using transmissibility $= 2.7/5.5$ (remember that [transmissibility = $R_0$/ infectious period](../episodes/simulating-transmission.md#the-basic-reproduction-number-r_0)), infectiousness rate $1/= 4$ and the recovery rate $= 1/5.5$ we run the model with` intervention = list(contacts = close_schools)` as follows : -```r +``` r output_school <- model_default_cpp( population = uk_population, transmissibility = 2.7 / 5.5, @@ -155,11 +155,15 @@ output_school <- model_default_cpp( ) ``` +``` error +Error in model_default_cpp(population = uk_population, transmissibility = 2.7/5.5, : could not find function "model_default_cpp" +``` + To be able to see the effect of our intervention, we also run the model where there is no intervention, combine the two outputs into one data frame and then plot the output. Here we plot the total number of infectious individuals in all age groups using `ggplot2::stat_summary()`: -```r +``` r # run baseline simulation with no intervention output_baseline <- model_default_cpp( population = uk_population, @@ -168,12 +172,38 @@ output_baseline <- model_default_cpp( recovery_rate = 1.0 / 5.5, time_end = 300, increment = 1.0 ) +``` +``` error +Error in model_default_cpp(population = uk_population, transmissibility = 2.7/5.5, : could not find function "model_default_cpp" +``` + +``` r # create intervention_type column for plotting output_school$intervention_type <- "school closure" +``` + +``` error +Error: object 'output_school' not found +``` + +``` r output_baseline$intervention_type <- "baseline" +``` + +``` error +Error: object 'output_baseline' not found +``` + +``` r output <- rbind(output_school, output_baseline) +``` + +``` error +Error in eval(expr, envir, enclos): object 'output_school' not found +``` +``` r ggplot(data = output[output$compartment == "infectious", ]) + aes( x = time, @@ -214,7 +244,9 @@ ggplot(data = output[output$compartment == "infectious", ]) + ) ``` - +``` error +Error in eval(expr, envir, enclos): object 'output' not found +``` We see that with the intervention in place, the infection still spreads through the population, though the peak number of infectious individuals is smaller than the baseline with no intervention in place (solid line). @@ -228,7 +260,7 @@ We expect that mask wearing will reduce an individual's infectiousness. As we ar We create an intervention object with `type = rate` and `reduction = 0.161`. Using parameters adapted from [Li et al. 2020](https://doi.org/10.1371/journal.pone.0237691) we have proportion wearing masks = coverage $\times$ availability = $0.54 \times 0.525 = 0.2835$, proportion reduction in transmissibility = $0.575$. Therefore, $\theta = 0.2835 \times 0.575 = 0.163$. We assume that the mask wearing mandate starts at day 40 and is in place for 200 days. -```r +``` r mask_mandate <- intervention( name = "mask mandate", type = "rate", @@ -241,7 +273,7 @@ mask_mandate <- intervention( To implement this intervention on the parameter $\beta$, we specify `intervention = list(beta = mask_mandate)`. -```r +``` r output_masks <- model_default_cpp( population = uk_population, transmissibility = 2.7 / 5.5, @@ -252,14 +284,38 @@ output_masks <- model_default_cpp( ) ``` +``` error +Error in model_default_cpp(population = uk_population, transmissibility = 2.7/5.5, : could not find function "model_default_cpp" +``` + -```r +``` r # create intervention_type column for plotting output_masks$intervention_type <- "mask mandate" +``` + +``` error +Error: object 'output_masks' not found +``` + +``` r output_baseline$intervention_type <- "baseline" +``` + +``` error +Error: object 'output_baseline' not found +``` + +``` r output <- rbind(output_masks, output_baseline) +``` +``` error +Error in eval(expr, envir, enclos): object 'output_masks' not found +``` + +``` r ggplot(data = output[output$compartment == "infectious", ]) + aes( x = time, @@ -300,7 +356,9 @@ ggplot(data = output[output$compartment == "infectious", ]) + ) ``` - +``` error +Error in eval(expr, envir, enclos): object 'output' not found +``` ::::::::::::::::::::::::::::::::::::: callout ### Intervention types @@ -340,7 +398,7 @@ To explore the effect of vaccination we need to create a vaccination object to p Here we will assume all age groups are vaccinated at the same rate 0.01 and that the vaccination program starts on day 40 and is in place for 150 days. -```r +``` r # prepare a vaccination object vaccinate <- vaccination( name = "vaccinate all", @@ -353,7 +411,7 @@ vaccinate <- vaccination( We pass our vaccination object using `vaccination = vaccinate`: -```r +``` r output_vaccinate <- model_default_cpp( population = uk_population, transmissibility = 2.7 / 5.5, @@ -364,6 +422,10 @@ output_vaccinate <- model_default_cpp( ) ``` +``` error +Error in model_default_cpp(population = uk_population, transmissibility = 2.7/5.5, : could not find function "model_default_cpp" +``` + ::::::::::::::::::::::::::::::::::::: challenge @@ -377,11 +439,24 @@ Plot the three interventions vaccination, school closure and mask mandate and th ## Output -```r +``` r # create intervention_type column for plotting output_vaccinate$intervention_type <- "vaccination" +``` + +``` error +Error: object 'output_vaccinate' not found +``` + +``` r output <- rbind(output_school, output_masks, output_vaccinate, output_baseline) +``` +``` error +Error in eval(expr, envir, enclos): object 'output_school' not found +``` + +``` r ggplot(data = output[output$compartment == "infectious", ]) + aes( x = time, @@ -406,7 +481,9 @@ ggplot(data = output[output$compartment == "infectious", ]) + ) ``` - +``` error +Error in eval(expr, envir, enclos): object 'output' not found +``` From the plot we see that the peak number of total number of infectious individuals when vaccination is in place is much lower compared to school closures and mask wearing interventions. diff --git a/quantify-transmissibility.md b/quantify-transmissibility.md index ed8c29b7..fe5c0c23 100644 --- a/quantify-transmissibility.md +++ b/quantify-transmissibility.md @@ -68,7 +68,7 @@ In Bayesian inference, we use prior knowledge (prior distributions) with data (i The first step is to load the `{EpiNow2}` package : -```r +``` r library(EpiNow2) ``` @@ -85,11 +85,11 @@ This tutorial illustrates the usage of `epinow()` to estimate the time-varying r To illustrate the functions of `EpiNow2` we will use outbreak data of the start of the COVID-19 pandemic from the United Kingdom. The data are available in the R package `{incidence2}`. -```r +``` r head(incidence2::covidregionaldataUK) ``` -```{.output} +``` output date region region_code cases_new cases_total deaths_new 1 2020-01-30 East Midlands E12000004 NA NA NA 2 2020-01-30 East of England E12000006 NA NA NA @@ -119,7 +119,7 @@ To use the data, we must format the data to have two columns: + `confirm` : number of confirmed cases on that date. -```r +``` r cases <- aggregate( cases_new ~ date, data = incidence2::covidregionaldataUK[, c("date", "cases_new")], @@ -169,17 +169,34 @@ The distribution of incubation period can usually be obtained from the literatur We specify a (fixed) gamma distribution with mean $\mu = 4$ and standard deviation $\sigma^2= 2$ (shape = $4$, scale = $1$) using the function `dist_spec()` as follows: -```r +``` r incubation_period_fixed <- dist_spec( mean = 4, sd = 2, max = 20, distribution = "gamma" ) -incubation_period_fixed ``` -```{.output} +``` warning +Warning: `dist_spec()` was deprecated in EpiNow2 1.5.0. +ℹ Please use distribution functions such as `Gamma()` or `Lognormal()` instead. +ℹ The function will become internal only in the next version. +This warning is displayed once every 8 hours. +Call `lifecycle::last_lifecycle_warnings()` to see where this warning was +generated. +``` + +``` r +incubation_period_fixed +``` - Fixed distribution with PMF [0.019 0.12 0.21 0.21 0.17 0.11 0.069 0.039 0.021 0.011 0.0054 0.0026 0.0012 0.00058 0.00026 0.00012 5.3e-05 2.3e-05 1e-05 4.3e-06] +``` output +- gamma distribution (max: 20): + shape: + - fixed value: + 4 + rate: + - fixed value: + 1 ``` The argument `max` is the maximum value the distribution can take, in this example 20 days. @@ -210,18 +227,40 @@ $$\mbox{Normal}(\sigma^2,\sigma_{\sigma^2}^2).$$ We specify this using `dist_spec` with the additional arguments `mean_sd` ($\sigma_{\mu}^2$) and `sd_sd` ($\sigma_{\sigma^2}^2$). -```r +``` r incubation_period_variable <- dist_spec( mean = 4, sd = 2, mean_sd = 0.5, sd_sd = 0.5, max = 20, distribution = "gamma" ) -incubation_period_variable ``` -```{.output} +``` warning +Warning in new_dist_spec(params, "gamma"): Uncertain gamma distribution +specified in terms of parameters that are not the "natural" parameters of the +distribution (shape, rate). Converting using a crude and very approximate +method that is likely to produce biased results. If possible, it is preferable +to specify the distribution directly in terms of the natural parameters. +``` + +``` r +incubation_period_variable +``` - Uncertain gamma distribution with (untruncated) mean 4 (SD 0.5) and SD 2 (SD 0.5) +``` output +- gamma distribution (max: 20): + shape: + - normal distribution: + mean: + 4 + sd: + 0.61 + rate: + - normal distribution: + mean: + 1 + sd: + 0.31 ``` @@ -235,7 +274,7 @@ When specifying a distribution, it is useful to visualise the probability densit If we want to assume that the mean reporting delay is 2 days (with a standard deviation of 1 day), the log normal distribution will look like: -```r +``` r log_mean <- convert_to_logmean(2, 1) log_sd <- convert_to_logsd(2, 1) x <- seq(from = 0, to = 10, length = 1000) @@ -254,7 +293,7 @@ ggplot(df) + Using the mean and standard deviation for the log normal distribution, we can specify a fixed or variable distribution using `dist_spec()` as before: -```r +``` r reporting_delay_variable <- dist_spec( mean = log_mean, sd = log_sd, mean_sd = 0.5, sd_sd = 0.5, @@ -265,7 +304,7 @@ reporting_delay_variable <- dist_spec( If data is available on the time between symptom onset and reporting, we can use the function `estimate_delay()` to estimate a log normal distribution from a vector of delays. The code below illustrates how to use `estimate_delay()` with synthetic delay data. -```r +``` r delay_data <- rlnorm(500, log(5), 1) # synthetic delay data reporting_delay <- estimate_delay( delay_data, @@ -281,7 +320,7 @@ We also must specify a distribution for the generation time. Here we will use a -```r +``` r generation_time_variable <- dist_spec( mean = 3.6, sd = 3.1, mean_sd = 0.5, sd_sd = 0.5, @@ -298,7 +337,7 @@ There are numerous other inputs that can be passed to `epinow()`, see `EpiNow2:: One optional input is to specify a log normal prior for the effective reproduction number $R_t$ at the start of the outbreak. We specify a mean and standard deviation as arguments of `prior` within `rt_opts()`: -```r +``` r rt_log_mean <- convert_to_logmean(2, 1) rt_log_sd <- convert_to_logsd(2, 1) rt <- rt_opts(prior = list(mean = rt_log_mean, sd = rt_log_sd)) @@ -312,7 +351,7 @@ The Bayesian inference is performed using MCMC methods with the program [Stan](h To reduce computation time, we can run chains in parallel. To do this, we must set the number of cores to be used. By default, 4 MCMC chains are run (see `stan_opts()$chains`), so we can set an equal number of cores to be used in parallel as follows: -```r +``` r withr::local_options(list(mc.cores = 4)) ``` @@ -325,7 +364,7 @@ To find the maximum number of available cores on your machine, use `parallel::de *Note : in the code below fixed distributions are used instead of variable. This is to speed up computation time. It is generally recommended to use variable distributions that account for additional uncertainty.* -```r +``` r reported_cases <- cases[1:90, ] estimates <- epinow( reported_cases = reported_cases, @@ -335,12 +374,30 @@ estimates <- epinow( ) ``` -```{.output} -WARN [2024-07-02 01:18:38] epinow: There were 1 divergent transitions after warmup. See +``` warning +Warning: The `reported_cases` argument of `epinow()` is deprecated as of EpiNow2 1.5.0. +ℹ Please use the `data` argument instead. +ℹ The argument will be removed completely in the next version. +This warning is displayed once every 8 hours. +Call `lifecycle::last_lifecycle_warnings()` to see where this warning was +generated. +``` + +``` output +WARN [2024-07-02 01:54:21] epinow: There were 1994 divergent transitions after warmup. See https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup to find out why this is a problem and how to eliminate them. - -WARN [2024-07-02 01:18:38] epinow: Examine the pairs() plot to diagnose sampling problems +WARN [2024-07-02 01:54:21] epinow: Examine the pairs() plot to diagnose sampling problems - +WARN [2024-07-02 01:54:22] epinow: The largest R-hat is NA, indicating chains have not mixed. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#r-hat - +WARN [2024-07-02 01:54:24] epinow: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#bulk-ess - +WARN [2024-07-02 01:54:26] epinow: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable. +Running the chains for more iterations may help. See +https://mc-stan.org/misc/warnings.html#tail-ess - ``` ### Results @@ -348,7 +405,7 @@ WARN [2024-07-02 01:18:38] epinow: Examine the pairs() plot to diagnose sampling We can extract and visualise estimates of the effective reproduction number through time: -```r +``` r estimates$plots$R ``` @@ -358,7 +415,7 @@ The uncertainty in the estimates increases through time. This is because estimat We can also visualise the growth rate estimate through time: -```r +``` r estimates$plots$growth_rate ``` @@ -367,29 +424,29 @@ estimates$plots$growth_rate To extract a summary of the key transmission metrics at the *latest date* in the data: -```r +``` r summary(estimates) ``` -```{.output} - measure estimate - -1: New confirmed cases by infection date 7228 (4056 -- 12526) -2: Expected change in daily cases Likely decreasing -3: Effective reproduction no. 0.89 (0.57 -- 1.3) -4: Rate of growth -0.015 (-0.065 -- 0.038) -5: Doubling/halving time (days) -46 (18 -- -11) +``` output + measure estimate + +1: New infections per day 5646 (590 -- 41311) +2: Expected change in daily reports Likely decreasing +3: Effective reproduction no. 0.91 (0.69 -- 1.2) +4: Rate of growth -0.078 (-0.23 -- 0.14) +5: Doubling/halving time (days) -8.8 (5.1 -- -3.1) ``` As these estimates are based on partial data, they have a wide uncertainty interval. -+ From the summary of our analysis we see that the expected change in daily cases is Likely decreasing with the estimated new confirmed cases 7228 (4056 -- 12526). ++ From the summary of our analysis we see that the expected change in daily cases is with the estimated new confirmed cases . -+ The effective reproduction number $R_t$ estimate (on the last date of the data) is 0.89 (0.57 -- 1.3). ++ The effective reproduction number $R_t$ estimate (on the last date of the data) is 0.91 (0.69 -- 1.2). -+ The exponential growth rate of case numbers is -0.015 (-0.065 -- 0.038). ++ The exponential growth rate of case numbers is -0.078 (-0.23 -- 0.14). -+ The doubling time (the time taken for case numbers to double) is -46 (18 -- -11). ++ The doubling time (the time taken for case numbers to double) is -8.8 (5.1 -- -3.1). ::::::::::::::::::::::::::::::::::::: callout ### `Expected change in daily cases` @@ -420,7 +477,7 @@ The outbreak data of the start of the COVID-19 pandemic from the United Kingdom + `confirm` : number of confirmed cases for a region on a given date. -```r +``` r regional_cases <- incidence2::covidregionaldataUK[, c("date", "cases_new", "region")] colnames(regional_cases) <- c("date", "confirm", "region") @@ -432,7 +489,7 @@ regional_cases <- regional_cases[which(regional_cases$date %in% dates), ] head(regional_cases) ``` -```{.output} +``` output date confirm region 1 2020-01-30 NA East Midlands 2 2020-01-30 NA East of England @@ -445,7 +502,7 @@ head(regional_cases) To find regional estimates, we use the same inputs as `epinow()` to the function `regional_epinow()`: -```r +``` r estimates_regional <- regional_epinow( reported_cases = regional_cases, generation_time = generation_time_opts(generation_time_fixed), @@ -454,73 +511,74 @@ estimates_regional <- regional_epinow( ) ``` -```{.output} -INFO [2024-07-02 01:18:43] Producing following optional outputs: regions, summary, samples, plots, latest -INFO [2024-07-02 01:18:43] Reporting estimates using data up to: 2020-04-28 -INFO [2024-07-02 01:18:43] No target directory specified so returning output -INFO [2024-07-02 01:18:43] Producing estimates for: East Midlands, East of England, England, London, North East, North West, Northern Ireland, Scotland, South East, South West, Wales, West Midlands, Yorkshire and The Humber -INFO [2024-07-02 01:18:43] Regions excluded: none -INFO [2024-07-02 02:02:28] Completed regional estimates -INFO [2024-07-02 02:02:28] Regions with estimates: 13 -INFO [2024-07-02 02:02:28] Regions with runtime errors: 0 -INFO [2024-07-02 02:02:28] Producing summary -INFO [2024-07-02 02:02:28] No summary directory specified so returning summary output -INFO [2024-07-02 02:02:29] No target directory specified so returning timings +``` warning +Warning: The `reported_cases` argument of `regional_epinow()` is deprecated as of +EpiNow2 1.5.0. +ℹ Please use the `data` argument instead. +ℹ The argument will be removed completely in the next version. +This warning is displayed once every 8 hours. +Call `lifecycle::last_lifecycle_warnings()` to see where this warning was +generated. ``` -```r +``` output +INFO [2024-07-02 01:54:29] Producing following optional outputs: regions, summary, samples, plots, latest +INFO [2024-07-02 01:54:29] Reporting estimates using data up to: 2020-04-28 +INFO [2024-07-02 01:54:29] No target directory specified so returning output +INFO [2024-07-02 01:54:29] Producing estimates for: East Midlands, East of England, England, London, North East, North West, Northern Ireland, Scotland, South East, South West, Wales, West Midlands, Yorkshire and The Humber +INFO [2024-07-02 01:54:29] Regions excluded: none +Error in FUN(X[[i]], ...) : + Stan does not support NA (in shifted_cases) in data +Error in FUN(X[[i]], ...) : + Stan does not support NA (in shifted_cases) in data +Error in FUN(X[[i]], ...) : + Stan does not support NA (in shifted_cases) in data +Error in FUN(X[[i]], ...) : + Stan does not support NA (in shifted_cases) in data +Error in FUN(X[[i]], ...) : + Stan does not support NA (in shifted_cases) in data +Error in FUN(X[[i]], ...) : + Stan does not support NA (in shifted_cases) in data +Error in FUN(X[[i]], ...) : + Stan does not support NA (in shifted_cases) in data +Error in FUN(X[[i]], ...) : + Stan does not support NA (in shifted_cases) in data +Error in FUN(X[[i]], ...) : + Stan does not support NA (in shifted_cases) in data +INFO [2024-07-02 02:21:19] Completed regional estimates +INFO [2024-07-02 02:21:19] Regions with estimates: 4 +INFO [2024-07-02 02:21:19] Regions with runtime errors: 9 +INFO [2024-07-02 02:21:19] Producing summary +INFO [2024-07-02 02:21:19] No summary directory specified so returning summary output +INFO [2024-07-02 02:21:19] No target directory specified so returning timings +``` + +``` r estimates_regional$summary$summarised_results$table ``` -```{.output} - Region New confirmed cases by infection date - - 1: East Midlands 347 (212 -- 549) - 2: East of England 534 (328 -- 853) - 3: England 3597 (2177 -- 5608) - 4: London 297 (195 -- 452) - 5: North East 253 (145 -- 436) - 6: North West 551 (318 -- 883) - 7: Northern Ireland 42 (23 -- 82) - 8: Scotland 284 (153 -- 531) - 9: South East 595 (354 -- 1004) -10: South West 414 (292 -- 592) -11: Wales 94 (66 -- 137) -12: West Midlands 272 (142 -- 495) -13: Yorkshire and The Humber 483 (272 -- 792) - Expected change in daily cases Effective reproduction no. - - 1: Likely increasing 1.2 (0.85 -- 1.6) - 2: Likely increasing 1.2 (0.82 -- 1.6) - 3: Likely decreasing 0.92 (0.64 -- 1.3) - 4: Likely decreasing 0.79 (0.57 -- 1.1) - 5: Likely decreasing 0.91 (0.6 -- 1.3) - 6: Likely decreasing 0.86 (0.57 -- 1.2) - 7: Likely decreasing 0.63 (0.38 -- 1) - 8: Likely decreasing 0.91 (0.58 -- 1.4) - 9: Stable 1 (0.67 -- 1.4) -10: Increasing 1.4 (1.1 -- 1.7) -11: Decreasing 0.57 (0.42 -- 0.75) -12: Likely decreasing 0.71 (0.42 -- 1.1) -13: Stable 1 (0.68 -- 1.4) - Rate of growth Doubling/halving time (days) - - 1: 0.024 (-0.02 -- 0.068) 29 (10 -- -34) - 2: 0.022 (-0.024 -- 0.065) 32 (11 -- -28) - 3: -0.011 (-0.052 -- 0.032) -63 (22 -- -13) - 4: -0.029 (-0.064 -- 0.0091) -24 (76 -- -11) - 5: -0.012 (-0.059 -- 0.037) -60 (19 -- -12) - 6: -0.018 (-0.065 -- 0.024) -38 (29 -- -11) - 7: -0.053 (-0.1 -- 0.0048) -13 (150 -- -6.9) - 8: -0.012 (-0.063 -- 0.042) -56 (16 -- -11) - 9: -0.00047 (-0.047 -- 0.051) -1500 (14 -- -15) -10: 0.046 (0.012 -- 0.083) 15 (8.4 -- 60) -11: -0.065 (-0.092 -- -0.034) -11 (-20 -- -7.5) -12: -0.041 (-0.092 -- 0.014) -17 (51 -- -7.5) -13: 0.0037 (-0.046 -- 0.052) 190 (13 -- -15) -``` - -```r +``` output + Region New infections per day + +1: England 1147 (8 -- 13672) +2: London 70 (9 -- 695) +3: South East 794 (125 -- 10455) +4: Yorkshire and The Humber 417 (3 -- 18022) + Expected change in daily reports Effective reproduction no. + +1: Stable 0.98 (0.35 -- 1.2) +2: Likely decreasing 0.9 (0.68 -- 1.1) +3: Stable 1 (0.91 -- 1.2) +4: Stable 0.99 (0.56 -- 1.5) + Rate of growth Doubling/halving time (days) + +1: -0.017 (-0.37 -- 0.17) -40 (4.1 -- -1.9) +2: -0.11 (-0.39 -- 0.055) -6.4 (13 -- -1.8) +3: -0.0021 (-0.092 -- 0.22) -330 (3.1 -- -7.5) +4: -0.0068 (-0.57 -- 0.4) -100 (1.7 -- -1.2) +``` + +``` r estimates_regional$summary$plots$R ``` diff --git a/read-delays.md b/read-delays.md index a2223974..a7c4b50d 100644 --- a/read-delays.md +++ b/read-delays.md @@ -84,7 +84,7 @@ generation_time <- dist_spec( Let's explore how we can access this and other time delays using `{epiparameter}`. We'll use the pipe `%>%` to connect some of their functions, so let's also call to the `{tidyverse}` package: -```r +``` r library(epiparameter) library(EpiNow2) library(tidyverse) @@ -186,25 +186,59 @@ First, let's assume that the data set `example_confirmed` has COVID-19 observed Let's start by looking at how many parameters we have in the epidemiological distributions database (`epidist_db`) for the `disease` named `covid`-19: -```r +``` r epiparameter::epidist_db( disease = "covid" ) ``` -```{.output} +``` output Returning 27 results that match the criteria (22 are parameterised). Use subset to filter by entry variables or single_epidist to return a single entry. -To retrieve the short citation for each use the 'get_citation' function +To retrieve the citation for each use the 'get_citation' function ``` -```{.output} -List of objects - Number of entries in library: 27 - Number of studies in library: 10 - Number of diseases: 1 - Number of delay distributions: 27 - Number of offspring distributions: 0 +``` output +# List of 27 objects +Number of diseases: 1 +❯ COVID-19 +Number of epi distributions: 5 +❯ hospitalisation to death ❯ incubation period ❯ onset to death ❯ onset to hospitalisation ❯ serial interval +[[1]] +Disease: COVID-19 +Pathogen: SARS-CoV-2 +Epi Distribution: incubation period +Study: Men K, Wang X, Yihao L, Zhang G, Hu J, Gao Y, Han H (2020). "Estimate +the incubation period of coronavirus 2019 (COVID-19)." _medRxiv_. +doi:10.1101/2020.02.24.20027474 +. +Parameters: + +[[2]] +Disease: COVID-19 +Pathogen: SARS-CoV-2 +Epi Distribution: incubation period +Study: Rai B, Shukla A, Dwivedi L (2022). "Incubation period for COVID-19: a +systematic review and meta-analysis." _Zeitschrift fur +Gesundheitswissenschaften_. doi:10.1007/s10389-021-01478-1 +. +Parameters: + +[[3]] +Disease: COVID-19 +Pathogen: SARS-CoV-2 +Epi Distribution: incubation period +Study: Alene M, Yismaw L, Assemie M, Ketema D, Gietaneh W, Birhan T (2021). +"Serial interval and incubation period of COVID-19: a systematic review +and meta-analysis." _BMC Infectious Diseases_. +doi:10.1186/s12879-021-05950-x +. +Parameters: + +# ℹ 24 more elements +# ℹ Use `print(n = ...)` to see more elements. +# ℹ Use `parameter_tbl()` to see a summary table of the parameters. +# ℹ Explore database online at: https://epiverse-trace.github.io/epiparameter/dev/articles/database.html ``` From the `{epiparameter}` package, we can use the `epidist_db()` function to ask for any `disease` and also for a specific epidemiological distribution (`epi_dist`). @@ -212,19 +246,19 @@ From the `{epiparameter}` package, we can use the `epidist_db()` function to ask Let's ask now how many parameters we have in the epidemiological distributions database (`epidist_db`) with the generation time using the string `generation`: -```r +``` r epiparameter::epidist_db( epi_dist = "generation" ) ``` -```{.output} +``` output Returning 1 results that match the criteria (1 are parameterised). Use subset to filter by entry variables or single_epidist to return a single entry. -To retrieve the short citation for each use the 'get_citation' function +To retrieve the citation for each use the 'get_citation' function ``` -```{.output} +``` output Disease: Influenza Pathogen: Influenza-A-H1N1 Epi Distribution: generation time @@ -242,20 +276,25 @@ Parameters: Currently, in the library of epidemiological parameters, we have one `generation` time entry for Influenza. Considering the abovementioned considerations, we can look at the `serial` intervals for `COVID`-19. -```r +``` r epiparameter::epidist_db( disease = "COVID", epi_dist = "serial" ) ``` -```{.output} +``` output Returning 4 results that match the criteria (3 are parameterised). Use subset to filter by entry variables or single_epidist to return a single entry. -To retrieve the short citation for each use the 'get_citation' function +To retrieve the citation for each use the 'get_citation' function ``` -```{.output} +``` output +# List of 4 objects +Number of diseases: 1 +❯ COVID-19 +Number of epi distributions: 1 +❯ serial interval [[1]] Disease: COVID-19 Pathogen: SARS-CoV-2 @@ -307,8 +346,8 @@ Parameters: mean: 4.600 sd: 4.400 -attr(,"class") -[1] "multi_epidist" +# ℹ Use `parameter_tbl()` to see a summary table of the parameters. +# ℹ Explore database online at: https://epiverse-trace.github.io/epiparameter/dev/articles/database.html ``` ::::::::::::::::: callout @@ -322,7 +361,7 @@ attr(,"class") We get more than one epidemiological delay. To summarize this view and get the column names from the underlying parameter dataset, we can add the `epiparameter::list_distributions()` function to the previous code using the pipe `%>%`: -```r +``` r epiparameter::epidist_db( disease = "covid", epi_dist = "serial" @@ -330,18 +369,8 @@ epiparameter::epidist_db( epiparameter::list_distributions() ``` -```{.output} -Returning 4 results that match the criteria (3 are parameterised). -Use subset to filter by entry variables or single_epidist to return a single entry. -To retrieve the short citation for each use the 'get_citation' function -``` - -```{.output} - disease epi_distribution prob_distribution author year -1 COVID-19 serial interval Muluneh .... 2021 -2 COVID-19 serial interval lnorm Hiroshi .... 2020 -3 COVID-19 serial interval weibull Hiroshi .... 2020 -4 COVID-19 serial interval norm Lin Yang.... 2020 +``` error +Error: 'list_distributions' is not an exported object from 'namespace:epiparameter' ``` ::::::::::::::::::::::::::::::::: challenge @@ -371,7 +400,7 @@ The `{epiparameter}` combo of `epidist_db()` plus `list_distributions()` list al ::::::::::::::::: solution -```r +``` r # 16 delays distributions epiparameter::epidist_db( disease = "ebola" @@ -396,7 +425,7 @@ Now, from the output of `epiparameter::epidist_db()`, What is an [offspring dist The `epiparameter::epidist_db()` function works as a filtering or subset function. Let's use the `author` argument to filter `Hiroshi Nishiura` parameters: -```r +``` r epiparameter::epidist_db( disease = "covid", epi_dist = "serial", @@ -405,22 +434,14 @@ epiparameter::epidist_db( epiparameter::list_distributions() ``` -```{.output} -Returning 2 results that match the criteria (2 are parameterised). -Use subset to filter by entry variables or single_epidist to return a single entry. -To retrieve the short citation for each use the 'get_citation' function -``` - -```{.output} - disease epi_distribution prob_distribution author year -1 COVID-19 serial interval lnorm Hiroshi .... 2020 -2 COVID-19 serial interval weibull Hiroshi .... 2020 +``` error +Error: 'list_distributions' is not an exported object from 'namespace:epiparameter' ``` We still get more than one epidemiological parameter. We can set the `single_epidist` argument to `TRUE` to only one: -```r +``` r epiparameter::epidist_db( disease = "covid", epi_dist = "serial", @@ -429,15 +450,15 @@ epiparameter::epidist_db( ) ``` -```{.output} +``` output Using Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel coronavirus (COVID-19) infections." _International Journal of Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060 .. -To retrieve the short citation use the 'get_citation' function +To retrieve the citation use the 'get_citation' function ``` -```{.output} +``` output Disease: COVID-19 Pathogen: SARS-CoV-2 Epi Distribution: serial interval @@ -470,7 +491,7 @@ Now, we have an epidemiological parameter we can reuse! We can replace the numbe Let's assign this `` class object to the `covid_serialint` object. -```r +``` r covid_serialint <- epiparameter::epidist_db( disease = "covid", @@ -480,19 +501,19 @@ covid_serialint <- ) ``` -```{.output} +``` output Using Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel coronavirus (COVID-19) infections." _International Journal of Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060 .. -To retrieve the short citation use the 'get_citation' function +To retrieve the citation use the 'get_citation' function ``` -```r +``` r covid_serialint ``` -```{.output} +``` output Disease: COVID-19 Pathogen: SARS-CoV-2 Epi Distribution: serial interval @@ -528,7 +549,7 @@ The `{epiparameter}` combo of `epidist_db()` plus `list_distributions()` list al This is a `` class object: -```r +``` r epiparameter::epidist_db( disease = "ebola", epi_dist = "incubation" @@ -541,7 +562,7 @@ epiparameter::epidist_db( ::::::::::::::::: solution -```r +``` r # the distribution with the highest sample size has a gamma distribution epiparameter::epidist_db( disease = "ebola", @@ -561,12 +582,12 @@ To access the `sample_size`, review an [issue reported in the GitHub repository] We can get the `mean` and standard deviation (`sd`) from this `` diving into the `summary_stats` object: -```r +``` r # get the mean covid_serialint$summary_stats$mean ``` -```{.output} +``` output [1] 4.7 ``` @@ -602,7 +623,7 @@ Use the `str()` to display the structure of the `` R object. :::::::::: solution -```r +``` r # get the sd covid_serialint$summary_stats$sd @@ -619,11 +640,11 @@ covid_serialint$metadata$sample_size An interesting element is the `method_assess` nested entry, which refers to the methods used by the study authors to assess for bias while estimating the serial interval distribution. -```r +``` r covid_serialint$method_assess ``` -```{.output} +``` output $censored [1] TRUE @@ -658,7 +679,7 @@ An informative delay measures the time from symptom onset to recovery or death. ::::::::::::::::: solution -```r +``` r # one way to get the list of all the available parameters epidist_db(disease = "all") %>% list_distributions() %>% @@ -688,7 +709,7 @@ ebola_severity$summary_stats$mean_ci_limits The following output has four entries with different content in the **probability distribution** (`prob_distribution`) column: -```r +``` r distribution <- epiparameter::epidist_db( disease = "covid", @@ -696,29 +717,25 @@ distribution <- ) ``` -```{.output} +``` output Returning 4 results that match the criteria (3 are parameterised). Use subset to filter by entry variables or single_epidist to return a single entry. -To retrieve the short citation for each use the 'get_citation' function +To retrieve the citation for each use the 'get_citation' function ``` -```r +``` r distribution %>% list_distributions() ``` -```{.output} - disease epi_distribution prob_distribution author year -1 COVID-19 serial interval Muluneh .... 2021 -2 COVID-19 serial interval lnorm Hiroshi .... 2020 -3 COVID-19 serial interval weibull Hiroshi .... 2020 -4 COVID-19 serial interval norm Lin Yang.... 2020 +``` error +Error in list_distributions(.): could not find function "list_distributions" ``` Entries with a missing value (``) in the `prob_distribution` column are *non-parameterised* entries. They have summary statistics but no probability distribution. Compare these two outputs: -```r +``` r distribution[[1]]$summary_stats distribution[[1]]$prob_dist ``` @@ -730,7 +747,7 @@ distribution[[1]]$prob_dist As detailed in `?is_parameterised`, a parameterised distribution is the entry that has a probability distribution associated with it provided by an `inference_method` as shown in `metadata`: -```r +``` r distribution[[1]]$metadata$inference_method distribution[[2]]$metadata$inference_method distribution[[4]]$metadata$inference_method @@ -741,17 +758,13 @@ distribution[[4]]$metadata$inference_method In the `epiparameter::list_distributions()` output, we can also find different types of probability distributions (e.g., Log-normal, Weibull, Normal). -```r +``` r distribution %>% list_distributions() ``` -```{.output} - disease epi_distribution prob_distribution author year -1 COVID-19 serial interval Muluneh .... 2021 -2 COVID-19 serial interval lnorm Hiroshi .... 2020 -3 COVID-19 serial interval weibull Hiroshi .... 2020 -4 COVID-19 serial interval norm Lin Yang.... 2020 +``` error +Error in list_distributions(.): could not find function "list_distributions" ``` In `{epiparameter}`, you will mostly find **continuous** distributions like these. You can visualize any of them with the `plot()` function and access to: @@ -760,7 +773,7 @@ In `{epiparameter}`, you will mostly find **continuous** distributions like thes - the *Cumulative Distribution Function (CDF)*. -```r +``` r plot(distribution[[2]]) ``` @@ -769,7 +782,7 @@ plot(distribution[[2]]) With the `day_range` argument, you can change the length or number of days in the `x` axis. Explore what it look like: -```r +``` r plot(distribution[[2]], day_range = 0:20) ``` @@ -819,48 +832,48 @@ If you need it, read in detail about the [R probability functions for the normal If you look at `?stats::Distributions`, each type of distribution has a unique set of functions. However, `{epiparameter}` gives you the same four functions to access each of the values above for any `` object you want! -```r +``` r # plot this to have a visual reference plot(covid_serialint, day_range = 0:20) ``` -```r +``` r # the density value at quantile value of 10 (days) density(covid_serialint, at = 10) ``` -```{.output} +``` output [1] 0.01911607 ``` -```r +``` r # the cumulative probability at quantile value of 10 (days) cdf(covid_serialint, q = 10) ``` -```{.output} +``` output [1] 0.9466605 ``` -```r +``` r # the quantile value (day) at a cumulative probability of 60% quantile(covid_serialint, p = 0.6) ``` -```{.output} +``` output [1] 4.618906 ``` -```r +``` r # generate 10 random values (days) given # the distribution family and its parameters generate(covid_serialint, times = 10) ``` -```{.output} - [1] 3.699438 1.079265 3.468895 7.358368 3.553574 7.719362 1.557383 1.134973 - [9] 9.779812 6.748984 +``` output + [1] 2.566130 2.840252 2.675253 5.179135 4.549618 5.261451 2.861703 1.690031 + [9] 4.533478 4.735248 ``` ::::::::: instructor @@ -893,12 +906,12 @@ In Figure 5 from the [R probability functions for the normal distribution](https ::::::::::::::::: solution -```r +``` r plot(covid_serialint) ``` -```r +``` r cdf(covid_serialint, q = 2) cdf(covid_serialint, q = 6) ``` @@ -920,7 +933,7 @@ If we exchange the question between days and cumulative probability to: - When considering secondary cases, how many days following the symptom onset of primary cases can we expect 55% of symptom onset to occur? -```r +``` r quantile(covid_serialint, p = 0.55) ``` @@ -945,14 +958,14 @@ We can use the set of distribution functions for a _continuous_ distribution (as When we `epiparameter::discretise()` the continuous distribution we get a **discrete**(-ized) distribution: -```r +``` r covid_serialint_discrete <- epiparameter::discretise(covid_serialint) covid_serialint_discrete ``` -```{.output} +``` output Disease: COVID-19 Pathogen: SARS-CoV-2 Epi Distribution: serial interval @@ -975,7 +988,7 @@ Distribution: discrete lnorm While for a **continuous** distribution, we plot the *Probability Density Function (PDF)*, for a **discrete** distribution, we plot the *Probability Mass Function (PMF)*: -```r +``` r # continuous plot(covid_serialint) @@ -986,7 +999,7 @@ plot(covid_serialint_discrete) To finally get a `max` value, let's access the quantile value of the 99.9th percentile or `0.999` probability of the distribution with the `prob_dist$q` notation, similarly to how we access the `summary_stats` values. -```r +``` r covid_serialint_discrete_max <- covid_serialint_discrete$prob_dist$q(p = 0.999) ``` @@ -1008,7 +1021,7 @@ What delay distribution measures the time between infection and the onset of sym The probability function for `` **discrete** distributions differ from the *continuous* ones! -```r +``` r # plot to have a visual reference plot(covid_serialint_discrete, day_range = 0:20) @@ -1030,7 +1043,7 @@ covid_serialint_discrete$prob_dist$r(10) ::::::::::::::::: solution -```r +``` r covid_incubation <- epiparameter::epidist_db( disease = "covid", @@ -1056,7 +1069,7 @@ Now, _Is this result expected in epidemiological terms?_ From a maximum value with `$prob_dist$q()`, we can create a sequence of quantile values as a numeric vector and map density values for each: -```r +``` r # create a discrete distribution visualization # from a maximum value from the distribution covid_serialint_discrete$prob_dist$q(0.999) %>% @@ -1091,7 +1104,7 @@ covid_serialint_discrete$prob_dist$q(0.999) %>% Now we can plug everything into the `EpiNow2::dist_spec()` function! -```r +``` r serial_interval_covid <- dist_spec( mean = covid_serialint$summary_stats$mean, @@ -1099,13 +1112,29 @@ serial_interval_covid <- max = covid_serialint_discrete_max, distribution = "lognormal" ) +``` -serial_interval_covid +``` warning +Warning: `dist_spec()` was deprecated in EpiNow2 1.5.0. +ℹ Please use distribution functions such as `Gamma()` or `Lognormal()` instead. +ℹ The function will become internal only in the next version. +This warning is displayed once every 8 hours. +Call `lifecycle::last_lifecycle_warnings()` to see where this warning was +generated. ``` -```{.output} +``` r +serial_interval_covid +``` - Fixed distribution with PMF [0.18 0.11 0.08 0.066 0.057 0.05 0.045 0.041 0.037 0.034 0.032 0.03 0.028 0.027 0.025 0.024 0.023 0.022 0.021 0.02 0.019 0.019 0.018] +``` output +- lognormal distribution (max: 23): + meanlog: + - fixed value: + 4.7 + sdlog: + - fixed value: + 2.9 ``` :::::::::: callout @@ -1119,7 +1148,7 @@ Using the serial interval instead of the generation time is an alternative that Let's replace the `generation_time` input we used for `EpiNow2::epinow()`. -```r +``` r epinow_estimates <- epinow( # cases reported_cases = example_confirmed[1:60], @@ -1165,7 +1194,7 @@ Key functions we applied in this episode are: -```r +``` r # read data # e.g.: if path to file is data/raw-data/ebola_cases.csv then: ebola_confirmed <- @@ -1261,7 +1290,7 @@ epinow_estimates <- epinow( ::::::::::::::::: solution -```r +``` r covid_incubation <- epiparameter::epidist_db( disease = "covid", epi_dist = "incubation", @@ -1334,7 +1363,7 @@ We can use two complementary delay distributions to estimate the $R_t$ at time $ -```r +``` r # read data # e.g.: if path to file is data/raw-data/ebola_cases.csv then: ebola_confirmed <- @@ -1418,7 +1447,7 @@ How to get the mean and standard deviation from a generation time with median an -```r +``` r # What parameters are available for Influenza? epidist_db(disease = "influenza") %>% list_distributions() %>% @@ -1507,11 +1536,11 @@ However, the methods to accurately estimate delays like the generation interval We can identify what entries in the `{epiparameter}` library assessed for these biases in their methodology with the `method_assess` nested entry: -```r +``` r covid_serialint$method_assess ``` -```{.output} +``` output $censored [1] TRUE @@ -1577,7 +1606,7 @@ Key functions: In this solution we use `purrr::pluck()` to extract elements within the `summary_stats` object which is of class `list`. -```r +``` r # pre-symptomatic transmission epidist_db( disease = "influenza", @@ -1614,7 +1643,7 @@ epidist_db( ``` -```r +``` r # pre-symptomatic transmission epidist_db( disease = "covid", diff --git a/simulating-transmission.md b/simulating-transmission.md index ec182352..987922b5 100644 --- a/simulating-transmission.md +++ b/simulating-transmission.md @@ -41,7 +41,14 @@ Learners should familiarise themselves with following concept dependencies befor Mathematical models are useful tools for generating future trajectories of disease spread. In this tutorial, we will use the R package `{epidemics}` to generate disease trajectories of an influenza strain with pandemic potential. By the end of this tutorial, you will be able to generate the trajectory below showing the number of infectious individuals in different age categories over time. - + +``` error +Error in model_default_cpp(population = uk_population, transmissibility = 1.46/7, : could not find function "model_default_cpp" +``` + +``` error +Error in eval(expr, envir, enclos): object 'output_plot' not found +``` :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: instructor @@ -144,7 +151,7 @@ contact_matrix ## Output -```{.output} +``` output contact.age.group [,1] [,2] [,3] [0,20) 7.883663 2.794154 1.565665 @@ -172,7 +179,7 @@ The initial conditions are the proportion of individuals in each disease state $ The initial conditions in the first age category are $S(0)=1-\frac{1}{1,000,000}$, $E(0) =0$, $I(0)=\frac{1}{1,000,000}$, $R(0)=0$. This is specified as a vector as follows: -```r +``` r initial_i <- 1e-6 initial_conditions_inf <- c( S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0 @@ -182,7 +189,7 @@ initial_conditions_inf <- c( For the age categories that are free from infection, the initial conditions are $S(0)=1$, $E(0) =0$, $I(0)=0$, $R(0)=0$. We specify this as follows, -```r +``` r initial_conditions_free <- c( S = 1, E = 0, I = 0, R = 0, V = 0 ) @@ -191,7 +198,7 @@ initial_conditions_free <- c( We combine the three initial conditions vectors into one matrix, -```r +``` r # combine the initial conditions initial_conditions <- rbind( initial_conditions_inf, # age group 1 @@ -204,7 +211,7 @@ rownames(initial_conditions) <- rownames(contact_matrix) initial_conditions ``` -```{.output} +``` output S E I R V [0,20) 0.999999 0 1e-06 0 0 [20,40) 1.000000 0 0e+00 0 0 @@ -218,13 +225,13 @@ initial_conditions The population object requires a vector containing the demographic structure of the population. The demographic vector must be a named vector containing the number of individuals in each age group of our given population. In this example, we can extract the demographic information from the `contact_data` object that we obtained using the `socialmixr` package. -```r +``` r demography_vector <- contact_data$demography$population names(demography_vector) <- rownames(contact_matrix) demography_vector ``` -```{.output} +``` output [0,20) [20,40) 40+ 14799290 16526302 28961159 ``` @@ -232,7 +239,7 @@ demography_vector To create our population object, we call the function `population()` specifying a name, the contact matrix, the demography vector and the initial conditions. -```r +``` r uk_population <- population( name = "UK", contact_matrix = contact_matrix, @@ -292,13 +299,13 @@ In `epidemics`, the [ODE solver](https://www.boost.org/doc/libs/1_82_0/libs/nume Now we are ready to run our model. Let's load the `epidemics` package : -```r +``` r library(epidemics) ``` Then we specify `time_end=600` to run the model for 600 days. -```r +``` r output <- model_default_cpp( population = uk_population, transmissibility = 1.46 / 7.0, @@ -306,17 +313,18 @@ output <- model_default_cpp( recovery_rate = 1.0 / 7.0, time_end = 600, increment = 1.0 ) +``` + +``` error +Error in model_default_cpp(population = uk_population, transmissibility = 1.46/7, : could not find function "model_default_cpp" +``` + +``` r head(output) ``` -```{.output} - time demography_group compartment value -1 0 [0,20) susceptible 14799275 -2 0 [20,40) susceptible 16526302 -3 0 40+ susceptible 28961159 -4 0 [0,20) exposed 0 -5 0 [20,40) exposed 0 -6 0 40+ exposed 0 +``` error +Error in eval(expr, envir, enclos): object 'output' not found ``` *Note : This model also has the functionality to include vaccination and tracks the number of vaccinated individuals through time. Even though we have not specified any vaccination, there is still a vaccinated compartment in the output (containing no individuals). We will cover the use of vaccination in future tutorials.* @@ -324,7 +332,7 @@ head(output) Our model output consists of the number of individuals in each compartment in each age group through time. We can visualise the infectious individuals only (those in the $I$ class) through time. -```r +``` r filter(output_plot, compartment %in% c("exposed", "infectious")) %>% ggplot( aes( @@ -363,7 +371,9 @@ filter(output_plot, compartment %in% c("exposed", "infectious")) %>% ) ``` - +``` error +Error in eval(expr, envir, enclos): object 'output_plot' not found +``` ::::::::::::::::::::::::::::::::::::: callout @@ -384,14 +394,14 @@ We ran our model with $R_0= 1.5$. However, we believe that $R_0$ follows a norma 1. Obtain 100 samples from the from a normal distribution -```r +``` r R0_vec <- rnorm(100, 1.5, 0.05) ``` 2. Run the model 100 times with $R_0$ equal to a different sample each time -```r +``` r output_samples <- Map( R0_vec, seq_along(R0_vec), @@ -410,16 +420,26 @@ output_samples <- Map( output } ) +``` +``` error +Error in model_default_cpp(population = uk_population, transmissibility = x/7, : could not find function "model_default_cpp" +``` + +``` r # combine to prepare for plotting output_samples <- bind_rows(output_samples) ``` +``` error +Error in eval(expr, envir, enclos): object 'output_samples' not found +``` + 3. Calculate the mean and 95% quantiles of number of infectious individuals across each model simulation and visualise output -```r +``` r ggplot( output_samples[output_samples$compartment == "infectious", ], aes(time, value) @@ -447,7 +467,9 @@ ggplot( ) ``` - +``` error +Error in eval(expr, envir, enclos): object 'output_samples' not found +``` Deciding which parameters to include uncertainty in depends on a few factors: how well informed a parameter value is e.g. consistency of estimates from the literature; how sensitive model outputs are to parameter value changes; and the purpose of the modelling task. See [McCabe et al. 2021](https://doi.org/10.1016%2Fj.epidem.2021.100520) to learn about different types of uncertainty in infectious disease modelling.