-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathplot_comparePopSizes.py
148 lines (119 loc) · 7.94 KB
/
plot_comparePopSizes.py
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
139
140
141
142
143
144
145
146
147
148
import numpy as np
import matplotlib.pyplot as plt
import pickle
from scenarioPlotter import samples_needed_each_test_plotter, diagnostic_probability_plotter, \
time_since_infection_plotter, case_seasonality_plotter, num_cases_legend_plotter, diagnostic_legend_plotter
# initialize parameters that describe the system
initialize_years = 2
test_detection_probs = [([0.0]*1 + [0.01*np.exp(np.log(70)/3)**y for y in range(4)] + [0.7]*10
+ [0.7*np.exp(np.log(0.00007)/30)**y for y in range(31)] + [0]*250),
([0.0]*1 + [0.01*np.exp(np.log(80)/2)**y for y in range(3)] + [0.8]*15
+ [0.8*np.exp(np.log(0.00008)/150)**y for y in range(151)] + [0]*150),
([0.0]*2 + [0.01*np.exp(np.log(90)/10)**y for y in range(11)] + [0.9]*20
+ [0.9*np.exp(np.log(0.00009)/300)**y for y in range(301)])]
diagnostic_names = ['RDT (pretend)', 'uRDT (pretend)', 'serology (pretend)']
false_pos_prob = [0.0]*len(test_detection_probs)
num_sims_each = 200
confidence_level = 0.95
# Lists of the three seasonality/sampling scenarios to explore
seasonal_scalar_list = [[1] * 365, [(1 + 0.5 * np.sin(2 * np.pi / 365 * y - 166)) for y in range(1, 366)],
[(1 + 0.5 * np.sin(2 * np.pi / 365 * y - 166)) for y in range(1, 366)]]
surveillance_days_list = [[250], [250], [90]]
# Create plot based on loaded output from previous run (if this parameter set has not already been saved, run it
# from run_diagnosticTestComparisons_plotPanel.py
pop_size = 500
sim_output_list = pickle.load(open("sim_output_list_pop%i_CI%i.p" % (pop_size, round(confidence_level * 100)), "rb"))
sim_num_samples_needed_list = pickle.load(open("sim_num_samples_needed_list_pop%i_CI%i.p" % (pop_size, round(confidence_level * 100)), "rb"))
pop_size = 1000
num_case_in_year_list = [round(pop_size * y) for y in [0.03, 0.06, 0.09, 0.22, 0.35, 0.48, 0.61, 0.74, 0.87, 1, 1.25]]
sim_output_list_1000 = pickle.load(
open("sim_output_list_pop%i_CI%i.p" % (pop_size, round(confidence_level * 100)), "rb"))
sim_num_samples_needed_list_1000 = pickle.load(
open("sim_num_samples_needed_list_pop%i_CI%i.p" % (pop_size, round(confidence_level * 100)), "rb"))
pop_size = 500
num_case_in_year_list = [round(pop_size * y) for y in [0.03, 0.06, 0.09, 0.22, 0.35, 0.48, 0.61, 0.74, 0.87, 1, 1.25]]
sim_output_list_500 = pickle.load(
open("sim_output_list_pop%i_CI%i.p" % (pop_size, round(confidence_level * 100)), "rb"))
sim_num_samples_needed_list_500 = pickle.load(
open("sim_num_samples_needed_list_pop%i_CI%i.p" % (pop_size, round(confidence_level * 100)), "rb"))
# Create multi-panel plot
# - Each columns of the new plot panel should correspond to a different assumption about seasonality of cases.
# The first row shows the expected number of cases through the year. The second row show the proportion of the
# population that has had their most recent infection within a certain number of days. The final row shows the
# number of samples needed to achieve a particular certainty that disease would be detected if circulation were
# occurring.
# - The last column (before the data) shows the legend for the different transmission intensity scenarios and the
# plot of the diagnostic test detection probabilities
# Set up the axes with gridspec
fig = plt.figure(figsize=(14, 10.5))
grid = plt.GridSpec(6, 4, hspace=1, wspace=0.6)
# iterate through three seasonality/sampling scenarios
for scenario in range(len(surveillance_days_list)):
# values for this seasonality / surveillance scenario
seasonal_scalar = seasonal_scalar_list[scenario]
surveillance_days = surveillance_days_list[scenario]
# Plotting
# first plotted row: seasonality in cases
ax_seasonality = fig.add_subplot(grid[0:2, scenario])
# plot expected number of cases on each day of year (seasonality of transmission)
case_seasonality_plotter(pop_size=pop_size, seasonal_scalar=seasonal_scalar,
num_case_in_year_list=num_case_in_year_list, surveillance_days=surveillance_days,
ax=ax_seasonality)
# second plotted row: 500 individuals
pop_size = 500
sim_output = sim_output_list_500[scenario]
sim_num_samples_needed = sim_num_samples_needed_list_500[scenario]
num_case_in_year_list = [round(pop_size * y) for y in
[0.03, 0.06, 0.09, 0.22, 0.35, 0.48, 0.61, 0.74, 0.87, 1, 1.25]]
ax_samples_needed_500 = fig.add_subplot(grid[2:4, scenario])
# Get the mean values as well as the 95CI for plotting lines and shaded CI regions
# save results in lines_for_plot
values_for_plot = [[[None]*len(num_case_in_year_list) for n in range(3)]
for m in range(len(test_detection_probs))]
for i1 in range(len(sim_num_samples_needed)):
for i2 in range(len(sim_num_samples_needed[i1])):
# mean value
values_for_plot[i1][0][i2] = np.percentile(sim_num_samples_needed[i1][i2], 50)
# upper 95% CI
values_for_plot[i1][1][i2] = min(pop_size, np.percentile(sim_num_samples_needed[i1][i2], 97.5))
# lower 95% CI
values_for_plot[i1][2][i2] = np.percentile(sim_num_samples_needed[i1][i2], 2.5)
# Call the plotting function to create the mean lines and 95CI areas for each of the diagnostic tests
samples_needed_each_test_plotter(x_values=num_case_in_year_list, values_for_plot=values_for_plot, pop_size=pop_size,
diagnostic_names=diagnostic_names, confidence_level=confidence_level,
num_case_in_year_list=num_case_in_year_list, ax=ax_samples_needed_500)
# third plotted row: 1000 individuals
pop_size = 1000
sim_output = sim_output_list_1000[scenario]
sim_num_samples_needed = sim_num_samples_needed_list_1000[scenario]
num_case_in_year_list = [round(pop_size * y) for y in
[0.03, 0.06, 0.09, 0.22, 0.35, 0.48, 0.61, 0.74, 0.87, 1, 1.25]]
ax_samples_needed_1000 = fig.add_subplot(grid[4:6, scenario])
# Get the mean values as well as the 95CI for plotting lines and shaded CI regions
# save results in lines_for_plot
values_for_plot = [[[None]*len(num_case_in_year_list) for n in range(3)]
for m in range(len(test_detection_probs))]
for i1 in range(len(sim_num_samples_needed)):
for i2 in range(len(sim_num_samples_needed[i1])):
# mean value
values_for_plot[i1][0][i2] = np.percentile(sim_num_samples_needed[i1][i2], 50)
# upper 95% CI
values_for_plot[i1][1][i2] = min(pop_size, np.percentile(sim_num_samples_needed[i1][i2], 97.5))
# lower 95% CI
values_for_plot[i1][2][i2] = np.percentile(sim_num_samples_needed[i1][i2], 2.5)
# Call the plotting function to create the mean lines and 95CI areas for each of the diagnostic tests
samples_needed_each_test_plotter(x_values=num_case_in_year_list, values_for_plot=values_for_plot, pop_size=pop_size,
diagnostic_names=diagnostic_names, confidence_level=confidence_level,
num_case_in_year_list=num_case_in_year_list, ax=ax_samples_needed_1000)
# Add legends for different transmission intensities and diagnostic tests.
ax_legend_intensity = fig.add_subplot(grid[0:3, -1])
num_cases_legend_plotter(num_case_in_year_list=[0.03, 0.06, 0.09, 0.22, 0.35, 0.48, 0.61, 0.74, 0.87, 1, 1.25],
ax=ax_legend_intensity)
ax_legend_test = fig.add_subplot(grid[3, -1])
diagnostic_legend_plotter(diagnostic_names, ax_legend_test)
# plot diagnostic tests
ax_diagnostic_test = fig.add_subplot(grid[4:6, -1])
diagnostic_probability_plotter(test_detection_probs=test_detection_probs, diagnostic_names=diagnostic_names)
fig.suptitle('Population size = 500 vs 1000; Confidence level = %.2f' % round(confidence_level, 2))
plt.show()
fig.savefig('NumSamplesNeeded_pop500vs1000_CI%i.pdf' % round(confidence_level * 100))