Universal inference for variance components

TL;DR

This paper develops finite-sample valid confidence intervals for variance components, especially near boundaries, with applications in genetics and mixed models, offering faster algorithms than existing methods.

Contribution

It introduces novel universal inference methods for variance components near boundaries, improving accuracy and computational efficiency over existing approaches.

Findings

Confidence intervals are valid in finite samples for boundary cases.

Proposed algorithms are significantly faster than naive universal inference implementations.

Methods are demonstrated effective through simulations and real data examples.

Abstract

We consider universal inference in variance components models, focusing on settings where the parameter is near or at the boundary of the parameter set. Two cases, which are not handled by existing state-of-the-art methods, are of particular interest: (i) inference on a variance component when other variance components are near or at the boundary, and (ii) inference on near-unity proportions of variability, that is, one variance component divided by the sum of all variance components. Case (i) is relevant, for example, for the construction of componentwise confidence intervals, as often used by practitioners. Case (ii) is particularly relevant when making inferences about heritability in modern genetics. For both cases, we show how to construct confidence intervals that are uniformly valid in finite samples. We propose algorithms which, by exploiting the structure of variance components models, lead to substantially faster computing than naive implementations of universal inference. The usefulness of the proposed methods is illustrated by simulations and a data example with crossed random effects, which are known to be complicated for conventional inference procedures.

Figures

Tables

Table 1. Table 1: Maximum likelihood estimates of variance components and proportions of variability for the resistor data. Standard errors based on observed information are in parentheses.

	resistor	operator	error/total variance
${\sigma }^2$	0.00 (5.91)	50.0 (42.5)	40.3 (7.74)
$h^2$	0.00 (0.0654)	0.554 (0.217)	90.3 (43.4)

Table 2. Table 2: Confidence intervals for square-root variance components for the resistor data.

	Wald	profile	bootstrap	split LRT
${\sigma }_1$	(0.00, 3.40)	(0.00, 2.80)	(0.00, 2.72)	(0.00, 5.43)
${\sigma }_2$	(0.00, 11.5)	(3.55, 21.3)	(1.16, 12.3)	(2.38, 78.8)

Equations

$$Y\sim \mathrm{N}(0,{\sigma }_{K_1}^2+\cdots +{\sigma }_M^2{K}_M+{\sigma }_{M+1}^2{I}_n),$$

$$Y=Z_1{U}_1+\cdots +Z_M{U}_M+E,$$

$$h_m^2=\frac{{\sigma }_m^2}{{\sum }_{m=1}^{M+1}{\sigma }_m^2},m\in \{1,\dots ,M\}.$$

$$\Sigma =\Sigma (\theta )={\tau }^2\left\{h_1^2{K}_1+\cdots +h_M^2{K}_M+\left(1-\sum _{m=1}^M{h}_m^2\right)I_n\right\}\in {\mathbb{R}}^{n\times n}.$$

$$\frac{{\sup }_{\theta \in \Theta }{\mathcal{L}}_Y(\theta )}{{\sup }_{\theta \in {\Theta }_0}{\mathcal{L}}_Y(\theta )}.$$

$$T_n=\frac{{\mathcal{L}}_{Y_{(0)}|Y_{(1)}}({\hat{\theta }}_1)}{{\mathcal{L}}_{Y_{(0)}|Y_{(1)}}({\hat{\theta }}_0)}.$$

$$\Psi (h^2)=\sum _{m=1}^M{h}_m^2{K}_m+\left(1-\sum _{m=1}^M{h}_m^2\right)I_n,h^2=(h_1^2,\dots ,h_M^2{)}^{\mathrm{T}}.$$

$$l_Y(\theta )$$

$$=-\frac{1}{2}\left\{n\log ({\tau }^2)+\log \left|\Psi (h^2)\right|+\frac{1}{{\tau }^2}{Y}^{\mathrm{T}}\Psi (h^2{)}^{-1}Y\right\}.$$

$$Y_{(0)}\mid Y_{(1)}\sim \mathrm{N}\left({\Sigma }_{(01)}{\Sigma }_{(11)}^{-1}{Y}_{(1)},{\Sigma }_{(00)}-{\Sigma }_{(01)}{\Sigma }_{(11)}^{-1}{\Sigma }_{(10)}\right),$$

$$=-\frac{1}{2}\log |{\Sigma }_{(00)}-{\Sigma }_{(01)}{\Sigma }_{(11)}^{-1}{\Sigma }_{(10)}|$$

$$-\frac{1}{2}\{(Y_{(0)}-{\Sigma }_{(01)}{\Sigma }_{(11)}^{-1}{Y}_{(1)}{)}^{\mathrm{T}}({\Sigma }_{(00)}-{\Sigma }_{(01)}{\Sigma }_{(11)}^{-1}{\Sigma }_{(10)}{)}^{-1}$$

$$(Y_{(0)}-{\Sigma }_{(01)}{\Sigma }_{(11)}^{-1}{Y}_{(1)})\}.$$

$${\Sigma }_{(ij)}=\sum _{m=1}^M{h}_m^2{K}_{(ij)}^m+\mathbb{I}(i=j)\left(1-\sum _{m=1}^M{h}_m^2\right)I_{n_{(0)}},i,j\in \{0,1\},$$

$$\frac{\partial l_Y(\theta )}{\partial {\theta }_j}=-\frac{1}{2}\mathrm{tr}\left[\left\{\Sigma (\theta {)}^{-1}-\Sigma (\theta {)}^{-1}Y{Y}^{\mathrm{T}}\Sigma (\theta {)}^{-1}\right\}\frac{\partial \Sigma (\theta )}{\partial {\theta }_j}\right],$$

$$-\frac{1}{2}\mathrm{tr}\left[\Sigma (\theta {)}^{-1}\left\{I_n-{\tau }^{-2}Y{Y}^{\mathrm{T}}\Psi (h{)}^{-1}\right\}\Psi (h^2)\right]$$

$$=-\frac{1}{2}\mathrm{tr}\left[{\tau }^{-2}\left\{I_n-{\tau }^{-2}Y{Y}^{\mathrm{T}}\Psi (h{)}^{-1}\right\}\right],$$

$${\tilde{\tau }}_{(1)}^2(h^2)=n_{(1)}^{-1}{Y}_{(1)}^{\mathrm{T}}{\Psi }_{(11)}^{-1}(h^2)Y_{(1)}.$$

$${\frac{\partial l_{Y_{(1)}}(h^2,{\tau }^2)}{\partial h_m^2}|}_{h^2,{\tilde{\tau }}_{(1)}^2(h^2)}+{\frac{\partial l_{Y_{(1)}}(h^2,{\tau }^2)}{\partial {\tau }^2}\frac{\partial {\tau }^2(h^2)}{\partial h_m^2}|}_{h^2,{\tilde{\tau }}_{(1)}^2(h^2)}$$

$$={\frac{\partial l_{Y_{(1)}}(h^2,{\tau }^2)}{\partial h_m^2}|}_{h^2,{\tilde{\tau }}_{(1)}^2(h^2)}+0,$$

$$-\frac{1}{2}\mathrm{tr}\left[\left\{{\Sigma }_{(11)}^{-1}(\theta )-{\Sigma }_{(11)}^{-1}(\theta )Y_{(1)}{Y}_{(1)}^{\mathrm{T}}{\Sigma }_{(11)}^{-1}(\theta )\right\}(K_{(11)}^m-I_{n_{(1)}})\right].$$

