Estou tentando usar o modelo de efeitos mistos lineares do Python statsmodels para ajustar um modelo que possui duas interceptações aleatórias, por exemplo dois grupos. Não consigo descobrir como inicializar o modelo para que eu possa fazer isso.

Aqui está o exemplo. Tenho dados parecidos com os seguintes (extraídos deaqui):

subject gender  scenario    attitude    frequency
F1  F   1   pol 213.3
F1  F   1   inf 204.5
F1  F   2   pol 285.1
F1  F   2   inf 259.7
F1  F   3   pol 203.9
F1  F   3   inf 286.9
F1  F   4   pol 250.8
F1  F   4   inf 276.8

Quero criar um modelo linear de efeitos mistos com dois efeitos aleatórios - um para o grupo de sujeitos e outro para o grupo de cenários. Estou tentando fazer isso:

import statsmodels.api as sm
model = sm.MixedLM.from_formula("frequency ~ attitude + gender", data, groups=data[['subject', 'scenario']])
result = model.fit()
print result.summary()

Eu continuo recebendo esse erro:

LinAlgError: Singular matrix

Funciona bem em R. Quando usolme4 em R com a renderização baseada em fórmula, ele se encaixa perfeitamente:

politeness.model = lmer(frequency ~ attitude + gender + 
        (1|subject)  + (1|scenario), data=politeness)

Não entendo por que isso está acontecendo. Funciona quando uso qualquer efeito / grupo aleatório, por exemplo

model = sm.MixedLM.from_formula("frequency ~ attitude + gender", data, groups=data['subject'])

Então eu recebo:

                 Mixed Linear Model Regression Results
===============================================================
Model:                MixedLM   Dependent Variable:   frequency
No. Observations:     83        Method:               REML     
No. Groups:           6         Scale:                850.9456 
Min. group size:      13        Likelihood:           -393.3720
Max. group size:      14        Converged:            Yes      
Mean group size:      13.8                                     
---------------------------------------------------------------
                 Coef.   Std.Err.   z    P>|z|  [0.025   0.975]
---------------------------------------------------------------
Intercept        256.785   15.226 16.864 0.000  226.942 286.629
attitude[T.pol]  -19.415    6.407 -3.030 0.002  -31.972  -6.858
gender[T.M]     -108.325   21.064 -5.143 0.000 -149.610 -67.041
Intercept RE     603.948   23.995                              
===============================================================

Como alternativa, se eu fizer:

model = sm.MixedLM.from_formula("frequency ~ attitude + gender", data, groups=data['scenario'])

Este é o resultado que recebo:

              Mixed Linear Model Regression Results
================================================================
Model:               MixedLM    Dependent Variable:    frequency
No. Observations:    83         Method:                REML     
No. Groups:          7          Scale:                 1110.3788
Min. group size:     11         Likelihood:            -402.5003
Max. group size:     12         Converged:             Yes      
Mean group size:     11.9                                       
----------------------------------------------------------------
                 Coef.   Std.Err.    z    P>|z|  [0.025   0.975]
----------------------------------------------------------------
Intercept        256.892    8.120  31.637 0.000  240.977 272.807
attitude[T.pol]  -19.807    7.319  -2.706 0.007  -34.153  -5.462
gender[T.M]     -108.603    7.319 -14.838 0.000 -122.948 -94.257
Intercept RE     182.718    5.502                               
================================================================

Eu não tenho ideia do que está acontecendo. Sinto que estou perdendo algo fundamental nas estatísticas do problema.