In order to verify the above results, we repeated the simulation study from
Section 2.5 but with Gaussian covariates, considering only the naive estimator
$\hat{\tau}^{2}$ and the OOE estimator $T_{oracle}.$ For the low-dimensional
case, we fixed $p=3$ and considered
$\beta_{j}^{2}=\frac{{\tau_{\bf{B}}^{2}}}{2}$ for $j=1,\dots,2$, and
$\beta_{3}^{2}={{\tau^{2}-\tau_{\bf{B}}^{2}}}.$ Table 4 suggests that the OOE
estimator achieves similar reduction in variance as claimed in theory, namely
that for the low-dimensional case, the reduction is about 66% and when $n=p$
the reduction is about 40%.

Table 4: Summary statistics for the naive and the OOE estimators; $n=300$;
$\tau^{2}=\sigma^{2}=1.$ Mean, mean square error (MSE) and percentage change
from the naive estimator (in terms of MSE) are shown. Simulation standard
errors are shown in parenthesis. The table results were computed over $100$
simulated datasets for each setting.

$p$ | $\tau^{2}_{\bf{B}}$ | Estimator | Mean | MSE | % Change
---|---|---|---|---|---
300 | 5% | naive | 1 (0.015) | 0.068 (0.106) | 0
300 | 5% | OOE | 1.01 (0.012) | 0.042 (0.07) | -38.24
300 | 35% | naive | 1.01 (0.015) | 0.068 (0.095) | 0
300 | 35% | OOE | 1.02 (0.011) | 0.039 (0.059) | -42.65
300 | 65% | naive | 1.01 (0.015) | 0.069 (0.101) | 0
300 | 65% | OOE | 1.01 (0.011) | 0.038 (0.056) | -44.93
300 | 95% | naive | 1.01 (0.015) | 0.072 (0.11) | 0
300 | 95% | OOE | 1.01 (0.012) | 0.04 (0.058) | -44.44
3 | 5% | naive | 0.99 (0.011) | 0.033 (0.05) | 0
3 | 5% | OOE | 1.01 (0.007) | 0.014 (0.022) | -57.58
3 | 35% | naive | 1.01 (0.012) | 0.043 (0.069) | 0
3 | 35% | OOE | 1.01 (0.007) | 0.016 (0.025) | -62.79
3 | 65% | naive | 1.01 (0.013) | 0.052 (0.087) | 0
3 | 65% | OOE | 1.01 (0.008) | 0.018 (0.029) | -65.38
3 | 95% | naive | 1.01 (0.013) | 0.049 (0.083) | 0
3 | 95% | OOE | 1.01 (0.007) | 0.017 (0.029) | -65.31

Proof of Proposition 2.3:
Write,

$\displaystyle\text{Var}\left(T\right)$
$\displaystyle=\text{Var}\left[{{{\hat{\tau}}^{2}}-2\sum\limits_{j=1}^{p}{\sum\limits_{j^{\prime}=1}^{p}{{{\hat{\psi}}_{jj^{\prime}}}}}}\right]$
$\displaystyle=\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-4{\text{Cov}}\left({{{\hat{\tau}}^{2}},\sum\limits_{j=1}^{p}{\sum\limits_{j^{\prime}=1}^{p}{{{\hat{\psi}}_{jj^{\prime}}}}}}\right)+4\text{Var}\left({\sum\limits_{j=1}^{p}{\sum\limits_{j^{\prime}=1}^{p}{{{\hat{\psi}}_{jj^{\prime}}}}}}\right).$
(2.50)

We start with calculating the middle term. Let $p_{n}(k)\equiv
n(n-1)(n-2)\cdots(n-k).$ Write,

$\displaystyle{\text{Cov}}\big{(}{{{\hat{\tau}}^{2}},\sum\limits_{j=1}^{p}{\sum\limits_{j^{\prime}=1}^{p}{{{\hat{\psi}}_{jj^{\prime}}}}}}\big{)}$
$\displaystyle={\text{Cov}}\left({\frac{1}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{{W_{{i_{1}}j}}{W_{{i_{2}}j}}}},\frac{1}{{n\left({n-1}\right)\left({n-2}\right)}}\sum\limits_{j,j^{\prime}}{\sum\limits_{{i_{1}}\neq{i_{2}}\neq{i_{3}}}{{W_{{i_{1}}j}}{W_{{i_{2}}j^{\prime}}}\left[{{X_{{i_{3}}j}}{X_{{i_{3}}j^{\prime}}}-E\left({{X_{{i_{3}}j}}{X_{{i_{3}}j^{\prime}}}}\right)}\right]}}}\right)$
$\displaystyle=C_{n}\sum\limits_{I}{\sum\limits_{J}{{\text{Cov}}\left({{W_{{i_{1}}{j_{1}}}}{W_{{i_{2}}{j_{1}}}},{W_{{i_{3}}{j_{2}}}}{W_{{i_{4}}{j_{3}}}}\left[{{X_{{i_{5}}{j_{2}}}}{X_{{i_{5}}{j_{3}}}}-E\left({{X_{{i_{5}}{j_{2}}}}{X_{{i_{5}}{j_{3}}}}}\right)}\right]}\right)}},$
(2.51)

where $C_{n}\equiv\frac{1}{p_{n}(1)\cdot p_{n}(2)},$ $I$ is the set of all
quintuples of indices $(i_{1},i_{2},i_{3},i_{4},i_{5})$ such that $i_{1}\neq
i_{2}$ and $i_{3}\neq i_{4}\neq i_{5},$ and $J$ is the set of all triples of
indices $(j_{1},j_{2},j_{3})$. For the set $I$, there are $\binom{2}{1}\cdot
3=6$ different cases to consider when one of $\\{i_{1},i_{2}\\}$ is equal to
one of $\\{i_{3},i_{4},i_{5}\\}$, and an additional $\binom{2}{2}\cdot 3!=6$
cases to consider when two of $\\{i_{1},i_{2}\\}$ are equal to two of
$\\{i_{3},i_{4},i_{5}\\}$. Similarly, for the set $J$ there are three cases to
consider when only two indices of $\\{j_{1},j_{2},j_{3}\\}$ are equal to each
other, (e.g., $j_{1}=j_{2}\neq j_{3}$) ; one case to consider when no pair of
indices is equal to each other and; one case to consider when all three
indices are equal. Thus, there are total of $(6+6)\times(3+1+1)=60$ cases to
consider. Here we demonstrate only one such case. Let
${I_{1}}=\left\\{{\left({{i_{1}},\dots,{i_{5}}}\right):{i_{1}}={i_{5}}\neq{i_{2}}\neq{i_{3}}\,\neq\,{i_{4}}}\right\\}$
and
${J_{1}}=\left\\{{\left({{j_{1}},{j_{2}},{j_{3}}}\right):{j_{1}}={j_{2}}={j_{3}}}\right\\}.$
Write,

$\begin{split}&C_{n}\sum\limits_{I_{1}}{\sum\limits_{J_{1}}{{\text{Cov}}\left({{W_{{i_{1}}{j_{1}}}}{W_{{i_{2}}{j_{1}}}},{W_{{i_{3}}{j_{2}}}}{W_{{i_{4}}{j_{3}}}}\left[{{X_{{i_{5}}{j_{2}}}}{X_{{i_{5}}{j_{3}}}}-E\left({{X_{{i_{5}}{j_{2}}}}{X_{{i_{5}}{j_{3}}}}}\right)}\right]}\right)}}\\\
&=C_{n}\sum\limits_{I_{1}}{\sum\limits_{j=1}^{p}{\text{Cov}\left({{W_{{i_{1}}{j}}}{W_{{i_{2}}{j}}},{W_{{i_{3}}{j}}}{W_{{i_{4}}{j}}}\left[{X_{{i_{1}}{j}}^{2}-1}\right]}\right)}}\\\
&=C_{n}\sum\limits_{I_{1}}{\sum\limits_{j=1}^{p}{E({{W_{{i_{2}}{j}}})E({W_{{i_{3}}{j}}})E({W_{{i_{4}}{j}}})}E\left({{W_{{i_{1}}{j}}}\left[{X_{{i_{1}}{j}}^{2}-1}\right]}\right)}}\\\
&=C_{n}\sum\limits_{I_{1}}{\sum\limits_{j=1}^{p}\beta_{{j}}^{3}E\left({{W_{{i}{j}}}\left[{X_{{i}{j}}^{2}-1}\right]}\right)}.\end{split}$
(2.52)

Now, notice that

$\begin{split}E\left[{{W_{ij}}\left({X_{ij}^{2}-1}\right)}\right]&=E\left[{{X_{ij}}{Y_{i}}\left({X_{ij}^{2}-1}\right)}\right]\\\
&=E\left[{X_{ij}^{3}\left({{\beta^{T}}X+{\varepsilon_{i}}}\right)}\right]-{\beta_{j}}\\\
&={\beta_{j}}E\left({X_{ij}^{4}}\right)-{\beta_{j}}\\\
&=\beta_{j}[E(X_{ij}^{4})-1].\end{split}$ (2.53)

Rewrite (2.52) to get

$\displaystyle
C_{n}\sum\limits_{I_{1}}{\sum\limits_{j=1}^{p}{\beta_{{j}}^{3}E\left[{W_{i{j}}\left({X_{i{j}}^{2}-1}\right)}\right]}}$
$\displaystyle=C_{n}\sum\limits_{I_{1}}{\sum\limits_{j=1}^{p}{\beta_{j}^{3}\left(\beta_{j}[E(X_{ij}^{4})-1]\right)}}$
$\displaystyle=\frac{p_{n}(3)}{p_{n}(1)\cdot
p_{n}(2)}\sum\limits_{j=1}^{p}{\beta_{j}^{4}}[E(X_{ij}^{4})-1]$
$\displaystyle=\frac{(n-3)}{n(n-1)}\sum\limits_{j=1}^{p}{\beta_{j}^{4}}[E(X_{ij}^{4})-1]$
$\displaystyle=\frac{1}{n}\sum\limits_{j=1}^{p}{\beta_{j}^{4}}[E(X_{ij}^{4})-1]+O(n^{-2}),$

where we used (2.53) to justify the first equality. By the same type of
calculation, one can compute the covariance in (2.51) over all 60 and obtain
that

${\text{Cov}}\big{(}{{{\hat{\tau}}^{2}},\sum\limits_{j=1}^{p}{\sum\limits_{j^{\prime}=1}^{p}{{{\hat{\psi}}_{jj^{\prime}}}}}}\big{)}=\frac{2}{n}\left\\{{\sum\limits_{j=1}^{p}{\beta_{j}^{4}\left[{E\left({X_{1j}^{4}-1}\right)}\right]+2\sum\limits_{j\neq
j^{\prime}}{\beta_{j}^{2}\beta_{j^{\prime}}^{2}}}}\right\\}+O\left({{n^{-2}}}\right).$
(2.54)

We now move to calculate the last term of (Supplementary Material). Recall
that