$${\Sigma }_{(0\mid 1)}(\theta )={\Sigma }_{(00)}(\theta )-{\Sigma }_{(01)}(\theta ){\Sigma }_{(11)}^{-1}(\theta ){\Sigma }_{(10)}(\theta ).$$

$${\tilde{\tau }}_{(0\mid 1)}^2$$

$$=n_{(0)}^{-1}(Y_{(0)}-{\Psi }_{(01)}{\Psi }_{(11)}^{-1}{Y}_{(1)}{)}^{\mathrm{T}}{\Psi }_{(0\mid 1)}^{-1}(Y_{(0)}-{\Psi }_{(01)}{\Psi }_{(11)}^{-1}{Y}_{(1)}).$$

$$l_{Y_{(0)}\mid Y_{(1)}}(\theta )=l_Y(\theta )-l_{Y_{(1)}}(\theta ).$$

$${\check{h}}_{(0\mid 1)}^2\in \underset{h^2\in {\mathbb{R}}^M}{\mathrm{argmax}}{l}_{Y_{(0)}\mid Y_{(1)}}\{h^2,{\tilde{\tau }}_{(0\mid 1)}^2(h^2)\}.$$

$${\check{h}}_{(1)}^2\in \underset{h^2\in {\mathbb{R}}^M}{\mathrm{argmax}}{l}_{Y_{(1)}}\{h^2,{\tilde{\tau }}_{(1)}^2(h^2)\}.$$

$$\Sigma (\theta )={\tau }^2\{h_M^2{K}_M+(1-h_M^2)I_n\}={\tau }^{2}O\{h_M^2\Lambda +(1-h_M^2)I_n\}{O}^{\mathrm{T}}.$$

$$-\frac{1}{2}\sum _{k:Y_k\in Y_{(0)}}\left\{\log ({\tau }^2)+\log (h_M^2{\lambda }_k+1-h_M^2)\frac{Y_k^2}{{\tau }^2(h_M^2{\lambda }_k+1-h_M^2)}\right\},$$

$$O^{\mathrm{T}}{K}_{m}O={\Lambda }_m,m\in \{1,\dots ,M\},$$

$$-\frac{1}{2}\sum _{k:Y_k\in Y_{(i)}}\left\{\log \left(\sum _{m=1}^M{\sigma }_m^2{\lambda }_{mk}+{\sigma }_{M+1}^2\right)+\frac{Y_k^2}{{\sum }_{m=1}^M{\sigma }_m^2{\lambda }_{mk}+{\sigma }_{M+1}^2}\right\},$$

$$=-\frac{1}{2}\sum _{k:Y_k\in Y_{(i)}}\left\{\frac{{\lambda }_{mk}}{{\sum }_{r=1}^M{\sigma }_r^2{\lambda }_{rk}+{\sigma }_{M+1}^2}-\frac{{\lambda }_{mk}{Y}_k^2}{{\left({\sum }_{r=1}^M{\sigma }_r^2{\lambda }_{rk}+{\sigma }_{M+1}^2\right)}^2}\right\},$$

$$l_{Y_{(i)}}({\sigma }^2)=-\frac{1}{2}\sum _{k:Y_k\in Y_{(i)}}\left\{\log \left({\sigma }_{(k)}^2{\lambda }_{(k)}+{\sigma }_{M+1}^2\right)+\frac{Y_k^2}{{\sigma }_{(k)}^2{\lambda }_{(k)}+{\sigma }_{M+1}^2}\right\},$$

Full text

1 Introduction

Variance components models are used routinely in a wide variety of scientific applications. Often times, multiple sources of variation are present, in which case practitioners want to understand the degree of total variation attributable to each source. In epidemiology, for example, researchers want to understand how a trait’s variability is affected by additive genetic effects and the environment (Heckerman et al.,, 2016), or by additive genetic effects and gene-environment interactions (Pazokitoroudi et al.,, 2024). Similarly, in statistical genetics, genetic effects can be partitioned by chromosome (Yang et al.,, 2011), partitioned into purely additive effects versus genetic interaction effects (Bloom et al.,, 2015; Vitezica et al.,, 2013), partitioned into additive and dominance effects of SNP markers (Da et al.,, 2014), among many other partitioning approaches (e.g., see Runcie and Crawford,, 2019).

More generally, hierarchical and multilevel mixed models aim to quantify the degree of a random variable’s variability that can be attributed to distinct sources of variation (Goldstein,, 2011; Kreft and De Leeuw,, 1998; Lee and Nelder,, 1996; Rasbash and Goldstein,, 1994). For example, the variability in student test scores may be attributed to variation arising from classroom effects, school-level effects, or community-wide effects.

One widely used model for capturing multiple sources of variation is the variance components model. A variance components model assumes a vector $Y\in\mathbb{R}^{n}$ satisfies, for some known symmetric and positive semi-definite $K_{1},\dots,K_{M}$, $M$ a positive integer,

[TABLE]

where $\sigma^{2}_{1},\dots,\sigma^{2}_{M}$ are the variance components. Here, we assume $\mathsf{E}(Y)=0$ for simplicity but later allow $\mathsf{E}(Y)=X\beta$ for known $X\in\mathbb{R}^{n\times p}$ and unknown $\beta\in\mathbb{R}^{p}$. The distribution in (1) sometimes results from a random effects model:

[TABLE]

where, independently for each $m\in\{1,\dots,M\}$, random effects satisfy $U_{m}\sim\mathrm{N}(0,\sigma^{2}_{m}I_{q_{m}})$ for some $q_{m}\in\{1,\dots,m\}$, and $Z_{m}\in\mathbb{R}^{n\times q_{m}}$. The error term $E\sim\mathrm{N}(0,\sigma^{2}_{M+1}I_{n})$ is independent of the random effects. Then, for every $m\leq M$, $K_{m}=Z_{m}Z_{m}^{\mathrm{\scriptscriptstyle T}}$, so that the rank of $K_{m}$ is at most $q_{m}\leq n$. To avoid degenerate distributions, we will typically assume the error variance is nonnegative, $\sigma^{2}_{M+1}>0$.

Our focus is hypothesis tests, or inferences more generally, for the variance components and the proportions

[TABLE]

The parameter $h_{m}^{2}$ is often interpreted as the proportion of variability attributable to sources encoded by $K_{m}$, $m\in\{1,\dots,M\}$. In statistical genetics, for example, $K_{m}$ can encode the genetic similarly of individuals’ $m$th chromosome. Then, $h^{2}_{m}$ is the proportion of variance in the outcome explained by the $m$th chromosome’s SNP genotypes (Yang et al.,, 2011). To facilitate inference on the $h^{2}_{m}$ we consider a reparameterization of (1) in terms of $\theta=(h^{2}_{1},\dots,h^{2}_{M},\tau^{2})^{{\mathrm{\scriptscriptstyle T}}}$, where $\tau^{2}=\sum_{m=1}^{M+1}\sigma^{2}_{m}$. Then, the parameter set $\Theta\subseteq\mathbb{R}^{M+1}$ is the set of $\theta$ such that $\sum_{m=1}^{M}h_{m}^{2}<1$, $h^{2}_{m}\geq 0$ for all $m$, and $\tau^{2}>0$. With these parameters $Y$ is multivariate normal with mean zero and covariance matrix

[TABLE]

