Guest User

Untitled

a guest
Oct 26th, 2020
76
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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×