Gibbs sampling for Bayesian estimation of disease prevalence

Questions
1

Hi,

I am a complete beginner with pymc3 so apologies for the basic questions.

I would like to implement the method described in:

Joseph, Lawrence, Theresa W. Gyorkos, and Louis Coupal. “Bayesian Estimation of Disease Prevalence and the Parameters of Diagnostic Tests in the Absence of a Gold Standard.” American Journal of Epidemiology 141, no. 3 (February 1, 1995): 263–72. https://doi.org/10.1093/oxfordjournals.aje.a117428.
PDF

I am not really clear how implement the gibbs sampling proposed there in PyMC3 (ie how to declare conditional distributions A1-A5) or whether it is better to use the full likelihood (unnumbered equation top left P265). I am interested initially in the case of a single diagnostic test and unknown ground truth.

I first tried directly, which works when alpha and beta (specificity and sensitivity of the diagnostic test) are specified, but doesn’t seem to converge well when beta has a broad prior.

from pymc3.distributions.dist_math import binomln, logpow, bound

def likelihood(alpha,beta,prev):
    def logp_(k,N):
        return bound(
            binomln(N,k) + logpow(alpha*(1-prev)+beta*prev,k) + logpow(1-alpha*(1-prev)-beta*prev,N-k),
            0 <= k, k<= N,
            0 <= prev, prev<=1,
            0 <= alpha, alpha<=1,
            0 <= beta, beta<=1)
    return logp_
        
with pm.Model() as model:
    #pi = pm.Uniform('prev',lower=0,upper=1)
    prev = pm.Beta('prev',1,1)
    #C = 1 - 0.05 # specifcity
    alpha = pm.Beta('spec',alpha=6,beta=96)
    beta = pm.Beta('sens',1,1)

    k = pm.DensityDist('k', likelihood(alpha,beta,prev), observed = {'k': 50, 'N': 100})
    trace = pm.sample()

I suppose I need to do the multi-step Gibbs sampling with the latent variables as proposed in the paper, but I’m not sure how to implement that in pymc3.

I tried this but the results look the same as the above:

a = 50;
b = 50;
N = a+b;
with pm.Model() as prev_model:
    #pi = pm.Uniform('prev',lower=0,upper=1)
    pi = pm.Beta('prev',1,1)
    #C = 1 - 0.05 # specifcity
    C = pm.Beta('spec',alpha=6,beta=96)
    S = pm.Beta('sens',1,1)
    
    # eq A1
    Y1 = pm.Binomial('Y1',a, (pi*S)/(pi*S + (1-pi)*(1-C)))
    Y2 = pm.Binomial('Y2',b, (pi*(1-S))/(pi*(1-S) + (1-pi)*C))
    
    pm.sample()
1 Like