Skip to content
Snippets Groups Projects
Commit 04fdb04e authored by Peter Jentsch's avatar Peter Jentsch
Browse files

initial model fits done, bounded rationality too low or immunity loss too high

parent af952914
No related branches found
No related tags found
No related merge requests found
Showing
with 192 additions and 205 deletions
......@@ -38,6 +38,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
......
......@@ -2,65 +2,4 @@ using CovidAlertVaccinationModel
using OnlineStats
using Plots
using StatsBase
const parameters = (
sim_length = 600,
num_households = 5000,
I_0_fraction = 0.002,
base_transmission_probability = 0.001,
recovery_rate = 1/7,
immunization_loss_prob = 0.0055, #mean time of 6 months
π_base = -0.3,
η = 0.0,
κ = 0.8,
ω = 0.00005,
ρ_y = 0.0,
ρ_m = 0.0,
ρ_o = 0.0,
ω_en = 0.00,
ρ_en = [0.0,0.0,0.0],
γ = 0.0,
β = 10.0,
notification_parameter = 0.001,
vaccinator_prob = 0.2,
app_user_fraction = 0.0,
notification_threshold = 17,
immunizing = true,
immunization_delay = 14,
immunization_begin_day = 30,
infection_introduction_day = 100
)
# seasonal_transmission_dist = CovidAlertVaccinationModel.fit_epi_parameters(parameters,0.073) ##seasonal
# outbreak_transmission_dist = CovidAlertVaccinationModel.fit_epi_parameters(parameters,0.241) ##outbreak
# function plot_max_posterior(fname,particles)
# samples = 5
# base_transmission = mode(particles.P.particles)
# p_tuple_without_vac = merge(parameters,
# (
# sim_length = 150,
# immunization_begin_day = 0,
# infection_introduction_day = 1,
# immunizing = false,
# )
# )
# new_params = merge(p_tuple_without_vac, (base_transmission_probability = base_transmission,))
# out = mean_solve(samples, new_params ,DebugRecorder)
# p = plot_model(nothing,[nothing],[out],new_params.infection_introduction_day,new_params.immunization_begin_day)
# savefig(p,"$fname.pdf")
# hist = fit(Histogram,particles.P.particles; nbins = 25)
# p = plot(hist;legend = false)
# savefig(p,"$(fname)_posterior.pdf")
# end
# plot_max_posterior("seasonal", seasonal_transmission_dist)
# plot_max_posterior("outbreak", outbreak_transmission_dist)
a = CovidAlertVaccinationModel.fit_behavioural_parameters(parameters,nothing)
function solve_and_plot_parameters()
out = mean_solve(samples, parameters ,DebugRecorder)
p = plot_model(nothing,[nothing],[out],parameters.infection_introduction_day,parameters.immunization_begin_day)
savefig(p,"timeseries.pdf")
return out
end
\ No newline at end of file
CovidAlertVaccinationModel.fit_parameters(CovidAlertVaccinationModel.get_parameters())
\ No newline at end of file
File added
File added
File added
File added
......@@ -2,41 +2,13 @@ using CovidAlertVaccinationModel
using OnlineStats
using Plots
const samples = 1
const parameters = (
sim_length = 600,
num_households = 5000,
I_0_fraction = 0.0,
base_transmission_probability = 0.001,
recovery_rate = 1/7,
immunization_loss_prob = 0.0055, #mean time of 6 months
π_base = -0.8,
η = 0.0,
κ = 0.0,
ω = 0.00005,
ρ_y = 0.0,
ρ_m = 0.0,
ρ_o = 0.0,
ω_en = 0.00,
ρ_en = [0.0,0.0,0.0],
γ = 0.0,
β = 10.0,
notification_parameter = 0.001,
vaccinator_prob = 0.2,
app_user_fraction = 0.2,
notification_threshold = 17,
immunizing = true,
immunization_delay = 14,
immunization_begin_day = 30,
infection_introduction_day = 100
)
function solve_and_plot_parameters()
out = mean_solve(samples, parameters ,DebugRecorder)
p = plot_model(nothing,[nothing],[out],parameters.infection_introduction_day,parameters.immunization_begin_day)
p = CovidAlertVaccinationModel.get_parameters()
out = mean_solve(samples, p,DebugRecorder)
p = plot_model(nothing,[nothing],[out],p.infection_introduction_day,p.immunization_begin_day)
savefig(p,"timeseries.pdf")
return out
end
out = solve_and_plot_parameters()
println(mean.(out.new_ymo_immunization.Y))
println(mean.(out.new_ymo_immunization.M))
println(mean.(out.new_ymo_immunization.O))
\ No newline at end of file
println(sum.(eachrow(mean.(out.new_ymo_immunization))))
\ No newline at end of file
function get_parameters()
params = (
sim_length = 600,
sim_length = 120,
num_households = 5000,
I_0_fraction = 0.002,
base_transmission_probability = 0.001,
recovery_rate = 1/7,
immunization_loss_prob = 0.0055, #mean time of 6 months
π_base = -0.3,
immunization_loss_prob = 0.0001, #mean time of 6 months
π_base_y = -5.0,
π_base_m = -5.0,
π_base_o = 1.0,
η = 0.0,
κ = 0.8,
κ = 0.0,
ω = 0.00005,
ρ_y = 0.0,
ρ_m = 0.0,
ρ_o = 0.0,
ω_en = 0.00,
ρ_en = [0.0,0.0,0.0],
ρ_en = [0.0,0.0,0.0],
γ = 0.0,
β = 10.0,
β = 5.0,
notification_parameter = 0.001,
vaccinator_prob = 0.2,
app_user_fraction = 0.2,
app_user_fraction = 0.0,
notification_threshold = 2,
immunizing = true,
immunization_delay = 14,
immunization_begin_day = 30,
infection_introduction_day = 100
immunization_begin_day = 60,
infection_introduction_day = 120
)
return params
end
......
......@@ -108,10 +108,29 @@ function OnlineStats.fit!(stat_type::R,pt) where {T,OS<:OnlineStat{T}, R<:Abstra
end
end
end
using Printf
const ts_colors = cgrad(:PuBu_9)
function plot_model(varname,univariate_series, output_list::Vector{T},infection_begin,vac_begin) where T<:DebugRecorder
plts = [plot() for i=1:7]
sim_length = length(output_list[1].recorded_status_totals.S)
ts_list(data) = [
(data.recorded_status_totals.S, "Susceptible over time"),
(data.recorded_status_totals.I, "Infected over time"),
(data.recorded_status_totals.R, "Recovered over time"),
(data.recorded_status_totals.V, "Vaccinated over time"),
(data.total_vaccinators, "No. vaccinators over time"),
(data.mean_time_since_last_notification, "Mean time since last notification"),
(data.daily_cases,"Daily (incident) cases"),
(data.new_ymo_immunization.Y, "new Y vaccinations each day"),
(data.new_ymo_immunization.M, "new M vaccinations each day"),
(data.new_ymo_immunization.O, "new O vaccinations each day"),
]
l = length(ts_list(output_list[1]))
plts = [plot() for i=1:l]
for (i,(p,data)) in enumerate(zip(univariate_series, output_list))
# display(p[varname])
if !isnothing(varname)
......@@ -120,65 +139,19 @@ function plot_model(varname,univariate_series, output_list::Vector{T},infection_
else
labelname = nothing
end
# display(p_val)
plot!(plts[1], mean.(data.recorded_status_totals.S); ribbon = std.(data.recorded_status_totals.S),
label = labelname, line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "Susceptible over time",
legend=:topright
)
plot!(plts[2], mean.(data.recorded_status_totals.I); ribbon = std.(data.recorded_status_totals.I),
label = labelname, line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "Infected over time",
legend=false,
)
plot!(plts[3], mean.(data.recorded_status_totals.R); ribbon = std.(data.recorded_status_totals.R),
label = labelname, line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "Recovered over time",
legend=false
)
plot!(plts[4], mean.(data.recorded_status_totals.V); ribbon = std.(data.recorded_status_totals.V),
label = labelname, line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "Vaccinated over time",
legend=false
)
plot!(plts[5], mean.(data.total_vaccinators); ribbon = std.(data.total_vaccinators),
label = labelname, line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = "No. vaccinators over time",
legend=false
)
plot!(plts[6], mean.(data.mean_time_since_last_notification); ribbon = std.(data.mean_time_since_last_notification),
label = labelname, line_z = i, color=:blues,
ylabel = "Days",
colorbar = false,
title = "Mean time since last notification",
legend=false
)
# ic,exp_growth_rate = exp_growth_rate_estimation(mean.(data.daily_cases[infection_begin,:]))
plot!(plts[7], mean.(data.daily_cases); ribbon = std.(data.daily_cases),
label = labelname, line_z = i, color=:blues,
ylabel = "Days",
colorbar = false,
title = "Daily (incident) cases",
legend= false
)
for (plt,(ts,title)) in zip(plts,ts_list(data))
plot!(plt, mean.(ts); ribbon = std.(ts),
label = labelname, line_z = i, color=:blues,
ylabel = "no. of people",
colorbar = false,
title = title,
)
end
end
plot!(plts[begin]; legends = :topright)
for p in plts
vline!(p,[infection_begin]; label = "infection begin", line =:dot)
vline!(p,[vac_begin]; label = "vaccination begin",line = :dot)
end
return plot(plts...;layout = (7,1),size=(800,2900),leftmargin = 8Plots.mm)
return plot(plts...;layout = (l,1),size=(800,400*l),leftmargin = 12Plots.mm)
end
......@@ -2,7 +2,16 @@ using KissABC
using CurveFit
using Loess
const target_cumulative_vac_proportion = 0.33
const vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
const ymo_vac = @SVector [0.255,0.278,0.602]
function solve_w_parameters(default_p, p_names, new_p_list)
new_params = merge(default_p, NamedTuple{p_names}(ntuple(i -> new_p_list[i],length(p_names))))
out = DebugRecorder(0,default_p.sim_length)
model = abm(new_params,out)
return out, model
end
@views function exp_growth_rate_estimation(incident_cases)
diff(l) = [l[i] - l[i-1] for i=2:length(l)]
l = length(incident_cases)
......@@ -14,13 +23,6 @@ using Loess
return ic,exp_growth_rate
end
function vaccination_data()
ymo_total_vaccinated = [0.279,0.38,0.7]
total_vaccinated = 0.33
end
function fit_epi_parameters(p_tuple, target_growth_rate)
p_tuple_without_vac = merge(p_tuple,
(
......@@ -30,71 +32,99 @@ function fit_epi_parameters(p_tuple, target_growth_rate)
immunizing = false,
)
)
samples = 1
p_names = (:base_transmission_probability,)
priors = Factored(Uniform(0.0001,0.005))
function cost(p)
new_params = merge(p_tuple_without_vac, (base_transmission_probability = p[1],))
out = mean_solve(samples, new_params ,DebugRecorder)
incident_cases = mean.(out.daily_cases)
output,model = solve_w_parameters(p_tuple_without_vac,p_names ,p)
incident_cases = output.daily_cases
_,exp_growth_rate = exp_growth_rate_estimation(incident_cases)
return (target_growth_rate - exp_growth_rate)^2
end
out = ABCDE(priors,cost,0.01; verbose=true, nparticles=100,generations=20, parallel = true) #this one has better NaN handling
return out
out = ABCDE(priors,cost,0.01; verbose=true, nparticles= 200,generations=20, parallel = true) #this one has better NaN handling
return NamedTuple{p_names}((out.P.particles,))
end
function fit_behavioural_parameters(p_tuple, target_growth_rate)
function fit_pre_inf_behavioural_parameters(p_tuple)
samples = 1
p_names = (:π_base,:ρ_y,:ρ_m,:ρ_o)
priors = Factored(Uniform(-2.0,0.0),Uniform(0.0,1.0),Uniform(0.0,1.0),Uniform(0.0,1.0))
#simulation begins in august
#30 days for opinion dynamics to stabilize, then immunization begins in september,
#infection is introduced at the end of november
sim_length = 120
p_names = (:π_base_y,:π_base_m,:π_base_o)
priors = Factored(
Uniform(-10.0,2.0),
Uniform(-10.0,2.0),
Uniform(-10.0,2.0)
)
#simulation begins in july
#60 days for opinion dynamics to stabilize, then immunization begins in september,
#infection is not considered
sim_length = 180
p_tuple_adjust = merge(p_tuple,
(
sim_length = sim_length,
I_0_fraction = 0.000,
immunization_begin_day =60,
infection_introduction_day = 91,
infection_introduction_day = 180,
immunizing = true,
)
)
target_cumulative_vac_proportion = 0.33
vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
ymo_vac = @SVector [0.255,0.278,0.602]
function cost(p)
new_params = merge(p_tuple_adjust, NamedTuple{p_names}(ntuple(i -> p[i],length(p_names))))
out = DebugRecorder(0,sim_length)
model = abm(new_params,out)
function cost(p)
output,model = solve_w_parameters(p_tuple_adjust, p_names,p)
target_cumulative_vaccinations = target_cumulative_vac_proportion*model.nodes
target_ymo_vac = ymo_vac .* sum(vaccination_data[1:4]) .* target_cumulative_vaccinations
ymo_vaccination_ts = out.new_ymo_immunization
ymo_vaccination_ts = output.new_ymo_immunization
total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
# display(target_ymo_vac)
# display(total_preinfection_vaccination)
return sum((total_preinfection_vaccination .- target_ymo_vac).^2)
end
#smc(priors,cost; verbose = true, nparticles = 100, parallel = true)#
out = ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=100, parallel = true) #this one has better NaN handling
return out
#smc(priors,cost; verbose = true, nparticles = 100, parallel = true)# smc(priors,cost; verbose = true, nparticles =300, parallel = true)#
out = smc(priors,cost; verbose = true, nparticles = 800, parallel = true)#
return NamedTuple{p_names}(ntuple(i -> out.P[i].particles,length(p_names)))
end
function fit_post_inf_behavioural_parameters(p_tuple)
samples = 1
p_names = (:ω,)
priors = Factored(Uniform(0.0,0.1))
#simulation begins in july
#60 days for opinion dynamics to stabilize, then immunization begins in september,
#infection introduced at beginning of december
sim_length = 300
p_tuple_adjust = merge(p_tuple,
(
sim_length = sim_length,
I_0_fraction = 0.002,
immunization_begin_day =60,
infection_introduction_day = 180,
immunizing = true,
)
)
target_cumulative_vac_proportion = 0.33
vaccination_data = @SVector [0.0,0.043,0.385,0.424,0.115,0.03,0.005] #by month starting in august
ymo_vac = @SVector [0.255,0.278,0.602]
function cost(p)
output,model = solve_w_parameters(p_tuple_adjust, p_names,p)
target_cumulative_vaccinations = target_cumulative_vac_proportion*model.nodes
target_ymo_vac = ymo_vac .* sum(vaccination_data[5:end]) .* target_cumulative_vaccinations
ymo_vaccination_ts = output.new_ymo_immunization
total_vaccination = sum.(eachrow(ymo_vaccination_ts[:,180:end]))
return sum((total_vaccination .- target_ymo_vac).^2)
end
out =smc(priors,cost; verbose = true, nparticles = 200, parallel = true)# ABCDE(priors,cost,1e6; verbose=true, nparticles=300,generations=1000, parallel = true) #this one has better NaN handling
return NamedTuple{p_names}((out.P.particles,))
end
function plot_behavioural_fit(particles,p_tuple)
p_names = (:π_base, :ρ_y,:ρ_m,:ρ_o)
p_names = (:π_base_y,:π_base_m,:π_base_o)
sim_length = 210
samples = 1
p_tuple_adjust = merge(p_tuple,
(
(
sim_length = sim_length,
I_0_fraction = 0.000,
immunization_begin_day =60,
infection_introduction_day = 91,
immunizing = true,
)
)
......@@ -108,10 +138,84 @@ function plot_behavioural_fit(particles,p_tuple)
ymo_vaccination_ts = mean.(out.new_ymo_immunization)
total_preinfection_vaccination = sum.(eachrow(ymo_vaccination_ts))
display(total_preinfection_vaccination)
p = [plot(),plot(),plot()]
p = plot_model(nothing,[nothing],[out],new_params.infection_introduction_day,new_params.immunization_begin_day)
savefig(p,"behaviour_fit.pdf")
return out
end
# outbreak_transmission_dist = CovidAlertVaccinationModel.fit_epi_parameters(default_parameters,0.241) ##outbreak
# serialize(joinpath(PACKAGE_FOLDER,"abm_parameter_fits","outbreak_inf_parameters.dat"),outbreak_transmission_dist)
# plot_max_posterior("outbreak", outbreak_transmission_dist,default_parameters)
function fit_parameters(default_parameters)
seasonal_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","outbreak_inf_parameters.dat")
pre_inf_behaviour_parameters_path =joinpath(PACKAGE_FOLDER,"abm_parameter_fits","pre_inf_behaviour_parameters.dat")
post_inf_behaviour_parameters_path = joinpath(PACKAGE_FOLDER,"abm_parameter_fits","post_inf_behaviour_parameters.dat")
seasonal_transmission_dist = CovidAlertVaccinationModel.fit_epi_parameters(default_parameters,0.073) ##seasonal
plot_max_posterior("seasonal", seasonal_transmission_dist,default_parameters)
pre_inf_behaviour_parameters = CovidAlertVaccinationModel.fit_pre_inf_behavioural_parameters(default_parameters)
fitted_parameters = map(mode,((;seasonal_transmission_dist...,pre_inf_behaviour_parameters...)))
fitted_parameter_tuple = merge(default_parameters,fitted_parameters)
serialize(seasonal_parameters_path,seasonal_transmission_dist)
serialize(pre_inf_behaviour_parameters_path, pre_inf_behaviour_parameters)
post_inf_behaviour_parameters = fit_post_inf_behavioural_parameters(fitted_parameter_tuple)
serialize(post_inf_behaviour_parameters_path, post_inf_behaviour_parameters)
fitted_parameters_with_post_inf_behaviour = (;
deserialize(seasonal_parameters_path)...,
deserialize(pre_inf_behaviour_parameters_path)...,
deserialize(post_inf_behaviour_parameters_path)...
)
display(map(mode,fitted_parameters_with_post_inf_behaviour))
return plot_fitting_posteriors("post_inf_fitting",fitted_parameters_with_post_inf_behaviour,default_parameters)
end
function plot_max_posterior(fname,particles,parameters)
samples = 5
base_transmission = mode(particles.base_transmission_probability)
p_tuple_without_vac = merge(parameters,
(
sim_length = 150,
immunization_begin_day = 0,
infection_introduction_day = 1,
immunizing = false,
)
)
new_params = merge(p_tuple_without_vac, (base_transmission_probability = base_transmission,))
out = mean_solve(samples, new_params ,DebugRecorder)
p = plot_model(nothing,[nothing],[out],new_params.infection_introduction_day,new_params.immunization_begin_day)
savefig(p,"$fname.pdf")
hist = StatsBase.fit(Histogram,particles.base_transmission_probability; nbins = 25)
p = plot(hist;legend = false)
savefig(p,"$(fname)_posterior.pdf")
end
#once we have vaccine parameters, with seasonal, look at cumulative infections
#
function plot_fitting_posteriors(fname,particles_tuple,parameters)
p_tuple_adjust = merge(parameters,
(
sim_length = 500,
I_0_fraction = 0.002,
immunization_begin_day =60,
infection_introduction_day = 180,
immunizing = true,
)
)
out = mean_solve(5, merge(p_tuple_adjust,map(mode,particles_tuple)) ,DebugRecorder)
p = plot_model(nothing,[nothing],[out],parameters.infection_introduction_day,parameters.immunization_begin_day)
savefig(p, "$fname.pdf")
plts = [plot() for i in 1:length(particles_tuple)]
for (plt,(k,v)) in zip(plts,pairs(particles_tuple))
hist = StatsBase.fit(Histogram,v; nbins = 40)
plot!(plt,hist;legend = false,xlabel = k)
end
p = plot(plts...; size = (1000,600))
savefig(p,"$(fname)_posteriors.pdf")
return out
end
......@@ -146,11 +146,9 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
vac_payoff = 0
num_soc_nbrs = 0
random_soc_network = sample(Random.default_rng(Threads.threadid()), soc_network.graph_list[t])
# display(random_soc_network)
if !isempty(neighbors(random_soc_network.g,i))
random_neighbour = sample(Random.default_rng(Threads.threadid()), neighbors(random_soc_network.g,i))
if u_vac[random_neighbour] == u_vac[i]
vac_payoff += π_base[Int(demographics[i])] + total_infections*ω
if app_user[i] && time_of_last_alert[app_user_list[i]]>=0
......@@ -169,7 +167,9 @@ Base.@propagate_inbounds @views function update_vaccination_opinion_state!(t,mod
end
end
end
modelsol.daily_vaccinators = count(==(true),u_vac) #could maybe make this more efficient
end
......
......@@ -24,7 +24,6 @@ using KissABC
using CSV
import Pandas: read_csv
using DataFrames
using VectorizedRNG
using StaticArrays
import LightGraphs.neighbors
......
No preview for this file type
File added
No preview for this file type
No preview for this file type
File added
File added
No preview for this file type
No preview for this file type
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment