This is the deformed Ferrara–Sabharwal metric appearing in [CDS17, Corollary
15 (p. 97)] upon taking

$\displaystyle\rho^{\mathrm{CDS}}=r^{2}-2k,\qquad
w^{\mathrm{CDS}}_{i}=\frac{(-1)^{\delta_{i,0}}}{\sqrt{2}}\,\overline{w}_{i},\qquad
t^{\mathrm{CDS}}=4t.$ (5.31)

This was originally derived by physicists in [CFG89, FS90, RSV06] as a string-
theoretic moduli space metric with perturbative quantum corrections.

When $k=0$, this is a left invariant metric on (a $\mathbb{Z}$-quotient of)
the non-compact Wolf space
$\mathrm{SU}(n+1,2)/\mathrm{S}(\mathrm{U}(n+1)\times\mathrm{U}(2))$.

###  5.3 Heisenberg action

While there is always a local basis of $\nabla$-parallel sections of $E$, a
global basis may fail to exist because of the presence of monodromy. We can
however get around this issue by considering the universal cover of $M$ and
pulling back $E$ along the covering map. This gives us a VPHS whose Gauß–Manin
connection has no monodromy. We will show that in the absence of monodromy,
the supergravity c-map admits the action of the Heisenberg group.

In particular, we would like to obtain the Killing vector fields on $N_{2k}$
that generate the Heisenberg action. In order to do this, we begin by
constructing vertical Killing vector fields on the _rigid_ c-map space. The
following proposition will show that the bundle map $\mathrm{vert}\colon
E_{\pi}\to T(E_{\pi})$ maps pullbacks of parallel sections of $E$ to Killing
vector fields of $E_{\pi}$.

######  Proposition 5.7.

Let $s\in\Gamma(E)$ be such that $\nabla s=0$. Then, the vertical vector field
$v_{s}\in\Gamma(T(E_{\pi}))$ corresponding to $s$ is Killing and
$\omega_{\mathrm{H}}$-Hamiltonian with Hamiltonian function

$f_{s}=Q(s,\Phi_{E}).$ (5.32)

Moreover, it preserves $Z$, $f_{\mathrm{H}}$, $I_{j}$, and
$\widehat{\omega}_{j}=\widehat{g}(I_{j}\cdot,\cdot)$, for $j=1,2,3$ as defined
in (5.10) and (5.5).

###### Proof.

First of all, notice that since $s$ is $\nabla$-parallel,
$\nabla(\nabla_{v_{s}}\Phi_{E})=\nabla s=0$. We now compute the Lie derivative
of $\widehat{g}$ along $v_{s}$, obtaining

$\displaystyle\begin{split}&\mathcal{L}_{v_{s}}\widehat{g}(X,Y)=\mathcal{L}_{v_{s}}\big{(}\widetilde{g}+\operatorname{Re}h_{\mathrm{G}}(\nabla\Phi_{E},\nabla\Phi_{E})\big{)}(X,Y)\\\
&\qquad=\mathcal{L}_{v_{s}}\big{(}\operatorname{Re}h_{\mathrm{G}}(\nabla\Phi_{E},\nabla\Phi_{E})\big{)}(X,Y)\\\
&\qquad=v_{s}\left(\operatorname{Re}h_{\mathrm{G}}(\nabla_{X}\Phi_{E},\nabla_{Y}\Phi_{E})\right)\\\
&\qquad\quad-\operatorname{Re}h_{\mathrm{G}}(\nabla_{[v_{s},X]}\Phi_{E},\nabla_{Y}\Phi_{E})-\operatorname{Re}h_{\mathrm{G}}(\nabla_{X}\Phi_{E},\nabla_{[v_{s},Y]}\Phi_{E})\\\
&\qquad=\operatorname{Re}(\nabla_{v_{s}}h_{\mathrm{G}})(\nabla_{X}\Phi_{E},\nabla_{Y}\Phi_{E})\\\
&\qquad\quad+\operatorname{Re}h_{\mathrm{G}}\left([\nabla_{v_{s}},\nabla_{X}]\Phi_{E}-\nabla_{[v_{s},X]}\Phi_{E},\nabla_{Y}\Phi_{E}\right)\\\
&\qquad\quad+\operatorname{Re}h_{\mathrm{G}}\left([\nabla_{v_{s}},\nabla_{Y}]\Phi_{E}-\nabla_{[v_{s},Y]}\Phi_{E},\nabla_{X}\Phi_{E}\right)\\\
&\qquad=\operatorname{Re}(\nabla_{v_{s}}h_{\mathrm{G}})(\nabla_{X}\Phi_{E},\nabla_{Y}\Phi_{E}).\end{split}$
(5.33)

The last expression vanishes because $\nabla_{v_{s}}h_{G}$ is to be
interpreted as $\nabla_{(\pi_{E})_{*}v_{s}}h_{G}$ by the definition of a
pullback connection, and $v_{s}$ is a vertical tangent vector on the total
space of the bundle $\pi_{E}\colon E\to M$.

We can now check that $f_{s}$ is a Hamiltonian function of $v_{s}$ by
computing

$\displaystyle\mathrm{d}f_{s}$
$\displaystyle=\nabla(Q(s,\Phi_{E}))=Q(s,\nabla\Phi_{E})=Q(\nabla_{v_{s}}\Phi_{E},\nabla\Phi_{E})=-\iota_{v_{s}}\omega_{\mathrm{H}}.$
(5.34)

The function $r$ is defined on the base, and therefore it is invariant along
vertical vector fields. In particular $\mathcal{L}_{v_{s}}f_{\mathrm{H}}=0$.

The vector field $Z$ is also invariant, as
$\widetilde{\varphi}=r^{-2}\widehat{g}(Z,\cdot)$ is

$\displaystyle\mathcal{L}_{v_{s}}\widetilde{\varphi}$
$\displaystyle=\mathrm{d}\iota_{v_{s}}\widetilde{\varphi}+\iota_{v_{s}}\mathrm{d}\widetilde{\varphi}=0-2\iota_{v_{s}}\widetilde{\omega}=0.$
(5.35)

The invariance of $I_{k}$ is equivalent to that of $\widehat{\omega}_{k}$.
Note that $\widehat{\omega}_{1}=2\widetilde{\omega}+\omega_{\mathrm{H}}$ is a
sum of invariant $2$-forms. Their invariance follows from the observation that
$\widetilde{\omega}$ is horizontal, and $v_{s}$ is
$\omega_{\mathrm{H}}$-Hamiltonian with $\omega_{\mathrm{H}}$ closed.

Meanwhile, for $\widehat{\omega}_{2}$, we compute

$\displaystyle\begin{split}\mathcal{L}_{v_{s}}\widehat{\omega}_{2}&=\mathrm{d}\iota_{v_{s}}\widehat{\omega}_{2}=\mathrm{d}\big{(}\widehat{g}(I_{2}v_{s},\cdot)\big{)}=-\mathrm{d}\Big{(}\widehat{g}\big{(}(\phi^{-1}(s))^{\nabla},\cdot\big{)}\Big{)}\\\
&=-\mathrm{d}\Big{(}\widetilde{g}\big{(}\phi^{-1}(s),\cdot\big{)}\Big{)}=-2\mathrm{Alt}\Big{(}\nabla^{\widetilde{g}}\widetilde{g}\big{(}\phi^{-1}(s),\cdot\big{)}\Big{)}\\\
&=-2\mathrm{Alt}\Big{(}\widetilde{g}\big{(}\nabla^{\widetilde{g}}(\phi^{-1}(s)),\cdot\big{)}\Big{)}=4\mathrm{Alt}\Big{(}\widetilde{g}\big{(}\operatorname{Re}\widetilde{\eta}_{\cdot}(\phi^{-1}(s)),\cdot\big{)}\Big{)}.\end{split}$
(5.36)

This quantity vanishes because $\operatorname{Re}\widetilde{\eta}$ is
symmetric. The invariance of $I_{3}=I_{1}I_{2}$ follows. ∎

Killing vector fields on a hyperkähler manifold satisfying the properties
resulting from Proposition 5.7 can be used to obtain Killing vector fields on
the corresponding quaternionic Käher manifold via the following proposition.

######  Proposition 5.8 ([CST21] Proposition 3.5 (p. 108)).

Let $X$ be a Killing vector field on $(E_{\pi},\widehat{g})$ which is
$\omega_{\mathrm{H}}$-Hamiltonian with Hamiltonian function $f_{X}$, and which
preserves the hyperkähler structure $I_{1}$, $I_{2}$, $I_{3}$ and the twist
data $Z$, $f_{\mathrm{H}}$. Then, the vector field

$X^{\mathrm{Q}}\colonequals\mathrm{tw}\left(X-\frac{f_{X}}{f_{\mathrm{H}}}Z\right)$
(5.37)

is Killing with respect to $g_{2k}$.

Let $s$ thus be a section of $E$ such that $\nabla s=0$, and let $v_{s}$ be
the corresponding vertical vector field on $E_{\pi}$. We obtain a Killing
vector field

$w_{s}\colonequals\mathrm{tw}\left(v_{s}-\frac{f_{s}}{f_{\mathrm{H}}}Z\right).$
(5.38)

Notice that $-\widetilde{\varphi}(v_{s})=0$, so by $Z$-invariance, $v_{s}$ is
the $(-\widetilde{\varphi})$-horizontal lift of some vector field on
$\mathbb{R}_{>2k}\times M$. Let $v^{\prime}_{s}$ be the
$\varphi_{k}$-horizontal lift of said vector field on $N_{2k}$. The twist of
$v_{s}$ is then $v^{\prime}_{s}-\beta(v_{s})Z_{k}$, and that of $Z$ is
$-f_{\mathrm{H}}Z_{k}$, where $Z_{k}$ is the fundamental vector field of the
$\mathrm{U}(1)$-bundle $N_{2k}\to E_{\mathrm{pr}_{M}}$. Therefore,

$\displaystyle\begin{split}w_{s}&=v^{\prime}_{s}-\beta(v_{s})Z_{k}+f_{s}Z_{k}=v^{\prime}_{s}+\frac{1}{2}\,Q(\Phi_{E},\nabla_{v_{s}}\Phi_{E})Z_{k}+Q(s,\Phi_{E})Z_{k}\\\
&=v^{\prime}_{s}+\frac{1}{2}\,Q(\Phi_{E},s)Z_{k}+Q(s,\Phi_{E})Z_{k}=v^{\prime}_{s}+\frac{1}{2}\,Q(s,\Phi_{E})Z_{k}.\end{split}$
(5.39)

Now, in order to see the Heisenberg action, we make use of the following
result relating commutators of Killing vector fields on the quaternionic
Kähler manifold to those on the hyperkähler manifold.

######  Proposition 5.9.

Let $X_{1},X_{2}$ be Killing vector fields on $(E_{\pi},\widehat{g})$ as in
Proposition 5.8 with Hamiltonian functions $f_{X_{1}},f_{X_{2}}$. Then the
corresponding Killing vector fields $X_{1}^{\mathrm{Q}}$ and
$X_{2}^{\mathrm{Q}}$ on $N_{2k}$ satisfy

$[X_{1}^{\mathrm{Q}},X_{2}^{\mathrm{Q}}]=\mathrm{tw}([X_{1},X_{2}])+\mathrm{tw}(\omega_{\mathrm{H}}(X_{1},X_{2}))Z_{k}.$
(5.40)

###### Proof.

See [CST21, Theorem 3.8 (p. 109)] combined with Proposition 5.8. ∎

######  Corollary 5.10.

The vertical Killing vector fields $w_{s}$ associated to $\nabla$-parallel
sections $s$ of $E$ and the Killing vector field $Z_{k}$ satisfy the
Heisenberg algebra relations

$\displaystyle[w_{s},Z_{k}]=0,$
$\displaystyle[w_{s_{1}},w_{s_{2}}]=Q(s_{1},s_{2})Z_{k}.$ (5.41)

######  Remark 5.11.

The vertical Killing vector field $w_{s}$ associated to a $\nabla$-parallel
section $s$ of $E$ can be readily integrated to yield an isometry $\psi_{s}$
of $N_{2k}$, given by

$\psi_{s}(p,x_{p},u_{p})=\left(p,x_{p}+s_{p},u_{p}\exp\bigg{(}\frac{1}{2}\,Q_{p}(s_{p},x_{p})\bigg{)}\right).$
(5.42)

Here, $p$ is a point on $M\times\mathbb{R}_{>2k}$, while $x_{p}$ and $u_{p}$
are points in the fibres over $p$ of the respective pullbacks of $E\rightarrow
M$ and $S_{k}\rightarrow M$ to $M\times\mathbb{R}_{>2k}$.

