The conditionally studentized test for high-dimensional parametric regressions

Write the underlying regression function as $m(\boldsymbol{x})=E(Y|\boldsymbol{X}=\boldsymbol{x})$ and the error as $\varepsilon=Y-m(\boldsymbol{X})$. To study the power performance, we also consider a sequence of alternative hypotheses:

$$H_{1n}:Y=g(\boldsymbol{\theta}_{0},\boldsymbol{X})+\delta_{n}l(\boldsymbol{X})+\varepsilon.$$

(2.1)

When $\delta_{n}\to 0$ as $n\to\infty$, (2.1) corresponds to a sequence of local alternatives. When $\delta_{n}$ is a fixed constant, (2.1) reduces to the global alternative model (1.2).

Set

$\displaystyle\boldsymbol{\theta}^{*}=\underset{\boldsymbol{\theta}\in\boldsymbol{\Theta}}{\arg\min}E\left\{Y-g(\boldsymbol{\theta},\boldsymbol{X})\right\}^{2}=\underset{\boldsymbol{\theta}\in\boldsymbol{\Theta}}{\arg\min}E\left\{m(\boldsymbol{X})-g(\boldsymbol{\theta},\boldsymbol{X})\right\}^{2}.$

(2.2)

Under the null hypothesis in (1.1) and the regularity condition 1 specified in Section 6, $\boldsymbol{\theta}^{*}=\boldsymbol{\theta}_{0}$. Under the alternatives, $\boldsymbol{\theta}^{*}$ typically depends on the distribution of $\boldsymbol{X}$. To save space, redefine $l(\boldsymbol{X})=m(\boldsymbol{X})-g(\boldsymbol{\theta}^{*},\boldsymbol{X})$ under the global alternative hypothesis with fixed $\delta_{n}=1$, while we still write $l(\boldsymbol{X})=m(\boldsymbol{X})-g(\boldsymbol{\theta}_{0},\boldsymbol{X})$ under the local alternative hypothesis.

Subsequently, we proceed to present additional notations that have been employed in this paper. Denote $\dot{g}(\boldsymbol{\theta},\boldsymbol{X})=\partial g(\boldsymbol{\theta},\boldsymbol{X})/\partial\boldsymbol{\theta}$ and $\ddot{g}(\boldsymbol{\theta},\boldsymbol{X})=\partial\dot{g}(\boldsymbol{\theta},\boldsymbol{X})/\partial\boldsymbol{\theta}^{\top}$. Under the null hypothesis and the local alternative hypothesis, we define $\boldsymbol{\Sigma}=E\left\{\dot{g}(\boldsymbol{\theta}_{0},\boldsymbol{X})\dot{g}(\boldsymbol{\theta}_{0},\boldsymbol{X})^{\top}\right\}$. Under the global alternatives, without confusion, we define $\boldsymbol{\Sigma}=E\left\{\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X})\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X})^{\top}\right\}-E\left[\left\{m(\boldsymbol{X})-g(\boldsymbol{\theta}^{*},\boldsymbol{X})\right\}\ddot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X})\right]$ and $\boldsymbol{\Sigma}_{*}=E\left\{\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X})\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X})^{\top}\right\}$. Use $\|\cdot\|$ to represent the $L_{2}$ norm. Write the conditional expectation $E(A|B=b)$ as $E_{B}(A)$ for any random variable $A$ and random variable/vector $B$.

The least squares estimator of $\boldsymbol{\theta}$ is defined by

$\displaystyle\boldsymbol{\hat{\theta}}=\underset{\boldsymbol{\theta}\in\boldsymbol{\Theta}}{\arg\min}\sum_{i=1}^{n}\left\{Y_{i}-g(\boldsymbol{\theta},\boldsymbol{X}_{i})\right\}^{2}.$

(2.3)

Under the regularity conditions 1-5 in Section 6, the convergence rate and asymptotically linear representation of $\boldsymbol{\hat{\theta}}$ can be derived, which are important for studying the asymptotic properties of the test statistic under the null and alternative hypotheses. We will state them as three lemmas in Supplementary Material, and the first two lemmas are Theorems 1 and 2 and the third lemma is an extension of Theorem 4 in Tan and Zhu (2022).

