Linear Mixed Effects Models

In [1]:
%matplotlib inline

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools
/usr/lib/python3/dist-packages/scipy/_lib/_util.py:18: RuntimeWarning: invalid value encountered in multiply
  out = np.ones(shape, dtype=bool) * value
In [2]:
%load_ext rpy2.ipython
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-2-f4f84eb1b14e> in <module>()
----> 1 get_ipython().magic('load_ext rpy2.ipython')

/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py in magic(self, arg_s)
   2158         magic_name, _, magic_arg_s = arg_s.partition(' ')
   2159         magic_name = magic_name.lstrip(prefilter.ESC_MAGIC)
-> 2160         return self.run_line_magic(magic_name, magic_arg_s)
   2161 
   2162     #-------------------------------------------------------------------------

/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line)
   2079                 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
   2080             with self.builtin_trap:
-> 2081                 result = fn(*args,**kwargs)
   2082             return result
   2083 

<decorator-gen-63> in load_ext(self, module_str)

/usr/lib/python3/dist-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
    186     # but it's overkill for just that one bit of state.
    187     def magic_deco(arg):
--> 188         call = lambda f, *a, **k: f(*a, **k)
    189 
    190         if callable(arg):

/usr/lib/python3/dist-packages/IPython/core/magics/extension.py in load_ext(self, module_str)
     35         if not module_str:
     36             raise UsageError('Missing module name.')
---> 37         res = self.shell.extension_manager.load_extension(module_str)
     38 
     39         if res == 'already loaded':

/usr/lib/python3/dist-packages/IPython/core/extensions.py in load_extension(self, module_str)
     81             if module_str not in sys.modules:
     82                 with prepended_to_syspath(self.ipython_extension_dir):
---> 83                     __import__(module_str)
     84             mod = sys.modules[module_str]
     85             if self._call_load_ipython_extension(mod):

ModuleNotFoundError: No module named 'rpy2'
In [3]:
%R library(lme4)
UsageError: Line magic function `%R` not found.

Comparing R lmer to Statsmodels MixedLM

The Statsmodels imputation of linear mixed models (MixedLM) closely follows the approach outlined in Lindstrom and Bates (JASA 1988). This is also the approach followed in the R package LME4. Other packages such as Stata, SAS, etc. should also be consistent with this approach, as the basic techniques in this area are mostly mature.

Here we show how linear mixed models can be fit using the MixedLM procedure in Statsmodels. Results from R (LME4) are included for comparison.

Here are our import statements:

Growth curves of pigs

These are longitudinal data from a factorial experiment. The outcome variable is the weight of each pig, and the only predictor variable we will use here is "time". First we fit a model that expresses the mean weight as a linear function of time, with a random intercept for each pig. The model is specified using formulas. Since the random effects structure is not specified, the default random effects structure (a random intercept for each group) is automatically used.

In [4]:
data = sm.datasets.get_rdataset('dietox', 'geepack', cache=True).data
md = smf.mixedlm("Weight ~ Time", data, groups=data["Pig"])
mdf = md.fit()
print(mdf.summary())
         Mixed Linear Model Regression Results
========================================================
Model:            MixedLM Dependent Variable: Weight    
No. Observations: 861     Method:             REML      
No. Groups:       72      Scale:              11.3669   
Min. group size:  11      Likelihood:         -2404.7753
Max. group size:  12      Converged:          Yes       
Mean group size:  12.0                                  
--------------------------------------------------------
             Coef.  Std.Err.    z    P>|z| [0.025 0.975]
--------------------------------------------------------
Intercept    15.724    0.788  19.952 0.000 14.179 17.268
Time          6.943    0.033 207.939 0.000  6.877  7.008
groups RE    40.394    2.149                            
========================================================

/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/regression/mixed_linear_model.py:2366: RuntimeWarning: invalid value encountered in multiply
  sdf = np.nan * np.ones((self.k_fe + self.k_re2 + self.k_vc, 6))
/usr/lib/python3/dist-packages/scipy/_lib/_util.py:18: RuntimeWarning: invalid value encountered in multiply
  out = np.ones(shape, dtype=bool) * value

Here is the same model fit in R using LMER:

In [5]:
%%R 
data(dietox, package='geepack')
UsageError: Cell magic `%%R` not found.
In [6]:
%R print(summary(lmer('Weight ~ Time + (1|Pig)', data=dietox)))
UsageError: Line magic function `%R` not found.

Note that in the Statsmodels summary of results, the fixed effects and random effects parameter estimates are shown in a single table. The random effect for animal is labeled "Intercept RE" in the Statmodels output above. In the LME4 output, this effect is the pig intercept under the random effects section.