We can now explicitly see that this gives an action of the Heisenberg group:

$\begin{split}\psi_{s_{1}}\circ\psi_{s_{2}}(p,x_{p},u_{p})&=\psi_{s_{1}+s_{2}}\bigg{(}p,x_{p},u_{p}\exp\bigg{(}\frac{1}{2}\,Q(s_{1},s_{2})_{p}\bigg{)}\bigg{)}.\end{split}$
(5.43)

###  5.4 Functoriality

The supergravity c-map is a natural construction. This means that isomorphisms
of PSK manifolds lift to isometries of quaternionic Kähler manifolds in a
natural way. In this subsection, we describe this lifting of isomorphisms
explicitly.

Consider the PSK isomorphism $(\widetilde{\psi},\psi)$ from
$(\pi\colon\widetilde{M}\rightarrow
M,\widetilde{g},I,\widetilde{\omega},\widetilde{\nabla},\xi)$ to
$(\pi^{\prime}\colon\widetilde{M^{\prime}}\rightarrow
M^{\prime},\widetilde{g}^{\prime},I^{\prime},\widetilde{\omega}^{\prime},\widetilde{\nabla}^{\prime},\xi^{\prime})$.
In the following discussion, any construction associated to the PSK manifold
$M$ will have a primed counterpart associated to the PSK manifold
$M^{\prime}$. Using the structure of
$L\colonequals\widetilde{M}\times_{\mathbb{C}^{\times}}\mathbb{C}$ and
$L^{\prime}$ as associated fibre bundles and the isomorphism
$\widetilde{\psi}\colon\widetilde{M}\rightarrow\widetilde{M}^{\prime}$, we can
produce a complex vector bundle isomorphism
$\psi^{L}\colonequals\widetilde{\psi}\times_{\mathbb{C}^{\times}}\mathrm{id}_{\mathbb{C}}\colon
L\rightarrow L^{\prime}$ covering $\psi$. This can be regarded as an extension
of $\widetilde{\psi}$. Using $\psi^{L}$, we can then construct the maps
$\psi^{E_{\mathbb{C}}}\colon E_{\mathbb{C}}\rightarrow
E^{\prime}_{\mathbb{C}}$ and $\psi_{\mathrm{pr}}^{E}\colon
E\times\mathbb{R}_{>2k}\rightarrow E^{\prime}\times\mathbb{R}_{>2k}$ via

$\begin{split}\psi^{E_{\mathbb{C}}}&=\psi^{L}+\psi^{L}\otimes(\psi^{1,0}_{*})_{\psi^{-1}}+\overline{\psi}^{L}\otimes(\psi^{0,1}_{*})_{\psi^{-1}}+\overline{\psi}^{L},\\\
\psi_{\mathrm{pr}}^{E}&=\left.\psi^{E_{\mathbb{C}}}\right|_{E}\times\mathrm{id}_{\mathbb{R}_{>2k}}.\end{split}$
(5.44)

Since $\widetilde{\psi}$ is compatible with the Kähler structures on
$\widetilde{M}$ and $\widetilde{M}^{\prime}$, the map $\psi^{L}$ is compatible
with the Hermitian forms $h_{L}$ and $h_{L^{\prime}}$, and with the Chern
connections $\chi$ and $\chi^{\prime}$ (whose imaginary parts are
$\widetilde{\varphi}$ and $\widetilde{\varphi}^{\prime}$ respectively). Since
$\widetilde{\psi}_{*}$ maps $\zeta$ to $\zeta^{\prime}$, the map $\psi^{L}$
sends the unit circle bundle $S\subset L$ to $S^{\prime}\subset L^{\prime}$.
So, we obtain a map $\psi^{S}_{k}\colon S_{k}\rightarrow S^{\prime}_{k}$
pulling back the principal $\mathrm{U}(1)$-connection $\varphi^{\prime}_{k}$
to $\varphi_{k}$ via the codifferential $(\psi^{S}_{k})^{*}$.

By Theorem 5.4, the quaternionic Kähler manifold $N_{2k}$ obtained from $M$
via the supergravity c-map is the fibred product of the bundles
$E\times\mathbb{R}_{>2k}\rightarrow M$ and $S_{k}\rightarrow M$. An analogous
statement holds for $N^{\prime}_{2k}$. Therefore, the maps
$\psi_{\mathrm{pr}}^{E}$ and $\psi^{S}_{k}$ (both of which cover $\psi$)
canonically induce a map $\psi^{N}\colon N_{2k}\rightarrow N^{\prime}_{2k}$
covering $\psi$.

As the map $\widetilde{\psi}$ relates $\widetilde{\omega}$ and
$\widetilde{\nabla}$ to $\widetilde{\omega}^{\prime}$ and
$\widetilde{\nabla}^{\prime}$, the induced map $\psi_{\mathrm{pr}}^{E}$
relates the (pullbacks of the) bilinear form $Q$ and the connection $\nabla$
on the bundle $E\times\mathbb{R}_{>2k}\rightarrow M\times\mathbb{R}_{>2k}$, to
$Q^{\prime}$ and $\nabla^{\prime}$ on
$E^{\prime}\times\mathbb{R}_{>2k}\rightarrow
M^{\prime}\times\mathbb{R}_{>2k}$. Similarly, the Kähler structure of $M$ and
the coordinate $r$ are related to the corresponding counterparts on
$M^{\prime}$. Therefore, the map $\psi^{N}$ is an isometry between the
quaternionic Kähler manifolds $(N_{2k},g_{2k})$ and
$(N^{\prime}_{2k},g^{\prime}_{2k})$. As a result of naturality, the assignment
$(\psi,\widetilde{\psi})\mapsto\psi^{N}$ is a functor. This provides us with
the desired functorial lift of PSK isomorphisms.

In particular, we can consider $1$-parameter subgroups of PSK automorphisms
and take the differential of the action at the identity to get vector fields
that act as infinitesimal automorphisms. Lifts of infinitesimal PSK
automorphisms to infinitesimal isometries (i.e. Killing vector fields) of the
supergravity c-map metric were constructed in a different way in [CST21,
Theorem 3.15 (p. 112)]. Unsurprisingly, the above functorial construction
reproduces the same result.

######  Proposition 5.12.

Let $\widetilde{M}\rightarrow M$ be a PSK manifold and let $\widetilde{X}$ and
$X$ be vector fields on $\widetilde{M}$ and $M$ respectively, which generate a
$1$-paramater family of PSK automorphisms $(\widetilde{\psi}_{t},\psi_{t})$ on
$\widetilde{M}\rightarrow M$. Then, the complete lift $\widetilde{X}^{T}$ on
$T\widetilde{M}\cong E_{\pi}$ is $\omega_{\mathrm{H}}$-Hamiltonian with
Hamiltonian function

$f^{T}_{X}\colonequals
f_{\widetilde{X}^{T}}=-\Big{(}f_{\mathrm{H}}+k\Big{)}\widetilde{\varphi}(X)-\frac{1}{2}\,Q(\Phi_{E},\nabla_{X}\Phi_{E}).$
(5.45)

Moreover, the functorial lift $\psi_{t}^{N}$ of
$(\widetilde{\psi}_{t},\psi_{t})$ is generated by the Killing vector field
$(\widetilde{X}^{T})^{\mathrm{Q}}$ on $N_{2k}$ arising from
$\widetilde{X}^{T}$ and the choice of $\omega_{\mathrm{H}}$-Hamiltonian
function $f^{T}_{X}$, in accordance with (5.37).

###### Proof.

In order to prove the first part, our strategy will be to use Lemma 4.4. We
already know $\omega_{\mathrm{H}}$ is non-degenerate. We need to show that
$\widetilde{X}^{T}$, interpreted as vector field tangent to $E_{\pi}$,
preserves

$\displaystyle
f_{\mathrm{H}}=-\frac{1}{2}(r^{2}+2k),\qquad\varphi=-\widetilde{\varphi},\qquad\beta=-\frac{1}{2}\,Q(\Phi_{E},\nabla\Phi_{E}).$
(5.46)

For this purpose, we will describe the flow of $\widetilde{X}^{T}$ on
$E_{\pi}$ and show that it preserves these objects.

Consider first a PSK automorphism $(\widetilde{\psi},\psi)$. By the universal
property of the pullback we get a map
$\psi^{(E_{\mathbb{C}})_{\pi}}\colon(E_{\mathbb{C}})_{\pi}\to E_{\pi}$, which
in turn restricts to a map $\psi^{E_{\pi}}\colon E_{\pi}\to E_{\pi}$, as
$\psi^{E_{\mathbb{C}}}$ preserves the real structure. Given $Y\in
T\widetilde{M}_{>2k}$, we determine the element of $E_{\pi}$ corresponding to
$\psi_{*}Y$ to be

$\begin{split}\phi(\widetilde{\psi}_{*}Y)&=\chi(\widetilde{\psi}_{*}Y)\Phi_{L}+\Phi_{L}\otimes\pi_{*}\widetilde{\psi}_{*}Y^{1,0}+\overline{\Phi}_{L}\otimes\pi_{*}\widetilde{\psi}_{*}Y^{0,1}+\overline{\chi}(\psi_{*}Y)\overline{\Phi}_{L}\\\
&=\chi(Y)\Phi_{L}+\Phi_{L}\otimes\psi_{*}\pi_{*}Y^{1,0}+\overline{\Phi}_{L}\otimes\psi_{*}\pi_{*}Y^{0,1}+\overline{\chi}(Y)\overline{\Phi}_{L}.\end{split}$
(5.47)

Here, we used the fact that $\chi$ and $\overline{\chi}$ are invariant under
the action of the diffeomorphism $\widetilde{\psi}$, and that
$\widetilde{\psi}$ is a map of the bundle $\pi\colon\widetilde{M}\rightarrow
M$ covering $\psi$. We know that the tautological section $\Phi_{L}$ is
preserved by the map (induced by) $\psi^{L}$. Therefore,
$\phi(\widetilde{\psi}_{*}Y)$ is actually equal to
$\psi^{(E_{\mathbb{C}})_{\pi}}(\phi(Y))$, implying the commutativity relation
$\psi^{E_{\pi}}\circ\phi=\phi\circ\widetilde{\psi}_{*}$. The vector field
$\widetilde{X}^{T}$, interpreted in $E_{\pi}$, can be written as

$\phi_{*}(\widetilde{X}^{T}{}_{Y})=\left.\frac{\mathrm{d}}{\mathrm{d}t}\phi\left((\widetilde{\psi}_{t})_{*}Y\right)\right|_{t=0}=\left.\frac{\mathrm{d}}{\mathrm{d}t}(\psi_{t})^{E_{\pi}}\left(\phi(Y)\right)\right|_{t=0}.$
(5.48)

It follows that $(\psi_{t})^{E_{\pi}}$ is the $1$-parameter family of
automorphisms induced by $\widetilde{X}^{T}$ in $E_{\pi}$. Since
$\widetilde{\psi}_{t}$ preserves $r^{2}=-\widetilde{g}(\xi,\xi)$,
$\widetilde{\varphi}=\operatorname{Im}(\chi)$, $\widetilde{\omega}$, and
$\widetilde{\nabla}$, the induced lift $\psi_{t}^{E_{\pi}}$ preserves the
objects in (5.46), which are then $\widetilde{X}^{T}$-invariant.

By combining (4.30) and (5.37), we deduce that
$(\widetilde{X}^{T})^{\mathrm{Q}}$ is the pushforward of $\widetilde{X}^{T}$
via the quotient map $T\widetilde{M}_{>2k}\to N_{2k}$. Therefore, in order to
prove the last part, it is enough to show that given some PSK automorphism
$(\widetilde{\psi},\psi)$, the maps induced by $\widetilde{\psi}_{*}\colon
T\widetilde{M}_{>2k}\to T\widetilde{M}_{>2k}$ on $E\times\mathbb{R}_{>2k}$ and
$S$ via the diagram (5.14), are $\psi^{E}_{\mathrm{pr}}$ and $\psi^{S}$
respectively. From the fact that $\widetilde{\psi}$ preserves $r^{2}$, we
infer that it restricts to $\psi^{S}$ on $S$ and it descends to
$\psi\times\mathrm{id}_{\mathbb{R}_{>2k}}$ on $M\times\mathbb{R}_{>2k}$.
Notice that in the diagram (5.13), for every object, we have defined an
associated automorphism, and that they all commute with the undashed arrows in
the diagram. It follows that they also commute with the dashed one, implying
that $\widetilde{\psi}_{*}$, descends to $\psi^{E}_{\mathrm{pr}}$ on
$E\times\mathbb{R}_{>2k}$. ∎