Inference on the $h^{2}_{m}$ is complicated in general because one or more of them are often close to zero. That is, the true parameter is often at or near the boundary of the parameter set. When $M=1$, so that there is a single $h^{2}_{m}$, there are methods based on inverting score test-statistics (Zhang et al.,, 2025) and simulation-based methods (Crainiceanu and Ruppert,, 2004; Schweiger et al.,, 2018, 2016). However, the supporting theory for the simulation-based methods is not applicable when $M>1$, and the score-based confidence intervals have non-nominal coverage probability when there are nuisance parameters near the boundary. For example, a score-based confidence interval for $h^{2}_{1}$ often has non-nominal coverage probability if another $h^{2}_{m}$, $m\neq 1$, is close to zero or one. Figure 1 illustrates this issue in a setting where $M=2$; the parameter of interest is $h_{1}^{2}$, whose true value is zero; and $h_{2}^{2}$ is a nuisance parameter whose true value is on the horizontal axis. Clearly, the coverage probability for the confidence interval for $h_{1}^{2}$ is affected by the true value of the nuisance parameter, especially when it is near one. Somewhat informally, the issue is that test-statistics for $h_{1}^{2}$ depend on a constrained maximum likelihood estimator of $h_{2}^{2}$, and that estimator behaves irregularly near the boundary.

The issues are easier to deal with when there are no nuisance parameters as, then, score-based test-statistics can be evaluated at the null hypothesis parameter vector, no estimation needed (Zhang et al.,, 2025; Ekvall and Bottai,, 2025). However, even without nuisance parameters, state-of-the art methods are unreliable when a $h_{m}^{2}$ is near unity. The main reason is that points where $\sum_{m=1}^{M}h_{m}^{2}=1$ are often hard boundary points, in the sense that the likelihood cannot be extended to such points, while points where a $h^{2}_{m}=0$ are soft boundary points (Elkantassi et al.,, 2023).

Here, we consider methods for componentwise inference for variance components in the presence of nuisance parameters. The methods are based on universal inference, and in particular split likelihood ratio tests (Wasserman et al.,, 2020). Consequently, the proposed tests and confidence intervals are uniformly valid in finite samples, regardless of how close to a soft or hard boundary the parameter is. Thus, even when there are no nuisance parameters, the methods studied here can be preferable to common ones, which are often motivated by asymptotic theory. We describe universal inference and its application to our setting in more detail in Section 2. Briefly, the method relies on data splitting, where one set of data is used to estimate parameters and the other to carry out the test. In our setting, where there is a potentially complicated dependence structure, more care is needed when randomly splitting the data than in settings with independent and identically distributed observations. A main concern with universal inference is a lack of power compared to classical parametric methods, which we address by using a randomized version of the split likelihood ratio test (Ramdas and Manole,, 2023).

There are several computational challenges with applying universal inference to (1), especially when $n$ is large. To address such issues we develop efficient algorithms for several special cases of interest. For example, in a setting with crossed random effects, which are known to complicate both computation and theory, we decrease the required time by several orders of magnitude compared to a naive implementation.

2 Randomized Split Likelihood Ratio Test

Let $Y\in\mathbb{R}^{n}$ be a random vector with density $f_{\theta^{*}}(y)$, $\theta^{*}\in\Theta$ being the true value of the parameter. Let $\mathcal{L}_{Y}(\theta)=f_{\theta}(Y)$ be the (random) likelihood function; that is, the density of $Y$, evaluated at $Y$, considered as a function of $\theta$ with domain $\Theta$. In our setting the density and likelihood correspond to (1), but universal inference applies more generally. Suppose we wish to test test $H_{0}:\theta^{*}\in\Theta_{0}\subseteq\Theta$ versus $H_{A}:\theta^{*}\in\Theta\setminus\Theta_{0}$. Recall that a conventional likelihood ratio test (LRT) can be based on the statistic

[TABLE]

Under well-known regularity conditions, two times the logarithm of this statistic has an asymptotic chi-square distribution, which can be used to construct tests with asymptotically correct size. However, such regularity conditions do not hold in our setting, so we instead consider the split LRT.

Let $Y_{(0)}\in\mathbb{R}^{n_{(0)}}$ and $Y_{(1)}\in\mathbb{R}^{n_{(1)}}$ be a random partition of $Y$, where $n_{(0)}+n_{(1)}=n$. That is, every element of $Y$ is in exactly one of the $Y_{(i)},i\in\{0,1\}$. The randomization is done independently of $Y$, but it need not be uniform on the set of possible partitions of given sizes. Let $\mathcal{L}_{Y_{(1)}}(\theta)$ be the likelihood based on $Y_{(1)}$ and $\mathcal{L}_{Y_{(0)}|Y_{(1)}}(\theta)$ likelihood based on the conditional distribution of $Y_{(0)}$ given $Y_{(1)}$. Let also $\hat{\theta}_{1}\in\operatorname*{arg\,max}_{\theta\in\Theta}\mathcal{L}_{Y_{(1)}}(\theta)$ and $\hat{\theta}_{0}\in\operatorname*{arg\,max}_{\theta\in\Theta_{0}}\mathcal{L}_{Y_{(0)}|Y_{(1)}}(\theta)$, assuming they exist. Then the split likelihood-ratio statistic is

[TABLE]

The test that rejects $H_{0}$ if $T_{n}>1/\alpha$ is valid at level $\alpha\in(0,1)$ (Wasserman et al.,, 2020); that is, the size of the test is at most $\alpha$, and hence confidence regions obtained by inverting the test have coverage probability at least $1-\alpha$. In fact, the test remains valid, and has greater power, if rejection is instead based on comparing to a uniform random variable (Ramdas and Manole,, 2023). We state these facts formally in the following known result and provide a proof in the Appendix for completeness. It will be important for later to note that the same result and proof holds if, in (4), $\hat{\theta}_{1}$ is replaced by any other estimator based on $Y_{(1)}$ only. Similarly, $\hat{\theta}_{0}$ can be replaced by any $\tilde{\theta}$ such that $\mathcal{L}_{Y_{(0)}\mid Y_{(1)}}(\tilde{\theta})\geq\mathcal{L}_{Y_{(0)}\mid Y_{(1)}}(\hat{\theta}_{0})$.

Theorem 1**.**

The split likelihood ratio test that rejects if $T_{n}>1/\alpha$ is a uniformly valid level $\alpha$ test, i.e., for any $\Theta_{0}\subseteq\Theta$ and $\theta^{*}\in\Theta_{0}$, it holds that $\mathrm{P}_{\theta^{*}}(T_{n}>1/\alpha)\leq\alpha$. Moreover, the randomized split likelihood ratio test that rejects if $T_{n}>U/\alpha$, with $U\sim\mathrm{U}(0,1)$ independently of $T_{n}$, is also uniformly valid, and it is more powerful than the split likelihood ratio test.

2.1 Application to Variance Components

Assume (1) with the parameterization in (3) and let

[TABLE]

Then, ignoring additive terms not depending on $\theta$,

[TABLE]

Because the partitioning is independent of $Y$, $Y_{(0)}$ and $Y_{(1)}$ are jointly multivariate normal. Thus, the marginal log-likelihood for $Y_{(1)}$ is similar to that for $Y$, with $Y_{(1)}$ and $\Sigma_{(11)}$ in place of $Y$ and $\Sigma$, respectively. Also,

[TABLE]

