Multi-Multilevel Modeling

Questions
7

Thanks @OriolAbril for all the details! At the moment I have non-ragged data (in terms of equal replicates per group), so I’ll focus on the first half first, but want to make sure I can get your second demo to work too. The definitions for mu and sigma are right now, but how do I handle defining the likelihood now?

This works:

group_idx, groups = pd.factorize(df_data["group"], sort=True)
replicate_idx, replicates = pd.factorize(df_data["replicate"], sort=True)

coords = {
    "group": groups,
    "obs_id": np.arange(df_data.shape[0]),
    "replicate":replicates,
}

model = pm.Model(coords=coords)

with model:
    BoundHalfNormal = pm.Bound(pm.HalfNormal, lower=0)
    
    mu_a = pm.Normal("μ_a", mu=0, sigma=10, dims="group")
    sigma_a = BoundHalfNormal("σ_a", sigma=10, dims="group")
    
    sigma_b = BoundHalfNormal("σ_b", sigma=10, dims="group")

    mu = pm.Normal("μ", mu=mu_a, sigma=sigma_a, dims=("replicate","group"))
    sigma = BoundHalfNormal("σ", sigma=sigma_b, dims=("replicate", "group"))

    obs = pm.Normal(
        "obs",
        mu=mu[replicate_idx, group_idx],
        sigma=sigma[replicate_idx, group_idx],
        observed=df_data["height"],
        dims="obs_id",
    )
    idata = pm.sample(draws=2_000, tune=1_000, chains=4, return_inferencedata=True)

image
image1024×764 25.1 KB

It’s now running (very slowly!) Does this otherwise look okay? And to be clear, the raggedness pertains to number of groups-replicates, vs number of observations per group-replicate?

2 Likes