######  Remark 5.13.

The integrality of $k$ is crucial in our functorial lifting of general PSK
isomorphisms. For instance, if $k=\frac{1}{2}$, by Remark 4.3, the twist of
$E_{\pi}$ can be identified with a double cover of itself, and in general,
automorphisms cannot be lifted to a double cover in a canonical way. In
[CST21], the integrality of $k$ was not assumed, but this is nevertheless
consistent with the above, as infinitesimal automorphisms can still be lifted
for general $k$, in a canonical way.

######  Remark 5.14.

The action of the automorphism group $\mathrm{Aut}_{\mathrm{PSK}}(M)$ of the
PSK manifold $\widetilde{M}\rightarrow M$ on the space of parallel sections
$s$ of $E$ induces an action on the Heisenberg isometries $\psi_{s}$. This
action is compatible with the conjugation by the functorial lift of PSK
automorphisms. Therefore, the subgroup of isometries generated by the
functorial lifts of PSK automorphisms and the Heisenberg isometries is a
semidirect product

$\mathrm{Aut}_{\mathrm{PSK}}(M)\ltimes\mathrm{Heis}_{\mathrm{rk}(E)+1}\leq\mathrm{Isom}(N_{2k}).$
(5.49)

This was shown for $\mathbf{H}_{\mathbb{C}}^{n}$ in a different manner in
[CRT21, §3]. For the $n=0$ case, it was shown in [CST21, Theorem 4.5 (p. 116)]
that this semidirect product is in fact the full group of isometries. An
advantage of our functorial approach is that it makes the action of the
semidirect product on $N_{2k}$ transparent and explicit in general.

##

## References

  * [ACDM15] D. V. Alekseevsky, V. Cortés, M. Dyckmanns, and T. Mohaupt.  Quaternionic Kähler metrics associated with special Kähler manifolds.  J. Geom. Phys., 92:271–287, 2015.
  * [ACDV03] D. V. Alekseevsky, V. Cortés, C. Devchand, and A. Van Proeyen.  Flows on quaternionic-Kähler and very special real manifolds.  Comm. Math. Phys., 238(3):525–543, 2003.
  * [BCOV94] M. Bershadsky, S. Cecotti, H. Ooguri, and C. Vafa.  Kodaira-Spencer theory of gravity and exact results for quantum string amplitudes.  Comm. Math. Phys., 165(2):311–427, 1994.
  * [BWW20] J. Buczyński, J. A. Wiśniewski, and A. Weber.  Algebraic torus actions on contact manifolds.  arXiv:1802.05002, 2020.
  * [CDF+97] A. Ceresole, R. D’Auria, S. Ferrara, W. Lerche, J. Louis, and T. Regge.  Picard-Fuchs equations, special geometry and target space duality.  In Mirror symmetry, II, volume 1 of AMS/IP Stud. Adv. Math., pages 281–353. Amer. Math. Soc., Providence, RI, 1997.
  * [CDJL18] V. Cortés, M. Dyckmanns, M. Jüngling, and D. Lindemann.  A class of cubic hypersurfaces and quaternionic Kähler manifolds of co-homogeneity one.  arXiv:1701.7882, 2018.
  * [CDS17] V. Cortés, M. Dyckmanns, and S. Suhr.  Completeness of projective special Kähler and quaternionic Kähler manifolds.  In Special metrics and group actions in geometry, volume 23 of Springer INdAM Ser., pages 81–106. Springer, Cham, 2017.
  * [CFG89] S. Cecotti, S. Ferrara, and L. Girardello.  Geometry of type II superstrings and the moduli of superconformal field theories.  Internat. J. Modern Phys. A, 4(10):2475–2529, 1989.
  * [CHM12] V. Cortés, X. Han, and T. Mohaupt.  Completeness in supergravity constructions.  Comm. Math. Phys., 311(1):191–213, 2012.
  * [Cor98] V. Cortés.  On hyper-Kähler manifolds associated to Lagrangian Kähler submanifolds of $T^{*}{\bf C}^{n}$.  Trans. Amer. Math. Soc., 350(8):3193–3205, 1998.
  * [CRT21] V. Cortés, M. Röser, and D. Thung.  Complete quaternionic kähler manifolds with finite volume ends.  arXiv:2105.00727, 2021.
  * [CST21] V. Cortés, A. Saha, and D. Thung.  Symmetries of quaternionic kähler manifolds with ${S}^{1}$-symmetry.  Trans. London Math. Soc., 8(1):95–119, 2021.
  * [DS99] A. Dancer and A. Swann.  Quaternionic Kähler manifolds of cohomogeneity one.  Internat. J. Math., 10(5):541–570, 1999.
  * [FK82] Th. Friedrich and H. Kurke.  Compact four-dimensional self-dual Einstein manifolds with positive scalar curvature.  Math. Nachr., 106:271–299, 1982.
  * [Fre99] D. S. Freed.  Special Kähler manifolds.  Comm. Math. Phys., 203(1):31–52, 1999.
  * [FS90] S. Ferrara and S. Sabharwal.  Quaternionic manifolds for type ${\rm II}$ superstring vacua of Calabi-Yau spaces.  Nuclear Phys. B, 332(2):317–332, 1990.
  * [Hay08] A. Haydys.  HyperKähler and quaternionic Kähler manifolds with $S^{1}$-symmetries.  J. Geom. Phys., 58(3):293–306, 2008.
  * [HHP10] C. Hertling, L. Hoevenaars, and H. Posthuma.  Frobenius manifolds, projective special geometry and Hitchin systems.  J. Reine Angew. Math., 649:117–165, 2010.
  * [Hit81] N. J. Hitchin.  Kählerian twistor spaces.  Proc. London Math. Soc. (3), 43(1):133–150, 1981.
  * [KMS93] I. Kolář, P. W. Michor, and J. Slovák.  Natural operations in differential geometry.  Springer-Verlag, Berlin, 1993.
  * [LS94] C. LeBrun and S. M. Salamon.  Strong rigidity of positive quaternion-Kähler manifolds.  Invent. Math., 118(1):109–132, 1994.
  * [Man21] M. Mantegazza.  Construction of projective special kähler manifolds.  Annali di Matematica Pura ed Applicata (1923 -), apr 2021.
  * [MS15] O. Macia and A. Swann.  Twist geometry of the c-map.  Comm. Math. Phys., 336(3):1329–1357, 2015.
  * [MS20] O. Macia and A. Swann.  The c-map on groups.  Classical Quantum Gravity, 37(1):015015, 17, 2020.
  * [PS91] Y. S. Poon and S. M. Salamon.  Quaternionic Kähler $8$-manifolds with positive scalar curvature.  J. Differential Geom., 33(2):363–378, 1991.
  * [RSV06] D. Robles-Llana, F. Saueressig, and S. Vandoren.  String loop corrected hypermultiplet moduli spaces.  J. High Energy Phys., (3):081, 35, 2006.
  * [Sal89] S. M. Salamon.  Riemannian geometry and holonomy groups, volume 201 of Pitman Research Notes in Mathematics Series.  Longman Scientific & Technical, Harlow; copublished in the United States with John Wiley & Sons, Inc., New York, 1989.
  * [Sim88] C. T. Simpson.  Constructing variations of Hodge structure using Yang-Mills theory and applications to uniformization.  J. Amer. Math. Soc., 1(4):867–918, 1988.
  * [Spi79] M. Spivak.  A comprehensive introduction to differential geometry. Vol. I.  Publish or Perish, Inc., Wilmington, Del., second edition, 1979.
  * [Swa10] A. Swann.  Twisting Hermitian and hypercomplex geometries.  Duke Math. J., 155(2):403–431, 2010.
  * [Tu17] L. W. Tu.  Differential geometry, volume 275 of Graduate Texts in Mathematics.  Springer, Cham, 2017.  Connections, curvature, and characteristic classes.
  * [Voi07] C. Voisin.  Hodge theory and complex algebraic geometry. I, volume 76 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, Cambridge, english edition, 2007.  Translated from the French by Leila Schneps.
  * [YI66] K. Yano and S. Ishihara.  Differential geometry in tangent bundle.  Kōdai Math. Sem. Rep., 18:271–292, 1966.


#  Estimating The Proportion of Signal Variables Under
Arbitrary Covariance Dependence

X. Jessie Jeng
Department of Statistics, North Carolina State University Address for
correspondence: Department of Statistics, North Carolina State University,
2311 Stinson Dr., Raleigh, NC 27695-8203, USA. Email<EMAIL_ADDRESS>
###### Abstract

Estimating the proportion of signals hidden in a large amount of noise
variables is of interest in many scientific inquires. In this paper, we
consider realistic but theoretically challenging settings with arbitrary
covariance dependence between variables. We define mean absolute correlation
(MAC) to measure the overall dependence level and investigate a family of
estimators for their performances in the full range of MAC. We explicit the
joint effect of MAC dependence and signal sparsity on the performances of the
family of estimators and discover that no single estimator in the family is
most powerful under different MAC dependence levels. Informed by the
theoretical insight, we propose a new estimator to better adapt to arbitrary
covariance dependence. The proposed method compares favorably to several
existing methods in extensive finite-sample settings with strong to weak
covariance dependence and real dependence structures from genetic association
studies.

Keywords: Dependence adaptivity; High-dimension data; Lower bound estimator;
Sparse signal

##  1 Introduction

We consider the problem of estimating the proportion of information bearing
signals that are sparsely located in a large amount of noise variables. This
problem is of interest in many scientific inquiries. For example, estimation
of signal proportion is required by multiple testing methods to calculate
local false discovery rate (Efron, 2007), to derive q-value (Storey, 2003),
and to improve power (Storey, 2002; Finner and Gontscharuk, 2009). Moreover,
in many multi-stage studies, estimation of signal proportion can assist
efficient pre-screening and sample size calculation (Cai and Sun, 2017). A
recent line of research, which focuses on retaining a high proportion of
signals through efficient false negative control, also replies on the
estimation of signal proportion as a benchmark for signal inclusion (Jeng et
al., 2016, 2019).

Although estimation of signal proportion is widely requested, methodology
development has met two major challenges. First, signals of different sparsity
levels often call for different estimation methods, while signal sparsity
levels are unknown a priori. Secondly, the large set of variables under
investigation may have complex dependence structures. There exist a number of
rigorously developed methods. Most of them, however, assume independence
between variables (Genovese and Wasserman, 2004; Meinshausen and Rice, 2006;
Jin and Cai, 2007a) and sparsity levels within a certain range (Cai et al.,
2007; Jin, 2008). An extensive review of the existing methods can be found in
Chen (2019). More recent developments extend the study to consider specific
dependence structures. For example, Jeng et al. (2019) studies the problem
assuming block-diagonal covariance structures for the variables; Jeng and Chen
(2019) considers the problem in linear regression and imposes certain
dependence and sparsity conditions to facilitate accurate precision matrix
estimation and bias mitigation. There lacks a method to consistently estimate
signal proportion under arbitrary covariance dependence when signal sparsity
is unknown and possibly falls in a wide range. Such an universal estimator can
have far-reaching impact in real applications.

In this paper, we define mean absolute correlation (MAC) to measure the
overall covariance dependence level and investigate a family of estimators for
their performances in the full extent of MAC. We explicate the joint effect of
MAC dependence and signal sparsity on the performances of the family of
estimators and discover that the most powerful estimator under independence is
no longer the best candidate under arbitrary covariance dependence. Moreover,
it shows that no single estimator in the family is most powerful under
different MAC dependence levels. Informed by the theoretical insight, we
propose a new estimator to better adapt to arbitrary covariance dependence and
signal sparsity.

The new method is compared with several popular methods in extensive
simulation examples including strong to weak covariance dependence between
variables and real dependence structures from genetic association studies. It
shows that although the winner of the several existing methods changes over
different settings, performance of the new method is either comparable to or
better than the performance of that winner in each setting. We apply the new
method to analyze two real datasets. The first dataset is from an expression
quantitative trait loci (eQTL) study with 8637 candidate single-nucleotide
polymorphisms (SNPs), for which the overall dependence in terms of the MAC
level is rather weak. The second dataset is from a classical association study
with microarray data, where 4088 candidate genes possess much stronger overall
dependence with a high MAC level. We compare the estimates of the new method
with those of the existing methods and validate the results through
simulations under the same dependence structures. Our method exhibits better
accuracy for signals of different signal sparsity levels under the real
dependence structures.

