MixedModelsSmallSample

Small-sample corrections for linear mixed models fitted by MixedModels.jl.

Linear Mixed Model

A linear mixed model expresses the $n$-vector of responses as

\[y = X\beta + Zu + \varepsilon, \quad u \sim \mathcal{N}(0, G(\theta)), \quad \varepsilon \sim \mathcal{N}(0, \sigma^2 I_n),\]

where $X$ is the $n \times p$ fixed-effects design matrix, $\beta$ is the $p$-vector of fixed-effects coefficients, $Z$ is the random-effects design matrix, $u$ is the vector of random effects with covariance $G(\theta)$, and $\varepsilon$ is the residual error.

The marginal variance of $y$ is

\[V(\theta) = ZG(\theta)Z^\top + \sigma^2 I_n.\]

The Small-Sample Problem

The generalized least squares estimator $\hat{\beta}$ has covariance

\[\Phi(\theta) = (X^\top V(\theta)^{-1} X)^{-1}.\]

Standard inference replaces $\theta$ with its REML estimate $\hat{\theta}$, yielding $\hat{\Phi} = \Phi(\hat{\theta})$, and uses asymptotic approximations (e.g., treating $t$-statistics as normally distributed). For small samples, this ignores:

1. Bias in $\hat{\Phi}$ due to estimating $\theta$
2. Uncertainty in $\hat{\theta}$ affecting the reference distribution

This package implements two corrections:

• Kenward-Roger (KenwardRoger): Applies a bias correction to $\hat{\Phi}$ and uses a scaled F-distribution with approximate degrees of freedom.
• Satterthwaite (Satterthwaite): Uses moment-matching to approximate degrees of freedom for the t/F reference distribution.

References

• Kenward, M. G. and Roger, J. H. (1997). "Small sample inference for fixed effects from restricted maximum likelihood." Biometrics, 53(3), 983-997.
• Satterthwaite, F. E. (1946). "An Approximate Distribution of Estimates of Variance Components." Biometrics Bulletin, 2(6), 110-114.
• Fai, A. H.-T. and Cornelius, P. L. (1996). "Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments." Journal of Statistical Computation and Simulation, 54(4), 363-378.

See Also

• small_sample_adjust: Apply KR or SW adjustment to a fitted model
• small_sample_ftest: Perform an adjusted F-test for contrasts

small_sample_adjust(m::MixedModel, method) -> LinearMixedModelAdjusted

Apply small-sample corrections to a fitted linear mixed model.

This is the main entry point for computing adjusted inference. The function:

1. Validates that m is unweighted; KenwardRoger additionally requires a REML fit
2. Computes $W = I(\theta)^{-1}$, the asymptotic covariance of $\hat{\theta}$ (see vcov_varpar)
3. For KenwardRoger: computes adjusted covariance $\Phi_A$

(see compute_adjusted_covariance)

1. Computes denominator degrees of freedom $\nu_k$ for each fixed effect (see compute_dof)

Arguments

• m::MixedModel: A fitted LinearMixedModel from MixedModels.jl (must use REML for KR; ML or REML for SW)
• method: Either KenwardRoger or Satterthwaite

Returns

• LinearMixedModelAdjusted

The returned object can be displayed to show adjusted coefficient tables with corrected standard errors (KR only) and degrees of freedom.

Example

using MixedModels, MixedModelsSmallSample
m = fit(MixedModel, @formula(y ~ x + (1|g)), data; REML=true)
kr = small_sample_adjust(m, KenwardRoger())
sw = small_sample_adjust(m, Satterthwaite())

See Also

• small_sample_ftest: Perform F-tests on contrasts

small_sample_ftest(m::LinearMixedModelAdjusted, L) -> (ν, F*)

Compute an adjusted F-test for the null hypothesis $L^\top\beta = 0$.

Returns $(ν, F^*)$ where $(ν, λ)$ is obtained from compute_dof and $F^* = λF$ with $F = β^\top Mβ/q$ and $M = L(L^\top\Phi L)^{-1}L^\top$ (using the unadjusted covariance $\Phi$).

Arguments

• m: An adjusted model (LinearMixedModelAdjusted)
• L: Contrast matrix ($p \times q$) where $p$ = number of fixed effects

Returns

