Example
To illustrate the trend fitting functionality of i2extras we will use the simulated Ebola Virus Disease (EVD) outbreak data from the outbreaks package.
Loading relevant packages and data
library(outbreaks)
library(i2extras)
#> Loading required package: incidence2
#> Loading required package: grates
library(ggplot2)
raw_dat <- ebola_sim_clean$linelistFor this example we will restrict ourselves to the first 20 weeks of data:
dat <- incidence(
raw_dat,
date_index = "date_of_onset",
interval = "week",
groups = "gender"
)[1:20, ]
dat
#> # incidence: 20 x 4
#> # count vars: date_of_onset
#> # groups: gender
#> date_index gender count_variable count
#> * <isowk> <fct> <chr> <int>
#> 1 2014-W15 f date_of_onset 1
#> 2 2014-W16 m date_of_onset 1
#> 3 2014-W17 f date_of_onset 4
#> 4 2014-W17 m date_of_onset 1
#> 5 2014-W18 f date_of_onset 4
#> 6 2014-W19 f date_of_onset 9
#> 7 2014-W19 m date_of_onset 3
#> 8 2014-W20 f date_of_onset 7
#> 9 2014-W20 m date_of_onset 10
#> 10 2014-W21 f date_of_onset 8
#> 11 2014-W21 m date_of_onset 7
#> 12 2014-W22 f date_of_onset 9
#> 13 2014-W22 m date_of_onset 10
#> 14 2014-W23 f date_of_onset 13
#> 15 2014-W23 m date_of_onset 10
#> 16 2014-W24 f date_of_onset 7
#> 17 2014-W24 m date_of_onset 14
#> 18 2014-W25 f date_of_onset 15
#> 19 2014-W25 m date_of_onset 15
#> 20 2014-W26 f date_of_onset 12
plot(dat, angle = 45, border_colour = "white")
Modeling incidence
We can use fit_curve() to fit the data with either a
poisson or negative binomial regression model. The output from this will
be a nested object with class incidence2_fit which has
methods available for both automatic plotting and the calculation of
growth (decay) rates and doubling (halving) times.
out <- fit_curve(dat, model = "poisson", alpha = 0.05)
out
#> # A tibble: 2 × 9
#> count_varia…¹ gender data model estimates fitti…² fitti…³ predi…⁴ predi…⁵
#> <chr> <fct> <list<t> <lis> <list> <list> <list> <list> <list>
#> 1 date_of_onset f [11 × 2] <glm> <trndng_p> <NULL> <NULL> <NULL> <NULL>
#> 2 date_of_onset m [9 × 2] <glm> <trndng_p> <NULL> <NULL> <NULL> <NULL>
#> # … with abbreviated variable names ¹count_variable, ²fitting_warning,
#> # ³fitting_error, ⁴prediction_warning, ⁵prediction_error
plot(out, angle = 45, border_colour = "white")
growth_rate(out)
#> # A tibble: 2 × 10
#> count_varia…¹ gender model r r_lower r_upper growt…² time time_…³ time_…⁴
#> <chr> <fct> <lis> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 date_of_onset f <glm> 0.137 0.0698 0.206 doubli… 5.07 3.36 9.93
#> 2 date_of_onset m <glm> 0.240 0.146 0.341 doubli… 2.89 2.03 4.75
#> # … with abbreviated variable names ¹count_variable, ²growth_or_decay,
#> # ³time_lower, ⁴time_upperTo unnest the data we can use unnest() (a function that
has been reexported from the tidyr package.
unnest(out, estimates)
#> # A tibble: 20 × 15
#> count_v…¹ gender data model count date_…² estim…³ lower…⁴ upper…⁵ lower…⁶
#> <chr> <fct> <list<t> <lis> <int> <isowk> <dbl> <dbl> <dbl> <dbl>
#> 1 date_of_… f [11 × 2] <glm> 1 2014-W… 3.27 1.90 5.61 0
#> 2 date_of_… f [11 × 2] <glm> 4 2014-W… 4.30 2.83 6.52 0
#> 3 date_of_… f [11 × 2] <glm> 4 2014-W… 4.93 3.44 7.06 0
#> 4 date_of_… f [11 × 2] <glm> 9 2014-W… 5.65 4.15 7.68 1
#> 5 date_of_… f [11 × 2] <glm> 7 2014-W… 6.47 4.99 8.40 1
#> 6 date_of_… f [11 × 2] <glm> 8 2014-W… 7.42 5.92 9.31 2
#> 7 date_of_… f [11 × 2] <glm> 9 2014-W… 8.51 6.90 10.5 2
#> 8 date_of_… f [11 × 2] <glm> 13 2014-W… 9.75 7.88 12.1 3
#> 9 date_of_… f [11 × 2] <glm> 7 2014-W… 11.2 8.82 14.2 4
#> 10 date_of_… f [11 × 2] <glm> 15 2014-W… 12.8 9.72 16.9 4
#> 11 date_of_… f [11 × 2] <glm> 12 2014-W… 14.7 10.6 20.4 5
#> 12 date_of_… m [9 × 2] <glm> 1 2014-W… 2.01 1.02 3.95 0
#> 13 date_of_… m [9 × 2] <glm> 1 2014-W… 2.56 1.43 4.59 0
#> 14 date_of_… m [9 × 2] <glm> 3 2014-W… 4.13 2.73 6.24 0
#> 15 date_of_… m [9 × 2] <glm> 10 2014-W… 5.25 3.75 7.36 1
#> 16 date_of_… m [9 × 2] <glm> 7 2014-W… 6.67 5.07 8.79 1
#> 17 date_of_… m [9 × 2] <glm> 10 2014-W… 8.48 6.69 10.8 2
#> 18 date_of_… m [9 × 2] <glm> 10 2014-W… 10.8 8.50 13.7 3
#> 19 date_of_… m [9 × 2] <glm> 14 2014-W… 13.7 10.4 18.0 5
#> 20 date_of_… m [9 × 2] <glm> 15 2014-W… 17.4 12.4 24.4 6
#> # … with 5 more variables: upper_pi <dbl>, fitting_warning <list>,
#> # fitting_error <list>, prediction_warning <list>, prediction_error <list>,
#> # and abbreviated variable names ¹count_variable, ²date_index, ³estimate,
#> # ⁴lower_ci, ⁵upper_ci, ⁶lower_pifit_curve() also works nicely with grouped
incidence2 objects. In this situation, we return a nested
tibble with some additional columns that indicate whether there has been
a warning or error during the fitting or prediction stages.
grouped_dat <- incidence(
raw_dat,
date_index = "date_of_onset",
interval = "week",
groups = "hospital"
)[1:120, ]
grouped_dat
#> # incidence: 120 x 4
#> # count vars: date_of_onset
#> # groups: hospital
#> date_index hospital count_variable count
#> * <isowk> <fct> <chr> <int>
#> 1 2014-W15 Military Hospital date_of_onset 1
#> 2 2014-W16 Connaught Hospital date_of_onset 1
#> 3 2014-W17 <NA> date_of_onset 2
#> 4 2014-W17 other date_of_onset 3
#> 5 2014-W18 <NA> date_of_onset 1
#> 6 2014-W18 Connaught Hospital date_of_onset 1
#> 7 2014-W18 Princess Christian Maternity Hospital (PCMH) date_of_onset 1
#> 8 2014-W18 Rokupa Hospital date_of_onset 1
#> 9 2014-W19 <NA> date_of_onset 1
#> 10 2014-W19 Connaught Hospital date_of_onset 3
#> # … with 110 more rows
out <- fit_curve(grouped_dat, model = "poisson", alpha = 0.05)
out
#> # A tibble: 6 × 9
#> count_vari…¹ hospi…² data model estimates fitti…³ fitti…⁴ predi…⁵ predi…⁶
#> <chr> <fct> <list<t> <lis> <list> <list> <list> <list> <list>
#> 1 date_of_ons… Connau… [22 × 2] <glm> <trndng_p> <NULL> <NULL> <NULL> <NULL>
#> 2 date_of_ons… Milita… [21 × 2] <glm> <trndng_p> <NULL> <NULL> <NULL> <NULL>
#> 3 date_of_ons… other [20 × 2] <glm> <trndng_p> <NULL> <NULL> <NULL> <NULL>
#> 4 date_of_ons… Prince… [17 × 2] <glm> <trndng_p> <NULL> <NULL> <NULL> <NULL>
#> 5 date_of_ons… Rokupa… [18 × 2] <glm> <trndng_p> <NULL> <NULL> <NULL> <NULL>
#> 6 date_of_ons… <NA> [22 × 2] <glm> <trndng_p> <NULL> <NULL> <NULL> <NULL>
#> # … with abbreviated variable names ¹count_variable, ²hospital,
#> # ³fitting_warning, ⁴fitting_error, ⁵prediction_warning, ⁶prediction_error
# plot with a prediction interval but not a confidence interval
plot(out, ci = FALSE, pi=TRUE, angle = 45, border_colour = "white")
growth_rate(out)
#> # A tibble: 6 × 10
#> count_vari…¹ hospi…² model r r_lower r_upper growt…³ time time_…⁴ time_…⁵
#> <chr> <fct> <lis> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 date_of_ons… Connau… <glm> 0.197 0.177 0.217 doubli… 3.53 3.20 3.92
#> 2 date_of_ons… Milita… <glm> 0.173 0.147 0.200 doubli… 4.00 3.46 4.71
#> 3 date_of_ons… other <glm> 0.170 0.141 0.200 doubli… 4.09 3.46 4.92
#> 4 date_of_ons… Prince… <glm> 0.142 0.101 0.188 doubli… 4.87 3.70 6.87
#> 5 date_of_ons… Rokupa… <glm> 0.178 0.133 0.228 doubli… 3.89 3.04 5.22
#> 6 date_of_ons… <NA> <glm> 0.184 0.164 0.205 doubli… 3.77 3.39 4.24
#> # … with abbreviated variable names ¹count_variable, ²hospital,
#> # ³growth_or_decay, ⁴time_lower, ⁵time_upperWe provide helper functions, is_ok(),
is_warning() and is_error() to help filter the
output as necessary.
out <- fit_curve(grouped_dat, model = "negbin", alpha = 0.05)
is_warning(out)
#> # A tibble: 5 × 7
#> count_variable hospital data model estimates fitti…¹ predi…²
#> <chr> <fct> <list<t> <list> <list> <list> <list>
#> 1 date_of_onset Connaught Hospital [22 × 2] <negbin> <trndng_p> <chr> <NULL>
#> 2 date_of_onset other [20 × 2] <negbin> <trndng_p> <chr> <NULL>
#> 3 date_of_onset Princess Christia… [17 × 2] <negbin> <trndng_p> <chr> <NULL>
#> 4 date_of_onset Rokupa Hospital [18 × 2] <negbin> <trndng_p> <chr> <NULL>
#> 5 date_of_onset <NA> [22 × 2] <negbin> <trndng_p> <chr> <NULL>
#> # … with abbreviated variable names ¹fitting_warning, ²prediction_warning
unnest(is_warning(out), fitting_warning)
#> # A tibble: 10 × 7
#> count_variable hospital data model estimates fitti…¹ predi…²
#> <chr> <fct> <list<t> <list> <list> <chr> <list>
#> 1 date_of_onset Connaught Hospit… [22 × 2] <negbin> <trndng_p> iterat… <NULL>
#> 2 date_of_onset Connaught Hospit… [22 × 2] <negbin> <trndng_p> iterat… <NULL>
#> 3 date_of_onset other [20 × 2] <negbin> <trndng_p> iterat… <NULL>
#> 4 date_of_onset other [20 × 2] <negbin> <trndng_p> iterat… <NULL>
#> 5 date_of_onset Princess Christi… [17 × 2] <negbin> <trndng_p> iterat… <NULL>
#> 6 date_of_onset Princess Christi… [17 × 2] <negbin> <trndng_p> iterat… <NULL>
#> 7 date_of_onset Rokupa Hospital [18 × 2] <negbin> <trndng_p> iterat… <NULL>
#> 8 date_of_onset Rokupa Hospital [18 × 2] <negbin> <trndng_p> iterat… <NULL>
#> 9 date_of_onset <NA> [22 × 2] <negbin> <trndng_p> iterat… <NULL>
#> 10 date_of_onset <NA> [22 × 2] <negbin> <trndng_p> iterat… <NULL>
#> # … with abbreviated variable names ¹fitting_warning, ²prediction_warningRolling average
We can add a rolling average, across current and previous intervals,
to an incidence2 object with the
add_rolling_average() function:
ra <- add_rolling_average(grouped_dat, n = 2L) # group observations with the 2 prior
plot(ra, border_colour = "white", angle = 45) +
geom_line(aes(x = date_index, y = rolling_average))
#> Warning: Removed 1 row containing missing values (`geom_line()`).