where $\Sigma_{(ij)}=\Sigma_{(ij)}(\theta)=\mathsf{E}_{\theta}(Y_{(i)}Y_{(j)}^{\mathrm{\scriptscriptstyle T}})-\mathsf{E}_{\theta}(Y_{(i)})\mathsf{E}_{\theta}(Y_{(j)})^{\mathrm{\scriptscriptstyle T}}$, $i,j\in\{0,1\}$. Thus, the conditional log-likelihood $l_{Y_{(0)}\mid Y_{(1)}}(\theta)$ is

[TABLE]

We can write the $\Sigma_{(ij)}$ in terms of the $K_{m}$ as follows. For each $m\in\{1,\dots,M\}$, let $K^{m}_{(ij)}$ be the matrix obtained from $K_{m}$ by keeping only the rows with indices corresponding to observations in $Y_{(i)}$ and columns with indices corresponding to observations in $Y_{(j)}$. Then

[TABLE]

where $\mathbb{I}(\cdot)$ is an indicator function and $n_{(0)}$ the number of elements in $Y_{(0)}$.

Several challenges with implementing universal inference are evident from (5). In general, finding the maximizers $\hat{\theta}_{0}$ and $\hat{\theta}_{1}$ is nontrivial and requires numerical optimization. To find $\hat{\theta}_{1}$, we use off-the-shelf, gradient-based methods to maximize $l_{Y_{(1)}}(\theta)$. To that end, note that for any element of $\theta$, say $\theta_{j}$,

[TABLE]

and similarly for $l_{Y_{(1)}}(\theta)$ but with $Y_{(1)}$ and $\Sigma_{(11)}$ in place of $Y$ and $\Sigma$, respectively. When $j=M+1$, so that $\theta_{j}=\tau^{2}$, (6) is

[TABLE]

which vanishes if evaluated at $\tau^{2}=n^{-1}Y^{\mathrm{\scriptscriptstyle T}}\Psi^{-1}(h^{2})Y$. It is routine to show this stationary point is in fact a global partial maximizer. Thus, an algorithm for finding $\hat{\theta}_{1}$ can alternate between updating optimization variables corresponding to $\tau^{2}$ and $h^{2}$, with the former update being available in closed form. By contrast, first order conditions for $h^{2}$ based on (6) cannot in general be solved analytically. Thus, we update $h^{2}$ using a gradient-based step. Equivalently, we use gradient-based methods to maximize the profile log-likelihood $h^{2}\mapsto l_{Y_{(1)}}\{h^{2},\tilde{\tau}^{2}_{(1)}(h^{2})\}$, where

[TABLE]

The derivatives of the profile log-likelihood, $\partial l_{Y_{(1)}}\{h^{2},\tilde{\tau}^{2}(h)\}/\partial h_{m}^{2}$, are

[TABLE]

where the last equality is due to $\tilde{\tau}_{(1)}^{2}(h^{2})$ satisfying the first order condition for an interior partial maximizer; these calculations can be formalized (Milgrom and Segal,, 2002). Thus, updating $\tau^{2}$ and then updating $h^{2}$ using the gradient of the log-likelihood, is equivalent to updating $h^{2}$ using the gradient of the profile log-likelihood. Using (6), $\partial l_{Y_{(1)}}(h^{2},\tau^{2})/\partial h_{m}^{2}$ is

[TABLE]

Similar arguments apply to the problem of finding $\hat{\theta}_{0}$. Let

[TABLE]

Note that $\Sigma_{(0\mid 1)}(\theta)=\tau^{2}\Psi_{(0\mid 1)}(h^{2})$, where $\Psi_{(0\mid 1)}$ is defined like $\Sigma_{(0\mid 1)}$ but replacing every $\Sigma$ by $\Psi$. Thus, omitting the $h^{2}$ argument for simplicity, we have the partial maximizer

[TABLE]

To apply gradient-based methods to the profile, conditional log-likelihood $h^{2}\mapsto l_{Y_{(0)}\mid Y_{(1)}}\{h^{2},\tau^{2}_{(0\mid 1)}(h^{2})\}$, observe

[TABLE]

Derivatives of the two terms on the right-hand side can be obtained as in (6) and (7).

In general, the objective functions we have discussed are nonconvex, and standard algorithms can be computationally expensive due to matrix decompositions needed to deal with the inverses, scaling approximately as $n^{3}$. Additionally, even when the parameter is identifiable in the distribution for $Y$, it can be unidentifiable, or nearly so, in the conditional distribution of $Y_{(0)}$ given $Y_{(1)}$, as we will see examples of later.

The constraint $\sum_{m=1}^{M}h_{m}^{2}<1$ is also nontrivial in general. An exception is the case where $M=1$ and $K^{1}_{(11)}$ is singular, because in that case any $h^{2}_{1}\leq 0$ would lead to a $\Sigma_{(11)}(\theta)$ that is not positive definite, and hence an undefined or vanishing likelihood. Thus, for that case, the log-determinant term in the multivariate normal log-likelihood acts as a barrier. For the other cases, implementation is made easier by the fact that Theorem 1 continues to hold if we ignore the constraint when finding $\hat{\theta}_{0}$ and $\hat{\theta}_{1}$. Specifically, we may replace $\hat{\theta}_{0}$ by $\check{\theta}_{0}=\{\check{h}_{(0\mid 1)}^{2},\tilde{\tau}^{2}_{(0\mid 1)}(\check{h}_{(0\mid 1)}^{2})\}^{\mathrm{\scriptscriptstyle T}}$, where

[TABLE]

As noted before Theorem 1, validity is retained since $\mathcal{L}_{Y_{(0)}\mid Y_{(1)}}(\check{\theta}_{0})\geq\mathcal{L}_{Y_{(0)}\mid Y_{(1)}}(\hat{\theta}_{0})$. Similarly, we can replace $\hat{\theta}_{1}$ by $\check{\theta}_{1}=\{\check{h}_{(1)}^{2},\tilde{\tau}^{2}_{(1)}(\check{h}^{2})\}^{\mathrm{\scriptscriptstyle T}}$, where

[TABLE]

Validity is retained since $\check{\theta}_{1}$ is a function of $Y_{(1)}$ only. Both $\check{\theta}_{0}$ and $\check{\theta}_{1}$ can be obtained using unconstrained optimization.

In the following section, we discuss special cases of interest where computation can be made more efficient by using additional structure.

3 Testing Variance Components

3.1 Boundary Points

Suppose we wish to test whether all but one variance component are zero; without loss of generality, $\sigma^{2}_{M}$ can be non-zero under the null hypothesis. Equivalently, $h^{2}_{M}$ can be non-zero under the null hypothesis. Thus, in the parameterization given by (3), $\Theta_{0}$ is the set of $\theta=(h_{1}^{2},\dots,h_{M}^{2},\tau^{2})^{\mathrm{\scriptscriptstyle T}}$ such that $h_{1}^{2}=\cdots=h_{M-1}^{2}=0$, $h_{M}^{2}<1$, and $\tau^{2}>0$. This setting is perhaps most natural when $M=2$, in which case it corresponds to testing one variance component with the other unconstrained. Even that special case is challenging for existing methods when the unconstrained parameter is near the boundary.

The structure of $\Theta_{0}$ enables substantial computational gains compared to a naive implementation. In particular, let $K_{M}=O\Lambda O^{\mathrm{\scriptscriptstyle T}}$ by eigendecomposition and note that for any $\theta\in\Theta_{0}$,

[TABLE]