The rest of the paper is organized as follows. Section 2 first introduces a
family of lower bound estimators and develops a general result on their
estimation consistency under arbitrary dependence. Then, for specific
estimators in the family, the joint effects of signal sparsity and covariance
dependence are explicated. Consequently, a new and more powerful estimator is
developed under arbitrary covariance dependence. Section 3 compares the
proposed method with existing methods in simulation examples. Section 4
applies the new and existing methods to real genetic association studies.
Section 5 concludes the work with further discussions. Technical proofs are
provided in Appendix.

##  2 Method and Theory

Denote $I_{0}$ and $I_{1}$ as the sets of indices for signal and noise
variables, respectively. We consider the marginal distribution of $p$
variables as

$X_{j}\sim F_{0}\cdot 1\\{j\in I_{0}\\}+F_{1}\cdot 1\\{j\in I_{1}\\},\qquad
j=1,\ldots,p,$

where $F_{0}$ and $F_{1}$ are the null and signal distribution, respectively.
Define the signal proportion

$\pi=|I_{1}|/p.$

We assume that $F_{0}$ is continuous and known a priori. All the other
components are unknown. Our goal is to estimate the signal proportion $\pi$
without the need to specify $F_{1}$ or to identify which variables are
signals.

###  2.1 A family of estimators

Meinshausen and Rice (2006) introduced a family of proportion estimators that
are built upon the empirical process of $p$-values under independence. This
family of estimators have been proved to provide lower bound estimates for the
true proportion $\pi$. Members in the family are indexed by the choice of a
bounding function, and it has been shown that the estimator with bounding
function $\delta(u)=\sqrt{u(1-u)}$ has the best overall performance for
signals of different sparsity levels. This conclusion, unfortunately, does not
hold anymore under arbitrary dependence. When variables are arbitrarily
dependent, the limiting distribution of the empirical process of $p$-values is
generally unknown and may not even have an analytic expression. Moreover, the
dependence effect mingles with signal sparsity to influence the performances
of different estimators in the family. These difficulties substantially
complicate the estimation problem and motivate us to develop new techniques to
study the family of estimators and to come up with a new powerful estimator
under arbitrary covariance dependence.

For presentation simplicity, we perform inverse normal transformation as
$Z_{j}=\Phi^{-1}(F_{0}(X_{j}))$, where $\Phi^{-1}$ is the inverse of the
cumulative distribution function of a standard Normal distribution. Then, we
have

$Z_{j}\sim\Phi\cdot 1\\{j\in I_{0}\\}+G\cdot 1\\{j\in I_{1}\\},\qquad
j=1,\ldots,p,$ (1)

where $G$ denotes the signal distribution after inverse normal transformation,
which remains unknown.

Next, we construct a modified family of estimators to accommodate dependence
among $Z_{j}$. Let

$\bar{W}_{p}(t)=p^{-1}\sum_{j=1}^{p}1\\{|W_{j}|>t\\},$

where $(W_{1},\ldots,W_{p})$ follow the joint null distribution of
$(Z_{1},\ldots,Z_{p})$. Denote $\bar{\Phi}(t)=1-\Phi(t)$. For a given function
$\delta(t)$ that is strictly positive on $(0,\infty)$, define

$V_{p,\delta}=\sup_{t>0}\frac{|\bar{W}_{p}(t)-2\bar{\Phi}(t)|}{\delta\left(t\right)}.$
(2)

Apparently, $V_{p,\delta}$ varies with the choice of $\delta(t)$, and
$\delta(t)$ is called a bounding function. For a given $\delta(t)$ and a
control level $\alpha$, we define the corresponding bounding sequence
$c_{p,\delta}$ as a function of $p,\delta$ and $\alpha$ that satisfies the
following properties:

(a) $pc_{p,\delta}>p_{0}c_{p_{0},\delta}$, where $p_{0}=|I_{0}|$, and

(b) $P(V_{p,\delta}>c_{p,\delta})<\alpha$ for all $p$.
From the above properties, it can be seen that $c_{p,\delta}$ is an upper
bound of $V_{p,\delta}$, while $V_{p,\delta}$ relies on the joint null
distribution. Both $V_{p,\delta}$ and $c_{p,\delta}$ carry the information of
dependence among the variables. Also, the absolute sign in the numerator of
(2) stabilizes $V_{p,\delta}$ as $\bar{W}_{p}(t)-2\bar{\Phi}(t)$ may not be
asymptotically symmetric anymore under arbitrary dependence.

Then, given the observed $Z_{j}$, a family of estimators indexed by $\delta$
are constructed as

$\hat{\pi}_{\delta}=\sup_{t>0}{\frac{\bar{F}_{p}(t)-2\bar{\Phi}(t)-c_{p,\delta}\delta\left(t\right)}{1-2\bar{\Phi}(t)}},$
(3)

where

$\bar{F}_{p}(t)=p^{-1}\sum_{j=1}^{p}1\\{|Z_{j}|>t\\}.$

It can be seen that the $\hat{\pi}_{\delta}$ family of estimators hinge on the
choices of the bounding function $\delta(t)$ and the corresponding bounding
sequence $c_{p,\delta}$. As the bounding sequence $c_{p,\delta}$ implies a
normal range of $\bar{F}_{p}(u)$ if all variables are noise, the information
carried by the observed $\bar{F}_{p}(u)$ that exceeds the normal range
presents evidence for the existence of signals. This version of lower bound
family is for signals of two-sided effects. Minor changes to accommodate one-
sided signal effect is straightforward. Details for the numerical
implementation of $\hat{\pi}_{\delta}$ can be found at the end of Section 2.3.

As shown in Meinshausen and Rice (2006), for signals of different sparsity
levels, the performances of the estimators in the original family can be very
different under independence. More specifically, the signal proportion can be
re-parameterized as $\pi=p^{-\gamma}$, $\gamma\in(0,1)$, with
$\gamma\in(0,1/2)$ representing the relatively dense case and
$\gamma\in[1/2,1)$ representing the more sparse case. Consistency of different
estimators in the original family was proved under independence for
$\gamma\in(0,1/2)$ and $\gamma\in[1/2,1)$ separately. Here, we consider the
estimation problem in the more challenging setting with arbitrary dependence
and the modified family of estimators to accommodate the dependence.

First, we impose a monotonicity constraint on the bounding function
$\delta(t)$, which allows us to present the following result under general
dependence for signals of sparsity levels in the full range of
$\gamma\in(0,1)$.

###### Theorem 2.1

Consider model (1). For a given bounding function $\delta(t)$, if there exist
a bounding sequence $c_{p,\delta}$ satisfying the properties in (a) and (b),
then

$P(\hat{\pi}_{\delta}<\pi)\geq 1-\alpha.$ (4)

On the other hand, for $\pi$ satisfying $0<\pi\ll 1$, if $\delta(t)$ is non-
increasing with respect to $t$, and $G=G_{p}$ such that $G_{p}(\tau)\to 0$ or
$G_{p}(-\tau)\to 1$ for some $\tau$ such that $\tau\gg 1$ and
$\delta(\tau)\ll\pi/c_{p,\delta}$, then

$P(\hat{\pi}_{\delta}>(1-\epsilon)\pi)\to 1$ (5)

for any constant $\epsilon>0$.

The above theorem says that the lower bound property of $\hat{\pi}_{\delta}$
in (4) holds as long as a bounding sequence $c_{p,\delta}$ satisfying (a) and
(b) can be found. On the other hand, the upper bound property of
$\hat{\pi}_{\delta}$ in (5) holds under certain conditions on the signal
distribution $G$, which essentially says that the signal effect, either
positive ($G<\Phi$) or negative ($G>\Phi$), is strong enough. When both lower
bound and upper bound conditions are satisfied for a given $\delta(t)$ and a
degenerating $\alpha$ such that $\alpha=\alpha_{p}\to 0$, $\hat{\pi}_{\delta}$
consistently estimate the true signal proportion, i.e., for any constant
$\epsilon>0$,

$P((1-\epsilon)\pi\leq\hat{\pi}_{\delta}<\pi)\to 1.$

The above results hold for arbitrarily dependent variables and for sparse
signals that are not of a fixed proportion of all the variables. Such sparsity
condition covers the full range of $\gamma\in(0,1)$.

###  2.2 Joint effect of MAC dependence and sparsity

As individual estimators in the $\hat{\pi}_{\delta}$ family hinge on the
choice of the bounding function, we further study their consistency with
specific $\delta(t)$ functions. We focus on $\delta(t)$ of the form
$\delta(t)=[\bar{\Phi}(t)]^{\theta}$, $\theta\in[0,1]$, so that the
corresponding $\hat{\pi}_{\delta}$ are closely related to several existing
estimators whose consistency has been studied under independence in literature
(please see Meinshausen and Rice (2006) and the references therein). Valuable
insight can be obtained by comparing the estimators’ performances under
independence and dependence, which helps us construct a new and more powerful
estimator under arbitrary covariance dependence.

In order to explicate the effect of dependence on $\hat{\pi}_{\delta}$, we
assume that

$(Z_{1},\ldots,Z_{p})\sim N_{p}(\mu,\Sigma),$ (6)

where $\mu$ is a $p$-dimensional sparse vector with $\mu_{j}=A\cdot 1\\{j\in
I_{1}\\}$, $A\neq 0$, and $\Sigma$ is an arbitrary correlation matrix, i.e.,
$\Sigma_{ij}=Corr(Z_{i},Z_{j})$. For presentation simplicity, assume $A>0$. We
define the Mean Absolute Correlation (MAC) to calibrate the covariance
dependence as

$\bar{\rho}_{\Sigma}=\sum_{i=1}^{p}\sum_{j=1}^{p}|\Sigma_{ij}|/p^{2}.$ (7)

A larger value of $\bar{\rho}_{\Sigma}$ indicates stronger overall dependence.

Moreover, we employ a discretization technique from Arias-Castro et al. (2011)
and Jeng and Chen (2019) as follows. Define $\mathbb{T}=[1,\sqrt{5\log
p}]\cap\mathbb{N}$ and the discretized version of $V_{p,\delta}$ as

$V^{*}_{p,\delta}=\max_{t\in\mathbb{T}}\frac{|\bar{W}_{p}(t)-2\bar{\Phi}(t)|}{\delta(t)}.$
(8)

Denote $c^{*}_{p,\delta}$ as the bounding sequence based on
$V^{*}_{p,\delta}$, and define the corresponding proportion estimator as

$\hat{\pi}^{*}_{\delta}=\max_{t\in\mathbb{T}}\frac{\bar{F}_{p}(t)-2\bar{\Phi}(t)-c^{*}_{p,\delta}\delta(t)}{1-2\bar{\Phi}(t)}.$
(9)

Next, we explicates how the MAC level interacts with signal sparsity and
signal intensity to influence the consistency of $\hat{\pi}^{*}_{\delta}$ with
$\delta(t)=[\bar{\Phi}(t)]^{\theta}$ and reveals very different results for
$\theta\in[0,1/2]$ and $\theta\in(1/2,1]$.

###### Theorem 2.2

Consider model (6). Let $\delta(t)=[\bar{\Phi}(t)]^{\theta}$ with
$\theta\in[0,1/2]$. Then, there exists a bounding sequence

$c_{p,\delta}^{*}=O\left(\sqrt{\bar{\rho}_{\Sigma}(\log
p)^{\theta+1/2}}\right),$

that satisfies properties (a) and (b), and the corresponding estimator
$\hat{\pi}^{*}_{\delta}$ satisfies $P(\hat{\pi}^{*}_{\delta}<\pi)\geq
1-\alpha$. Moreover, for $\pi$ satisfying $0<\pi\ll 1$, if $A=A_{p}$ such that
$A_{p}\gg 1$ and

$A_{p}-\bar{\Phi}^{-1}\left({\pi^{1/\theta}\over\bar{\rho}_{\Sigma}^{1/(2\theta)}(\log
p)^{(\theta+1/2)/(2\theta)}}\right)\to\infty,$ (10)

then $P(\hat{\pi}^{*}_{\delta}>(1-\epsilon)\pi)\to 1$ for any constant
$\epsilon>0$.

The above theorem says that for $\delta(t)=[\bar{\Phi}(t)]^{\theta}$ with
$\theta\in[0,1/2]$, we can find a bounding sequence $c_{p,\delta}^{*}$, whose
order increases with respect to $\bar{\rho}_{\Sigma},p$ and $\theta$,
respectively. The second part of the theorem provides a condition on the
signal distribution for the consistency of $\hat{\pi}^{*}_{\delta}$. This
condition is easier to be satisfied for less sparse signals (larger $\pi$) or
less dependent variables (smaller $\bar{\rho}_{\Sigma}$). Note that
$\bar{\Phi}^{-1}(\pi^{1/\theta}/(\bar{\rho}_{\Sigma}^{1/(2\theta)}(\log
p)^{(\theta+1/2)/(2\theta)}))$ is well-defined only for
$\pi<\sqrt{\bar{\rho}_{\Sigma}(\log p)^{\theta+1/2}}$. In the case
$\pi\geq\sqrt{\bar{\rho}_{\Sigma}(\log p)^{\theta+1/2}}$, condition on $A_{p}$
is simply $A_{p}\gg 1$. For the special case with independent variables,
$\bar{\rho}_{\Sigma}=1/p$ and(10) degenerates to

