Hi, I’m cross-posting from stackoverflow since I just realised I might be better off asking here. https://stackoverflow.com/questions/56928544/how-can-i-calculate-risk-ratio-using-pymc3
I’d like to model the risk ratio for data in a survey, but I’m having difficulty framing the problem to pymc3. Risk ratio is the ratio of probability of some outcome given an exposure, to the probability without exposure.
So as I understand it, risk ratio is: (incidence in exposed group) / (incidence in notExposed group)
I’ve tried posing the problem as shown in the sample code below:
import pymc3 exposed_num_observations = 25 exposed_observed = 11 notExposed_num_observations = 853 notExposed_observed = 202 with pm.Model() as rr_mdoel: # Prior p_exposed = pm.Uniform("p_exposed", 0 ,1) p_notExposed = pm.Uniform("p_notExposed", 0 ,1) sd = pm.Uniform("sd", 0,5) # Likelihood exposed = pm.Binomial("exposed", p=p_exposed, n=exposed_num_observations, observed=exposed_observed) notExposed = pm.Binomial("notExposed", p=p_notExposed, n=notExposed_num_observations, observed=notExposed_observed) risk_ratio = pm.Normal("risk_ratio", mu=(p_exposed/p_notExposed), sd=sd) # Inference Run/ Markov chain Monte Carlo trace = pm.sample(10000, chains=2, tune=1000)
I expected the
risk_ratio to always be above 0 (because you cant have a negative risk ratio), and I expected
sd to be fit to some value. The results (as shown by
az.plot_posterior(trace, round_to=3) for example) show some negative values for
sd is still uniform.
Any tips here are greatly appreciated!