On Testing Marginal versus Conditional Independence

TL;DR

This paper investigates the statistical challenge of distinguishing between marginal and conditional independence in Gaussian models, proposing a robust model selection method with uniform error guarantees even at the boundary of detectability.

Contribution

It introduces a novel approach to test and select between marginal and conditional independence models, addressing non-nested models and local alternatives with a new envelope distribution technique.

Findings

Envelope distributions provide uniform error control.

The proposed method achieves near-optimal power.

Model selection is robust even at indistinguishability boundaries.

Abstract

We consider testing marginal independence versus conditional independence in a trivariate Gaussian setting. The two models are non-nested and their intersection is a union of two marginal independences. We consider two sequences of such models, one from each type of independence, that are closest to each other in the Kullback-Leibler sense as they approach the intersection. They become indistinguishable if the signal strength, as measured by the product of two correlation parameters, decreases faster than the standard parametric rate. Under local alternatives at such rate, we show that the asymptotic distribution of the likelihood ratio depends on where and how the local alternatives approach the intersection. To deal with this non-uniformity, we study a class of "envelope" distributions by taking pointwise suprema over asymptotic cumulative distribution functions. We show that these envelope distributions are well-behaved and lead to model selection procedures with rate-free uniform error guarantees and near-optimal power. To control the error even when the two models are indistinguishable, rather than insist on a dichotomous choice, the proposed procedure will choose either or both models.

Tables

Table 1. Table 1: Envelope quantiles − F ¯ ρ − 1 ​ ( α ) superscript subscript ¯ 𝐹 𝜌 1 𝛼 -\bar{F}_{\rho}^{-1}(\alpha) , Monte Carlo standard errors ≤ 0.01 absent 0.01 \leq 0.01

$\rho$	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
$\alpha =0.05$	2.71	2.71	2.68	2.65	2.58	2.48	2.42	2.39	2.58	2.90	3.19
$\alpha =0.01$	5.41	5.41	5.40	5.34	5.27	5.21	5.11	5.05	5.02	5.40	5.97

Equations

$${\mathcal{M}}_0:X_1\text{scalebox1.07⊥⊥}{X}_2\text{versus}{\mathcal{M}}_1:X_1\text{scalebox1.07⊥⊥}{X}_2\mid X_3.$$

$${\mathcal{M}}_0:\{{\sigma }_{12}=0\},{\mathcal{M}}_1:\{{\sigma }_{12}{\sigma }_{33}={\sigma }_{13}{\sigma }_{23}\}$$

$${\mathcal{M}}_{\text{sing}}:\{{\sigma }_{12}={\sigma }_{13}={\sigma }_{23}=0\},$$

$${\ell }_n(\Sigma )=\frac{n}{2}(-\log |\Sigma |-\mathrm{Tr}(S_n{\Sigma }^{-1})),$$

$${\lambda }_n^{(i)}:=2\left(\underset{\Sigma \in \Theta }{\sup }{\ell }_n(\Sigma )-\underset{\Sigma \in {\Theta }_i}{\sup }{\ell }_n(\Sigma )\right)\ge 0$$

$$\begin{array}{rl}{\ell }_n^{\text{sat}}& :=\underset{\Sigma \in \Theta }{\sup }{\ell }_n(\Sigma )\\ & =-\frac{n}{2}\left(\log \left(s_{11}{s}_{22}{s}_{33}+2s_{12}{s}_{23}{s}_{13}-s_{11}{s}_{23}^2-s_{13}^2{s}_{22}-s_{12}^2{s}_{33}\right)+3\right).\end{array}$$

$$\begin{array}{rl}{\lambda }_n^{(0:1)}& :={\lambda }_n^{(1)}-{\lambda }_n^{(0)}\\ & =2\left(\underset{\Sigma \in {\Theta }_0}{\sup }{\ell }_n(\Sigma )-\underset{\Sigma \in {\Theta }_1}{\sup }{\ell }_n(\Sigma )\right)=2\left({\ell }_n({\hat{\Sigma }}_n^{(0)})-{\ell }_n({\hat{\Sigma }}_n^{(1)})\right),\end{array}$$

$$\begin{array}{rl}p(X_1,X_2,X_3)& =p(X_1)p(X_2)p(X_3\mid X_1,X_2)\\ & =\mathcal{N}(X_1;0,{\sigma }_{11})\mathcal{N}(X_2;0,{\sigma }_{22})\mathcal{N}(X_3;{\beta }_{32\cdot 1}{X}_1+{\beta }_{31\cdot 2}{X}_2,{\sigma }_{33\cdot 12}).\end{array}$$

$${\hat{\sigma }}_{11}^{(0)}=s_{11},{\hat{\sigma }}_{22}^{(0)}=s_{22},$$

$${\hat{\beta }}_{32\cdot 1}^{(0)}=\frac{s_{22}{s}_{13}-s_{12}{s}_{23}}{s_{11}{s}_{22}-s_{12}^2},{\hat{\beta }}_{31\cdot 2}^{(0)}=\frac{s_{11}{s}_{23}-s_{12}{s}_{13}}{s_{11}{s}_{22}-s_{12}^2},$$

$${\hat{\sigma }}_{33\cdot 12}^{(0)}=s_{33}-\frac{s_{22}{s}_{13}^2-2s_{12}{s}_{23}{s}_{13}+s_{11}{s}_{23}^2}{s_{11}{s}_{22}-s_{12}^2}.$$

$${\hat{\sigma }}_{13}^{(0)}=\frac{s_{11}\left(s_{22}{s}_{13}-s_{12}{s}_{23}\right)}{s_{11}{s}_{22}-s_{12}^2},{\hat{\sigma }}_{23}^{(0)}=\frac{s_{22}\left(s_{11}{s}_{23}-s_{12}{s}_{13}\right)}{s_{11}{s}_{22}-s_{12}^2}$$

$${\hat{\sigma }}_{33}^{(0)}=s_{33}-\frac{2s_{12}\left(s_{12}{s}_{13}-s_{11}{s}_{23}\right)\left(s_{12}{s}_{23}-s_{13}{s}_{22}\right)}{\left(s_{11}{s}_{22}-s_{12}^2\right){}^2}.$$

$${\ell }_n^{(0)}=-\frac{n}{2}\left[\log \left(s_{11}{s}_{22}\left(\frac{s_{22}{\hat{s}}_{13}^2-2s_{12}{s}_{23}{s}_{13}+s_{11}{s}_{23}^2}{s_{12}^2-s_{11}{s}_{22}}+s_{33}\right)\right)+3\right].$$

$${\hat{\sigma }}_{11}^{(1)}=s_{11},{\hat{\sigma }}_{22}^{(1)}=s_{22},{\hat{\sigma }}_{33}^{(1)}=s_{33},{\hat{\sigma }}_{13}^{(1)}=s_{13},{\hat{\sigma }}_{23}^{(1)}=s_{23},$$

$${\ell }_n^{(1)}=-\frac{n}{2}\left[\log \left(\frac{\left(s_{13}^2-s_{11}{s}_{33}\right)\left(s_{23}^2-s_{22}{s}_{33}\right)}{s_{33}}\right)+3\right].$$

$$\begin{gathered} {\lambda }_n^{(0:1)}=2({\ell }_n^{(0)}-{\ell }_n^{(1)})=n\log \left(\frac{\left(s_{13}^2-s_{11}{s}_{33}\right)\left(s_{23}^2-s_{22}{s}_{33}\right)}{s_{33}}\right)- \\ n\log \left(s_{11}{s}_{22}\left(\frac{s_{22}{s}_{13}^2-2s_{12}{s}_{23}{s}_{13}+s_{11}{s}_{23}^2}{s_{12}^2-s_{11}{s}_{22}}+s_{33}\right)\right). \end{gathered}$$

$$H^2(P_n^n,Q_n^n)=1-{\left\{1-H^2(P_n,Q_n)\right\}}^n.$$

$$H^2(P_n^n,Q_n^n)\le {\mathrm{d}}_{\text{TV}}(P_n^n,Q_n^n)\le H(P_n^n,Q_n^n){\left\{2-H^2(P_n^n,Q_n^n)\right\}}^{1/2}.$$

