-
Notifications
You must be signed in to change notification settings - Fork 2
/
setup_model.jl
138 lines (109 loc) · 4.2 KB
/
setup_model.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
## UK group sizes
N_uk = 67.22e6
prop_ovr18 = 0.787
prop_men = 0.5
prop_msm = 0.034 #https://www.ons.gov.uk/peoplepopulationandcommunity/culturalidentity/sexuality/bulletins/sexualidentityuk/2020
prop_sexual_active = 1 - 0.154 #A dynamic power-law sexual network model of gonorrhoea outbreaks
N_msm = round(Int64, N_uk * prop_men * prop_ovr18 * prop_msm * prop_sexual_active) #~760k
## Incubation period
# Using this best-fitting distribution, the mean incubation period was estimated
# to be 8.5 days (95% credible intervals (CrI): 6.6–10.9 days),
# with the 5th percentile of 4.2 days and the 95th percentile of 17.3 days (Table 2)
d_incubation = Gamma(6.77, 1 / 0.77)#Fit from rerunning
negbin_std = fill(0.0, 8)
for r = 1:8
p = r / mean(d_incubation)
negbin_std[r] = std(NegativeBinomial(r, p))
end
# plt_incfit = bar(negbin_std,
# title="Discrete time model vs data-driven model for incubation",
# lab="",
# xticks=1:8,
# xlabel="Number of stages",
# ylabel="Std. incubation (days)",
# size=(800, 600), dpi=250,
# legend=:top,
# left_margin=5mm,
# guidefont=16,
# tickfont=13,
# titlefont=17,
# legendfont=12,
# right_margin=5mm)
# hline!(plt_incfit, [std(d_incubation)], lab="std. of data-driven model")
#Optimal choice is 4 stages with effective rate to match the mean
p_incubation = 4 / mean(d_incubation)
α_incubation_eff = -log(1 - p_incubation)
prob_bstfit = [zeros(3); [pdf(NegativeBinomial(4, p_incubation), k) for k = 0:30]]
# bar!(
# plt_incfit,
# prob_bstfit,
# lab="",
# color = :green,
# xlabel = "Incubation period (days)",
# ylabel = "Probability",
# inset = bbox(0.65, 0.15, 0.275, 0.3, :top, :left),
# subplot = 2,
# grid = nothing,
# title = "Best fit incubation distribution",
# yguidefont = 15
# )
# display(plt_incfit)
# savefig(plt_incfit, "plots/incubation_fit.png")
## Set up equipotential groups
n_grp = 10
α_scaling = 0.82
x_min = 1.0
x_max = 3650.0
#Mean rate of yearly contacts over all MSMS
X̄ = (α_scaling / (α_scaling - 1)) * (x_min^(1 - α_scaling) - x_max^(1 - α_scaling)) / (x_min^(-α_scaling) - x_max^(-α_scaling))
#Calculation
C = (x_min^(-α_scaling) - x_max^(-α_scaling)) * X̄ / n_grp * ((1 - α_scaling) / α_scaling)
xs = zeros(n_grp + 1)
xs[1] = x_min
for k = 2:n_grp
xs[k] = (xs[k-1]^(1 - α_scaling) + C)^(1 / (1 - α_scaling))
end
xs[end] = x_max
#Percentages in each group
ps = map(x -> (x_min^(-α_scaling) - x^(-α_scaling)) / (x_min^(-α_scaling) - x_max^(-α_scaling)), xs) |> diff
#Mean daily contact rates within groups
xs_pairs = [(xs[i], xs[i+1]) for i = 1:(length(xs)-1)]
mean_daily_cnts = map(x -> (α_scaling / (α_scaling - 1)) * (x[1]^(1 - α_scaling) - x[2]^(1 - α_scaling)) / (x[1]^(-α_scaling) - x[2]^(-α_scaling)), xs_pairs) .|> x -> x / 365.25
##Plot sexual contact groups
# plt_ps = bar(ps,
# yscale=:log10,
# title="Proportion MSM in each group",
# xticks=1:10,
# ylabel="Proportion",
# xlabel="Sexual activity group",
# lab="")
# plt_μs = bar(mean_daily_cnts,
# yscale=:log10,
# title="Mean daily contact rates in each group",
# xticks=1:10,
# ylabel="Rate (days)",
# xlabel="Sexual activity group",
# lab="")
# hline!(plt_μs, [1 / 31], lab="Vac. threshold", lw=3, legend=:topleft)
# plt = plot(plt_ps, plt_μs,
# size=(1000, 400),
# bottom_margin=5mm, left_margin=5mm)
# display(plt)
# savefig(plt, "plots/sexual_activity_groups.png")
## Set up for ABC
ingroup = 0.99
n_cliques = 50
ts = wks .|> d -> d - Date(2021, 12, 31) .|> t -> t.value
# wkly_vaccinations = [zeros(12); 1000; 2000; fill(5000, 23)] * 1.5
wkly_vaccinations = [[zeros(12); 1000; 2000; fill(5000, 4)] * 1.675
fill(650, 20)
]
wkly_vaccinations_ceased = [copy(wkly_vaccinations)[1:length(wks)+1]; fill(0, 52)]
constants = [N_uk, N_msm, ps, mean_daily_cnts, ingroup, ts, α_incubation_eff, n_cliques, wkly_vaccinations, 0.8, 204] #Constant values passed to the MPX model
## Check model runs
# _p = [0.01, 0.5, 20, 0.2, 0.5, 6, 1.5, 130, 0.7, 0.3, 0.5, 0.5]
# err, pred = MonkeypoxUK.mpx_sim_function_chp(_p, constants, mpxv_wkly[1:9, :])
# plt = plot(pred, color=[1 2])
# scatter!(plt, mpxv_wkly, color=[1 2])
# display(plt)
# print(err)