Linear mixed-effect model analysis via python

Question

I want to perform linear mixed-effect analyses for my research. I am trying to understand and compare the effect of 3 different intervention models on the outcome. I have 2 measures for the outcome, pre- and post-analysis. I worked with a code on Python, but when I compared it with the SPSS results, it was quite different. I was wondering whether I used the correct code.

Here are my variables DV: ACC IV: Group (3), Time(2) Covariates: sex, age, education years.

Sex, Group and time are categorical variables.

I used the code below

model = smf.mixedlm(
    "ACC ~ C(Time) * C(Group) + Age + C(Sex) + Education years", 
    data, 
    groups=data["ID"]
)
result = model.fit()

What do you think? Should I change something?

Answer 1

Your model:

smf.mixedlm(
    "ACC ~ C(Time) * C(Group) + Age + C(Sex) + Education years", 
    data, 
    groups=data["ID"]
)

...has the following features:

1. Fixed effects for Time and Group, and the interaction between them.

2. Fixed effects for Sex and Education years. However please note that Education years will cause a problem because spaces are not allowed within variable names - thus you should change the variable name to something like Education_years

3. Random intercepts for ID. This is appropriate when observations are grouped and therefore are likely to be non-independent - that is, observations for one ID are likely to be more similar to each other than observations for another ID.

The categorical variables (Group, Time, Sex) are specified using C(), which tells statsmodels to treat them as factors. This is consistent with SPSS, which also treats such variables as categorical by default. I am not very experienced with SPSS but If you want to ensure consistency between the two, double-check that Python might handle categorical variables differently in terms of reference categories and dummy coding. By default, statsmodels uses "reference" coding, where one level of a categorical variable is treated as the baseline (this is also known as "treatment" coding in some software). In SPSS, you can also set how categories are treated (e.g., reference or deviation coding).

Since this is a longitudinal study, we would expect there to be serial correlation between observations within subjects. However, this model does not account for it. Rather, it assumes that the correlation between consecutive observations is the same as the correlation between observations that are further apart. One way to handle this is by specifying random slopes for Time. Additionally, several different covariance structures should be considered, including AR(1), Toeplitz, Ante-dependence, and Unstructured. However, to the best of my knowledge, these options are not available in the mixedlm implementation from statsmodels in Python. Further details about covariance structures in mixed models can be found in these threads:

Beyond AR(1) as a covariance structure for mixed models with repeated measures

How to handle the residual covariance structure in a mixed model with repeated measurements?