$A_{p}-\bar{\Phi}^{-1}\left(\pi^{1/\theta}p^{1/2\theta}/(\log
p)^{(\theta+1/2)/(2\theta)}\right)\to\infty,$

which agrees with the sufficient and necessary condition for the consistency
of $\hat{\pi}_{\delta}$ under independence for relatively sparse signals. The
comparison can be made by adopting the same parameterization as in Theorem 3
of Meinshausen and Rice (2006) with $\pi=p^{-\gamma}$, $\gamma\in[1/2,1)$,
$\nu=\theta$, and $\kappa=2$.

Results in Theorem 2.2 explicate how the performance of
$\hat{\pi}^{*}_{\delta}$ with
$\delta(t)=[\bar{\Phi}(t)]^{\theta},\theta\in[0,1/2]$ deteriorates as the MAC
dependence gets stronger. These results, however, cannot be extended to
$\hat{\pi}^{*}_{\delta}$ with $\theta\in(1/2,1]$. For the latter, we present
the following theorem.

###### Theorem 2.3

Consider model (6). Let $\delta(t)=[\bar{\Phi}(t)]^{\theta}$ with
$\theta\in(1/2,1]$. Then, there exists a bounding sequence
$c_{p,\delta}^{*}=O(\sqrt{\log p})$, that satisfies the properties (a) and
(b), and the corresponding estimator $\hat{\pi}^{*}_{\delta}$ satisfies
$P(\hat{\pi}^{*}_{\delta}<\pi)\geq 1-\alpha$. Moreover, for $\pi$ satisfying
$0<\pi\ll 1$, if $A=A_{p}$ such that $A_{p}\gg 1$ and

$A_{p}-\bar{\Phi}^{-1}\left({\pi^{1/\theta}\over(\log
p)^{1/(2\theta)}}\right)\to\infty.$ (11)

then $P(\hat{\pi}^{*}_{\delta}>(1-\epsilon)\pi)\to 1$ for any constant
$\epsilon>0$.

This theorem shows that for $\delta(t)=[\bar{\Phi}(t)]^{\theta}$ with
$\theta\in(1/2,1]$, we can find a bounding sequence $c_{p,\delta}^{*}$, whose
order does not involve the MAC dependence level $\bar{\rho}_{\Sigma}$.
Moreover, the intensity condition in (11) does not involve
$\bar{\rho}_{\Sigma}$, which means that $\hat{\pi}^{*}_{\delta}$ with
$\theta\in(1/2,1]$ is consistent under (11), no matter how strong the MAC
dependence is.

Although the above analyses are for the discretized version of
$\hat{\pi}_{\delta}$, they provides important insight on the performance of
the $\hat{\pi}_{\delta}$ family under dependence, which is also supported by
extensive simulation studies in Section 3.

###  2.3 A new estimator for dependent variables

Study in the previous section reveals very different effects of the MAC
dependence on different estimators in the $\hat{\pi}_{\delta}$ family. When
MAC dependence is relatively weak, estimators with $\theta\in[0,1/2]$ may be
more powerful as condition (10) is less restrictive for smaller
$\bar{\rho}_{\Sigma}$. On the other hand, when MAC dependence gets stronger,
estimators with $\theta\in(1/2,1]$ could be more powerful as condition (11) is
free of $\bar{\rho}_{\Sigma}$. Motivated by these findings, we propose to
construct a new estimator of the form

$\hat{\pi}_{adap}=\max\\{\hat{\pi}_{\delta},\delta\in\Delta\\},$

where $\Delta$ is the set of $\delta(t)$ functions that render the most
powerful estimators in different dependence scenarios. Note that such
$\delta(t)$ functions in $\Delta$ all result in conservative estimators with
the lower bound property as stated in Theorem 2.1. Therefore, the new
estimator $\hat{\pi}_{adap}$ naturally inherit the lower bound property.

The power of $\hat{\pi}_{adap}$ depends on the candidate $\delta(t)$ functions
in $\Delta$. We consider $\delta(t)=[\bar{\Phi}(t)]^{\theta}$ with specific
$\theta$ values. Informed by Theorem 2.2, which is for all the estimators with
$\theta\in[0,1/2]$, we see that condition (10) is less stringent with larger
$\theta$. Therefore, we expect the estimator with $\theta=1/2$ to be the most
powerful candidate in this group. On the other hand, based on the results in
Theorem 2.3, we expect the most powerful estimator in the group of estimators
with $\theta\in(1/2,1]$ to have $\theta=1$. The above analysis helps us narrow
down to two candidates that are most powerful in their own $\theta$ groups.
Comparing these two candidates, we see that $\theta=1/2$ can result in a more
(or less) powerful estimator than $\theta=1$ when MAC dependence is relatively
weak (or strong) with $\bar{\rho}_{\Sigma}\leq\pi/\sqrt{\log p}$ (or
$\bar{\rho}_{\Sigma}>\pi/\sqrt{\log p}$). Since none of the candidates
dominates the other under arbitrary covariance dependence, our new estimator
is constructed as

$\hat{\pi}_{adap}=\max\\{\hat{\pi}_{0.5},\hat{\pi}_{1}\\},$ (12)

where $\hat{\pi}_{0.5}$ denotes $\hat{\pi}_{\delta}$ with
$\delta=[\bar{\Phi}(t)]^{1/2}$, and $\hat{\pi}_{1}$ denotes
$\hat{\pi}_{\delta}$ with $\delta=\bar{\Phi}(t)$. $\hat{\pi}_{adap}$ is likely
to be comparable to $\hat{\pi}_{0.5}$ and outperform $\hat{\pi}_{1}$ under
relatively weak dependence, and be comparable to $\hat{\pi}_{1}$ and
outperform $\hat{\pi}_{0.5}$ when covariance dependence is strong.

Numerical Implementation. We conclude this section with additional notes on
the numerical implementation of $\hat{\pi}_{adap}$. Specifically, we simulate
$(w_{1},\ldots,w_{p})$ following the joint null distribution of $Z_{j}$. When
the joint null distribution is unknown in real applications,
$(w_{1},\ldots,w_{p})$ can often be simulated non-parametrically. For example,
when $(Z_{1},\ldots,Z_{p})$ are a set of test statistics for associations
between a set of explanatory variables and a response variable, a common
practice to simulate $(w_{1},\ldots,w_{p})$ is by randomly shuffling only the
sample of the response variable to remove the potential associations. More
details for such permutation approaches can be found in Westfall and Young
(1993). Then, following (2), we simulate $V_{p,0.5}$ and $V_{p,1}$
corresponding to $\theta=1/2$ and $1$ by taking the maximums over
$t=w_{1},\ldots,w_{p}$. Next, we set $\alpha=0.1$ and generate $c_{p,0.5}$ and
$c_{p,1}$ as the $(1-\alpha)$th quantile of 1000 replicates of $V_{p,0.5}$ and
$V_{p,1}$, respectively. The simulated $c_{p,0.5}$ and $c_{p,1}$ are
implemented to calculate $\hat{\pi}_{0.5}$ and $\hat{\pi}_{1}$ as in (3) by
taking the maximums over $t=z_{1},\ldots z_{p}$, where $z_{1},\ldots z_{p}$
are the observed variables after inverse normal transformation. Finally, the
new estimator $\hat{\pi}_{adap}$ is calculated by (12).

##  3 Simulation Study

In the following simulation examples, we consider six dependence structures:
(a)-(d) are commonly observed correlation structures in literature and (e)-(f)
are real correlation structures from genetic association studies. In all the
examples, $\Sigma_{ii}=1,i=1,\ldots,p$.

  1. (a)

Autoregressive. $\Sigma_{ij}=r^{|i-j|}$ and $r=0.9$.

  2. (b)

Equal correlation. $\Sigma_{ij}=0.5$ for $i\neq j$.

  3. (c)

Block correlation. $\Sigma$ has square diagonal blocks. The off-diagonal
elements in the blocks are 0.5, and the elements outside the blocks are zero.

  4. (d)

Sparse correlation. $\Sigma$ has nonzero elements randomly located. The data
generation process is similar to Model 3 in Cai et al. (2013). Let
$\Sigma^{*}=(\sigma_{ij})$, where $\sigma_{ii}=1$, $\sigma_{ij}=0.9*$
Bernoulli$(1,0.1)$ for $i<j$ and $\sigma_{ji}=\sigma_{ij}$. Then
$\Sigma=I^{1/2}(\Sigma^{*}+\delta I)/(1+\delta)I^{1/2}$, where
$\delta=|\lambda_{min}(\Sigma^{*})|+0.05$.

  5. (e)

SNP correlation. $\Sigma$ is the sample correlation matrix of the real SNP
data on Chromosome 21 from 90 individuals in the International HapMap project.

  6. (f)

Gene correlation. $\Sigma$ is the sample correlation matrix of the real gene
expression data from 71 individuals in a riboflavin production study.

We generate test statistics $Z_{1},\ldots,Z_{p}\sim
N((\mu_{1},\ldots,\mu_{p}),\Sigma)$ and set $p=2000$ for cases (a)-(d) above.
Case (e) has $p=8657$, which is the number of SNPs in the dataset, and case
(f) has $p=4088$, which is the number of genes in the dataset. Additional
details of the datasets can be found in Section 4. The block size in case (c)
is set as $400\times 400$. $(\mu_{1},\ldots,\mu_{p})$ is a sparse vector with
randomly located non-zero elements. We consider both relatively sparse signals
with $\pi=0.02$ and more dense signals with $\pi=0.1$.

###  3.1 MAC dependence effect on bounding sequences

We first calculate the MAC levels as defined in (7) for the dependence
structures in (a)-(f) above and report the realized values of $c_{p,0.5}$ and
$c_{p,1}$, which are generated by the procedure described at the end of
Section 2.3. Recall that a larger value of $\bar{\rho}_{\Sigma}$ indicates
stronger overall covariance dependence. It can be seen in Table 1 that
$\bar{\rho}_{\Sigma}$ is fairly small for cases (a) and (d), moderately small
for cases (c) and (e), and fairly large for cases (b) and (f). Moreover,
$c_{p,0.5}$ seems to vary positively with $\bar{\rho}_{\Sigma}$, whereas
$c_{p,1}$ does not show such tendency. The numerical results are demonstrated
more clearly in Figure 1, which seem to agree with the theoretical results in
Theorem 2.2 and 2.3 about the MAC dependence effect on $c_{p,0.5}$ and
$c_{p,1}$.

Table 1: MAC levels and realized values of bounding sequences. | Autocorr | Equal corr | Block corr | Sparse corr | SNP corr | Gene corr
---|---|---|---|---|---|---
$\bar{\rho}_{\Sigma}$ | 0.0095 | 0.5003 | 0.1003 | 0.0042 | 0.0869 | 0.3353
$c_{p,0.5}$ | 0.178 | 0.87 | 0.397 | 0.099 | 0.222 | 0.706
$c_{p,1}$ | 8.46 | 4.39 | 5.58 | 6.79 | 13.6 | 6.42
Figure 1:  Trends of MAC levels and bounding sequences.

###  3.2 Comparison with existing methods

We compare the new estimator $\hat{\pi}_{adap}$ with $\hat{\pi}_{0.5}$ and
$\hat{\pi}_{1}$ from the estimator family as well as two other popular
methods. Note that $\hat{\pi}_{0.5}$ is closely related to the optimal
estimator in the original family of Meinshausen and Rice (2006) under
independence, and $\hat{\pi}_{1}$ has been discovered to have good power for
relatively sparse signals under independence. The two other existing methods
are $\hat{\pi}_{GW}$ developed in Genovese and Wasserman (2004) and
$\hat{\pi}_{JC}$ developed in Jin and Cai (2007b). These two existing methods
have been studied for relatively dense signals under independence.

Besides various dependence structures in (a)-(f), we consider sparse and
relatively dense signals with $\pi=0.02$ and $0.1$, respectively, and varying
signal intensity with non-zero $\mu_{j}=3,4,5,6$. The comparisons are
organized into two sets of examples.