In the literature, several existing tests have a similar structure, with a weight function $W_{n}(\cdot,\cdot)$ depending on the sample size $n$ such that a test statistic can be written as, before standardizing, $\sum_{i=1}^{n}\sum_{j\neq i}\hat{e}_{i}\hat{e}_{j}W_{n}(\boldsymbol{X}_{i},\boldsymbol{X}_{j})/n(n-1)$ where $\hat{e}_{i}=Y_{i}-g(\boldsymbol{\hat{\theta}},\boldsymbol{X}_{i})$ is the residual. Some classic tests are listed as follows. Bierens (1982) proposed the integrated conditional moment (ICM) test based on the weight function $\exp(-\|\boldsymbol{X}_{1}-\boldsymbol{X}_{2}\|^{2}/2)$; The tests proposed by Escanciano (2009) discussed the case with general weight functions; Li et al. (2019) used the weight function $1/\sqrt{\|\boldsymbol{X}_{1}-\boldsymbol{X}_{2}\|^{2}+1}$ induced by an idea bridging local and global smoothing tests. Other tests can be written or approximately written as U-statistics including those local smoothing tests proposed by Härdle and Mammen (1993), Zheng (1996), and those global smoothing tests suggested by Stute et al. (1998b, a), Tan and Zhu (2019) and Tan and Zhu (2022). However, these tests do not apply to the cases with diverging dimensions $p$ without dimension reduction structures.

We now construct a novel test by combining the sample-splitting and the conditional studentization approaches. The following observations induce the test construction. Let $e=Y-g(\boldsymbol{\theta}_{0},\boldsymbol{X})$. Under the null hypothesis, $m(\boldsymbol{X})=g(\boldsymbol{\theta}_{0},\boldsymbol{X})$ and $e=\varepsilon$. Note that $E(e|\boldsymbol{X})=0$ under the null hypothesis, that is, $E\left\{eE(e|\boldsymbol{X})f(\boldsymbol{X})\right\}=E\left\{E^{2}(e|\boldsymbol{X})f(\boldsymbol{X})\right\}=0$, otherwise greater than zero under the alternatives. In a more general presentation, when we use an approximation $E_{\boldsymbol{X}}\{eW_{n}(\boldsymbol{X},\boldsymbol{X}_{2})\}$ with, say, a kernel function $K\left\{(\boldsymbol{X}-\boldsymbol{X}_{2})/h\right\}/h^{q}$ in lieu of $W_{n}(\boldsymbol{X},\boldsymbol{X}_{2})$ to $E(e|\boldsymbol{X})f(\boldsymbol{X})$, this quantity can be approximated by $E\left\{e_{1}e_{2}W_{n}(\boldsymbol{X}_{1},\boldsymbol{X}_{2})\right\}$. Simply, we abbreviate $W_{n}(\boldsymbol{X}_{i},\boldsymbol{X}_{j})$ as $w_{ij}$. Thus, in general, this quantity can be estimated by $\frac{1}{n}\sum_{i=1}^{n}\hat{e}_{i}\left(\frac{1}{n-1}\sum_{j=1,j\not=i}^{n}\hat{e}_{j}w_{ij}\right)=\frac{1}{n(n-1)}\sum_{i=1}^{n}\sum_{j=1,j\not=i}^{n}\hat{e}_{i}\hat{e}_{j}w_{ij}$ that can be used to define a non-standardized statistic:

$\displaystyle U_{n}=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\hat{e}_{i}\left(\frac{1}{\sqrt{n-1}}\sum_{j=1,j\not=i}^{n}\hat{e}_{j}w_{ij}\right)=\frac{1}{\sqrt{n(n-1)}}\sum_{i=1}^{n}\sum_{j=1,j\not=i}^{n}\hat{e}_{i}\hat{e}_{j}w_{ij}.$