$\hat{\psi}_{jj^{\prime}}=\frac{1}{{n\left({n-1}\right)\left({n-2}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}\neq{i_{3}}}{{W_{{i_{1}}j}}{W_{{i_{2}}j^{\prime}}}\left[{{X_{{i_{3}}j}}{X_{{i_{3}}j^{\prime}}}-E\left({{X_{{i_{3}}j}}{X_{{i_{3}}j^{\prime}}}}\right)}\right]}.$

Therefore,

$\displaystyle\text{Var}\big{(}{\sum\limits_{j=1}^{p}{\sum\limits_{j^{\prime}=1}^{p}{{{\hat{\psi}}_{jj^{\prime}}}}}}\big{)}=\sum\limits_{J}{{\text{Cov}}\left({{{\hat{\psi}}_{{j_{1}}{j_{2}}}},{{\hat{\psi}}_{{j_{3}}{j_{4}}}}}\right)}$
(2.55)
$\displaystyle=\frac{1}{{{{\left[{n\left({n-1}\right)\left({n-2}\right)}\right]}^{2}}}}\sum\limits_{J}{{\text{Cov}}\left({\sum\limits_{{i_{1}}\neq{i_{2}}\neq{i_{3}}}{{W_{{i_{1}}{j_{1}}}}{W_{{i_{2}}{j_{2}}}}{X_{{i_{3}}{j_{1}}}}{X_{{i_{3}}{j_{2}}}}},\sum\limits_{{i_{1}}\neq{i_{2}}\neq{i_{3}}}{{W_{{i_{1}}{j_{3}}}}{W_{{i_{2}}{j_{4}}}}{X_{{i_{3}}{j_{3}}}}{X_{{i_{3}}{j_{4}}}}}}\right)}$
$\displaystyle=p_{n}^{-2}(2)\sum\limits_{J}{\sum\limits_{I}{{\text{Cov}}\left({{W_{{i_{1}}{j_{1}}}}{W_{{i_{2}}{j_{2}}}}{X_{{i_{3}}{j_{1}}}}{X_{{i_{3}}{j_{2}}}},{W_{{i_{4}}{j_{3}}}}{W_{{i_{5}}{j_{4}}}}{X_{{i_{6}}{j_{3}}}}{X_{{i_{6}}{j_{4}}}}}\right)}\,,}$

where $J$ is now defined to be the set of all quadruples
$(j_{1},j_{2},j_{3},j_{4}),$ and $I$ is now defined to be the set of all
sextuples $(i_{1},...,i_{6})$ such that $i_{1}\neq i_{2}\neq i_{3}$ and
$i_{4}\neq i_{5}\neq i_{6}.$ For the set $I$, there are three different cases
to consider: (1) when one of $\\{i_{1},i_{2},i_{3}\\}$ is equal to one of
$\\{i_{4},i_{5},i_{6}\\}$; (2) when two of $\\{i_{1},i_{2},i_{3}\\}$ are equal
to two of $\\{i_{4},i_{5},i_{6}\\};$ and (3) when $\\{i_{1},i_{2},i_{3}\\}$
are equal to $\\{i_{4},i_{5},i_{6}\\}.$ There are $\binom{3}{1}\cdot 3=9$
options for the first case, $\binom{3}{2}\cdot 3!=18$ for the second case, and
$\binom{3}{3}\cdot 3!=6$ options for the third case. For the set $J$, there
are five different cases to consider: (1) when there is only one pair of equal
indices (e.g., $j_{1}=j_{2}\neq j_{3}\neq j_{4}$); (2) when there are two
pairs of equal indices (e.g., $j_{1}=j_{2}\neq j_{3}=j_{4}$); (3) when only
three indices are equal (e.g., $j_{1}=j_{2}=j_{3}\neq j_{4}$); (4) when all
four indices are equal and; (5) all four indices are different from each
other. Note that there are $\binom{4}{2}=6$ combinations for the first case,
$\binom{4}{2}=6$ for the second case, $\binom{4}{3}=4$ combinations for the
third case, and a single combination for each of the last two cases. Thus,
there are total of $(9+18+6)\times(6+6+4+1+1)=594.$ Again we demonstrate only
one such calculation. Let
${I_{2}}=\left\\{{\left({{i_{1}},...,{i_{6}}}\right):{i_{1}}={i_{4}},{i_{2}}={i_{5}},{i_{3}}={i_{6}}}\right\\}$
and
${J_{2}}=\left\\{{\left({{j_{1}},{j_{2}},{j_{3}},{j_{4}}}\right):{j_{1}}={j_{3}}\neq{j_{2}}={j_{4}}}\right\\}.$
In the view of (2.55),

$\displaystyle
p_{n}^{-2}(2)\sum\limits_{{J_{2}}}{\sum\limits_{{I_{2}}}{\text{Cov}\left({{W_{{i_{1}}{j_{1}}}}{W_{{i_{2}}{j_{2}}}}{X_{{i_{3}}{j_{1}}}}{X_{{i_{3}}{j_{2}}}},{W_{{i_{4}}{j_{3}}}}{W_{{i_{5}}{j_{4}}}}{X_{{i_{6}}{j_{3}}}}{X_{{i_{6}}{j_{4}}}}}\right)}}=$
$\displaystyle=p_{n}^{-2}(2)\sum\limits_{{J_{2}}}{\sum\limits_{{I_{2}}}{\text{Cov}\left({{W_{{i_{1}}{j_{1}}}}{W_{{i_{2}}{j_{2}}}}{X_{{i_{3}}{j_{1}}}}{X_{{i_{3}}{j_{2}}}},{W_{{i_{1}}{j_{1}}}}{W_{{i_{2}}{j_{2}}}}{X_{{i_{3}}{j_{1}}}}{X_{{i_{3}}{j_{2}}}}}\right)}}$
$\displaystyle=p_{n}^{-2}(2)\sum\limits_{{J_{2}}}{\sum\limits_{{I_{2}}}{E\left({W_{{i_{1}}{j_{1}}}^{2}}\right)E\left({W_{{i_{2}}{j_{2}}}^{2}}\right)E\left({X_{{i_{3}}{j_{1}}}^{2}}\right)E\left({X_{{i_{3}}{j_{2}}}^{2}}\right)}}$
$\displaystyle=p_{n}^{-2}(2)\sum\limits_{{J_{2}}}{\sum\limits_{{I_{2}}}{\left({\sigma_{Y}^{2}+\beta_{{j_{1}}}^{2}\left\\{{E\left({X_{i{j_{1}}}^{4}}\right)-1}\right\\}}\right)\left({\sigma_{Y}^{2}+\beta_{{j_{2}}}^{2}\left\\{{E\left({X_{i{j_{2}}}^{4}}\right)-1}\right\\}}\right)}}$
$\displaystyle\leq
p_{n}^{-2}(2)\sum\limits_{{J_{2}}}{\sum\limits_{{I_{2}}}{\left({\sigma_{Y}^{2}+\beta_{{j_{1}}}^{2}\left({C-1}\right)}\right)\left({\sigma_{Y}^{2}+\beta_{{j_{2}}}^{2}\left({C-1}\right)}\right)}}$
$\displaystyle=p_{n}^{-1}(2)\sum\limits_{{j_{1}}\neq{j_{2}}}{\left[{\sigma_{Y}^{4}+\sigma_{Y}^{2}\left({C-1}\right)\left({\beta_{{j_{1}}}^{2}+\beta_{{j_{2}}}^{2}}\right)+{{\left({C-1}\right)}^{2}}\beta_{{j_{1}}}^{2}\beta_{{j_{2}}}^{2}}\right]}$
$\displaystyle=p_{n}^{-1}(2)\left[{p\left({p-1}\right)\sigma_{Y}^{4}+\sigma_{Y}^{2}\left({C-1}\right)\sum\limits_{{j_{1}}\neq{j_{2}}}{\left({\beta_{{j_{1}}}^{2}+\beta_{{j_{2}}}^{2}}\right)}+{{\left({C-1}\right)}^{2}}\sum\limits_{{j_{1}}\neq{j_{2}}}{\beta_{{j_{1}}}^{2}\beta_{{j_{2}}}^{2}}}\right]$
$\displaystyle\leq
p_{n}^{-1}(2)\left[{p\left({p-1}\right)\sigma_{Y}^{4}+\sigma_{Y}^{2}\left({C-1}\right)\left({2p{\tau^{2}}}\right)+{{\left({C-1}\right)}^{2}}{\tau^{4}}}\right],$

where the fourth equality we use
$E\left({W_{ij}^{2}}\right)=\sigma_{Y}^{2}+\beta_{j}^{2}[E(X_{ij}^{2})-1]$,
which is given by (2.28), and in the fifth equality we used the assumption
that $E(X_{ij}^{4})\leq C$ for some positive $C.$ Since we assume $p/n=O(1),$
the above expression can be further simplified to
$\frac{{{p^{2}}{\sigma_{Y}^{4}}}}{{{n^{3}}}}~{}+~{}O\left({{n^{-2}}}\right).$

By the same type of calculation, one can compute the covariance in (2.55) over
all 594 cases and obtain that

$\text{Var}\big{(}{\sum\limits_{j=1}^{p}{\sum\limits_{j^{\prime}=1}^{p}{{{\hat{\psi}}_{jj^{\prime}}}}}}\big{)}=\frac{1}{n}\left\\{{\sum\limits_{j=1}^{p}{\beta_{j}^{4}\left[{E\left({X_{1j}^{4}-1}\right)}\right]+2\sum\limits_{j\neq
j^{\prime}}{\beta_{j}^{2}\beta_{j^{\prime}}^{2}}}}\right\\}+\frac{{{2p^{2}}{\sigma_{Y}^{4}}}}{{{n^{3}}}}~{}+~{}O\left({{n^{-2}}}\right).$
(2.56)

Lastly, plug-in (2.54) and (2.56) into (Supplementary Material) to get

$\displaystyle\text{Var}\left(T\right)$
$\displaystyle=\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-4{\text{Cov}}\left({{{\hat{\tau}}^{2}},\sum\limits_{j=1}^{p}{\sum\limits_{j^{\prime}=1}^{p}{{{\hat{\psi}}_{jj^{\prime}}}}}}\right)+4\text{Var}\left({\sum\limits_{j=1}^{p}{\sum\limits_{j^{\prime}=1}^{p}{{{\hat{\psi}}_{jj^{\prime}}}}}}\right)$
$\displaystyle=\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-4\left(\frac{2}{n}\left\\{{\sum\limits_{j=1}^{p}{\beta_{j}^{4}\left[{E\left({X_{1j}^{4}-1}\right)}\right]+2\sum\limits_{j\neq
j^{\prime}}{\beta_{j}^{2}\beta_{j^{\prime}}^{2}}}}\right\\}\right)$
$\displaystyle+4\left({\frac{1}{n}\left\\{{\sum\limits_{j=1}^{p}{\beta_{j}^{4}\left[{E\left({X_{1j}^{4}-1}\right)}\right]+2\sum\limits_{j\neq
j^{\prime}}{\beta_{j}^{2}\beta_{j^{\prime}}^{2}}}}\right\\}+\frac{{{2p^{2}}{\sigma_{Y}^{4}}}}{{{n^{3}}}}}\right)+O\left({{n^{-2}}}\right)$
$\displaystyle=\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-\frac{4}{n}\left\\{{\sum\limits_{j=1}^{p}{\beta_{j}^{4}\left[{E\left({X_{1j}^{4}-1}\right)}\right]+2\sum\limits_{j\neq
j^{\prime}}{\beta_{j}^{2}\beta_{j^{\prime}}^{2}}}}\right\\}+\frac{{{8p^{2}}{\sigma_{Y}^{4}}}}{{{n^{3}}}}+O({n^{-2}})$
$\displaystyle=\text{Var}\left({{T_{oracle}}}\right)+\frac{{4{p^{2}}{\sigma_{Y}^{4}}}}{{{n^{3}}}}+O\left({{n^{-2}}}\right),$

where the last equality holds by (2.12). ∎

######  Remark 2.5.

Calculations for equation 2.16:
Write,

$\displaystyle\text{Cov}\left({{{\hat{\tau}}^{2}},{g_{n}}}\right)$
$\displaystyle=\text{Cov}\left({\frac{1}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{{W_{{i_{1}}j}}{W_{{i_{2}}j}}}},\frac{1}{n}\sum\limits_{i=1}^{n}{{g_{i}}}}\right)$
$\displaystyle=\frac{1}{{{n^{2}}\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{\sum\limits_{i=1}^{n}{E\left({{W_{{i_{1}}j}}{W_{{i_{2}}j}}{g_{i}}}\right)}}}$
$\displaystyle=\frac{2}{{{n^{2}}\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{E\left({{W_{{i_{1}}j}}{g_{{i_{1}}}}}\right)E\left({{W_{{i_{2}}j}}}\right)}}$
$\displaystyle=\frac{2}{{{n^{2}}\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{E\left({{W_{{i_{1}}j}}{g_{{i_{1}}}}}\right){\beta_{j}}}}$
$\displaystyle=\frac{2}{n}\sum\limits_{j=1}^{p}{E\left({{S_{ij}}}\right){\beta_{j}}},$

where $S_{ij}\equiv W_{ij}g_{i}$. Also notice that
$\text{Var}\left({{g_{n}}}\right)=\text{Var}\left({\frac{1}{n}\sum\limits_{i=1}^{n}{{g_{i}}}}\right)=\frac{{\text{Var}\left({{g_{i}}}\right)}}{n}.$
Thus, by (2.8) we get

${c^{*}}=\frac{{{\text{Cov}}\left({{{\hat{\tau}}^{2}},{g_{n}}}\right)}}{{\text{Var}\left({{g_{n}}}\right)}}=\frac{{2\sum\limits_{j=1}^{p}{E\left({{S_{ij}}}\right){\beta_{j}}}}}{{\text{Var}\left({{g_{i}}}\right)}}.$

######  Remark 2.6.

Calculations for Example 2.3:
In order to calculate $\text{Var}(T_{c^{*}})$ we need to calculate the
numerator and denominator of (2.17). Consider first $\theta_{j}\equiv
E(S_{ij}).$ Write,

$\begin{array}[]{l}{\theta_{j}}\equiv
E\left({{S_{ij}}}\right)=E\left({{X_{ij}}{Y_{i}}{g_{i}}}\right)=E\left({{X_{ij}}\left({{\beta^{T}}{X_{i}}+\varepsilon_{i}}\right){g_{i}}}\right)=\\\
E\left({{X_{ij}}\left({\sum\limits_{m=1}^{p}{{\beta_{m}}{X_{im}}}+\varepsilon_{i}}\right)\sum\limits_{k<k^{\prime}}{{X_{ik}}{X_{ik^{\prime}}}}}\right)=\sum\limits_{m=1}^{p}{}\sum\limits_{k<k^{\prime}}{{\beta_{m}}E\left({{X_{ij}}{X_{im}}{X_{ik}}{X_{ik^{\prime}}}}\right)}\end{array}$

where in the last equality we used the assumption that $E(\varepsilon|X)=0$.
Since the columns of $\bf{X}$ are independent, the summation is not zero (up
to permutations) when $j=k$ and $m=k^{\prime}$. In this case we have

${\theta_{j}}=\sum\limits_{m=1}^{p}{\sum\limits_{k<k^{\prime}}{{\beta_{m}}E\left({{X_{ij}}{X_{im}}{X_{ik}}{X_{ik^{\prime}}}}\right)}}=\sum\limits_{m\neq
j}^{p}{{\beta_{m}}E\left({X_{ij}^{2}X_{im}^{2}}\right)}=\sum\limits_{m\neq
j}^{p}{{\beta_{m}}E\left({X_{ij}^{2}}\right)E\left({X_{im}^{2}}\right)}=\sum\limits_{m\neq
j}^{p}{{\beta_{m}}}.$

Notice that in the fourth equality we used the assumption that
$E(X_{ij}^{2})=1$ for all $j=1,...,p.$ Thus,

$\sum\limits_{j=1}^{p}{{\beta_{j}}E\left({{S_{ij}}}\right)=}\sum\limits_{j=1}^{p}{{\beta_{j}}\sum\limits_{m\neq
j}^{p}{{\beta_{m}}}=}\sum\limits_{j=1}^{p}{{\beta_{j}}\left({\sum\limits_{m=1}^{p}{{\beta_{m}}}-{\beta_{j}}}\right)={{\left({\sum\limits_{j=1}^{p}{{\beta_{j}}}}\right)}^{2}}-\sum\limits_{j=1}^{p}{\beta_{j}^{2}}=}{\left({\sum\limits_{j=1}^{p}{{\beta_{j}}}}\right)^{2}}-{\tau^{2}}.$
(2.57)

plug-in $\tau^{2}=1$ and $\beta_{j}=\frac{1}{\sqrt{p}}$ to get the numerator
of (2.17):

${\left[{2\sum\limits_{j=1}^{p}{{\beta_{j}}}E\left({{S_{ij}}}\right)}\right]^{2}}=4{\left[{{{\left({\sum\limits_{j=1}^{p}{{\beta_{j}}}}\right)}^{2}}-{\tau^{2}}}\right]^{2}}=4{\left[{{{\left({p\frac{1}{{\sqrt{p}}}}\right)}^{2}}-1}\right]^{2}}=4{\left({{p}-1}\right)^{2}}.$

Consider now the denominator of (2.17). Write,

$\text{Var}\left({{g_{i}}}\right)=E\left({g_{i}^{2}}\right)=E\left[{{{\left({\sum\limits_{j<j^{\prime}}{{X_{ij}}{X_{ij^{\prime}}}}}\right)}^{2}}}\right]=\sum\limits_{{j_{1}}<{j_{2}}}{\sum\limits_{{j_{3}}<{j_{4}}}{E\left({{X_{i{j_{1}}}}{X_{i{j_{2}}}}{X_{i{j_{3}}}}{X_{i{j_{4}}}}}\right)}}.$

Since we assume that the columns of X are independent, the summation is not
zero when ${j_{1}}={j_{3}}$ and ${j_{2}}={j_{4}}.$ Thus,

$\text{Var}\left({{g_{i}}}\right)=\sum\limits_{{j_{1}}<{j_{2}}}{E\left({X_{i{j_{1}}}^{2}X_{i{j_{2}}}^{2}}\right)=\sum\limits_{{j_{1}}<{j_{2}}}{E\left({X_{i{j_{1}}}^{2}}\right)E\left({X_{i{j_{2}}}^{2}}\right)}}=p\left({p-1}\right)/2.$
(2.58)

Notice that we used the assumption that since we assume that
${\bf{\Sigma}}={\bf{I}}$ in the last equality. Now, recall by (2.49) that
$\text{Var}\left({{{\hat{\tau}}^{2}}}\right)=\frac{{20}}{n}+O\left({\frac{1}{{{n^{2}}}}}\right).$
Therefore, we have

$\text{Var}\left({{T_{{c^{*}}}}}\right)=\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-\frac{{{{\left[{2\sum\limits_{j=1}^{p}{{\beta_{j}}E\left({{S_{ij}}}\right)}}\right]}^{2}}}}{{n\text{Var}\left({{g_{i}}}\right)}}=\frac{{20}}{n}+O\left({\frac{1}{{{n^{2}}}}}\right)-\frac{{4{{\left({p-1}\right)}^{2}}}}{{n\cdot\left[{p\left({p-1}\right)/2}\right]}}=\frac{{12}}{n}+O\left({\frac{1}{{{n^{2}}}}}\right),$
(2.59)

where we used the assumption that $n=p$ in the last equality.

Proof of Proposition 2.4:
We need to prove that
$\sqrt{n}\left[{{T_{{c^{*}}}}-{{T}_{{\hat{c}^{*}}}}}\right]\overset{p}{\rightarrow}0.$
Write,

$\sqrt{n}\left[{{T_{c^{*}}}-T_{\hat{c}^{*}}}\right]=\sqrt{n}\left[{{{\hat{\tau}}^{2}}-{c^{*}}{g_{n}}-\left({{{\hat{\tau}}^{2}}-{{\hat{c}}^{*}}{g_{n}}}\right)}\right]={\sqrt{n}{g_{n}}}{\left({{{\hat{c}}^{*}}-{c^{*}}}\right)}.$

By Markov and Cauchy-Schwarz inequalities, it is enough to show that

$\displaystyle
P\left\\{{\left|{\sqrt{n}{g_{n}}\left({{\hat{c}^{*}}-{{c}^{*}}}\right)}\right|>\varepsilon}\right\\}$
$\displaystyle\leq\frac{{E\left\\{{\left|{\sqrt{n}{g_{n}}\left({{\hat{c}^{*}}-{{c}^{*}}}\right)}\right|}\right\\}}}{\varepsilon}\leq\frac{{\sqrt{nE\left({g_{n}^{2}}\right)E\left[{{{\left({{\hat{c}^{*}}-{{c}^{*}}}\right)}^{2}}}\right]}}}{\varepsilon}$
$\displaystyle=\frac{\sqrt{\text{Var}(g)\text{Var}(\hat{c}^{*})}}{\varepsilon}\rightarrow
0.$ (2.60)

Notice that by (2.18) we have

$\text{Var}\left({{g}}\right)\text{Var}\left({{{\hat{c}}^{*}}}\right)=\frac{{\text{Var}\left(U\right)}}{{\text{Var}\left({{g}}\right)}},$
(2.61)

where
$U\equiv\frac{2}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{}\sum\limits_{j=1}^{p}{{W_{{i_{1}}j}}W_{i_{2}j}g_{i_{2}}}.$
The variance of $g$ is

$\displaystyle\text{Var}({{g}})$
$\displaystyle=E({g^{2}})=E[{{{({\sum\limits_{j<j^{\prime}}{{X_{ij}}{X_{ij^{\prime}}}}})^{2}}}}]=\sum\limits_{{j_{1}}<{j_{2}}}{\sum\limits_{{j_{3}}<{j_{4}}}{E\left({{X_{i{j_{1}}}}{X_{i{j_{2}}}}{X_{i{j_{3}}}}{X_{i{j_{4}}}}}\right)}}$
$\displaystyle=\sum\limits_{{j_{1}}<{j_{2}}}{E\left({X_{i{j_{1}}}^{2}X_{i{j_{2}}}^{2}}\right)=\sum\limits_{{j_{1}}<{j_{2}}}{E\left({X_{i{j_{1}}}^{2}}\right)E\left({X_{i{j_{2}}}^{2}}\right)}}=p\left({p-1}\right)/2,$
(2.62)

where the equation above holds since we assume that ${\bf{\Sigma}}={\bf{I}}$
and that the columns of X are independent. Hence, by (Supplementary Material)
- (Supplementary Material) it enough to prove
$\frac{\text{Var}(U)}{p^{2}}\rightarrow 0.$

The variance of $U$ is

$\displaystyle\text{Var}(U)$
$\displaystyle=\text{Var}\left[{\frac{2}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{{W_{{i_{1}}j}}{W_{{i_{2}}j}}g_{i_{2}}}}}\right]$
$\displaystyle=\frac{4}{{{n^{2}}{{\left({n-1}\right)}^{2}}}}\sum\limits_{j,j^{\prime}}^{p}{\sum\limits_{{i_{1}}\neq{i_{2}},{i_{3}}\neq{i_{4}}}{\text{Cov}\left[{{W_{{i_{1}}j}}{W_{{i_{2}}j}}g_{i_{2}},{W_{{i_{3}}j^{\prime}.}}{W_{{i_{4}}j^{\prime}}}g_{i_{4}}}\right]}}$
(2.63)

The covariance in (Supplementary Material) is different from zero in the two
following cases:

  1. 1.

When $\left\\{{{i_{1}},{i_{2}}}\right\\}$ equals to
$\left\\{{{i_{3}},{i_{4}}}\right\\}.$

  2. 2.

When one of $\left\\{{{i_{1}},{i_{2}}}\right\\}$ equals to
$\left\\{{{i_{3}},{i_{4}}}\right\\}$ while the other is different.

The first condition includes two different sub-cases and each of those
consists $n(n-1)$ quadruples $(i_{1},i_{2},i_{3},i_{4})$ that satisfy the
condition. Similarly, the second condition above includes four different sub-
cases and each of those consists $n(n-1)(n-2)$ quadruples that satisfy the
condition.

We now calculate the covariance for all these six sub-cases.
(1) The covariance when ${i_{1}}={i_{3}},{i_{2}}={i_{4}}$ is

$\displaystyle{\delta_{1}}$
$\displaystyle\equiv\text{Cov}\left[{{W_{j}}{{\tilde{W}}_{j}}\tilde{g},{W_{j^{\prime}}}{{\tilde{W}_{j^{\prime}}}}\tilde{g}}\right]=E\left({{W_{j}}{W_{j^{\prime}}}}\right)E\left[{{{\tilde{W}}_{j}}{{\tilde{W}_{j^{\prime}}}}\tilde{g}^{2}}\right]-E\left({{W_{j}}}\right)E\left[{{{\tilde{W}}_{j}}\tilde{g}}\right]E\left({{{\tilde{W}_{j^{\prime}}}}}\right)E\left[{{{\tilde{W}_{j^{\prime}}}}\tilde{g}}\right]$
$\displaystyle=E\left({{W_{j}}{W_{j^{\prime}}}}\right)E\left[{{{\tilde{W}}_{j}}{{\tilde{W}_{j^{\prime}}}}\tilde{g}^{2}}\right]-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{j}}{\theta_{j^{\prime}}},$
(2.64)

where $\tilde{W}$ and $\tilde{g}$ are independent copies of $W$ and $g$
respectively.
(2) The covariance when and ${i_{1}}={i_{4}},{i_{2}}={i_{3}}$ is

$\displaystyle{\delta_{2}}$
$\displaystyle\equiv\text{Cov}\left[{{W_{j}}{{\tilde{W}}_{j}}\tilde{g},{{\tilde{W}_{j^{\prime}}}}{W_{j^{\prime}}}g}\right]=E\left[{{W_{j}}{W_{j^{\prime}}}g}\right]E\left[{{{\tilde{W}}_{j}}{{\tilde{W}_{j^{\prime}}}}\tilde{g}}\right]-E\left(W\right)E\left[{{{\tilde{W}}_{j}}\tilde{g}}\right]E\left({{{\tilde{W}_{j^{\prime}}}}}\right)E\left[{{W_{j^{\prime}}}g}\right]$
$\displaystyle={\left\\{{E\left[{{W_{j}}{W_{j^{\prime}}}g}\right]}\right\\}^{2}}-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{j}}{\theta_{j^{\prime}}}.$
(2.65)

(3) The covariance when $\,{i_{1}}={i_{3}},{i_{2}}\neq{i_{4}}$ is

${\delta_{3}}\equiv\text{Cov}\left[{{W_{j}}{{\tilde{W}}_{j}}\tilde{g},{W_{j^{\prime}}}{{\tilde{\tilde{W}}}_{j^{\prime}}}\tilde{\tilde{g}}}\right]=E\left({{W_{j}}{W_{j^{\prime}}}}\right){\theta_{j}}{\theta_{j^{\prime}}}-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{j}}{\theta_{j^{\prime}}},$
(2.66)

where $\tilde{\tilde{W}}$ and $\tilde{\tilde{g}}$ are another independent
copies of $W$ and $g$ respectively.
(4) The covariance when $i_{1}=i_{4},i_{2}\neq i_{3}$ is

${\delta_{4}}\equiv\text{Cov}\left[{{W_{j}}{{\tilde{W}}_{j}}\tilde{g},{{\tilde{\tilde{W}}}_{j^{\prime}}}{W_{j^{\prime}}}g}\right]={\theta_{j}}{\beta_{j^{\prime}}}E\left[{{W_{j}}{W_{j^{\prime}}}g}\right]-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{j}}{\theta_{j^{\prime}}}.$
(2.67)

(5) The covariance when $i_{2}=i_{3},i_{1}\neq i_{4}$ is similar to
$\delta_{4}$, i.e.,

${\delta_{5}}\equiv\text{Cov}\left[{{W_{j}}{{\tilde{W}}_{j}}\tilde{g},{{\tilde{W}_{j^{\prime}}}}{{\tilde{\tilde{W}}}_{j^{\prime}}}\tilde{\tilde{g}}}\right]={\beta_{j}}{\theta_{j^{\prime}}}E\left[{{{\tilde{W}}_{j}}{{\tilde{W}_{j^{\prime}}}}\tilde{g}}\right]-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{j}}{\theta_{j^{\prime}}}.$
(2.68)

(6) The covariance when $i_{2}=i_{4},i_{1}\neq i_{3}$ is

${\delta_{6}}\equiv\text{Cov}\left[{{W_{j}}{{\tilde{W}}_{j}}\tilde{g},{{\tilde{\tilde{W}}}_{j^{\prime}}}{{\tilde{W}_{j^{\prime}}}}\tilde{g}}\right]={\beta_{j}}{\beta_{j^{\prime}}}E\left[{{{\tilde{W}}_{j}}{{\tilde{W}_{j^{\prime}}}}\tilde{g}^{2}}\right]-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{j}}{\theta_{j^{\prime}}}$
(2.69)

Thus, plugging-in (Supplementary Material) - (2.69) into (Supplementary
Material) gives

$\text{Var}\left(U\right)=4\sum\limits_{j,j^{\prime}}{\left\\{{\frac{1}{{n\left({n-1}\right)}}\left({{\delta_{1}}+{\delta_{2}}}\right)+\frac{{\left({n-2}\right)}}{{n\left({n-1}\right)}}\left({{\delta_{3}}+{\delta_{4}}+{\delta_{5}}+{\delta_{6}}}\right)}\right\\}}.$
(2.70)

Recall that we wish to show that $\frac{\text{Var}(U)}{p^{2}}\rightarrow 0.$
Since we assume that $n/p=O(1),$ it is enough to show that
$\sum\limits_{j,j^{\prime}}{\frac{{\left({{\delta_{1}}+{\delta_{2}}}\right)}}{{{n^{4}}}}}\to
0$ and
$\sum\limits_{j,j^{\prime}}{\frac{{\left({{\delta_{3}}+{\delta_{4}}+{\delta_{5}}+{\delta_{6}}}\right)}}{{{n^{3}}}}}\to
0$. Careful calculations, which are not presented here, show that under the
linear model when the covariates are independent

$\displaystyle\sum_{j,j^{\prime}}\delta_{1}$ $\displaystyle\leq
C^{2}\tau^{2}\sigma^{2}p^{2}\frac{p-1}{2},~{}\sum_{j,j^{\prime}}\delta_{2}\leq
C^{2}\tau^{3}p^{3}\quad\text{and},$
$\displaystyle\sum_{j,j^{\prime}}\delta_{3}$ $\displaystyle\leq
p^{2}\tau^{2},\sum_{j,j^{\prime}}\delta_{4}\leq
C\tau^{4}p^{2},~{}\sum_{j,j^{\prime}}\delta_{5}\leq
C\tau^{2}(\tau^{2}+\sigma^{2})\frac{p(p-1)}{2},$

where $C$ is a bound on $E(X_{ij}^{8})$ for all $j.$ It follows that
$\frac{\text{Var}(U)}{p^{2}}\to 0$ because $n/p=O(1)$, which completes the
proof of the proposition.

######  Remark 2.7.

We now calculate the asymptotic improvement of $T_{\bf{B}}$ over the naive
estimator. For simplicity, consider the case when $\tau^{2}=\sigma^{2}=1.$
Recall the variance of $\hat{\tau}^{2}$ and $T_{\bf{B}}$ given in (2.6) and
(2.21), respectively. Write,

$\displaystyle\mathop{\lim}\limits_{n,p\to\infty}\frac{{\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-\text{Var}\left({{T_{\bf{B}}}}\right)}}{{\text{Var}\left({{{\hat{\tau}}^{2}}}\right)}}$
$\displaystyle=\mathop{\lim}\limits_{n,p\to\infty}\frac{{8{\tau_{\bf{B}}^{4}}/n}}{{\frac{4}{n}\left[{\frac{{\left({n-2}\right)}}{{\left({n-1}\right)}}\left({\sigma_{Y}^{2}{\tau^{2}}+{\tau^{4}}}\right)+\frac{1}{{2\left({n-1}\right)}}\left({p{\sigma_{Y}^{4}}+4\sigma_{Y}^{2}{\tau^{2}}+3{\tau^{4}}}\right)}\right]}}$
$\displaystyle=\frac{{2{\tau_{\bf{B}}^{4}}}}{{3{\tau^{4}}+\frac{{4p{\tau^{4}}+4\sigma_{Y}^{2}{\tau^{2}}+3{\tau^{4}}}}{{2n}}}}=\frac{0.5}{{3+2\frac{p}{n}}},$

where we used (2.6) in the first equality, and the fact that
$\sigma_{Y}^{2}=2\tau^{2}=2$ in the second equality. Now, notice that when
$p=n$ and $\tau_{\bf{B}}^{2}=0.5$ then the reduction is
$\frac{0.5}{{3+2}}=10\%$ and when $p/n$ converges to zero, the reduction is
$16\%.$

######  Remark 2.8.

Calculations for Example 2.4:
Consider the first scenario where $\beta_{j}^{2}=\frac{1}{p}.$ Recall that we
assume that the set B is a fixed set of indices such that
$\left|{\bf{B}}\right|\ll p$. Therefore, we have
$\tau_{\bf{B}}^{2}=\sum\limits_{j\in{\bf{B}}}{\beta_{j}^{2}=O\left({\frac{1}{p}}\right)}$.
Now, by (2.22) we have
$\text{Var}(T_{\textbf{B}})=\text{Var}(\hat{\tau}^{2})-\frac{8}{n}\tau^{2}_{\textbf{B}}+O(n^{-2})$
and by Remark 2.6 we have $\text{Var}(\hat{\tau}^{2})=\frac{20}{n}+O(n^{-2}).$
Using the assumption that $n=p$ we can conclude that
$\text{Var}\left({{T_{\bf{B}}}}\right)=\frac{{20}}{n}+O\left({\frac{1}{{{n^{2}}}}}\right).$
Hence, in this scenario, $T_{\textbf{B}}$ and the naive estimator have the
same asymptotic variance. In contrast, recall that in Example 2.3 we showed
that the asymptotic variance of $T_{c^{*}}$ is $40\%$ lower than the variance
of the naive estimator.

