(When) Does GLMM provide better predictions than logistic regression?

Question

I am trying to test statsmodels GLMM vs logistic regression (by either statsmodels or scikit-learn - see the code with a toy example below). I understand the utility of mixed models in parameter estimation, and that this model might be more truthful in terms of Akaike or some other criterion. However, I am not sure whether and when I should expect improvement in prediction accuracy, when testing on data coming from unknown groups.

Here is the code with a toy example:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

import statsmodels.api as sm
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

from scipy.special import expit

seed = 2024

"""
Comparing logistic regression with GLM
"""

#Generating mixed data
n_smpls = 1000
n_feats = 5
n_groups = 100

groups = np.repeat(np.arange(n_groups), n_smpls//n_groups)

rng = np.random.default_rng(seed=2024)
X_train = np.random.normal(size=(n_smpls, n_feats))

offs = np.random.normal(size=n_groups, scale=0.5)
offsets = np.array([offs[g] for g in groups])

x0 = 0.; b = 10.
y_train = (np.random.random_sample(n_smpls) < expit(b * (X_train[:,0]-x0-offsets))).astype(np.float64) 

#Generating test data
m_smpls = 100
xmax = 1.
xx = np.linspace(-xmax, xmax, m_smpls)
X_test = np.random.normal(size=(m_smpls, n_feats))
X_test[:,0] = xx
prob_test = expit(b * (X_test[:,0] - x0))
y_test = (np.random.random_sample(m_smpls) < prob_test).astype(np.float64) 


#--------------------------------------------------------------------------
#scaling data
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

preds = []
probs = [prob_test]
#--------------------------------------------------------------------------
#sklearn logistic regression
clf = LogisticRegression(random_state=seed, penalty='l2')
clf.fit(X_train, y_train)

y_pred = clf.predict(X_test)
preds.append(y_pred)

prob = clf.predict_proba(X_test)[:,1]
probs.append(prob)

#--------------------------------------------------------------------------
#Logistic regression by Statsmodels
clf = sm.GLM(y_train, X_train_scaled, family=sm.families.Binomial()).fit_regularized(L1_wt=.9, alpha=.01)

prob = clf.predict(X_test_scaled)
probs.append(prob)

y_pred = (prob > 0.5).astype(float)
preds.append(y_pred)

#--------------------------------------------------------------------------
#GLMM with statsmodels
Z_train = np.zeros((n_smpls, n_groups))
Z_train[np.arange(n_smpls), groups] = 1
glmm = sm.BinomialBayesMixedGLM(y_train, X_train_scaled, Z_train, np.arange(n_groups))
clf = glmm.fit_vb(scale_fe=True)

prob = clf.predict(X_test_scaled)
probs.append(prob)

y_pred = (prob > 0.5).astype(float)
preds.append(y_pred)

#--------------------------------------------------------------------------
#Analyzing results
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
for name, prob in zip(['test', 'lr', 'glm', 'glmm'], probs):
    ax.plot(xx, prob, label=name)
ax.legend()
ax.set_xlabel(r'$X_0$')
ax.set_ylabel('prob')
ax.set_title('Predicted probability')
plt.tight_layout()
plt.show()

The resulting figure is:

We see here that in terms of predicting the probability GLMM does better, but it produces the same intercept, and no improvement in the classification accuracy: one could evaluate the numbers of true/false positive and the related metrics. E.g., for true positive/negative rates one gets

[np.sum((y_pred==1)&(y_test==1))/np.sum(y_test==1) for y_pred in preds]
[0.9782608695652174, 0.9782608695652174, 0.9565217391304348]

with the results

[0.9782608695652174, 0.9782608695652174, 0.9565217391304348]
[0.9074074074074074, 0.9259259259259259, 0.9259259259259259]

Answer 1

There is no good answer to this. to some extent L2 regularised logistic regression is a subset of a GLMM model. As such you should be able to replicate any L2 regularised logistic regression by GLMM.

in your statsmodels code you are using elastic net regularisation, in sklearn code, you are using L2.

GLMM estimates the regularisation parameters from the training data, whereas in a typical ML application, you estimate the parameter using crossvalidation to maximise test prediction score.

Just as what value of the regularisation parameteer or whether L1 or L2 regularisation performs better depends on your data, so will whether using regularised glm or GLMM.