There has been a lot of debate about whether the standard errors for random effect variance and covariance parameters are useful. In LME4, these standard errors are not displayed, because the authors of the package believe they are not very informative. While there is good reason to question their utility, we elected to include the standard errors in the summary table, but do not show the corresponding Wald confidence intervals.

Next we fit a model with two random effects for each animal: a random intercept, and a random slope (with respect to time). This means that each pig may have a different baseline weight, as well as growing at a different rate. The formula specifies that "Time" is a covariate with a random coefficient. By default, formulas always include an intercept (which could be suppressed here using "0 + Time" as the formula).

In [7]:
md = smf.mixedlm("Weight ~ Time", data, groups=data["Pig"], re_formula="~Time")
mdf = md.fit()
print(mdf.summary())
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/regression/mixed_linear_model.py:2001: ConvergenceWarning: Gradient optimization failed.
  warnings.warn(msg, ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/regression/mixed_linear_model.py:2366: RuntimeWarning: invalid value encountered in multiply
  sdf = np.nan * np.ones((self.k_fe + self.k_re2 + self.k_vc, 6))
/usr/lib/python3/dist-packages/scipy/_lib/_util.py:18: RuntimeWarning: invalid value encountered in multiply
  out = np.ones(shape, dtype=bool) * value
              Mixed Linear Model Regression Results
=================================================================
Model:               MixedLM    Dependent Variable:    Weight    
No. Observations:    861        Method:                REML      
No. Groups:          72         Scale:                 5.7891    
Min. group size:     11         Likelihood:            -2220.3890
Max. group size:     12         Converged:             No        
Mean group size:     12.0                                        
-----------------------------------------------------------------
                       Coef.  Std.Err.   z    P>|z| [0.025 0.975]
-----------------------------------------------------------------
Intercept              15.739    0.672 23.438 0.000 14.423 17.055
Time                    6.939    0.085 81.326 0.000  6.772  7.106
Intercept RE           30.266    4.271                           
Intercept RE x Time RE  0.746    0.304                           
Time RE                 0.483    0.046                           
=================================================================

Here is the same model fit using LMER in R:

In [8]:
%R print(summary(lmer("Weight ~ Time + (1 + Time | Pig)", data=dietox)))
UsageError: Line magic function `%R` not found.

The random intercept and random slope are only weakly correlated $(0.294 / \sqrt{19.493 * 0.416} \approx 0.1)$. So next we fit a model in which the two random effects are constrained to be uncorrelated:

In [9]:
.294 / (19.493 * .416)**.5
Out[9]:
0.10324316832591753
In [10]:
md = smf.mixedlm("Weight ~ Time", data, groups=data["Pig"],
                  re_formula="~Time")
free = sm.regression.mixed_linear_model.MixedLMParams.from_components(np.ones(2), 
                                                                      np.eye(2))

mdf = md.fit(free=free)
print(mdf.summary())
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
              Mixed Linear Model Regression Results
=================================================================
Model:               MixedLM    Dependent Variable:    Weight    
No. Observations:    861        Method:                REML      
No. Groups:          72         Scale:                 5.8015    
Min. group size:     11         Likelihood:            -2220.0996
Max. group size:     12         Converged:             No        
Mean group size:     12.0                                        
-----------------------------------------------------------------
                       Coef.  Std.Err.   z    P>|z| [0.025 0.975]
-----------------------------------------------------------------
Intercept              15.739    0.672 23.416 0.000 14.421 17.056
Time                    6.939    0.084 83.012 0.000  6.775  7.103
Intercept RE           30.322    4.025                           
Intercept RE x Time RE  0.000    0.000                           
Time RE                 0.462    0.040                           
=================================================================

/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/base/model.py:496: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals
  "Check mle_retvals", ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/regression/mixed_linear_model.py:2001: ConvergenceWarning: Gradient optimization failed.
  warnings.warn(msg, ConvergenceWarning)
/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/regression/mixed_linear_model.py:2366: RuntimeWarning: invalid value encountered in multiply
  sdf = np.nan * np.ones((self.k_fe + self.k_re2 + self.k_vc, 6))
/usr/lib/python3/dist-packages/scipy/_lib/_util.py:18: RuntimeWarning: invalid value encountered in multiply
  out = np.ones(shape, dtype=bool) * value

The likelihood drops by 0.3 when we fix the correlation parameter to 0. Comparing 2 x 0.3 = 0.6 to the chi^2 1 df reference distribution suggests that the data are very consistent with a model in which this parameter is equal to 0.

Here is the same model fit using LMER in R (note that here R is reporting the REML criterion instead of the likelihood, where the REML criterion is twice the log likeihood):

In [11]:
%R print(summary(lmer("Weight ~ Time + (1 | Pig) + (0 + Time | Pig)", data=dietox)))
UsageError: Line magic function `%R` not found.

Sitka growth data

This is one of the example data sets provided in the LMER R library. The outcome variable is the size of the tree, and the covariate used here is a time value. The data are grouped by tree.

In [12]:
data = sm.datasets.get_rdataset("Sitka", "MASS", cache=True).data
endog = data["size"]
data["Intercept"] = 1
exog = data[["Intercept", "Time"]]
---------------------------------------------------------------------------
gaierror                                  Traceback (most recent call last)
/usr/lib/python3.7/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
   1316                 h.request(req.get_method(), req.selector, req.data, headers,
-> 1317                           encode_chunked=req.has_header('Transfer-encoding'))
   1318             except OSError as err: # timeout error

/usr/lib/python3.7/http/client.py in request(self, method, url, body, headers, encode_chunked)
   1243         """Send a complete request to the server."""
-> 1244         self._send_request(method, url, body, headers, encode_chunked)
   1245 

/usr/lib/python3.7/http/client.py in _send_request(self, method, url, body, headers, encode_chunked)
   1289             body = _encode(body, 'body')
-> 1290         self.endheaders(body, encode_chunked=encode_chunked)
   1291 

/usr/lib/python3.7/http/client.py in endheaders(self, message_body, encode_chunked)
   1238             raise CannotSendHeader()
-> 1239         self._send_output(message_body, encode_chunked=encode_chunked)
   1240 

/usr/lib/python3.7/http/client.py in _send_output(self, message_body, encode_chunked)
   1025         del self._buffer[:]
-> 1026         self.send(msg)
   1027 

/usr/lib/python3.7/http/client.py in send(self, data)
    965             if self.auto_open:
--> 966                 self.connect()
    967             else:

/usr/lib/python3.7/http/client.py in connect(self)
   1398 
-> 1399             super().connect()
   1400 

/usr/lib/python3.7/http/client.py in connect(self)
    937         self.sock = self._create_connection(
--> 938             (self.host,self.port), self.timeout, self.source_address)
    939         self.sock.setsockopt(socket.IPPROTO_TCP, socket.TCP_NODELAY, 1)

/usr/lib/python3.7/socket.py in create_connection(address, timeout, source_address)
    706     err = None
--> 707     for res in getaddrinfo(host, port, 0, SOCK_STREAM):
    708         af, socktype, proto, canonname, sa = res

/usr/lib/python3.7/socket.py in getaddrinfo(host, port, family, type, proto, flags)
    747     addrlist = []
--> 748     for res in _socket.getaddrinfo(host, port, family, type, proto, flags):
    749         af, socktype, proto, canonname, sa = res

gaierror: [Errno -3] Temporary failure in name resolution

During handling of the above exception, another exception occurred:

URLError                                  Traceback (most recent call last)
<ipython-input-12-567a79faedbd> in <module>()
----> 1 data = sm.datasets.get_rdataset("Sitka", "MASS", cache=True).data
      2 endog = data["size"]
      3 data["Intercept"] = 1
      4 exog = data[["Intercept", "Time"]]

/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/datasets/utils.py in get_rdataset(dataname, package, cache)
    288                      "master/doc/"+package+"/rst/")
    289     cache = _get_cache(cache)
