Unexpected behavior while implementing OpenBUGS code in PyMC3

Questions
1

Hello!
So i am trying to implement a certain model written in BUGS into PyMC3.

The OpenBUGS code goes something like this:

model{

lambdac ~ dgamma(ac,bc)
lambdav ~ dgamma(av,bv)

sc ~ dnorm(expsc, tausc)
sv ~ dnorm(expsv, tausv)

tausc <- 1/(nv*varsc)
tausv <- 1/(nv*varsv)

expsc <- nc*(1- (1-exp(-lambdac*D))/(lambdac*D)) / lambdac
varsc <- (2*exp(-lambdac*D) + 4*exp(-lambdac*D)/(lambdac*D) -
pow((1-(1-exp(-lambdac*D))/(lambdac*D)),2))/pow(lambdac,2)

expsv <- nv*(1- (1-exp(-lambdav*D))/(lambdav*D))/lambdav
varsv <- (2*exp(-lambdav*D) + 4*exp(-lambdav*D)/(lambdav*D) -
pow((1-(1-exp(-lambdav*D))/(lambdav*D)),2))/pow(lambdav,2)

cases ~ dpois(meancases)
meancases <- sc*lambdac + sv*lambdav

xv ~ dbin(theta, cases)
theta <- sv*lambdav/(sc*lambdac + sv*lambdav)

VE <- 1 - lambdav/lambdac

}

and they have taken the following values: av=1, bc=0.01917808, ac= 2.428571, bv=0.01917808, nc=17511, nv=17411, sc=2222,sv=2214, cases=170, xv=8, D=0.29.
As an alternate, we also have, ac=1, bc=2222, av=0.7, bv=2214 (rest everything same)

My implementation is the following:

with covid_model:
    lambda_v = pm.Gamma("lambda_v", alpha=a_v, beta=b_v)
    lambda_c = pm.Gamma("lambda_c", alpha=a_c, beta=b_c)

    exp_sv = nv*(1 - (1 - tt.exp(-lambda_v*D))/(lambda_v*D))/lambda_v
    var_sv = (2*tt.exp(-lambda_v*D) + (4*tt.exp(-lambda_v*D)/lambda_v*D) - tt.sqr(1 + (tt.exp(-lambda_v*D) - 1)/lambda_v*D))/tt.sqr(lambda_v)
    tau_sv = 1/(nv*var_sv)

    exp_sc = nc*(1 - (1 - tt.exp(-lambda_c*D))/(lambda_c*D))/lambda_c
    var_sc = (2*tt.exp(-lambda_c*D) + (4*tt.exp(-lambda_c*D)/lambda_c*D) - tt.sqr(1 + (tt.exp(-lambda_c*D) - 1)/lambda_c*D))/tt.sqr(lambda_v)
    tau_sc = 1/(nc*var_sc)

    sv = pm.Normal("sv", mu=exp_sv, tau=tau_sv)
    sc = pm.Normal("sc", mu=exp_sc, tau=tau_sc)

    mean_case = sc*lambda_c + sv*lambda_v
    cases = pm.Poisson("cases", mu=mean_case, observed=n_cases)
    theta_v = sv*lambda_v/mean_case
    theta_c = sc*lambda_c/mean_case
    xv = pm.Binomial("xv", n=cases, p=theta_v, observed=x_v)
    xc = pm.Binomial("xc", n=cases, p=theta_c, observed=x_c)
    ve = pm.Deterministic("ve", (1 - (lambda_v/lambda_c)))

With the first set of data, I get an error regarding the initial values of the parameters.
With the alternate set of data, i dont get any such error, but the model fails to converge spectacularly.

image
image883×213 16.9 KB

What i need is the ve parameter. it should be a real number between 1 and 0 (its actually the vaccine efficacy). I have had worse r-hat values with NUTS, which had over 18000 divergences in 50000 draws. This output is after using Metropolis as the step method, no divergences, but the results dont make any sense.

As i am new to PyMC3 and Probabilistic programming in general, what exactly am i doing wrong?

(The BUGS code can be found in this paper, with other details: Improving Bayesian estimation of Vaccine Efficacy

2

in the line that calculates var_sc, you have lambda_v right at the end when it should be lambda_c ? Not sure if that’s sufficient to explain all the issues.

3

oh yes figured out that typo, forgot to update it here

4

Figured it out myself. Not to mention there is a typo at the end of s_c line, where it should be lambda_c, sv and sc are observed random variables so I needed to pass the observed values 2214 to sv and 2222 to sc. That solves all of it