The first set of examples demonstrate the performances of $\hat{\pi}_{adap}$,
$\hat{\pi}_{0.5}$, and $\hat{\pi}_{1}$ as they all possess the lower bound
property. Recall the theoretical results in Section 2.2 that $\hat{\pi}_{0.5}$
may outperform $\hat{\pi}_{1}$ when dependence is weak enough or signals are
less sparse, and $\hat{\pi}_{1}$ may perform better in the other scenarios. We
observe such tendencies in Figure 2-7. Specifically, the autocorrelation case
(Figure 2) has small $\bar{\rho}_{\Sigma}=0.0095$. It shows that
$\hat{\pi}_{0.5}$ has comparable results as those of $\hat{\pi}_{1}$ for small
$\pi=0.02$, and outperforms $\hat{\pi}_{1}$ for larger $\pi=0.1$. The equal
correlation case (Figure 3) has the largest $\bar{\rho}_{\Sigma}=0.5$. It
shows that $\hat{\pi}_{1}$ outperforms $\hat{\pi}_{0.5}$ for both $\pi=0.02$
and $0.1$. The block diagonal case (Figure 4) has moderate
$\bar{\rho}_{\Sigma}=0.1$. It shows that $\hat{\pi}_{1}$ outperforms
$\hat{\pi}_{0.5}$ for small $\pi$, and is comparable to $\hat{\pi}_{0.5}$ for
larger $\pi$. The sparse correlation case (Figure 5) has the smallest
$\bar{\rho}_{\Sigma}=0.0042$, we see that $\hat{\pi}_{0.5}$ outperforms
$\hat{\pi}_{1}$ for both $\pi=0.02$ and $0.1$. The SNP correlation case
(Figure 6) has $\bar{\rho}_{\Sigma}=0.0869$, which is moderately small. We see
that $\hat{\pi}_{1}$ is slightly better for small $\pi$, and $\hat{\pi}_{0.5}$
is better for larger $\pi$. The gene correlation case (Figure 7) has
$\bar{\rho}_{\Sigma}=0.3353$, which is fairly large. It shows that
$\hat{\pi}_{1}$ outperforms $\hat{\pi}_{0.5}$ for both small and larger $\pi$.
In all these examples, the new estimator $\hat{\pi}_{adap}$ always matches the
winner of $\hat{\pi}_{0.5}$ and $\hat{\pi}_{1}$ and exhibits better adaptivity
to different dependence structures and signal sparsity levels.

Figure 2: Comparison under autocorrelation. “est_0.05”, “est_1”, and
“est_adapt” represent $\hat{\pi}_{0.5}$, $\hat{\pi}_{1}$, and
$\hat{\pi}_{adap}$, respectively. The top row has $\pi=0.02$, and the bottom
row has $\pi=0.1$. The true $\pi$ values are highlighted by the red horizontal
lines.

Figure 3: Comparison under equal correlation. Notations and symbols are the
same as in Figure 2.

Figure 4: Comparison under block dependence. Notations and symbols are the
same as in Figure 2.

Figure 5: Comparison under sparse dependence. Notations and symbols are the
same as in Figure 2.

Figure 6: Comparison under SNP dependence. Notations and symbols are the same
as in Figure 2.

Figure 7: Comparison under gene dependence. Notations and symbols are the same
as in Figure 2.

The second set of examples compares the performance of $\hat{\pi}_{adap}$ with
those of $\hat{\pi}_{GW}$ from Genovese and Wasserman (2004) and
$\hat{\pi}_{JC}$ from Jin and Cai (2007b). Table 2 and 3 show that for the
cases (a)-(f), $\hat{\pi}_{GW}$ and $\hat{\pi}_{JC}$ tend to over-estimate the
true $\pi$ when dependence is strong (such as in (b) Equal correlation and (f)
Gene correlation) or signals are sparse (such as for $\pi=0.02$). On the other
hand, the performance of $\hat{\pi}_{adap}$ seems to be generally more
accurate and stable over different dependence structures and $\pi$ values.

Table 2: Mean values and standard deviations (in brackets) of $\hat{\pi}_{adap}$, $\hat{\pi}_{GW}$, and $\hat{\pi}_{JC}$ when signals are relatively spare with $\pi=0.02$. Dependence | Method | $\mu=3$ | $\mu=4$ | $\mu=5$ | $\mu=6$
---|---|---|---|---|---
Autocorr | $\hat{\pi}_{adap}$ | 0.020(0.003) | 0.020(0.003) | 0.020(0.002) | 0.021(0.002)
| $\hat{\pi}_{GW}$ | 0.014(0.024) | 0.015(0.024) | 0.014(0.024) | 0.014(0.024)
| $\hat{\pi}_{JC}$ | 0.050(0.045) | 0.050(0.045) | 0.050(0.045) | 0.051(0.044)
Equal corr | $\hat{\pi}_{adap}$ | 0.034(0.051) | 0.036(0.050) | 0.038(0.049) | 0.038(0.049)
| $\hat{\pi}_{GW}$ | 0.194(0.274) | 0.194(0.273) | 0.194(0.273) | 0.194(0.274)
| $\hat{\pi}_{JC}$ | 0.295(0.402) | 0.295(0.402) | 0.295(0.401) | 0.29690.402
Block corr | $\hat{\pi}_{adap}$ | 0.012(0.019) | 0.018(0.017) | 0.023(0.016) | 0.024(0.016)
| $\hat{\pi}_{GW}$ | 0.060(0.091) | 0.060(0.091) | 0.059(0.091) | 0.060(0.091)
| $\hat{\pi}_{JC}$ | 0.115(0.143) | 0.115(0.142) | 0.116(0.141) | 0.118(0.140)
Sparse corr | $\hat{\pi}_{adap}$ | 0.009(0.003) | 0.015(0.022) | 0.018(0.001) | 0.020(0.001)
| $\hat{\pi}_{GW}$ | 0.003(0.006) | 0.003(0.007) | 0.003(0.006) | 0.003(0.007)
| $\hat{\pi}_{JC}$ | 0.030(0.021) | 0.031(0.020) | 0.032(0.020) | 0.033(0.020)
SNP corr | $\hat{\pi}_{adap}$ | 0.034 (0.051) | 0.036(0.050) | 0.038(0.049) | 0.038(0.049)
| $\hat{\pi}_{GW}$ | 0.194(0.274) | 0.194(0.273) | 0.194(0.273) | 0.194(0.274)
| $\hat{\pi}_{JC}$ | 0.295(0.402) | 0.295(0.402) | 0.295(0.401) | 0.296(0.402)
Gene corr | $\hat{\pi}_{adap}$ | 0.019(0.038) | 0.026(0.037) | 0.030(0.035) | 0.032(0.035)
| $\hat{\pi}_{GW}$ | 0.151(0.187) | 0.151(0.188) | 0.151(0.188) | 0.151(0.187)
| $\hat{\pi}_{JC}$ | 0.257(0.308) | 0.256(0.307) | 0.257(0.307) | 0.258(0.306)
Table 3: Mean values and standard deviations (in brackets) of $\hat{\pi}_{adap}$, $\hat{\pi}_{GW}$, and $\hat{\pi}_{JC}$ when signals are relatively dense with $\pi=0.1$. Dependence | Method | $\mu=3$ | $\mu=4$ | $\mu=5$ | $\mu=6$
---|---|---|---|---|---
Autocorr | $\hat{\pi}_{adap}$ | 0.063(0.011) | 0.083(0.007) | 0.094(0.004) | 0.100(0.002)
| $\hat{\pi}_{GW}$ | 0.078(0.029) | 0.086(0.026) | 0.088(0.023) | 0.088(0.024)
| $\hat{\pi}_{JC}$ | 0.119(0.046) | 0.124(0.043) | 0.127(0.043) | 0.125(0.041)
Equal corr | $\hat{\pi}_{adap}$ | 0.067(0.053) | 0.098(0.042) | 0.110(0.038) | 0.113(0.036)
| $\hat{\pi}_{GW}$ | 0.233(0.263) | 0.246(0.256) | 0.250(0.254) | 0.250(0.253)
| $\hat{\pi}_{JC}$ | 0.333(0.394) | 0.336(0.377) | 0.344(0.374) | 0.349(0.373)
Block diag | $\hat{\pi}_{adap}$ | 0.052(0.026) | 0.082(0.018) | 0.097(0.013) | 0.102(0.012)
| $\hat{\pi}_{GW}$ | 0.114(0.094) | 0.126(0.089) | 0.128(0.085) | 0.129(0.086)
| $\hat{\pi}_{JC}$ | 0.173(0.144) | 0.179(0.135) | 0.183(0.131) | 0.184(0.131)
Sparse corr | $\hat{\pi}_{adap}$ | 0.073(0.006) | 0.088(0.004) | 0.096(0.001) | 0.099(0.001)
| $\hat{\pi}_{GW}$ | 0.069(0.012) | 0.075(0.010) | 0.077(0.009) | 0.077(0.009)
| $\hat{\pi}_{JC}$ | 0.106(0.019) | 0.108(0.019) | 0.110(0.019) | 0.110(0.018)
SNP corr | $\hat{\pi}_{adap}$ | 0.057(0.013) | 0.080(0.007) | 0.093(0.004) | 0.098(0.002)
| $\hat{\pi}_{GW}$ | 0.100(0.043) | 0.107(0.038) | 0.110(0.036) | 0.111(0.035)
| $\hat{\pi}_{JC}$ | 0.129(0.070) | 0.133(0.066) | 0.137(0.064) | 0.138(0.063)
Gene corr | $\hat{\pi}_{adap}$ | 0.055(0.038) | 0.087(0.033) | 0.102(0.028) | 0.107(0.026)
| $\hat{\pi}_{GW}$ | 0.202(0.181) | 0.213(0.177) | 0.216(0.173) | 0.217(0.172)
| $\hat{\pi}_{JC}$ | 0.301(0.302) | 0.304(0.292) | 0.312(0.287) | 0.314(0.284)

##  4 Real Application

