Merge branch 'main' into 43-do-inference-on-generated-data-using-this…

…-in-conjunction-with-the-already-existing-rw-latent-process

EpiAware/src/EpiAware.jl

@@ -31,18 +31,10 @@ using Distributions,
     QuadGK
 
 # Exported utilities
-export scan,
-    create_discrete_pmf,
-    growth_rate_to_reproductive_ratio,
-    generate_observation_kernel,
-    default_rw_priors,
-    default_delay_obs_priors,
-    neg_MGF,
-    dneg_MGF_dr,
-    R_to_r
+export create_discrete_pmf, default_rw_priors, default_delay_obs_priors
 
 # Exported types
-export EpiData, Renewal, ExpGrowthRate, DirectInfections, AbstractEpiModel
+export EpiData, Renewal, ExpGrowthRate, DirectInfections
 
 # Exported Turing model constructors
 export make_epi_inference_model, random_walk, delay_observations

EpiAware/src/utilities.jl

@@ -173,7 +173,7 @@ end
 
 
 """
-    growth_rate_to_reproductive_ratio(r, w)
+    r_to_R(r, w)
 
 Compute the reproductive ratio given exponential growth rate `r`
     and discretized generation interval `w`.
@@ -185,7 +185,7 @@ Compute the reproductive ratio given exponential growth rate `r`
 # Returns
 - The reproductive ratio.
 """
-function growth_rate_to_reproductive_ratio(r, w::AbstractVector)
+function r_to_R(r, w::AbstractVector)
     return 1 / neg_MGF(r, w::AbstractVector)
 end
 

EpiAware/test/predictive_checking/fast_approx_for_r.jl