Consequently, upon replacing $Y$ and $K_{m}$ by $O^{\mathrm{\scriptscriptstyle T}}Y$ and $O^{\mathrm{\scriptscriptstyle T}}K_{m}O$, respectively, $m\in\{1,\dots,M\}$, we may assume without loss of generality that $K_{M}=\Lambda$. With this assumption, $\Sigma(\theta)$ is diagonal for $\theta\in\Theta_{0}$, which simplifies computation of $\hat{\theta}_{0}$.

Specifically, when $K_{M}$ is diagonal, $\Sigma_{(01)}(\theta)$ is a matrix of zeros for $\theta\in\Theta_{0}$. Thus, for such $\theta$, $l_{Y_{(0)}\mid Y_{(1)}}(\theta)=l_{Y_{(0)}}(\theta)$, which equals

[TABLE]

where $\lambda_{k}$ is the $k$th element of $\Lambda$ and, with a slight abuse of notation, $Y_{k}\in Y_{(0)}$ means the $k$th element of $Y$ is an element of $Y_{(0)}$. Thus, finding $\hat{\theta}_{0}$ reduces to a one-dimensional optimization problem over the interval $[0,1)$ with an easy-to-compute derivative. This problem is substantially simpler than maximizing (5) in general, without diagonalization.

3.2 Shared Eigenvectors

In some settings the eigenvectors of $K_{m}$ can be chosen not to depend on $m$; that is, one can find an orthogonal $O$ such that

[TABLE]

where $\Lambda_{m}$ is a diagonal matrix with the eigenvalues of $K_{m}$ on the diagonal.

One such example is when the $K_{m}$ model variation in orthogonal directions, or more formally, $K_{m}K_{\ell}=0$ for every $m\neq\ell$. Equivalently, by symmetry, every column of $K_{m}$ is orthogonal to every column of $K_{\ell}$. This claim follows from well-known facts about commuting, symmetric matrices; for completeness we give a more direct and constructive proof in the Appendix.

Theorem 2**.**

If symmetric matrices $K_{1},\dots,K_{M}$ satisfy $K_{m}K_{\ell}=0$ for every $m\neq\ell$, then they satisfy (8).

Assuming (8), $O^{\mathrm{\scriptscriptstyle T}}Y$ has the multivariate normal distribution in (1), but with $\Lambda_{m}$ in place of $K_{m}$, $m\in\{1,\dots,M\}$. Thus, replacing $Y$ by $O^{\mathrm{\scriptscriptstyle T}}Y$ if needed, we may assume every $K_{m}=\Lambda_{m}$ is diagonal. We make this assumption for the remainder of the section. To set us up for a motivating example in the next section, we use the parameterization with $\sigma^{2}=(\sigma_{1}^{2},\dots,\sigma_{M+1}^{2})^{\mathrm{\scriptscriptstyle T}}\in\Omega=[0,\infty)^{M}\times(0,\infty)$ and with some abuse of notation write $\Sigma(\sigma^{2})$ for the covariance matrix of the multivariate normal distribution in (1).

Since $\Sigma(\sigma^{2})$ is diagonal for every $\sigma^{2}\in\Omega$, $l_{Y_{(i)}}(\sigma^{2})=$ is, for $i\in\{0,1\}$,

[TABLE]

where $\lambda_{mk}$ is the $k$th diagonal element of $\Lambda_{m}$, $m\in\{1,\dots,M\}$. Moreover, $Y_{(0)}$ and $Y_{(1)}$ are independent, and hence $l_{Y_{(0)}\mid Y_{(1)}}=l_{Y_{(0)}}$, which simplifies finding $\hat{\theta}_{0}$ compared to the general case. The derivatives $\partial l_{Y_{(i)}}(\sigma^{2})/\partial\sigma^{2}_{m}$ needed for gradient-based methods are, for $i\in\{0,1\}$,

[TABLE]

where $m\in\{1,\dots,M+1\}$ and $\lambda_{(M+1)k}=1$ for all $k$.

When $K_{m}K_{\ell}=0$ further simplifications are possible: not only may we assume $K_{m}=\Lambda_{m}$, but after doing so it holds that for each $k$, there is at most one $m\in\{1,\dots,M\}$ for which $\lambda_{mk}\neq 0$. Let $\lambda_{(k)}$ be that $\lambda_{mk}$, with $\lambda_{(k)}=0$ if no such $m$ exists, and let $\sigma^{2}_{(k)}$ be the corresponding $\sigma^{2}_{m}$. Then, for $i\in\{0,1\}$

[TABLE]

and similarly for its derivatives.

3.3 Crossed Random Effects

Crossed random effects are known to be challenging both for theory and computation (Jiang,, 2013; Ekvall and Jones,, 2020; Papaspiliopoulos et al.,, 2020; Ghosh et al.,, 2022; Lyu et al.,, 2024; Jiang et al.,, 2024; Jiang,, 2025; Ekvall and Bottai,, 2025). However, more efficient computing is possible by using a connection to the setting with shared eigenvectors, in the sense of (8). To introduce the setting, suppose momentarily that responses are naturally organized as a matrix $(Y_{ij})$ with $n_{1}$ rows and $n_{2}$ columns, with observations in the same row or column potentially dependent. A simple model with two crossed random effects is

[TABLE]

where $U_{1i}\sim\mathrm{N}(0,\sigma_{1}^{2})$, $U_{2j}\sim\mathrm{N}(0,\sigma^{2}_{2})$, and $E_{ij}\sim\mathrm{N}(0,\sigma^{2}_{3})$, independently for all $i\in\{1,\dots,n_{1}\}$ and $j\in\{1,\dots,n_{2}\}$. The $U_{1i}$ can be interpreted as row effects and the $U_{2j}$ as column effects.

More generally, suppose there are $M$ crossed random effects, with $n_{m}$ observations along the $m$th dimension, $m\in\{1,\dots,M\}$. Define the index set

[TABLE]

and suppose that for ${(j)}\in\mathcal{J}$,

[TABLE]

where $U_{mj_{m}}\sim\mathrm{N}(0,\sigma^{2}_{m})$ and $E_{{(j)}}\sim\mathrm{N}(0,\sigma^{2}_{M+1})$ are independent for all $m\in\{1,2,\dots,M\}$ and $j_{m}\in\{1,2,\dots,n_{m}\}$. Following Ekvall and Bottai, (2025), we can write this model as $Y=ZU+E$ by letting

[TABLE]

$m\in\{1,\dots,M\}$. Accordingly,

[TABLE]

Now, to explore the connection to (8), define the projection matrices $P_{m}=1_{n_{m}}1^{T}_{n_{m}}/n_{m}$, and $R_{m}^{I}=P_{1}\otimes\cdots\otimes P_{m-1}\otimes I_{n_{m}}\otimes P_{m+1}\otimes\cdots\otimes P_{M}$. Let also $w_{m}=\prod_{k\neq m}n_{k}$. Then the covariance matrix of $Y$ is

[TABLE]

Define $Q_{m}=I_{n_{m}}-P_{m}$, $R_{m}^{Q}=P_{1}\otimes\cdots\otimes P_{m-1}\otimes Q_{m}\otimes P_{m+1}\otimes\cdots\otimes P_{M}$, and $R^{P}=P_{1}\otimes\cdots\otimes P_{M}$. Then $R_{m}^{I}=R^{P}+R_{m}^{Q}$ and, consequently,

[TABLE]