(2.4)

Different tests uses different weight functions $w_{ij}$. Examples include the kernel function $K\left\{(\boldsymbol{X}_{i}-\boldsymbol{X}_{j})/h\right\}/h^{q}$ with a bandwidth $h$ (see, e.g., Zheng (1996)), exponential function $\exp(-\|\boldsymbol{X}_{i}-\boldsymbol{X}_{j}\|^{2}/2)$ (see e.g., Bierens (1982)), and the weight function $1/\sqrt{\|\boldsymbol{X}_{i}-\boldsymbol{X}_{j}\|^{2}+1}$ used by Li et al. (2019). These tests often have no tractable limiting null distributions.

We now modify $U_{n}$ to avoid duplicated use of the samples. Estimate $\boldsymbol{\theta}$ by two parts of samples respectively. Define

$$\boldsymbol{\hat{\theta}}_{1}=\underset{\boldsymbol{\theta}\in\boldsymbol{\Theta}}{\arg\min}\sum_{i=1}^{n_{1}}\left\{Y_{i}-g(\boldsymbol{\theta},\boldsymbol{X}_{i})\right\}^{2}\text{ and }\boldsymbol{\hat{\theta}}_{2}=\underset{\boldsymbol{\theta}\in\boldsymbol{\Theta}}{\arg\min}\sum_{j=n_{1}+1}^{n}\left\{Y_{j}-g(\boldsymbol{\theta},\boldsymbol{X}_{j})\right\}^{2},$$

where the samples of size $n$ are divided into two disjoint parts $\mathcal{N}_{1}=\{(\boldsymbol{X}_{i},Y_{i})\}_{i=1}^{n_{1}}\text{ and }\mathcal{N}_{2}=\{(\boldsymbol{X}_{j},Y_{j})\}_{j=n_{1}+1}^{n}$ of sizes $n_{1}$ and $n_{2}$ satisfying $n=n_{1}+n_{2}$. All the results in Section 2 hold when $n$ is replaced by $n_{1}$ or $n_{2}$. Define a modified test statistic as

$$U_{\mathcal{N}_{1},\mathcal{N}_{2}}=\frac{1}{\sqrt{n_{1}}}\sum_{i=1}^{n_{1}}\hat{e}_{i}\left(\frac{1}{\sqrt{n_{2}}}\sum_{j=n_{1}+1}^{n}\hat{e}_{j}w_{ij}\right)=\frac{1}{\sqrt{n_{1}n_{2}}}\sum_{i=1}^{n_{1}}\hat{e}_{i}\sum_{j=n_{1}+1}^{n}\hat{e}_{j}w_{ij},$$

(2.5)

where $\hat{e}_{i}=Y_{i}-g(\boldsymbol{\hat{\theta}}_{1},\boldsymbol{X}_{i}),\ i=1,\dots,n_{1}$ and $\hat{e}_{j}=Y_{j}-g(\boldsymbol{\hat{\theta}}_{2},\boldsymbol{X}_{j}),\ j=n_{1}+1,\dots,n.$ Again, this test statistic usually has no significant difference from the previous $U_{n}$ as they have similar asymptotic behaviors. However, this seemingly minor modification plays a vital role in constructing a conditionally studentized test with a normal weak limit under the null hypothesis. The key ingredient in the construction is that in the two independent sums the residuals are not used duplicately, but the weight function $w_{ij}$ links the two sums.

To see how to define the final statistic, we give the decomposition of $U_{\mathcal{N}_{1},\mathcal{N}_{2}}$. Under the null hypothesis, we have the followings:

	$\displaystyle U_{\mathcal{N}_{1},\mathcal{N}_{2}}$	$\displaystyle=\frac{1}{\sqrt{n_{1}n_{2}}}\sum_{i=1}^{n_{1}}\sum_{j=n_{1}+1}^{n}\hat{e}_{i}\hat{e}_{j}w_{ij}$
		$\displaystyle=\frac{1}{\sqrt{n_{1}}}\sum_{i=1}^{n_{1}}e_{i}\frac{1}{\sqrt{n_{2}}}\sum_{j=n_{1}+1}^{n}\hat{e}_{j}\left[w_{ij}-\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X}_{i})^{\top}\boldsymbol{\Sigma}^{-1}E_{\boldsymbol{X}_{j}}\left\{\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X}_{1})w_{1j}\right\}\right]+o_{p}(1)$
		$\displaystyle=:\frac{1}{\sqrt{n_{1}}}\sum_{i=1}^{n_{1}}e_{i}\tilde{w}_{i}(\mathcal{N}_{2})+o_{p}(1),$		(2.6)

