GLMM Fundamentals

This page covers the statistical theory behind the GLMM tab. See that page for usage instructions.

Hierarchical Data and the Independence Problem

Most statistical models assume observations are independent. In practice, data often has a hierarchical structure where observations are grouped:

Students nested within schools
Patients nested within hospitals
Repeated measurements within subjects

Observations in the same group share group-level factors (school policies, hospital resources, individual physiology) and tend to be more similar to each other. Analyzing such data with GLM while ignoring group structure is a model misspecification: the data-generating process has hierarchical structure that the model fails to represent.

This misspecification harms the estimation results in two ways.

Underestimated standard errors: A GLM that ignores group structure computes standard errors under the assumption of zero covariance between observations. Because observations within the same group are actually positively correlated, standard errors are underestimated for the intercept and for coefficients of predictors that vary between groups, so 95% confidence intervals do not actually achieve 95% coverage. The degree of underestimation depends on whether a predictor's variation lies between or within groups (see the design effect in the ICC section).

Biased fixed effects: When group-level confounders exist, the fixed effect estimates suffer omitted variable bias. For example, if school policy influences both study hours and test scores, a regression ignoring school differences cannot correctly estimate the effect of study hours. Random intercepts model mean-level differences between groups, but their estimation assumes that the random intercepts are uncorrelated with the predictors. When unmeasured group-level confounders are correlated with predictors, this assumption fails and the fixed effect estimates remain biased. Measured confounders must be included as fixed effects.

Mixed models resolve the underestimation of standard errors by explicitly modeling within-group correlation as random effects. They also estimate between-group and within-group variability separately, quantifying how much of the variability stems from group differences; the ICC described below summarizes this decomposition.

Fixed Effects and Random Effects

The "mixed" in mixed model refers to containing both fixed and random effects.

Fixed effects are systematic predictor–response relationships assumed to be common to all groups, modeled with unknown constant coefficients. For example, "how many points does test score increase per additional hour of study" is a fixed effect.

Random effects are a modeling device for representing unobserved group-level variability within the model. Each school's average score differs due to policies, student demographics, and other factors. This variability is not generated by a random mechanism — it stems from concrete causes — but measuring and including all those factors as predictors is rarely practical. Random effects approximate this "unobserved group difference" using a probability distribution.

If we write the difference for group $j$ as $u_j$ , the assumption $u_j \sim N(0, \sigma_u^2)$ is not a claim that school differences are generated from a normal distribution. It is a modeling assumption that enables estimation of the overall magnitude of group variability $\sigma_u^2$ and produces shrinkage estimates (see below).

The practical difference between the two shows in how parameters are counted. Treating the groups as dummy variables with fixed effects adds one free parameter per group, and each group's effect becomes an estimation target in its own right. With random effects, the only variance component estimated is the single parameter $\sigma_u^2$ ; each group's $u_j$ is not a free parameter but a predicted value based on that variance. The fixed/random distinction is not about whether groups were "randomly sampled" but about which of these two forms the group differences take (Gelman, 2005).

Random Intercept Model

The simplest mixed model is the random intercept model, which gives each group its own baseline:

g(\mu_i) = x_i'\beta + u_{j[i]}, \quad u_j \sim N(0, \sigma_u^2)

where $g$ is the link function and $x_i'\beta$ is the fixed-effect linear predictor. The subscript $j[i]$ denotes the group that observation $i$ belongs to, and $u_{j[i]}$ is that group's random intercept. The statement $u_j \sim N(0, \sigma_u^2)$ means that the group-level values $u_1, \dots, u_J$ each independently follow the same normal distribution.