--> 290     data, from_cache = _get_data(data_base_url, dataname, cache)
    291     data = read_csv(data, index_col=0)
    292     data = _maybe_reset_index(data)

/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/datasets/utils.py in _get_data(base_url, dataname, cache, extension)
    219     url = base_url + (dataname + ".%s") % extension
    220     try:
--> 221         data, from_cache = _urlopen_cached(url, cache)
    222     except HTTPError as err:
    223         if '404' in str(err):

/build/statsmodels-0.8.0/.pybuild/cpython3_3.7_statsmodels/build/statsmodels/datasets/utils.py in _urlopen_cached(url, cache)
    210     # not using the cache or didn't find it in cache
    211     if not from_cache:
--> 212         data = urlopen(url).read()
    213         if cache is not None:  # then put it in the cache
    214             _cache_it(data, cache_path)

/usr/lib/python3.7/urllib/request.py in urlopen(url, data, timeout, cafile, capath, cadefault, context)
    220     else:
    221         opener = _opener
--> 222     return opener.open(url, data, timeout)
    223 
    224 def install_opener(opener):

/usr/lib/python3.7/urllib/request.py in open(self, fullurl, data, timeout)
    523             req = meth(req)
    524 
--> 525         response = self._open(req, data)
    526 
    527         # post-process response