• ν: Denominator degrees of freedom
• F*: Scaled F-statistic

Under $H_0$, $F^* \sim F(q, \nu)$ approximately.

Example

kr = small_sample_adjust(model, KenwardRoger())
L = [0 1 0; 0 0 1]'  # Test β₂ = β₃ = 0
ν, Fstar = small_sample_ftest(kr, L)
pvalue = 1 - cdf(FDist(size(L,2), ν), Fstar)

small_sample_ftest(m::MixedModel, L; method=KenwardRoger()) -> (ν, F*)

Convenience wrapper that first applies small_sample_adjust then computes the F-test.

Equivalent to:

m_adj = small_sample_adjust(m, method)
small_sample_ftest(m_adj, L)

Arguments

• m: A fitted MixedModel (must use REML for KR; ML or REML for SW)
• L: Contrast matrix
• method: KenwardRoger or Satterthwaite

vcov_varpar(m::MixedModel; fim=ObservedFIM(), parameterization=Unstructured())

Compute W = I(θ)^{-1}, the asymptotic covariance of the fitted variance parameters.

The Fisher information I(θ) is computed using either ObservedFIM or ExpectedFIM.

KenwardRoger(; fim=ObservedFIM(), parameterization=Unstructured())

Kenward-Roger method for small_sample_adjust and small_sample_ftest.

See also Satterthwaite for an alternative method.

Satterthwaite(; fim=ObservedFIM(), parameterization=Unstructured())

Satterthwaite method for small_sample_adjust and small_sample_ftest.

See also KenwardRoger for an alternative method.

ObservedFIM

Use the observed (Hessian-based) Fisher Information Matrix for the variance parameters.

Let $V$, $V^{-1}$, $dV_i$ and $d^2V_{ij}$ come from VarianceDecomposition. Define

\[\Phi = (X^\top V^{-1} X)^{-1},\qquad P = V^{-1} - V^{-1}X\,\Phi\,X^\top V^{-1},\qquad r = y - X\hat{\beta}.\]

Then the observed information is

\[I^{\text{obs}}_{ij}(\theta) = -\frac{1}{2}\operatorname{tr}\left(P dV_i P dV_j\right) + r^{\top} V^{-1} dV_i P dV_j V^{-1} r + \frac{1}{2}\operatorname{tr}\left(P d^2V_{ij}\right) - \frac{1}{2} r^{\top} V^{-1} d^2V_{ij} V^{-1} r.\]

See Also

• ExpectedFIM: Alternative using expected values
• vcov_varpar: Computes the inverse of the FIM

ExpectedFIM

Use the expected Fisher Information Matrix for the variance parameters.

Let $V^{-1}$, $dV_i$, $P_i$ and $Q_{ij}$ come from VarianceDecomposition, and define $\Phi = (X^\top V^{-1} X)^{-1}$.

Then the expected information used by vcov_varpar is

\[I^{\text{exp}}_{ij}(\theta) = \frac{1}{2}\operatorname{tr}\left(V^{-1} dV_i V^{-1} dV_j\right) - \operatorname{tr}(\Phi\,Q_{ij}) + \frac{1}{2}\operatorname{tr}(\Phi\,P_i\,\Phi\,P_j).\]

See Also

• ObservedFIM: Alternative using observed Hessian
• vcov_varpar: Computes the inverse of the FIM

Unstructured

Unstructured (variance component) parameterization.

In this parameterization, the marginal variance is linear in the variance parameters. Writing the model as $V(\theta) = \sigma^2 I_n + \sum_b Z_b G_b(\theta) Z_b^\top$, Unstructured parameterizes each block covariance $G_b(\theta)$ by its unique (co)variance components.

Gradient Structure

Because $V$ is linear in $\theta$, the derivatives are trivial: each $\partial V/\partial\theta_i$ is constant and $\partial^2 V/(\partial\theta_i\,\partial\theta_j) = 0$.

This implies the $R_{ij}$ terms used by the Kenward-Roger covariance correction vanish.

See Also

• FactorAnalytic: Alternative Cholesky-based parameterization

FactorAnalytic

Factor analytic (Cholesky) parameterization.

For each random-effects block $b$, decompose the block covariance as $G_b = L_b L_b^\top$ with $L_b$ lower triangular. The variance parameters are $\sigma^2$ and the free elements of the $L_b$.