Now note, for any $m\neq\ell$, by properties of Kronecker products, $P_{m}Q_{m}=0$ and $R^{P}R_{m}^{Q}=R_{m}^{Q}R_{\ell}^{Q}=0$. Thus, by Theorem 2, there is an orthogonal $O$ such that $O^{\mathrm{\scriptscriptstyle T}}R_{m}^{Q}O=D^{Q}_{m}$ and $O^{\mathrm{\scriptscriptstyle T}}R^{P}O=D^{P}$, for diagonal $D_{m}^{Q}$, $m\in\{1,\dots,M\}$, and $D^{P}$. Thus, upon replacing $Y$ by $O^{\mathrm{\scriptscriptstyle T}}Y$, we may assume $\Sigma(\sigma^{2})$ is diagonal. Specifically,

[TABLE]

where $\Lambda_{m}=w_{m}(D_{m}^{Q}+D^{P})$. Thus, inference can be based on (9). Note, however, that $\Lambda_{m}\Lambda_{\ell}=w_{m}w_{\ell}D^{P}\neq 0$ for $m\neq\ell$, so the further simplifications discussed following (9) are not applicable.

The columns in the matrix $O$ can be computed relatively cheaply using that eigenvectors of $R_{m}^{Q}$, $m\in\{1,\dots,M\}$, are Kronecker products of eigenvectors of the matrices $P_{1},\dots,P_{m-1}$, $Q_{m}$, $P_{m+1},\dots,P_{M}$. Specifically, there are $n_{m}-1$ orthonormal eigenvectors corresponding to the eigenvalue one. Each of these is in the form $1_{m_{-}}\otimes v\otimes 1_{m_{+}}/\sqrt{w_{m}}$, where $v$ is an eigenvector of $Q_{m}$ corresponding to the eigenvalue one, and $1_{m_{-}}$ and $1_{m_{+}}$ are vectors of ones of lengths $m_{-}=\sum_{k=1}^{m-1}n_{k}$ and $m_{+}=\sum_{k=m+1}^{M}n_{k}$, respectively. Additionally, one column of $O$ can be the eigenvector $1_{n}/\sqrt{n}$ of $R^{P}$. Thus, we have $1+\sum_{m=1}^{M}(n_{m}-1)$ columns of $O$; the remaining can be obtained by completing the orthonormal basis in any way.

3.4 Approximate Diagonalization

In some settings (8) may not hold exactly but approximately. That is, for an orthogonal $O$, $\|\Lambda_{m}-O^{\mathrm{\scriptscriptstyle T}}K_{m}O\|$ is small for every $m$, and hence, intuitively, $O^{\mathrm{\scriptscriptstyle T}}Y$ has a covariance matrix that is close to diagonal. Similarly, $\|K_{m}K_{\ell}\|$ can be small for every $m\neq\ell$. In either case, it may be useful to replace $K_{m}$ by an approximation $\tilde{K}_{m}$, $m\in\{1,\dots,M\}$, such that the $\tilde{K}_{m}$ are jointly diagonalizable as in (8). To formalize, let us consider an example that we will examine in simulations.

Suppose $M=2$ for simplicity, and that the $q_{m}<n$ greatest eigenvalues of $K_{m}$ are much larger than the trailing $n-q_{m}$, $m\in\{1,\dots,M\}$. Suppose also $K_{1}=\Lambda_{1}$ is diagonal, with the diagonal elements of $\Lambda_{1}$ sorted in decreasing order. Suppose also that $K_{2}=O_{2}\Lambda_{2}O_{2}^{\mathrm{\scriptscriptstyle T}}$, with the eigenvalues of $\Lambda_{2}$ sorted in decreasing order, and $O_{2}=\operatorname{bdiag}(I_{q_{2}},O_{(2)})$ for an orthogonal $O_{(2)}$ that is not the identity matrix. For simplicity we assume that the diagonal elements of each are distinct. More specifically, the $q_{m}$ first elements of $\Lambda_{m}$ are evenly spaced on, say, $(a_{2},a_{3})$, and for some constant $c>0$, the trailing $n-q_{m}$ eigenvalues are evenly spaced numbers on $(0,a_{1})$ divided by a constant $c\geq 1$. The larger $c$ is, the better $K_{m}$ is approximated by the $\tilde{K}_{m}$ that sets small eigenvalues to zero. Notably, for those $\tilde{K}_{m}$, (8) holds with $O=I_{n}$.

4 Simulations

Fig.2 shows, in a setting with $M=2$ variance components, Monte Carlo estimates of the coverage probabilities for the score-based confidence interval discussed in the introduction and the proposed randomized split LRT confidence interval. The confidence intervals are for $h_{1}^{2}$, with $h^{2}_{2}$ a nuisance parameter. In the simulations, the matrices $K_{1}$ and $K_{2}$ were set to diagonal matrices with the eigenvalues of autoregressive correlation matrices with correlation parameters 0.95 and 0.5, respectively. That is, $K_{1}$ was diagonal with the eigenvalues of the $n\times n$ matrix $(0.95^{|i-j|})$, and similarly for $K_{2}$. We set $n=300$, $h_{2}^{2}=0$, and $\tau^{2}=1$. The value of $h_{1}^{2}$ is on the horizontal axis. Monte Carlo estimates, and the corresponding confidence bands, are based on 10,000 replications.

Fig. 2 indicates that, when the nuisance parameter $h_{2}^{2}$ is near one, the score-based confidence interval for $h_{1}^{2}$ is invalid. The distortion is substantial, with coverage as low as 0.88 for extreme parameter values. By contrast, the proposed confidence interval, while conservative, is everywhere valid, as guaranteed by theory. The actual coverage probability of the proposed interval is around 0.975 regardless of the value of the nuisance parameter. For some values of the nuisance parameter this is substantially higher than the coverage probability of the score-based interval, but when $h_{2}^{2}$ is near zero the two intervals have similar coverage probabilities.

Fig. 3 shows Mote Carlo estimates of rejection probabilities for the split likelihood ratio test in three different scenarios, all with $M=2$ variance components. For all three scenarios, $K_{1}$ and $K_{2}$ were constructed as described in Sec. 3.4, with $a_{1}=a_{2}=5$, $a_{3}=10$, and $c=100$. Thus, the $q_{m}$th eigenvalue of $K_{m}$ in decreasing order is $5$ while the $(q_{m}+1)$th is $5/100$. To ensure a non-identity $O_{(2)}$, we drew it, before starting the simulations, uniformly on the Stiefel manifold as the left singular vectors of a $(n-q_{2})\times(n-q_{2})$ matrix with independent standard normal entries. We set $q_{1}=100\neq q_{2}=120$ to ensure identifiability. The true $h_{1}^{2}$, $h_{2}^{2}$, and $\tau^{2}$ were, respectively, 0, 0, and 1; the null hypothesis value of $h_{1}^{2}$ is on the horizontal axis.

The left plot in Fig. 3 shows the proposed test is conservative–the rejection probability at $h_{1}^{2}=0$ is below the nominal level. As the null hypothesis value moves away from the truth, the rejection probability increases monotonically, as expected.

In the middle plot of Fig. 3, the proposed test is implemented with approximations $\tilde{K}_{1}$ and $\tilde{K}_{2}$ in place of the true $K_{1}$ and $K_{2}$ used to generate the data. Specifically, as discussed in Sec. 3.4, $\tilde{K}_{m}$ sets the $n-q_{m}$ smallest eigenvalues to zero, which leads to jointly diagonalizable $\tilde{K}_{1}$ and $\tilde{K}_{2}$. Thus, the test is much faster to implement than when using the true $K_{1}$ and $K_{2}$. Nevertheless, both size and power appear to be only minimally affected; the rejection probability curve is similar to that in the left plot, for which the true $K_{1}$ and $K_{2}$ were used.