/usr/lib/python3.7/urllib/request.py in _open(self, req, data)
    541         protocol = req.type
    542         result = self._call_chain(self.handle_open, protocol, protocol +
--> 543                                   '_open', req)
    544         if result:
    545             return result

/usr/lib/python3.7/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
    501         for handler in handlers:
    502             func = getattr(handler, meth_name)
--> 503             result = func(*args)
    504             if result is not None:
    505                 return result

/usr/lib/python3.7/urllib/request.py in https_open(self, req)
   1358         def https_open(self, req):
   1359             return self.do_open(http.client.HTTPSConnection, req,
-> 1360                 context=self._context, check_hostname=self._check_hostname)
   1361 
   1362         https_request = AbstractHTTPHandler.do_request_

/usr/lib/python3.7/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
   1317                           encode_chunked=req.has_header('Transfer-encoding'))
   1318             except OSError as err: # timeout error
-> 1319                 raise URLError(err)
   1320             r = h.getresponse()
   1321         except:

URLError: <urlopen error [Errno -3] Temporary failure in name resolution>

Here is the statsmodels LME fit for a basic model with a random intercept. We are passing the endog and exog data directly to the LME init function as arrays. Also note that endog_re is specified explicitly in argument 4 as a random intercept (although this would also be the default if it were not specified).

In [13]:
md = sm.MixedLM(endog, exog, groups=data["tree"], exog_re=exog["Intercept"])
mdf = md.fit()
print(mdf.summary())
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-13-e6b2b322447f> in <module>()
----> 1 md = sm.MixedLM(endog, exog, groups=data["tree"], exog_re=exog["Intercept"])
      2 mdf = md.fit()
      3 print(mdf.summary())

NameError: name 'endog' is not defined

Here is the same model fit in R using LMER:

In [14]:
%%R
data(Sitka, package="MASS")
print(summary(lmer("size ~ Time + (1 | tree)", data=Sitka)))
UsageError: Cell magic `%%R` not found.

We can now try to add a random slope. We start with R this time. From the code and output below we see that the REML estimate of the variance of the random slope is nearly zero.

In [15]:
%R print(summary(lmer("size ~ Time + (1 + Time | tree)", data=Sitka)))
UsageError: Line magic function `%R` not found.

If we run this in statsmodels LME with defaults, we see that the variance estimate is indeed very small, which leads to a warning about the solution being on the boundary of the parameter space. The regression slopes agree very well with R, but the likelihood value is much higher than that returned by R.

In [16]:
exog_re = exog.copy()
md = sm.MixedLM(endog, exog, data["tree"], exog_re)
mdf = md.fit()
print(mdf.summary())
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-16-cdffaa9ce318> in <module>()
----> 1 exog_re = exog.copy()
      2 md = sm.MixedLM(endog, exog, data["tree"], exog_re)
      3 mdf = md.fit()
      4 print(mdf.summary())

NameError: name 'exog' is not defined

We can further explore the random effects struture by constructing plots of the profile likelihoods. We start with the random intercept, generating a plot of the profile likelihood from 0.1 units below to 0.1 units above the MLE. Since each optimization inside the profile likelihood generates a warning (due to the random slope variance being close to zero), we turn off the warnings here.

In [17]:
import warnings

with warnings.catch_warnings():
    warnings.filterwarnings("ignore")
    likev = mdf.profile_re(0, 're', dist_low=0.1, dist_high=0.1)

Here is a plot of the profile likelihood function. We multiply the log-likelihood difference by 2 to obtain the usual $\chi^2$ reference distribution with 1 degree of freedom.

In [18]:
import matplotlib.pyplot as plt
In [19]:
plt.figure(figsize=(10,8))
plt.plot(likev[:,0], 2*likev[:,1])
plt.xlabel("Variance of random slope", size=17)
plt.ylabel("-2 times profile log likelihood", size=17)
Out[19]:
Text(0, 0.5, '-2 times profile log likelihood')

Here is a plot of the profile likelihood function. The profile likelihood plot shows that the MLE of the random slope variance parameter is a very small positive number, and that there is low uncertainty in this estimate.

In [20]:
re = mdf.cov_re.iloc[1, 1]
likev = mdf.profile_re(1, 're', dist_low=.5*re, dist_high=0.8*re)

plt.figure(figsize=(10, 8))
plt.plot(likev[:,0], 2*likev[:,1])
plt.xlabel("Variance of random slope", size=17)
plt.ylabel("-2 times profile log likelihood", size=17)
Out[20]:
Text(0, 0.5, '-2 times profile log likelihood')
In [ ]: