How to properly index multilevel model with repeated measures data

Questions version agnostic
5

@jessegrabowski Yes, I am also quite concerned about taking the mean of a. It was the only way I could get it to not throw out an error of p being multidimensional. Do you have a reference I could check out on how to index into long-form data? The data I have is already in long-form each observation is a row. When I was implementing the pymc code I did initially try to only index my level 0 parameters like b_gender only by dims=“gender” but when I try to construct the a variable that has dims = “patients” it would complain that it could not broadcast variables.

edit:
I think I was able to make some progress thanks to your suggestions. I am now able to sample from the following model:

with pm.Model(coords = coords) as multilevel_model:
    # Prior
    sigma_a = pm.HalfCauchy("sigma_a", 2)

    # Patient Baseline
    gamma = pm.Normal("gamma_0", 0, 1)
    gamma_gender = pm.Normal("gamma_gender", 0, 1, dims = "gender")
    gamma_unit = pm.Normal("gamma_unit", 0, 1, dims = "unit")
    gamma_care = pm.Normal("gamma_care", 0, 1, dims = "care")
    gamma_age = pm.Normal("b_age", 0, 1)

    mu_a = pm.Deterministic(
        "mu_a",
        var = gamma + 
        gamma_gender[gender_idx] + 
        gamma_unit[unit_idx] + 
        gamma_care[care_idx] + 
        gamma_age*age
    )

    # Unexplained Patient Variation
    epsilon_a = pm.Normal("epsilon_a", 0, 1)
    a = pm.Deterministic("a", var=mu_a + sigma_a * epsilon_a)
    
    # Common slopes priors
    b_time = pm.Normal("b_time", 0, 1)
    b_hematocrit = pm.Normal("b_hematocrit", 0, 1)
    b_min_bp = pm.Normal("b_min_bp", 0, 1)
    b_max_bp = pm.Normal("b_max_bp", 0, 1)
    b_min_temp = pm.Normal("b_min_temp", 0, 1)
    b_max_temp = pm.Normal("b_max_temp", 0, 1)
    b_min_rsp = pm.Normal("b_min_rsp", 0, 1)
    b_max_rsp = pm.Normal("b_max_rsp", 0, 1)
    b_min_pls = pm.Normal("b_min_pls", 0, 1)
    b_max_pls = pm.Normal("b_max_pls", 0, 1)
    b_min_sgr = pm.Normal("b_min_sgr", 0, 1)
    b_max_sgr = pm.Normal("b_max_sgr", 0, 1)
    b_min_pn = pm.Normal("b_min_pn", 0, 1)
    b_max_pn = pm.Normal("b_max_pn", 0, 1)

# Expectation
    p = pm.Deterministic(
        "p", 
        pm.math.invlogit(
            a[patient_idx] + 
            b_time*time +
            b_hematocrit*hematocrit +
            b_min_bp*min_BP +
            b_max_bp*max_BP +
            b_min_temp*min_temp +
            b_max_temp*max_temp +
            b_min_rsp*min_rsp +
            b_max_rsp*max_rsp +
            b_min_pls*min_pls + 
            b_max_pls*max_pls +
            b_min_sgr*min_sgr +
            b_max_sgr*max_sgr +
            b_min_pn*min_pn +
            b_max_pn*max_pn
        ),
        dims="obs"
    )
    likelihood = pm.Bernoulli("likelihood", p = p, observed = transferred, dims="obs")

multilevel_model-1
multilevel_model-11920×421 47.7 KB

The only piece that I am still struggling with is the variable “epsilon_a” which represents the unobserved variation between patients. I am unable to set the dims of “epsilon_a” to be the number of patients in the data set. I get broadcasting errors if I do.