 # Untitled

a guest
Oct 26th, 2020
250
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. using Turing, StatsPlots, Distributions, Printf
2.
3. @model function antibody_and_pcr_tests(
4.             results,n,deaths,population
5.             )
6.     # Model for Streeck et al
7.     #Prior for the ifr
8.     ifr ~ Uniform(0.0,0.05)
9.     # The population is grouped into 4 partitions
10.     # α - people who have antibodies but are not currently infected
11.     # β - people who have antibodies and are currently infected
12.     # γ - people who are currently infected but do not have antibodies
13.     # δ - people who are neither infected, nor have antibodies
14.     #Uniform priors
15.     subpopulations ~ Dirichlet([1,1,1,1])
16.     α = subpopulations
17.     β = subpopulations
18.     γ = subpopulations
19.     δ = subpopulations
20.     #
21.     #How many people test positive for antibody test from subpopulations α,β,γ,δ = (1,2,3,4) respectively
22.     #There is some uncertainty associated to these variables, which we attempt to model
23.     #These first variables are sensitivities
24.     # Streeck et al claim that sensitivity is 90.9%, but no number of samples given
25.     # So we assume it's something like 200 and of these 180 are positive and 20 negative
26.     antibody_test_1 ~ Beta(1+180,1+20)
27.     #For people that are currently infected but have antibodies we assume they have less antibodies, so smaller positivity rate.
28.     #Probably also more uncertainty with these numbers
29.     antibody_test_2 ~ Beta(1+80,1+20)
30.     #These people are infected but without substantial antibodies. Let's assume the false positive rate is slightly higher here
31.     #Than completely healthy people
32.     antibody_test_3 ~ Beta(1+20,1+80)
33.     #For these, Streeck et. all claim that 1656 samples of the manufacturer had 15 false positives
34.     #Then they themselves get 2 false positives for 109 tests
35.     #So in total, we know of 17 false positives and  1748 negatives
36.     antibody_test_4 ~ Beta(1+17,1+1748)
37.     #
38.     #Same numbers for the PCR test
39.     #
40.     #We know that the PCR test tests quite long even after infection, so let's allow for around 5% "false positives" amongst those with antibodies
41.     #of people who were infected but are not longer infected
42.     pcr_test_1 ~ Beta(1+5,1+95)
43.     #Let's say PCR has around 70% sensitivity, based on 500 tests.
44.     pcr_test_2 ~ Beta(1+350,1+150)
45.     pcr_test_3 ~ Beta(1+350,1+150)
46.     #slightly less than 1% false positives amongst complete negatives due to e.g. contamination etc seems reasonable
47.     #But we're not so sure, so we act like we've seen 1 false positive test after doing 150 tests
48.     pcr_test_4 ~ Beta(1+1,1+149)
49.     #
50.     #
51.     #Probability of getting a positive antibody test
52.     prob_antibody = α*antibody_test_1 + β*antibody_test_2 + γ*antibody_test_3 + δ*antibody_test_4
53.     #Probability of getting a positive PCR test
54.     prob_pcr = α*pcr_test_1 + β*pcr_test_2 + γ*pcr_test_3 + δ*pcr_test_4
55.     #Probability that both are positive (we assume a certain statistical independence here that is probably not warranted)
56.     prob_both = α*antibody_test_1*pcr_test_1 + β*antibody_test_2*pcr_test_2 + γ*antibody_test_3*pcr_test_3 + δ*antibody_test_4*pcr_test_4
57.     #
58.     categoryprobs = [
59.                     prob_antibody - prob_both,
60.                     prob_both,
61.                     prob_pcr - prob_both,
62.                     1 - prob_pcr - prob_antibody + prob_both
63.                     ]
64.     results ~ Multinomial(n,categoryprobs)
65.     #
66.     total_infected = α + β + γ
67.     #
68.     deaths ~ Binomial(population,ifr*total_infected)
69.     #
70. end
71.
72.
73. nsamples = 919 #Total number of samples
74. #We ignore the two old PCR tests that weren't serology positive, as we do account for serology false-negatives already
75. results = [102,10,23, nsamples - 102 - 23 - 10]
76. #Total population and deaths in the region
77. population = 12597
78. deaths = 7
79. #
80. iterations = 500000
81. chain = sample(
82.     antibody_and_pcr_tests(
83.         results,nsamples,deaths,population
84.         ),
85.     HMC(0.03,10), iterations, progress=true);
86.
87.
88. total_positive_res = chain["subpopulations"] .+ chain["subpopulations"] .+ chain["subpopulations"]
89.
90. fig1 = density(total_positive_res,xlim=(0.05,0.2),label="",title="Infection rate posterior",color=:black)
91. q1t = quantile(total_positive_res[:],0.025)
92. q2t = quantile(total_positive_res[:],0.975)
93. vline!([q1t], label=(@sprintf "2.5%% quantile = %.2f%%" q1t*100),color=:red,legend=:topleft)
94. vline!([q2t], label=(@sprintf "97.5%% quantile = %.2f%%" q2t*100),color=:red,legend=:topleft)
95.
96. ifr = (collect(chain[:ifr].data)[:])
97. fig2 = density(ifr,xlim=(0.0,0.01),label="",title="IFR posterior",color=:black)
98. q1 = quantile(ifr, 0.025)
99. q2 = quantile(ifr,0.975)
100. vline!([q1], label=(@sprintf "2.5%% quantile = %.2f%%" q1*100),color=:red,legend=:topright)
101. vline!([q2], label=(@sprintf "97.5%% quantile = %.2f%%" q2*100),color=:red,legend=:topright)
102.
103. fig3 = plot(fig1,fig2)
104.