Consider now the second scenario where $\hat{\tau}_{\bf{B}}^{2}={\tau^{2}}=1$.
By (2.22) we have

$\text{Var}\left({{T_{\bf{B}}}}\right)={\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-\frac{8}{n}\tau_{\bf{B}}^{4}}+O(n^{-2})=\frac{12}{n}+O(n^{-2}).$

Hence, in this scenario the asymptotic variance of $T_{\textbf{B}}$ is $40\%$
smaller than the variance of the naive estimator. Consider now
$\text{Var}(T_{c^{*}}).$ By Cauchy–Schwarz inequality
${\left({\sum\limits_{j\in{\bf{B}}}{{\beta_{j}}}}\right)^{2}}\leq\sum\limits_{j\in{\bf{B}}}{\beta_{j}^{2}\cdot\left|{\bf{B}}\right|}=\tau^{2}_{\textbf{B}}\left|{\bf{B}}\right|=O(1)$,
where the last equality holds since we assume that
${\bf{B}}\subset\left\\{{1,...,p}\right\\}$ be a fixed set of some indices
such that $\left|{\bf{B}}\right|\ll p.$ Now, By (2.57) we have

$\sum\limits_{j=1}^{p}{{\beta_{j}}}{\theta_{j}}={\left({\sum\limits_{j=1}^{p}{{\beta_{j}}}}\right)^{2}}-{\tau^{2}}={\left({\sum\limits_{j\in{\bf{B}}}^{p}{{\beta_{j}}}+\overbrace{\sum\limits_{j\notin{\bf{B}}}^{p}{{\beta_{j}}}}^{0}}\right)^{2}}-\tau_{\bf{B}}^{2}\leq\left|{\bf{B}}\right|-1=O\left(1\right).$

Now, recall that
$\text{Var}\left({{{\hat{\tau}}^{2}}}\right)=\frac{{20}}{n}+O\left({\frac{1}{{{n^{2}}}}}\right)$
and $\text{Var}(g_{i})=p(p-1)/2$ by (2.49) and (2.58) respectively. Therefore,
we have

$\text{Var}\left({{T_{{c^{*}}}}}\right)=\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-\frac{{{{\left[{2\sum\limits_{j=1}^{p}{{\beta_{j}}\theta_{j}}}\right]}^{2}}}}{{n\text{Var}\left({{g_{i}}}\right)}}=\frac{{20}}{n}+O\left({\frac{1}{{{n^{2}}}}}\right)-O\left(\frac{1}{np^{2}}\right)=\frac{{20}}{n}+O\left({\frac{1}{{{n^{2}}}}}\right),$

Hence, in this scenario, $T_{c^{*}}$ and the naive estimator have the same
asymptotic variance.

Lastly, recall that in Example 2.2 we already showed that, asymptotically, the
variance of $T_{oracle}$ (i.e., the optimal oracle estimator) is $40\%$ lower
than the naive variance (without any assumptions about the structure of the
coefficient vector $\beta$).

Proof of Proposition 2.5:
In order to prove that
$\sqrt{n}\left({{T_{\gamma}}-{T_{\bf{B}}}}\right)\overset{p}{\rightarrow}0,$
it is enough to show that

$\displaystyle
E\left\\{{\sqrt{n}\left({{T_{\gamma}}-{T_{\bf{B}}}}\right)}\right\\}\to 0,$
(2.71)
$\displaystyle\text{Var}\left\\{{\sqrt{n}\left({{T_{\gamma}}-{T_{\bf{B}}}}\right)}\right\\}\to
0.$ (2.72)

We start with the first equation. Let $A$ denote the event that the selection
algorithm $\gamma$ perfectly identifies the set of large coefficients, i.e.,
$A=\left\\{{{{\bf{B}}_{\gamma}}={\bf{B}}}\right\\}.$ Let $p_{A}\equiv P(A)$
denote the probability that $A$ occurs, and let $\mathbbm{1}_{A}$ denote the
indicator of $A.$ Notice that $E(T_{\textbf{B}})=\tau^{2}$ and
${T_{\gamma}}\mathbbm{1}_{A}={T_{\rm B}}\mathbbm{1}_{A}.$ Thus,

$\displaystyle
E\left\\{{\sqrt{n}\left({{T_{\gamma}}-{T_{\bf{B}}}}\right)}\right\\}=\sqrt{n}\left[{E\left({{T_{\gamma}}}\right)-{\tau^{2}}}\right]$
$\displaystyle=\sqrt{n}\left(E\left[{{T_{\gamma}}(1-\mathbbm{1}_{A})}\right]+E\left[{{T_{\gamma}}\mathbbm{1}_{A}}\right]-{\tau^{2}}\right)$
$\displaystyle=\sqrt{n}E\left[{{T_{\gamma}}(1-\mathbbm{1}_{A})}\right]+\sqrt{n}\left[E\left({{T_{\bf{B}}}\mathbbm{1}_{A}}\right)-{\tau^{2}}\right],$
(2.73)

where the last equality holds since ${T_{\gamma}}\mathbbm{1}_{A}={T_{\rm
B}}\mathbbm{1}_{A}.$ For the convenience of notation, let $C$ be an upper
bound of the maximum over all first four moments of ${T_{\gamma}}$ and
$T_{\bf{B}},$ and consider the first term of (2.73). By the Cauchy–Schwarz
inequality,

$\sqrt{n}E\left[{{T_{\gamma}}\left({1-\mathbbm{1}_{A})}\right)}\right]\leq\sqrt{n}{\left\\{{E\left[{T_{\gamma}^{2}}\right]}\right\\}^{1/2}}{\left\\{{E\left[{{{\left({1-\mathbbm{1}_{A}}\right)}^{2}}}\right]}\right\\}^{1/2}}\leq\sqrt{n}{C^{1/2}}{\left\\{{1-{p_{A}}}\right\\}^{1/2}}\underset{n\to\infty}{\rightarrow}0,$
(2.74)

where the last inequality holds since
$\mathop{\lim}\limits_{n\to\infty}n\left({1-{p_{A}}}\right)^{1/2}=0$ by
assumption. We now consider the second term of (2.73). Write,

$\sqrt{n}\left[{E\left({{T_{\bf{B}}}\mathbbm{1}_{A}}\right)-{\tau^{2}}}\right]=\sqrt{n}E\left({{T_{\bf{B}}}\mathbbm{1}_{A}-{T_{\bf{B}}}}\right)=-\sqrt{n}E\left[{{T_{\bf{B}}}\left({1-\mathbbm{1}_{A}}\right)}\right],$
(2.75)

and notice that by the same type of argument as in (2.74) we have
$\sqrt{n}E\left[{{T_{\bf{B}}}\left({1-\mathbbm{1}_{A}}\right)}\right]\underset{n\to\infty}{\rightarrow}~{}0.$
This completes the proof of (2.71).

We now move to show that
$\text{Var}\left\\{{\sqrt{n}\left({{T_{\gamma}}-{T_{\bf{B}}}}\right)}\right\\}\to
0$. Write,

$\displaystyle\text{Var}\left\\{{\sqrt{n}\left({{T_{\gamma}}-{T_{\bf{B}}}}\right)}\right\\}=n\text{Var}\left({{T_{\gamma}}-{T_{\bf{B}}}}\right)$
$\displaystyle=n\left[{\text{Var}\left({{T_{\gamma}}}\right)+\text{Var}\left({{T_{\bf{B}}}}\right)-2{\text{Cov}}\left({{T_{\gamma}},{T_{\bf{B}}}}\right)}\right]$
$\displaystyle=n\left\\{{E\left({T_{\gamma}^{2}}\right)-{{\left[{E\left({{T_{\gamma}}}\right)}\right]}^{2}}+E\left({T_{\bf{B}}^{2}}\right)-{\tau^{4}}-2\left[{E\left({{T_{\gamma}}{T_{\bf{B}}}}\right)-E\left({{T_{\gamma}}}\right){\tau^{2}}}\right]}\right\\}$
$\displaystyle=n\left\\{{E\left({T_{\gamma}^{2}}\right)-E\left({{T_{\gamma}}{T_{\bf{B}}}}\right)+E\left({T_{\bf{B}}^{2}}\right)-E\left({{T_{\gamma}}{T_{\bf{B}}}}\right)+E\left({{T_{\gamma}}}\right)\left[{{\tau^{2}}-E\left({{T_{\gamma}}}\right)}\right]-{\tau^{2}}\left[{{\tau^{2}}-E\left({{T_{\gamma}}}\right)}\right]}\right\\}$
$\displaystyle=n\left\\{{\underbrace{E\left({T_{\gamma}^{2}}\right)-E\left({{T_{\gamma}}{T_{\bf{B}}}}\right)}_{{\theta_{1}}}+\underbrace{E\left({T_{\bf{B}}^{2}}\right)-E\left({{T_{\gamma}}{T_{\bf{B}}}}\right)}_{{\theta_{2}}}-\underbrace{{{\left[{{\tau^{2}}-E\left({{T_{\gamma}}}\right)}\right]}}}_{{\theta_{3}}}}^{2}\right\\}.$

Thus, it is enough to show that $n\theta_{1}\rightarrow 0$,
$n\theta_{2}\rightarrow 0$ and $n\theta_{3}\rightarrow 0.$

We start with showing that $n\theta_{1}\rightarrow 0.$ Notice that
$T_{\textbf{B}}^{2}\mathbbm{1}_{A}=T_{\textbf{B}}T_{\gamma}\mathbbm{1}_{A}=T_{\gamma}^{2}\mathbbm{1}_{A}$
Thus,

