Guest User

Untitled

a guest
Oct 26th, 2020
50
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[1]
  17.     β = subpopulations[2]
  18.     γ = subpopulations[3]
  19.     δ = subpopulations[4]
  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[1]"] .+ chain["subpopulations[2]"] .+ chain["subpopulations[3]"]
  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.  
RAW Paste Data