4.7 Clustered data: mixed models

When data come in clusters (pupils in schools, measurements of the same person, years of the same firm), cases in a cluster are alike and carry less information than their number suggests. A mixed model gives each cluster its own level, but pulls the levels towards the average, more so for clusters with little data.

import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from matplotlib.axes import Axes

from statly import style

style.use()
warnings.filterwarnings("ignore")  # mixed models warn about their numerical optimisation

grunfeld = sm.datasets.grunfeld.load_pandas().data  # investment of 11 US firms, 1935 to 1954
grunfeld["log_invest"] = np.log(grunfeld.invest)
grunfeld["log_value"] = np.log(grunfeld.value)
FIRMS = grunfeld.firm.unique()
FIRM_COLORS = dict(zip(FIRMS, plt.cm.tab20(np.linspace(0, 1, len(FIRMS)))))

4.7.1 Clusters

Eleven large US firms, each followed for 20 years: how much they invested, and what the firm was worth on the market. Do firms invest more when they are worth more? There are 220 rows, but they are not 220 independent observations. A firm’s investment this year is much like its investment last year, and General Motors is always much bigger than Diamond Match.

pooled = smf.ols("log_invest ~ log_value", grunfeld).fit()
grunfeld["pooled_residual"] = pooled.resid
order = grunfeld.groupby("firm").pooled_residual.median().sort_values().index


def plot_firms(ax: Axes) -> None:
    """All firm-years coloured by firm, with the pooled line."""
    for firm, part in grunfeld.groupby("firm"):
        ax.scatter(part.log_value, part.log_invest, s=10, color=FIRM_COLORS[firm], alpha=0.8)
    grid = np.array([grunfeld.log_value.min(), grunfeld.log_value.max()])
    ax.plot(
        grid, pooled.params.Intercept + pooled.params.log_value * grid, color=style.INK, lw=2
    )
    ax.set(title="Eleven firms", xlabel="log market value", ylabel="log investment")


def plot_residuals_by_firm(ax: Axes) -> None:
    """Residuals of the pooled line, one column of dots per firm."""
    for i, firm in enumerate(order):
        part = grunfeld[grunfeld.firm == firm]
        ax.scatter(
            np.full(len(part), i),
            part.pooled_residual,
            s=10,
            color=FIRM_COLORS[firm],
            alpha=0.8,
        )
    ax.axhline(0, color=style.INK, lw=1)
    ax.set_xticks(range(len(order)), order, rotation=60, ha="right", fontsize=8)
    ax.set(title="Residuals by firm", ylabel="residual")


fig, axes = style.subplots(2, height=3.6, width=10)
plot_firms(axes[0])
plot_residuals_by_firm(axes[1])
fig.tight_layout()
plt.show()

Figure 4.7.1: Left: log investment against log market value, 20 years for each of 11 firms (one colour per firm), with one line through all 220 points. Right: the residuals of that line, grouped by firm. Each firm’s residuals sit mostly above or mostly below zero: the firms differ in a way the line does not capture, and the 20 years of one firm are not 20 independent pieces of information.

The share of the variation that lies between clusters is the intraclass correlation (ICC): the correlation between two observations from the same cluster. Here it is about 0.77. With an ICC that high, 20 years of one firm are worth only a little more than one year.

4.7.2 Three ways to handle clusters

• Complete pooling ignores the firms: one line through everything. Its coefficient may be fine, but its standard errors treat 220 rows as independent and are far too small.
• No pooling fits each firm separately, or gives each firm its own intercept (fixed effects: one indicator per firm, Chapter 4.2). It respects the firms but learns nothing about one firm from the others, and cannot say anything about a new firm.
• Partial pooling is the middle way. A mixed model gives each firm its own intercept, but treats these intercepts as drawn from a common normal distribution, estimated from the data. Firms with little or noisy data are pulled more towards the average: they borrow strength from the others.

The random intercept model writes this as

$$y_{ij} = \beta_0 + u_j + \beta_1 x_{ij} + \varepsilon_{ij}, \qquad u_j \sim N(0, \tau^2), \quad \varepsilon_{ij} \sim N(0, \sigma^2)$$

for year i of firm j. $\beta_0$ and $\beta_1$ are fixed effects, the same for every firm; $u_j$ is the random effect of firm j, its own shift. The ICC is $\tau^2 / (\tau^2 + \sigma^2)$.

mixed = smf.mixedlm("log_invest ~ log_value", grunfeld, groups=grunfeld.firm).fit()
centre = grunfeld.log_value.mean()
fixed = smf.ols(
    "log_invest ~ log_value + C(firm)", grunfeld
).fit()  # each firm its own level, one slope
own_level = {
    firm: fixed.predict(pd.DataFrame({"firm": [firm], "log_value": [centre]})).iloc[0]
    for firm in FIRMS
}
pooled_level = {
    firm: mixed.params.Intercept
    + mixed.random_effects[firm].iloc[0]
    + mixed.params.log_value * centre
    for firm in FIRMS
}
average_level = mixed.params.Intercept + mixed.params.log_value * centre


def plot_shrinkage(ax: Axes) -> None:
    """Separate-fit and mixed-model levels per firm, joined by arrows."""
    firms = sorted(FIRMS, key=own_level.get)
    for i, firm in enumerate(firms):
        ax.annotate(
            "",
            xy=(pooled_level[firm], i),
            xytext=(own_level[firm], i),
            arrowprops=dict(arrowstyle="->", color=style.INK, lw=1),
        )
        ax.plot(own_level[firm], i, "o", ms=8, mfc="none", mec=FIRM_COLORS[firm], mew=1.6)
        ax.plot(pooled_level[firm], i, "o", ms=7, color=FIRM_COLORS[firm])
    ax.axvline(average_level, color=style.INK, ls="--", lw=1)
    ax.set_yticks(range(len(firms)), firms, fontsize=8.5)
    ax.set(xlabel="log investment at the average market value")