where $E_{\boldsymbol{X}_{j}}\left\{\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X}_{1})w_{1j}\right\}$ is the conditional expectation given $\boldsymbol{X}_{j}$, and $\boldsymbol{\theta}^{*}$ is defined in (2.2). The detailed justification can be found in Corollary 1 of Supplementary Material. The decomposition is only about the residuals $\hat{e}_{i}$’s in the first part of the samples. Given $\mathcal{N}_{2}$, $e_{1}\tilde{w}_{1}(\mathcal{N}_{2}),\ldots,e_{n_{1}}\tilde{w}_{n_{1}}(\mathcal{N}_{2})$ are conditionally independent and identically distributed random variables. Intuitively, under the null hypothesis, the following random sequence would have a normal weak limit by using the Central Limit Theorem conditionally:

$$V_{n}:=\frac{\frac{1}{\sqrt{n_{1}}}\sum_{i=1}^{n_{1}}e_{i}\tilde{w}_{i}(\mathcal{N}_{2})}{\sqrt{\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\left\{e_{i}\tilde{w}_{i}(\mathcal{N}_{2})-\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}e_{i}\tilde{w}_{i}(\mathcal{N}_{2})\right\}^{2}}}$$

(2.7)

that is a conditionally studentized version of $\sum_{i=1}^{n_{1}}e_{i}\tilde{w}_{i}(\mathcal{N}_{2})$. To define the final test statistic, we use $\sum_{i=1}^{n_{1}}\sum_{j=n_{1}+1}^{n}\hat{e}_{i}\hat{e}_{j}w_{ij}/\sqrt{n_{1}n_{2}}$ in lieu of $\sum_{i=1}^{n_{1}}e_{i}\tilde{w}_{i}(\mathcal{N}_{2})/\sqrt{n_{1}}$ and for the denominator that is a conditional standard deviation, we use $\hat{e}_{i}$, $\hat{\boldsymbol{\Sigma}}$ and $\boldsymbol{\hat{\theta}}$ to replace $e_{i}$, $\boldsymbol{\Sigma}$ and $\boldsymbol{\theta}^{*}$:

$$\hat{V}_{n}:=\frac{\frac{1}{\sqrt{n_{1}n_{2}}}\sum_{i=1}^{n_{1}}\sum_{j=n_{1}+1}^{n}\hat{e}_{i}\hat{e}_{j}w_{ij}}{\sqrt{\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\left\{\hat{e}_{i}\tilde{\tilde{w}}_{i}(\mathcal{N}_{2})-\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\hat{e}_{i}\tilde{\tilde{w}}_{i}(\mathcal{N}_{2})\right\}^{2}}},$$

(2.8)

where

$$\tilde{\tilde{w}}_{i}(\mathcal{N}_{2})=\frac{1}{\sqrt{n_{2}}}\sum_{j=n_{1}+1}^{n}\hat{e}_{j}\left[w_{ij}-\dot{g}\left(\boldsymbol{\hat{\theta}},\boldsymbol{X}_{i}\right)^{\top}\hat{\boldsymbol{\Sigma}}^{-1}\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\left\{\dot{g}\left(\boldsymbol{\hat{\theta}},\boldsymbol{X}_{i}\right)w_{ij}\right\}\right],$$