The marginal variance is

\[V = \sigma^2 I_n + \sum_b Z_b G_b Z_b^\top = \sigma^2 I_n + \sum_b Z_b L_b L_b^\top Z_b^\top.\]

For MixedModels.jl models, the internal random-effects factor is stored as a relative factor $\lambda_b$; we use the scaled factor $L_b = \sigma\,\lambda_b$.

Derivatives of V

The derivative of $V$ with respect to a Cholesky element $L_{b,ij}$ uses the chain rule through $G_b = L_b L_b^\top$:

\[\frac{\partial V}{\partial L_{ij}} = Z_b \frac{\partial G_b}{\partial L_{ij}} Z_b^\top, \qquad \frac{\partial G_b}{\partial L_{ij}} = E_{ij} L_b^\top + L_b E_{ij}^\top,\]

where $E_{ij}$ is the elementary matrix with a 1 at position $(i,j)$.

The second derivatives are non-zero. For indices within the same block,

\[\frac{\partial^2 G_b}{\partial L_{i_1 j_1}\,\partial L_{i_2 j_2}} = \delta_{j_1 j_2}\,(E_{i_1 i_2} + E_{i_2 i_1}),\]

and hence

\[\frac{\partial^2 V}{\partial L_{i_1 j_1}\,\partial L_{i_2 j_2}} = Z_b\,\frac{\partial^2 G_b}{\partial L_{i_1 j_1}\,\partial L_{i_2 j_2}}\,Z_b^\top.\]

Unlike Unstructured, this parameterization has non-zero second derivatives $\partial^2 V / \partial L_{ij} \partial L_{kl}$, which must be included in the Kenward-Roger bias correction.

See Also

• Unstructured: Alternative with zero second derivatives

LinearMixedModelAdjusted{T<:LinearMixedModel, O<:AbstractLinearMixedModelAdjustment}

A wrapper around a LinearMixedModel that holds the small-sample adjusted results.

Fields

• m: The original LinearMixedModel.
• adj: The adjustment option used (KenwardRoger or Satterthwaite).
• varcovar_adjusted: The adjusted variance-covariance matrix of the fixed effects (corrected for KR, unadjusted for SW).
• W: The asymptotic covariance matrix of the variance parameters.
• P: Derived matrix (P_i).
• Q: Derived matrix (Q_{ij}) (optional, nothing for Satterthwaite).
• R: Derived matrix (R_{ij}) (optional, nothing for Satterthwaite).
• R_factor: Scaling factor for (R_{ij}) (0.0 for Satterthwaite).
• ν: The denominator degrees of freedom for each fixed effect coefficient.

VarianceDecomposition

Cached matrices derived from the marginal covariance V(θ) and its derivatives.

Let θs denote the variance parameters, and write dVs[i] = ∂V/∂θ_i and d2Vs[i,j] = ∂²V/(∂θ_i∂θ_j).

This struct stores V, Vinv, θs, dVs, d2Vs, and the derived matrices

\[P_i = -X^\top V^{-1} dV_i V^{-1} X, \qquad Q_{ij} = X^\top V^{-1} dV_i V^{-1} dV_j V^{-1} X, \qquad R_{ij} = X^\top V^{-1} d^2V_{ij} V^{-1} X,\]

and R_factor. R_factor is a parameterization-dependent scalar applied to R_{ij} in the Kenward-Roger covariance bias correction (see compute_adjusted_covariance).

The definition of θ (and thus the rest of the calculated quantities here) depends on the chosen parameterization; see Unstructured and FactorAnalytic.

compute_adjusted_covariance(Φ, vd, W, method) -> Matrix

Compute the (possibly adjusted) covariance matrix of fixed effects.

For KenwardRoger, applies the bias correction:

\[\Phi_A = \hat{\Phi} + 2\hat{\Phi}\,B\,\hat{\Phi},\]

where the bias matrix $B$ accounts for the variability in estimating $\theta$:

\[B = \sum_{i,j} W_{ij} \bigl(Q_{ij} - P_i\hat{\Phi}P_j + R_{\text{factor}}\cdot R_{ij}\bigr).\]

