Как правильно указать скрещенные эффекты в смешанной линейной модели в statsmodels?

Question

У меня вопрос о скрещенных эффектах линейных смешанных моделей в statsmodels. В частности, я вижу два способа представления своих данных, и я не уверен, какой из них подходит. Любой совет приветствуется!

Мои данные, как показано ниже. Я sh, чтобы установить, предсказывает ли объективное качество книги («хорошо» или «плохо») оценку, присвоенную книге. Заголовки вложены в переменную качества, но заголовки и оценщики пересекаются. (Это фальшивые данные, поэтому я не беспокоюсь о сходстве моделей.)

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

Один из способов анализа этих данных - это взять качество в качестве моей независимой переменной, оценить в качестве моей зависимой переменной, оценить в качестве моей группировки. переменная, и использовать компоненты дисперсии для захвата скрещенных эффектов на заголовок. Это дает:

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()

Сводка модели:

           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                                     
===========================================================

Мне кажется, что это правильный подход. Это дает мне p-значение и коэффициент для моего IV и учитывает скрещенные эффекты.

Однако в других местах мне советовали, что подобные перекрестные эффекты следует указывать, рассматривая набор данных как одну группу и указав вариацию полностью, используя компоненты дисперсии. Таким образом:

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()

Выход:

        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                                    
=====================================================

Эта модель не предлагает мне никакого p-значения, разных коэффициентов и разных статистических тестов модели. Теперь я либо ошибаюсь в использовании обеих моделей, либо ошибаюсь в использовании одной из них. Может кто-нибудь посоветовать мне, в чем дело? Благодаря.

Answer 1

Вы сравниваете две совершенно разные модели, которые дают разную интерпретацию.

Если вас интересует, как качество влияет на школу, а другие ковариаты ставятся в качестве случайного перехвата, первая модель верна. Я быстро проверил и нашел этот пост , в котором упоминается перекрестный эффект, и если вы запустите модель со всеми вашими случайными перехватами в компоненте дисперсии, вы получите примерно тот же результат:

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           

Вопрос в том, появляется ли «рейтинг» в этой таблице. V C также предназначен для определения более сложных моделей, например, с переменным наклоном, поэтому, если это чистый случайный перехват, я думаю, вы можете использовать комбинацию group и V C.

Во второй модели вы моделируете счет только с перехватом, независимым от качества, что не имеет смысла, если вы заинтересованы в качестве.

Наконец, мы можем проверить результаты в R со следующим:

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

Более или менее похож на вашу первую модель, но из-за ваших данных не легко получить оценку ваших случайных эффектов.