We apply the proposed method to two real datasets. The first dataset is from
an eQTL study with the goal to identify SNPs that potentially govern the
expression of gene CCT8 on chromosome 21. This gene has been found to be
relevant to Down Syndrome (Bradic et al., 2011; Fan et al., 2012). We obtain
the SNP data of unaffected subjects from the International HapMap project
(http://zzz.bwh.harvard.edu/plink/res.shtml#hapmap) and the gene expression
data from ftp://ftp.sanger.ac.uk/pub/genevar/. Our data includes 90 samples
from Asian population ($45$ Japanese in Tokyo, Japan (JPT), and $45$ Han
Chinese in Beijing). We consider SNPs without missing values, which results in
8657 candidate SNPs.

Test statistics for the associations between each SNP and the expression level
of CCT8 are derived by marginal linear regression as in Bradic et al. (2011)
and Fan et al. (2012). Histogram of the test statistics is presented in Figure
8, where the long and thin right tail indicates possibly a small proportion of
signals with positive signal effects. The correlation matrix of the test
statistics, which is the same as the correlation matrix of the SNPs, has the
MAC level of $\bar{\rho}_{\Sigma}=0.087$. The heatmap of the correlation
matrix of the first 50 SNPs is illustrated in Figure 9.

Figure 8: Histogram of test statistics for eQTL analysis Figure 9: Heatmap of
the absolute value of correlations for 50 SNPs in Hapmap data.

We apply the proposed estimator $\hat{\pi}_{adap}$ and the existing estimators
$\hat{\pi}_{GW}$ and $\hat{\pi}_{JC}$ to the dataset. The results are
$\hat{\pi}_{adap}=0.0016$, $\hat{\pi}_{GW}=0.0068$, and
$\hat{\pi}_{JC}=0.0186$. Consequently, the estimated numbers of relevant SNPs
are $14$, $60$, and $161$ by the three methods, respectively. The dependence
structure of the real data has been adopted in simulation study to help with
justifying the real data results. The MAC level and the realized values of the
bounding sequences $c_{p,0.5}$ and $c_{p,1}$ have been demonstrated in Table
1, case (e). Simulation results that are most relevant are presented in the
top row of Figure 6 and the sixth row of Table 2. The real application results
here seem to be consistent with the findings in simulation, where
$\hat{\pi}_{adap}$ is much smaller than $\hat{\pi}_{GW}$ and $\hat{\pi}_{JC}$
and is closer to the true $\pi$.

The second real application example has microarray data from a study on
riboflavin (vitamin B2) production in bacillus subtilis. This dataset is
available at https://www.annualreviews.org/doi/suppl/10.1146/annurev-
statistics-022513-115545 and has been studied in Buhlmann et al. (2014). The
dataset includes the expression levels of 4088 genes and the logarithm of
riboflavin production rate of $71$ individuals. Marginal regression
coefficients are used as test statistics for associations between genes and
riboflavin production. The histogram of the test statistics is presented in
Figure 10, which suggests a larger signal proportion than that in the first
real data example. Figure 11 shows the heatmap of the correlation matrix of
the first 50 genes, which indicates a more complicated dependence structure.
The MAC level of the genes is $\bar{\rho}_{\Sigma}=0.335$, which is fairly
large. The realized value of the bounding sequences $c_{p,0.5}$ and $c_{p,1}$
have been demonstrated in Table 1, case (f).

Figure 10: Histogram of test statistics for gene expression association Figure
11: Heatmap of the absolute value of correlations for 50 gene expressions in
riboflavin production study.

In this example, we have $\hat{\pi}_{adap}=0.064$, $\hat{\pi}_{GW}=0.258$, and
$\hat{\pi}_{JC}=0.374$, which corresponds to $261$, $1053$, and $1530$
relevant genes. This real dependence structure has been adopted in simulation
study. The most relevant simulation results are presented in the bottom row of
Figure 7 and the last row of Table 3, which help us justify the real
application results. We can see that the results here seem to be consistent
with the findings in simulation, where $\hat{\pi}_{adap}$ is much smaller than
$\hat{\pi}_{GW}$ and $\hat{\pi}_{JC}$, and is closer to the true $\pi$.

##  5 Conclusion and Discussion

Estimating the proportion of sparse signals is notoriously difficult when
dealing with large-scale data with complex dependence structures. In this
paper, we define the MAC level to measure arbitrary covariance dependence and
explicate the joint effect of MAC dependence and signal sparsity on a family
of estimators. We find that no single estimator in the family is most powerful
under different MAC dependence levels. We identify candidate estimators that
are most powerful in different dependence scenarios and develop a new
estimator that better adapts to arbitrary covariance dependence. The new
estimator inherits the lower bound property of the family and provides a
conservative estimate under very general conditions. This property is valuable
in real applications as it requires no conditions on the unknown signals.
Moreover, the new estimator is more powerful than representative members in
the estimator family and compares favorably to other popular methods in
extensive numerical examples including weak to strong covariance dependence
and real dependence structures from genetic associations studies. By pushing
the frontier of high-dimensional sparse inference to better accommodate the
complex data structures in real applications, we expect the impact of the
proposed research to be far-reaching.

## Appendix

This section presents the proofs of Theorem 2.1, Theorem 2.2, and Theorem 2.3.
The symbol $C$ denotes a genetic, finite constant whose values can be
different at different occurrences.

###  5.1 Proof of Theorem 2.1

We first show $P(\hat{\pi}_{\delta}<\pi)\geq 1-\alpha$. Let $Z_{j}^{0}=Z_{j}$
for $j\in I_{0}$ and $Z_{j}^{1}=Z_{j}$ for $j\in I_{1}$. Denote
$p_{0}=|I_{0}|$ and $s=|I_{1}|$. Then

$\displaystyle\bar{F}_{p}(t)$ $\displaystyle=$ $\displaystyle p^{-1}\sum_{j\in
I_{0}}1_{\\{|Z_{j}^{0}|>t\\}}+p^{-1}\sum_{j\in
I_{1}}1_{\\{|Z_{j}^{1}|>t\\}}\leq p^{-1}\sum_{j\in
I_{0}}1_{\\{|Z_{j}^{0}|>t\\}}+p^{-1}s$ $\displaystyle=$
$\displaystyle(1-\pi)p_{0}^{-1}\sum_{j\in I_{0}}1_{\\{|Z_{j}^{0}|>t\\}}+\pi.$

Consequently,

$\displaystyle P(\hat{\pi}_{\delta}>\pi)$ $\displaystyle\leq$ $\displaystyle
P\left(\sup_{t>0}\left\\{(1-\pi)\left(p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|Z_{j}^{0}|>t\\}}-2\bar{\Phi}(t)\right)-c_{p,\delta}\delta(t)\right\\}>0\right)$
$\displaystyle\leq$ $\displaystyle
P\left(\sup_{t>0}\left\\{p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|Z_{j}^{0}|>t\\}}-2\bar{\Phi}(t)-c_{p_{0},\delta}\delta(t)\right\\}>0\right)$
$\displaystyle\leq$ $\displaystyle
P\left(V_{p_{0}}>c_{p_{0},\delta}\right)=\alpha,$

where the second and last inequalities are by properties (a) and (b) of
$c_{p,\delta}$, respectively. The claim $P(\hat{\pi}_{\delta}<\pi)\geq
1-\alpha$ follows.

Next, we show $P(\hat{\pi}_{\delta}>(1-\epsilon)\pi)\to 1$. Because
$\hat{\pi}_{\delta}>\bar{F}_{p}(t)-2\bar{\Phi}(t)-c_{p,\delta}\delta\left(t\right)$
for any $t>0$ and

$\bar{F}_{p}(t)={\frac{1-\pi}{p_{0}}}\sum_{j\in
I_{0}}1_{\\{|Z_{j}^{0}|>t\\}}+{\frac{\pi}{s}}\sum_{j\in
I_{1}}1_{\\{|Z_{j}^{1}|>t\\}},$

then

$\displaystyle\frac{\hat{\pi}_{\delta}}{\pi}-1$ $\displaystyle>$
$\displaystyle-\frac{1}{\pi}c_{p,\delta}\delta(t)+\frac{1-\pi}{\pi}\left(p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|Z_{j}^{0}|>t\\}}-2\bar{\Phi}(t)\right)$ (13) $\displaystyle+$
$\displaystyle\left({s}^{-1}\sum\nolimits_{j\in
I_{1}}1_{\\{|Z_{j}^{1}|>t\\}}-1\right)-2\bar{\Phi}(t)$

for any $t>0$. Now, set $t$ in (13) at $\tau$ such that $\tau\gg 1$ and
$\delta(\tau)\ll\pi/c_{p,\delta}$. We will show that each term on the right
hand side of (13) at $t=\tau$ is of $o_{p}(1)$.

First, by the condition $\delta(\tau)\ll\pi/c_{p,\delta}$, we have the first
term $A_{1}=-c_{p,\delta}\delta(\tau)/\pi=o(1)$.

Consider the second term
$A_{2}=\pi^{-1}(1-\pi)\left(p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|Z_{j}^{0}|>\tau\\}}-2\bar{\Phi}(\tau)\right)$ in (13). The
following lemma is proved in Section 5.2.

###### Lemma 5.1

For any $\tau$ such that $\tau\gg 1$ and $\delta(\tau)\ll\pi/c_{p,\delta}$, we
have

$\pi^{-1}\left(p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|Z_{j}^{0}|>\tau\\}}-2\bar{\Phi}(\tau)\right)=o_{p}(1).$

Lemma 5.1 and the condition $\pi=o(1)$ are enough to show $A_{2}=o_{p}(1)$.

For the third term $A_{3}={s}^{-1}\sum\nolimits_{j\in
I_{1}}1_{\\{|Z_{j}^{1}|>\tau\\}}-1$ in (13), we have

$\displaystyle P(|A_{3}|>a)$ $\displaystyle=$ $\displaystyle
P(1-{s}^{-1}\sum\nolimits_{j\in I_{1}}1_{\\{|Z_{j}^{1}|>\tau\\}}>a)$
$\displaystyle\leq$ $\displaystyle a^{-1}(1-P(|Z_{j}^{1}|>\tau))$
$\displaystyle=$ $\displaystyle a^{-1}\left(G(\tau)-G(-\tau)\right)=o(1)$

for any fixed $a>0$, where the third step is by $Z_{j}^{1}\sim G$ for $j\in
I_{1}$, and the last the step is by the condition $G(\tau)\to 0$ or
$G(-\tau)\to 1$.

Last but not least, the forth term in (13): $A_{4}=-2\bar{\Phi}(\tau)=o(1)$
given $\tau\gg 1$.

Summarizing the above gives the desired result
$P(\hat{\pi}_{\delta}/\pi>1-\epsilon)\to 1$.

###  5.2 Proof of Lemma 5.1

Recall the definitions of $V_{p,\delta}$ in (2) and

$P\left(\sup_{t>0}{\frac{|p^{-1}\sum_{j=1}^{p}1_{\\{|W_{j}|>t\\}}-2\bar{\Phi}(t)|}{\delta\left(t\right)}}>c_{p,\delta}\right)<\alpha.$

where $W_{1},\ldots,W_{p}\sim N_{p}(0,\Sigma)$. Then, for $t=\tau$,

$P\left({\frac{|p^{-1}\sum_{j=1}^{p}1_{\\{|W_{j}|>\tau\\}}-2\bar{\Phi}(\tau)|}{\pi}}>\frac{c_{p,\delta}\delta(\tau)}{\pi}\right)<\alpha.$

Given $c_{p,\delta}\delta(\tau)/\pi=o(1)$, we have

${\frac{|p^{-1}\sum_{j=1}^{p}1_{\\{|W_{j}|>\tau\\}}-2\bar{\Phi}(\tau)|}{\pi}}=o_{p}(1).$
(14)

Decompose the left hand side above as

$\displaystyle\pi^{-1}\left|p^{-1}\sum\nolimits_{j=1}^{p}1_{\\{|W_{j}|>\tau\\}}-2\bar{\Phi}(\tau)\right|$
$\displaystyle\geq$ $\displaystyle\pi^{-1}\left|p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|W_{j}|>\tau\\}}-2\bar{\Phi}(\tau)\right|$ $\displaystyle-$
$\displaystyle\pi^{-1}\left|p^{-1}\sum_{j=1}^{p}1_{\\{|W_{j}|>\tau\\}}-p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|W_{j}|>\tau\\}}\right|$

For the second term on the right hand side,

$\displaystyle\pi^{-1}\left|p^{-1}\sum\nolimits_{j=1}^{p}1_{\\{|W_{j}|>\tau\\}}-p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|W_{j}|>\tau\\}}\right|$ $\displaystyle=$
$\displaystyle\pi^{-1}\left|p^{-1}\sum\nolimits_{j\in
I_{1}}1_{\\{|W_{j}|>\tau\\}}-\pi p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|W_{j}|>\tau\\}}\right|$ (15) $\displaystyle\leq$ $\displaystyle
s^{-1}\sum\nolimits_{j\in
I_{1}}1_{\\{|W_{j}|>\tau\\}}+p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|W_{j}|>\tau\\}}.$

Since $\tau\gg 1$, both $s^{-1}\sum\nolimits_{j\in
I_{1}}1_{\\{|W_{j}|>\tau\\}}=o_{p}(1)$ and $p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|W_{j}|>\tau\\}}=o_{p}(1)$ by Markov’s inequality. Combining this
with (14) and (15) gives

$\pi^{-1}\left|p_{0}^{-1}\sum\nolimits_{j\in
I_{0}}1_{\\{|W_{j}|>\tau\\}}-2\bar{\Phi}(\tau)\right|=o_{p}(1).$

Now, because the joint distribution of $Z_{j}^{0},j\in I_{0}$, is the same as
the joint distribution of $W_{j},j\in I_{0}$, claim in Lemma 5.1 follows.

###  5.3 Proof of Theorem 2.2

First, we show that given
$\delta(t)=[\bar{\Phi}(t)]^{\theta},\theta\in[0,1/2]$,
$c^{*}_{p,\delta}=C\sqrt{\bar{\rho}_{\Sigma}(\log p)^{\theta+1/2}}$, with a
large enough constant $C$, satisfies properties (a)
$pc^{*}_{p,\delta}>p_{0}c^{*}_{p_{0},\delta}$ and (b)
$P(V^{*}_{p,\delta}>c^{*}_{p,\delta})<\alpha$ for all $p$.

Consider property (a). Define $\Sigma_{0}$ as the covariance matrix of
$W_{j},j\in I_{0}$ and

$\bar{\rho}_{\Sigma_{0}}=\sum_{i\in I_{0}}\sum_{j\in
I_{0}}|\Sigma_{ij}|/p^{2}_{0}.$

It can be shown that

$\bar{\rho}_{\Sigma}>{1\over p^{2}}\sum_{i\in I_{0}}\sum_{j\in
I_{0}}|\Sigma_{ij}|={(1-\pi)^{2}\over p_{0}^{2}}\sum_{i\in I_{0}}\sum_{j\in
I_{0}}|\Sigma_{ij}|=(1-\pi)^{2}\bar{\rho}_{\Sigma_{0}}.$

Then it follows that

$c^{*}_{p,\delta}>C(1-\pi)\sqrt{\bar{\rho}_{\Sigma_{0}}(\log
p)^{\theta+1/2}}>(1-\pi)c^{*}_{p_{0},\delta}=(p_{0}/p)c^{*}_{p_{0},\delta},$

