Hi @cowirihy,
This is interesting! I am by no means an expert but thought I could offer a different perspective.
First, I am not sure if the Dirichlet would work because if you constraint the sum of your event probabilities to 1 then you are telling the model to treat the events as being dependent.
I think an alternative would be similar to what Krushke from DBDA2 chapter 12 does and model the events as either the event happened or did not happen. That way you can model the joint event types using a Beta distribution. To check for dependence, you would compute the marginal event probabilities from the posterior and subtract the joint posterior probabilities from the product of the marginal posterior probabilities you computed. You can then define a region of practical equivalence (ROPE) to decide what probability mass would constitute dependence/independence.
Here is some code showing how this could be done with PyMC:
import arviz as az
import numpy as np
import pymc as pm
# From table in post
a0b0 = 33
a0b1 = 343
a0b2 = 30
a1b0 = 64
a1b1 = 360
a1b2 = 70
# Long format A event type, B event type, event count
d = np.concatenate(
(
np.stack([np.array([0, 0, 1]) for _ in range(a0b0)]),
np.stack([np.array([0, 0, 0]) for _ in range(900 - a0b0)]), # these are the failures
np.stack([np.array([0, 1, 1]) for _ in range(a0b1)]),
np.stack([np.array([0, 1, 0]) for _ in range(900 - a0b1)]),
np.stack([np.array([0, 2, 1]) for _ in range(a0b2)]),
np.stack([np.array([0, 2, 0]) for _ in range(900 - a0b2)]),
np.stack([np.array([1, 0, 1]) for _ in range(a1b0)]),
np.stack([np.array([1, 0, 0]) for _ in range(900 - a1b0)]),
np.stack([np.array([1, 1, 1]) for _ in range(a1b1)]),
np.stack([np.array([1, 1, 0]) for _ in range(900 - a1b1)]),
np.stack([np.array([1, 2, 1]) for _ in range(a1b2)]),
np.stack([np.array([1, 2, 0]) for _ in range(900 - a1b2)]),
)
)
a_event_uniques, a_event_idx = np.unique(d[:, 0], return_inverse=True)
b_event_uniques, b_event_idx = np.unique(d[:, 1], return_inverse=True)
y_event_counts = d[:, 2]
coords = {
"a_event": a_event_uniques,
"b_event": b_event_uniques,
"n_obs": np.arange(len(y_event_counts))
}
with pm.Model(coords=coords) as model:
covariate_a = pm.Data("covariate_a", a_event_idx, dims="n_obs")
covariate_b = pm.Data("covariate_b", b_event_idx, dims="n_obs")
joint_prob = pm.Beta("joint_prob", alpha=2, beta=2, dims=("a_event", "b_event"))
pm.Bernoulli("likelihood", p=joint_prob[covariate_a, covariate_b], observed=y_event_counts, dims="n_obs")
with model:
idata = pm.sample()
# Get the marginals
a0_marginal = idata.posterior['joint_prob'].sel(a_event=0).mean(('b_event'))
b0_marginal = idata.posterior['joint_prob'].sel(b_event=0).mean(('a_event'))
# Get difference between product of marginals and joint (should be close to zero if independent)
difference_between_joint_and_prod_marginal = idata.posterior['joint_prob'].sel(a_event=0, b_event=0) - (a0_marginal * b0_marginal)
# Let's say that they are independent if most of the posterior mass is inside 0% to 5%.
az.plot_posterior(difference_between_joint_and_prod_marginal, rope=[0, 0.05])
You’ll get something like this:
Here you can see for the particular events A0 and B0 with a ROPE spanning 0 to 0.05 (which may be too wide) we would say that the events are independent because almost all of the posterior mass lies in the region.
