Estoy tratando de usar el modelo de efectos mixtos lineales de Python statsmodels para ajustar un modelo que tiene dos intercepciones aleatorias, p. dos grupos. No puedo entender cómo inicializar el modelo para poder hacer esto.

Aquí está el ejemplo. Tengo datos que se parecen a lo siguiente (tomado deaquí):

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

Quiero hacer un modelo de efectos lineales mixtos con dos efectos aleatorios: uno para el grupo de sujetos y otro para el grupo de escenarios. Estoy tratando de hacer esto:

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

Sigo recibiendo este error:

LinAlgError: Singular matrix

Funciona bien en R. Cuando usolme4 en R con la representación basada en fórmulas, encaja perfectamente:

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

No entiendo por qué está sucediendo esto. Funciona cuando uso cualquier efecto / grupo aleatorio, p.

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

Entonces me sale:

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

Alternativamente, si lo hago:

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

Este es el resultado que obtengo:

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

No tengo idea de lo que está pasando. Siento que me falta algo fundamental en las estadísticas del problema.