$\displaystyle n{\theta_{1}}$
$\displaystyle=n\left\\{{E\left({T_{\gamma}^{2}}\right)-E\left({{T_{\gamma}}{T_{\bf{B}}}}\right)}\right\\}$
$\displaystyle=n\left\\{{E\left({T_{\gamma}^{2}}\right)-E\left[{{T_{\gamma}}{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]-E\left({{T_{\gamma}}{T_{\bf{B}}}{\mathbbm{1}_{A}}}\right)}\right\\}$
$\displaystyle=n\left\\{{E\left({T_{\gamma}^{2}}\right)-E\left[{{T_{\gamma}}{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]-E\left({T_{\gamma}^{2}{\mathbbm{1}_{A}}}\right)}\right\\}$
$\displaystyle=n\left\\{{E\left[{T_{\gamma}^{2}\left({1-{\mathbbm{1}_{A}}}\right)}\right]-E\left[{{T_{\gamma}}{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]}\right\\}.$

Now, notice that
$n(E\left[{T_{\gamma}^{2}\left({1-{\mathbbm{1}_{A}}}\right)}\right])\rightarrow
0$ by similar arguments as in (2.74), with a slight modification of using the
existence of the fourth moments of $T_{\gamma}$ and $T_{\bf{B}}$, rather than
the second moments. Also, by Cauchy–Schwarz inequality we have,

$\displaystyle
nE\left[{{T_{\gamma}}{T_{\bf{B}}}\left({1-{1_{A}}}\right)}\right]$
$\displaystyle\leq
n{\left\\{{E\left({T_{\gamma}^{2}T_{\bf{B}}^{2}}\right)}\right\\}^{1/2}}{\left\\{{E\left[{{{\left({1-{1_{A}}}\right)}^{2}}}\right]}\right\\}^{1/2}}$
$\displaystyle\leq
n{\left\\{{E\left({T_{\gamma}^{4}}\right)E\left({T_{\bf{B}}^{4}}\right)}\right\\}^{1/4}}{\left\\{{1-{p_{A}}}\right\\}^{1/2}}$
$\displaystyle\leq n{C^{1/2}}{\left\\{{1-{p_{A}}}\right\\}^{1/2}}\to 0,$

where $C$ is an upper bound of the maximum over all first four moments of
${T_{\gamma}}$ and $T_{\bf{B}}.$ Therefore, $n\theta_{1}\rightarrow 0.$

Consider now $n\theta_{2}.$ Write,

$\displaystyle n{\theta_{2}}$
$\displaystyle=n\left\\{{E\left({T_{\bf{B}}^{2}}\right)-E\left({{T_{\gamma}}{T_{\bf{B}}}}\right)}\right\\}$
$\displaystyle=n\left\\{{E\left({T_{\bf{B}}^{2}}\right)-E\left[{{T_{\gamma}}{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]-E\left({{T_{\gamma}}{T_{\bf{B}}}{\mathbbm{1}_{A}}}\right)}\right\\}$
$\displaystyle=n\left\\{{E\left({T_{\bf{B}}^{2}}\right)-E\left[{{T_{\gamma}}{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]-E\left({T_{\bf{B}}^{2}{\mathbbm{1}_{A}}}\right)}\right\\}$
$\displaystyle=n\left\\{{E\left[{T_{\bf{B}}^{2}\left({1-{\mathbbm{1}_{A}}}\right)}\right]-E\left[{{T_{\gamma}}{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]}\right\\}\rightarrow
0,$

and notice that the last equation follows by similar arguments.

Consider now $n\theta_{3}.$ Write,

$\displaystyle n\theta_{3}$
$\displaystyle=n\left[{\tau^{2}-E\left({{T_{\gamma}}}\right)}\right]$
$\displaystyle=n\left[{E\left({{T_{\bf{B}}}}\right)-E\left({{T_{\gamma}}}\right)}\right]$
$\displaystyle=n\left[{E\left[{{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)+{T_{\bf{B}}}{\mathbbm{1}_{A}}}\right]-E\left({{T_{\gamma}}}\right)}\right]$
$\displaystyle=n\left\\{{E\left[{{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]+E\left({{T_{\bf{B}}}{\mathbbm{1}_{A}}-{T_{\gamma}}}\right)}\right\\}$
$\displaystyle=n\left\\{{E\left[{{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]+E\left({{T_{\gamma}}{\mathbbm{1}_{A}}-{T_{\gamma}}}\right)}\right\\}$
$\displaystyle=n\left\\{{E\left[{{T_{\bf{B}}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]-E\left[{{T_{\gamma}}\left({1-{\mathbbm{1}_{A}}}\right)}\right]}\right\\}\rightarrow
0,$

where the last equation follows by similar arguments as in (2.74). This
completes the proof of (2.72) and we conclude that
$\sqrt{n}\left({{T_{\gamma}}-{T_{\bf{B}}}}\right)\overset{p}{\rightarrow}0.$ ∎

Proof of Proposition 2.6:
We wish to prove that

$\displaystyle
n\left[{\widehat{\text{Var}\left({{{\hat{\tau}}^{2}}}\right)}-\text{Var}\left({{\hat{\tau}^{2}}}\right)}\right]\overset{p}{\rightarrow}0.$
(2.76)

Recall by (2.6) that

${\text{Var}}\left({{{\hat{\tau}}^{2}}}\right)=\frac{{4\left({n-2}\right)}}{{n\left({n-1}\right)}}\left[{{\beta^{T}}{\bf{A}}\beta-{{\left\|\beta\right\|}^{4}}}\right]+\frac{2}{{n\left({n-1}\right)}}\left[{\left\|{\bf{A}}\right\|_{F}^{2}-{{\left\|\beta\right\|}^{4}}}\right].$

Now, when we assume standard Gaussian covariates, one can verify that
${\beta^{T}}{\bf{A}}\beta-{\left\|\beta\right\|^{4}}=\sigma_{Y}^{2}{\tau^{2}}+{\tau^{4}}$
and
$\left\|{\bf{A}}\right\|_{F}^{2}-{\left\|\beta\right\|^{4}}=p{\sigma_{Y}^{4}}+4\sigma_{Y}^{2}{\tau^{2}}+3{\tau^{4}},$
where $\sigma_{Y}^{2}=\sigma^{2}+\tau^{2}.$ Thus, in this case we can write

$\text{Var}\left({{{\hat{\tau}}^{2}}}\right)=\frac{4}{n}\left[{\frac{{\left({n-2}\right)}}{{\left({n-1}\right)}}\left[{\sigma_{Y}^{2}{\tau^{2}}+{\tau^{4}}}\right]+\frac{1}{{2\left({n-1}\right)}}\left({p{\sigma_{Y}^{4}}+4\sigma_{Y}^{2}{\tau^{2}}+3{\tau^{4}}}\right)}\right].$
(2.77)

In order to prove that (2.76) holds, it is enough to prove the consistency of
$\hat{\tau}^{2}$ and $\hat{\sigma}_{Y}^{2}.$ Consistency of the sample
variance $\hat{\sigma}_{Y}^{2}$ is a standard result, and since
$\hat{\tau}^{2}$ is an unbiased estimator, it is enough to show that its
variance converges to zero as $n\to\infty.$ Since we assume
$\hat{\tau}^{2}+\sigma^{2}=O(1)$ and $p/n=O(1),$ we have by (2.48) that
$\text{Var}(\hat{\tau}^{2})\underset{n\to\infty}{\rightarrow}0$, and (2.76)
follows.

Proof of Proposition 2.7:
We now move to prove that

$n\left[{\widehat{\text{Var}\left({T_{\gamma}}\right)}-\text{Var}\left({{T_{\gamma}}}\right)}\right]\overset{p}{\rightarrow}0,$
(2.78)

Recall that by Proposition 2.5 we have
$\mathop{\lim}\limits_{n\to\infty}n\left[{\text{Var}\left({{T_{\bf{B}}}}\right)-{\text{Var}}\left({{T_{\gamma}}}\right)}\right]=0.$
Hence, it is enough to show that

$n\left[{\widehat{\text{Var}\left({T_{\gamma}}\right)}-\text{Var}\left(T_{\bf{B}}\right)}\right]\overset{p}{\rightarrow}~{}0.$

Since we assume
${X_{i}}\mathop{\sim}\limits^{i.i.d}N\left({\bf{0},\bf{I}}\right)$ then by
(2.21) we have
$\text{Var}\left({{T_{\bf{B}}}}\right)={\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-\frac{8}{n}\tau_{\bf{B}}^{4}}+O(n^{-2}).$
Recall that by definition we have
$\widehat{\text{Var}\left({{T_{\gamma}}}\right)}={\widehat{\text{Var}\left({{{\hat{\tau}}^{2}}}\right)}-\frac{8}{n}\hat{\tau}_{{{\bf{B}}_{\gamma}}}^{4}}.$
Also recall that $\widehat{\text{Var}({\hat{\tau}}^{2}})$ is consistent by
Proposition 2.6. Thus, it is enough to prove that
$\hat{\tau}_{{{\bf{B}}_{\gamma}}}^{2}-\tau_{\bf{B}}^{2}\overset{p}{\rightarrow}0.$
Now, since we assumed that
$n\left[{P\left(\left\\{{{{\bf{B}}_{\gamma}}\neq{\bf{B}}}\right\\}\right)}\right]^{1/2}\xrightarrow[n\rightarrow\infty]{}0$
then clearly
$P\left({{{\bf{B}}_{\gamma}}={\bf{B}}}\right)\xrightarrow[n\rightarrow\infty]{}1.$
Thus, it is enough to show that
$\hat{\tau}_{{{\bf{B}}}}^{2}-\tau_{\bf{B}}^{2}\overset{p}{\rightarrow}0.$
Recall that $E(\hat{\beta}_{j}^{2})=\beta_{j}^{2}$ for $j=1,...,p$ and notice
that $\text{Var}(\hat{\beta}_{j}^{2})\underset{n\to\infty}{\rightarrow}0$ by
similar arguments that were used to derive (2.6). Hence, we have
$\hat{\beta}_{j}^{2}-\beta_{j}^{2}\overset{p}{\rightarrow}0.$ Since we assumed
that ${\bf{B}}$ is finite, we have

$\hat{\tau}_{\bf{B}}^{2}-\tau_{\bf{B}}^{2}=\sum\limits_{j\in{\bf{B}}}{\left({\hat{\beta}_{j}^{2}-\beta_{j}^{2}}\right)}\overset{p}{\rightarrow}0,$

and (2.78) follows.

######  Remark 2.9.

We use the the following simple selection algorithm $\gamma$:

Input: A dataset $\left({{{\bf{X}}_{n\times p}},{{\bf{Y}}_{n\times
1}}}\right)$.

  1. 1.

Calculate $\hat{\beta}_{1}^{2},...,\hat{\beta}_{p}^{2}$ where
$\hat{\beta}_{j}^{2}$ is given in (2.4) for $j=1,...,p.$

  2. 2.

Calculate the differences
${\lambda_{j}}=\hat{\beta}_{\left(j\right)}^{2}-\hat{\beta}_{\left({j-1}\right)}^{2}$
for $j=2,\ldots,p$ where
$\hat{\beta}_{\left(1\right)}^{2}<\hat{\beta}_{\left(2\right)}^{2}<...<\hat{\beta}_{\left(p\right)}^{2}$
denotes the order statistics.

  3. 3.

Select the covariates
${{\bf{B}}_{\gamma}}=\left\\{{j:\hat{\beta}_{\left(j\right)}^{2}>\hat{\beta}_{\left({{j^{*}}}\right)}^{2}}\right\\}$,
where ${j^{*}}=\mathop{\arg\max}\limits_{j}{\lambda_{j}}$.

Result: Return $\bf{B_{\gamma}}$.

Algorithm 3 Covariate selection $\gamma$

The algorithm above finds the largest gap between the ordered estimated
squared coefficients and then uses this gap as a threshold to select a set of
coefficients $\bf{B_{\gamma}}\subset\left\\{{1,...,p}\right\\}.$ The algorithm
works well in scenarios where a relatively large gap truly separates between
larger coefficients and the smaller coefficients of the vector $\beta$.

######  Remark 2.10.

The following algorithm is used to construct the Selection estimator that
improves an initial estimator of $\tau^{2}$, as presented in Table 3.

Input: A dataset $\left({{{\bf{X}}_{n\times p}},{Y_{n\times 1}}}\right)$, an
estimation procedure $\tilde{\tau}^{2},$ and a covariate-selection procedure
$\delta$.

  1. 1.

Apply the procedure $\delta$ to the dataset $\left({{{\bf{X}}_{n\times
p}},{Y_{n\times 1}}}\right)$ to obtain ${{\bf{B}}_{\delta}}.$

  2. 2.

Apply the procedure $\tilde{\tau}^{2}$ to the dataset
$\left({{{\bf{X}}_{n\times p}},{Y_{n\times 1}}}\right)$.

  3. 3.

Calculate the zero-estimator
${Z_{h}}=\frac{1}{n}\sum\limits_{i=1}^{n}{h\left({{X_{i}}}\right)}$, where
$h\left({{X_{i}}}\right)=\sum\limits_{j<j^{\prime}\in{{\bf{B}}_{\delta}}}{{X_{ij}}{X_{ij^{\prime}}}}.$

  4. 4.

Bootstrap step:

     * •

Sample $n$ observations at random from $\left({{{\bf{X}}_{n\times
p}},{Y_{n\times 1}}}\right)$, with replacement, to obtain a bootstrap dataset.

     * •

Repeat steps 2 and 3 based on the bootstrap dataset.

The bootstrap step is repeated $M$ times in order to produce
$(\tilde{\tau}^{2})^{*1},...,(\tilde{\tau}^{2})^{*M}$ and
$Z_{h}^{*1},...,Z_{h}^{*M}.$

  5. 5.

Approximate the coefficient
$\tilde{c}_{h}^{*}=\frac{{\widehat{\text{Cov}\left({\tilde{\tau}^{2},{Z_{h}}}\right)}}}{{\text{Var}\left({{Z_{h}}}\right)}}$
where $\widehat{{\text{Cov}}\left(\cdot\right)}$ denotes the empirical
covariance from the bootstrap samples.

Output: Return the empirical estimator
$T_{\tilde{h}}\equiv\tilde{\tau}^{2}-\tilde{c}_{h}^{*}Z_{h}.$

Algorithm 4 Empirical Estimator

The paper “A zero-estimator approach for estimating the signal level in a
high-dimensional model-free setting” is submitted to the Journal of
Statistical Planning and Inference.

##  3 A zero-estimator approach for estimating the signal level in a high-
dimensional model-free setting

###  3.1 Introduction

In many regression settings, an important goal is to estimate the signal and
noise levels, i.e., to quantify the amount of variance in the response
variable that can be explained by a set of covariates, versus how much of the
variation is left unexplained. When the covariates’ dimension is low and a
linear regression model is assumed, the ordinary least squares method can be
used to find a consistent estimator for the signal level. However, in a high-
dimensional setting, the least squares method breaks down and it becomes more
challenging to develop good estimators without further assumptions. In recent
years, several methods have been proposed for estimating the signal level
under the assumption that the regression coefficient vector $\beta$ is sparse
(Fan et al. 2012, Sun and Zhang 2012, Chatterjee and Jafarov 2015, Verzelen et
al. 2018, Tony Cai and Guo 2020). Other widely-used methods assume some
probabilistic structure on $\beta$ (e.g., $\beta$ is Gaussian) and use maximum
likelihood to derive consistent estimators of the signal level (Yang et al.
2010, Bonnet et al. 2015). These methods have been extensively studied in the
literature of _random-effect models_ where $\beta$ is treated as random; see
de Los Campos et al. (2015) and references therein. However, methods that rely
on the assumption that $\beta$ is either sparse or highly structured may not
perform well when these assumptions fail to hold. For example, a known problem
in genetics is the problem of missing heritability (de Los Campos et al. 2015,
Zhu and Zhou 2020). Heritability is defined as the fraction of the observed
outcome (phenotype) that is explained by genetic factors. The term “missing
heritability” is traditionally used to describe the gap between heritability
estimates from genome-wide-association-studies (GWAS) and the corresponding
estimates from family studies. To explain the gap, it has been suggested that
some phenotypes are explained by a numerous number of genetic factors that
their individual effect is too small to detect, but their collective effect is
significant (Yang et al. 2010, Young 2022). In such a setting, methods that
rely on the sparsity assumption may fail to provide accurate estimates.

Rather than assuming sparsity or other structural assumptions on $\beta$, a
different approach for estimating the signal level in a high-dimensional
setting is to assume some or complete knowledge about the covariate
distribution. This can be justified, for example, in the semi-supervised
setting when one has access to a large amount of unlabeled (covariate) data
without the corresponding labels (responses). When the covariates are assumed
independent Gaussian, Dicker (2014) proposed estimators based on the method-
of-moments and Janson et al. (2017) used convex optimization techniques. In
both methods, the Gaussian assumption was used to show consistency and
asymptotic-normality, and it is not clear how robust these methods are when
the assumptions are violated. Dropping the Gaussian independent covariate
assumption, Livne et al. (2021) proposed a consistent estimator under the
assumption that the first two moments of covariates are known. More recently,
Chen (2022) proposed an estimator that is consistent and asymptotically-normal
when the covariates are independent and the entries of $\beta$ are small and
dense.

All of the estimators that we reviewed above were developed under the
assumption that the linear model is true, which can be unrealistic in many
situations. In this work, we focus our attention on the _model-free_ setting,
i.e., no assumptions are made about the relationship between the covariates
and the response. Under this setting, Kong and Valiant (2018) proposed a
consistent estimator under some assumptions on the covariance matrix. In this
paper we follow the two-stage approach presented in Livne et al. (2021), where
an initial estimator is first suggested and then a zero-estimator is used to
reduce its variance. Our initial estimator is the same as in Kong and Valiant
(2018) and Livne et al. (2021), and the zero estimators we use are tailored to
the model-free framework. Furthermore, we provide a general algorithm that, in
principle, can improve any initial estimator and we also demonstrate the
usefulness of the algorithm for several initial estimators.

The rest of this work is organized as follows. In Section 3.2, we discuss the
parameters of interest in a model-free setting under the assumption that the
first two moments of the covariates are known. In Section 3.3, we present our
initial estimators and prove that they are consistent under some minimal
assumptions. In Section 3.4, we use the zero-estimator approach to construct
two improved estimators and then study some theoretical properties of the
improved estimators. Simulation results are given in Section 2.5. Section 3.6
demonstrates how the zero-estimator approach can be generalized to other
estimators. A discussion is given in Section 2.7. The proofs are provided in
the Appendix.

###  3.2 Preliminaries

Let $X\in{{\mathbb{R}}^{p}}$ be a random vector of covariates and let
$Y\in\mathbb{R}$ be the response. The conditional mean $E(Y|X)$ is the _best
predictor_ in the sense that it minimizes the mean squared error
$E\\{[Y-g(X)]^{2}\\}$ over all measurable functions $g(X)$ (see Hansen 2022,
p. $25$). However, the functional form of $E(Y|X)$ is typically unknown and
difficult to estimate, especially in a high-dimensional setting. Consequently,
we can define the best linear approximation to $E(Y|X)$.

######  Definition 3.1.

Assume that both $E(Y^{2})$, $E({{{\left\|X\right\|}^{2}}})$ exist and that
the covariance matrix of $X$, denoted by ${\Sigma_{X}}$, is invertible. Then,
the _best linear predictor_ , $\alpha+\beta^{T}X$, is defined by the unique
$\alpha$ and $\beta$ that minimize the mean squared error

$(\alpha,\beta)=\mathop{\rm
argmin}_{a\in\mathbb{R},b\in\mathbb{R}^{p}}E(Y-a-b^{T}X)^{2},$

and, by Hansen (2022), pp. 34-36, is given by

$\beta={\Sigma_{X}}^{-1}\\{E(XY)-E(X)E(Y)\\}\quad\text{and}\quad\alpha=E(Y)-\beta^{T}E(X).$
(3.1)

The best linear predictor is essentially the population version of the OLS
method. Notice that $\alpha,\beta$ also satisfy $(\alpha,\beta)=\mathop{\rm
argmin}_{a\in\mathbb{R},b\in\mathbb{R}^{p}}E[E(Y|X)-a-b^{T}X]^{2}$. It is a
model-free quantity, i.e., no specific assumptions are made about the
relationship between $X$ and $Y$. In particular, we do not assume that
$E(Y|X)$ is linear in $X$. If $E(Y|X)$ happens to be linear, say
$E(Y|X)=\tilde{\alpha}+\tilde{\beta}^{T}X$, then the best linear predictor
parameters $(\alpha,\beta)$ coincide with the model parameters
$(\tilde{\alpha},\tilde{\beta}).$ However, when $E(Y|X)$ is not linear, the
parameters $\alpha$ and $\beta$ are still meaningful: they describe the
overall direction of the association between $X$ and $Y$ (Buja et al. 2019).
Hence, Definition 3.1 is useful since in most cases we have no reason to
believe that $E(Y|X)$ is indeed linear in $X$.

We now wish to decompose the variance of $Y$ into signal and noise levels. Let
$\sigma_{Y}^{2}$ denote the variance of $Y$ and define the residual
$\varepsilon\equiv Y-\alpha-\beta^{T}X$. Notice that both $E(X\varepsilon)=0$
and $E(\varepsilon)=0$ by construction. Write

$\sigma_{Y}^{2}=\text{Var}(\alpha+X^{T}\beta+\varepsilon)=\beta^{T}\text{Var}(X)\beta+\text{Var}(\varepsilon)=\beta^{T}{\Sigma}_{X}\beta+\sigma^{2},$
(3.2)

where $\sigma^{2}\equiv\text{Var}(\varepsilon)$ and
${\Sigma}_{X}\equiv\text{Var}(X)$. Here, the _signal_ level
$\tau^{2}\equiv\beta^{T}{\Sigma}_{X}\beta$ can be thought of as the total
variance explained by the best linear function of the covariates. The _noise_
level $\sigma^{2}$ is the variance left unexplained. Notice that the
parameters $\tau^{2}$, $\sigma^{2}$ and $\sigma_{Y}^{2}$ depend on $p$ but
this is suppressed in the notation.

A common starting point for many works of regression problems is to use strong
assumptions about $E(Y|X)$, and minimal assumptions, if any, about the
covariate $X$. In this work, we take the opposite approach: we make no
assumptions about $E(Y|X)$ but assume instead that the distribution of $X$ is
known. This can be justified, for example, in a semi-supervised setting where,
in addition to the labeled data, we also have access to a large amount of
unlabeled data; see, for example, the work of Zhang and Bradic (2019) who
study estimation of the mean and variance of $Y$ in a semi-supervised, model-
free setting. A sensitivity analysis of the known-covariate distribution
assumption is carried out in a simulation study in Section 2.5, where the
distribution of $X$ is estimated based on an unlabeled dataset. Note that the
assumption of known-covariate distribution has already been presented and
discussed in the context of high-dimensional regression as in Candes et al.
(2018), Berrett et al. (2020) and Wang and Janson (2020). Our goal here is to
develop good estimators for $\tau^{2}$ and $\sigma^{2}$ in a model-free
setting under the assumption that the distribution of the covariates is known.

Let $(X_{1},Y_{1}),...,(X_{n},Y_{n})$ be i.i.d. observations drawn from an
unknown distribution where $X_{i}\in\mathbb{R}^{p}$ and $Y_{i}\in\mathbb{R}$.
Let $(X,Y)$ denote a generic observation from the sample. We assume that
$E(X)\equiv\mu$ is known and also that the variance matrix ${\Sigma}_{X}$ is
known and invertible. Linear transformations do not affect the signal and
noise levels. Thus, we can apply the transformation
$X\mapsto{\Sigma}_{X}^{-1/2}(X-~{}\mu)$ and assume w.l.o.g. that

$E(X)=0\quad\text{and}\quad{\Sigma}_{X}=\textbf{I}.$ (3.3)

By (3.2), $\sigma_{Y}^{2}=\|\beta\|^{2}+\sigma^{2}$, which means that in order
to estimate $\sigma^{2}$ it is enough to estimate both $\sigma_{Y}^{2}$ and
$\tau^{2}\equiv\|\beta\|^{2}$. The variance term $\sigma_{Y}^{2}$ can be
easily estimated from the sample. Hence, the main challenge is to derive an
estimator for $\tau^{2}.$

###  3.3 Initial Estimators

In this section, we present our initial estimators for the signal and noise
levels, $\tau^{2}$ and $\sigma^{2}$. Interestingly, when the linear model is
true and $X$ is assumed to be constant, no consistent estimator of
$\sigma^{2}$ exists (Azriel 2019). However, when $X$ is random, a consistent
estimator does exist if ${\Sigma}$ is known (Verzelen et al. 2018). The
current work goes one step further as we generalize this result without
assuming linearity. Indeed, Proposition 3.2 below demonstrates that by knowing
the first and second moments of $X$, it is possible, under some mild
assumptions, to construct consistent estimators of $\tau^{2}$ and $\sigma^{2}$
in a high-dimensional setting without assuming that $E(Y|X)$ is linear. The
estimator we use below was suggested by Kong and Valiant (2018) who provided
an upper bound on the variance. Our analysis is more general and we discuss
sufficient conditions for consistency.

Let $W_{ij}\equiv X_{ij}Y_{i}$ for $i=1,...,n$, and $j=1,...,p$. Notice that

$E\left({{W_{ij}}}\right)=E\left({{X_{ij}}{Y_{i}}}\right)=E\left[{{X_{ij}}\left({{\alpha+\beta^{T}}{X_{i}}+{\varepsilon_{i}}}\right)}\right]={\beta_{j}},$
(3.4)

where in the last equality we used (3.3) and the orthogonality between $X$ and
$\varepsilon.$ Now, since

$\\{E(W_{ij})\\}^{2}=E(W_{ij}^{2})-\text{Var}(W_{ij}),$

a natural unbiased estimator for $\beta_{j}^{2}$ is

${\hat{\beta}_{j}^{2}}\equiv\frac{1}{n}\sum\limits_{i=1}^{n}W_{ij}^{2}-\frac{1}{n-1}\sum\limits_{i=1}^{n}(W_{ij}-\overline{W}_{j})^{2}=\binom{n}{2}^{-1}\sum_{i_{1}<i_{2}}^{n}W_{i_{1}j}W_{i_{2}j},$
(3.5)

where $\overline{W}_{j}=\frac{1}{n}\sum_{i=1}^{n}W_{ij}$. Thus, an unbiased
estimator of $\tau^{2}$ is given by

${\hat{\tau}^{2}}=\sum\limits_{j=1}^{p}{\hat{\beta}_{j}^{2}=\binom{n}{2}^{-1}\sum\limits_{{i_{1}}<{i_{2}}}{W_{{i_{1}}}^{T}{W_{{i_{2}}}}}},$
(3.6)

where ${W_{i}}=\left({{W_{i1}},...,{W_{ip}}}\right)^{T}$. We call
$\hat{\tau}^{2}$ the _naive estimator_. Notice that $\hat{\tau}^{2}$ is a
U-statistic with the kernel $h\left({{W_{1}},{W_{2}}}\right)=W_{1}^{T}{W_{2}}$
and thus its variance can be calculated directly by using U-statistic
properties (see van der Vaart (2000), Theorem 12.3).

Let ${\zeta_{1}}\equiv{\beta^{T}}{\bf{A}}\beta-{\left\|\beta\right\|^{4}}$,
${\zeta_{2}}\equiv\left\|{\bf{A}}\right\|_{F}^{2}-{\left\|\beta\right\|^{4}}$,
${\bf{A}}\equiv E\left({{W}W^{T}}\right)$ and $\|\mathbf{A}\|_{F}$ denotes the
Frobenius norm of $\bf{A}.$ The following proposition calculates the variance
of $\hat{\tau}^{2}$.

######  Proposition 3.1.

Assuming that $\zeta_{1}$ and $\zeta_{2}$ are finite, the variance of the
naive estimator $\hat{\tau}^{2}$ is given by

$\text{Var}\left({{{\hat{\tau}}^{2}}}\right)=\frac{{4\left({n-2}\right)}}{{n\left({n-1}\right)}}{\zeta_{1}}+\frac{2}{{n\left({n-1}\right)}}{\zeta_{2}}.$
(3.7)

The next proposition shows that the naive estimator $\hat{\tau}^{2}$ is
consistent under some assumptions.

######  Proposition 3.2.

Assume that
$\frac{\|\textbf{A}\|^{2}_{F}}{n^{2}}\xrightarrow{n\rightarrow\infty}0,$ and
that $\tau^{2}$ is bounded. Then,

$\hat{\tau}^{2}-\tau^{2}\overset{p}{\rightarrow}~{}0.$

Similarly, an estimator for the noise level $\sigma^{2}$ can be obtained by

${\hat{\sigma}^{2}}\equiv\hat{\sigma}_{Y}^{2}-{\hat{\tau}^{2}},$ (3.8)

where
$\hat{\sigma}_{Y}^{2}\equiv\frac{1}{{n-1}}\sum\limits_{i=1}^{n}{{{\left({{Y_{i}}-\bar{Y}}\right)}^{2}}}$
is the standard unbiased estimator of $\sigma_{Y}^{2}.$ Let ${\mu_{4}}\equiv
E\left[{{{\left({Y-\alpha}\right)}^{4}}}\right]$,
$\pi={\left({{\pi_{1}},...,{\pi_{p}}}\right)^{T}}$ where ${\pi_{j}}\equiv
E\left[{{{\left({Y_{1}-\alpha}\right)}^{2}}{W_{1j}}}\right]$. The variance of
$\hat{\sigma}^{2}$ is given by the following proposition.

######  Proposition 3.3.

The variance of $\hat{\sigma}^{2}$ is

$\displaystyle\text{Var}\left({{{\hat{\sigma}}^{2}}}\right)$
$\displaystyle=\left[{\frac{1}{n}{\mu_{4}}-\frac{{\left({n-3}\right)}}{{n\left({n-1}\right)}}\sigma_{Y}^{4}}\right]+\text{Var}(\hat{\tau}^{2})$
$\displaystyle-\frac{4}{n}\left({{\pi^{T}}\beta-{\tau^{2}}\sigma_{Y}^{2}}\right)+\frac{4}{{n\left({n-1}\right)}}\sum\limits_{j=1}^{p}{{{\left\\{{E\left[{{W_{1j}}\left({Y_{1}-\alpha}\right)}\right]}\right\\}}^{2}}},$
(3.9)

where $\text{Var}(\hat{\tau}^{2})$ is given in (2.6).

The following result is a corollary of Propositions 2 and 3.

######  Corollary 3.1.

Assume that $\mu_{4}$ and $\tau^{2}$ are bounded and that $\frac{\|{\bf
A}\|_{F}^{2}}{n^{2}}\rightarrow 0.$ Then,

$\hat{\sigma}^{2}-\sigma^{2}\overset{p}{\rightarrow}~{}0.$

The condition $\frac{\|{\bf A}\|_{F}^{2}}{n^{2}}\rightarrow 0$ holds in
various settings. For example, it can be shown to hold when
$p/n^{2}\rightarrow 0$ and $Y$ is bounded. For more examples and details, see
Remark 3.1 in the Appendix.

###  3.4 Reducing Variance Using a Zero-Estimator

In this section, we study how the naive estimator $\hat{\tau}^{2}$, and
consequently $\hat{\sigma}^{2}$, can be improved by using the assumption that
the distribution of $X$ is known. We use zero-estimators to construct an
improved unbiased estimator of $\tau^{2}$. This is also known as the method of
_control variables_ from the Monte-Carlo literature; see, e.g., Glynn and
Szechtman (2002), and Lavenberg and Welch (1981). Here, a zero-estimator is
defined as a statistic $Z$ such that $E(Z)=0$. For a given zero-estimator $Z$
and a constant $c$, we define a new estimator $\hat{\tau}^{2}(Z,c)$ as

$\hat{\tau}^{2}(Z,c)=\hat{\tau}^{2}-cZ.$ (3.10)

For a fixed $c$, notice that $\hat{\tau}^{2}(Z,c)$ is an unbiased estimator
for $\tau^{2}.$ Also notice that for every function $f$ of the covariates
$X_{1},...,X_{n}$, one can always define a zero-estimator

${Z}={f\left({{X_{1},\dots,X_{n}}}\right)-E\left[{f\left({{X_{1},\dots,X_{n}}}\right)}\right]}.$

This is possible since we assume that the distribution of the covariates is
known and hence $E[f\left({{X_{1},\dots,X_{n}}}\right)]$ is known. The
variance of $\hat{\tau}^{2}(Z,c)$ is

$\text{Var}[\hat{\tau}^{2}(Z,c)]=\text{Var}\left({{{\hat{\tau}}^{2}}}\right)+{c^{2}}\text{Var}\left({{Z}}\right)-2c\,\text{Cov}\left({{{\hat{\tau}}^{2}},{Z}}\right).$

Minimizing the variance with respect to $c$ yields the minimizer

${c^{*}}=\frac{{\text{Cov}\left({{{\hat{\tau}}^{2}},{Z}}\right)}}{{\text{Var}\left({{Z}}\right)}}.$
(3.11)

Hence, the corresponding oracle-estimator is

$\hat{\tau}^{2}(Z,c^{*})=\hat{\tau}^{2}-c^{*}Z.$

We use the term _oracle_ since the optimal coefficient $c^{*}$ is an unknown
quantity. The variance of the above oracle-estimator is

$\displaystyle\text{Var}[\hat{\tau}^{2}(Z,c^{*})]$
$\displaystyle=\text{Var}\left({{{\hat{\tau}}^{2}}}\right)-\frac{{\left[{\text{Cov}\left({{{\hat{\tau}}^{2}},{Z}}\right)}\right]^{2}}}{{\text{Var}\left({{Z}}\right)}}$
$\displaystyle=\left({1-\rho_{{{\hat{\tau}}^{2}},{Z}}^{2}}\right)\text{Var}\left({{{\hat{\tau}}^{2}}}\right),$
(3.12)

where $\rho_{{{\hat{\tau}}^{2}},{Z}}$ is the correlation coefficient between
$\hat{\tau}^{2}$ and $Z$. The term ${1-\rho_{{{\hat{\tau}}^{2}},{Z}}^{2}}$ is
the factor by which $\text{Var}(\hat{\tau}^{2})$ could be reduced if the
optimal coefficient $c^{*}$ was known. Thus, the more correlation there is
between the zero-estimator $Z$ and the naive estimator $\hat{\tau}^{2}$, the
greater the reduction in variance.

There are two challenges to be addressed with the above approach. First, one
should find a simple zero-estimator $Z$ which is correlated with the naive
estimator $\hat{\tau}^{2}$. Second, the optimal coefficient $c^{*}$ is an
unknown quantity and therefore needs to be estimated.

To address the first challenge, we propose the following zero-estimator

$Z_{g}\equiv{\bar{g}_{n}}=\frac{1}{n}\sum\limits_{i=1}^{n}{{g}(X_{i})},$

where ${g}(X_{i})=\sum\limits_{j<j^{\prime}}{{X_{ij}}{X_{ij^{\prime}}}}.$ In
Remark 3.2 in the Appendix, we show that the optimal coefficient, with respect
to $Z_{g}$, is

${c_{g}^{*}}=\frac{{{2\beta^{T}}\theta_{g}}}{{\text{Var}[{{g(X)}}]}},$ (3.13)

where ${\theta_{g}}=E\left[{Wg\left(X\right)}\right].$ Notice that
$\text{Var}[{g(X)}]$ is a known quantity since the distribution of $g(X)$ is
assumed to be known. Hence, the corresponding oracle-estimator is

$\displaystyle
T_{g}\equiv\hat{\tau}^{2}(Z_{g},c_{g}^{*})=\hat{\tau}^{2}-c_{g}^{*}Z_{g}.$
(3.14)

To address the second challenge, i.e., to estimate the optimal coefficient
$c_{g}^{*}$, we suggest the following unbiased U-statistic estimator

$\hat{c}_{g}^{*}=\frac{\binom{n}{2}^{-1}{\sum\limits_{{i_{1}}\neq{i_{2}}}{W_{{i_{1}}}^{T}{W_{{i_{2}}}}g\left({{X_{{i_{2}}}}}\right)}}}{{\text{Var}\left[{g\left(X\right)}\right]}}.$
(3.15)

Thus, the corresponding improved estimator is

$T_{\hat{g}}\equiv\hat{\tau}^{2}(Z_{g},\hat{c}_{g}^{*})={\hat{\tau}^{2}}-{\hat{c}_{g}^{*}}Z_{g}.$
(3.16)

Using the zero-estimator $Z_{g}$ has a potential drawback. It uses _all_ the
$p$ covariates of the vector $X$ regardless of the sparsity level in the data,
which can result in some additional variability due to unnecessary estimation.
Intuitively, when the sparsity level is high, i.e., only a small number of
covariates plays an important role in explaining the response $Y$, it is
inefficient to use a zero-estimator that incorporates all the $p$ covariates.
In such a setting, it is reasonable to modify the zero-estimator $Z_{g}$ such
that only a small set of covariates will be included, preferably the
covariates that capture a significant part of the signal level $\tau^{2}.$
Selecting such a set of covariates can be difficult and one may use a
covariate-selection procedure for this purpose.

We call $\delta$ a covariate _selection procedure_ if for every dataset it
chooses a subset of indices ${{\bf{S}}_{\delta}}$ from
$\left\\{{1,...,p}\right\\}$. Different covariate-selection methods exist in
the literature (see Oda et al. 2020 and references therein) but these are not
a primary focus of this work. For a given selection procedure $\delta$ we
modify the estimator $T_{\hat{g}}$ such that only the indices in
${\bf{S}}_{\delta}$ will be included in its zero-estimator term. This modified
estimator, which is based on a given selection procedure $\delta$, is
presented in the algorithm below.

Input: A dataset $\left({{{\bf{X}}_{n\times p}},{Y_{n\times 1}}}\right)$ and a
selection procedure $\delta$.

  1. 1.

Calculate the naive estimator
${\hat{\tau}^{2}}=\binom{n}{2}^{-1}\sum\limits_{{i_{1}}<{i_{2}}}{W_{{i_{1}}}^{T}{W_{{i_{2}}}}},$

  2. 2.

Apply procedure $\delta$ to $\left({{{\bf{X}}_{n\times p}},{Y_{n\times
1}}}\right)$ to construct the set ${{\bf{S}}_{\delta}}.$

Output: Return the estimator
$T_{\hat{h}}\equiv{\hat{\tau}^{2}}-\hat{c}_{h}^{*}{Z_{h}},$ (3.17) where
$\hat{c}_{h}^{*}=\frac{{\binom{n}{2}^{-1}\sum\limits_{{i_{1}}\neq{i_{2}}}{W_{{i_{1}}}^{T}{W_{{i_{2}}}}h\left({{X_{{i_{2}}}}}\right)}}}{{\text{Var}\left[{h\left({{X}}\right)}\right]}}$,
${Z_{h}}=\frac{1}{n}\sum\limits_{i=1}^{n}{{h}\left({{X_{i}}}\right)}$ and
$h\left({{X_{i}}}\right)=\sum\limits_{j<j^{\prime}\in{\bf{S}}_{\delta}}{{X_{ij}}{X_{ij^{\prime}}}}$.

Algorithm 5 Proposed estimator for $\tau^{2}.$

Notice that the estimator $T_{\hat{g}}$ defined in (3.16) is a special case of
the estimator $T_{\hat{h}}$ defined in Algorithm 5, when
${\bf{S}}_{\delta}=\\{1,\dots p\\}$, i.e., when $\delta$ selects all the $p$
covariates.

Recall that in this work we treat $p$ as a function of $n$, i.e., $p\equiv
p_{n}$ but this is suppressed in the notation. Let
$\Theta\equiv\Theta_{n}\subseteq\\{1,\ldots,p\\}$ be a deterministic sequence
of subsets. In order to analyze the estimator $T_{\hat{h}}$ we define a
stability property, which is given next.

######  Definition 3.2.

A selection procedure $\delta$ is _stable_ if there exists a deterministic
sequence of subsets $\Theta$ such that

$\lim_{n\to\infty}n[P({\bf S}_{\delta}\neq\Theta)]^{1/2}=0.$ (3.18)

Definition 3.2 states that a selection procedure $\delta$ is stable if it is
asymptotically close to a deterministic procedure at a suitable rate. The
convergence rate of many practical selection procedures is exponential, which
is much faster than is required for the condition to hold. For example, the
lasso algorithm asymptotically selects the support of $\beta$ at an
exponential rate under some assumptions (see Hastie et al. 2015, Theorem
11.3). Notice also that the stability condition holds trivially when
${\bf{S}}_{\delta}=\\{1,\dots p\\}$, i.e., when $\delta$ selects all the $p$
covariates for all $n$.

Define the oracle-estimator
${T_{h}}\equiv{\hat{\tau}^{2}}\left({{Z_{h}},c_{h}^{*}}\right)$, where
$c_{h}^{*}=\frac{{{\beta^{T}}{\theta_{h}}}}{{\text{Var}\left[{h\left(X\right)}\right]}}$
and ${\theta_{h}}=E\left[{Wh\left(X\right)}\right]$. Let

$\displaystyle f\left({{X_{i}}}\right)$
$\displaystyle\equiv\sum\limits_{j<j^{\prime}\in\Theta}{{X_{ij}}{X_{ij^{\prime}}}},\quad{T_{f}}\equiv{\hat{\tau}^{2}}\left({{Z_{f}},c_{f}^{*}}\right),\quad{T_{\hat{f}}}\equiv{\hat{\tau}^{2}}\left({{Z_{f}},\hat{c}_{f}^{*}}\right),\quad
c_{f}^{*}\equiv\frac{{2{\beta^{T}}{\theta_{f}}}}{{\text{Var}\left[{f\left(X\right)}\right]}},$
$\displaystyle\hat{c}_{f}^{*}$
$\displaystyle\equiv\frac{{\frac{2}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{W_{{i_{1}}}^{T}{W_{{i_{2}}}}f\left({{X_{{i_{2}}}}}\right)}}}{{\text{Var}\left[{f\left(X\right)}\right]}},\quad
b\equiv Wf(X),\quad\theta_{b}\equiv
E(b),\quad{Z_{f}}\equiv\frac{1}{n}\sum\limits_{i=1}^{n}{f\left({{X_{i}}}\right)},$
$\displaystyle{\bf{B}}$
$\displaystyle=E\left[{W{W^{T}}f\left(X\right)}\right],\quad{\bf{C}}=E\left[{W{W^{T}}f^{2}{{\left(X\right)}}}\right].$
(3.19)

We now prove that the proposed estimator in Algorithm 5 is asymptotically
equivalent to its oracle version $T_{h}$ under some conditions.

######  Proposition 3.4.

Assume that the selection procedure $\delta$ is stable with respect to
$\Theta$. Assume also that $n/p$, $\|\beta\|^{2}$,
$\frac{{{{\left\|\theta_{b}\right\|}^{2}}}}{{\text{Var}[f(X)]}}$ and
$\frac{{E\left({{{\left\|b\right\|}^{2}}}\right)}}{{n\text{Var}[f(X)]}}$ are
bounded, and $\frac{{\left\|{\bf{A}}\right\|_{F}^{2}}}{{{n^{2}}}}\to 0$,
$\frac{{\left\|{\bf{B}}\right\|_{F}^{2}}}{{{n^{2}\text{Var}[f(X)]}}}\to~{}0$,
and
$\frac{{\left\|{\bf{C}}\right\|_{F}^{2}}}{{\\{n\text{Var}[f(X)]\\}^{2}}}\to
0.$ In addition, assume that the first four moments of $T_{h}$, $T_{\hat{h}}$,
$T_{f}$, and $T_{\hat{f}}$ are bounded. Then,

$\sqrt{n}\left[T_{\hat{h}}-T_{h}\right]\overset{p}{\rightarrow}~{}0.$ (3.20)

Our proof of Proposition 3.4 shows a slightly stronger result: the proposed
estimator $T_{\hat{h}}$ is also asymptotically equivalent to $T_{f}$, the
oracle-estimator that originally knows the set of indices $\Theta$.

We now discuss the assumptions of Proposition 3.4. In Remark 3.1, several
sufficient conditions implying that
$\frac{{\left\|{\bf{A}}\right\|_{F}^{2}}}{{{n^{2}}}}\to 0$ were presented .
Similarly, in Remark 3.3 we show that if the covariates $X_{ij}$, for
$j=1,...,p$, $i=1,...,n,$ and the response $Y$ are bounded, then so is
$\frac{{E\left({{{\left\|b\right\|}^{2}}}\right)}}{{n\text{Var}[f(X)]}}$. It
is also shown, that if in addition $|\Theta|$ is bounded and
$\text{Var}[f(X)]$ is bounded away from zero, then
$\frac{{\left\|{\bf{B}}\right\|_{F}^{2}}}{{{n^{2}\text{Var}[f(X)]}}}\to 0$ and
$\frac{{\left\|{\bf{C}}\right\|_{F}^{2}}}{{{\\{n\text{Var}[f(X)]\\}^{2}}}}\to
0.$ Proposition 4 in Livne et al. (2021) shows that
$\frac{{\left\|{\bf{B}}\right\|_{F}^{2}}}{{{n^{2}\text{Var}[f(X)]}}}\to 0$ and
$\frac{{\left\|{\bf{C}}\right\|_{F}^{2}}}{{{\\{n\text{Var}[f(X)]\\}^{2}}}}\to
0$ hold also when $|\Theta|$ is unbounded, but with additional conditions on
linearity and independence of the covariates. It is also shown there that
under those assumptions,
$\frac{{{{\left\|\theta_{b}\right\|}^{2}}}}{{\text{Var}[f(X)]}}$ is bounded.
In simulations, which are not presented here, we observed that these
conditions also hold for various non-linear models.

Notice that the zero-estimator $Z_{g}$, which is not based on any covariate-
selection procedure, is just a special case of $Z_{h}$ when $\delta$ selects
all the $p$ covariates, i.e., ${\bf{S}}_{\delta}=\\{1,\dots p\\}$. Hence, if
the conditions of Proposition 3.4 hold, then
$\sqrt{n}\left[T_{\hat{g}}-T_{g}\right]\overset{p}{\rightarrow}~{}0,$ where
$T_{\hat{g}}$ and $T_{g}$ are given in (3.16) and (3.14), respectively.

###  3.5 Simulations Results

In this section, we illustrate the performance of the proposed estimators
using simulations. Specifically, we compare the naive estimator
$\hat{\tau}^{2}$ and the improved estimators $T_{\hat{g}}$ and $T_{\hat{h}}$
which are defined in (2.5), (3.16), and Algorithm 5, respectively. The code
for reproducing the results of this section and the next section (2.6) is
available at https://t.ly/dwJg.

For demonstration purposes, we consider a setting in which $K$ entries of the
vector $\beta$ are relatively large (in absolute value), and all other entries
are small. The proportion of the signal in those $K$ entries is defined as the
_sparsity level_ of the vector $\beta$. Next, we study different sparsity
levels by defining the following non-linear model,

${Y_{i}}={\gamma_{L}}\sum\limits_{j\in\Theta}{\left[{{X_{ij}}+\sin\left({{X_{ij}}}\right)}\right]}+{\gamma_{S}}\sum\limits_{j\notin\Theta}{\left[{{X_{ij}}+\sin\left({{X_{ij}}}\right)}\right]}+{\xi_{i}},\quad\quad
i=1,\dots,n,$ (3.21)

where
${\gamma_{L}}\equiv{\left\\{{\frac{{\eta{\tau^{2}}}}{{k{{\left({1+E\left[{X\sin\left(X\right)}\right]}\right)}^{2}}}}}\right\\}^{1/2}}$
,
${\gamma_{S}}\equiv{\left\\{{\frac{{{\tau^{2}}\left({1-\eta}\right)}}{{\left({p-k}\right){{\left({1+E\left[{X\sin\left(X\right)}\right]}\right)}^{2}}}}}\right\\}^{1/2}}$,
and $\Theta$ is the set of the largest $K$ entries of the vector $\beta.$ The
model has two parameters, $\tau^{2}$ and $\eta$, that vary across the
different simulation scenarios. The covariates were generated from the
centered exponential distribution, i.e.,
$X_{ij}\overset{iid}{\sim}\text{Exp}(1)-1,$ $i=1,\dots,n$, $j=1,\dots,p$. The
noise level $\xi_{i}$ was generated from the standard normal distribution. One
can verify that under the above model
$\beta_{j}^{2}=\frac{{\eta{\tau^{2}}}}{K}$ for $j\in\Theta$, and that
$\beta_{j}^{2}=\frac{{{\tau^{2}}\left({1-\eta}\right)}}{{\left({p-K}\right)}}\,$
for $j\notin\Theta.$ Define
$\tau_{\Theta}^{2}\equiv\sum\limits_{j\in\Theta}{\beta_{j}^{2}}$. From the
above definitions, it follows that $\eta=~{}\tau_{\Theta}^{2}/{\tau^{2}}$. The
parameter $\eta$ is the proportion of signal that is captured by the set
$\Theta,$ which is the sparsity level as defined above. The case of full
sparsity, where the entire signal level $\tau^{2}$ comes only from the set
$\Theta$, corresponds to $\eta=1$, and is not assumed here.

We fix $n=p=300$ and $K=6$. For each combination of the parameters
$\tau^{2}\in\\{1,2\\}$ and $\eta\in\\{{0.1,0.3,0.5,0.7,0.9}\\},$ we generated
100 independent datasets from model (3.21) and estimated $\tau^{2}$ using the
different estimators. The covariate-selection procedure $\delta$ that was used
in the estimator $T_{\hat{h}}$ is defined in Remark 2.9 in the Appendix.

Figure 4 plots the RMSE of each estimator as a function of the sparsity level
$\eta$ and the signal level $\tau^{2}.$ It is demonstrated that the estimators
$T_{\hat{g}}$ and $T_{\hat{h}}$ improve (i.e., lower or equal RMSE) the naive
estimator in all settings. The improved estimators are complementary to each
other, i.e., for small values of $\eta$ the estimator $T_{\hat{g}}$ performs
better than $T_{\hat{h}},$ and the opposite occurs for large values of $\eta.$
This is expected since when the sparsity level $\eta$ is small, the
improvement of $T_{\hat{h}}$ is smaller as it ignores much of the signal that
lies outside of the set $\Theta.$ On the other hand, when a large portion of
the signal $\tau^{2}$ is captured by only the few covariates in $\Theta$, it
is sufficient to make use of only these covariates in the zero-estimator term,
and the improvement of $T_{\hat{h}}$ is greater.

Table 1 shows the RMSE, bias, standard error, and the relative improvement,
for the different estimators. It can be observed that the degree of
improvements depends on the sparsity level of the data $\eta$. For example,
when $\tau^{2}=1$ and sparsity level is low ($\eta=0.1$), the estimator
$T_{\hat{g}}$ improves the naive estimator by $11\%$, while the estimator
$T_{\hat{h}}$ presents a similar performance to the naive estimator. On the
other hand, when the sparsity level is high ($\eta=0.9$), the estimator
$T_{\hat{h}}$ improves the naive by $11\%$, while $T_{\hat{g}}$ presents a
similar performance to the naive estimator, as expected. Notice that when
$\tau^{2}=2$ these improvements are even more substantial.

Figure 4:  Root mean square error (RMSE) for the proposed estimators. The
x-axis stands for the sparsity level $\eta$.

Table 5: Summary statistics for the proposed estimators; $n=p=300.$ Bias,
standard error (SE), root mean square error (RMSE) and percentage change from
the naive estimator (in terms of RMSE) are shown. The table results were
computed over $100$ simulated datasets for each setting. An estimate for the
standard deviation of RMSE ($\hat{\sigma}_{RMSE}$) was calculated using the
delta method.

$\eta$ | $\tau^{2}$ | Estimator | Bias | SE | RMSE | % Change | $\hat{\sigma}_{RMSE}$
---|---|---|---|---|---|---|---
10% | 1 | $\hat{\tau}^{2}$ | -0.03 | 0.274 | 0.274 | 0.00 | 0.018
10% | 1 | $T_{\hat{g}}$ | 0.04 | 0.242 | 0.244 | -10.95 | 0.016
10% | 1 | $T_{\hat{h}}$ | -0.03 | 0.274 | 0.274 | 0.00 | 0.017
30% | 1 | $\hat{\tau}^{2}$ | -0.02 | 0.284 | 0.283 | 0.00 | 0.019
30% | 1 | $T_{\hat{g}}$ | 0.05 | 0.255 | 0.258 | -8.83 | 0.018
30% | 1 | $T_{\hat{h}}$ | 0.00 | 0.282 | 0.28 | -1.06 | 0.017
50% | 1 | $\hat{\tau}^{2}$ | 0.00 | 0.286 | 0.284 | 0.00 | 0.021
50% | 1 | $T_{\hat{g}}$ | 0.05 | 0.264 | 0.268 | -5.63 | 0.020
50% | 1 | $T_{\hat{h}}$ | 0.02 | 0.277 | 0.276 | -2.82 | 0.017
70% | 1 | $\hat{\tau}^{2}$ | 0.01 | 0.285 | 0.284 | 0.00 | 0.022
70% | 1 | $T_{\hat{g}}$ | 0.05 | 0.272 | 0.275 | -3.17 | 0.021
70% | 1 | $T_{\hat{h}}$ | 0.04 | 0.265 | 0.267 | -5.99 | 0.017
90% | 1 | $\hat{\tau}^{2}$ | 0.03 | 0.281 | 0.282 | 0.00 | 0.021
90% | 1 | $T_{\hat{g}}$ | 0.05 | 0.276 | 0.279 | -1.06 | 0.020
90% | 1 | $T_{\hat{h}}$ | 0.07 | 0.242 | 0.25 | -11.35 | 0.015
10% | 2 | $\hat{\tau}^{2}$ | -0.06 | 0.458 | 0.46 | 0.00 | 0.030
10% | 2 | $T_{\hat{g}}$ | 0.08 | 0.379 | 0.386 | -16.09 | 0.024
10% | 2 | $T_{\hat{h}}$ | -0.05 | 0.457 | 0.458 | -0.43 | 0.029
30% | 2 | $\hat{\tau}^{2}$ | -0.03 | 0.476 | 0.475 | 0.00 | 0.033
30% | 2 | $T_{\hat{g}}$ | 0.09 | 0.408 | 0.416 | -12.42 | 0.029
30% | 2 | $T_{\hat{h}}$ | -0.01 | 0.469 | 0.466 | -1.89 | 0.029
50% | 2 | $\hat{\tau}^{2}$ | -0.01 | 0.481 | 0.478 | 0.00 | 0.036
50% | 2 | $T_{\hat{g}}$ | 0.09 | 0.431 | 0.439 | -8.16 | 0.033
50% | 2 | $T_{\hat{h}}$ | 0.03 | 0.458 | 0.457 | -4.39 | 0.028
70% | 2 | $\hat{\tau}^{2}$ | 0.02 | 0.477 | 0.475 | 0.00 | 0.038
70% | 2 | $T_{\hat{g}}$ | 0.09 | 0.448 | 0.454 | -4.42 | 0.035
70% | 2 | $T_{\hat{h}}$ | 0.08 | 0.429 | 0.433 | -8.84 | 0.026
90% | 2 | $\hat{\tau}^{2}$ | 0.05 | 0.468 | 0.468 | 0.00 | 0.035
90% | 2 | $T_{\hat{g}}$ | 0.08 | 0.456 | 0.462 | -1.28 | 0.034
90% | 2 | $T_{\hat{h}}$ | 0.12 | 0.376 | 0.392 | -16.24 | 0.023

In the following we consider the case that the distribution of the covariates
is only partially known. It is assumed that a large amount of unlabeled data
is available and therefore the distribution of the covariates can be estimated
by using this data. Specifically, we assume that in addition to the label data
$(X_{1},Y_{1}),...,(X_{n},Y_{n}),$ we have access to additional $N$ i.i.d
observations $X_{n+1},...,X_{n+N}.$ However, the responses
$Y_{n+1},...,Y_{n+N}$ are missing. Additionally, a specific correlation
structure is assumed. Namely, we consider the case that $\Sigma_{X}$ is a
$b$-banded matrix, i.e., $[\Sigma_{X}]_{jj^{\prime}}=0$ for all
$|j-j^{\prime}|>b,$ where $b$ is referred as the bandwidth of the matrix
diagonal. This assumption is common in GWAS; see e.g., Bickel and Levina
(2008) and Cai et al. (2010). Unlike in the previous simulation setting, here
$\mu$ and ${\Sigma_{X}}$ are estimated using the unlabeled data, namely,
$\hat{\mu}\equiv\frac{1}{N}\sum_{i={n+1}}^{n+N}X_{i}$ and

$[{\hat{\Sigma}_{X}}]_{jj^{\prime}}\equiv\left\\{\begin{array}[]{cc}\frac{1}{N}\sum_{i=n+1}^{n+N}(X_{i,j}-\hat{\mu}_{j})(X_{i,j^{\prime}}-\hat{\mu}_{j^{\prime}})&|j-j^{\prime}|\leq
b\\\ 0&|j-j^{\prime}|>b\end{array}\right..$

We then apply the linear transformation,
$X\mapsto{\hat{\Sigma}_{X}}^{-1/2}(X-\hat{\mu}),$ which corresponds the
transformation shown in Section 3.2. We fixed $b=5$ and repeated the
simulation study above for different values of $N.$ Notice that under the
above setting, the true $\Sigma_{X}$ is the identity, and in particular it is
indeed a $5$-banded matrix.

Table 6 is similar to Table 1 but includes different values of $N$ rather than
different values of $\eta.$ For brevity, we present only the scenario of
$\tau^{2}=2$ and $\eta=0.9$ but the results for other scenarios are similar.
It can be observed from the table that for large values of $N$ the results are
fairly similar to those in Table 1.

Table 6: Summary statistics similar to Table 1; $n=p=300;$ $\eta=0.9;$
$\tau^{2}=2.$

$N$ | Estimator | Bias | SE | RMSE | % Change | $\hat{\sigma}_{RMSE}$
---|---|---|---|---|---|---
5000 | $\hat{\tau}^{2}$ | 0.01 | 0.478 | 0.476 | 0 | 0.038
5000 | $T_{\hat{g}}$ | 0.05 | 0.466 | 0.466 | -2.1 | 0.036
5000 | $T_{\hat{h}}$ | 0.1 | 0.38 | 0.39 | -18.07 | 0.023
10000 | $\hat{\tau}^{2}$ | 0.1 | 0.457 | 0.466 | 0 | 0.033
10000 | $T_{\hat{g}}$ | 0.14 | 0.447 | 0.465 | -0.21 | 0.033
10000 | $T_{\hat{h}}$ | 0.1 | 0.378 | 0.388 | -16.74 | 0.023
20000 | $\hat{\tau}^{2}$ | 0.06 | 0.467 | 0.469 | 0 | 0.035
20000 | $T_{\hat{g}}$ | 0.1 | 0.456 | 0.463 | -1.28 | 0.034
20000 | $T_{\hat{h}}$ | 0.12 | 0.377 | 0.392 | -16.42 | 0.023
$\infty$ | $\hat{\tau}^{2}$ | 0.05 | 0.468 | 0.468 | 0.00 | 0.035
$\infty$ | $T_{\hat{g}}$ | 0.08 | 0.456 | 0.462 | -1.28 | 0.034
$\infty$ | $T_{\hat{h}}$ | 0.12 | 0.376 | 0.392 | -16.24 | 0.023

###  3.6 Generalization to Other Estimators

The suggested methodology in this paper is not limited to improving only the
naive estimator, but can also be generalized to other estimators. As before,
the key idea is to use a zero-estimator that is correlated with an initial
estimator of $\tau^{2}$ in order to reduce its variance. Unlike the naive
estimator $\hat{\tau}^{2}$, which has a closed-form expression, other
estimators, such as the EigenPrism estimator (Janson et al. 2017), are
computed numerically by solving a convex optimization problem. For a given
zero-estimator, this makes the task of estimating the optimal-coefficient
$c^{*}$ more challenging than before. To overcome this challenge, we
approximate the optimal coefficient $c^{*}$ using bootstrap samples. This is
described in the following algorithm.

Input: A dataset $\left({{{\bf{X}}_{n\times p}},{Y_{n\times 1}}}\right)$, an
estimation procedure $\tilde{\tau}^{2},$ and a covariate-selection procedure
$\delta$.

  1. 1.

Apply the procedure $\delta$ to the dataset $\left({{{\bf{X}}_{n\times
p}},{Y_{n\times 1}}}\right)$ to obtain ${{\bf{S}}_{\delta}}.$

  2. 2.

Apply the procedure $\tilde{\tau}^{2}$ to the dataset
$\left({{{\bf{X}}_{n\times p}},{Y_{n\times 1}}}\right)$.

  3. 3.

Calculate the zero-estimator
${Z_{h}}=\frac{1}{n}\sum\limits_{i=1}^{n}{h\left({{X_{i}}}\right)}$, where
$h\left({{X_{i}}}\right)=\sum\limits_{j<j^{\prime}\in{{\bf{S}}_{\delta}}}{{X_{ij}}{X_{ij^{\prime}}}}.$

  4. 4.

Bootstrap step:

     * •

Sample $n$ observations at random from $\left({{{\bf{X}}_{n\times
p}},{Y_{n\times 1}}}\right)$, with replacement, to obtain a bootstrap dataset.

     * •

Repeat steps 2 and 3 based on the bootstrap dataset.

The bootstrap step is repeated $M$ times in order to produce
$(\tilde{\tau}^{2})^{*1},...,(\tilde{\tau}^{2})^{*M}$ and
$Z_{h}^{*1},...,Z_{h}^{*M}.$

  5. 5.

Approximate the coefficient
$\tilde{c}_{h}^{*}=\frac{{\widehat{\text{Cov}\left({\tilde{\tau}^{2},{Z_{h}}}\right)}}}{{\text{Var}\left({{Z_{h}}}\right)}}$
where $\widehat{{\text{Cov}}\left(\cdot\right)}$ denotes the empirical
covariance from the bootstrap samples.

Output: Return the empirical estimator
$T_{\tilde{h}}\equiv\tilde{\tau}^{2}-\tilde{c}_{h}^{*}Z_{h}.$

Algorithm 6 Empirical Estimators

In the special case when $\delta$ selects all the $p$ covariates, i.e.,
${\bf{S}}_{\delta}=\\{1,\dots p\\}$, we use the notations $Z_{g}$ and
$\tilde{c}_{g}^{*}$ rather than $Z_{h}$ and $\tilde{c}_{h}^{*}$, respectively,
i.e., $T_{\tilde{g}}\equiv\tilde{\tau}^{2}-\tilde{c}_{g}^{*}Z_{g}.$

We illustrate the improvement obtained by Algorithm 6 by choosing
$\tilde{\tau}^{2}$ to be the EigenPrism procedure (Janson et al. 2017), but
other estimators can be used as well. We consider the same setting as in
Section 2.5. The number of bootstrap samples is $M=100.$

The simulation results appear in Table 7 and Figure 5. Both estimators
$T_{\tilde{h}}$ and $T_{\tilde{g}}$ show an improvement over the EigenPrism
estimator $\tilde{\tau}^{2}.$ The results here are fairly similar to the
results shown for the naive estimator in Section 2.5, with just a smaller
degree of improvement. As before, the improved estimators $T_{\tilde{h}}$ and
$T_{\tilde{g}}$ are complementary to each other, i.e., for small values of
$\eta$ the estimator $T_{\tilde{g}}$ performs better than $T_{\tilde{h}},$ and
the opposite occurs for large values of $\eta.$

Figure 5:  Root mean square error (RMSE) for the proposed estimators. The
x-axis stands for the sparsity level $\eta$.  Table 7: Summary statistics for
the EigenPrism-related estimators; $n=p=300.$ Biases, standard errors (SE) and
the root mean square errors (RMSE) of the different estimators, computed over
$100$ independent datasets for each setting. The relative improvement over the
EigenPrism estimator $\tilde{\tau}^{2}$ (in terms of RMSE) is also shown. An
estimate for the standard deviation of RMSE ($\hat{\sigma}_{RMSE}$) was
calculated using the delta method.

$\eta$ | $\tau^{2}$ | Estimator | Bias | SE | RMSE | % Change | $\hat{\sigma}_{RMSE}$
---|---|---|---|---|---|---|---
10% | 1 | $\tilde{\tau}^{2}$ | 0.01 | 0.204 | 0.204 | 0.00 | 0.014
10% | 1 | $T_{\tilde{g}}$ | 0.01 | 0.19 | 0.19 | -6.86 | 0.012
10% | 1 | $T_{\tilde{h}}$ | 0.01 | 0.204 | 0.203 | -0.49 | 0.014
30% | 1 | $\tilde{\tau}^{2}$ | 0.01 | 0.202 | 0.201 | 0.00 | 0.014
30% | 1 | $T_{\tilde{g}}$ | 0.01 | 0.195 | 0.194 | -3.48 | 0.014
30% | 1 | $T_{\tilde{h}}$ | 0.01 | 0.201 | 0.201 | 0.00 | 0.014
50% | 1 | $\tilde{\tau}^{2}$ | 0.00 | 0.206 | 0.205 | 0.00 | 0.014
50% | 1 | $T_{\tilde{g}}$ | 0.00 | 0.202 | 0.201 | -1.95 | 0.014
50% | 1 | $T_{\tilde{h}}$ | 0.01 | 0.202 | 0.201 | -1.95 | 0.015
70% | 1 | $\tilde{\tau}^{2}$ | 0.00 | 0.212 | 0.211 | 0.00 | 0.015
70% | 1 | $T_{\tilde{g}}$ | 0.00 | 0.209 | 0.208 | -1.42 | 0.015
70% | 1 | $T_{\tilde{h}}$ | 0.00 | 0.203 | 0.202 | -4.27 | 0.015
90% | 1 | $\tilde{\tau}^{2}$ | -0.01 | 0.219 | 0.218 | 0.00 | 0.016
90% | 1 | $T_{\tilde{g}}$ | -0.01 | 0.218 | 0.218 | 0.00 | 0.016
90% | 1 | $T_{\tilde{h}}$ | 0.00 | 0.204 | 0.203 | -6.88 | 0.015
10% | 2 | $\tilde{\tau}^{2}$ | 0.00 | 0.291 | 0.29 | 0.00 | 0.020
10% | 2 | $T_{\tilde{g}}$ | 0.02 | 0.251 | 0.25 | -13.79 | 0.016
10% | 2 | $T_{\tilde{h}}$ | 0.01 | 0.291 | 0.29 | 0.00 | 0.020
30% | 2 | $\tilde{\tau}^{2}$ | 0.03 | 0.279 | 0.279 | 0.00 | 0.019
30% | 2 | $T_{\tilde{g}}$ | 0.03 | 0.258 | 0.258 | -7.53 | 0.018
30% | 2 | $T_{\tilde{h}}$ | 0.03 | 0.276 | 0.276 | -1.08 | 0.019
50% | 2 | $\tilde{\tau}^{2}$ | 0.02 | 0.285 | 0.284 | 0.00 | 0.019
50% | 2 | $T_{\tilde{g}}$ | 0.01 | 0.273 | 0.272 | -4.23 | 0.019
50% | 2 | $T_{\tilde{h}}$ | 0.02 | 0.274 | 0.273 | -3.87 | 0.020
70% | 2 | $\tilde{\tau}^{2}$ | 0.01 | 0.294 | 0.293 | 0.00 | 0.020
70% | 2 | $T_{\tilde{g}}$ | 0.00 | 0.289 | 0.287 | -2.05 | 0.019
70% | 2 | $T_{\tilde{h}}$ | 0.01 | 0.271 | 0.27 | -7.85 | 0.019
90% | 2 | $\tilde{\tau}^{2}$ | -0.01 | 0.308 | 0.307 | 0.00 | 0.022
90% | 2 | $T_{\tilde{g}}$ | -0.01 | 0.307 | 0.305 | -0.65 | 0.022
90% | 2 | $T_{\tilde{h}}$ | 0.00 | 0.269 | 0.268 | -12.7 | 0.019

###  3.7 Discussion and Future Work

In this work, we proposed a zero-estimator approach for improving estimation
of the signal and noise levels explained by a set of covariates in a high-
dimensional regression setting when the covariate distribution is known. We
presented theoretical properties of the naive estimator $\hat{\tau}^{2},$ and
the proposed improved estimators $T_{\hat{h}}$ and $T_{\hat{g}}$. In a
simulation study, we demonstrated that the zero-estimator approach leads to a
significant reduction in the RMSE. Our method does not rely on sparsity
assumptions of the regression coefficient vector, normality of the covariates,
or linearity of $E(Y|X)$. The goal in this work is to estimate the signal
coming from the best linear function of the covariates, which is a model-free
quantity. Our simulations demonstrate that our approach can be generalized to
improve other estimators as well.

We suggest the following directions for future work. One natural extension is
to relax the assumption of known covariate distribution to allow for a more
general setting. This may be studied under the semi-supervised setting where
one has access to a large amount of unlabeled data ($N\gg n$) in order to
obtain theoretical results as a function of $N.$ Another possible future
research might be to extend the proposed approach to generalized linear models
(GLM) such as logistic and Poisson regression, or survival models.

## Supplementary Material

###### Proof of Proposition 3.1.

Let ${W_{i}}={\left({{W_{i1}},...,{W_{ip}}}\right)^{T}}$ and notice that
$\hat{\tau}^{2}=\frac{1}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}^{n}\sum\limits_{j=1}^{p}{{W_{{i_{1}}j}}{W_{{i_{2}}j}}}$
is a U-statistic of order 2 with the kernel
$h_{\tau}\left({{W_{1}},{W_{2}}}\right)=W_{1}^{T}{W_{2}}=\sum\limits_{j=1}^{p}{{W_{1j}}{W_{2j}}}.$

By Theorem 12.3 in van der Vaart (2000),

${\text{Var}}\left({{{\hat{\tau}}^{2}}}\right)=\frac{{4\left({n-2}\right)}}{{n\left({n-1}\right)}}{\zeta_{1}}+\frac{2}{{n\left({n-1}\right)}}{\zeta_{2}},$
(3.22)

where
${\zeta_{1}}\equiv{\text{Cov}}\left[{h\left({{W_{1}},{W_{2}}}\right),h\left({{W_{1}},{{\widetilde{W}_{2}}}}\right)}\right]$
and
${\zeta_{2}}\equiv{\text{Var}}\left[{h\left({{W_{1}},{W_{2}}}\right)}\right],$
and where $\widetilde{W}_{2}$ is an independent copy of $W_{2}$. Define the
$p\times p$ matrix ${\bf{A}}=E\left({{{W}}{W}^{T}}\right)$ and notice that

$\displaystyle{\zeta_{1}}$
$\displaystyle\equiv\text{Cov}\left[{h\left({{{W}_{1}},{{W}_{2}}}\right),h\left({{{W}_{1}},{{\widetilde{W}}_{2}}}\right)}\right]$
$\displaystyle=\sum\limits_{j,j^{\prime}}^{p}{\text{Cov}\left({{W_{1j}}{W_{2j}},{W_{1j^{\prime}}}{{\widetilde{W}}_{2j^{\prime}}}}\right)}=\sum\limits_{j,j^{\prime}}^{p}{\left({{\beta_{j}}{\beta_{j^{\prime}}}E\left[{{W_{1j}}{W_{1j^{\prime}}}}\right]-\beta_{j}^{2}\beta_{j^{\prime}}^{2}}\right)}$
$\displaystyle={\beta^{T}}{\bf{A}}\beta-{\left\|\beta\right\|^{4}},$

and

$\displaystyle{\zeta_{2}}$
$\displaystyle\equiv{\text{Cov}}\left[{h\left({{{W}_{1}},{{W}_{2}}}\right),h\left({{{W}_{1}},{{W}_{2}}}\right)}\right]$
$\displaystyle=\sum\limits_{j,j^{\prime}}{{\text{Cov}}\left({{W_{1j}}{W_{2j}},{W_{1j^{\prime}}}{W_{2j^{\prime}}}}\right)}=\sum\limits_{j,j^{\prime}}\left({{{\left({E\left[{{W_{1j}}{W_{1j^{\prime}}}}\right]}\right)}^{2}}-\beta_{j}^{2}\beta_{j^{\prime}}^{2}}\right)$
$\displaystyle=\left\|{\bf{A}}\right\|_{F}^{2}-{\left\|\beta\right\|^{4}},$

where $\|\textbf{A}\|_{F}^{2}$ is the Frobenius norm of $\textbf{A}.$ Thus,
(2.6) follows from (3.22). ∎

###### Proof of Proposition 2.2.

Since $\hat{\tau}^{2}$ is an unbiased estimator of $\tau^{2},$ then it is
enough to prove that
$\text{Var}[\hat{\tau}^{2}]\xrightarrow{n\rightarrow\infty}0$. By (2.6) it is
enough to require that
$\frac{\beta^{T}\textbf{A}\beta}{n}\xrightarrow{n\rightarrow\infty}0$ and
$\frac{\|\textbf{A}\|^{2}_{F}}{n^{2}}\xrightarrow{n\rightarrow\infty}0$. The
latter is assumed and we now show that the former also holds true.
Let $\lambda_{1{\bf A}}\geq...\geq\lambda_{p{\bf A}}$ be the eigenvalues of A
and notice that A is symmetric. Thus,

$\frac{1}{n^{2}}(\lambda_{1{\bf
A}})^{2}\leq\frac{1}{n^{2}}\sum_{j=1}^{p}\lambda_{j{\bf
A}}^{2}=\frac{1}{n^{2}}tr(\textbf{A}^{2})=\frac{1}{n^{2}}\|\textbf{A}\|_{F}^{2}.$

Since we assume that
$\frac{\|\textbf{A}\|^{2}_{F}}{n^{2}}\xrightarrow{n\rightarrow\infty}0$, we
can conclude that $\frac{\lambda_{1{\bf
A}}}{n}\xrightarrow{n\rightarrow\infty}0$. Now, recall that the maximum of the
quadratic form ${a^{T}}{\bf{A}}a$ satisfies ${\lambda_{1{\bf
A}}}=\mathop{\max}\limits_{a}\,\frac{{{a^{T}}{\bf{A}}a}}{{{{\left\|a\right\|}^{2}}}}.$
Hence,

$\frac{1}{n}\beta^{T}\textbf{A}\beta\equiv\frac{1}{n}\|\beta\|^{2}[(\frac{\beta}{\|\beta\|})^{T}\textbf{A}\frac{\beta}{\|\beta\|}]\leq\frac{1}{n}\|\beta\|^{2}\lambda_{1}=\frac{\lambda_{1}}{n}\tau^{2}\xrightarrow{n\rightarrow\infty}0,$

where the last limit follows from the assumption that $\tau^{2}=O(1),$ and
from the fact that $\frac{\lambda_{1}}{n}\xrightarrow{n\rightarrow\infty}0.$
We conclude that
$\text{Var}[\hat{\tau}^{2}]\xrightarrow{n\rightarrow\infty}0$. ∎

###### Proof of Proposition 3.3.

Recall that (3.8) states that
${\hat{\sigma}^{2}}\equiv\hat{\sigma}_{Y}^{2}-{\hat{\tau}^{2}}$. Thus,

$\displaystyle\text{Var}\left({{{\hat{\sigma}}^{2}}}\right)=\text{Var}\left({{{\hat{\tau}}^{2}}}\right)+\text{Var}\left({\hat{\sigma}_{Y}^{2}}\right)-2\text{Cov}\left(\hat{\sigma}_{Y}^{2},{\hat{\tau}}^{2}\right).$
(3.23)

The variance of $\hat{\tau}^{2}$ is given in (2.6). By standard U-statistic
calculations (see e.g., Example 1.8 in Bose and Chatterjee 2018), the variance
of $\hat{\sigma}_{Y}^{2}$ is

$\displaystyle\text{Var}(\hat{\sigma}_{Y}^{2})=\frac{4(n-2)}{n(n-1)}\psi_{1}+\frac{2}{n(n-1)}\psi_{2},$

where $\psi_{1}\equiv\frac{\mu_{4}-\sigma_{Y}^{4}}{4}$ and
$\psi_{2}\equiv\frac{\mu_{4}+\sigma_{Y}^{4}}{2}.$ This can be further
simplified to obtain

$\displaystyle\text{Var}\left({\hat{\sigma}_{Y}^{2}}\right)=\frac{1}{n}{\mu_{4}}-\frac{{\left({n-3}\right)}}{{n\left({n-1}\right)}}\sigma_{Y}^{4}.$
(3.24)

We now calculate the covariance between and $\hat{\sigma}_{Y}^{2}$ and
$\hat{\tau}^{2}.$ Write,

$\displaystyle\text{Cov}\left({\hat{\sigma}_{Y}^{2},{{\hat{\tau}}^{2}}}\right)=\text{Cov}\left[{\frac{1}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{{{\left({{Y_{{i_{1}}}}-{Y_{{i_{2}}}}}\right)}^{2}}/2,\frac{1}{{n\left({n-1}\right)}}\sum\limits_{j=1}^{p}{\sum\limits_{{i_{1}}\neq{i_{2}}}{{W_{{i_{1}}j}}{W_{{i_{2}}j}}}}}}\right]$
$\displaystyle=\frac{1}{{2{n^{2}}{{\left({n-1}\right)}^{2}}}}\sum\limits_{j=1}^{p}{\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{{i_{3}}\neq{i_{4}}}{\text{Cov}\left[{{W_{{i_{3}}j}}{W_{{i_{4}}j}}},{{\left({{Y_{{i_{1}}}}-{Y_{{i_{2}}}}}\right)}^{2}}\right]}}}.$

The covariance above is different from zero either when one of
$\\{i_{1},i_{2}\\}$ is equal to one of $\\{i_{3},i_{4}\\}$, or when
$\\{i_{1},i_{2}\\}$ is equal to $\\{i_{3},i_{4}\\}$. There are $4n(n-1)(n-2)$
quadruples $(i_{1},i_{2},i_{3},i_{4})$ for the first case, and $2n(n-1)$ for
the second case. Therefore,

$\displaystyle\text{Cov}[\hat{\sigma}_{Y}^{2},\hat{\tau}^{2}]=\frac{2(n-2)}{n(n-1)}\sum_{j=1}^{p}\text{Cov}[W_{1j}{W_{2j}},(Y-Y_{3})^{2}]+\frac{1}{n(n-1)}\sum_{j=1}^{p}\text{Cov}[W_{1j}W_{2j},(Y_{1}-{Y_{2}})^{2}].$
(3.25)

The first covariance term in (3.25) is

$\displaystyle\text{Cov}[W_{1j}W_{2j},(Y_{1}-Y_{3})^{2}]=E\Big{(}[W_{1j}W_{2j}-E(W_{1j})E(W_{2j})][(Y_{1}-Y_{3})^{2}-2\sigma_{Y}^{2}]\Big{)}$
$\displaystyle=E\Big{(}[W_{1j}W_{2j}-\beta_{j}^{2}][(Y_{1}-\alpha-(Y_{3}-\alpha))^{2}-2\sigma_{Y}^{2}]\Big{)}$
$\displaystyle=E\Big{(}[W_{1j}W_{2j}-\beta_{j}^{2}][\\{(Y_{1}-\alpha)^{2}-\sigma_{Y}^{2}\\}-2(Y_{1}-\alpha)(Y_{3}-\alpha)+(Y_{3}-\alpha)^{2}-\sigma_{Y}^{2}]\Big{)}$
$\displaystyle=E\Big{(}[W_{1j}W_{2j}-\beta_{j}^{2}][\\{(Y_{1}-\alpha)^{2}-\sigma_{Y}^{2}\\}]\Big{)}$
$\displaystyle=E\left[{{W_{1j}}{{\left({Y_{1}-\alpha}\right)}^{2}}}\right]\beta_{j}^{2}-\beta_{j}^{2}\sigma_{Y}^{2}.$

The second covariance term of (3.25) is

$\displaystyle\text{Cov}\left[{W_{1j}W_{2j},{{\left(Y_{1}-Y_{2}\right)}^{2}}}\right]=E\left\\{{\left[{W_{1j}W_{2j}-E\left({W_{1j}W_{2j}}\right)}\right]\left[{{{\left(Y_{1}-Y_{2}\right)}^{2}}-E\left\\{{{{\left(Y_{1}-Y_{2}\right)}^{2}}}\right\\}}\right]}\right\\}$
$\displaystyle=E\left\\{{\left({W_{1j}W_{2j}-\beta_{j}^{2}}\right)\left[{{{\left(Y_{1}-Y_{2}\right)}^{2}}-2\sigma_{Y}^{2}}\right]}\right\\}$
$\displaystyle=E\left\\{{\left({W_{1j}W_{2j}-\beta_{j}^{2}}\right)\left[{{{\left\\{{\left({Y_{1}-\alpha}\right)-\left({Y_{2}-\alpha}\right)}\right\\}}^{2}}-2\sigma_{Y}^{2}}\right]}\right\\}$
$\displaystyle=E\left\\{{\left({W_{1j}W_{2j}-\beta_{j}^{2}}\right)\left[{\left\\{{{{\left({Y_{1}-\alpha}\right)}^{2}}-\sigma_{Y}^{2}}\right\\}-2\left({Y_{1}-\alpha}\right)\left({Y_{2}-\alpha}\right)+\left\\{{{{\left({Y_{2}-\alpha}\right)}^{2}}-\sigma_{Y}^{2}}\right\\}}\right]}\right\\}$
$\displaystyle=2E\left\\{{\left[{{{\left({Y_{1}-\alpha}\right)}^{2}}-\sigma_{Y}^{2}}\right]\left({W_{1j}W_{2j}-\beta_{j}^{2}}\right)}\right\\}-2E\left[{\left({W_{1j}W_{2j}-\beta_{j}^{2}}\right)\left({Y_{1}-\alpha}\right)\left({Y_{2}-\alpha}\right)}\right]$
$\displaystyle=2\left\\{{E\left[{{{\left({Y_{1}-\alpha}\right)}^{2}}{W_{1j}}}\right]{\beta_{j}}-\sigma_{Y}^{2}\beta_{j}^{2}}\right\\}-2E\left[{{W_{1j}}\left({Y_{1}-\alpha}\right)W_{2j}\left({Y_{2}-\alpha}\right)}\right]$
$\displaystyle=2E\left[{{{\left({Y_{1}-\alpha}\right)}^{2}}{W_{1j}}}\right]{\beta_{j}}-2\sigma_{Y}^{2}\beta_{j}^{2}-2{\left\\{{E\left[{{W_{1j}}\left({Y_{1}-\alpha}\right)}\right]}\right\\}^{2}},$

Hence, by (3.25) we have

$\displaystyle\text{Cov}\left({\hat{\sigma}_{Y}^{2},{{\hat{\tau}}^{2}}}\right)$
$\displaystyle=\frac{{2\left({n-2}\right)}}{{n\left({n-1}\right)}}\sum\limits_{j=1}^{p}{\left\\{{E\left[{W_{1j}{{\left({Y_{1}-\alpha}\right)}^{2}}}\right]{\beta_{j}}-\beta_{j}^{2}\sigma_{Y}^{2}}\right\\}}$
$\displaystyle+\frac{1}{{n\left({n-1}\right)}}\sum\limits_{j=1}^{p}{\left\\{{2E\left[{{{\left({Y_{1}-\alpha}\right)}^{2}}W_{1j}}\right]{\beta_{j}}-2\sigma_{Y}^{2}\beta_{j}^{2}-2{{\left\\{{E\left[{W_{1j}\left({Y_{1}-\alpha}\right)}\right]}\right\\}}^{2}}}\right\\}},$

which can be further simplified as

$\displaystyle\text{Cov}\left({\hat{\sigma}_{Y}^{2},{{\hat{\tau}}^{2}}}\right)=\frac{2}{n}\left({{\pi^{T}}\beta-{\tau^{2}}\sigma_{Y}^{2}}\right)-\frac{2}{{n\left({n-1}\right)}}\sum\limits_{j=1}^{p}{{{\left\\{{E\left[{W_{1j}\left({Y_{1}-\alpha}\right)}\right]}\right\\}}^{2}}},$
(3.26)

where $\pi={\left({{\pi_{1}},...,{\pi_{p}}}\right)^{T}}$ and ${\pi_{j}}\equiv
E\left[{{{\left({Y_{1}-\alpha}\right)}^{2}}W_{1j}}\right]$.
Plugging (3.24) and (3.26) into (3.23) leads to (3.3). ∎

###### Proof of Corollary 3.1.

Since $\hat{\sigma}^{2}$ is an unbiased estimator of $\sigma^{2},$ it is
enough to prove that
$\text{Var}(\hat{\sigma}^{2})\xrightarrow{n\rightarrow\infty}~{}0.$ Recall
that $\hat{\sigma}^{2}=\hat{\sigma}^{2}_{Y}-\hat{\tau}^{2}$. It follows that
$\text{Var}\left({{{\hat{\sigma}}^{2}}}\right)\leq
2\text{Var}\left({\hat{\sigma}_{Y}^{2}}\right)+2\text{Var}\left({{{\hat{\tau}}^{2}}}\right)$.
Thus, it is enough to prove that
$\text{Var}(\hat{\sigma}_{Y}^{2})\xrightarrow{n\rightarrow\infty}~{}0$ and
$\text{Var}(\hat{\tau}^{2})\xrightarrow{n\rightarrow\infty}~{}0.$ The former
requires, by (3.24), the assumption that $\mu_{4}$ is bounded and the latter
holds true by Proposition 2.2. ∎

######  Remark 3.1.

The condition $\frac{\|{\bf A}\|_{F}^{2}}{n^{2}}\rightarrow 0$ holds in the
homoskedastic linear model with the additional assumption that the columns of
$\bf{X}$ are independent (Livne et al. 2021, Proposition 2). We now show two
more examples where the condition $\frac{\|{\bf
A}\|_{F}^{2}}{n^{2}}\rightarrow 0$ holds without assuming linearity.

1) We show that if $p/n^{2}\rightarrow 0$ and $Y^{2}\leq C$ for some constant
$C$, then $\frac{\|{\bf A}\|_{F}^{2}}{n^{2}}\rightarrow 0$. For
$a\in\mathbb{R}^{p}$ we have,

$\displaystyle{a^{T}}{\bf{A}}a$
$\displaystyle={a^{T}}E\left({{\bf{W}}{{\bf{W}}^{T}}}\right)a=E\left({{a^{T}}{\bf{W}}{{\bf{W}}^{T}}a}\right)=E\left[{{{\left({{a^{T}}{\bf{W}}}\right)}^{2}}}\right]$
$\displaystyle=E\left[{{{\left(Y{\sum\limits_{j=1}^{p}{{a_{j}}{X_{ij}}}}\right)}^{2}}}\right]\leq{C}E\left[{{{\left({\sum\limits_{j=1}^{p}{{a_{j}}{X_{ij}}}}\right)}^{2}}}\right]$
$\displaystyle={C}\left\\{{\sum\limits_{j=1}^{p}{a_{j}^{2}\underbrace{E\left({X_{ij}^{2}}\right)}_{1}}}\right\\}+CE\left({\sum\limits_{j\neq
j^{\prime}}{{a_{j}}{a_{j^{\prime}}}}\underbrace{E\left({{X_{ij}}{X_{j^{\prime}}}}\right)}_{0}}\right)$
$\displaystyle={C}{\left\|a\right\|^{2}},$ (3.27)

where the last equality follows from ${\Sigma}={\bf I}.$ Now, let
$\lambda_{1{\bf A}}\geq\lambda_{2}^{\bf A}\geq,...,\geq\lambda_{p{\bf A}}^{\bf
A}$ be the eigenvalues of $\bf{A},$ and recall that the extrema of the
quadratic form ${a^{T}}{\bf{A}}a$ satisfies ${\lambda_{1{\bf
A}}}=\mathop{\max}\limits_{a}\,\frac{{{a^{T}}{\bf{A}}a}}{{{{\left\|a\right\|}^{2}}}},$
and hence by (3.1) we have $\lambda_{1}\leq C.$ Now, since $p/n^{2}\rightarrow
0$ by assumption, it follows that

$\frac{{\left\|{\bf{A}}\right\|_{F}^{2}}}{{{n^{2}}}}=\frac{{trace\left({{{\bf{A}}^{2}}}\right)}}{{{n^{2}}}}=\frac{{\sum\limits_{j=1}^{p}{\lambda_{j{\bf
A}}^{2}}}}{{{n^{2}}}}\leq\frac{{p\lambda_{1{\bf
A}}^{2}}}{{{n^{2}}}}\leq\frac{{p{C^{2}}}}{{{n^{2}}}}\to 0.$ (3.28)

2) We show that if $p/n^{2}\rightarrow 0$, $E(Y^{4}),E(X_{ij}^{4})\leq C$ for
$j=1,...,p$ and $C\geq 1$, and the columns of $\bf{X}$ are independent, then
$\frac{\|{\bf A}\|_{F}^{2}}{n^{2}}\rightarrow 0$. For $a\in\mathbb{R}^{p}$ we
have by Cauchy–Schwarz,

$\displaystyle{a^{T}}{\bf{A}}a$
$\displaystyle=E\left[{Y_{i}^{2}{{\left({\sum\limits_{j=1}^{p}{{a_{j}}{X_{ij}}}}\right)}^{2}}}\right]\leq{\left[{E\left({Y_{i}^{4}}\right)}\right]^{1/2}}{\left\\{{E\left[{{{\left({\sum\limits_{j=1}^{p}{{a_{j}}{X_{ij}}}}\right)}^{4}}}\right]}\right\\}^{1/2}}$
$\displaystyle\leq{C^{1/2}}{\left\\{{\sum\limits_{{j_{1}}{j_{2}}{j_{3}}{j_{4}}}{{a_{{j_{1}}}}{a_{{j_{2}}}}{a_{{j_{3}}}}{a_{{j_{4}}}}E\left({{X_{i{j_{1}}}}{X_{i{j_{2}}}}{X_{i{j_{3}}}}{X_{i{j_{4}}}}}\right)}}\right\\}^{1/2}}$
$\displaystyle\leq{C^{1/2}}{\left\\{{\sum\limits_{j=1}^{p}{a_{j}^{4}E\left({X_{ij}^{4}}\right)+\sum\limits_{j\neq
j^{\prime}}{a_{j}^{2}a_{j^{\prime}}^{2}E\left({X_{ij}^{2}}\right)E\left({X_{ij^{\prime}}^{2}}\right)}}}\right\\}^{1/2}}$
$\displaystyle\leq{C^{1/2}}{\left\\{{C\sum\limits_{j=1}^{p}{a_{j}^{4}+C\sum\limits_{j\neq
j^{\prime}}{a_{j}^{2}a_{j^{\prime}}^{2}}}}\right\\}^{1/2}}$ $\displaystyle\leq
C{\left\\{{\sum\limits_{j=1}^{p}{a_{j}^{4}+\sum\limits_{j\neq
j^{\prime}}{a_{j}^{2}a_{j^{\prime}}^{2}}}}\right\\}^{1/2}}=C{\left\|a\right\|^{2}}.$

Notice that since that the columns of $\bf{X}$ are independent, the
expectation
$E\left({{X_{i{j_{1}}}}{X_{i{j_{2}}}}{X_{i{j_{3}}}}{X_{i{j_{4}}}}}\right)$ is
not zero (up to permutation) when $j_{1}=j_{2}$ and $j_{3}=j_{4}$ or when
$j_{1}=j_{2}=j_{3}=j_{4}.$ Also notice we obtained the same result as in
(3.1), and hence $\frac{\|{\bf A}\|_{F}^{2}}{n^{2}}\rightarrow 0$ follows by
the same arguments as in the previous example.

######  Remark 3.2.

_Calculations for Equation 3.13_:
Write,

$\displaystyle\text{Cov}\left({{{\hat{\tau}}^{2}},Z_{g}}\right)$
$\displaystyle=\text{Cov}\left({\frac{1}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{{W_{{i_{1}}j}}{W_{{i_{2}}j}}}},\frac{1}{n}\sum\limits_{i=1}^{n}{g(X_{i})}}\right)$
$\displaystyle=\frac{1}{{{n^{2}}\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{\sum\limits_{i=1}^{n}{E\left({{W_{{i_{1}}j}}{W_{{i_{2}}j}}g(X_{i})}\right)}}}$
$\displaystyle=\frac{2}{{{n^{2}}\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{E\left({{W_{{i_{1}}j}}{g_{{i_{1}}}}}\right)E\left({{W_{{i_{2}}j}}}\right)}}$
$\displaystyle=\frac{2}{{{n^{2}}\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{E\left({{W_{{i_{1}}j}}{g_{{i_{1}}}}}\right){\beta_{j}}}}$
$\displaystyle=\frac{2}{n}\beta^{T}\theta_{g},$

where ${\theta_{g}}=E\left[{Wg\left(X\right)}\right]$. Also notice that
$\text{Var}\left(Z_{g}\right)=\text{Var}\left({\frac{1}{n}\sum\limits_{i=1}^{n}{g(X_{i})}}\right)=\frac{{\text{Var}\left({g(X)}\right)}}{n}.$
Thus, by (3.11) we get

${c_{g}^{*}}=\frac{{{\text{Cov}}\left({{{\hat{\tau}}^{2}},Z_{g}}\right)}}{{\text{Var}\left({Z_{g}}\right)}}=\frac{2\beta^{T}\theta_{g}}{{\text{Var}\left({g(X)}\right)}}.$

###### Proof of Proposition 3.4.

We wish to prove that
$\sqrt{n}\left[T_{h}-T_{\hat{h}}\right]\overset{p}{\rightarrow}~{}0.$ Write,

$\sqrt{n}\left({{T_{h}}-{T_{\hat{h}}}}\right)=\sqrt{n}\left[{\left({{T_{h}}-{T_{f}}}\right)+\left({{T_{f}}-{T_{\hat{f}}}}\right)+\left({{T_{\hat{f}}}-{T_{\hat{h}}}}\right)}\right].$

Thus, we need to show that

$\displaystyle\sqrt{n}\left({T_{h}-T_{f}}\right)\overset{p}{\rightarrow}0,$
(3.29)
$\displaystyle\sqrt{n}\left({T_{\hat{f}}-T_{\hat{h}}}\right)\overset{p}{\rightarrow}~{}0,\text{and}$
(3.30)
$\displaystyle\sqrt{n}\left({T_{f}-T_{\hat{f}}}\right)\overset{p}{\rightarrow}0,$
(3.31)

The proofs of (3.29) and (3.30) are essentially the same as the proof of
Proposition 5 in Livne et al. (2021). This is true since
$h(X)\mathbbm{1}_{A}=f(X)\mathbbm{1}_{A},$ where $A$ denotes the event that
the selection procedure $\delta$ perfectly selects the set $\Theta$, i.e.,
$A\equiv\left\\{{{{\bf{S}}_{\delta}}=\Theta}\right\\},$ and $\mathbbm{1}_{A}$
denotes the indicator of $A.$

We now wish to prove (3.31). Write,

$\sqrt{n}\left[T_{\hat{f}}-T_{f}\right]=\sqrt{n}\left[{{{\hat{\tau}}^{2}}-{\hat{c}_{f}^{*}}Z_{f}-\left({{{\hat{\tau}}^{2}}-{c_{f}^{*}}Z_{f}}\right)}\right]={\sqrt{n}}Z_{f}{\left({{{c_{f}^{*}}}-{\hat{c}_{f}^{*}}}\right)}.$

By Markov and Cauchy-Schwarz inequalities, it is enough to show that

$P\left\\{{\left|{\sqrt{n}Z_{f}\left({\hat{c}_{f}^{*}-c_{f}^{*}}\right)}\right|>\varepsilon}\right\\}\leq\frac{{E\left\\{{\left|{\sqrt{n}Z_{f}\left({\hat{c}_{f}^{*}-c_{f}^{*}}\right)}\right|}\right\\}}}{\varepsilon}\leq\frac{{\sqrt{nE\left({Z_{f}^{2}}\right)E\left[{{{\left({\hat{c}_{f}^{*}-c_{f}^{*}}\right)}^{2}}}\right]}}}{\varepsilon}\rightarrow
0.$

Notice that $E(Z_{f}^{2})=\text{Var}(Z_{f})=\frac{\text{Var}[f(X)]}{n}$ and
$E[(\hat{c}_{f}^{*}-c_{f}^{*})^{2}]=\text{Var}(\hat{c}_{f}^{*}),$ where by
(3.4) we have

$\hat{c}_{f}^{*}=\frac{{\frac{2}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{W_{{i_{1}}}^{T}{W_{{i_{2}}}}f\left({{X_{{i_{2}}}}}\right)}}}{{\text{Var}\left[{f\left(X\right)}\right]}}\equiv\frac{U}{\text{Var}[f(X)]}.$

Hence, it enough to show that

$\frac{{\text{Var}\left(U\right)}}{{\text{Var}[f(X)]}}\rightarrow 0.$ (3.32)

The variance of $U$ is

$\displaystyle\text{Var}(U)$
$\displaystyle=\text{Var}\left[{\frac{2}{{n\left({n-1}\right)}}\sum\limits_{{i_{1}}\neq{i_{2}}}{\sum\limits_{j=1}^{p}{{W_{{i_{1}}j}}{W_{{i_{2}}j}}f(X_{i_{2}})}}}\right]$
$\displaystyle=\frac{4}{{{n^{2}}{{\left({n-1}\right)}^{2}}}}\sum\limits_{j,j^{\prime}}^{p}{\sum\limits_{{i_{1}}\neq{i_{2}},{i_{3}}\neq{i_{4}}}{\text{Cov}\left[{{W_{{i_{1}}j}}{W_{{i_{2}}j}}f(X_{i_{2}}),{W_{{i_{3}}j^{\prime}}}{W_{{i_{4}}j^{\prime}}}f(X_{i_{4}})}\right].}}$
(3.33)

The covariance in (Supplementary Material) is different from zero in the
following two cases:

  1. 1.

When $\left\\{{{i_{1}},{i_{2}}}\right\\}$ is equal to
$\left\\{{{i_{3}},{i_{4}}}\right\\}.$

  2. 2.

When one of $\left\\{{{i_{1}},{i_{2}}}\right\\}$ equals to
$\left\\{{{i_{3}},{i_{4}}}\right\\}$ while the other is different.

The first condition includes two different sub-cases and each of those
consists of $n(n-1)$ quadruples $(i_{1},i_{2},i_{3},i_{4})$ that satisfy the
condition. Similarly, the second condition above includes four different sub-
cases and each of those consists of $n(n-1)(n-2)$ quadruples that satisfy the
condition.

We now calculate the covariance for all these six sub-cases.
(1) The covariance when ${i_{1}}={i_{3}},{i_{2}}={i_{4}}$ is

$\displaystyle\text{Cov}\left[{{W_{i_{1}j}}{W_{i_{2}j}}f(X_{i_{2}}),{W_{i_{1}j^{\prime}}}{{W_{i_{2}j^{\prime}}}}f(X_{i_{2}})}\right]$
$\displaystyle=E\left({{W_{i_{1}j}}{W_{i_{1}j^{\prime}}}}\right)E\left[{{W_{i_{2}j}}{{W_{i_{2}j^{\prime}}}}f^{2}(X_{i_{2}})}\right]$
$\displaystyle-E({{W_{i_{1}j}}})E\left[{{W_{i_{2}j}}f(X_{i_{2}})}\right]E\left({{{W_{i_{1}j^{\prime}}}}}\right)E\left[{{{W_{i_{2}j^{\prime}}}}f(X_{i_{2}})}\right]$
$\displaystyle=E\left({{W_{i_{1}j}}{W_{i_{1}j^{\prime}}}}\right)E\left[{{W_{i_{2}j}}{{W_{i_{2}j^{\prime}}}}f^{2}(X_{i_{2}})}\right]-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{fj}}{\theta_{fj^{\prime}}},$

where recall that $b\equiv Wf(X).$ Thus, we define

${\delta_{1}}\equiv\sum\limits_{j,j^{\prime}}{\left\\{{E\left({{W_{i_{1}j}}{W_{i_{1}j^{\prime}}}}\right)E\left[{{W_{i_{2}j}}{W_{i_{2}j^{\prime}}}{f^{2}(X_{i_{2}})}}\right]-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{fj}}{\theta_{fj^{\prime}}}}\right\\}}=E\left({{b^{T}}{\bf{A}}b}\right)-{\left({{\beta^{T}}\theta_{b}}\right)^{2}}.$
(3.34)

(2) The covariance when ${i_{1}}={i_{4}},{i_{2}}={i_{3}}$ is

$\displaystyle\text{Cov}\left[{{W_{i_{1}j}}{W_{i_{2}}}f(X_{i_{2}}),{{W_{i_{2}j^{\prime}}}}{W_{i_{1}j^{\prime}}}f(X_{i_{1}})}\right]$
$\displaystyle=E\left[{{W_{i_{1}j}}{W_{i_{1}j^{\prime}}}f(X_{i_{1}})}\right]E\left[{{W_{i_{2}j}}{{W_{i_{2}j^{\prime}}}}f(X_{i_{2}})}\right]$
$\displaystyle-E(W_{i_{1}j})E\left[{{{W_{i_{2}j}}}f(X_{i_{2}})}\right]E\left({{{W_{i_{2}j^{\prime}}}}}\right)E\left[{{W_{i_{1}j^{\prime}}}f(X_{i_{1}})}\right]$
$\displaystyle={\left\\{{E\left[{{W_{ij}}{W_{ij^{\prime}}}f(X_{i})}\right]}\right\\}^{2}}-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{fj}}{\theta_{fj^{\prime}}}.$

Thus, we define

${\delta_{2}}\equiv\sum\limits_{j,j^{\prime}}{\left\\{{{{\left[{E\left({{W_{ij}}{W_{ij^{\prime}}}f\left(X_{i}\right)}\right)}\right]}^{2}}-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{fj}}{\theta_{fj^{\prime}}}}\right\\}}=\left\|{\bf{B}}\right\|_{F}^{2}-{\left({{\beta^{T}}\theta_{b}}\right)^{2}},$
(3.35)

where ${\bf{B}}\equiv E\left({W{W^{T}}f\left(X\right)}\right).$
(3) Similarly, when $\,{i_{1}}={i_{3}},{i_{2}}\neq{i_{4}}$ we have

${\delta_{3}}\equiv\sum\limits_{j,j^{\prime}}{\left\\{{\left[{E\left({{W_{ij}}{W_{ij^{\prime}}}}\right){\theta_{fj}}{\theta_{fj^{\prime}}}}\right]-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{fj}}{\theta_{fj^{\prime}}}}\right\\}}={\theta_{b}^{T}}{\bf{A}}\theta_{b}-{\left({{\beta^{T}}\theta_{b}}\right)^{2}.}$
(3.36)

(4)-(5) When $i_{1}=i_{4},i_{2}\neq i_{3}$ or when $i_{2}=i_{3},i_{1}\neq
i_{4}$ we have

${\delta_{4}}=\delta_{5}\equiv\sum\limits_{j,j^{\prime}}{\left\\{{{\theta_{fj}}{\beta_{j^{\prime}}}E\left[{{W_{ij}}{W_{ij^{\prime}}}f\left(X_{i}\right)}\right]-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{fj}}{\theta_{fj^{\prime}}}}\right\\}}={\beta^{T}}{\bf{B}}\theta_{b}-{\left({{\beta^{T}}\theta_{b}}\right)^{2}.}$
(3.37)

(6) When $i_{2}=i_{4},i_{1}\neq i_{3}$,

${\delta_{6}}\equiv\sum\limits_{j,j^{\prime}}{\left\\{{{\beta_{j}}{\beta_{j^{\prime}}}E\left[{{W_{ij}}{W_{ij^{\prime}}}{f^{2}(X_{i})}}\right]-{\beta_{j}}{\beta_{j^{\prime}}}{\theta_{fj}}{\theta_{fj^{\prime}}}}\right\\}}={\beta^{T}}{\bf{C}}\beta-{\left({{\beta^{T}}\theta_{b}}\right)^{2}},$
(3.38)

where ${\bf{C}}=E\left({W{W^{T}}f^{2}{{\left(X\right)}}}\right).$ Thus,
plugging-in (3.34) - (3.38) into (Supplementary Material) gives

$\text{Var}\left(U\right)=4{\left\\{{\frac{1}{{n\left({n-1}\right)}}\left({{\delta_{1}}+{\delta_{2}}}\right)+\frac{{\left({n-2}\right)}}{{n\left({n-1}\right)}}\left({{\delta_{3}}+{\delta_{4}}+{\delta_{5}}+{\delta_{6}}}\right)}\right\\}}.$
(3.39)

Recall that we wish to show that
$\frac{\text{Var}(U)}{\text{Var}[f(X)]}\rightarrow 0.$ Thus, it is enough to
show that

${\frac{{\left({{\delta_{1}}+{\delta_{2}}}\right)}}{{{n^{2}\text{Var}[f(X)]}}}}\to
0,$ (3.40)

and

${\frac{{\left({{\delta_{3}}+{\delta_{4}}+{\delta_{5}}+{\delta_{6}}}\right)}}{{{n\text{Var}[f(X)]}}}}\to
0.$ (3.41)

Consider $\delta_{1}.$ For any square matrix ${\bf M},$ we denote
${\lambda_{1{\bf{M}}}}$ to be the largest eigenvalue of ${\bf M}.$ Write,

$\frac{{{\delta_{1}}}}{{{n^{2}}\text{Var}\left[{f\left(X\right)}\right]}}\leq\frac{{E\left({{b^{T}}{\bf{A}}b}\right)}}{{{n^{2}}\text{Var}[f(X)]}}\leq\frac{{{\lambda_{1{\bf{A}}}}}}{n}\frac{{E\left({{{\left\|b\right\|}^{2}}}\right)}}{{n\text{Var}[f(X)]}}\leq\sqrt{\frac{{\left\|{\bf{A}}\right\|_{F}^{2}}}{{{n^{2}}}}}\frac{{E\left({{{\left\|b\right\|}^{2}}}\right)}}{{n\text{Var}[f(X)]}}\to
0,$

where the last inequality holds since we assume that
$\frac{{\left\|{\bf{A}}\right\|_{F}^{2}}}{{{n^{2}}}}\to 0$ and
$\frac{{E\left({{{\left\|b\right\|}^{2}}}\right)}}{{n\text{Var}[f(X)]}}$ is
bounded. Similarly,

$\frac{{{\delta_{2}}}}{{{n^{2}}\text{Var}[f(X)]}}\leq\frac{{\left\|{\bf{B}}\right\|_{F}^{2}}}{{{n^{2}}\text{Var}[f(X)]}}\to
0,$

where $\frac{{\left\|{\bf{B}}\right\|_{F}^{2}}}{{{n^{2}\text{Var}[f(X)]}}}\to
0$ by assumption.
Consider now $\delta_{3}.$ Write,

$\frac{{{\delta_{3}}}}{{n\text{Var}[f(X)]}}\leq\frac{{{\theta_{b}^{T}}{\bf{A}}\theta_{b}}}{{n\text{Var}[f(X)]}}\leq\frac{{{\lambda_{1{\bf{A}}}}}}{n}\frac{{{{\left\|\theta_{b}\right\|}^{2}}}}{{\text{Var}[f(X)]}}\leq\sqrt{\frac{{\left\|{\bf{A}}\right\|_{F}^{2}}}{{{n^{2}}}}}\frac{{{{\left\|\theta_{b}\right\|}^{2}}}}{{\text{Var}[f(X)]}}\to
0,$

where the above holds since we assume that
$\frac{\|\theta_{b}\|^{2}}{\text{Var}[f(X)]}$ is bounded and
$\frac{{\left\|{\bf{A}}\right\|_{F}^{2}}}{{{n^{2}}}}\to 0.$

Recall that $\delta_{4}=\delta_{5}.$ By the Cauchy–Schwarz inequality and the
inequality of arithmetic means we have

$\frac{{{\delta_{4}}}}{{n\text{Var}[f(X)]}}\leq\frac{{{\beta^{T}}{\bf{B}}\theta_{b}}}{{n\text{Var}[f(X)]}}\leq\frac{{\left|{{\lambda_{1{\bf{B}}}}}\right|}}{n}\frac{{\left\|\beta\right\|\left\|\theta_{b}\right\|}}{{\text{Var}[f(X)]}}\leq\sqrt{\frac{{\left\|{\bf{B}}\right\|_{F}^{2}}}{{{n^{2}}}}}\frac{{\left({{{\left\|\beta\right\|}^{2}}+{{\left\|\theta_{b}\right\|}^{2}}}\right)}}{{2n\text{Var}[f(X)]}}\rightarrow
0,$ (3.42)

where the expression in (3.42) converges to zero by similar same arguments as
shown above.
Similarly for $\delta_{6},$

$\frac{{{\delta_{6}}}}{{n\text{Var}[f(X)]}}\leq\frac{{{\beta^{T}}{\bf{C}}\beta}}{{n\text{Var}[f(X)]}}\leq\frac{{{\lambda_{1{\bf{C}}}}}}{n}\frac{{{{\left\|\beta\right\|}^{2}}}}{{\text{Var}[f(X)]}}\leq\sqrt{\frac{{\left\|{\bf{C}}\right\|_{F}^{2}}}{{{\\{n\text{Var}[f(X)]\\}^{2}}}}}\|\beta\|^{2}\to
0,$

where
$\frac{{\left\|{\bf{C}}\right\|_{F}^{2}}}{{\\{n\text{Var}[f(X)]\\}^{2}}}\to 0$
and $\tau^{2}\equiv\|\beta\|^{2}=O(1)$ by assumptions. Hence, (3.32) follows.
This completes the proof that
$\sqrt{n}\left[T_{\hat{h}}-T_{h}\right]\overset{p}{\rightarrow}~{}0.$ ∎

######  Remark 3.3.

We first show that if $Y^{2}$ and $X_{ij}^{2}$ are bounded for all
$j=1,...,p$, $i=1,...,n$ and $p/n=O(1)$, then
$\frac{{E\left({{{\left\|b\right\|}^{2}}}\right)}}{{n\text{Var}[f(X)]}}$ is
bounded. Let $C$ be the upper bound of the maximum of $Y^{2}$ and
$X_{ij}^{2}$, for $j=1,...,p$ and $i=1,...,n.$ Then,

$\frac{{E\left({{{\left\|b\right\|}^{2}}}\right)}}{{n\text{Var}\left[{f\left(X\right)}\right]}}=\frac{{\sum\limits_{j=1}^{p}{E\left({b_{j}^{2}}\right)}}}{{nE\left[{{f^{2}}\left(X\right)}\right]}}=\frac{{\sum\limits_{j=1}^{p}{E\left[{X_{1j}^{2}{Y^{2}}{f^{2}}\left({{X}}\right)}\right]}}}{{nE\left[{{f^{2}}\left(X\right)}\right]}}\leq\frac{{{C^{2}}pE\left[{{f^{2}}\left(X\right)}\right]}}{{nE\left[{{f^{2}}\left(X\right)}\right]}}=\frac{{{C^{2}}p}}{n}=O(1),$

where notice that we used the assumption that $p/n=O(1).$

We now show that under the same assumptions as above, together with the
assumptions that $\Theta$ is bounded and $\text{Var}[f(X)]\geq c>0$, then
$\frac{{\left\|{\bf{B}}\right\|_{F}^{2}}}{{{n^{2}\text{Var}[f(X)]}}}\to 0.$
Recall that
$f\left({{X}}\right)=\sum\limits_{j<j^{\prime}}{{X_{ij}}{X_{ij^{\prime}}}}.$
Notice that when $\Theta$ is bounded, and when the covariates $X_{ij}$,
$i=1,...,n$, $j=1,...,p$ are bounded, then so is $f(X)$. Let $C$ be the upper
bound of $|f(X)|$. Similarly to (3.1), for $a\in\mathbb{R}^{p}$, we have

$\displaystyle{a^{T}}{\bf{B}}a$
$\displaystyle=E\left[{{a^{T}}W{W^{T}}af\left(X\right)}\right]\leq
E\left[{\left|{f\left(X\right)}\right|{{\left({{a^{T}}W}\right)}^{2}}}\right]\leq
CE\left[{{{\left({{a^{T}}XY}\right)}^{2}}}\right]$
$\displaystyle=CE\left[{{Y^{2}}{{\left({\sum\limits_{j=1}^{p}{{a_{j}}{X_{j}}}}\right)}^{2}}}\right]\leq{C^{2}}E\left[{{{\left({\sum\limits_{j=1}^{p}{{a_{j}}{X_{j}}}}\right)}^{2}}}\right]$
$\displaystyle={C^{2}}\sum\limits_{j\neq
j^{\prime}}{{a_{j}}{a_{j^{\prime}}}\underbrace{E\left({{X_{j}}{X_{j^{\prime}}}}\right)}_{0}}+{C^{2}}\sum\limits_{j=1}^{p}{a_{j}^{2}E\left({X_{j}^{2}}\right)}={C^{2}}{\left\|a\right\|^{2}}.$
(3.43)

It follows that
${\lambda_{1{\bf{B}}}}\leq\frac{{{C^{2}}{{\left\|a\right\|}^{2}}}}{{{{\left\|a\right\|}^{2}}}}={C^{2}},$
and by a similar argument as in (3.28) we conclude that

$\frac{{\left\|{\bf{B}}\right\|_{F}^{2}}}{{{n^{2}\text{Var}[f(X)]}}}\leq\frac{{p\lambda_{1{\bf{B}}}^{2}}}{{{n^{2}\text{Var}[f(X)]}}}\leq\frac{pC^{4}}{{{n^{2}c}}}\to
0,$

where recall we assume that $n/p=O(1)$. A similar argument can be used to show
that under the above conditions,
$\frac{{\left\|{\bf{C}}\right\|_{F}^{2}}}{{{\\{n\text{Var}[f(X)]\\}^{2}}}}\to
0.$

######  Remark 3.4.

We use the the following simple selection procedure $\delta$:

Input: A dataset $\left({{{\bf{X}}_{n\times p}},{{\bf{Y}}_{n\times
1}}}\right)$.

  1. 1.

Calculate $\hat{\beta}_{1}^{2},...,\hat{\beta}_{p}^{2}$ where
$\hat{\beta}_{j}^{2}$ is given in (2.4) for $j=1,...,p.$

  2. 2.

Calculate the differences
${\lambda_{j}}=\hat{\beta}_{\left(j\right)}^{2}-\hat{\beta}_{\left({j-1}\right)}^{2}$
for $j=2,\ldots,p$ where
$\hat{\beta}_{\left(1\right)}^{2}<\hat{\beta}_{\left(2\right)}^{2}<...<\hat{\beta}_{\left(p\right)}^{2}$
denotes the order statistics.

  3. 3.

Select the covariates
${{\bf{S}}_{\delta}}=\left\\{{j:\hat{\beta}_{\left(j\right)}^{2}>\hat{\beta}_{\left({{j^{*}}}\right)}^{2}}\right\\}$,
where ${j^{*}}=\mathop{\arg\max}\limits_{j}{\lambda_{j}}$.

Output: Return ${{\bf{S}}_{\delta}}$.

Algorithm 7 Covariate selection $\delta$

The algorithm above finds the largest gap between the ordered estimated
squared coefficients and then uses this gap as a threshold to select a set of
coefficients ${{\bf{S}}_{\delta}}\subset\left\\{{1,...,p}\right\\}.$ The
algorithm works well in scenarios where a relatively large gap truly separates
between larger coefficients and the smaller coefficients of the vector
$\beta$.

##  4 Real Data Application

In this chapter, we use the methods proposed in the previous chapters to
estimate the signal level $\tau^{2}$ on different datasets. In real-world data
(on opposed to synthetic data), it is common to encounter undesired issues
such as error of measurements, outliers, highly-skewed covariates, collinear
covariates, violation of the i.i.d. assumption, etc. In practice, these issues
might impact the performance of the proposed estimators. Moreover, the
theoretical results provided is this thesis are based on the underlying
assumption that the covariates’ distribution is known, while in practice one
needs to estimate it by using the unlabeled data. Our aim now in this chapter
is to investigate how well our proposed estimators perform on real datasets.

We provide below a brief description of the four real datasets and the
synthetic data we used in our simulations in this chapter. Data pre-possessing
steps that were applied are also described.

  1. 1.

The Blog Feedback dataset, which is publicly available at the UCI machine
learning repository (https://archive.ics.uci.edu/ml/datasets/BlogFeedback),
contains data from Hungarian blog sites. The covariates include blog posts
characteristics such as: the number of comments in different time intervals
relative to the publication of the blog post, and the length of the blog post.
The response variable is the number of comments that a given blog post
received within the first 24 hours after publication. A log transformation was
applied to the response variable. The dataset includes $N=122,993$ rows and
281 columns. After we removed collinear columns the resulting number of
covariates was $p=248.$

  2. 2.

The House Price dataset, which is publicly available at Kaggle
(https://www.kaggle.com/datasets/egebozoglu/house-price-linear-regression),
contains different house characteristics such as: number of bedrooms, square
feet of different areas in the house, the number of years since renovation.
The response variable is the house price. The dataset includes $N=12,613$ rows
and 21 columns. We also included all the pairwise interaction terms. The
resulting number of covariates after removing collinear columns is $p=138.$

  3. 3.

The Life Expectancy dataset, which is publicly available from Kaggle
(https://www.kaggle.com/datasets/kumarajarshi/life-expectancy-who), consists
of health, economic and social factors that affect life expectancy. The data
was collected by the World Health Organization (WHO) between the years 2000 to
2015 for 193 countries. The covariates include information such as
immunization level for Polio, the gross domestic product (GDP) and alcohol
consumption level. The response variable is life expectancy in a given country
in a given year. The dataset includes $N=2,928$ rows and 22 columns. We also
included all the pairwise interaction terms. The resulting number of
covariates after removing collinear columns is $p=190.$

  4. 4.

The Superconductivty dataset, which is publicly available at the UCI machine
learning repository
(https://archive.ics.uci.edu/ml/datasets/superconductivty+data), contains
characteristics of superconducting materials (materials that conduct current
with zero resistance). The covariates are based on properties such as atomic
mass and radius, electron affinity and fusion heat. The response variable is
the superconducting critical temperature. A log transformation was applied to
the response variable. The dataset includes $N=21,263$ rows and $81$ columns.
We also included the pairwise interactions terms of the top ten covariates
that have the largest absolute value of standardized coefficient (i.e., the
”t-value”) from running a linear regression model. The resulting number of
covariates after removing collinear columns is $p=126.$

  5. 5.

A synthetic dataset was simulated similarly to the simulation study described
in Chapter 2.5 with $p=200,\tau^{2}=2$ and $\eta=0.5.$

Recall that the parameter of interest is $\tau^{2}=\beta^{T}\Sigma_{X}\beta$,
where $\beta$ is best linear predictor coefficient vector from Chapter 3.1 and
$\Sigma_{X}$ is the covariance matrix of $X.$ In each dataset, we approximate
the parameter $\tau^{2}$ by $\hat{\beta}^{T}\hat{\Sigma}_{X}\hat{\beta}$ where
$\hat{\beta}$ and $\hat{\Sigma}_{X}$ are the least squares estimator and the
plug-in estimator over the entire dataset, respectively. We consider this
approximation as the true value of $\tau^{2}$.

In order to asses the performance of the different estimators, we run 100
independent simulations for each dataset. In each simulation run, we randomly
take subsamples of $n=p$ or $n=p/2$ observations
$(X_{1},Y_{1}),...,(X_{n},Y_{n})$ from the $N$ observations that exist in the
dataset. We treat these $n$ observations as the labeled data while the other
$N-n$ observations are considered the unlabeled data, i.e., the response
variable for this part was omitted. We then apply three initial estimators:
the Naive estimator from Chapter 2.2.2, the EigenPrism procedure (Janson et
al. 2017) and the PSI estimator (Taylor and Tibshirani 2018); see Chapter 2.5.
For each initial estimator we compute the two additional proposed estimators
that aim to improve the initial ones. We denote them as the Single and the
Selection estimators, where the LASSO (Hastie et al. 2015, Chapter 2.2) was
used as the covariate-selection procedure in this chapter. For more details
about the above estimators see Chapter 3.6, where the Single and the Selection
estimators are denoted there by $T_{\hat{g}}$ and $T_{\hat{h}}$, respectively.
We calculated the mean and the MSE from the $100$ independent runs by
$\frac{1}{{100}}\sum\limits_{k=1}^{100}{\hat{\tau}_{k}^{2}}$ and
$\frac{1}{{100}}\sum\limits_{k=1}^{100}{{{\left({\hat{\tau}_{k}^{2}-{{\tau}^{2}}}\right)}^{2}}}$,
respectively, where $\hat{\tau}^{2}_{k}$ is an estimate of $\tau^{2}$ from the
$k^{th}$ run.

Tables 8-12 show the mean, the MSE and the relative change in the MSE for the
different estimators based on the simulation runs. Standard errors are given
in parenthesis. Table 13 shows the approximated correlations between the
different zero-estimators and the initial estimators. In order the approximate
these correlations, we first randomly drew 300 independent subsamples of size
$n=p$ from each dataset. Second, for each subsample, we calculated the initial
estimators of the signal level and the different zero-estimators. We then
calculated the empirical correlations, which are shown in Table 13.

Important points to notice are:

  * •

For the Naive estimator, we see that both of the proposed estimators
successfully improve the MSE in all datasets and scenarios.

  * •

For the EigenPrism estimator, the Selection estimator improves the MSE in all
datasets, while the Single estimator improves it in the synthetic dataset, the
House dataset and the Superconductivty dataset, but not in the Blog Feedback
and the Life Expectancy datasets. In the latter two, the single estimator in
fact resulted in a slightly larger MSE in some scenarios.

  * •

For the PSI estimator, the performance of the proposed estimators generally
did not improve the MSE, and in some cases the MSE was higher compared to the
original PSI. Notice, however, that in all the datasets PSI yielded
underbiased estimates since the estimator relies on sparsity assumptions,
which are likely to not hold in these datasets; see Chapter 2.5.

  * •

In general, we noticed that the larger the correlation between a given initial
estimator and a zero-estimator, the better we were able to reduce the
variance, as expected from the theory; see (3.4). For example, Table 13 shows
that in the Blog Feedback dataset, the zero-estimator that is based on a
selection procedure, is more correlated with the naive and EigenPrism
estimators than the zero-estimator that is not based on a selection procedure.
This aligns with the simulation results shown in Table 8, which demonstrate
that the improvement of the Selection estimator was generally better than the
improvement of the Single estimator for this dataset.

To a practitioner who wishes to apply the proposed estimators, we recommend
keeping in mind that, in general, the Selection estimator tends to work better
when the data is sparse, while the Single estimator tends to work better when
the data is dense. Thus, in theory, a practitioner who knows the level of
sparsity in advance could benefit from choosing the more appropriate zero-
estimator. Note that, as shown in the simulations of this chapter, the
proposed estimators do not guarantee to improve any given initial estimator
(e.g., the PSI estimator) and there is a risk it will increase the MSE.
However, the more correlation exists between a zero-estimator and an initial
estimator, the more likely it is to create an improvement. Hence, in order to
implement our methods one needs to find a zero estimator that is highly
correlated with the initial estimator. The question of how to find a highly
correlated zero-estimator is not discussed here and remains open.

Table 8: The Blog Feedback dataset.
Summary statistics for the different estimators. Mean, mean square error (MSE)
and percentage change from the initial estimator (in terms of MSE) are shown.
Simulation standard errors are shown in parenthesis. The table results were
computed over $100$ subsamples for each setting; $p=248$, $\tau^{2}=0.609.$

n | p | Estimator | Mean | MSE | % Change
---|---|---|---|---|---
124 | 248 | Naive | 0.6195 (0.066) | 0.427 (0.107) | 0
124 | 248 | Selection | 0.5631 (0.058) | 0.3304 (0.069) | -22.623
124 | 248 | Single | 0.6322 (0.064) | 0.4088 (0.107) | -4.262
248 | 248 | Naive | 0.6215 (0.034) | 0.1122 (0.022) | 0
248 | 248 | Selection | 0.5909 (0.031) | 0.0934 (0.019) | -16.756
248 | 248 | Single | 0.6225 (0.033) | 0.1109 (0.022) | -1.159