# Specifying Mixed Effects Model Formulas

I have a dataset of experiment groups, different tests, schools, and students (note students are in the dataset multiple times due to multiple tests so they're not iid either) and want to specify a mixed effects model formula.

I would think it makes sense to have the fixed effects be the experiment group, and the random effects be tests/schools/students

Would this mixed effects formula be correct?

score ~ group + (1|test) + (1|test:school) + (1|school:student)


And how would I specify that same formula using Statsmodels? Right now I'm trying the following, but I don't totally understand how a mixed effects formula translates into statsmodels, and I think I'm missing the mark here:

vc = {'school': '0 + C(school)', 'student': '0 + C(student)' }
re_formula = '1'

model = smf.mixedlm('score ~ group', test_scores_df, re_formula = re_formula, vc_formula=vc, groups = test_scores_df['test'])


There are ~20000 total observations, with ~ 2000 unique tests, ~ 50 unique schools, and ~ 5000 unique students.

As you can see there are quite a lot of gaps in the data so some students are only taking some tests, some have data for multiple different tests, etc.

Students are nested in schools, but tests I believe are crossed with respect to both students and schools as Kerby Shedden specified below, however I believe they are partially crossed, since not all students/students receive the same tests.

So for Test A maybe students 10001, 10002, 10003, 10004 and schools 151, 152, 153 got this test A, and Test B maybe students 10002, 10003, 10004, 10005 and schools 152, 153, 154 got the test. Students are definitely nested within schools though. Would this still be specified the same as a full cross?

• Can you edit the question to add more detail. How many tests, schools and student do you have ? Plesse describe the nesting/hierarchy. Jan 14, 2020 at 19:34
• Thanks Robert I just updated it Jan 14, 2020 at 21:11
• It's still not clear. Are you saying it's a 3-level model, or a crossed 2-level model. Please be explicit and say what is nested in what. Jan 14, 2020 at 21:40
• I've just updated it with hopefully more details, sorry I'm a little lost on this. Hopefully the newest update is a little more clear Jan 21, 2020 at 16:32

A mixed effects model with random intercepts for students nested in schools, but also with random intercepts for tests, since this is a crossed factor, should be appropriate here. I don't know statsmodels, but using the standard formula used with R packages such as lme4 the model would look like this:

score ~ group + (1 | test) +  (1 | school) + (1| school:student)

• Since there's within student correlation, does the nested (1| school:student) account for this or would I need an additional random effect (1| student) be needed? I.e. score ~ group + (1 | test) + (1 | school) + (1| school:student) + (1|student) Jan 21, 2020 at 19:02
• @JDoe you would only need (1 | student) if student and school are crossed, otherwise the interaction should handle it. It's a bit less clear how test is related to school and student. If they are not nested then the formula I have should OK though ideally I would like to see the data. Have you seen this ? Jan 21, 2020 at 19:28
• Thank you this answered all of my questions. I found the link to be very helpful as well, appreciate all of your help! I think it's about time I go ahead and pick up a textbook on this topic Jan 21, 2020 at 19:43
• @JDoe you are very welcome. There are lots of free online resources for mixed models. Dimitris Rizopoulous (on here a lot) has a whole course online Jan 21, 2020 at 20:18

I will guess that students are nested in schools (i.e. each student only appears in one school), but tests are crossed with respect to both students and schools. This means that each test is taken by several different students, each student takes several different tests, each test is administered in several different schools, and each school administers several different tests.

This leads to a difficult situation that Statsmodels may not be able to handle very efficiently. If you want to try, the code below should work in theory, but may take a very long time to run with your sample size.

test_scores_df["groups"] = 1 # Put everyone in the same group.
vcf = {"student": "0 + C(student)", "test": "0 + C(test)", "school": "0 + C(school)"}
model = smf.mixedlm("score ~ 1", re_formula="1", vc_formula=vcf, groups="groups", data=test_scores_df)
result = model.fit()

• You're spot on with your assumptions on the structure. Could you explain why you chose that variance component formula, and also why you decided to place them all within the same higher level group? Isn't it failing to represent the nested structure / shouldn't school be the higher level group? Jan 21, 2020 at 16:10
• Putting everything in the same group allows the most general specification of variance/covariance structures. The only downside is that it is slower and won't scale to very large data sets. Observations within different top-level groups are independent of each other. In this example, there are "test effects", and the same tests are used in different schools (that is my understanding at least). In other words, "tests' and "schools" are crossed. If you had schools as the top level group, then students in different schools who take the same test would be independent. Jan 23, 2020 at 1:16