Intuitively, there is a single regression line shared by all groups ( $x_i'\beta$ ), and each group's line is shifted up or down by $u_j$ . The normal distribution $N(0, \sigma_u^2)$ assumption is a tool for summarizing the magnitude of these shifts with a single parameter $\sigma_u^2$ . The actual value of each $u_j$ is inferred from the data (see Predicting Random Effects and Shrinkage).

For the Gaussian family, this is a linear mixed model (LMM):

y_i = x_i'\beta + u_{j[i]} + \varepsilon_i, \quad u_j \sim N(0, \sigma_u^2), \quad \varepsilon_i \sim N(0, \sigma_e^2)

$\sigma_u^2$ is the between-group variance; $\sigma_e^2$ is the within-group (individual-level) variance. The residual variance $\sigma_e^2$ is assumed common to all groups: the random intercept model represents group differences in mean level and does not model structures where the magnitude of within-group variability differs by group.

In this equation, the similarity of observations within a group appears as covariance. Two observations in the same group share $u_j$ , so their covariance is $\sigma_u^2$ ; observations in different groups have covariance 0. Since each observation has variance $\sigma_u^2 + \sigma_e^2$ , the correlation between two observations in the same group is $\sigma_u^2 / (\sigma_u^2 + \sigma_e^2)$ — the quantity measured by the ICC described below.

When the link function is not the identity, the coefficients $\beta$ have a conditional interpretation: they express the effect of moving a predictor within the same group, on the link scale. Because $g^{-1}$ is nonlinear, the effect on the population-averaged response, with the random effects averaged out, generally differs from this whenever $\sigma_u^2 > 0$ . For example, an odds ratio obtained from a Binomial + logit model is read as a within-group comparison, not as a ratio of population-averaged odds.

Parameter Estimation

Estimation in mixed models is more complex than in standard regression because fixed effects $\beta$ and variance components ( $\sigma_u^2$ , $\sigma_e^2$ ) must be estimated simultaneously.

REML (Restricted Maximum Likelihood)

Maximum likelihood (ML) tends to underestimate variance components because the likelihood does not account for the degrees of freedom consumed by estimating $\beta$ . This is the same phenomenon that makes the sample variance with divisor $n$ an underestimate, which the usual $n-1$ divisor corrects.

To avoid this bias, REML estimates the variance components from quantities unaffected by $\beta$ . Viewing the response vector $y$ as a point in $n$ -dimensional space, the fixed effects $X\beta$ can only move within the $p$ -dimensional subspace spanned by the columns of the design matrix $X$ . Projecting $y$ onto the orthogonal complement of that subspace yields residuals whose $X\beta$ component vanishes exactly, so their distribution depends only on the variance components, not on $\beta$ . REML maximizes the likelihood of these residuals. Since no degrees of freedom are spent estimating $\beta$ , the variance estimation is effectively based on $n - p$ dimensions of data, which corrects the bias. The sample variance with divisor $n - 1$ is the special case where $X$ contains only an intercept ( $p = 1$ ).

MIDAS maximizes the REML likelihood by profiling, which reduces the search to one dimension. For each value of the variance ratio $\theta^2 = \sigma_u^2/\sigma_e^2$ , the optimal $\beta$ and $\sigma_e^2$ under that ratio have closed forms, so the only quantity searched numerically is $\log\theta$ , maximized by golden-section search.

Laplace Approximation and PIRLS

For the Gaussian family with the identity link, the random-effect integral is analytically tractable, so REML suffices. For all other combinations (Binomial, Poisson, Gamma, and Gaussian + log), integrating out the random effects:

L(\beta, \sigma_u^2) = \int \prod_i f(y_i \mid \beta, u) \cdot f_u(u) \, du

has no closed-form solution. Here $f$ is the density of each observation under the chosen distribution family (a probability mass function for discrete families such as Binomial and Poisson), and $f_u$ is the density of the random effects, i.e. the $N(0, \sigma_u^2)$ density.

The Laplace approximation replaces this integral with a second-order expansion around the mode of the integrand. The approximation may lose accuracy when group sizes are small or when events are rare in Binomial data.

MIDAS's estimation algorithm follows the nested two-loop structure described in Bates et al. (2015); the choice of optimizer for the outer loop is specific to MIDAS.

Outer loop: Optimizes the log of the relative covariance parameter $\theta$ via golden-section search
Inner loop: PIRLS (Penalized IRLS) simultaneously estimates $(\beta, u)$

PIRLS extends GLM's IRLS with a random-effect penalty and minimizes the following penalized deviance:

D_p = \sum_i d(y_i, \mu_i) + \sum_j \frac{u_j^2}{\theta^2}

$\theta$ is the relative covariance parameter optimized by the outer loop; it expresses the standard deviation of the random effects relative to the residual scale: $\theta = \sigma_u/\sigma_e$ for Gaussian, $\theta = \sigma_u$ for Poisson/Binomial, and $\theta = \sigma_u/\sqrt\phi$ for Gamma, where $\phi$ is the dispersion parameter. The first term measures fit to data; the second penalizes $u_j$ values that deviate too far from the normal prior. This penalty pulls estimates for data-sparse groups toward the overall mean (see shrinkage below).

Boundary Estimates (Singular Fit)

The estimate of the variance $\sigma_u^2$ has a lower bound of 0. In either estimation path, a fit whose optimum lands on this bound is called a singular fit. It occurs when the between-group variation in the data is negligible relative to the residual variation.

A boundary estimate changes how the results should be read. $\hat\sigma_u^2 = 0$ is not the conclusion that no group differences exist; it means the data could not distinguish between-group variation from residual variation. Quantities that depend on $\hat\sigma_u^2$ — the ICC and the predicted random effects — are also nearly 0. Moreover, asymptotic approximations for the variance components assume that the true parameter lies in the interior of the parameter space, and this assumption fails at the boundary. See the GLMM tab page for how this state is detected and handled.

AIC/BIC Limitations

Since the REML projection depends on the fixed-effect structure, REML-based AIC/BIC cannot be compared across models with different fixed effects. Gaussian + identity is the only combination MIDAS estimates via REML, so it is the only one subject to this constraint. All other combinations derive AIC/BIC from a maximum-likelihood quantity (the Laplace-approximated marginal log-likelihood), so models with different fixed effects can be compared.

In either path, however, comparisons are valid only within the same family and link. Changing the family or link mixes factors unrelated to model fit into the AIC/BIC difference: a different estimation path changes the basis of the log-likelihood (REML vs. maximum likelihood), a different family changes what the log-likelihood measures (probability mass for discrete families, density for continuous ones) and how the dispersion parameter is handled, and a different link changes how the Laplace approximation error behaves.

Predicting Random Effects and Shrinkage

The fitted model yields not only the fixed-effect coefficients but also a value of $u_j$ for each group. These values are not computed from each group's data alone. Because $u_j$ is assigned the probability distribution $N(0, \sigma_u^2)$ , its value is inferred by combining two sources of information: the data and the assumed distribution. This combination produces the shrinkage described in this section.

This difference from the fixed effects also shows in terminology. The fixed effects $\beta$ are unknown constants, so they are "estimated"; the $u_j$ are treated as random variables within the model, so they are "predicted." This is a formal distinction reflecting the fact that for $\beta$ only the finiteness of the data contributes uncertainty, while for $u_j$ the assumed distribution also acts as an information source — not a claim that $u_j$ was truly generated randomly.

For Gaussian models (LMM), random effects are predicted with the BLUP (Best Linear Unbiased Predictor). The BLUP is the linear predictor that minimizes the mean squared prediction error $E[(\hat u_j - u_j)^2]$ subject to the constraint that unbiasedness with respect to $u_j$ — $E[\hat u_j] = E[u_j] = 0$ — holds for every value of the unknown fixed effects $\beta$ . It derives from the mixed model equations of Henderson (1975) and takes the form of a shrinkage estimator:

\hat{u}_j = \frac{n_j}{n_j + \sigma_e^2 / \sigma_u^2} \times \bar{r}_j

where $n_j$ is the group size and $\bar{r}_j$ is the mean residual for group $j$ (the part not explained by fixed effects).

The coefficient $n_j / (n_j + \sigma_e^2/\sigma_u^2)$ ranges from 0 to 1, approaching 1 as group size increases:

Large groups: Enough data to trust the group-specific estimate, so use it nearly as-is
Small groups: Limited data, so pull the estimate toward the overall mean (zero)

This is also called "borrowing strength." The phrase comes from Tukey, who used it for this idea of shrinkage in his election-night forecasting work in the 1960s (Brillinger, 2002). Data-sparse groups borrow information from other groups to stabilize their estimates. Using group-specific estimates directly would have high variance; using only the overall mean would ignore group characteristics. BLUP optimally balances this bias-variance tradeoff.

Each predicted value comes with an assessment of its uncertainty. Writing the shrinkage coefficient as $c_j = n_j / (n_j + \sigma_e^2/\sigma_u^2)$ , the standard error of $\hat u_j$ is the square root of the variance of $u_j$ conditional on the observed data, $\sigma_u^2 (1 - c_j)$ , with the estimated variance components treated as the true values. As a group's data becomes sparse, $c_j$ approaches 0, so the prediction shrinks to 0 and the standard error approaches the distribution's standard deviation $\sigma_u$ ; as data accumulates, the standard error tends to 0.

For non-Gaussian families, random effects are estimated as conditional modes instead: the value of $u_j$ that maximizes the density of its distribution conditional on the observed responses. If the assumed distribution $N(0, \sigma_u^2)$ is viewed as a prior, this is the mode of the posterior distribution. Conditional modes are obtained as the PIRLS solution and exhibit shrinkage toward the overall mean just like the Gaussian BLUP, but they are not linear predictors, so strictly speaking they are not BLUPs. Their uncertainty is assessed from a normal approximation around the conditional mode.

ICC (Intraclass Correlation Coefficient)

ICC measures the share of unexplained variance — the variance remaining after the fixed effects are accounted for — attributable to between-group differences:

\text{ICC} = \frac{\sigma_u^2}{\sigma_u^2 + \sigma_e^2}

$\sigma_u^2$ and $\sigma_e^2$ are the variance components of the fitted model. In a model with fixed effects $X\beta$ , both are residual-side quantities after removing the variance explained by the fixed effects.

The ICC is also the correlation between two observations in the same group. As derived in the random intercept model section, that correlation is $\sigma_u^2 / (\sigma_u^2 + \sigma_e^2)$ — the same quantity as this definition. This property is what the name intraclass correlation coefficient refers to.

This form of the ICC is called the conditional ICC, and it is what MIDAS GLMM computes. Its counterpart, the unconditional ICC, is computed from an intercept-only model with no predictors. Because the conditional ICC depends on the fixed-effect specification, adding or removing predictors changes its value, so it generally does not equal the unconditional ICC.

Non-Gaussian Families

For non-Gaussian families, $\sigma_e^2$ is not directly estimated from data. For the ICC to be a variance decomposition, $\sigma_u^2$ and $\sigma_e^2$ must be components that split the variability of the same quantity on the same scale. A $\sigma_e^2$ satisfying this condition exists only when it can be derived theoretically from the link function and the distributional assumptions.

Binomial models with logit or probit links have a threshold model (Goldstein et al. 2002): an unobserved continuous variable $y^* = X\beta + u + e$ generates the binary response as $y = \mathbf{1}(y^* > 0)$ . The distribution of $e$ — and hence its variance — follows deductively from the choice of link function. With the logit link, $e$ follows a logistic distribution with variance $\pi^2/3$ . With the probit link, $e$ follows $N(0, 1)$ with variance $1$ . In this framework, $\sigma_u^2 + \sigma_e^2$ represents the total variance of the latent variable, and ICC is a genuine variance decomposition.

Family + Link	Residual Variance	Basis
Binomial + Logit	$\pi^2/3 \approx 3.29$	Logistic distribution (threshold model)
Binomial + Probit	$1$	Standard normal distribution (threshold model)

The Binomial ICC is a value on the scale of the latent variable $y^*$ . The within-group correlation of the observed binary responses is a different quantity: it is smaller than the latent-scale value because dichotomization discards information, and it also depends on the overall event probability.

ICC is a valid variance decomposition only for Gaussian + identity (where $\sigma_e^2$ is estimated directly via REML) and the two Binomial combinations above. For other family+link combinations, ICC is not defined:

Poisson (all links): Poisson distributions have $\operatorname{Var}(Y \mid \mu) = \mu$ , so the variance is fully determined by the mean and there is no independent residual variance parameter. Some conventions place $1$ in the ICC denominator, but $\sigma_u^2 / (\sigma_u^2 + 1)$ is a monotone transformation of $\sigma_u^2$ , not a variance decomposition.
Gamma (all links): The profiled dispersion $\phi$ is a conditional variance parameter, not a latent-scale residual variance.
Gaussian + log: $\sigma_u^2$ is on the link scale while $\sigma_e^2$ is estimated on the response scale. There is no theoretical residual variance on the link scale to reconcile the two.

The Design Effect and Choosing Between GLM and Mixed Models

How much a GLM that ignores group structure understates the standard errors can be gauged from the ICC and the group sizes. Observations in the same group resemble each other and carry less information than independent observations, so the variance of a coefficient estimator is larger than it would be with independent data of the same size. This ratio is the design effect (DEFF), expressed as $\text{DEFF} = 1 + (\bar n - 1) \times \text{ICC}$ , where $\bar n$ is the average group size. The formula is exact when group sizes are equal and approximate when they are unequal. Because standard errors computed under the independence assumption do not include this inflation, the true standard error is roughly $\sqrt{\text{DEFF}}$ times the reported one.

The inflation differs by predictor. For group means and coefficients of predictors that vary at the group level this guide applies as is; for coefficients of predictors that vary within groups the inflation is smaller, and for predictors with little between-group variation it is nearly absent.

This estimate informs the choice between GLM and a mixed model. When the ICC is small and groups are not large, DEFF stays close to 1, and ignoring the group structure with a GLM yields nearly identical results. When DEFF is well above 1, a mixed model is necessary.

Fixed Effect Inference and the Normal Approximation

Fixed-effect confidence intervals are Wald-based, constructed as $\hat\beta_j \pm z_{1-\alpha/2} \cdot \text{SE}(\hat\beta_j)$ . $\text{SE}(\hat\beta_j)$ is the square root of the $j$ -th diagonal entry of $(X'\hat V^{-1}X)^{-1}$ , with the estimated variance components plugged into $\hat V$ . For Gaussian + identity, $\hat V = \hat\sigma_u^2 ZZ' + \hat\sigma_e^2 I$ is the covariance matrix of the response, where $Z$ is the indicator matrix of group membership. For all other combinations, $\hat V = \hat W^{-1} + \hat\theta^2 ZZ'$ is the approximate covariance of the working response built from the PIRLS weights at convergence, further multiplied by $\hat\phi$ for families that estimate a dispersion parameter. This construction relies on $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ approximately following the standard normal distribution.

In ordinary regression (OLS), the error variance $\sigma^2$ is estimated, so $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ follows a t-distribution with $n - p$ degrees of freedom. In mixed models, there are multiple variance components ( $\sigma_u^2$ , $\sigma_e^2$ ), and the corresponding degrees of freedom are not uniquely determined. The normal approximation treats the estimated variance components as known.

In LMM, the normality assumptions on $u_j$ and $\varepsilon_i$ ensure that $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ follows the standard normal exactly when the variance components are known. In non-Gaussian GLMM, the normal distribution is an asymptotic approximation. In both cases, replacing the true variance components with estimates introduces additional uncertainty, making the distribution of $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ heavier-tailed than the standard normal.

When the number of groups is large, the variance component estimates are stable and this approximation is adequate. With few groups, confidence intervals tend to be too narrow.

Random Slopes and Crossed Random Effects

MIDAS GLMM supports random intercept models only. Since a random intercept captures group differences in baseline level and nothing more, this section describes the two extensions beyond that scope, as material for judging whether a random intercept model suffices for your data.

When the effect of a predictor also varies by group, a random slope model is needed:

g(\mu_i) = x_i'\beta + u_{0j[i]} + u_{1j[i]} x_{1i}

where $u_{0j}$ is the random intercept and $u_{1j}$ is the random slope, jointly following a multivariate normal distribution.

If the slopes do vary between groups and the data is nevertheless fit with a random intercept model, that variation enters the error term. The within-group covariance then actually depends on the predictor values, but the random intercept model replaces it with the constant $\sigma_u^2$ . This leaves the covariance structure misspecified, so the standard error of that predictor's coefficient becomes unreliable.

When observations belong to multiple grouping variables that are not nested in each other, crossed random effects are used. For example, if students belong to both a school and a region, a random effect for each grouping variable enters the linear predictor:

g(\mu_i) = x_i'\beta + u_{j[i]} + v_{k[i]}, \quad u_j \sim N(0, \sigma_u^2), \quad v_k \sim N(0, \sigma_v^2)

where $j[i]$ and $k[i]$ index the school and the region that observation $i$ belongs to, and the two random effects are assumed independent of each other.

A random intercept model has only one grouping variable, so with crossed data the variation from the other grouping variable remains in the residuals. That remaining variation produces correlation within the groups of that variable, and ignoring it understates standard errors in the same way described in Hierarchical Data and the Independence Problem.

References

Gelman, A. (2005). Analysis of variance—why it is more important than ever. The Annals of Statistics, 33(1), 1-53. https://www.jstor.org/stable/3448650
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1-48. https://www.jstatsoft.org/v67/i01/
Henderson, C. R. (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics, 31(2), 423-447. https://www.jstor.org/stable/2529430
Brillinger, D. R. (2002). John W. Tukey: His life and professional contributions. The Annals of Statistics, 30(6), 1535-1575. https://doi.org/10.1214/aos/1043351246
Goldstein, H., Browne, W., & Rasbash, J. (2002). Partitioning variation in multilevel models. Understanding Statistics, 1(4), 223-231. https://doi.org/10.1207/S15328031US0104_02

Also available as a Markdown file.