Hi, I am completely new to Bayesian statistics and pymc3 and have run into trouble.
I am trying to measure the change in the number of deaths in hospitals before and after the intervention. I decided to mostly follow BEST procedure (Bayesian Estimation Supersedes the T-Test — PyMC3 3.10.0 documentation) to evaluate the effect. However, I am not sure if my model is correct.
My questions are:
- is my model adequate to the situation?
- can we use the difference of means in this case? (the resulting posterior distributions don’t look normal to me)
- what are other ways to evaluate the pre and post difference?
Here is my code:
import pymc3 as pm pre_mu = pre.mean() pre_sd = pre.std() * 2 post_mu = post.mean() post_sd = post.std() * 2 σ_low = 0 σ_high = 10 with pm.Model() as model: pre_mean = pm.Normal("pre_mean", mu=pre_mu, sd=pre_sd) post_mean = pm.Normal("post_mean", mu=post_mu, sd=post_sd) pre_std = pm.Uniform("pre_std", lower=σ_low, upper=σ_high) post_std = pm.Uniform("post_std", lower=σ_low, upper=σ_high) diff_of_means = pm.Deterministic("difference of means", pre_mean - post_mean) diff_of_stds = pm.Deterministic("difference of stds", pre_std - post_std) effect_size = pm.Deterministic( "effect size", diff_of_means / np.sqrt((pre_std ** 2 + post_std ** 2) / 2) ) group1 = pm.NegativeBinomial("pre", mu=pre_mu, alpha=pre_sd, observed=pre_nzk) group2 = pm.NegativeBinomial("post", mu=post_mu, alpha=post_std, observed=post_nzk) trace = pm.sample(5000)