and $\hat{\boldsymbol{\Sigma}}=\frac{1}{n}\sum_{i=1}^{n}\dot{g}(\boldsymbol{\hat{\theta}},\boldsymbol{X}_{i})\dot{g}(\boldsymbol{\hat{\theta}},\boldsymbol{X}_{i})^{\top}$. We will prove that $\boldsymbol{\hat{\theta}}$, $\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\{\dot{g}(\boldsymbol{\hat{\theta}},\boldsymbol{X}_{i})w_{ij}\}$ and $\hat{\boldsymbol{\Sigma}}$ are the consistent estimators of $\boldsymbol{\theta}^{*}$, $E_{\boldsymbol{X}_{j}}\left\{\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X}_{1})w_{1j}\right\}$ and $\boldsymbol{\Sigma}_{*}$, respectively such that the consistency of this estimated conditional standard deviation holds. Note that in $\tilde{\tilde{w}}_{i}(\mathcal{N}_{2})$, we use the full data-based estimator $\boldsymbol{\hat{\theta}}$ instead of $\boldsymbol{\hat{\theta}}_{2}$. We find that asymptotically, there is no difference by using either estimator, but using $\boldsymbol{\hat{\theta}}$ can make a faster convergence rate of the estimator to $\boldsymbol{\theta}^{*}$. In Corollary 2 of Supplementary Material, we will show that under the null hypothesis and the local alternative hypothesis in (2.1) with $np^{3}\delta_{n}^{4}\to 0$,

$$\frac{\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\left\{e_{i}\tilde{w}_{i}(\mathcal{N}_{2})-\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}e_{i}\tilde{w}_{i}(\mathcal{N}_{2})\right\}^{2}}{\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\left\{\hat{e}_{i}\tilde{\tilde{w}}_{i}(\mathcal{N}_{2})-\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}\hat{e}_{i}\tilde{\tilde{w}}_{i}(\mathcal{N}_{2})\right\}^{2}}\stackrel{{\scriptstyle p}}{{\to}}1,$$

(2.9)

where $``\stackrel{{\scriptstyle p}}{{\to}}"$ stands for convergence in probability. While, under the global alternative hypothesis in (2.1) with fixed $\delta_{n}$, $\tilde{\tilde{w}}_{i}(\mathcal{N}_{2})$ will be proved to be a consistent estimator of $\tilde{w}_{i}^{(0)}(\mathcal{N}_{2}):=\sum_{j=n_{1}+1}^{n}\hat{e}_{j}\left[w_{ij}-\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X}_{i})^{\top}\boldsymbol{\Sigma}_{*}^{-1}E_{\boldsymbol{X}_{j}}\left\{\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X}_{1})w_{1j}\right\}\right]/\sqrt{n_{2}}$. Note that when a kernel function $K\left\{(\boldsymbol{X}-\boldsymbol{X}_{2})/h\right\}/h^{q}$ is used as the weight function $W_{n}(\boldsymbol{X},\boldsymbol{X}_{2})$, it can go to infinity as the sample size goes to infinity. But this is not a problem in our construction as the studentized test is scale-invariant, and the quantity $h^{q}$ in the weight $w_{ij}$ is eliminated from the numerator and denominator. Thus, we can consider the weight function without this quantity.

Under the null hypothesis, the conditionally studentized version $V_{n}$ of $\sum_{i=1}^{n_{1}}e_{i}\tilde{w}_{i}(\mathcal{N}_{2})$ has a normal weak limit under some regularity conditions (see Corollary 3 in Supplementary Material for details). We can prove the asymptotic equivalence between the numerator(denominator) of $\hat{V}_{n}$ and $V_{n}$. The result is stated as follows.

Therefore, we can compute the critical values easily.

