I am trying to fit a multilevel model on some longitudinal data. In this data I have ~300k observations on ~ 1400 patients. The goal is to have a multilevel model where level 0 models a patients baseline risk and level 1 models additional risk based on the repeated measures variables.

Below is my model code and an graph image of the 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", "patients"])
gamma_unit = pm.Normal("gamma_unit", 0, 1, dims = ["unit", "patients"])
gamma_care = pm.Normal("gamma_care", 0, 1, dims = ["care", "patients"])
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,
dims = ["obs", "patients"]
)
# Unexplained Patient Variation
epsilon_a = pm.Normal("epsilon_a", 0, 1, dims = ["obs", "patients"])
a = pm.Deterministic("a", var=mu_a + sigma_a * epsilon_a, dims = ["obs", "patients"])
# 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(
pm.math.mean(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")
```

I keep getting shape error messages when I try to sample this model.

Any tips/advice would really be appreciated.