and property (a) is verified.

Next consider property (b). By Chebyshev’s inequality and direct calculation,

$\displaystyle P(V_{p,\delta}^{*}>c_{p,\delta}^{*})$ $\displaystyle\leq$
$\displaystyle\left(c_{p,\delta}^{*}\right)^{-2}\mathsf{Var}(V_{p,\delta}^{*})\leq\left(c_{p,\delta}^{*}\right)^{-2}E([V_{p,\delta}^{*}]^{2})$
$\displaystyle=$
$\displaystyle\left(c_{p,\delta}^{*}\right)^{-2}E\left[\max_{t\in\mathbb{T}}\left(\frac{|\bar{W}_{p}(t)-2\bar{\Phi}(t)|}{[\bar{\Phi}(t)]^{\theta}}\right)^{2}\right]$

Let $A(t)=[\bar{\Phi}(t)]^{-2\theta}(\bar{W}_{p}(t)-2\bar{\Phi}(t))^{2}$. It
can be shown that

$\displaystyle E\left(\max_{t\in\mathbb{T}}A(t)\right)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}P(\max_{t\in\mathbb{T}}A(t)>c)dc\leq\int_{0}^{\infty}\sum_{t\in\mathbb{T}}P(A(t)>c)dc$
$\displaystyle=$ $\displaystyle\sum_{t\in\mathbb{T}}E[A(t)]\leq C\sqrt{\log
p}\cdot\max_{t\in\mathbb{T}}E[A(t)]$ $\displaystyle=$ $\displaystyle
C\sqrt{\log
p}\cdot\max_{t\in\mathbb{T}}\left\\{[\bar{\Phi}(t)]^{-2\theta}\mathsf{Var}\left(\bar{W}_{p}(t)\right)\right\\}$

The following lemma provides the order of $\mathsf{Var}(\bar{W}_{p}(t))$.

###### Lemma 5.2

For $W_{1},\ldots,W_{p}\sim N_{p}(0,\Sigma)$ and $\bar{\rho}_{\Sigma}$ in (7),

$\mathsf{Var}\left(\bar{W}_{p}(t)\right)=O\left(\bar{\rho}_{\Sigma}\cdot
e^{-t^{2}/2}\right).$ (16)

Therefore,

$\displaystyle[\bar{\Phi}(t)]^{-2\theta}\cdot\mathsf{Var}\left(\bar{W}_{p}(t)\right)$
$\displaystyle\leq$ $\displaystyle
C[\bar{\Phi}(t)]^{-2\theta}\cdot\bar{\rho}_{\Sigma}\cdot e^{-t^{2}/2}$
$\displaystyle\leq$ $\displaystyle C\left({t\over
e^{-t^{2}/2}}\right)^{2\theta}\cdot\bar{\rho}_{\Sigma}\cdot e^{-t^{2}/2}\leq
C(\log p)^{\theta}\cdot\bar{\rho}_{\Sigma}\cdot e^{(\theta-1/2)t^{2}}\leq
C\bar{\rho}_{\Sigma}\cdot(\log p)^{\theta}$

where the first step above is by Lemma 5.2, the second step is by Mill’s
ratio, the third step is by $t\in\mathbb{T}$, and the last step is by
$\theta\in[0,1/2]$. Combining the above, we have

$P(V_{p,\delta}^{*}>c_{p,\delta}^{*})\leq
C\left(c_{p,\delta}^{*}\right)^{-2}\cdot\bar{\rho}_{\Sigma}\cdot(\log
p)^{\theta+1/2},$

and $c^{*}_{p,\delta}=O(\sqrt{\bar{\rho}_{\Sigma}(\log p)^{\theta+1/2}})$
follows.

Next, we demonstrate the upper bound property of $\hat{\pi}_{\delta}^{*}$.
Denote

$B_{p}=\bar{\Phi}^{-1}\left({\pi^{1/\theta}\over\bar{\rho}_{\Sigma}^{1/(2\theta)}(\log
p)^{(\theta+1/2)/(2\theta)}}\right).$

Let $\tau=(A_{p}+B_{p})/2$. Then, $A_{p}\gg 1$ and condition (10) imply that
$\tau\gg 1$ and $\tau-B_{p}\to\infty$, which further imply

$[\bar{\Phi}(\tau)]^{\theta}\ll\pi/c^{*}_{p,\delta}.$

On the other hand,

$G_{p}(\tau)=\Phi(\tau-A_{p})=\Phi(-(A_{p}-B_{p})/2)=o(1),$

where the last step is by condition (10). The rest is straightforward by
applying Theorem 2.1.

###  5.4 Proof of Lemma 5.2

$\mathsf{Var}\left(\bar{W}_{p}(t)\right)=p^{-2}\sum_{j=1}^{p}Var(1_{\\{|W_{j}|>t\\}})+p^{-2}\sum_{i\neq
j}Cov(1_{\\{|W_{i}|>t\\}},1_{\\{|W_{j}|>t\\}}).$

By Mill’s ratio,

$\displaystyle p^{-2}\sum_{j=1}^{p}Var(1_{\\{|W_{j}|>t\\}})\leq
p^{-1}2\bar{\Phi}(t)(1-2\bar{\Phi}(t))\leq Cp^{-1}e^{-t^{2}/2}.$

For $p^{-2}\sum_{i\neq j}Cov(1_{\\{|W_{i}|>t\\}},1_{\\{|W_{j}|>t\\}})$, we
have

$Cov(1_{\\{|W_{i}|>t\\}},1_{\\{|W_{j}|>t\\}})=4\int_{-\infty}^{t}\int_{-\infty}^{t}f(x,y)dxdy-4\int_{-\infty}^{t}\phi(x)dx\int_{-\infty}^{t}\phi(y)dy\leq
C|\Sigma_{ij}|e^{-t^{2}/2},$

where the last step follows from Corollary 2.1 in Li and Shao (2002).
Combining the above with the definition of $\bar{\rho}_{\Sigma}$ results in
(16).

###  5.5 Proof of Theorem 2.3

First, it is easy to see that $c^{*}_{p,\delta}=C\sqrt{\log p}$ satisfies
property (a) $pc_{p,\delta}>p_{0}c_{p_{0},\delta}$.

For property (b), by Markov’s inequality,

$P(V_{p,\delta}^{*}>c_{p,\delta}^{*})\leq\left(c_{p,\delta}^{*}\right)^{-1}\mathsf{E}\left(\max_{t\in\mathbb{T}}\frac{|\bar{W}_{p}(t)-2\bar{\Phi}(t)|}{[\bar{\Phi}(t)]^{\theta}}\right).$

Let $B(t)=[\bar{\Phi}(t)]^{-\theta}|\bar{W}_{p}(t)-2\bar{\Phi}(t)|$, and by
the similar arguments as in Section 5.3, we have

$\mathsf{E}[\max_{t\in\mathbb{T}}B(t)]\leq C\sqrt{\log
p}\cdot\max_{t\in\mathbb{T}}\mathsf{E}[B(t)].$

Further,
$\mathsf{E}[B(t)]\leq[\bar{\Phi}(t)]^{-\theta}(\mathsf{E}[\bar{W}_{p}(t)]+2\bar{\Phi}(t))=4[\bar{\Phi}(t)]^{1-\theta}\leq
4$ for $\theta\in(1/2,1]$.

Summing up the above, we have

$P(V_{p,\delta}^{*}>c_{p,\delta}^{*})\leq
C\left(c_{p,\delta}^{*}\right)^{-1}\sqrt{\log p}<\alpha,$

where the last step is by $c^{*}_{p,\delta}=C_{0}\sqrt{\log p}$ with a large
enough constant $C_{0}$.

Next, we demonstrate the upper bound property of $\hat{\pi}_{\delta}^{*}$ with
$\delta(t)=[\bar{\Phi}(t)]^{\theta},\theta\in(1/2,1]$. Similar arguments as in
the proof of Theorem 2.2 for the upper bound can be applied with condition
(10) replaced by condition (11). We omit the details to save space.

## References

  * Arias-Castro et al. (2011) Arias-Castro, E., E. J. Candès, and Y. Plan (2011).  Global testing under sparse alternatives: Anova, multiple comparisons and the higher criticism.  The Annals of Statistics 39(5), 2533–2556.
  * Bradic et al. (2011) Bradic, J., J. Fan, and W. Wang (2011).  Penalized composite quasi-likelihood for ultrahigh dimensional variable selection.  Journal of the Royal Statistical Society: Series B 73(3), 325–349.
  * Buhlmann et al. (2014) Buhlmann, P., K. Kalisch, and L. Meier (2014).  High-dimensional statistics with a view toward applications in biology.  Annu Rev Stat Appl. 1, 255–278.
  * Cai et al. (2013) Cai, T., W. Liu, and Y. Xia (2013).  Two-sample covariance matrix testing and support recovery in high-dimensional and sparse settings.  Journal of the American Statistical Association 108(501), 265–277.
  * Cai and Sun (2017) Cai, T. and W. Sun (2017).  Optimal screening and discovery of sparse signals with applications to multistage high-throughput studies.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 79(1), 197.
  * Cai et al. (2007) Cai, T. T., J. Jin, M. G. Low, et al. (2007).  Estimation and confidence sets for sparse normal mixtures.  The Annals of Statistics 35(6), 2421–2449.
  * Chen (2019) Chen, X. (2019).  Uniformly consistently estimating the proportion of false null hypotheses via lebesgue–stieltjes integral equations.  Journal of Multivariate Analysis 173, 724–744.
  * Efron (2007) Efron, B. (2007).  Size, power and false discovery rates.  The Annals of Statistics 35(4), 1351–1377.
  * Fan et al. (2012) Fan, J., X. Han, and W. Gu (2012).  Estimating false discovery proportion under arbitrary covariance dependence.  Journal of the American Statistical Association 107(499), 1019–1035.
  * Finner and Gontscharuk (2009) Finner, H. and V. Gontscharuk (2009).  Controlling the familywise error rate with plug-in estimator for the proportion of true null hypotheses.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71(5), 1031–1048.
  * Genovese and Wasserman (2004) Genovese, C. and L. Wasserman (2004).  A stochastic process approach to false discovery control.  Annals of Statistics 32, 1035–1061.
  * Jeng and Chen (2019) Jeng, X. J. and X. Chen (2019).  Variable selection via adaptive false negative control in linear regression.  Electronic Journal of Statistics 13(2), 5306–5333.
  * Jeng et al. (2016) Jeng, X. J., Z. J. Daye, W. Lu, and J.-Y. Tzeng (2016).  Rare variants association analysis in large-scale sequencing studies at the single locus level.  PLoS computational biology 12(6), e1004993.
  * Jeng and Hu (2020) Jeng, X. J. and Y. Hu (2020).  Dual control of testing errors in high-dimensional data analysis.  arXiv preprint arXiv:2006.15667.
  * Jeng et al. (2019) Jeng, X. J., T. Zhang, and J.-Y. Tzeng (2019).  Efficient signal inclusion with genomic applications.  Journal of the American Statistical Association 114(528), 1787–1799.
  * Jin (2008) Jin, J. (2008).  Proportion of non-zero normal means: universal oracle equivalences and uniformly consistent estimators.  Journal of the Royal Statistical Society: Series B 70(3), 461–493.
  * Jin and Cai (2007a) Jin, J. and T. T. Cai (2007a).  Estimating the null and the proportion of nonnull effects in large-scale multiple comparisons.  Journal of the American Statistical Association 102(478), 495–506.
  * Jin and Cai (2007b) Jin, J. and T. T. Cai (2007b).  Estimating the null and the proportion of nonnull effects in large-scale multiple comparisons.  Journal of the American Statistical Association 102(478), 495–506.
  * Li and Shao (2002) Li, W. V. and Q.-M. Shao (2002).  A normal comparison inequality and its applications.  Probability Theory and Related Fields 122(4), 494–508.
  * Meinshausen and Rice (2006) Meinshausen, N. and J. Rice (2006).  Estimating the proportion of false null hypotheses among a large number of independently tested hypotheses.  The Annals of Statistics 34(1), 373–393.
  * Storey (2002) Storey, J. D. (2002).  A direct approach to false discovery rates.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 64(3), 479–498.
  * Storey (2003) Storey, J. D. (2003).  The positive false discovery rate: a bayesian interpretation and the q-value.  The Annals of Statistics 31(6), 2013–2035.
  * Westfall and Young (1993) Westfall, P. H. and S. S. Young (1993).  Resampling-based multiple testing: Examples and methods for p-value adjustment, Volume 279.  John Wiley & Sons.