In the right plot of Fig. 3, the true $K_{1}$ and $K_{2}$ are used, but $\hat{\theta}_{0}$ and $\hat{\theta}_{1}$ are replaced by, respectively, the estimators $\check{\theta}_{0}$ and $\check{\theta}_{1}$ which ignore the constraints on $h^{2}$. Because these estimators are unconstrained, they are simpler to compute using off-the-shelf solvers; see Sec. 2.1. The rejection probability curve shows the test retains validity, as guaranteed by the arguments given before Theorem 1. However, the power of the test is clearly lower than that of using constrained estimators. This is intuitive as the unconstrained estimators effectively ignore information about the parameter $h^{2}$.

Fig.4 shows computing times for different implementations of the proposed split LRT statistic. We consider a setting where the $K_{m}$ are jointly diagonalizable as in (8). Specifically, we generated $K_{m}$, $m\in\{1,\dots,M\}$ as follows. First, we eigendecomposed an $n\times n$ autoregressive covariance matrix $A=(A_{ij})=(0.5^{|i-j|})=O\Lambda O^{\mathrm{\scriptscriptstyle T}}$. Then, $\Lambda_{m}$ was created by setting $\lfloor n/(M+1)\rfloor$ of its diagonal entries equal to the corresponding entries of $\Lambda$, and the remaining entries to zero. The indices for the non-zero entries for each $\Lambda_{m}$ were randomly sampled in such a way that they were different for each $\Lambda_{m}$, so that $\Lambda_{m}\Lambda_{\ell}=0$ for $\ell\neq m$, and $\sum_{m=1}^{M}\Lambda_{m}$ equals $\Lambda$ with $n-M\lfloor n/(M+1)\rfloor$ entries set to zero. Finally, we set $K_{m}=O\Lambda_{m}O^{\mathrm{\scriptscriptstyle T}}$, $m\in\{1,\dots,M\}$, so $K_{m}K_{\ell}=0$ for $\ell\neq m$. When there are $M=2$ components (left plot), data are generated with $h^{2}=(0,0.2)^{\mathrm{\scriptscriptstyle T}}$ and the null hypothesis is $h_{1}^{2}=0$. When $M=3$, data are generated with $h^{2}=(0,0,0.2)^{\mathrm{\scriptscriptstyle T}}$ and the null hypothesis is $h_{1}^{2}=h^{2}_{2}=0$.

The three methods in Fig. 4 include a naive method that implements the split likelihood ratio test as described in Sec. 2.1, without using the fact that the $K_{m}$ can be jointly diagonalized. The second considered method uses this fact only when calculating $\hat{\theta}_{0}$, emulating a setting where diagonalization is possible under the null hypothesis, but not in general, as in Sec. 3.1. The third method uses joint diagonalization both when computing $\hat{\theta}_{0}$ and $\hat{\theta}_{1}$. For small sample sizes, the methods are all fast and hence no large differences are seen in computing times. In contrast, when $n$ is in the thousands, using diagonalization leads to substantially faster computing. For example, with $M=2$ components and $n=3000$, the naive method takes on average about 900s, while the method that uses diagonalization under the null only takes about 270s on average, and the method which fully uses diagonalization takes about 80s on average. The plots are consistent with the facts that evaluating the likelihood and its derivatives takes $O(n)$ operations when using diagonalization and up to $O(n^{3})$ operations otherwise.

The simulation results, along with those of the data example in the next section, can be reproduced using code at https://github.com/koekvall/univ_vc_suppl.

5 Data Example

To illustrate the proposed methods, we apply them to a well-known dataset (Hicks and Turner,, 1999, Problem 6.18). The data consist of resistance measurements (in milliohms) obtained through a fully randomized design. Ten resistors and three operators were randomly selected. Each operator independently measured the resistance of each resistor twice, resulting in a total of 60 observations. Both resistor and operator effects are treated as random in the analysis. This can be motivated, for example, by thinking of the operators and resistors in the experiment as drawn from larger populations of potential operators and resistors, respectively. The random effects model correlations between measurement from the same operator, and measurements from the same resistor. Because every operator measures every resistor, the random effects are crossed. Specifically, suppose the $k$th measurement by operator $j$ on resistor $i$ satisfies

[TABLE]

where $\beta\in\mathbb{R}$ is the population mean, $U_{1i}\sim\mathrm{N}(0,\sigma^{2}_{1})$ is the resistor random effect, and $U_{2j}\sim\mathrm{N}(0,\sigma^{2}_{2})$ is the operator random effect. For simplicity, we estimate $\beta$ using the sample mean and apply our method to the centered responses $Y_{ijk}=\tilde{Y}_{ijk}-\sum_{i,j,k}\tilde{Y}_{ijk}/60$. That is, we fit the model

[TABLE]

where $Y$ is obtained by stacking the $Y_{ijk}$.

Maximum likelihood estimates of $\sigma^{2}$, computed with the lme4 package (Bates et al.,, 2015), are in Table 1. The table also includes the corresponding estimates of $h^{2}$ and $\tau^{2}$ (“total variance”). The estimate of $\sigma_{1}^{2}$, the variance of the resistor random effect, is zero. The estimate of the operator random effect $\sigma^{2}_{2}$ is about 50, which leads to an estimate of $h_{2}^{2}$ of 0.554. That is, it is estimated that about 55% of the total variability is due to the operator random effect.

Table 2 shows confidence intervals for the random effect standard deviations $\sigma_{1}$ and $\sigma_{2}$; confidence intervals for the variances can be obtained by squaring the endpoints. The table includes three standard methods, two of which are included by default in lme4, namely profile likelihood and bootstrap-based intervals; and the proposed method based on the split LRT. The Wald interval is not available by default in lme4 since it is known to be unreliable, but we nevertheless compute it for comparison. All methods give confidence intervals for $\sigma_{1}$ that include zero. The confidence intervals for $\sigma_{2}$ are relatively large, likely due to the small number of operators (three). However, only one of the intervals for $\sigma_{2}$, the Wald interval, includes zero. For both parameters, the proposed interval is substantially wider than the other three. However, it is also the only one that is known to be valid.

To compute confidence intervals with the proposed method, given the small number of observations, we use a $k$-fold method that can reduce the variability introduced by random data splits (Wasserman et al.,, 2020). We use $k=4$ folds to strike a balance between reducing randomness and retaining sufficient data within each fold. In each iteration, one fold is used to compute $\hat{\theta}_{1}$ and the three other folds are used to compute $\hat{\theta}_{0}$. Then the average of the four test-statistics are compared to the threshold $U/\alpha$ (c.f. 4). Figure 6 shows graphs of the resulting test-statistics for a range of $\sigma_{1}$ and $\sigma_{2}$. The confidence intervals contain the points where the corresponding graphs are below the critical value, which is drawn as a horizontal line. The particular realization of $U$ was 0.742, so the critical value is $0.742/0.05\approx 15$. For comparison, recall the non-randomized test has critical value $1/\alpha=20$.

We can find $p$-values for the null hypotheses $\sigma_{m}^{2}=0$, $m\in\{1,2\}$, by finding the smallest $\alpha\in(0,1)$ for which the tests would reject, with the convention that the $p$-value is one if no such $\alpha$ exists. For the randomized split LRT, the $p$-value for $\sigma_{1}=0$ is one while the $p$-value for $\sigma_{2}=0$ is zero. For comparison we used the exactRLRT function from the R package RLRsim, which provides likelihood ratio tests for boundary points based on simulation. The $p$-values for $\sigma_{1}=0$ and $\sigma_{2}=0$ are 0.452 and 0.00, respectively.

