How to have multiple groups in Python statsmodels linear mixed effects model?
<p>I am trying to use the Python statsmodels linear mixed effects model to fit a model that has two random intercepts, e.g. two groups.  I cannot figure out how to initialize the model so that I can do this. </p>

<p>Here's the example. I have data that looks like the following (taken from <a href="http://www.bodowinter.com/tutorial/politeness_data.csv" rel="nofollow">here</a>): </p>

<pre><code>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
</code></pre>

<p>I want to make a linear mixed effects model with two random effects -- one for the subject group and one for the scenario group.  I am trying to do this: </p>

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

<p>I keep getting this error: </p>

<pre><code>LinAlgError: Singular matrix
</code></pre>

<p>It works fine in R.  When I use <code>lme4</code> in R with the formula-based rendering it fits just fine: </p>

<pre><code>politeness.model = lmer(frequency ~ attitude + gender + 
        (1|subject)  + (1|scenario), data=politeness)
</code></pre>

<p>I don't understand why this is happening. It works when I use any one random effect/group, e.g.</p>

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

<p>Then I get: </p>

<pre><code>                 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&gt;|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                              
===============================================================
</code></pre>

<p>Alternatively, if I do: </p>

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

<p>This is the result I get: </p>

<pre><code>              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&gt;|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                               
================================================================
</code></pre>

<p>I have no idea what's going on. I feel like I am missing something foundational in the statistics of the problem. </p>

<p>You are trying to fit a model with <em>crossed random effects</em>, i.e., you want to allow for consistent variation among subjects across scenarios as well as consistent variation among scenarios across subjects. You <em>can</em> use multiple random-effects terms in statsmodels, but they must be nested. Fitting crossed (as opposed to nested) random effects requires more sophisticated algorithms, and indeed the <a href="http://statsmodels.sourceforge.net/devel/mixed_linear.html" rel="nofollow">statsmodels documentation</a> says (as of 25 Aug 2016, emphasis added):</p>

<blockquote>
  <p>Some limitations of the current implementation are that it does not support structure more complex on the residual errors (they are always homoscedastic), and <strong>it does not support crossed random effects</strong>. We hope to implement these features for the next release.</p>
</blockquote>

<p>As far as I can see, your choices are (1) fall back to a nested model (i.e. fit the model as though either scenario is nested within subject or <em>vice versa</em> - or try both and see if the difference matters); (2) fall back to <code>lme4</code>, either within R or via <a href="http://rpy2.bitbucket.org/" rel="nofollow">rpy2</a>.</p>

<p>As always, you're entitled to a full refund of the money you paid to use statsmodels ...</p>

python statsmodels lme4 mixed-models