Thank you for your response! My most recent problem is a two connected hierarchical models, small dataset here:
I want to test if there is a difference in scientific ability between the Jefferson school and the Washington school as measured by both of the science tests (NOTE: in the full dataset there are more than 2 schools). In each of the 2 schools each of 4 students takes each of 2 tests 3 different times. I am able to estimate one beta for each school-test pair and pool into two betas one level above the test level:
But this model doesn’t cluster scores at the StudentID level and that clustering is what I am struggling with. I want to estimate a model like this:
I think I’m encoding the hierarchy correctly but something is going wrong with how I’m giving PyMC the indicators for SchoolID and StudentID and I don’t understand why. When I try to sample this model it gives a ValueError: Incompatible Elemwise input shapes.
How does this clustering work with PyMC? In the final model call I’m telling PyMC to cluster by school and student by subsetting a pm.Normal object like:
mod_test1 = pm.Normal("mod_cluster_test1",
mu = beta_test1_rep[school, student],
sigma = 2,
observed = df.science_test1,
dims = "match")
But I suspect this is incorrect, or at least not flexible enough. It’s only this one model that’s going wrong because I’ve done this double subsetting before and it’s worked just fine. Perhaps I’ve encoded something else wrong in this model?