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$ . The REML log-likelihood can have more than one local maximum in $\log\theta$ , so MIDAS evaluates it at evenly spaced points over the search range and runs golden-section search on the interval between the two neighbors of the point with the highest value. The reason for this procedure, and the case in which it still misses the maximum, are explained for the outer loop in Laplace Approximation and PIRLS.

Laplace Approximation and PIRLS

For Gaussian + identity, the integral over the random effects can be computed analytically, so MIDAS estimates that combination by REML as described in the previous section. For all other combinations (Binomial, Poisson, Gamma, and Gaussian + log), the marginal likelihood obtained by integrating over the random effects

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

has no closed form. 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. $\mu_i$ is the mean of observation $i$ , determined by $\beta$ and $u$ , and $\phi$ is the dispersion parameter of the family. MIDAS estimates $\phi$ for Gaussian + log and Gamma, and fixes $\phi = 1$ for Poisson and Binomial.

The Laplace approximation replaces the log of the integrand with a second-order expansion around its maximum, which gives the integrand the shape of a normal density, and computes this integral from that shape. The approximation loses accuracy when group sizes are small and when events are rare in Binomial data.

MIDAS writes the random-effect variance as $\sigma_u^2 = \theta^2\phi$ using the relative covariance parameter $\theta$ , and computes the Laplace-approximated marginal likelihood with the following expression. For Gaussian, $\phi = \sigma_e^2$ , so $\theta^2 = \sigma_u^2/\sigma_e^2$ . For Gaussian + log, however, $\sigma_u$ is on the scale of the linear predictor (the log scale), so unlike the $\theta$ of the REML section this $\theta$ depends on the units of the response.

-2 \log L \approx \sum_j \log(1 + \theta^2 C_j) + \frac{1}{\phi}\left(\sum_i d(y_i, \hat\mu_i) + \sum_j \frac{\hat u_j^2}{\theta^2}\right) + \sum_i a(y_i, \phi)

The terms $d$ and $a$ in the expression come from splitting the density of each observation as $-2 \log f(y \mid \mu, \phi) = d(y, \mu)/\phi + a(y, \phi)$ . $d(y, \mu)$ is the observation's contribution to the deviance, and $a(y, \phi)$ is a term that does not involve $\mu$ . With this split, the log of the integrand for group $j$ takes the form $-\{\sum_{i \in j} d(y_i, \mu_i) + u_j^2/\theta^2\}/(2\phi)$ plus terms that do not involve $u_j$ . The value $\hat u_j$ that maximizes the integrand is therefore determined by $\theta$ and $\beta$ alone. $\hat\mu_i$ is the mean of observation $i$ with $u_j = \hat u_j$ . $C_j$ is the sum, over the observations in group $j$ , of half the second derivative of $d$ with respect to the linear predictor, evaluated at $u_j = \hat u_j$ . The first term of the expression is the sum over groups of the log of the ratio between the variance $\sigma_u^2$ of the random-effect distribution and the variance $\phi/(C_j + 1/\theta^2)$ of the normal density that replaces the integrand.

MIDAS maximizes this approximate marginal likelihood over all of the fixed effects $\beta$ , $\theta$ , and $\phi$ . $\hat\phi$ is the value that minimizes the expression above over $\phi$ for given $\theta$ and $\beta$ . For Gaussian + log it is $\hat\phi = \{\sum_i d(y_i, \hat\mu_i) + \sum_j \hat u_j^2/\theta^2\}/n$ . For Gamma, MIDAS solves the corresponding equation numerically to obtain $\hat\phi$ . The reported log-likelihood, estimate of $\phi$ , and $\sigma_u^2 = \theta^2\phi$ all use this $\hat\phi$ . For Gaussian + log, $\hat\phi$ is the estimate of the residual variance $\sigma_e^2$ ; for Gamma, $\hat\phi$ is the coefficient of the conditional variance of the response $\phi\mu^2$ and is not a variance itself. Because the maximization is iterative, a converged fit gives a local maximum, which is not necessarily the global maximum. When the iterations do not converge, MIDAS either attaches a warning to the estimates or aborts the fit because the estimate cannot be determined. See Convergence Issues on the GLMM tab for which happens when.

