Simulation tests for growth rate estimators

This test data is based around a known time varying exponential growth rate with an initial epidemic seed size of 100. The data is simulated with weekly periodicity.

Locfit models

Simple incidence test with a poisson model

An incidence mode based on absolute counts:

data = sim_poisson_model()
data %>% dplyr::glimpse()
#> Rows: 105
#> Columns: 6
#> Groups: statistic [1]
#> $ time      <time_prd> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
#> $ growth    <dbl> 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, …
#> $ imports   <dbl> 100, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ rate      <dbl> 100.0000, 110.5171, 122.1403, 134.9859, 149.1825, 164.8721, …
#> $ count     <int> 94, 116, 124, 144, 137, 165, 203, 211, 229, 262, 294, 282, 3…
#> $ statistic <chr> "infections", "infections", "infections", "infections", "inf…
tmp = data %>% poisson_locfit_model(window=7, deg = 2)

plot_incidence(tmp, data)+ggplot2::geom_line(
  mapping=ggplot2::aes(x=as.Date(time),y=rate), data=data, colour="red",inherit.aes = FALSE)

Estimated absolute growth rate versus simulation (red)


plot_growth_rate(tmp)+
  ggplot2::geom_line(mapping=ggplot2::aes(x=as.Date(time),y=growth), data=data, colour="red",inherit.aes = FALSE)


plot_growth_phase(tmp)

Multinomial data

Multiple classes are simulated as 3 independent epidemics (‘variant1’, ‘variant2’ and ‘variant3’) with known growth rates and initial sample size resulting in 3 parallel time series. These are combined to give an overall epidemic and a proportional distribution of each ‘variant’ as a fraction of the whole. A relative growth rate is calculated based on set parameters.

data2 = sim_multinomial() %>% dplyr::group_by(class) %>% dplyr::glimpse()
#> Rows: 315
#> Columns: 10
#> Groups: class [3]
#> $ time            <time_prd> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,…
#> $ growth          <dbl> 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,…
#> $ imports         <dbl> 100, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ rate            <dbl> 100.0000, 110.5171, 122.1403, 134.9859, 149.1825, 164.…
#> $ count           <int> 118, 99, 129, 144, 156, 175, 198, 214, 227, 258, 293, …
#> $ statistic       <chr> "infections", "infections", "infections", "infections"…
#> $ class           <chr> "variant1", "variant1", "variant1", "variant1", "varia…
#> $ proportion      <dbl> 0.3333333, 0.3382826, 0.3420088, 0.3445125, 0.3458146,…
#> $ proportion.obs  <dbl> 0.3333333, 0.3390411, 0.3403694, 0.3453237, 0.3451327,…
#> $ relative.growth <dbl> 0.025000000, 0.019385523, 0.013833622, 0.008404115, 0.…

Poisson model

Firstly fitting the same incidence model in a groupwise fashion:

tmp2 = data2 %>% poisson_locfit_model(window=7, deg = 1)

plot_incidence(tmp2, data2)+scale_y_log1p()
#> Scale for y is already present.
#> Adding another scale for y, which will replace the existing scale.

And the absolute growth rates:

plot_growth_rate(modelled = tmp2)+
   ggplot2::geom_line(mapping=ggplot2::aes(x=as.Date(time),y=growth, colour=class), data=data2, inherit.aes = FALSE)+
   ggplot2::facet_wrap(dplyr::vars(class), ncol=1)

One versus others Binomial model

This looks at the proportions of the three variants and their growth rate relative to each other:

# This will reinterpret total to be the total of positives across all variants
data3 = data2 %>% 
  dplyr::group_by(time) %>% 
  dplyr::mutate(denom = sum(count)) %>%
  dplyr::group_by(class) %>%
  dplyr::glimpse()
#> Rows: 315
#> Columns: 11
#> Groups: class [3]
#> $ time            <time_prd> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,…
#> $ growth          <dbl> 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,…
#> $ imports         <dbl> 100, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ rate            <dbl> 100.0000, 110.5171, 122.1403, 134.9859, 149.1825, 164.…
#> $ count           <int> 118, 99, 129, 144, 156, 175, 198, 214, 227, 258, 293, …
#> $ statistic       <chr> "infections", "infections", "infections", "infections"…
#> $ class           <chr> "variant1", "variant1", "variant1", "variant1", "varia…
#> $ proportion      <dbl> 0.3333333, 0.3382826, 0.3420088, 0.3445125, 0.3458146,…
#> $ proportion.obs  <dbl> 0.3333333, 0.3390411, 0.3403694, 0.3453237, 0.3451327,…
#> $ relative.growth <dbl> 0.025000000, 0.019385523, 0.013833622, 0.008404115, 0.…
#> $ denom           <int> 354, 292, 379, 417, 452, 506, 575, 624, 667, 767, 880,…

Firstly proportions:


tmp3 = data3 %>% proportion_locfit_model(window=14, deg = 2)

plot_proportion(modelled = tmp3,raw = data3)+
  ggplot2::facet_wrap(dplyr::vars(class), ncol=1)

And secondly relative growth rate:



plot_growth_rate(modelled = tmp3)+
   ggplot2::geom_line(mapping=ggplot2::aes(x=as.Date(time),y=relative.growth, colour=class), data=data2, inherit.aes = FALSE)+
   ggplot2::facet_wrap(dplyr::vars(class), ncol=1)


plot_growth_phase(tmp3)

Multinomial model

The multinomial model gives us absolute proportions only (and no growth rates)

# we don't need to calculate the denominator as it is done automatically by the 
# multinomial model

tmp4 = data2 %>% multinomial_nnet_model()
#> # weights:  30 (18 variable)
#> initial  value 354566.130045 
#> iter  10 value 178654.226096
#> iter  20 value 175460.540475
#> final  value 173863.879898 
#> converged
plot_multinomial(tmp4)


# plot_multinomial(tmp3, events = event_test,normalise = TRUE)

GLM models

Poisson model

  • Spline currently only good for incidence

tmp5 = data %>% poisson_glm_model(window=7)
plot_incidence(tmp5,data)

Binomial model

Absolute proportions only


tmp6 = data3 %>% proportion_glm_model(window=14, deg = 2)
plot_proportion(tmp6,data3)