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(incidence2)
library(i2extras)
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"
)[1:20, ]
dat
#> An incidence object: 20 x 2
#> date range: [2014-W15] to [2014-W34]
#> cases: 783
#> interval: 1 (Monday) week
#>
#> date_index count
#> <yrwk> <int>
#> 1 2014-W15 1
#> 2 2014-W16 1
#> 3 2014-W17 5
#> 4 2014-W18 4
#> 5 2014-W19 12
#> 6 2014-W20 17
#> 7 2014-W21 15
#> 8 2014-W22 19
#> 9 2014-W23 23
#> 10 2014-W24 21
#> 11 2014-W25 30
#> 12 2014-W26 22
#> 13 2014-W27 34
#> 14 2014-W28 38
#> 15 2014-W29 61
#> 16 2014-W30 59
#> 17 2014-W31 80
#> 18 2014-W32 86
#> 19 2014-W33 116
#> 20 2014-W34 139
plot(dat)
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: 1 x 8
#> count_variable data model estimates fitting_warning fitting_error
#> <chr> <list<tibble[> <list> <list> <list> <list>
#> 1 count [20 × 2] <glm> <df [20 × … <NULL> <NULL>
#> # … with 2 more variables: prediction_warning <list>, prediction_error <list>
plot(out)
growth_rate(out)
#> # A tibble: 1 x 9
#> count_variable model r r_lower r_upper growth_or_decay time time_lower
#> <chr> <list> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
#> 1 count <glm> 0.184 0.168 0.200 doubling 3.77 3.46
#> # … with 1 more variable: time_upper <dbl>To unnest the data we can use unnest() (a function that has been reexported from the tidyr package.
unnest(out, estimates)
#> # A tibble: 20 x 14
#> count_variable data model count date_index estimate lower_ci upper_ci
#> <chr> <list<tibbl> <lis> <int> <yrwk> <dbl> <dbl> <dbl>
#> 1 count [20 × 2] <glm> 1 2014-W15 4.09 3.20 5.23
#> 2 count [20 × 2] <glm> 1 2014-W16 4.92 3.91 6.19
#> 3 count [20 × 2] <glm> 5 2014-W17 5.92 4.77 7.33
#> 4 count [20 × 2] <glm> 4 2014-W18 7.11 5.82 8.68
#> 5 count [20 × 2] <glm> 12 2014-W19 8.55 7.11 10.3
#> 6 count [20 × 2] <glm> 17 2014-W20 10.3 8.67 12.2
#> 7 count [20 × 2] <glm> 15 2014-W21 12.3 10.6 14.4
#> 8 count [20 × 2] <glm> 19 2014-W22 14.8 12.9 17.1
#> 9 count [20 × 2] <glm> 23 2014-W23 17.8 15.7 20.3
#> 10 count [20 × 2] <glm> 21 2014-W24 21.4 19.1 24.0
#> 11 count [20 × 2] <glm> 30 2014-W25 25.8 23.3 28.5
#> 12 count [20 × 2] <glm> 22 2014-W26 31.0 28.3 33.9
#> 13 count [20 × 2] <glm> 34 2014-W27 37.2 34.3 40.4
#> 14 count [20 × 2] <glm> 38 2014-W28 44.8 41.5 48.2
#> 15 count [20 × 2] <glm> 61 2014-W29 53.8 50.1 57.7
#> 16 count [20 × 2] <glm> 59 2014-W30 64.7 60.3 69.4
#> 17 count [20 × 2] <glm> 80 2014-W31 77.7 72.2 83.7
#> 18 count [20 × 2] <glm> 86 2014-W32 93.4 86.2 101.
#> 19 count [20 × 2] <glm> 116 2014-W33 112. 103. 123.
#> 20 count [20 × 2] <glm> 139 2014-W34 135. 122. 149.
#> # … with 6 more variables: lower_pi <int>, upper_pi <int>,
#> # fitting_warning <list>, fitting_error <list>, prediction_warning <list>,
#> # prediction_error <list>fit_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
#> An incidence object: 120 x 3
#> date range: [2014-W15] to [2014-W38]
#> cases: 1621
#> interval: 1 (Monday) week
#>
#> date_index hospital count
#> <yrwk> <fct> <int>
#> 1 2014-W15 Military Hospital 1
#> 2 2014-W16 Connaught Hospital 1
#> 3 2014-W17 <NA> 2
#> 4 2014-W17 other 3
#> 5 2014-W18 <NA> 1
#> 6 2014-W18 Connaught Hospital 1
#> 7 2014-W18 Princess Christian Maternity Hospital (PCMH) 1
#> 8 2014-W18 Rokupa Hospital 1
#> 9 2014-W19 <NA> 1
#> 10 2014-W19 Connaught Hospital 3
#> # … with 110 more rows
out <- fit_curve(grouped_dat, model = "poisson", alpha = 0.05)
out
#> # A tibble: 6 x 9
#> count_variable hospital data model estimates fitting_warning fitting_error
#> <chr> <fct> <list<t> <lis> <list> <list> <list>
#> 1 count Connaug… [22 × 2] <glm> <df [22 … <NULL> <NULL>
#> 2 count Militar… [21 × 2] <glm> <df [21 … <NULL> <NULL>
#> 3 count other [20 × 2] <glm> <df [20 … <NULL> <NULL>
#> 4 count Princes… [17 × 2] <glm> <df [17 … <NULL> <NULL>
#> 5 count Rokupa … [18 × 2] <glm> <df [18 … <NULL> <NULL>
#> 6 count <NA> [22 × 2] <glm> <df [22 … <NULL> <NULL>
#> # … with 2 more variables: prediction_warning <list>, prediction_error <list>
# plot with a prediction interval but not a confidence interval
plot(out, ci = FALSE, pi=TRUE)
growth_rate(out)
#> # A tibble: 6 x 10
#> count_variable hospital model r r_lower r_upper growth_or_decay time
#> <chr> <fct> <lis> <dbl> <dbl> <dbl> <chr> <dbl>
#> 1 count Connaught Ho… <glm> 0.197 0.177 0.217 doubling 3.53
#> 2 count Military Hos… <glm> 0.173 0.147 0.200 doubling 4.00
#> 3 count other <glm> 0.170 0.141 0.200 doubling 4.09
#> 4 count Princess Chr… <glm> 0.142 0.101 0.188 doubling 4.87
#> 5 count Rokupa Hospi… <glm> 0.178 0.133 0.228 doubling 3.89
#> 6 count <NA> <glm> 0.184 0.164 0.205 doubling 3.77
#> # … with 2 more variables: time_lower <dbl>, time_upper <dbl>We 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 x 7
#> count_variable hospital data model estimates fitting_warning
#> <chr> <fct> <list<tibb> <lis> <list> <list>
#> 1 count Connaught Hospital [22 × 2] <neg… <df [22 ×… <chr [2]>
#> 2 count other [20 × 2] <neg… <df [20 ×… <chr [2]>
#> 3 count Princess Christia… [17 × 2] <neg… <df [17 ×… <chr [2]>
#> 4 count Rokupa Hospital [18 × 2] <neg… <df [18 ×… <chr [2]>
#> 5 count <NA> [22 × 2] <neg… <df [22 ×… <chr [2]>
#> # … with 1 more variable: prediction_warning <list>
unnest(is_warning(out), fitting_warning)
#> # A tibble: 10 x 7
#> count_variable hospital data model estimates fitting_warning
#> <chr> <fct> <list<tibb> <list> <list> <chr>
#> 1 count Connaught Hospi… [22 × 2] <negb… <df [22 … iteration limit…
#> 2 count Connaught Hospi… [22 × 2] <negb… <df [22 … iteration limit…
#> 3 count other [20 × 2] <negb… <df [20 … iteration limit…
#> 4 count other [20 × 2] <negb… <df [20 … iteration limit…
#> 5 count Princess Christ… [17 × 2] <negb… <df [17 … iteration limit…
#> 6 count Princess Christ… [17 × 2] <negb… <df [17 … iteration limit…
#> 7 count Rokupa Hospital [18 × 2] <negb… <df [18 … iteration limit…
#> 8 count Rokupa Hospital [18 × 2] <negb… <df [18 … iteration limit…
#> 9 count <NA> [22 × 2] <negb… <df [22 … iteration limit…
#> 10 count <NA> [22 × 2] <negb… <df [22 … iteration limit…
#> # … with 1 more variable: prediction_warning <list>Rolling 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, before = 2) # group observations with the 2 prior
ra
#> # A tibble: 6 x 3
#> count_variable hospital rolling_average
#> <chr> <fct> <list>
#> 1 count Connaught Hospital <tibble [22 × 3]>
#> 2 count Military Hospital <tibble [21 × 3]>
#> 3 count other <tibble [20 × 3]>
#> 4 count Princess Christian Maternity Hospital (PCMH) <tibble [17 × 3]>
#> 5 count Rokupa Hospital <tibble [18 × 3]>
#> 6 count <NA> <tibble [22 × 3]>
unnest(ra, rolling_average)
#> # A tibble: 120 x 5
#> count_variable hospital date_index count rolling_average
#> <chr> <fct> <yrwk> <int> <dbl>
#> 1 count Connaught Hospital 2014-W16 1 NA
#> 2 count Connaught Hospital 2014-W18 1 NA
#> 3 count Connaught Hospital 2014-W19 3 1.67
#> 4 count Connaught Hospital 2014-W20 2 2
#> 5 count Connaught Hospital 2014-W21 5 3.33
#> 6 count Connaught Hospital 2014-W22 4 3.67
#> 7 count Connaught Hospital 2014-W23 6 5
#> 8 count Connaught Hospital 2014-W24 6 5.33
#> 9 count Connaught Hospital 2014-W25 12 8
#> 10 count Connaught Hospital 2014-W26 8 8.67
#> # … with 110 more rows
plot(ra, color = "white")
#> Warning: Removed 12 rows containing missing values (position_stack).
#> Warning: Removed 12 rows containing missing values (position_stack).