Which is the correct way to specify crossed effects in a mixed linear model in statsmodels?

Question

I have a question about crossed effects linear mixed models in statsmodels. Specifically, I can see two ways of representing my data and I'm not sure which is appropriate. Any advice appreciated!

My data are as below. I wish to ascertain whether the objective quality of a book ('good' or 'bad') predicts the score the book is assigned. Titles are nested within the quality variable, but titles and raters are crossed. (This is fake data so I'm not worried about models converging.)

rater   title   quality score
john    book_1  good    0.600833333
frank   book_2  bad     0.683020833
emma    book_3  good    0.653645833
john    book_4  bad     0.6528125
frank   book_5  good    0.6040625
emma    book_1  good    0.600833333
john    book_2  bad     0.522
frank   book_3  good    0.600833333
emma    book_4  bad     0.619464286
john    book_5  good    0.600833333
frank   book_1  good    0.57125
emma    book_2  bad     0.6296875
john    book_3  good    0.607205882
frank   book_4  bad     0.61203125
emma    book_5  good    0.600833333

One way to analyse this data is to take quality as my independent variable, score as my dependent variable, rater as my grouping variable, and use variance components to capture the crossed effects on title. This gives:

import statsmodels.api as sm
import statsmodels.fomula.api as smf

md = smf.mixedlm('score ~ quality', vc_formula = {"title":"0 + title"}, groups = data['rater'], data = data).fit().summary()

Model summary:

           Mixed Linear Model Regression Results
===========================================================
Model:              MixedLM   Dependent Variable:   score  
No. Observations:   15        Method:               REML   
No. Groups:         3         Scale:                0.0007 
Min. group size:    5         Log-Likelihood:       22.1997
Max. group size:    5         Converged:            Yes    
Mean group size:    5.0                                    
-----------------------------------------------------------
                Coef.  Std.Err.    z    P>|z| [0.025 0.975]
-----------------------------------------------------------
Intercept        0.620    0.001 841.098 0.000  0.618  0.621
quality[T.good] -0.015    0.013  -1.158 0.247 -0.041  0.011
title Var        0.001                                     
===========================================================

This intuitively seems to me to be the correct approach. It gives me a p-value and coefficient for my IV and accounts for the crossed effects.

However, I've been advised elsewhere that crossed effects like this should be specified by treating the dataset as one group and specifying variation entirely using variance components. Thus:

data['groups'] = 1

md = smf.mixedlm('score ~ 1', vc_formula = {"rater":"0 + rater", "title":"0 + title", "quality":"0 + quality"}, groups = data['groups'], data = data).fit().summary()

Yielding:

        Mixed Linear Model Regression Results
=====================================================
Model:            MixedLM Dependent Variable: score  
No. Observations: 15      Method:             REML   
No. Groups:       1       Scale:              0.0013 
Min. group size:  15      Log-Likelihood:     24.4023
Max. group size:  15      Converged:          No     
Mean group size:  15.0                               
-----------------------------------------------------
            Coef. Std.Err.   z    P>|z| [0.025 0.975]
-----------------------------------------------------
Intercept   0.612    0.017 35.011 0.000  0.577  0.646
quality Var 0.000                                    
rater Var   0.000    0.020                           
title Var   0.000                                    
=====================================================

This model offers me no p-value, different coefficients, and different model test statistics all around. Now, I'm either wrong in my use of both models, or I'm wrong in my use of one of them. Can anyone advise me which is the case? Thanks.

score 3 · Accepted Answer · answered Mar 17 '20 at 11:42

You are comparing two totally different models, which gives different interpretation.

If you are interested in how the quality affects the school, while putting the other covariates as a random intercept, the first model is correct. I quickly checked and found this post mentioning the cross effect, and if you run the model with all your random intercepts in the variance component, you get roughly the same result:

data['group'] = 1
md2 = smf.mixedlm('score ~ quality', 
                  vc_formula = {"title":"0 + title","rater":"0 + rater"}, 
                  groups = data['group'], data = data).fit().summary()

Model:  MixedLM Dependent Variable: score
No. Observations:   15  Method: REML
No. Groups: 1   Scale:  0.0014
Min. group size:    15  Log-Likelihood: 22.0951
Max. group size:    15  Converged:  No
Mean group size:    15.0        
Coef.   Std.Err.    z   P>|z|   [0.025  0.975]
Intercept   0.620   0.016   38.313  0.000   0.588   0.652
quality[T.good] -0.015  0.021   -0.736  0.462   -0.056  0.026
rater Var   0.000   0.030               
title Var   0.000

It's a matter of whether the 'rater' appears in this table. The VC is also for specifying more complex models e.g variable slope, so if it is pure random intercept, I think you can use a combination of group and VC.

In the second model, you are modeling score with only an intercept, independent of quality, which doesn't makes sense if you are interested in quality.

Lastly we can check with the results in R with the following:

df = structure(list(rater = structure(c(3L, 2L, 1L, 3L, 2L, 1L, 3L, 
2L, 1L, 3L, 2L, 1L, 3L, 2L, 1L), .Label = c("emma", "frank", 
"john"), class = "factor"), title = structure(c(1L, 2L, 3L, 4L, 
5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), .Label = c("book_1", 
"book_2", "book_3", "book_4", "book_5"), class = "factor"), quality = structure(c(2L, 
1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L), .Label = c("bad", 
"good"), class = "factor"), score = c(0.600833333, 0.683020833, 
0.653645833, 0.6528125, 0.6040625, 0.600833333, 0.522, 0.600833333, 
0.619464286, 0.600833333, 0.57125, 0.6296875, 0.607205882, 0.61203125, 
0.600833333)), class = "data.frame", row.names = c(NA, -15L))

library(lme4)
summary(lmer(score ~ quality + (1|rater) + (1|title),data=df))
boundary (singular) fit: see ?isSingular
Linear mixed model fit by REML ['lmerMod']
Formula: score ~ quality + (1 | rater) + (1 | title)
   Data: df

REML criterion at convergence: -44.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.60015 -0.09695 -0.09695  0.16712  1.67924 

Random effects:
 Groups   Name        Variance Std.Dev.
 title    (Intercept) 0.000000 0.00000 
 rater    (Intercept) 0.000000 0.00000 
 Residual             0.001416 0.03763 
Number of obs: 15, groups:  title, 5; rater, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.61984    0.01536  40.351
qualitygood -0.01535    0.01983  -0.774

Correlation of Fixed Effects:
            (Intr)
qualitygood -0.775
convergence code: 0
boundary (singular) fit: see ?isSingular

More or less similar to your first model, but because of your data, it's not easy to get an estimate of your random effects.

Thanks for this post, in which you went to a lot of effort to answer my question. It really clarifies the issue for me. Using the combination of group and VC seems optimal in my case so. I really appreciate this! — Lodore66, Mar 17 '20 at 13:00
By the way, you note that my first model (VC plus group) gives roughly the same results as specifying everything with just VC. This is true on most test statistics bar the p value (.247 vs .462). I appreciate that p values are a bit of a distraction in mixed models, but do you have any thoughts in which is the more reliable estimation of significance? — Lodore66, Mar 17 '20 at 13:11
Hey @Lodore66, glad it was useful for you. Which statsmodel version are you on? I am on 0.11.1 and when I ran your code, I got -0.015 for coefficient quality and 0.020 as its standard error.. which is very similar to that above — StupidWolf, Mar 17 '20 at 13:21
If you look at the standard error, it's not that big a difference, so until you see something where the standard error shows hugh difference, I will not be so worried about it — StupidWolf, Mar 17 '20 at 13:22
Thanks for getting back. I'm on version 0.11, but I also got those values too for coefficient and SE; the difference was on the p value? — Lodore66, Mar 17 '20 at 13:24
The p-value is calculated based on a t-statistic which is the coefficient over the standard error (https://stats.stackexchange.com/questions/337070/compute-standard-error-from-beta-p-value-sample-size-and-the-number-of-regres) , hence my guess is that the standard error is slightly different but not much — StupidWolf, Mar 17 '20 at 13:27

Which is the correct way to specify crossed effects in a mixed linear model in statsmodels?

1 Answers1