From 784e10761eab10e75332dbb70ab1ba18c7241131 Mon Sep 17 00:00:00 2001 From: cbernalz Date: Sun, 21 Jul 2024 14:42:04 -0700 Subject: [PATCH] 2024-07-21 update : added main model and added sim func. --- Project.toml | 20 +++ docs/Project.toml | 1 + docs/make.jl | 2 +- docs/src/license.md | 2 +- docs/src/tutorial.md | 51 ++++++- src/UCIWWEIHR.jl | 28 +++- src/bayes_eihr_model.jl | 190 ++++++++++++++++++++++++ src/distribution_functions.jl | 60 ++++++++ src/eihr_ode.jl | 38 +++++ src/generate_simulation_data_ww_eihr.jl | 95 ++++++++++++ test/runtests.jl | 2 +- 11 files changed, 482 insertions(+), 7 deletions(-) create mode 100644 src/bayes_eihr_model.jl create mode 100644 src/distribution_functions.jl create mode 100644 src/eihr_ode.jl create mode 100644 src/generate_simulation_data_ww_eihr.jl diff --git a/Project.toml b/Project.toml index e0b3d01..74c9e1d 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,26 @@ uuid = "176850df-a79a-485c-b76d-f2c16e00fafb" authors = ["Christian O. Bernal Zelaya"] version = "1.0.0-DEV" +[deps] +AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" +NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" +PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" + + [compat] julia = "1" diff --git a/docs/Project.toml b/docs/Project.toml index 6fddde0..68f248b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,4 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" UCIWWEIHR = "176850df-a79a-485c-b76d-f2c16e00fafb" diff --git a/docs/make.jl b/docs/make.jl index 4ddf10f..f453632 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,7 +11,7 @@ makedocs(; sitename = "UCIWWEIHR.jl", format = Documenter.HTML(; prettyurls = get(ENV, "CI", "false") == "true", - canonical = "https://cbernalz.github.io/UCIWWEIHR.jl", + repolink = "https://cbernalz.github.io/UCIWWEIHR.jl", ), pages = [ "HOME" => "index.md", diff --git a/docs/src/license.md b/docs/src/license.md index c49b699..6028fe8 100644 --- a/docs/src/license.md +++ b/docs/src/license.md @@ -1,6 +1,6 @@ # [UCIWWEIHR.jl Package License](@id license) -The UCIWWEIHR.jl package is licensed under the [MIT License](https://opensource.org/licenses/MIT). +The UCIWWEIHR.jl package is licensed under the MIT License. ## MIT License diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md index 3f22e1a..4dc2c1c 100644 --- a/docs/src/tutorial.md +++ b/docs/src/tutorial.md @@ -1,14 +1,59 @@ +```@setup tutorial +using Plots; gr() +Plots.reset_defaults() +``` # [Tutorials](@id tutorial) Welcome to the Tutorials page for the UCIWWEIHR.jl project. This section provides step-by-step guides and examples to help you get started with using the package and understanding its features. -## Getting Started +## 1. Getting Started -### Installation +### 1.1 Installation To install the UCIWWEIHR.jl package, open the Julia REPL and run the following command: -```julia +``` julia using Pkg Pkg.add("UCIWWEIHR") +``` + +## 2. Generating simulated data + +This package provides a way to also simulate data using the model specified in the future paper. The function called `generate_simulation_data_ww_eihr` can be used to generate synthetic data for a given number of samples and features. Here we provide a way to generate synthetic data for the default settings of `generate_simulation_data_ww_eihr`: + +``` @example tutorial +using UCIWWEIHR +using Plots +# Running simulation function with defaults +df = generate_simulation_data_ww_eihr() +first(df, 5) +``` + +Here we can make simple plots to visualize the data generated using the [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package. + +### 2.1. Concentration of Pathogen Genome in WW +```@example tutorial +plot(df.obstimes, df.log_ww_conc, + label=nothing, + xlabel="Obstimes", + ylabel="Conc. of Pathogen Genome in WW", + title="Plot of Conc. of Pathogen Genome in WW Over Time") +``` + +### 2.2. Hospitalizations +```@example tutorial +plot(df.obstimes, df.hosp, + label=nothing, + xlabel="Obstimes", + ylabel="Hosp", + title="Plot of Hosp Over Time") +``` + +### 2.3. Reproductive Number +```@example tutorial +plot(df.obstimes, df.rt, + label=nothing, + xlabel="Obstimes", + ylabel="RT", + title="Plot of RT Over Time") ``` \ No newline at end of file diff --git a/src/UCIWWEIHR.jl b/src/UCIWWEIHR.jl index 5e2711b..26c6c3d 100644 --- a/src/UCIWWEIHR.jl +++ b/src/UCIWWEIHR.jl @@ -1,8 +1,34 @@ module UCIWWEIHR -# Write your package code here. +using AxisArrays +using MCMCChains +using Optim +using LineSearches +using Random +using LinearAlgebra +using NaNMath +using LogExpFunctions +using PreallocationTools +using Distributions +using Turing +using Random +using ForwardDiff +using Logging +using CSV +using DataFrames +using DifferentialEquations +using Plots include("sum_values.jl") +include("generate_simulation_data_ww_eihr.jl") +include("bayes_eihr_model.jl") +include("distribution_functions.jl") +include("eihr_ode.jl") +export eihr_ode export sum_values +export generate_simulation_data_ww_eihr +export bayes_eihr_model +export NegativeBinomial2 +export GeneralizedTDist end \ No newline at end of file diff --git a/src/bayes_eihr_model.jl b/src/bayes_eihr_model.jl new file mode 100644 index 0000000..76003e2 --- /dev/null +++ b/src/bayes_eihr_model.jl @@ -0,0 +1,190 @@ +# bayesian non-constant Rt and hosp rate model eihr +# ------------------------------------------------- +""" + bayes_eihr_model(...) +This is the bayesian semi-parametric model for the wastewater EIHR compartmental model. +The defaults for this fuction will follow those of the default simulation in generate_simulation_data_ww_eihr.jl function. + +# Arguments +- `data_hosp::Array{Float64}`: An array of hospital data. +- `data_wastewater::Array{Float64}`: An array of pathogen genome concentration in localized wastewater data. +- `obstimes::Array{Float64}`: An array of timepoints for observed hosp/wastewater. +- `E_init_sd::Float64`: Standard deviation for the initial number of exposed individuals. +- `E_init_mean::Int64`: Mean for the initial number of exposed individuals. +- `I_init_sd::Float64`: Standard deviation for the initial number of infected individuals. +- `I_init_mean::Int64`: Mean for the initial number of infected individuals. +- `H_init_sd::Float64`: Standard deviation for the initial number of hospitalized individuals. +- `H_init_mean::Int64`: Mean for the initial number of hospitalized individuals. +- `gamma_sd::Float64`: Standard deviation for the rate of incubation. +- `log_gamma_mean::Float64`: Mean for the rate of incubation on log scale. +- `nu_sd::Float64`: Standard deviation for the rate of laeving the infected compartment. +- `log_nu_mean::Float64`: Mean for the rate of leaving the infected compartment on the log scale. +- `epsilon_sd::Float64`: Standard deviation for the rate of hospitalization recovery. +- `log_epsilon_mean::Float64`: Mean for the rate of hospitalization recovery on the log scale. +- `rho_gene_sd::Float64`: Standard deviation for the rho prior. +- `log_rho_gene_mean::Float64`: Mean for the row prior on log scale. +- `tau_sd::Float64`: Standard deviation for the scale/variation of the log scale data. +- `log_tau_mean::Float64`: Mean for the scale/variation of the log scale data on log scale itself. +- `df_shape::Float64`: Shape parameter for the gamma distribution. +- `df_scale::Float64`: Scale parameter for the gamma distribution. +- `sigma_hosp_sd::Float64`: Standard deviation for the negative binomial distribution for hospital data. +- `sigma_hosp_mean::Float64`: Mean for the negative binomial distribution for hospital data. +- `Rt_init_sd::Float64`: Standard deviation for the initial value of the time-varying reproduction number. +- `Rt_init_mean::Float64`: Mean for the initial value of the time-varying reproduction number. +- `sigma_Rt_sd::Float64`: Standard deviation for normal prior of log time-varying reproduction number standard deviation. +- `sigma_Rt_mean::Float64`: Mean for normal prior of log time-varying reproduction number standard deviation. +- `w_init_sd::Float64`: Standard deviation for the initial value of the time-varying hospitalization rate. +- `w_init_mean::Float64`: Mean for the initial value of the time-varying hospitalization rate. +- `sigma_w_sd::Float64`: Standard deviation for normal prior of log time-varying hospitalization rate standard devaiation. +- `sigma_w_mean::Float64`: Mean for normal prior of time-varying hospitalization rate standard devaiation. +- `param_change_times::Array{Float64}`: An array of timepoints where the parameters change. + +""" + +@model function bayes_eihr_model( + data_hosp::Array{Float64}, + data_wastewater::Array{Float64}, + obstimes::Array{Float64}, + E_init_sd::Float64, E_init_mean::Int64, + I_init_sd::Float64, I_init_mean::Int64, + H_init_sd::Float64, H_init_mean::Int64, + gamma_sd::Float64, log_gamma_mean::Float64, + nu_sd::Float64, log_nu_mean::Float64, + epsilon_sd::Float64, log_epsilon_mean::Float64, + rho_gene_sd::Float64, log_rho_gene_mean::Float64, + tau_sd::Float64, log_tau_mean::Float64, + df_shape::Float64, df_scale::Float64, + sigma_hosp_sd::Float64, sigma_hosp_mean::Float64, + Rt_init_sd::Float64, Rt_init_mean::Float64, + sigma_Rt_sd::Float64, sigma_Rt_mean::Float64, + w_init_sd::Float64, w_init_mean::Float64, + sigma_w_sd::Float64, sigma_w_mean::Float64, + param_change_times::Array{Float64} + ) + + + # Prelims + max_neg_bin_sigma = 1e10 + min_neg_bin_sigma = 1e-10 + + + # Calculate number of observed datapoints timepoints + l_obs = length(obstimes) + l_param_change_times = length(param_change_times) + + + # PRIORS----------------------------- + # Compartments + E_init_non_centered ~ Normal() + I_init_non_centered ~ Normal() + H_init_non_centered ~ Normal() + # Parameters for compartments + gamma_non_centered ~ Normal() + nu_non_centered ~ Normal() + epsilon_non_centered ~ Normal() + # Parameters for wastewater + rho_gene_non_centered ~ Normal() # gene detection rate + tau_non_centered ~ Normal() # standard deviation for log scale data + df ~ Gamma(df_shape, df_scale) + # Parameters for hospital + sigma_hosp_non_centered ~ Normal() + # Non-constant Rt + Rt_params_non_centered ~ MvNormal(zeros(l_param_change_times + 2), I) # +2 for sigma and init + sigma_Rt_non_centered = Rt_params_non_centered[1] + Rt_init_non_centered = Rt_params_non_centered[2] + log_Rt_steps_non_centered = Rt_params_non_centered[3:end] + # Non-constant Hosp Rate w + w_params_non_centered ~ MvNormal(zeros(l_param_change_times + 2), I) # +2 for sigma and init + sigma_w_non_centered = w_params_non_centered[1] + w_init_non_centered = w_params_non_centered[2] + log_w_steps_non_centered = w_params_non_centered[3:end] + + + # TRANSFORMATIONS-------------------- + # Compartments + E_init = E_init_non_centered * E_init_sd + E_init_mean + I_init = I_init_non_centered * I_init_sd + I_init_mean + H_init = H_init_non_centered * H_init_sd + H_init_mean + # Parameters for compartments + gamma = exp(gamma_non_centered * gamma_sd + log_gamma_mean) + nu = exp(nu_non_centered * nu_sd + log_nu_mean) + epsilon = exp(epsilon_non_centered * epsilon_sd + log_epsilon_mean) + # Parameters for wastewater + rho_gene = exp(rho_gene_non_centered * rho_gene_sd + log_rho_gene_mean) + tau = exp(tau_non_centered * tau_sd + log_tau_mean) + # Parameters for hospital + sigma_hosp = clamp.(sigma_hosp_non_centered * sigma_hosp_sd + sigma_hosp_mean, min_neg_bin_sigma, max_neg_bin_sigma) + # Non-constant Rt + Rt_init = exp(Rt_init_non_centered * Rt_init_sd + Rt_init_mean) + sigma_Rt = exp(sigma_Rt_non_centered * sigma_Rt_sd + sigma_Rt_mean) + alpha_t_no_init = exp.(log(Rt_init) .+ cumsum(log_Rt_steps_non_centered) * sigma_Rt) * nu + alpha_init = Rt_init * nu + alpha_t = vcat(alpha_init, alpha_t_no_init) + # Non-constant Hosp Rate w + w_init = exp(w_init_non_centered * w_init_sd + w_init_mean) + sigma_w = exp(sigma_w_non_centered * sigma_w_sd + sigma_w_mean) + w_no_init = exp.(log(w_init) .+ cumsum(log_w_steps_non_centered) * sigma_w) + w_t = vcat(w_init, w_no_init) + + + # ODE SETUP-------------------------- + prob = ODEProblem{true}(eihr_ode!, zeros(3), (0.0, obstimes[end]), ones(5)) + function param_affect_beta!(integrator) + ind_t = searchsortedfirst(param_change_times, integrator.t) # Find the index of param_change_times that contains the current timestep + integrator.p[1] = alpha_t_no_init[ind_t] # Replace alpha with a new value from alpha_t_no_init + integrator.p[4] = w_no_init[ind_t] # Replace w with a new value from w_no_init + end + param_callback = PresetTimeCallback(param_change_times, param_affect_beta!, save_positions=(false, false)) + u0 = [E_init, I_init, H_init] + p0 = [alpha_init, gamma, nu, w_init, epsilon] + extra_ode_precision = false + abstol = extra_ode_precision ? 1e-11 : 1e-9 + reltol = extra_ode_precision ? 1e-8 : 1e-6 + sol = solve(prob, Tsit5(); callback=param_callback, saveat=obstimes, save_start=true, + verbose=false, abstol=abstol, reltol=reltol, u0=u0, p=p0, tspan=(0.0, obstimes[end])) + # If the ODE solver fails, reject the sample by adding -Inf to the likelihood + if sol.retcode != :Success + println("An error occurred during ODE solution!!!") + Turing.@addlogprob! -Inf + return + end + sol_array = Array(sol) + I_comp_sol = clamp.(sol_array[2,2:end],1, 1e10) + + + # W-W means-------------------------- + # E - 1 // I - 2 // H - 3 // R - 4 + log_genes_mean = log.(I_comp_sol) .+ log(rho_gene) # first entry is the initial conditions, we want 2:end + # Likelihood calculations------------ + sol_hosp = clamp.(sol_array[3,2:end], 1, 1e10) + for i in 1:l_obs + data_wastewater[i] ~ GeneralizedTDist(log_genes_mean[i], tau, df) + data_hosp[i] ~ NegativeBinomial2(sol_hosp[i], sigma_hosp) + end + + + # Generated quantities + H_comp = sol_array[3, :] + rt_vals = alpha_t / nu + w_t = w_t + + return ( + E_init, + I_init, + H_init, + alpha_t = alpha_t, + gamma = gamma, + nu = nu, + w_t = w_t, + epsilon = epsilon, + rt_vals = rt_vals, + rho_gene = rho_gene, + tau = tau, + df = df, + sigma_hosp = sigma_hosp, + H = H_comp, + log_genes_mean = log_genes_mean + ) + + + end diff --git a/src/distribution_functions.jl b/src/distribution_functions.jl new file mode 100644 index 0000000..bb1388a --- /dev/null +++ b/src/distribution_functions.jl @@ -0,0 +1,60 @@ + +# credit to Damon Bayer for all of these helper functions +""" +Create a re-parametrized negative binomial distribution in terms of mean and overdispersion. + +# Arguments +- `μ`: Mean of the distribution. +- `ϕ`: Overdispersion parameter. + +# Returns +A `Distributions.NegativeBinomial` distribution object. + +""" +function NegativeBinomial2(μ, ϕ) + p = 1 / (1 + μ / ϕ) + + if p <= zero(p) + p = eps(zero(p)) + end + + if p >= one(p) + p = prevfloat(one(p)) + end + + r = ϕ + + if r <= zero(r) + r = eps(zero(r)) + end + + Distributions.NegativeBinomial(r, p) +end + +""" +Create a generalized t distribution with location and scale. + +# Arguments +- `μ`: Location parameter. +- `σ`: Scale parameter. +- `ν`: Degrees of freedom. + +# Returns +A generalized t distribution object. + +""" + +function GeneralizedTDist(μ, σ, ν) + if ν <= zero(ν) + ν = nextfloat(zero(ν)) + end + + if σ <= zero(σ) + σ = nextfloat(zero(σ)) + end + + Generalized_T = μ + σ*Distributions.TDist(ν) + + return(Generalized_T) + +end diff --git a/src/eihr_ode.jl b/src/eihr_ode.jl new file mode 100644 index 0000000..cf3d89b --- /dev/null +++ b/src/eihr_ode.jl @@ -0,0 +1,38 @@ +""" +eihr_ode!(du, u, p, t) + +Calculate the ordinary differential equations (ODEs) for the EIHR model. + +Parameters: +- `du`: Array{Float64,1} - The derivative of the state variables. +- `u`: Array{Float64,1} - The current state variables. +- `p`: Tuple{Float64,Float64,Float64,Float64,Float64} - The model parameters (alpha, gamma, nu, w, epsilon). +- `t`: Float64 - The current time. + +Returns: +- `du`: Array{Float64,1} - The derivative of the state variables. + +""" +function eihr_ode!(du, u, p, t) + # ODE for EIHR model with constant beta + (E, I, H) = u + (alpha, gamma, nu, w, epsilon) = p + # w - ratio of infected individuals going to hospital + + # -> E + exposed_in = alpha * I + # E -> I + progression = gamma * E + # I -> H + hospitalization = nu * w * I + # I -> R + non_hospitalized_recovery = nu * (1 - w) * I + # H -> R + hospitalized_recovery = epsilon * H + + @inbounds begin + du[1] = exposed_in - progression # E + du[2] = progression - (hospitalization + non_hospitalized_recovery) # I + du[3] = hospitalization - hospitalized_recovery # H + end +end \ No newline at end of file diff --git a/src/generate_simulation_data_ww_eihr.jl b/src/generate_simulation_data_ww_eihr.jl new file mode 100644 index 0000000..0b9a19a --- /dev/null +++ b/src/generate_simulation_data_ww_eihr.jl @@ -0,0 +1,95 @@ +""" +## Generating Simulation Data for Wastewater EIHR + +To generate simulation data for wastewater EIHR, you can use the `generate_simulation_data_ww_eihr` function defined in the `UCIWWEIHR.jl` package. This function allows you to customize various parameters for the simulation. + +### Function Signature + +```julia + +# Arguments +- time_points::Int64: Number of time points wanted for simulation. Default value is 150. +- seed::Int64: Seed for random number generation. Default value is 1. +- E_init::Int64: Initial number of exposed individuals. Default value is 200. +- I_init::Int64: Initial number of infected individuals. Default value is 100. +- H_init::Int64: Initial number of hospitalized individuals. Default value is 20. +- gamma::Float64: Rate of incubation. Default value is 1/4. +- nu::Float64: Rate of leaving the infected compartment. Default value is 1/7. +- epsilon::Float64: Rate of hospitalization recovery. Default value is 1/5. +- rho_gene::Float64: Contribution of infected individual's pathogen genome into wastewater. Default value is 0.011. +- tau::Float64: Scale/variation of the log concentration of pathogen genome in wastewater. Default value is 0.1. +- df::Float64: Degrees of freedom for generalized t distribution for log concentration of pathogen genome in wastewater. Default value is 29. +- sigma_hosp::Float64: Standard deviation for the negative binomial distribution for hospital data. Default value is 800. +- Rt_init::Float64: Initial value of the time-varying reproduction number. Default value is 1. +- sigma_Rt::Float64: Standard deviation for random walk of time-varying reproduction number. Default value is sqrt(0.02). +- w_init::Float64: Initial value of the time-varying hospitalization rate. Default value is 0.35. +- sigma_w::Float64: Standard deviation for random walk of time-varying hospitalization rate. Default value is sqrt(0.02). +""" +function generate_simulation_data_ww_eihr( + time_points::Int64=150, seed::Int64=1, + E_init::Int64=200, I_init::Int64=100, H_init::Int64=20, + gamma::Float64=1/4, nu::Float64=1/7, epsilon::Float64=1/5, + rho_gene::Float64=0.011, tau::Float64=0.1, df::Float64=29.0, + sigma_hosp::Float64=800.0, + Rt_init::Float64=1.0, sigma_Rt::Float64=sqrt(0.001), + w_init::Float64=0.35, sigma_w::Float64=sqrt(0.001), +) + + Random.seed!(seed) + + # Rt and W SETUP-------------------------- + Rt_t_no_init = Float64[] # Pre-defined vector + w_no_init = Float64[] # Pre-defined vector + log_Rt_t = log(Rt_init) + log_w_t = log(w_init) + for _ in 1:time_points + log_Rt_t = log_Rt_t + rand(Normal(0, sigma_Rt)) + log_w_t = log_w_t + rand(Normal(0, sigma_w)) + push!(Rt_t_no_init, exp(log_Rt_t)) + push!(w_no_init, exp(log_w_t)) + end + alpha_t_no_init = Rt_t_no_init * nu + alpha_init = Rt_init * nu + + # ODE SETUP-------------------------- + prob = ODEProblem{true}(eihr_ode!, zeros(3), (0.0, time_points), ones(5)) + function param_affect_beta!(integrator) + ind_t = searchsortedfirst(collect(1:time_points), integrator.t) # Find the index of collect(1:time_points) that contains the current timestep + integrator.p[1] = alpha_t_no_init[ind_t] # Replace alpha with a new value from alpha_t_no_init + integrator.p[4] = w_no_init[ind_t] # Replace w with a new value from w_no_init + end + param_callback = PresetTimeCallback(collect(1:time_points), param_affect_beta!, save_positions=(false, false)) + u0 = [E_init, I_init, H_init] + p0 = [alpha_init, gamma, nu, w_init, epsilon] + extra_ode_precision = false + abstol = extra_ode_precision ? 1e-11 : 1e-9 + reltol = extra_ode_precision ? 1e-8 : 1e-6 + sol = solve(prob, Tsit5(); callback=param_callback, saveat=collect(1:time_points), save_start=true, + verbose=false, abstol=abstol, reltol=reltol, u0=u0, p=p0, tspan=(0.0, time_points)) + if sol.retcode != :Success + Turing.@addlogprob! -Inf + return + end + sol_array = Array(sol) + I_comp_sol = clamp.(sol_array[2,2:end],1, 1e10) + H_comp_sol = clamp.(sol_array[3,2:end], 1, 1e10) + + # Log Gene SETUP-------------------------- + log_genes_mean = log.(I_comp_sol) .+ log(rho_gene) # first entry is the initial conditions, we want 2:end + data_wastewater = zeros(time_points) + data_hosp = zeros(time_points) + for t_i in 1:time_points + data_wastewater[t_i] = rand(GeneralizedTDist(log_genes_mean[t_i], tau, df)) + data_hosp[t_i] = rand(NegativeBinomial2(H_comp_sol[t_i], sigma_hosp)) + end + + df = DataFrame( + obstimes = 1:time_points, + log_ww_conc = data_wastewater, + hosp = data_hosp, + rt = Rt_t_no_init + ); + return df + + +end diff --git a/test/runtests.jl b/test/runtests.jl index 741a97c..509dade 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,4 +7,4 @@ using Test y = 2 @test UCIWWEIHR.sum_values(x,y) == 4 -end \ No newline at end of file +end