Consider the power performance under the alternative hypothesis in (2.1). The fixed non-zero constant $\delta_{n}$ corresponds to the global alternative hypothesis with $Y=m(\boldsymbol{X})+\varepsilon$, where $m(\boldsymbol{X})\neq g(\boldsymbol{\theta},\boldsymbol{X})$ for $\forall\boldsymbol{\theta}\in\boldsymbol{\Theta}$, and $\delta_{n}\to 0$ as $n\to\infty$ to local alternatives. Define $\boldsymbol{\Sigma}_{\varepsilon}=E\left\{\varepsilon^{2}\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X})\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X})^{\top}\right\},$ $\boldsymbol{\Sigma}_{l}=E\left\{l^{2}(\boldsymbol{X})\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X})\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X})^{\top}\right\}$ and $\boldsymbol{\Sigma}_{cov}^{*}=E\left\{\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X}_{1})\dot{g}(\boldsymbol{\theta}^{*},\boldsymbol{X}_{2})^{\top}w_{12}\right\}.$ We state the following results.

In this section, some numerical studies are conducted to examine the performance of our test proposed in Section 3. As Li et al. (2019) used, the weight function $1/\sqrt{\|\boldsymbol{X}_{1}-\boldsymbol{X}_{2}\|^{2}+1}$ has the merits that combine local smoothing and global smoothing tests. But it is theoretically flawed for large $p$ as it converges to a constant when $p$ goes to infinity. To remedy this defect, not only should our chosen weight function include this weight function, but also it contains another weight function $\sum_{k=1}^{p}K\left\{(X_{ik}-X_{jk})/h\right\}/h$, where $K(\cdot)$ denotes the kernel density function and $h$ is the bandwidth. The summation form ensures that it works in diverging dimension cases. As a result, the weight function, defined as $0.5\times\left[1/\sqrt{\|\boldsymbol{X}_{1}-\boldsymbol{X}_{2}\|^{2}+1}+\sum_{k=1}^{p}K\left\{(X_{ik}-X_{jk})/h\right\}/h\right]$, is a hybrid of two equally weighted functions.

To make the simulation results convincing, we compare with the test $AICM_{n}$ proposed by Tan and Zhu (2022) that can also handle diverging dimension cases and shows its advantages. It is worth noticing that the method in Janková et al. (2020) could also be applied to non-sparse models with random designs and dimensions under our constraints, Tan and Zhu (2022) made a comparison for Logistic model with three model settings, and $AICM_{n}$ outperformed the test in Janková et al. (2020) in two of the three model settings. We have also conducted the simulations under the same model settings and found that $COST$ worked similarly to $AICM_{n}$. Therefore, we put the results in Section 4.2, and the main text here only reports the numerical comparison with $AICM_{n}$. In addition, we also evaluate the performance of another related test, $DrCost_{n}$, which can be seen as a modified version of $Cost_{n}$. When the underlying model has a dimension reduction structure, $DrCost_{n}$ contains the dimension reduction step (see, Tan and Zhu (2022)) . We design three studies: models with dimension reduction structures; with dimension reduction structures and diverging dimensions; and without dimension reduction structures. The first two studies favor Tan and Zhu (2022)’s method, and the third study deals with general models. The predictor vectors $\boldsymbol{X}_{i}$’s are independently generated from multivariate normal distribution $\mathbf{N}(\boldsymbol{0},\boldsymbol{\Sigma})$. Here $\boldsymbol{\Sigma}=\boldsymbol{\Sigma}_{1}=I_{p}$ or $\boldsymbol{\Sigma}=\boldsymbol{\Sigma}_{2}=(0.5^{|i-j|})_{p\times p}$. The error $\varepsilon$’s are independently drawn from the standard normal distribution $\mathbf{N}(0,1)$. The simulation results are all based on $1,000$ replications. The results of $AICM_{n}$ are computed by the code provided by the authors of the paper Tan and Zhu (2022). To choose the bandwidth $h$, we consider it to be $c\cdot n^{-0.2}$, where $c$ is a constant taking the following five values: $0.5,0.8,1,1.2,1.5$. To check how robust the test is against the value $c$, we use the model $H_{11}$ as an example in Figure 1. It shows that our test performs robustly, and the size and power level can be similar. Therefore, we use the bandwidth with $c=1$.

