Hello,
I am trying to calculate a disproportionality measure using an Empirical Bayes method. A really good work has been done in SAS to achieve the goal and available here : Disproportionality-Measures-with-SAS/ebgm.sas at master · statmike/Disproportionality-Measures-with-SAS · GitHub
I am trying to replicate this method in Python using Pymc library. However I’m can’t get the MCMC to work to estimate the parameters. The distribution is assumed to be a mixture of two gamma distribution (priors) and the data (N11, and E11 observed) allow to generate the posterior. I get an output in python but it is not aligned with the SAS results.
In particular here is the SAS code I am trying to replicate in Python :
ods output PostSumInt = work.MCMC_PARMS;
proc mcmc data=&outds. seed=32259 nmc=10000 thin=5 nthread=6 propcov=quanew monitor=(Mix_p alpha1 alpha2 beta1 beta2);
parms Mix_p 0.3333 alpha1 .2 alpha2 2 beta1 .1 beta2 4;
/**Generic gamma priors*/
prior alpha1 beta1 ~ gamma(1, iscale=1);
prior alpha2 beta2 ~ gamma(1, iscale=1);
prior Mix_p ~ uniform(0,1);
const1=lgamma(alpha1+N11) - lgamma(alpha1) - lgamma(N11+1);
const2=lgamma(alpha2+N11) - lgamma(alpha2) - lgamma(N11+1);
LogL1=const1 - N11*log(1+beta1/E11) - alpha1*log(1+E11/beta1);
LogL2=const2 - N11*log(1+beta2/E11) - alpha2*log(1+E11/beta2);
llike = log(Mix_p*exp(LogL1) + (1-Mix_p)*exp(LogL2));
model N11 ~ general(llike);
run;
Here is my code in Python :
with pm.Model() as model2:
def likelihood(theta,E11,N11):
Mix_p,alpha1,alpha2,beta1,beta2 = theta
const1=at.gammaln(alpha1+N11) - at.gammaln(alpha1) - at.gammaln(N11+1)
const2=at.gammaln(alpha2+N11) - at.gammaln(alpha2) - at.gammaln(N11+1)
LogL1=const1 - N11*np.log(1+beta1/E11) - alpha1*np.log(1+E11/beta1)
LogL2=const2 - N11*np.log(1+beta2/E11) - alpha2*np.log(1+E11/beta2)
return(np.log(Mix_p*np.exp(LogL1) + (1-Mix_p)*np.exp(LogL2)))
Mix_p = pm.Uniform('Mix_p', 0, 1)
alpha1 = pm.Exponential('alpha1', lam=0.2)
beta1 = pm.Exponential('beta1', lam=0.1)
alpha2 = pm.Exponential('alpha2', lam=2)
beta2 = pm.Exponential('beta2', lam=4)
#g1 = pm.Gamma('gamma1',alpha = alpha1, beta = beta1)
#g2 = pm.Gamma('gamma2',alpha = alpha2, beta = beta2)
theta = (Mix_p,alpha1,alpha2,beta1,beta2)
like = pm.Potential('likelihood',likelihood(theta,E11,N11))
I am not very familiar with Pymc so I have tried to user other people experience on other feeds but I don’t seem to get the right result.
Thank you in advance
Victor