The maximization follows the nested two-loop structure described in Bates et al. (2015). The optimizer of the outer loop and the procedure that finds the fixed effects in the inner loop are specific to MIDAS.

Outer loop: Maximizes the approximate marginal likelihood over the log of the relative covariance parameter $\theta$ , first by evaluating it at points over the search range and then by golden-section search
Inner loop: With $\theta$ fixed, starts from the PIRLS (Penalized IRLS) solution and maximizes the approximate marginal likelihood over $\beta$

The outer loop evaluates the approximate marginal likelihood at several points before golden-section search because this likelihood, maximized over the parameters other than $\theta$ , can have more than one local maximum in $\theta$ . For some data both the boundary $\theta = 0$ and a point with $\theta > 0$ are local maxima, and the point with $\theta > 0$ has the higher likelihood. Golden-section search assumes a single maximum and discards one side of the search range by comparing the values at two points. When applied to the whole search range, it can therefore discard the side containing the maximum with $\theta > 0$ and stop at a boundary estimate. MIDAS evaluates the approximate marginal likelihood at evenly spaced points over the search range of $\log\theta$ and runs golden-section search on the interval between the two neighbors of the point with the highest value. Golden-section search narrows the interval starting from the point with the highest value, so the resulting likelihood is never lower than the value at that point. If the local maximum with the highest likelihood lies outside this interval, MIDAS misses it. When the interval contains more than one local maximum, the result is one of them.

PIRLS extends GLM's IRLS with a random-effect penalty and minimizes the following penalized deviance over $(\beta, u)$ :

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

The first term of $D_p$ measures fit to data, and the second is a penalty that grows as $u_j$ moves away from 0. The second term comes from the log of the density of the distribution $N(0, \sigma_u^2)$ of $u_j$ , and it is heavier when $\theta$ is smaller. This penalty pulls the predictions for data-sparse groups toward the overall mean (see shrinkage below).

The $\beta$ of the PIRLS solution generally does not coincide with the maximum of the approximate marginal likelihood. PIRLS minimizes $D_p$ alone, whereas the term $\sum_j \log(1 + \theta^2 C_j)$ of the approximate marginal likelihood also depends on $\beta$ through $C_j$ . MIDAS therefore starts from the PIRLS solution and maximizes the approximate marginal likelihood, including this term, over $\beta$ . Each time $\beta$ moves, it recomputes $\hat u_j$ as the value that minimizes $D_p$ over $u_j$ .

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. A boundary estimate occurs when the between-group variation in the data is small relative to the residual variation. It also occurs when the search over $\log\theta$ stops at a local maximum on the boundary $\theta = 0$ and misses a maximum with a higher likelihood away from the boundary. Laplace Approximation and PIRLS explains when the search misses a maximum.

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 is the state in which the estimate lies on the boundary. 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

Models with different fixed effects can be compared by AIC/BIC only for the combinations estimated with the Laplace approximation. 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 the Laplace-approximated marginal log-likelihood maximized over all parameters, including the fixed effects. This value approximates the marginal likelihood of the response $y$ itself rather than of projected residuals, and the definition of that likelihood does not change with the fixed-effect structure, so models with different fixed effects can be compared. The error of the Laplace approximation also changes with the fixed-effect structure, so a difference in AIC/BIC includes the difference in the approximation errors. This comparison is also limited to models fitted to the same observations. MIDAS drops rows in which any selected column is missing, so adding or removing a predictor with missing values changes the observations used, and the log-likelihoods can no longer 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 combinations other than Gaussian + identity, MIDAS uses conditional modes as the predictions of the random effects: 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. The conditional mode is the value of $u_j$ that minimizes $D_p$ at the estimated $\beta$ and $\theta$ . Conditional modes shrink toward the overall mean just like the BLUP, but they are not linear predictors, so they are not BLUPs.