fig, ax = plt.subplots(figsize=(8, 3.8))
plot_shrinkage(ax)
plt.show()

Figure 4.7.2: Each firm’s own level of investment at the average market value, estimated without pooling (open circles: fixed effects, one indicator per firm) and with partial pooling (filled dots: the mixed model). The arrows show the pull towards the overall average (dashed). With 20 years per firm, and firms that differ far more than years do, the pull is small; with few observations per cluster it would be large.

4.7.3 Why it matters: the standard errors

The slope of log investment on log market value is close to 1 whichever way it is estimated: firms invest roughly in proportion to their value. What changes is how sure we can be:

codes = grunfeld.firm.astype("category").cat.codes
cluster = pooled.get_robustcov_results("cluster", groups=codes)
ESTIMATES = [
    ("pooled, ordinary SE", pooled.params.log_value, pooled.bse.log_value),
    ("pooled, cluster-robust SE", cluster.params[1], cluster.bse[1]),
    ("firm fixed effects", fixed.params.log_value, fixed.bse.log_value),
    ("mixed model", mixed.params.log_value, mixed.bse.log_value),
]


def plot_estimates(ax: Axes) -> None:
    """Slope and 95% interval for each approach."""
    for i, (name, estimate, se) in enumerate(ESTIMATES):
        color = style.CORAL if i == 0 else style.BLUE
        ax.plot([estimate - 1.96 * se, estimate + 1.96 * se], [i, i], color=color, lw=3)
        ax.plot(estimate, i, "o", ms=7, color=color)
    ax.set_yticks(range(len(ESTIMATES)), [e[0] for e in ESTIMATES])
    ax.invert_yaxis()
    ax.set(xlabel="slope: log investment per log unit of market value")


fig, ax = plt.subplots(figsize=(8, 2.6))
plot_estimates(ax)
plt.show()

Figure 4.7.3: The slope of log investment on log market value with its 95% interval, estimated four ways. Complete pooling with ordinary standard errors treats the 220 rows as independent and gives an interval about a quarter as wide as the others. Cluster-robust standard errors, firm fixed effects and the mixed model all account for the firms.

The naive interval is wrong, not just a little optimistic. Every test of Part 3 and every regression in this part assumes independent cases; with clusters, use one of the other three.

• Cluster-robust standard errors keep the simple model and only fix the uncertainty. They need a fair number of clusters (dozens, ideally).
• Fixed effects remove everything that differs between firms and does not change over time, which protects against confounders at the firm level, but they cannot estimate the effect of anything that is constant within a firm.
• Mixed models describe the clusters, predict for new ones, and extend to random slopes and nested clusters (pupils in classes in schools).

4.7.4 Random slopes

A random intercept lets each firm sit higher or lower; a random slope also lets the effect of market value differ by firm, again pulled towards the common slope:

$$y_{ij} = (\beta_0 + u_{0j}) + (\beta_1 + u_{1j})\, x_{ij} + \varepsilon_{ij}$$

slopes_model = smf.mixedlm(
    "log_invest ~ log_value", grunfeld, groups=grunfeld.firm, re_formula="~log_value"
).fit()


def plot_random_slopes(ax: Axes) -> None:
    """The data and one fitted line per firm."""
    for firm, part in grunfeld.groupby("firm"):
        effect = slopes_model.random_effects[firm]
        a = slopes_model.params.Intercept + effect.iloc[0]
        b = slopes_model.params.log_value + effect.iloc[1]
        grid = np.array([part.log_value.min(), part.log_value.max()])
        ax.scatter(part.log_value, part.log_invest, s=8, color=FIRM_COLORS[firm], alpha=0.5)
        ax.plot(grid, a + b * grid, color=FIRM_COLORS[firm], lw=2)
    grid = np.array([grunfeld.log_value.min(), grunfeld.log_value.max()])
    ax.plot(
        grid,
        slopes_model.params.Intercept + slopes_model.params.log_value * grid,
        color=style.INK,
        lw=2.5,
        ls="--",
    )
    ax.set(xlabel="log market value", ylabel="log investment")


fig, ax = plt.subplots(figsize=(8, 3.6))
plot_random_slopes(ax)
plt.show()

Figure 4.7.4: Lines for each firm from a mixed model with random intercepts and random slopes. The firms differ mostly in level; their slopes are similar and pulled close to the common slope (black).

Mixed models are also the standard tool for repeated measurements (several measurements of the same person over time) and for designs with several layers of clustering. Generalised versions handle yes/no and count outcomes, combining this chapter with Chapters 4.4 and 4.6.

4.7.5 In one picture

• Clustered cases are alike; the intraclass correlation measures how much.
• Ignoring clusters gives far too small standard errors.
• No pooling (fixed effects) respects clusters but shares nothing; complete pooling shares everything; partial pooling (mixed models) shrinks each cluster towards the average.
• Random intercept model: $y_{ij} = \beta_0 + u_j + \beta_1 x_{ij} + \varepsilon_{ij}$; add random slopes when effects differ by cluster.
• Alternatives: cluster-robust standard errors with many clusters, fixed effects to remove cluster-level confounding.

Next: the same units measured again and again over time.