Here $W$ is the asymptotic covariance of $\hat{\theta}$, and $P_i, Q_{ij}, R_{ij}$ are defined in VarianceDecomposition. The R_factor is -1/4 for FactorAnalytic parameterization and 1 for Unstructured.

For Satterthwaite, returns $\Phi$ unchanged (no bias correction).

Arguments

• Φ: Unadjusted covariance matrix $(X^\top V^{-1} X)^{-1}$
• vd: VarianceDecomposition containing derivative matrices
• W: Asymptotic covariance of variance parameters
• method: Adjustment method (KenwardRoger or Satterthwaite)

compute_dof(m::LinearMixedModelAdjusted, L) -> (ν, λ)

Compute degrees of freedom and F-statistic scaling factor for contrast L.

Dispatches to the specific method based on m.adj:

• For KenwardRoger: computes scaled F-statistic with $λ ≠ 1$
• For Satterthwaite: returns $λ = 1$ (no scaling)

Arguments

• m: A LinearMixedModelAdjusted from small_sample_adjust
• L: Contrast matrix ($p × q$) where $p$ = number of fixed effects

Returns

• ν: Denominator degrees of freedom
• λ: F-statistic scaling factor (1.0 for Satterthwaite)

compute_dof(::KenwardRoger, m, L)

Kenward-Roger denominator degrees of freedom and F-scaling factor.

For testing $H_0: L^\top\beta = 0$ with $q$ constraints, the standard Wald F-statistic $F = \beta^\top M \beta / q$ (where $M = L(L^\top\Phi L)^{-1}L^\top$) is scaled to $F^* = \lambda F$. The pair $(\nu, \lambda)$ is chosen so that $F^*$ approximately follows an $F(q, \nu)$ distribution.

Algorithm

1. Compute moment quantities:

   \[A_1 = \sum_{i,j} W_{ij} \operatorname{tr}(M\Phi P_i \Phi) \operatorname{tr}(M\Phi P_j \Phi)\]

   \[A_2 = \sum_{i,j} W_{ij} \operatorname{tr}(M\Phi P_i \Phi\, M\Phi P_j \Phi)\]

2. Compute intermediate terms:

   \[B = \frac{A_1 + 6A_2}{2q}, \quad g = \frac{(q+1)A_1 - (q+4)A_2}{(q+2)A_2}\]

   \[c_1 = \frac{g}{3q + 2(1-g)}, \quad c_2 = \frac{q-g}{3q + 2(1-g)}, \quad c_3 = \frac{q+2-g}{3q + 2(1-g)}\]

3. Compute expectation and variance of the scaled statistic:

   \[E^* = \frac{1}{1 - A_2/q}, \quad V^* = \frac{2}{q} \cdot \frac{1 + c_1 B}{(1-c_2 B)^2 (1-c_3 B)}\]

4. Match to F-distribution moments:

   \[\rho = \frac{V^*}{2(E^*)^2}, \quad \nu = 4 + \frac{q+2}{q\rho - 1}, \quad \lambda = \frac{\nu}{E^*(\nu - 2)}\]

compute_dof(::Satterthwaite, m, L)

Satterthwaite denominator degrees of freedom for contrast L.

For a scalar contrast $c^\top\beta$, the variance is $v = c^\top\Phi\,c$. The Satterthwaite approximation matches the first two moments of $\hat{v}$ to a scaled chi-squared distribution, yielding:

\[\nu = \frac{2v^2}{(\nabla_\theta v)^\top W (\nabla_\theta v)},\]

where the gradient is $\nabla_\theta v = [-c^\top\Phi\,P_i\,\Phi\,c]_i$ and $W$ is the asymptotic covariance of $\hat{\theta}$.

Multivariate Contrasts

For a matrix $L$ (multiple constraints), the algorithm:

1. Computes the eigendecomposition of $L^\top\Phi\,L$
2. Transforms to orthogonal contrasts $\tilde{L} = L \cdot \text{eigenvectors}$
3. Computes $\nu_i$ for each orthogonal contrast
4. Combines using: $\nu = 2 E_Q / (E_Q - q)$ where $E_Q = \sum_i \nu_i/(\nu_i - 2)$

Reference

Satterthwaite, F. E. (1946). "An Approximate Distribution of Estimates of Variance Components." Biometrics Bulletin, 2(6), 110-114.