For comparisons, we consider four numerical studies. Studies 1 and 2 have multi-index model structures that favor $AICM_{n}$, and Study 3 does not have such a structure. As under some strong conditions on the regression function our method can handle the cases where the dimension $q$ of the predictor vector can be much higher than that $p$ of parameter vector, we consider Study 4 that has $q=p^{2}$ and $q=n$. As the limiting null distribution of $AICM_{n}$ is generally intractable, the wild bootstrap approximation is used to determine the critical values. For model $H_{11}$ in Study 1 and model $H_{21}$ in Study 2, the numerical results of $AICM_{n}$ are excerpted from Tan and Zhu (2022) to make the paper self-contained.

Study 1. Generate data from the following double-index and triple-index models:

	$\displaystyle H_{11}:Y=\boldsymbol{\beta}_{1}^{\top}\boldsymbol{X}+a\left(\boldsymbol{\beta}_{2}^{\top}\boldsymbol{X}\right)^{2}+\varepsilon,$
	$\displaystyle H_{12}:Y=X_{1}+\cos(2X_{2})+a\exp(3X_{2})+\varepsilon.$

We set $\boldsymbol{\beta}_{1}=(\underbrace{1,\cdots,1}_{q_{1}},0,\cdots,0)^{\top}/\sqrt{q_{1}}$ and $\boldsymbol{\beta}_{2}=(0,\cdots,0,\underbrace{1,\cdots,1}_{q_{1}})^{\top}/\sqrt{q_{1}}$ with $q_{1}=[q/2]_{-}$, where $[\quad]_{-}$ indicates the largest integer smaller than or equal to $q/2$. Here $X_{i}$ denotes the $i$-th component of $\boldsymbol{X}$. The first hypothetical model is single-index, and the second is double-index, containing the second index and third index in alternative models respectively. The empirical sizes and powers are reported in Tables 1-2 in Section 4.2.

Table 1 shows that $AICM_{n}$ works better for model $H_{11}$, and $DrCost_{n}$ with the dimension reduction step outperforms $Cost_{n}$, but slightly. When the sample size gets larger, our tests gradually work closer to $AICM_{n}$. However, for model $H_{12}$, the situation changes. The results of Table 2 suggest that $Cost_{n}$ outperforms $DrCost_{n}$ and $AICM_{n}$, although the model favors them. We checked the details and found that the structural dimension of the central subspace is underestimated to $1$ in this model and the residuals cannot be estimated well under the alternative hypothesis. This might be one of the reasons that $AICM_{n}$ and $DrCost_{n}$ do not work well.

Study 2. Generate data from the following multi-index models with $q=0.1n$ and $q=\sqrt{n}$ respectively:

	$\displaystyle H_{21}:Y=\boldsymbol{\beta}_{0}^{\top}\boldsymbol{X}+a\exp\left(\boldsymbol{\beta}_{0}^{\top}\boldsymbol{X}\right)+\varepsilon,$
	$\displaystyle H_{22}:Y=\boldsymbol{\beta}_{1}^{\top}\boldsymbol{X}+\exp\left(\boldsymbol{\beta}_{2}^{\top}\boldsymbol{X}\right)+a\exp\left(-\boldsymbol{\beta}_{0}^{\top}\boldsymbol{X}\right)+\varepsilon,$

where $\boldsymbol{\beta}_{0}=(1,1,\cdots,1)^{\top}/\sqrt{q}$ and other notations are the same as stated above. In this study, the dimension $q$ is large; thus, the regularity conditions fail to hold in theory. The empirical sizes and powers of Study 2 are presented in Tables 3-4 in Section 4.2.

The simulation results suggest that all competitors cannot maintain the significance level for model $H_{21}$ when $q$ diverges at the rate of $0.1n$. This means that the condition of dimensionality is violated too much. For model $H_{22}$ with $q=\sqrt{n}$ that also violates the condition, we can see that when the sample size is large, our test performs well in the significance level maintenance with relatively high power. At the same time, $AICM_{n}$ is liberal in general. The test $DrCost_{n}$ with the dimension reduction step works slightly worse than $Cost_{n}$. This suggests that the test may still be usable when the sample size is large.