Because the critical value for the randomized test in general is different each time it is implemented, $p$-values and widths of confidence intervals also vary from implementation to implementation. In Figure 7 we examine the distribution of the widths of the confidence intervals in this data example. Because the non-randomized test compares to $1/\alpha$ rather than $U/\alpha$, its width is the greatest value on the horizontal axis where the density has support, which corresponds to a realization $U=1$. For example, the right plot in Fig. 7 indicates the non-randomized CI for $\sigma_{2}$ has width of about 90, as also seen in Fig. 6. The distribution of the width of the randomized CI has a mode around 80, which happens to be about the length we observed in Table 2. Using Monte Carlo, we found that the randomized CI widths were on average $81.1\%$ and $63.7\%$ of the non-randomized ones, for $\sigma_{1}$ and $\sigma_{2}$ respectively.

Finally, since the proofs of validity of the proposed methods require a correctly specified likelihood, we consider diagnostic plots. With $G=\operatorname{bdiag}(\hat{\sigma}^{2}_{1}I_{10},\hat{\sigma}^{2}_{2}I_{3})$, we can predict the random random effects with the best linear unbiased predictions (BLUPs) (Robinson,, 1991), i.e.,

[TABLE]

Thus, we can predict $\mathsf{E}(Y\mid U)=ZU$ by $Z\hat{U}$, and $E=Y-ZU$ by $\hat{E}=Y-Z\hat{U}$. If the predictions are accurate and the model is correct, we expect $\hat{E}$ to be distributed approximately as $E\sim\mathrm{N}(0,\sigma_{3}^{2}I_{n})$; the quantile-quantile plot in Fig.5 supports this approximation. The right plot in Fig. 5 shows $\hat{E}$ plotted against the predictions $Z\hat{U}$; there are only three levels because $\hat{\sigma}^{2}_{1}=0$. The variability in the $\hat{E}_{i}$ is perhaps slightly larger for the smallest value of $\hat{Y}_{i}$ compared to the other two, but overall the assumption that $E_{i}$ is independent of $ZU$, with $E_{i}\sim\mathrm{N}(0,\psi_{3})$, appears serviceable.

6 Conclusions

The proposed methods lead to valid inference on variance components even in settings where existing methods fail. In particular, to the best of our knowledge, it is the first method to be uniformly valid in settings where heritability, or a proportion of variation more generally, is near unity. The main drawback of the proposed methods is that they are conservative, but that is not a problem if a test rejects or if the confidence interval is narrow enough to be useful. Whether a confidence interval is narrow enough to be useful has to be assessed on a case-by-case basis by practitioners.

In many settings of interest the eigenvectors of the covariance matrix of the response vector do not depend on the variance components. Then, the algorithms provided here lead to much faster computing than naive ones. Further improvements may be possible with more research. For example, it is unclear how different methods for splitting the data affect the methods. It would be of interest to understand, for example, whether diagonalization is best performed before or after splitting. Similarly, in settings with crossed random effects, it may be preferable to, instead of splitting uniformly at random, balance the randomization so that an equal number of levels of a given factor is present in both splits. Finally, there may be room for improvements in both size and power by using estimators other than maximum likelihood. We focused on maximum likelihood estimators because they are common and have good large sample properties, but there are certainly many settings where we expect better estimation is possible. For instance, penalized likelihood-based estimators can be preferable in settings with many parameters.

Appendix A Technical details

Proof of Theorem 1.

The arguments are available in the literature (Ramdas and Manole,, 2023; Wasserman et al.,, 2020, Theorems 1.2 and 3, respectively). Pick an arbitrary $\theta^{*}\in\Theta_{0}$. By Markov’s inequality,

[TABLE]

so the first claim follows if $\mathsf{E}_{\theta^{*}}\left(T_{n}\right)\leq 1$. To show the latter, note that since $\theta^{*}\in\Theta_{0}$ and $\hat{\theta}_{0}$ is a maximizer over that set,

[TABLE]

Thus, it suffices to show that the last expectation equals one. To that end, condition on $Y_{(1)}$ and note that, since $\hat{\theta}_{1}$ is measurable with respect to the $\sigma$-algebra generated by $Y_{(1)}$, with probability one,

[TABLE]

Here, $f_{\theta}(y_{(0)}\mid Y_{(1)})$ denotes the density of the conditional distribution of $Y_{(0)}$ given $Y_{(1)}$ under $\theta$, evaluated at a fixed $y_{(0)}$ and random $Y_{(1)}$. Since the conditional expectation equals one, so does the unconditional one, and that proves validity of the split likelihood ratio test.

The proof of validity of the randomized split likelihood ratio test is similar except for the first step. Specifically, instead of using Markov’s inequality we note that by independence of $U$ and $T_{n}$,

[TABLE]

Thus, validity again follows from $\mathsf{E}_{\theta^{*}}\left(T_{n}\right)\leq 1$.

Finally, the claim about power is immediate from the fact that the event $\{\alpha T_{n}>1\}$ is a subset of the event $\{\alpha T_{n}>U\}$ since $U$ has support $(0,1)$. ∎

Proof of Theorem 2.

First note that, by symmetry, each $K_{m}$ has an eigendecomposition $K_{m}=O_{m}\Lambda_{m}O_{m}^{\mathrm{\scriptscriptstyle T}}$, where $\Lambda_{m}=\operatorname{diag}(\lambda_{m1},\dots,\lambda_{mn})$. Since the conclusion is obvious if $K_{m}=0$ for all $m$, suppose not. Pick arbitrary $K_{m}$ and $K_{\ell}$, $m\neq\ell$, such that $\lambda_{\ell k}\neq 0$ is an element of $\Lambda_{\ell}$ and $o_{\ell k}$ the corresponding column of $O_{\ell}$. Thus, $K_{\ell}o_{\ell k}=\lambda_{\ell k}o_{\ell k}$. Left-multiplying by $K_{m}$ and using $K_{m}K_{\ell}=0$ leads to

[TABLE]

so $o_{\ell k}$ is an eigenvector of $K_{m}$, corresponding to the eigenvalue zero. Thus, $o_{\ell k}$ is orthogonal to every column of $O_{m}$ corresponding to a nonzero eigenvalue of $K_{m}$. Take now every vector that, for some $m\in\{1,\dots,m\}$, is a column of $O_{m}$ corresponding to a nonzero eigenvalue of $K_{m}$, say

[TABLE]

This set is orthonormal. Indeed, we already showed that any two such vectors from different $O_{m}$ are orthogonal, and if they are from the same $O_{m}$ they are orthogonal by construction. Because $\mathcal{O}\subseteq\mathbb{R}^{n}$ is orthonormal, it contains at most $q\leq n$ vectors. Let these vectors be the leading $q$ columns of $O$ (in any order), and take the remaining $n-q$ vectors to be any orthonormal basis for the orthogonal complement of the span of $\mathcal{O}$. Note the last $n-q$ vectors are in the null space of every $K_{m}$, $m\in\{1,\dots,M\}$ by construction of $\mathcal{O}$. Now (8) holds upon possibly reordering elements in $\Lambda_{m}$, $m\in\{1,\dots,M\}$. ∎