$$1-{\mathrm{d}}_{\text{TV}}(P_n^n,Q_n^n)\to \{\begin{array}{ll}0,& H^2(P_n,Q_n)=\omega (n^{-1})\\ 1,& H^2(P_n,Q_n)=o(n^{-1})\end{array}.$$

$$\begin{array}{rl}0<1-\{1-\exp (-2h){\}}^{1/2}& \le \underset{n\to \infty }{\mathrm{liminf}}\left\{1-{\mathrm{d}}_{\text{TV}}(P_n^n,Q_n^n)\right\}\\ & \le \underset{n\to \infty }{\mathrm{limsup}}\left\{1-{\mathrm{d}}_{\text{TV}}(P_n^n,Q_n^n)\right\}\le \exp (-h)<1.\end{array}$$

$$\begin{array}{rl}1-{\mathrm{d}}_{\text{TV}}(P_n^n,Q_n^n)\le 1-H^2(P_n^n,Q_n^n)& ={\left\{1-H^2(P_n,Q_n)\right\}}^n\\ & =\exp \left[n\log \left\{1-H^2(P_n,Q_n)\right\}\right].\end{array}$$

$$\begin{array}{rl}\underset{n\to \infty }{\mathrm{limsup}}\left\{1-{\mathrm{d}}_{\text{TV}}(P_n^n,Q_n^n)\right\}& \le \underset{n\to \infty }{\mathrm{limsup}}\exp \left[n\log \left\{1-H^2(P_n,Q_n)\right\}\right]\\ & =\{\begin{array}{ll}0,& H^2(P_n,Q_n)=\omega (n^{-1})\\ \exp (-h),& n{H}^2(P_n,Q_n)\to h>0\end{array}.\end{array}$$

$$\begin{array}{rl}\underset{n\to \infty }{\mathrm{liminf}}\left\{1-{\mathrm{d}}_{\text{TV}}(P_n^n,Q_n^n)\right\}& \ge \underset{n\to \infty }{\mathrm{liminf}}1-H(P_n^n,Q_n^n){\left\{2-H^2(P_n^n,Q_n^n)\right\}}^{1/2}\\ & =\{\begin{array}{ll}1,& H^2(P_n,Q_n)=o(n^{-1})\\ 1-\{1-\exp (-2h){\}}^{1/2},& n{H}^2(P_n,Q_n)\to h>0\end{array}.\end{array}$$

$$1-\exp (-h)\le \text{optimal asymptotic power}\le \alpha +\{1-\exp (-2h){\}}^{1/2}.$$

$$\begin{array}{rl}{\Sigma }_n^{(0)}& =\left(\begin{array}{ccc}{\sigma }_{11,n}& 0& {\rho }_{13,n}\sqrt{{\sigma }_{11,n}{\sigma }_{33,n}}\\ 0& {\sigma }_{22,n}& {\rho }_{23,n}\sqrt{{\sigma }_{11,n}{\sigma }_{33,n}}\\ {\rho }_{13,n}\sqrt{{\sigma }_{11,n}{\sigma }_{33,n}}& {\rho }_{23,n}\sqrt{{\sigma }_{11,n}{\sigma }_{33,n}}& {\sigma }_{33,n}\end{array}\right),\\ {\Sigma }_n^{(1)}& =\left(\begin{array}{ccc}{\sigma }_{11,n}& {\rho }_{13,n}{\rho }_{23,n}\sqrt{{\sigma }_{11,n}{\sigma }_{22,n}}& {\rho }_{13,n}\sqrt{{\sigma }_{11,n}{\sigma }_{33,n}}\\ {\rho }_{13,n}{\rho }_{23,n}\sqrt{{\sigma }_{11,n}{\sigma }_{22,n}}& {\sigma }_{22,n}& {\rho }_{23,n}\sqrt{{\sigma }_{22,n}{\sigma }_{33,n}}\\ {\rho }_{13,n}\sqrt{{\sigma }_{11,n}{\sigma }_{33,n}}& {\rho }_{23,n}\sqrt{{\sigma }_{22,n}{\sigma }_{33,n}}& {\sigma }_{33,n}\end{array}\right).\end{array}$$

$$H^2\left(P_{{\Sigma }_n^{(0)}},P_{{\Sigma }_n^{(1)}}\right)$$

$$=\{\begin{array}{ll}{\rho }_{13,n}^2{\rho }_{23,n}^2/8+O({\rho }_{13,n}^4+{\rho }_{23,n}^4),& {\rho }_{13,n},{\rho }_{23,n}\to 0\\ {\rho }_{23}^2(1-{\rho }_{23}^2{)}^{-1}{\rho }_{13,n}^2/8+O({\rho }_{13,n}^4),& {\rho }_{13,n}\to 0,{\rho }_{23,n}\to {\rho }_{23}\ne 0\\ {\rho }_{13}^2(1-{\rho }_{13}^2{)}^{-1}{\rho }_{23,n}^2/8+O({\rho }_{23,n}^4),& {\rho }_{23,n}\to 0,{\rho }_{13,n}\to {\rho }_{13}\ne 0\end{array}.$$

$$\{\begin{array}{ll}n=\omega \left(\frac{1}{{\rho }_{13}^2{\rho }_{23}^2}\right),& \text{for consistent model selection}\\ n\asymp \left(\frac{1}{{\rho }_{13}^2{\rho }_{23}^2}\right),& \text{for asymptotic total error ∈(0,1)}\end{array}.$$

$${\Sigma }^{\ast }=\left(\begin{array}{ccc}{\sigma }_{11}& 0& 0\\ 0& {\sigma }_{22}& \rho \sqrt{{\sigma }_{22}{\sigma }_{33}}\\ 0& \rho \sqrt{{\sigma }_{22}{\sigma }_{33}}& {\sigma }_{33}\end{array}\right)\in {\mathcal{M}}_0\cap {\mathcal{M}}_1\setminus {\mathcal{M}}_{\text{sing}}.$$

$${\Sigma }_n^{(0)}=\left(\begin{array}{ccc}{\sigma }_{11}& 0& \gamma \sqrt{{\sigma }_{11}{\sigma }_{33}}/\sqrt{n}\\ 0& {\sigma }_{22}& \rho \sqrt{{\sigma }_{22}{\sigma }_{33}}\\ \gamma \sqrt{{\sigma }_{11}{\sigma }_{33}}/\sqrt{n}& \rho \sqrt{{\sigma }_{22}{\sigma }_{33}}& {\sigma }_{33}\end{array}\right)\in {\mathcal{M}}_0\setminus {\mathcal{M}}_1,$$

Full text

model selection,

collider,

conditional independence,

likelihood ratio test,

confidence,

Gaussian graphical model,

journalname \usetkzobjall

and

1. 1 Introduction

2. 2 Maximum likelihood in a trivariate Gaussian model

3. 2.1 MLE within $\mathcal{M}_{0}$

4. 2.2 MLE within $\mathcal{M}_{1}$

5. 2.3 Likelihood ratio

6. 3 Optimal error

7. 4 Local asymptotics

8. 4.1 The weak-strong regime

9. 4.2 The weak-weak regime

10. 4.3 Limit experiments

11. 4.3.1 The weak-strong regime

12. 4.3.2 The weak-weak regime

13. 5 Envelope distributions

14. 5.1 The weak-weak regime

15. 5.2 The weak-strong regime

16. 6 Model selection procedures

17. 7 Simulations

18. 7.1 Local hypotheses

19. 7.2 Projected Wishart

20. 7.3 Conditional on covariates

21. 8 Real data example

22. 9 Discussion

1 Introduction

It is often of interest to test marginal or conditional independence for a set of random variables. For example, in the context of graphical modeling, the PC algorithm (Spirtes et al., 2000) for directed acyclic graph model selection determines the orientation of an unshielded triple $X-Z-Y$ based on whether the separating set of $X$ and $Y$ contains $Z$: if so, X\mathrel{\text{\scalebox{1.07}{\perp\mkern-10.0mu\perp}}}Y\mid Z and $Z$ is not a collider; if not, X\mathrel{\text{\scalebox{1.07}{\perp\mkern-10.0mu\perp}}}Y and the triple is oriented as $X\rightarrow Z\leftarrow Y$. The reader is referred to Dawid (1979); Lauritzen (1996); Koller et al. (2009) and Reichenbach (1956) for more discussion.

Here we consider the simplest case, namely testing X_{1}\mathrel{\text{\scalebox{1.07}{\perp\mkern-10.0mu\perp}}}X_{2} versus X_{1}\mathrel{\text{\scalebox{1.07}{\perp\mkern-10.0mu\perp}}}X_{2}\mid X_{3} in a trivariate Gaussian setting. For testing whether a specific marginal or conditional independence holds, it is common to use the correlation coefficient or partial correlation coefficient under Fisher’s $z$-transformation (Fisher, 1924) as the test statistic. Under independence, the transformed correlation coefficient is approximately distributed as a normal distribution with zero mean and variance determined by the sample size and the number of variables being conditioned on (Hotelling, 1953; Anderson, 1984). In this paper, however, we assume at least one type of independence holds (from prior knowledge or precursory inference) and we want to contrast the two types. To this end, we will use the likelihood ratio statistic, which often provides intuitively reasonable tests for composite hypotheses (Perlman and Wu, 1999), especially in terms of model selection.

Contributions

We briefly highlight our main contributions as follows. Firstly, we consider an important problem in non-nested model selection, which is in general less well-understood than the nested case. Secondly, we take an approach that is different from the usual Neyman-Pearson framework, in the sense that we treat the two models symmetrically and allow them to be both selected if the data does not significantly prefer one over the other. Thirdly, by introducing a new family of envelope distributions, we deal with non-uniform asymptotic laws of the likelihood ratio statistic. The model selection procedures we propose come with asymptotic guarantees that are applicable to all varieties of relations between the sample size and the signal strength; an assumption on the asymptotic rate is not required.

Notation

The following notation is used through the paper. $\mathbb{R}_{\text{PD}}^{n\times n}$ denotes $n\times n$ positive definite matrices. $\Theta$ denotes the parameter space and ${\mathcal{M}}$ denotes a model, which is subset of the parameter space. ${\mathcal{M}}_{1}\setminus{\mathcal{M}}_{2}$ denotes the set of parameters that belong to ${\mathcal{M}}_{1}$ but not belong to ${\mathcal{M}}_{2}$.

We use $P$ and $Q$ to denote measures, and similarly $P_{n}$ and $Q_{n}$ to denote sequences of measures. $\mu$ is reserved for the Lebesgue measure. Lower-case letters $p,q$ denote the densities of $P,Q$ with respect to $\mu$. $P_{n}^{n}$ denotes the $n$-sample product (tensorized) measure of $P_{n}$, namely the law of $X_{1},\dots,X_{n}\overset{\textsf{iid}}{\sim}P_{n}$. We write $P_{n}\overset{d}{\Rightarrow}P$ if $P_{n}$ converges (weakly) to $P$ in law. For $X_{n}(t)$ a stochastic process indexed by $t\in T$, we write $X_{n}\rightsquigarrow X$ if $X_{n}(t)$ converges weakly to $X(t)$.

For two sequences $\{a_{n}\}$ and $\{b_{n}\}$, we write $a_{n}=O(b_{n})$ if there exists a constant $c<\infty$ such that $a_{n}\leq cb_{n}$ for large enough $n$; $a_{n}=o(b_{n})$ and $b_{n}=\omega(a_{n})$ if $\lim_{n\rightarrow\infty}a_{n}/b_{n}=0$; $a_{n}\asymp b_{n}$ if $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$.

Also, we write $x\lesssim y$ if $x\leq cy$ for some constant $c>0$. We write $x\vee y=\max(x,y)$ and $x\wedge y=\min(x,y)$.

Setup

For $(X_{1},X_{2},X_{3})\sim\mathcal{N}\{0,\Sigma=(\sigma_{ij})\}$ with parameter space $\Theta$ being the set of $3\times 3$ real positive definite matrices $\mathbb{R}_{\text{PD}}^{3\times 3}$, we consider testing

[TABLE]

${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ are algebraic models (Drton and Sullivant, 2007) as represented by equality constraints

[TABLE]

imposed on $\Theta$. They are visualized in the correlation space (ignoring the variances) in Figure 1.

$\mathcal{M}_{0}$ and $\mathcal{M}_{1}$ are non-nested and they further intersect at the origin

[TABLE]

which is a singularity within ${\mathcal{M}}_{0}\cap{\mathcal{M}}_{1}$ that corresponds to diagonal covariances. At ${\mathcal{M}}_{\text{sing}}$ the likelihood ratio statistic is not regular in the sense that the tangent cones (linear approximations to the parameter space; see Bertsekas et al. (2003, Chap. 4.6)) of the two models coincide. As pointed out by Evans (2020), we will see that the equivalence of local geometry between the two models presents a challenge for model selection. It is also worth mentioning that, in the setting of nested model selection, the behavior of the likelihood ratio of testing $\mathcal{M}_{0}\cap\mathcal{M}_{1}$ against a saturated model, especially at the singularity, has been studied by Drton (2006); Drton and Sullivant (2007); Drton (2009).

Organization

The paper is organized as follows. In Section 2, we derive the maximum likelihood estimates under the two types of independence models, and obtain the loglikelihood ratio statistic in a closed form. In Section 3, we characterize the information-theoretic limit to distinguishing the two models, and outline two regimes on the boundary of distinguishability. Then in Section 4, we consider local alternative sequences in the two aforementioned regimes and establish the asymptotic distribution of the loglikelihood ratio. Section 4.3 provides a geometric perspective in terms of limit experiments. We then deal with non-uniformity issue of asymptotic distributions in Section 5 by introducing a family of envelope distributions. Next in Section 6 we propose model selection procedures with a uniform error guarantee. In Section 7, we compare the performance of several methods through simulation studies. We present a realistic example in Section 8 on inferring the American occupational structure. Finally some discussions are given in Section 9.

2 Maximum likelihood in a trivariate Gaussian model

The log-likelihood of a Gaussian graphical model under sample size $n$ (Lauritzen, 1996, Chap. 5) is

[TABLE]

where $S_{n}$ is the sample covariance computed with respect to mean zero (i.e., the scatter matrix divided by $n$). A model can be scored by its log-likelihood maximized within the model contrasted against the saturated model

[TABLE]

for $i=0,1$, which is the quantity considered in nested model selection. The saturated model attains maximal likelihood when $\Sigma=S_{n}$, yielding

[TABLE]

To contrast $\mathcal{M}_{0}$ and $\mathcal{M}_{1}$, we instead consider

[TABLE]

where $\hat{\Sigma}_{n}^{(0)}$ and $\hat{\Sigma}_{n}^{(1)}$ are MLEs within the two models. Intuitively, a positive value of $\lambda_{n}^{(0:1)}$ prefers $\mathcal{M}_{0}$, and a negative value prefers $\mathcal{M}_{1}$.

2.1 MLE within $\mathcal{M}_{0}$

By X_{1}\mathrel{\text{\scalebox{1.07}{\perp\mkern-10.0mu\perp}}}X_{2}, we can factorize the likelihood of $\mathcal{M}_{0}$ as

[TABLE]

where the parameters $\sigma_{11},\sigma_{22},\beta_{32\cdot 1},\beta_{31\cdot 2},\sigma_{33\cdot 12}$ are variation independent (Barndorff-Nielsen, 2014, Chap. 10.2). The MLEs for them are given by

[TABLE]

[TABLE]

and

[TABLE]

Mapping them back to the original parameters via relations $\beta_{32\cdot 1}=\sigma_{13}/\sigma_{11}$, $\beta_{31\cdot 2}=\sigma_{23}/\sigma_{22}$ and $\sigma_{33\cdot 12}=\sigma_{33}-\sigma_{13}^{2}/\sigma_{11}-\sigma_{23}^{2}/\sigma_{22}$, in addition to Eq. 7 we have the MLEs as

[TABLE]

and

[TABLE]

This derivation is essentially the same as executing the iterative conditional fitting algorithm of Chaudhuri et al. (2007) in the order of $X_{1},X_{2},X_{3}$. Plugging Eqs. 7, 8 and 9 into Eq. 4, we have the following closed-form expression of maximized log-likelihood of ${\mathcal{M}}_{0}$

[TABLE]

2.2 MLE within $\mathcal{M}_{1}$

The MLE within $\mathcal{M}_{1}$ in the covariance parametrization is simpler. By writing $\sigma_{12}=\sigma_{13}\sigma_{23}/\sigma_{33}$ and simplifying the score condition, we obtain

[TABLE]

all of which are their sample counterparts. Plugging into Eq. 4, we have

[TABLE]

2.3 Likelihood ratio

Finally, $\mathcal{M}_{0}$ and $\mathcal{M}_{1}$ are contrasted with

[TABLE]

3 Optimal error

We study the information-theoretic limit to distinguishing the two models. Specifically, consider two sequences of sampling distributions — one within ${\mathcal{M}}_{0}$ and the other within ${\mathcal{M}}_{1}$, as they approach the same limit in ${\mathcal{M}}_{0}\cap{\mathcal{M}}_{1}$. Let $P_{n}$ be the sequence in ${\mathcal{M}}_{0}$ under covariance $\Sigma_{n}^{(0)}\in{\mathcal{M}}_{0}\setminus{\mathcal{M}}_{1}$, and let $Q_{n}$ be the sequence in ${\mathcal{M}}_{1}$ under covariance $\Sigma_{n}^{(1)}\in{\mathcal{M}}_{1}\setminus{\mathcal{M}}_{0}$. Further, let $P_{n}^{n}$ and $Q_{n}^{n}$ be the product measures of $n$ independent copies of $P_{n}$ and $Q_{n}$ respectively.

The fundamental limit to distinguishing two distributions $P$ and $Q$ is characterized by their total variation distance $\operatorname*{\mathrm{d}_{\text{TV}}}(P,Q):=\sup_{A}\{P(A)-Q(A)\}$. We have the following classical result on testing two simple hypotheses, where the minimum total error is achieved by the likelihood ratio test.

Lemma 1** (Theorem 13.1.1 of Lehmann and Romano (2006)).**

For testing $H_{0}:X\sim P$ versus $H_{1}:X\sim Q$, the minimum sum of type-I and type-II errors is $1-\operatorname*{\mathrm{d}_{\text{TV}}}(P,Q)$.

The optimal error above does not permit a tractable formula. The analysis for a product measure is more tractable in terms of the Hellinger squared distance $H^{2}(P,Q):=(1/2)\int(p^{1/2}-q^{1/2})^{2}\mathop{}\!\mathrm{d}\mu$, for which it holds that

[TABLE]

The total variation is related to Hellinger by Le Cam’s inequality (Tsybakov, 2009, Lemma 2.3)

[TABLE]

Lemma 2** (see also Theorem 13.1.3 of Lehmann and Romano (2006)).**

It holds that

[TABLE]

And when $nH^{2}(P_{n},Q_{n})\rightarrow h>0$, it holds that

[TABLE]

Proof.

Using Eq. 14 and Eq. 15, we have

[TABLE]

It follows that

[TABLE]

Similarly, we also have

[TABLE]

The proof is finished by combining the previous two displays with the fact that $\liminf_{n}\left\{1-\operatorname*{\mathrm{d}_{\text{TV}}}(P_{n}^{n},Q_{n}^{n})\right\}\leq\limsup_{n}\left\{1-\operatorname*{\mathrm{d}_{\text{TV}}}(P_{n}^{n},Q_{n}^{n})\right\}$ and noting $\operatorname*{\mathrm{d}_{\text{TV}}}\in[0,1]$. ∎

Corollary 1**.**

Under $nH^{2}(P_{n},Q_{n})\rightarrow h>0$, the optimal power of an asymptotic $\alpha$-level procedure satisfies

[TABLE]

Proof.

This directly follows from Lemma 2 since $(1-\text{optimal power})+\text{type-I error}=1-\operatorname*{\mathrm{d}_{\text{TV}}}(P_{n}^{n},Q_{n}^{n})$ for type-I error asymptotically between 0 and $\alpha$, and then passing to the limit. ∎

By Lemma 2, the asymptotic error converges to zero (exponentially fast) if $P_{n}$ and $Q_{n}$ are separated by a distance that is decreasing more slowly than rate $n^{-1/2}$. For example, when $P_{n}=P$, $Q_{n}=Q$ are fixed distributions from which we observe $n$ independent samples, that is, when $P_{n}^{n}=P^{n}$ and $Q_{n}^{n}=Q^{n}$. The analysis above shows that the ability to differentiate $P_{n}$ and $Q_{n}$ based on $n$ samples depends on the distance between $P_{n}$ and $Q_{n}$. The consideration of $P_{n}$ and $Q_{n}$ as $n\rightarrow\infty$ is necessitated by the development of asymptotic results that are applicable in a specific analysis with a fixed $n$. In particular, here we want to investigate what happens when the sample size is small compared to the signal strength, or equivalently, when signal strength is weak under a given sample size. This is modeled by the regime that yields a non-trivial optimal error strictly between 0 and 1. By Lemma 2, we need to choose sequences $\Sigma_{n}^{(0)}$ and $\Sigma_{n}^{(1)}$ such that $H^{2}\left(P_{\Sigma_{n}^{(0)}},P_{\Sigma_{n}^{(1)}}\right)\asymp n^{-1}$. More specifically, we choose $\Sigma_{n}^{(1)}$ that is the most difficult to distinguish from $\Sigma_{n}^{(0)}$. That is, we choose $\Sigma_{n}^{(1)}$ to minimize $\operatorname*{\mathcal{D}_{\text{KL}}}\left(P_{\Sigma_{n}^{(0)}}\|P_{\Sigma_{n}^{(1)}}\right)$, i.e., $\Sigma_{n}^{(1)}$ is the MLE projection of $\Sigma_{n}^{(0)}$ in ${\mathcal{M}}_{1}$ by Eq. 11. The two sequences take the form of

[TABLE]

Both of them converge to $\Sigma^{\ast}\in{\mathcal{M}}_{0}\cap{\mathcal{M}}_{1}$ as $n\rightarrow\infty$. We assume the variances $\sigma_{ii,n}\rightarrow\sigma_{ii}>0$ for $i=1,2,3$. For $H^{2}\left(P_{\Sigma_{n}^{(0)}},P_{\Sigma_{n}^{(1)}}\right)\rightarrow 0$, it is necessary that either (or both) $\rho_{13,n}$ and $\rho_{23,n}$ converges to zero. The squared Hellinger distance is calculated as

[TABLE]

The calculation reveals that $H^{2}(P_{\Sigma_{n}^{(0)}},Q_{\Sigma_{n}^{(1)}})\asymp 1/n$ if and only if $\rho_{13,n}\rho_{23,n}\asymp n^{-1/2}$. This entails two distinct regimes.

The weak-strong regime

Between $\rho_{13,n}$ and $\rho_{23,n}$, one (the weak edge) converges to zero at $n^{-1/2}$ rate, and the other (the strong edge) converges to a non-zero limit $\rho\in(-1,1)$. The limiting model is on ${\mathcal{M}}_{0}\cap{\mathcal{M}}_{1}\setminus{\mathcal{M}}_{\text{sing}}$, namely one of the axes excluding the origin in Fig. 1.

The weak-weak regime

$\rho_{13,n},\rho_{23,n}\rightarrow 0$ and $\sqrt{n}\rho_{13,n}\rho_{23,n}\rightarrow\delta\neq 0$. The limiting model is on ${\mathcal{M}}_{\text{sing}}$, namely the origin in Fig. 1.

Remark 1**.**

The result can be rephrased as the sample size required to distinguish ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$. Consider distinguishing ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ in a Euclidean $m^{-1/2}$-neighborhood of $\Sigma^{\ast}\in{\mathcal{M}}_{0}\cap{\mathcal{M}}_{1}$ as $m\rightarrow\infty$. The sample size required is $m^{2}$ if $\Sigma^{\ast}\in{\mathcal{M}}_{\text{sing}}$, and $m$ if $\Sigma^{\ast}\notin{\mathcal{M}}_{\text{sing}}$. This phenomenon is described by Evans (2020) in terms of equivalence of local geometry. ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ are 1-equivalent at $\Sigma^{\ast}\in{\mathcal{M}}_{\text{sing}}$ in the sense that their tangent cones coincide; and they are 1-near-equivalent at $\Sigma^{\ast}\notin{\mathcal{M}}_{\text{sing}}$ in the sense that they have distinct tangent cones. See Evans (2020, Theorem 2.8).

Proposition 1**.**

In testing $\mathcal{M}_{0}$ versus $\mathcal{M}_{1}$, the sample complexity required is

[TABLE]

4 Local asymptotics

In this section, we analyze the asymptotic distribution of the log-likelihood ratio statistic $\lambda_{n}^{(0:1)}=2(\ell_{n}^{(0)}-\ell_{n}^{(1)})$ under the two regimes outlined earlier.

4.1 The weak-strong regime

Without loss of generality, we choose $\rho_{13,n}=\gamma/\sqrt{n}$ as the weak edge and $\rho_{23,n}\rightarrow\rho\neq 0$ as the strong edge. $\gamma\in\mathbb{R}$ characterizes the size of the local asymptotic, and is also referred to as a local parameter (van der Vaart, 2000, Chapter 9). We consider asymptotics under local sequences $\Sigma_{n}^{(0)}$ and $\Sigma_{n}^{(1)}$ approaching the limiting covariance

[TABLE]

We consider the following local alternatives of size $\gamma$ on the correlation scale. Again, $\Sigma_{n}^{(1)}$ is the KL-projection (i.e., MLE-projection) of $\Sigma_{n}^{(0)}$.

[TABLE]

[TABLE]

At the limit, both models are correct (intersection); see Figure 3. However, the sequence approaches the limit only on one of the models, and the size of violation of the other model is $|\gamma|n^{-1/2}$. To ensure positive definiteness, we require $|\rho|<1$.

Proposition 2**.**

Under local alternative $\Sigma_{n}^{(0)}$,

[TABLE]

and under local alternative $\Sigma_{n}^{(1)}$,

[TABLE]

where $Z_{1},Z_{2}$ are two independent standard normal variables.

We leave the proof to the next section, where we will present a geometric interpretation of the asymptotic distributions. Alternatively, the distribution can be derived by a change of measure with Le Cam’s third lemma; see van der Vaart (2000, Example 6.7).

Asymptotically the log-likelihood ratio statistic is distributed as a scaled difference of two independent non-central $\chi^{2}_{1}$ variables, with non-centralities scaled by $\gamma$ and weighted by $\rho$ differently, depending on the true model. Note that the distribution only depends on the absolute values of $\gamma$ and $\rho$. The asymptotic distributions under the two types of sequences (truths) are visualized in Fig. 2. We can see that the mean is positive under ${\mathcal{M}}_{0}\setminus{\mathcal{M}}_{1}$ and negative under ${\mathcal{M}}_{1}\setminus{\mathcal{M}}_{0}$. However, a pair of these distributions are not symmetric to each other in terms of shape. They are further separated apart (more easily distinguished) as $|\gamma|$ or $|\rho|$ becomes bigger, and only become identical (distributed as $\rho(Z_{1}^{2}-Z_{2}^{2})$) when $\gamma\rightarrow 0$.

Remark 2**.**

The models are locally asymptotically normal at $\Sigma^{\ast}\notin{\mathcal{M}}_{\text{sing}}$. By regularity, replacing the constant elements in Eqs. 21 and 22 with sequences in $n$ that converge to the corresponding limits does not alter the asymptotic distribution of $\lambda_{n}^{(0:1)}$.

4.2 The weak-weak regime

Now we study the asymptotic under $\rho_{13,n},\rho_{23,n}\rightarrow 0$ and $\sqrt{n}\rho_{13,n}\rho_{23,n}\rightarrow\delta$. The limiting covariance is $\Sigma^{\ast}=\operatorname{\mathsf{diag}}(\sigma_{11},\sigma_{22},\sigma_{33})\in{\mathcal{M}}_{\text{sing}}$, towards which we consider two local sequences

[TABLE]

and

[TABLE]

Proposition 3**.**

Given $\rho_{13,n}\rho_{23,n}=\delta n^{-1/2}+o(n^{-1/2})$ for $\delta\neq 0$ and $\rho_{13,n},\rho_{23,n}\rightarrow 0$. Under $\Sigma_{n}^{(i)}\in{\mathcal{M}}_{i}\setminus{\mathcal{M}}_{1-i}$ for $i=0,1$, we have

[TABLE]

The limit is a centered Gaussian shifted and then scaled by $\delta$. Plots for a few values of $\delta$ are given by Fig. 4.

Proof of Proposition 3.

By the coincidence of tangent cones of the two models in this regime, the distribution cannot be obtained from local asymptotic normality or contiguity relative to a tensorized static law. Instead, we perform a direct calculation. For convenience, we assume the form of (sub)sequences of $\rho_{13,n}$ and $\rho_{23,n}$ as

[TABLE]

for $a\in(0,1/2)$ and $\eta\tau=\delta$. We perform a manual change of measure by relating the law under $P_{\Sigma_{n}^{(i)}}$ to that under $P_{I}$, which is iid sampling of $\mathcal{N}(0,I)$. Under sample size $n$, suppose $\Omega_{n}$ is the sample covariance under $\mathcal{N}(0,I)$. Now suppose $S_{n}^{(i)}$ is the sample covariance under $P_{\Sigma_{n}^{(i)}}$ for $i=0,1$. Then it holds that

[TABLE]

for some $L_{n}^{(i)}$ such that $\Sigma_{n}^{(i)}=L_{n}^{(i)}L_{n}^{{}^{(i)}\intercal}$. Here we choose them as the Cholesky decompositions

[TABLE]

and

[TABLE]

By the central limit theorem, we have

[TABLE]

for $W$ a $3\times 3$ matrix of joint Gaussian variables whose covariance is determined by the Isserlis matrix. The asymptotic distribution of $\lambda_{n}^{(0:1)}$ can be obtained by substituting

[TABLE]

into the closed-form expression of Eq. 13 and simplifying. We have under $\Sigma_{n}^{(0)}$

[TABLE]

and under $\bm{\Sigma}_{n}^{(1)}$

[TABLE]

The result is immediate from $w_{12}\sim\mathcal{N}(0,1)$ and $\gamma\tau=\delta$. ∎

Remark 3**.**

The Gaussian asymptotic in Proposition 3 does not depend on how $\rho_{13,n}$ and $\rho_{23,n}$ approach zero individually. We verify it with simulations shown in Figure 5. We simulate under $n=10,000$ for $5,000$ replicates. We set $\rho_{13,n}=rn^{-a}$ and $\rho_{23,n}=tn^{-(1/2-a)}$ such that $\rho_{13,n}\rho_{23,n}=\delta n^{-1/2}$ for $\delta=rt$ under different values of $a$.

4.3 Limit experiments

We establish the equivalence of testing the two models local asymptotics to that of a limit experiment, which sheds light on the form of the asymptotic distribution. As we will see, the limit experiments are Gaussian location experiments and the problem is asymptotically equivalent to testing the location between two lines from a single normal observation. Further by weak convergence, $\lambda_{n}^{(0:1)}$ is asymptotically distributed as the likelihood ratio statistic arising from the limit experiment. The reader is referred to van der Vaart (2000, Chapter 7, 9 and 16) for more background.

4.3.1 The weak-strong regime

We characterize the limit experiment in the weak-strong regime.

Proposition 4**.**

The family of distributions $\{P_{\bm{\Sigma}^{\ast}+\operatorname{\mathsf{G}}h/\sqrt{n}}:h\in\mathbb{R}^{2}\}$ is locally asymptotically normal, where

[TABLE]

and $\operatorname{\mathsf{G}}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{3\times 3}$ is a linear operator

[TABLE]

Proof.

The Gaussian model $P_{\Sigma}$ is differentiable in quadratic mean at $\Sigma^{\ast}$. The result follows from van der Vaart (2000, 7.14 and 7.15). ∎

The limit experiment of a LAN (locally asymptotically normal) family is a normal location experiment.

Proposition 5**.**

The sequence of experiments indexed by the local parameter $h$ converges to the following normal location experiment

[TABLE]

where

[TABLE]

Proof.

$\{P_{\bm{\Sigma}^{\ast}+\operatorname{\mathsf{G}}h/\sqrt{n}}:h\in\mathbb{R}^{2}\}$ is LAN with non-singular Fisher information $I_{\Sigma}^{\ast}$, which is the conditional information matrix of $(\sigma_{12},\sigma_{13})$ under $P_{\Sigma}$ given $(\sigma_{11},\sigma_{22},\sigma_{33},\sigma_{23})$, corresponding to $(\mathsf{h}_{1},\mathsf{h}_{2})$. The result then follows from van der Vaart (2000, Corollary 9.5). ∎

The local sequences Eq. 21 and Eq. 22 can be identified as $\Sigma^{\ast}+\operatorname{\mathsf{G}}h/\sqrt{n}$ with $h$ taking value of

[TABLE]

respectively. Models ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ correspond to the set of $h_{0}$ and $h_{1}$ respectively as $\gamma$ varies in $\mathbb{R}$. That is, ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ are represented by local parameter spaces

[TABLE]

which consist of all limits of $\sqrt{n}\operatorname{\mathsf{G}}^{-1}(\Sigma_{n}^{(i)}-\Sigma^{\ast})$ for $i=0,1$ (see van der Vaart (2000, Chapter 7.4)). Note $H_{0}$ and $H_{1}$ are lines in $\mathbb{R}^{2}$ (affine) and they correspond to tangent cones from ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ at $\Sigma^{\ast}$ under Chernoff regularity; see also Drton (2009) and Geyer (1994).

Proposition 6**.**

Suppose $I_{\Sigma^{\ast}}^{-1}=LL^{\intercal}$. For $i=0,1$, under $P^{n}_{\Sigma^{\ast}+\operatorname{\mathsf{G}}h/\sqrt{n}}$ for $h=h_{i}$, it holds that $(-1)^{i}\lambda_{n}^{(0:1)}$ is asymptotically distributed as the likelihood ratio statistic of testing

[TABLE]

from a single observation $Z\sim\mathcal{N}(\mu=\bm{0},I_{2})$.

Proof.

Under $P^{n}_{\Sigma^{\ast}+\operatorname{\mathsf{G}}h/\sqrt{n}}$, by van der Vaart (2000, Theorem 16.7) $\lambda_{n}^{(0:1)}$ is asymptotically distributed as the log-likelihood ratio statistic for testing $H_{0}$ and $H_{1}$ based on a single sample from $\mathcal{N}(h,I_{\Sigma^{\ast}}^{-1})$. Note that the theorem still applies to our case even though $H_{0}$ and $H_{1}$ are non-nested, as its proof does not require the two models to be nested. That is, given $X\sim\mathcal{N}(m=\bm{0},I_{\Sigma^{\ast}}^{-1})$, we have

[TABLE]

which is equivalent to testing $m\in H_{0}-h$ versus $m\in H_{1}-h$ from $X$. Given $I_{\Sigma^{\ast}}^{-1}=LL^{\intercal}$, by rewriting $X=_{d}LZ$ for $Z\sim\mathcal{N}(\mu=\bm{0},I_{2})$, the testing problem is mapped to that from $Z$ by $L^{-1}$. Hence, this is further equivalent to testing

[TABLE]

from $Z$. Note that $H_{i}-h_{i}=H_{i}$ since $H_{i}$ is affine. ∎

Now we derive limit experiments based on Proposition 6. The Cholesky decomposition gives

[TABLE]

We have, when $h=h_{0}$

[TABLE]

and when $h=h_{1}$

[TABLE]

They are visualized in Figure 6. The limit experiments Eq. 38 and Eq. 39 are of the same type as they are both characterized by an angle and an intercept. The two have the same angle $\theta=\arcsin\rho$ and their intercepts are related by a factor of $1/\sqrt{1-\rho^{2}}$.

Now we prove the form of the asymptotic distributions in Proposition 2 from the limit experiment.

Proof of Proposition 2.

Since the limit experiments are of the same type, we only derive for local alternatives $\Sigma_{n}^{(1)}\in{\mathcal{M}}_{1}\setminus{\mathcal{M}}_{0}$. We set the coordinate system as in Fig. 7, where the bisector of angle $\angle BCA=\theta=\arcsin\rho$ is the $y$-axis. The standard Gaussian vector centered at $B$ is represented as $Z=(x,y)=(\gamma\sin(\theta/2)-Z_{1},Z_{2})$. By the limit experiment, we have $\lambda_{n}^{(0:1)}\overset{d}{\Rightarrow}d_{1}^{2}-d_{2}^{2}$. ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ are respectively represented by lines $y=\pm kx+a$ for $k=\cot(\theta/2)$ and $a=-\gamma\cos(\theta/2)$. We have

[TABLE]

where we used

[TABLE]

By a change of variables $(Z_{1},Z_{2})=_{d}((U+V)/\sqrt{2},(U-V)/\sqrt{2})$ for another pair of independent standard normals and using the fact

[TABLE]

upon simplifying we have

[TABLE]

∎

4.3.2 The weak-weak regime

The Gaussian limit in Proposition 3 shows that the limit experiment of the weak-weak regime is testing the location of a univariate normal between two points; see the last panel of Fig. 6.

Corollary 2**.**

Testing ${\mathcal{M}}_{0}$ versus ${\mathcal{M}}_{1}$ under $\sqrt{n}\rho_{12,n}\rho_{13,n}\rightarrow\delta$ for $\delta\neq 0$ with $\rho_{12,n},\rho_{13,n}\rightarrow 0$ is asymptotically equivalent to testing $H_{0}:\mu=0$ versus $H_{1}:\mu=\delta$ from a single observation $Z\sim\mathcal{N}(\mu,1)$.

5 Envelope distributions

Though it may at first appear otherwise, the asymptotic distributions as obtained in Proposition 2 and Proposition 3 are not directly applicable to forming decision rules. This is due to the non-uniformity of the asymptotics.

Firstly, the asymptotic depends on the regime: weak-strong versus weak-weak, namely where the local sequence converges to. And the law is discontinuous between the two regimes. That is, the law in the weak-strong regime (scaled difference of noncentral chi-squares) does not converge to that of the weak-weak regime (Gaussian) as $\rho\rightarrow 0$. Furthermore, a procedure that firstly estimates the regime and then uses the corresponding distribution to form decision boundary, is susceptible to irregularity issues. Additionally, it is difficult to judge if an edge is weak based on whether its confidence interval contains zero without further assumptions, as illustrated by the following example.

Example 1**.**

Suppose $X_{i}\overset{\textsf{iid}}{\sim}\mathcal{N}(\gamma/\sqrt{n},\sigma^{2})$ for $i=1,\cdots,n$. The usual $(1-\alpha)$-level confidence interval for the mean of $X$ is $\bar{X}_{n}\pm z_{\alpha/2}\hat{\sigma}_{n}/\sqrt{n}$. The probability that it contains zero is

[TABLE]

for $\gamma\neq 0$ and $Z\sim\mathcal{N}(0,1)$. A large enough $\gamma$ can be chosen to make this probability arbitrarily small.

Secondly, given the regime, the distribution depends on the value of a local parameter ($\gamma$ for strong-weak and $\delta$ for weak-weak), which determines how the local sequence converges. Due to the $\sqrt{n}$ factor, the standard error for its estimator does not vanish and in general the local parameter cannot be consistently estimated. The reader is referred to Berger and Boos (1994); Andrews (2001) for discussions in the literature on the treatment of asymptotic distributions involving nuisance parameters that are not point-identified. Here we take a different approach, presented as follows.

The non-uniformity of asymptotic distributions motivates us to seek a procedure that circumvents the inference on the regime and the local parameter. In this section, we study the “extremal” distributions arising from the asymptotic distributions as the local parameter varies in $\mathbb{R}$.

Definition 1**.**

Given a family of distribution functions $\{F_{h}:h\in\mathcal{H}\}$ on $\mathbb{R}$, define

[TABLE]

and

[TABLE]

We call $\bar{F}$ the envelope distribution of $\{F_{h}:h\in\mathcal{H}\}$ if $\bar{F}$ is a valid distribution function.

Lemma 3**.**

$\bar{F}^{\ast}(x)$* is left-continuous if every $F_{h}(x)$ for $h\in\mathcal{H}$ is continuous.*

Proof.

Fix any $x$ and $\delta>0$, for $\epsilon>0$ we have $|\bar{F}^{\ast}(x)-\bar{F}^{\ast}(x-\epsilon)|=\sup_{h}F_{h}(x)-\sup_{h}F_{h}(x-\epsilon)$. By definition of supremum, there exists $h^{\prime}\in\mathcal{H}$ such that $F_{h^{\prime}}(x)\geq\sup_{h}F_{h}(x)-\delta/2$. Hence, $|\bar{F}^{\ast}(x)-\bar{F}^{\ast}(x-\epsilon)|\leq\delta/2+F_{h^{\prime}}(x)-F_{h^{\prime}}(x-\epsilon)$. By continuity of $F_{h^{\prime}}$, choosing $\epsilon>0$ such that $F_{h^{\prime}}(x)-F_{h^{\prime}}(x-\epsilon)\leq\delta/2$ shows that $\bar{F}^{\ast}(x)$ is left-continuous. ∎

Lemma 4**.**

If $\bar{F}^{\ast}(x)\rightarrow 0$ as $x\rightarrow-\infty$, then $\bar{F}(x)$ is a valid distribution function.

Proof.

Given any $x\leq x^{\prime}$, $\sup_{h}F_{h}(x)\leq\sup_{h}F_{h}(x^{\prime})$ by monotonicity of every $F_{h}$. Since $\bar{F}^{\ast}$ is non-decreasing, by Folland (1999, Theorem 3.23), the set of points at which $\bar{F}^{\ast}$ is discontinuous is countable. By redefining the function value at these points to be their right limits, $\bar{F}$ is right continuous. Also, $\bar{F}(x)\geq\bar{F}^{\ast}(x)\rightarrow 1$ as $x\rightarrow+\infty$ since every $F_{h}(x)\rightarrow 1$. Finally, as $x\rightarrow-\infty$ if $\bar{F}^{\ast}(x)\rightarrow 0$ , then $\bar{F}(x)\rightarrow 0$. $\bar{F}$ is a distribution function. ∎

5.1 The weak-weak regime

Proposition 7**.**

Let $G_{\delta}=\{\mathcal{N}(\delta^{2},(2\delta)^{2}):\delta\in\mathbb{R}\}$ be the asymptotic distributions for the weak-weak regime under ${\mathcal{M}}_{0}\setminus{\mathcal{M}}_{1}$. The envelope of $\{G_{\delta}\}$ is an equal-probability mixture of $(-\chi_{1}^{2})$ and a point mass at zero, namely

[TABLE]

The corresponding envelope under ${\mathcal{M}}_{1}\setminus{\mathcal{M}}_{0}$ is distributed as its negation.

Proof.

It suffices to consider $\delta\geq 0$. Given any $x<0$,

[TABLE]

where $\delta^{\ast}=\sqrt{-x}$ is the maximizer; Given any $x\geq 0$, $\delta=0$ maximizes the probability to one. Hence, the envelope CDF is

[TABLE]

from which it follows that

[TABLE]

The envelope for ${\mathcal{M}}_{1}\setminus{\mathcal{M}}_{0}$ follows from symmetry. ∎

Since when ${\mathcal{M}}_{0}$ is true, the region for rejecting ${\mathcal{M}}_{0}$ should take the form $(-\infty,r)$ for some $r<0$, only the negative part of $\bar{G}$ is relevant for decision making. It follows from Proposition 7 that the negative part of $\bar{G}$ is distributed as $\chi_{1}^{2}$. The formation of the envelope is visualized in Fig. 8, which aligns with the behavior observed in Fig. 4, where as $\delta$ grows, the quantiles for $\alpha=0.05$ first moves outward for $\delta\in(0.5,1.64)$ and then moves inward for $\delta\in(1.64,\infty)$.

5.2 The weak-strong regime

Now we study the envelope distributions under the weak-strong regime. We first observe that the envelope distributions, if they exist, must be symmetric for Eq. 23 and Eq. 24, in the sense that they are distributed as the negation of each other. The symmetry holds because the two local parameters are related by a factor of $1/\sqrt{1-\rho^{2}}$ (see Fig. 6), and hence the suprema are taken over the same set of laws up to a difference in the sign. Fix $\rho$, let $\{F_{\rho,\gamma}:\gamma\in\mathbb{R}\}$ be the family of asymptotic distributions in the weak-strong regime under ${\mathcal{M}}_{0}\setminus{\mathcal{M}}_{1}$ as given in Eq. 23. Let $\bar{F}_{\rho}$ be its envelope distribution function.

Proposition 8**.**

$\bar{F}_{\rho}$* is a valid distribution function for $|\rho|\in(0,1]$.*

Proof.

Note since $\bar{F}_{\rho}=\bar{F}_{-\rho}$, it suffices to consider $\rho\in(0,1]$. First consider $\varphi_{\rho,\gamma}(x)$, the density function for $X^{2}-Y^{2}$ with $X\sim\mathcal{N}\left(\mu_{1}=\gamma\sqrt{\frac{1-\rho}{2}},1\right)$ and $Y\sim\mathcal{N}\left(\mu_{2}=\gamma\sqrt{\frac{1+\rho}{2}},1\right)$ for $\gamma\in\mathbb{R}$ and $\rho\in(0,1]$. Since $p(X^{2}-Y^{2}=v^{2},Y^{2}=t)=p(Y^{2}=t)p(X^{2}=t+v^{2})$, the density $\varphi_{\rho,\gamma}$ has the following integral representation from marginalization

[TABLE]

Recall that $\lesssim$ allows for a positive multiplicative constant. Using $\cosh(x)<\exp(x)$ for $x>0$, we have

[TABLE]

We note that

[TABLE]

by completing the square in $\gamma$. It then follows that

[TABLE]

where we used $2\sqrt{t(t+v^{2})}\leq 2t+v^{2}$ in the third line. $K_{\nu}(\cdot)$ is the modified Bessel function of the second kind, and has the following asymptotic expansion for $z>0$ (Abramowitz and Stegun, 1972, Page 378)

[TABLE]

Hence for large $v^{2}$, we have

[TABLE]

Recall that $\{F_{\rho,\gamma}:\gamma\in\mathbb{R}\}$ is the family of distributions for the RHS of Eq. 23. With $\gamma^{\prime}=\sqrt{1-\rho^{2}}\gamma$,

[TABLE]

It follows that the density function

[TABLE]

where the exponent $\frac{1+\rho-\sqrt{1-\rho^{2}}}{4\rho}\in(1/4,1/2]$. By Definition 1, we have

[TABLE]

and hence $\bar{F}^{\ast}_{\rho}(-v^{2})\rightarrow 0$ as $v\rightarrow\infty$. By Lemma 4, $\bar{F}_{\rho}$ is a distribution function for every $\rho\in(0,1]$. ∎

The following proposition shows that $F_{\rho,\gamma=0}$ constitutes the envelope for the positive part of $\bar{F}_{\rho}$.

Proposition 9**.**

The positive part of $\bar{F}_{\rho}$ for $|\rho|\in(0,1]$ is distributed as the positive part of $\rho(Z_{1}^{2}-Z_{2}^{2})$ for $Z_{1},Z_{2}\overset{\textsf{iid}}{\sim}\mathcal{N}(0,1)$.

Proof.

Fix $\rho\in(0,1]$ and $v^{2}\geq 0$, with $\gamma^{\prime}=\gamma\sqrt{1-\rho^{2}}$ it follows from Proposition 2 that

[TABLE]

where $\mu_{1}=\gamma\sqrt{(1+\rho)/2}$, $\mu_{2}=\gamma\sqrt{(1-\rho)/2}$. Since $F_{\rho,\gamma^{\prime}}$ is symmetric in $\gamma$, we show $\gamma=0$ maximizes $F_{\rho,\gamma^{\prime}}(v^{2})$ by showing that the probability on the RHS of Eq. 43 increases in $\gamma\in(0,\infty)$. The probability can be interpreted as the standard Gaussian measure of the hyperbolic set $\{(x,y):x^{2}-y^{2}\geq v^{2}/\rho\}$ with the Gaussian centered at $G=(\mu_{1},\mu_{2})=\gamma(\sqrt{(1+\rho)/2},\sqrt{(1-\rho)/2})$. This is visualized in Fig. 9, where $\gamma=\overline{OG}$, $\tan\phi=\sqrt{(1-\rho)/(1+\rho)}$ and the hyperbolic set consists of the area inside the two branches. As $\gamma$ increases from zero, the center moves away from the origin along the $V$ line. Let $U$ be the line perpendicular to $V$. The Gaussian measure has two independent standard normal projections $(U,V)$, which is a rotation of $(Z_{1},Z_{2})$. Now we show that for every $v>0$, by conditioning on $|V|=v$, the conditional probability of $U$ in the appropriate “section” of the hyperbolic set, denoted by probability $q(v)$, increases with $\gamma$.

Let $[A,B]$ and $[C,D]$ be the line segments that $V=-v$ and $V=v$ intersect the hyperbola respectively. By independence of $U$ and $V$, we have $q(v)=\Pr(U\in[A,B])+\Pr(U\in[C,D])$. Let $\hat{v}$ and $\bar{v}$ be the distance from $G$ to the tangent to the left and right branch of the hyperbola respectively, parallel to line $U$; see Fig. 9. There are three cases. (i) When $v\leq\bar{v}$ (the first panel of Fig. 9), as $\gamma$ increases, both $[A,B]$ and $[C,D]$ become bigger, and thus $q(v)$ increases. (ii) When $\bar{v}<v\leq\hat{v}$, $[A,B]$ is empty but $[C,D]$ becomes bigger, so $q(v)$ increases. (iii) When $v>\hat{v}$, as $\gamma$ increases (the second panel of Fig. 9), $[C,D]$ increases but $[A,B]$ decreases. Let $[E,F]$ be the segment symmetric to $[A,B]$ about the origin. We observe that, as $\gamma$ increases by an infinitesimal $\Delta\gamma$, the amount that $\Pr(U\in[A,B])$ decreases equals the amount that $\Pr(U\in[E,F])$ increases, which is smaller than the amount that $\Pr(U\in[C,D])$ increases. Hence, $q(v)$ still increases.

By the monotonicity for every value of $|V|$, we conclude that the total probability on the RHS of Eq. 43 increases in $\gamma$. Hence, $F_{\gamma,\rho}(v^{2})$ is maximized at $\gamma=0$ for every $v$, namely $\bar{F}_{\rho}=F_{\rho,\gamma=0}$. It follows that for $X\sim\bar{F}_{\gamma}$, $(X)_{+}=_{d}\rho(Z_{1}^{2}-Z_{2}^{2})_{+}$ for two independent standard normal variables $Z_{1},Z_{2}$. ∎

Corollary 3**.**

$\bar{F}_{\rho}(0)\equiv 1/2$.

Unfortunately, we do not have an analytic form of the distribution for the negative part of $\bar{F}_{\rho}$, which is the part relevant for decision making, except for $\rho\rightarrow 0$ and $\rho=1$.

Proposition 10** (Bessel envelope).**

$\bar{F}_{\rho=1}=_{d}Z_{1}^{2}-Z_{2}^{2}$* for $Z_{1},Z_{2}\overset{\textsf{iid}}{\sim}\mathcal{N}(0,1)$.*

Proof.

Under $\rho=1$, the CDF is

[TABLE]

Since the conditional probability is non-negative, it suffices to show that given any $x\in\mathbb{R}$, $\gamma=0$ maximizes $\Pr\left((Z_{1}+\gamma)^{2}\leq x+z_{2}^{2}\mid Z_{2}=z_{2}\right)=\Pr\left((Z_{1}+\gamma)^{2}\leq x+z_{2}^{2}\right)$ for all $z_{2}\in\mathbb{R}$. When $x+z_{2}^{2}\leq 0$, the conditional probability is zero and $\gamma=0$ is trivially a maximizer. When $x+z_{2}^{2}>0$, then $\Pr\left((Z_{1}+\gamma)^{2}\leq x+z_{2}^{2}\right)=\Phi(\sqrt{x+z_{2}^{2}}-\gamma)-\Phi(-\sqrt{x+z_{2}^{2}}-\gamma)$. Setting the derivative with respect to $\gamma$ to zero requires $\phi(-\sqrt{x+z_{2}^{2}}-\gamma)=\phi(\sqrt{x+z_{2}^{2}}-\gamma)$, to which $\gamma=0$ is the unique solution. Therefore, $\gamma=0$ is the unique maximizer of $F_{\gamma}(x)$ for all $x$. ∎

The distribution in Proposition 10 is a difference between two independent $\chi_{1}^{2}$ variables. The density, as plotted in Fig. 10, is

[TABLE]

where $K_{0}$ is a modified Bessel function of the second kind. It is referred to as a $K$-form Bessel distribution in the literature; see Johnson et al. (1995, Chapter 4.4), Bhattacharyya (1942) and Simon (2007, Page 25).

Proposition 11** (Continuity of envelope).**

$\bar{F}_{\rho}\overset{d}{\Rightarrow}\bar{G}$* as $\rho\rightarrow 0$, where $\bar{G}$ is the envelope distribution for the weak-weak regime given in Proposition 7.*

Proof.

Firstly, we note that $\bar{F}_{\rho}(0)=\bar{G}(0)=1/2$ and by Proposition 9 the non-negative part of $\bar{F}_{\rho}$ also converges to that of $\bar{G}$ as $\rho\rightarrow 0$, namely a point mass at zero. It remains to be shown that the negative part of $\bar{F}_{\rho}$ converges in law to the negative part of $\bar{G}$. It suffices to show for any $x\leq 0$

[TABLE]

as $\rho\rightarrow 0$. Given $\rho>0$, the maximized probability can be rewritten as

[TABLE]

where we define

[TABLE]

for $\rho\in[0,1)$ and $\delta\in\mathbb{R}$. Note that $\sup_{\delta}\Pr(X_{0}(\delta)\leq x)=\sup_{\delta}\Pr(\delta^{2}+2\delta Z\leq x)=\Pr(-Z^{2}\leq x)/2$ for $Z\sim\mathcal{N}(0,1)$ by Proposition 7. We are left to show $\sup_{\delta}\Pr\left(X_{\rho}(\delta)\leq x\right)\rightarrow\sup_{\delta}\Pr\left(X_{0}(\delta)\leq x\right)$ as $\rho\rightarrow 0$. Choose $x<M<\infty$ and define

[TABLE]

We observe that

[TABLE]

where the last step follows from weak convergence $\{Y_{\rho}(\delta):\delta\in\mathbb{R}\}\rightsquigarrow\{Y_{0}(\delta):\delta\in\mathbb{R}\}$ in $\ell^{\infty}(\mathbb{R})$ as $\rho\rightarrow 0$ for a bounded stochastic process; see van der Vaart (2000, Chap. 18). ∎

Perhaps surprisingly, Proposition 11 shows that the asymptotic envelope is continuous between the two regimes, which bridges the discontinuity of the asymptotic distributions of $\lambda_{n}^{(0:1)}$ as presented in Propositions 2 and 3. Therefore, taking the envelope resolves the non-uniformity issue in terms of both the regime and the local parameter. Now with this we extend the definition of the envelope $\bar{F}_{\rho}$ to $\rho\in[0,1]$ by writing $\bar{F}_{\rho=0}=\bar{G}$.

Figure 11 showcases two envelopes. In the absence of an analytic form for $\rho\in(0,1)$, the envelopes can be numerically simulated by taking the supremum over a grid of values for $\gamma$. We observe from simulations that there exists $\gamma^{\ast}(x,\rho)\in(0,\infty)$ such that $\pm\gamma^{\ast}$ uniquely maximizes $F_{\gamma,\rho}(x)$.

Finally, we conclude this section by noting the following envelope of envelopes. See Figure 12 for an illustration. This result will be used in the next section to form simple decision rules based on the Bessel distribution.

Proposition 12** (Envelope of envelopes).**

The negative part of the envelope of $\{\bar{F}_{\rho}:\rho\in[0,1]\}$ is distributed as the negative part of $\bar{F}_{\rho=1}$ (Bessel).

6 Model selection procedures

Since we are selecting between two non-nested models, we want to refrain from choosing one of them as the default (the null hypothesis). By treating ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ symmetrically, however, a procedure that takes output value in $\{{\mathcal{M}}_{0},{\mathcal{M}}_{1}\}$ cannot simultaneously control both types of error under a given tolerance. It can be seen from Figs. 2 and 5 that there are cases where the asymptotic distributions of $\lambda_{n}^{(0:1)}$ under $P_{\Sigma_{n}^{(0)}}$ and $P_{\Sigma_{n}^{(1)}}$ significantly overlap. In such cases, insisting on a dichotomous choice will inevitably result in a high probability of error under at least one model.

To deal with this possible indistinguishability, we opt for a procedure with three options: if two models can be sufficiently distinguished, it selects one of them; otherwise it refrains from commitment by selecting both models, formally denoted as the union ${\mathcal{M}}_{0}\cup{\mathcal{M}}_{1}$. It is worth stressing that we always assume at least one of the two models is true. By such a design, when the procedure does not output the union, we are ensured that the probability of choosing the wrong model is small, being controlled below a given tolerance $\alpha$. In contrast, in the usual hypothesis testing framework where supposedly ${\mathcal{M}}_{0}$ is the null and ${\mathcal{M}}_{1}$ is the alternative, one typically cannot simultaneously control both type-I and type-II error. In other words, our procedure selects model with “confidence”. Recently the same notion has been investigated by Lei (2014) in a classification setting; Robins et al. (2003) also allows a test to make no decision when faced with ambiguity. We formalize the concept as follows.

Suppose $(X_{1},X_{2},X_{3})\in\mathbb{R}^{n\times 3}$ consists of $n$ independent samples from $\mathcal{N}(\bm{0},\Sigma_{n})$, where $\Sigma_{n}\in{\mathcal{M}}_{0}\cup{\mathcal{M}}_{1}$ is allowed to change with $n$ and $\Sigma_{n}\rightarrow\Sigma^{\ast}$. The sequence $\Sigma_{n}$ models signal strength relative to the sample size. We consider a deterministic decision rule

[TABLE]

where sample covariance $S_{n}$ is the sufficient statistic. For a given sequence of $\Sigma_{n}\in{\mathcal{M}}_{i}\setminus{\mathcal{M}}_{1-i}$ with $\Sigma_{n}\rightarrow\Sigma^{\ast}\in{\mathcal{M}}_{i}\cup{\mathcal{M}}_{1-i}$, we define the asymptotic (type-I) error of $\phi_{n}$ as the large-sample probability of rejecting the true model, i.e.,

[TABLE]

where the probability is taken under $P_{\Sigma_{n}}^{n}$. Similarly, the asymptotic power is defined as

[TABLE]

We say that the error is uniformly controlled below a given size $\alpha>0$, if

[TABLE]

where for $i=0,1$ the supremum for $(\Sigma_{n}^{(i)})$ is taken over all converging sequences of $\Sigma_{n}^{(i)}$ within ${\mathcal{M}}_{i}\setminus{\mathcal{M}}_{1-i}$ (the limit can be in either ${\mathcal{M}}_{i}\setminus{\mathcal{M}}_{1-i}$ or ${\mathcal{M}}_{i}\cap{\mathcal{M}}_{1-i}$). In general, the power $p_{\text{pow}}((\Sigma_{n}))$ depends on the sequence considered and we do not seek power optimality or guarantee in a uniform sense. In the next section, we will compare the power of several proposed procedures to the theoretical optimal for $\Sigma_{n}$ considered in the two regimes of local asymptotics.

By construction, using the $\alpha$-quantile of the envelope as the decision boundary achieves uniform error control. Based on the envelope of envelopes, a simple uniform rule is

[TABLE]

To gain more power, since $\bar{F}_{\rho}$ is continuous in $\rho$ and $\rho$ can be consistently estimated (recall that $\rho=\rho_{\text{strong}}$ in the weak-strong regime, and $\rho=0$ in the weak-weak regime), an adaptive rule can be formed as

[TABLE]

where $\hat{\rho}_{n}=|\hat{\rho}_{13,n}|\vee|\hat{\rho}_{23,n}|$ is the MLE for $|\rho|$. If it is desired to report a $p$-value, consider a potentially conservative $p\text{-value}=\bar{F}_{\rho}(-|\lambda_{n}^{(0:1)}|)$. For $\rho=1$ and $\rho=\hat{\rho}_{n}$ respectively, the uniform rule and the adaptive rule can be then restated as

[TABLE]

The conservative $p$-value can be computed numerically by Monte Carlo and then taking the maximum over a grid of values for $\gamma$.

Theorem 1**.**

The adaptive rule $\phi^{\text{ada}}_{n}$ controls asymptotic error uniformly below $\alpha$ for $0<\alpha<1/2$.

Proof.

We show error guarantee when ${\mathcal{M}}_{0}$ is true. The same argument holds when ${\mathcal{M}}_{1}$ is true. It suffices to show for any converging sequence $\Sigma_{n}\in{\mathcal{M}}_{0}\setminus{\mathcal{M}}_{1}$,

[TABLE]

where the probability is measured under $P_{\Sigma_{n}}^{n}$. Suppose $\Sigma_{n}\rightarrow\Sigma^{\ast}\in{\mathcal{M}}_{0}\cup{\mathcal{M}}_{1}$. If $\Sigma^{\ast}\notin{\mathcal{M}}_{0}\cap{\mathcal{M}}_{1}$, then $\lambda_{n}^{(0:1)}$ is unbounded in probability towards $+\infty$. Hence $\Pr(\phi_{n}(S_{n})={\mathcal{M}}_{0})\rightarrow 1$ and $p_{\text{err}}(\Sigma_{n})=0$. In the following we prove the claim for $\Sigma^{\ast}\in{\mathcal{M}}_{0}\cap{\mathcal{M}}_{1}$. Suppose $\hat{\rho}_{ij,n}$, $\rho_{ij,n}$ and $\rho_{ij}$ respectively denote the corresponding correlation coefficient of $S_{n}$, $\Sigma_{n}$ and $\Sigma^{\ast}$. We have three cases depending on the rate at which $\Sigma_{n}$ converges.

When $|\rho_{13,n}\rho_{23,n}|\asymp 1/\sqrt{n}$, there are two regimes depending on $\Sigma^{\ast}$.

1. (a)

In the weak-strong regime, without loss of generality suppose $\sqrt{n}\rho_{13,n}\rightarrow\gamma\neq 0$ and $\rho_{23,n}\rightarrow\rho\neq 0$. By consistency $\hat{\rho}_{n}=|\hat{\rho}_{13,n}|\vee|\hat{\rho}_{23,n}|\rightarrow_{p}|\rho|$ and the definition of envelope, we have

[TABLE] 2. (b)

In the weak-weak regime, suppose $\sqrt{n}\rho_{13,n}\rho_{23,n}\rightarrow\delta\neq 0$. We have $\hat{\rho}_{n}=|\hat{\rho}_{13,n}|\vee|\hat{\rho}_{23,n}|=(|\rho_{13,n}|\vee|\rho_{23,n}|)+O_{p}(1/\sqrt{n})\rightarrow_{p}0$ since both $\rho_{13,n},\rho_{23,n}\rightarrow 0$. By Proposition 11, we have

[TABLE] 2. 2.

When $|\rho_{13,n}\rho_{23,n}|=o(1/\sqrt{n})$, we have $\lambda_{n}^{(0:1)}\rightarrow_{p}0$. Since $F^{-1}_{\rho}(\alpha)<c<0$ for $\alpha<1/2$, we have $\Pr(\phi_{n}(S_{n})={\mathcal{M}}_{0}\cup{\mathcal{M}}_{1})\rightarrow 1$. 3. 3.

When $|\rho_{13,n}\rho_{23,n}|=\omega(1/\sqrt{n})$, we have $\Pr(\lambda_{n}^{(0:1)}>c)\rightarrow 1$ for any constant $c$ and hence $\Pr(\phi_{n}(S_{n})={\mathcal{M}}_{0})\rightarrow 1$.

∎

As can be seen from the proof, the consistency of model selection based on the loglikelihood (or AIC/BIC since in this case ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ have the same dimensions) is a special case when $|\rho_{13,n}\rho_{23,n}|=\omega(n^{-1/2})$, i.e., under strong signal or large enough sample size. However, under $|\rho_{13,n}\rho_{23,n}|=O(n^{-1/2})$, simply choosing the model with the highest loglikelihood (or the lowest AIC/BIC) can lead to large errors, as we will illustrate in the next section. Note that Theorem 1 provides a “rate-free” guarantee, in the sense that it does not require any a priori assumption on the rate of signal strength relative to the sample size. The envelope of envelopes leads to the same guarantee for the uniform rule.

Corollary 4**.**

The decision rule $\phi_{n}^{\text{unif}}$ controls asymptotic error uniformly below $\alpha$ for $0<\alpha<1/2$.

Proof.

It follows from Theorem 1 and Proposition 12. ∎

The uniform rule can be easily applied by comparing the difference in log-likelihoods to a single number, e.g., 3.19 for $\alpha=0.05$ and 5.97 for $\alpha=0.01$. The adaptive rule can be implemented by numerically evaluating $\bar{F}_{\rho}^{-1}(\alpha)$ via Monte Carlo on a grid of $\rho$ and interpolating. Some values are plotted in Fig. 13 and tabulated in Table 1 based on $10^{7}$ samples. It is interesting to notice that $\bar{F}^{-1}_{\rho}(\alpha)$ is not monotonic in $\rho\in[0,1]$.

7 Simulations

In this section we conduct numerical simulations to assess the performance of the adaptive and uniform decision rules proposed in the previous section. In subsequent simulations we use $\alpha=0.05$. In addition to the two methods we propose, we also consider the following methods for comparison.

Naive

The naive procedure selects the model with a higher likelihood (or equivalently, a lower AIC/BIC, since ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ have the same dimensions), namely

[TABLE]

This is effectively choosing a single model based on AIC/BIC since the penalty terms cancel out as ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$ have the same dimension.

Interval Selection

This method is adapted from Drton and Perlman (2004). We construct $(1-\alpha)$-level non-simultaneous confidence intervals on correlation coefficients $\rho_{12}$ and $\rho_{12\cdot 3}$ with Fisher’s $z$-transform (Fisher, 1924). The decision rule is

[TABLE]

Note that the interval selection method controls asymptotic error below $\alpha$. For example, when ${\mathcal{M}}_{0}$ is true,

[TABLE]

We conduct numerical simulations in the following three settings.

7.1 Local hypotheses

We simulate under $\Sigma_{n}^{(0)}$ and $\Sigma_{n}^{(1)}$ (variances are set to unity) for the two regimes considered in Section 4. The power is compared to the theoretically optimal. Since exact values of the total variation distance are intractable, we plot bounds given by Eq. 16 in grey curves. We perform 4,000 replications for each point on the graphs.

See Figures 14 and 15 for the size and power in the weak-strong regime (Eqs. 21 and 22) under $n=1,000$. Smaller sample sizes $n=100,200,\dots$ generate very similar results. See Figure 16 for the size and power in the weak-weak regime (Eqs. 25 and 26), where we set $\rho_{13,n}=n^{-a/4}$, $\rho_{23,n}=n^{-1/2+a/4}$ and let $a$ vary. We observe that (i) the naive method does not control error at all; (ii) the other three methods control error uniformly even under relatively small $n$. We also observe that the relation “adaptive” $>$ “uniform” $>$ “interval” holds in general in terms of both size and power. By comparing to the grey curves, we regard the adaptive rule as achieving near-optimal power in these settings.

7.2 Projected Wishart

We generate a covariance matrix by firstly drawing $\tilde{\Sigma}$ from the Wishart distribution (with the scale matrix chosen as $\sigma_{ij}=(-1/2)^{|i-j|}$) and then projecting $\tilde{\Sigma}$ into ${\mathcal{M}}_{0}$ or ${\mathcal{M}}_{1}$ respectively by finding the MLE under each model. Then we perform model selection based on two sets of zero-mean Gaussian samples generated with the two projected covariances respectively. We vary the degrees of freedom for the Wishart distribution. See Figure 17 for the results.

7.3 Conditional on covariates

We consider the common regression setting where two types of independences are contrasted conditional on a set of covariates $X\in\mathbb{R}^{p}$. In other words, we want to select between \mathcal{M}_{0}:Y_{1}\mathrel{\text{\scalebox{1.07}{\perp\mkern-10.0mu\perp}}}Y_{2}\mid X and \mathcal{M}_{1}:Y_{1}\mathrel{\text{\scalebox{1.07}{\perp\mkern-10.0mu\perp}}}Y_{2}\mid Y_{3},X. We generate instances by

[TABLE]

where we use the previous projected Wishart to generate error covariance $\Sigma$ under ${\mathcal{M}}_{0}$ and ${\mathcal{M}}_{1}$. We perform model selection by firstly regressing $(Y_{1},Y_{2},Y_{3})$ onto $X$ with least squares and then apply the model selection procedures to the residual covariance. Covariates are randomly drawn from standard Gaussians and regression coefficients are generated from a $t$-distribution with $4$ degrees of freedom. We fix $n=1,000$ and vary the number of covariates $p$. The results are presented in Figure 18. We observe that the proposed procedure continues to maintain nominal size until $p$ is relatively large compared to $n$. The power performance, on the other hand, does not seem to vary much as $p$ grows.

8 Real data example

In this section we showcase an example of applying the method to edge orientation in learning a DAG. In studying the American occupational structure, Blau and Duncan (1967) measured the following covariates on $n=20,700$ subjects:

• $V$: father’s educational attainment,
• $X$: father’s occupational status,
• $U$: child’s educational attainment,
• $W$: status of child’s first job,
• $Y$: status of child’s occupation in 1962.

The data is summarized as the following correlation matrix of $(V,X,U,W,Y)$

[TABLE]

At level $\alpha=0.01$, the PC algorithm identifies the skeleton by $d$-separation, which only removes the edge between $V$ and $Y$ based on Y\mathrel{\text{\scalebox{1.07}{\perp\mkern-10.0mu\perp}}}V\mid U,X. This is because the PC algorithm tests for conditional independence given smaller conditioning sets first. By a common-sense temporal ordering $\{V,X\}<U<\{W,Y\}$ among the variables, edges can be oriented except for $X-V$ and $W-Y$; see Fig. 19. The edge $V-X$ does not involve a collider and the orientation is statistically unidentifiable.

However, the orientation of $W-Y$ raises the interesting question of testing

[TABLE]

We apply our method to the conditional correlation of $(V,W,Y)$ given $(U,X)$. We have $\lambda_{n}^{(0:1)}=3.72$ and $p\text{-value}=0.026$ under the envelope distribution $\bar{F}_{\hat{\rho}_{n}}$. Therefore, under $\alpha=0.01$ the adaptive procedure would choose ${\mathcal{M}}_{0}\cup{\mathcal{M}}_{1}$ and leave the orientation undetermined (the procedure would choose $Y\rightarrow W$ under $\alpha=0.05$). This example illustrates the potential ambiguity in model selection even under a large sample size. The reader is referred to Spirtes et al. (2000, Section 5.8.4) for another discussion of the same example.

9 Discussion

We have considered choosing between marginal independence and conditional independence in a Gaussian graphical model, assuming we know at least one of them is true. The loglikelihood ratio statistic converges to a tight law under a sequence of truths converging to the intersection of the two models at a certain rate. The asymptotic distribution is shown to be non-uniform as it depends on where and how the sequence converges. We address this non-uniformity issue by introducing a family of envelope distributions that are well-behaved and bring back the continuity of asymptotic laws, as indexed by a parameter that can be consistently estimated. Contrary to the usual Neyman–Pearson hypothesis testing, we treat the two models symmetrically and develop model selection rules that choose both models when they are indistinguishable under a given sample size. Such rules can be designed according to the quantiles of the envelope distributions to uniformly control the type-I error below a desired level. As noted before we believe that “rate-free” asymptotic guarantees that are uniform are more useful in practice, since they do not rely upon untestable assumptions regarding the sample size and the signal strength.

In this report we restricted ourselves to the Gaussian case. For testing conditional independence, some form of distributional assumption seems inevitable, since recent work of Shah and Peters (2020) shows that testing conditional independence without restricting the form of conditional independence is impossible in general.

Selection of non-nested models routinely relies on penalized scores based on loglikelihoods, such as the negated AIC and BIC. However, as we show, in the context of a weak signal relative to the sample size, simply choosing the model with the highest score can lead to considerable errors. To select models with “confidence”, one should also look at the “gaps” between the top scores. We believe that the method developed in this paper may be generalizable to a wider range of model selection problems.

Acknowledgements

RG thanks Michael Perlman for helpful discussions. TR thanks Robin Evans and Peter Spirtes. The research was supported by the U.S. Office of Naval Research.