Prevestim: true prevalence calculator

Bayesian estimation of true prevalence. See documentation for a complete explanation. The test parameters specificity and sensitivity are given by θ⁻ and θ⁺. The total number of survey participants is indicated by N, and the positives amongst them by N⁺. Enter those parameters below and then click on Run the calculator. The results are computed via a Bayesian probabilistic simulation (code shown below) written in WebPPL and displayed in section Results.

Results

Citation

Code

Parameter	Value
Crude estimate
Maximum a-posteriori estimate (MAP) of λ
Expected a-posteriori estimate (EAP) of λ
Credible interval level
Credible interval	–

/* INPUT */
var theta_plus  = globalStore.thetaPlus / 100   // Sensitivity
var theta_minus = globalStore.thetaMinus / 100  // Specificity
var n    = globalStore.n                        // Survey participants
var npos = globalStore.npos                     // Number of positives


/* PRIORS */
/* As default we will use a Beta prior with mean the crude estimate 
   (npos/n). We use the μ-ν-parameterization of the Beta distribution.
   The ν-parameter is the so-called "sample size," but not to be confused
   with n. Advisable is to leave it set to 100. */
var mu = npos/n // mean
var nu = 100    // "sample size"
var prior = Beta({
  a: mu*nu,
  b: (1 - mu)*nu
})


/* SIMULATION */
var prevalence_model = function()
{
  /* Draw prior seroprevalence
     We assume that all individuals have the same true
     underlying probability of seropositive, i.e. we
     do a marginal model with respect to e.g. age, sex, etc.

     Using notation from Stringhini*/
    var pi = sample(prior) // prior seroprevalence mean

    /* Condition on observerd data, adjust for test
       specificity and sensitivity*/
    observe(
    Binomial({
        p: pi*theta_plus + (1- pi)*(1 - theta_minus),
        n: n
    }), npos)

    /* Return the updated seroprevalence rate */
    return {pi}
}


/* INFERENCE */
/* Sequential Monte-Carlo estimation using 50,000 particles.
   The simulation will take several seconds. */
var PosteriorDistribution = Infer({method: 'SMC',
                       particles: 50000,
                       model: prevalence_model})

var posteriorPi = marginalize(PosteriorDistribution, "pi")
var posteriorPiMAP = posteriorPi.MAP().val

// From https://mhtess.github.io/bdappl/chapters/03-simpleModels.html
var getSamples = function(Dist, key){
  return _.map(_.map(PosteriorDistribution.samples, "value"), key)
}

var credibleInterval = function(mySamples, credMass){
  var sortedPts = sort(mySamples)
  var ciIdxInc = Math.ceil(credMass*sortedPts.length)
  var nCIs = sortedPts.length - ciIdxInc

  var ciWidth = map(function(i){
    sortedPts[i + ciIdxInc] - sortedPts[i]
  },_.range(nCIs))

  var i = _.indexOf(ciWidth, _.min(ciWidth))

  return [sortedPts[i], sortedPts[i+ciIdxInc]]
}

var posteriorSamples = getSamples(PosteriorDistribution, "pi")

// Credible Interval
// see https://easystats.github.io/bayestestR/articles/credible_interval.html
var masspercent = .95
var ci = credibleInterval(posteriorSamples, masspercent)

globalStore.mu = mu
globalStore.posteriorPiMAP = posteriorPiMAP
globalStore.posteriorPiEAP = expectation(posteriorPi)
globalStore.ciLevel = masspercent
globalStore.ciLow = ci[0]
globalStore.ciHigh = ci[1]

© 2020 Viktor Senderov, Jan Kudlicka, and Savina Stoitsova. This software is provided under the Apache License 2.0 and the authors bear no liability for its usage. This research has been financed under the EU Horizon 2020 research and innovation program, Marie Skłodowska-Curie grant agreement PhyPPL No 898120.