# Gibbs sampling for Bayesian estimation of disease prevalence

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