The standard error of a conditional mode is the square root of $\phi/(C_j + 1/\theta^2)$ , the variance of the normal approximation around the conditional mode. $C_j$ is the same quantity as the $C_j$ in the expression of the Laplace approximation. The estimated $\beta$ , $\theta$ , and $\phi$ are treated as known.

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

The construction of fixed-effect confidence intervals differs between LMM (Gaussian + identity) and all other GLMMs. For LMM, the uncertainty from estimating the variance components is accounted for by the Kenward-Roger method; for other GLMMs, the standard normal approximation is used.

LMM: the Kenward-Roger Method

LMM fixed-effect confidence intervals are constructed as $\hat\beta_j \pm t_{1-\alpha/2,\,\nu_j} \cdot \text{SE}(\hat\beta_j)$ . A naive standard error can be read off the diagonal of $(X'\hat V^{-1}X)^{-1}$ , where $\hat V = \hat\sigma_u^2 ZZ' + \hat\sigma_e^2 I$ and $Z$ is the indicator matrix of group membership, but this quantity omits the extra variability introduced by plugging estimates of the variance components $(\sigma_u^2, \sigma_e^2)$ into $\hat V$ , and therefore understates the standard error. The Kenward-Roger method takes the standard error from an adjusted covariance matrix that corrects this understatement, and derives a per-coefficient degrees of freedom $\nu_j$ from the estimation precision of the variance components (Kenward & Roger, 1997).

The degrees of freedom $\nu_j$ expresses how much information the coefficient is estimated from. Coefficients of predictors that vary only between groups have small degrees of freedom, because their effective sample size is the number of groups. In a balanced design with only group-level predictors, inference for those coefficients coincides with an ordinary regression on the group means, and the degrees of freedom equals its residual degrees of freedom (the number of groups minus the number of estimated group-level coefficients). Coefficients of predictors that vary within groups have large degrees of freedom, and their t-distribution approaches the standard normal. For coefficients with small degrees of freedom, the t quantile is larger and the confidence interval widens accordingly.

With few groups, the corrected intervals are distinctly wider than normal-approximation intervals. For example, the t quantile at 3 degrees of freedom, $t_{0.975,3} \approx 3.18$ , is about 1.6 times $z_{0.975} \approx 1.96$ . This width faithfully reflects how little information a handful of groups carries about the variance components; the narrower normal-approximation interval falls short of its nominal coverage.

This construction depends on the normality of the response, the normality of the random intercepts, independence across groups, and the correctness of the covariance structure $V = \sigma_u^2 ZZ' + \sigma_e^2 I$ . None of these assumptions can be relaxed. When they fail, for example when data that require a random slope are analyzed with a random intercept model, the standard errors themselves are unreliable with or without the correction (see Random Slopes and Crossed Random Effects).

When the variance component structure is degenerate and the adjustment cannot be computed, MIDAS shows a warning and constructs the intervals from the unadjusted standard errors and the standard normal approximation.

Non-Gaussian GLMM: the Standard Normal Approximation

For combinations other than Gaussian + identity, confidence intervals are 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 corresponding diagonal element of $(X'\hat V^{-1}X)^{-1}$ . $\hat V = \hat\phi(\hat W^{-1} + \hat\theta^2 ZZ')$ is the approximate covariance of the adjusted dependent variable of GLM's IRLS, with $\hat\phi = 1$ for Poisson and Binomial. $\hat W$ is the IRLS weight matrix evaluated at the means $\hat\mu$ computed from the fixed-effect estimates and the conditional modes. $\hat\phi$ is defined differently from the dispersion parameter of GLM (deviance/(n−p) for Gaussian, Pearson $\chi^2$ /(n−p) for Gamma), so the standard errors of GLMM and GLM do not coincide even when the between-group variance is nearly 0. This construction relies on $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ asymptotically following the standard normal distribution, and treats the estimated variance components as known.

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. $\hat\phi$ is also a maximum-likelihood estimate without a degrees-of-freedom correction, so with few observations it understates $\phi$ and makes the confidence intervals narrower.

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
Kenward, M. G., & Roger, J. H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53(3), 983-997. https://www.jstor.org/stable/2533558
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.