@@ -29,10 +29,17 @@ To rapidly improve the estimate for `r` we use Newton steps given by
 r_{n+1} = r_n - {\mathcal{R}_0 G(r) - 1\over \mathcal{R}_0 G'(r)}.
 ```
 
-=#
+### Test mode
+
+Run this script in Test environment mode. If not, run the following command:
 
+```julia
 using TestEnv
 TestEnv.activate()
+```
+=#
+
+
 using EpiAware
 using Distributions
 using StatsPlots
@@ -52,10 +59,10 @@ doubling_times = [1.0, 3.5, 7.0, 14.0]
 
 errors = mapreduce(hcat, doubling_times) do T_2
     true_r = log(2) / T_2 # 7 day doubling time
-    R0 = growth_rate_to_reproductive_ratio(true_r, w)
+    R0 = EpiAware.r_to_R(true_r, w)
 
     return map(idxs) do ns
-        @time r = R_to_r(R0, w, newton_steps = ns)
+        @time r = EpiAware.R_to_r(R0, w, newton_steps = ns)
         abs(r - true_r) + jitter
     end
 end

EpiAware/test/predictive_checking/toy_model_log_infs_RW.jl

@@ -48,17 +48,18 @@ r &\sim \text{Gamma}(3, 0.05/3).
 
 ## Load dependencies
 
-This script should be run from the root folder of `EpiAware` and with the active environment.
-
-NB: This script is intended to be run in the test environment.
+This script should be run from Test environment mode. If not, run the following command:
 
 ```julia
 using TestEnv # Run in Test environment mode
 TestEnv.activate()
 ```
+
 =#
+
 # using TestEnv # Run in Test environment mode
 # TestEnv.activate()
+
 using EpiAware
 using Turing
 using Distributions

EpiAware/test/test_utilities.jl

@@ -7,7 +7,7 @@
     xs = [1, 2, 3, 4, 5]
     expected_ys = [1, 3, 6, 10, 15]
     expected_carry = 15
-    ys, carry = scan(add, 0, xs)
+    ys, carry = EpiAware.scan(add, 0, xs)
     @test ys == expected_ys
     @test carry == expected_carry
 end
@@ -22,7 +22,7 @@ end
     expected_ys = [1, 2, 6, 24, 120]
     expected_carry = 120
 
-    ys, carry = scan(multiply, 1, xs)
+    ys, carry = EpiAware.scan(multiply, 1, xs)
     @test ys == expected_ys
     @test carry == expected_carry
 end
@@ -81,21 +81,21 @@ end
 
 end
 
-@testitem "Testing growth_rate_to_reproductive_ratio function" begin
+@testitem "Testing r_to_R function" begin
     #Test that zero exp growth rate imples R0 = 1
     @testset "Test case 1" begin
         r = 0
         w = ones(5) |> x -> x ./ sum(x)
         expected_ratio = 1
-        ratio = growth_rate_to_reproductive_ratio(r, w)
+        ratio = EpiAware.r_to_R(r, w)
         @test ratio ≈ expected_ratio atol = 1e-15
     end
 
     #Test MethodError when w is not a vector
     @testset "Test case 2" begin
         r = 0
         w = 1
-        @test_throws MethodError growth_rate_to_reproductive_ratio(r, w)
+        @test_throws MethodError EpiAware.r_to_R(r, w)
     end
 
 end
@@ -114,7 +114,7 @@ end
                 0 0 0.3 0.5 0.2
             ],
         )
-        K = generate_observation_kernel(delay_int, time_horizon)
+        K = EpiAware.generate_observation_kernel(delay_int, time_horizon)
         @test K == expected_K
     end
 
@@ -126,7 +126,7 @@ end
         r = 0.5
         w = [0.2, 0.3, 0.5]
         expected_result = 0.2 * exp(-0.5 * 1) + 0.3 * exp(-0.5 * 2) + 0.5 * exp(-0.5 * 3)
-        result = neg_MGF(r, w)
+        result = EpiAware.neg_MGF(r, w)
         @test result ≈ expected_result atol = 1e-15
     end
 
@@ -136,20 +136,21 @@ end
         w = [0.1, 0.2, 0.3, 0.4]
         expected_result =
             0.1 * exp(-0 * 1) + 0.2 * exp(-0 * 2) + 0.3 * exp(-0 * 3) + 0.4 * exp(-0 * 4)
-        result = neg_MGF(r, w)
+        result = EpiAware.neg_MGF(r, w)
         @test result ≈ expected_result atol = 1e-15
     end
 
 end
 
 @testitem "Testing dneg_MGF_dr function" begin
+
     # Test case 1: Testing with positive r and non-empty weight vector
     @testset "Test case 1" begin
         r = 0.5
         w = [0.2, 0.3, 0.5]
         expected_result =
             -(0.2 * 1 * exp(-0.5 * 1) + 0.3 * 2 * exp(-0.5 * 2) + 0.5 * 3 * exp(-0.5 * 3))
-        result = dneg_MGF_dr(r, w)
+        result = EpiAware.dneg_MGF_dr(r, w)
         @test result ≈ expected_result atol = 1e-15
     end
 
@@ -163,7 +164,7 @@ end
             0.3 * 3 * exp(-0 * 3) +
             0.4 * 4 * exp(-0 * 4)
         )
-        result = dneg_MGF_dr(r, w)
+        result = EpiAware.dneg_MGF_dr(r, w)
         @test result ≈ expected_result atol = 1e-15
     end
 

manuscript/index.qmd

@@ -49,27 +49,58 @@ We assume:
    - AR(1) process
    - Differenced AR(1) process
 
-### Simulation model
+### Simulations
 
-We use the generic model structure described above with a renewal process. To simulate noise in the infection process we assume additional Brownian noise for the effective reproduction number of XX.
+#### Observation-generating process
 
-### Simulations
+We use the generic model structure described above with a renewal process as it represents (in its equivalence to an SEIR compartmental model) domain understanding of a model that can capture known infectious dynamics.
+
+We simulate from the renewal process through the following procedure:
+1. Take a fixed timeseries of Rt for 160 days. See the next subsection for more description of these scenarios.
+2. Add noise to the fixed Rt estimates draws from a N(0, 0.1) with a fixed seed of `12345`.
+3. Simulate daily incidence starting from $I_0 = 10$ cases and a fixed generation interval. See the next subsection for more detail.
+4. The delay between infection and case ascertainment is represented as a convolution on the true incidence timeseries, as is standard in the literature **CITATIONS**. For any given infected person the delay between infection and ascertainment is distributed **SOME GAMMA/LOGNORMAL**; this is mapped to our discrete time forward simulations using double interval censoring of both the time of infection and the time of ascertainment **CITE SWP + OTHERS**.
+5. Simulate additional negative binomial observation noise on the delayed cases drawn with mean of the true cases and overdispersion of 10.
+
+We do not add a day-of-week effect.
+
+#### Generation intervals
+
+We use two generation intervals, corresponding to pathogens with long and short GIs. We use descretized, double-censored versions of the GI PMFs.
+1. *Short:* We use a Gamma(shape = 2, scale = 1), corresponding to a pathogen with a relatively short generation interval. Vaguely corresponds to flu A in Wallinga & Lipsitch, 2006
+3. *Medium:* We use a Gamma(shape = 2, scale = 5, corresponding to lots of
+2. *Long:* We use a Gamma(shape = 2, scale = 10), corresponding to a pathogen with a moderately long generation interval (Smallpox? I don't know if we need to ground this in anything real and if we do we could drop this down to 15 days and use varicella?)
 
-We test the following general scenarios:
-- Piecewise constant Rt in an epidemic setting
-      - Generation time:
-- An endemic setting with smoothly varying Rt
-- An outbreak setting with changes in Rt comparable to that observed due to susceptible depletion
-- A mixed outbreak setting with both smooth changes and piecewise changes in Rt
 
-We assume a delay distribution of ** motivated by **.
+We produce the simulations described in the next section for both of these GIs.
 
+#### Scenarios
+
+##### Reproduction number scenarios
+
+We test the following scenario:
+- Piecewise constant Rt
+    - 1.1 for two weeks
+    - 2 for two weeks
+    - 0.5 for two weeks
+    - 1.5 for two weeks
+    - 0.75 for two weeks
+    - 1.1 for six weeks
+    - sine curve centered at 1 with amplitude of 0.3 afterwards
+
+We simulate out of this scenario for the GIs described in the previous section.
+
+This scenario provides both sharp changes at the start of the timeseries and more gradual transitions towards the end. The rolling windows allow for exploration of both of these situations in a single case study. The longer fit to the entire timeseries tests the ability to flexibily handle both of these paradigms in a single fit.
+
+##### Inference scenarios
 We explore the following misspecification scenarios for the generation interval:
 
 - Correct
 - Too short
 - Too long
 
+For each simulated scenario we fit to 12 weeks of data or as much as possible if the scenario is shorter than 12 weeks.
+
 ### Case studies
 
 - [ ] 2014-2016 Sierra Leone Ebola virus disease outbreak