Simultaneous estimation of normal means with side information

TL;DR

This paper develops a data-driven method for estimating normal means using side information, improving accuracy in high-dimensional settings and demonstrating its application in gene expression data integration.

Contribution

It introduces an asymptotically optimal, data-driven decision rule for normal mean estimation leveraging side information, with applications in genomics.

Findings

The proposed rule asymptotically achieves the minimum possible risk.

It can outperform existing methods in numerical experiments.

Demonstrated effectiveness in integrating gene expression datasets.

Abstract

The integrative analysis of multiple datasets is an important strategy in data analysis. It is increasingly popular in genomics, which enjoys a wealth of publicly available datasets that can be compared, contrasted, and combined in order to extract novel scientific insights. This paper studies a stylized example of data integration for a classical statistical problem: leveraging side information to estimate a vector of normal means. This task is formulated as a compound decision problem, an oracle integrative decision rule is derived, and a data-driven estimate of this rule based on minimizing an unbiased estimate of its risk is proposed. The data-driven rule is shown to asymptotically achieve the minimum possible risk among all separable decision rules, and it can outperform existing methods in numerical properties. The proposed procedure leads naturally to an integrative high-dimensional classification procedure, which is illustrated by combining data from two independent gene expression profiling studies.

Tables

Table 1. Table 1: Average total squared errors for the classical normal means problem without side information. GMLEB: the g 𝑔 g -modeling method of Jiang et al. ( 2009 ) .

# nonzero	5				50				500
$\mu$	3	4	5	7	3	4	5	7	3	4	5	7
GMLEB	39	34	23	11	157	105	58	14	459	285	139	18
Proposed	37	32	21	11	158	110	56	14	460	289	133	21

Equations

$$R_n(\mathbf{\theta },\mathbf{\delta })=\frac{1}{n}\sum _{i=1}^{n}E[\{{\theta }_{i1}-{\delta }_i({\mathbf{X}}_{\cdot 1},{\mathbf{X}}_{\cdot 2}){\}}^2]$$

$$\mathcal{S}=\{\mathbf{\delta }({\mathbf{X}}_{\cdot 1},{\mathbf{X}}_{\cdot 2}):{\delta }_i({\mathbf{X}}_{\cdot 1},{\mathbf{X}}_{\cdot 2})=f(X_{i1},X_{i2})\},$$

$$p(x_1,x_2;t_1,t_2)=$$

$$p_i^0(x_1,x_2)=$$

$$f^{\star }(x_1,x_2)=\frac{{\sum }_{j=1}^n{\theta }_{j1}{p}_j^0(x_1,x_2)}{{\sum }_{j=1}^n{p}_j^0(x_1,x_2)}.$$

$$G_n(t_1,t_2)=\frac{1}{n}\sum _{i=1}^{n}I({\theta }_{i1}\le t_1,{\theta }_{i2}\le t_2),$$

$$f_{\rho }^{\star }(x_1,x_2)=x_1+\frac{{\sum }_{j=1}^n({\theta }_{j1}-x_1)p_j^0(x_1,x_2)}{\rho +{\sum }_{j=1}^n{p}_j^0(x_1,x_2)}$$

$$f^{\star }(x_1,x_2)=x_1+\frac{{\tilde{p}}^{\prime }(x_1,x_2)}{\tilde{p}(x_1,x_2)},$$

$$\tilde{p}(x_1,x_2)=\int p(x_1,x_2;t_1,t_2)d{G}_n(t_1,t_2)$$

$$f^{\star }(x_1,x_2)=x_1+\frac{\int (t_1-x_1)p(x_1,x_2;t_1,t_2)d{G}_n(t_1,t_2)}{\int p(x_1,x_2;t_1,t_2)d{G}_n(t_1,t_2)},$$

$$\underset{G}{\mathrm{argmax}}\prod _{i=1}^n\int p(X_{i1},X_{i2};t_1,t_2)dG(t_1,t_2).$$

$${\delta }_{\rho i}^{\mathbf{t}}(x_1,x_2)=x_1+\frac{{\sum }_{j=1}^n(t_{j1}-x_1)p(x_1,x_2;t_{j1},t_{j2})}{\rho +{\sum }_{j=1}^{n}p(x_1,x_2;t_{j1},t_{j2})},$$

$$\text{textsc}sure(\mathbf{t})$$

$$=$$

$$\frac{1}{n}\sum _{i=1}^n{\left\{\frac{{\sum }_j(t_{1j}-X_{i1})p(X_{i1},X_{i2};t_{j1},t_{j2})}{\rho +{\sum }_{j}p(X_{i1},X_{i2};t_{j1},t_{j2})}\right\}}^2+{\sigma }_1^2.$$

$${\ell }_n(\mathbf{t})=\frac{1}{n}\sum _{i=1}^n\{{\theta }_{i1}-{\delta }_{\rho i}^{\mathbf{t}}(X_{i1},X_{i2}){\}}^2$$

$$\mathcal{T}=\{\mathbf{t}:|t_{jd}|\le C{n}^{1/4-\eta },j=1,\dots ,n,d=1,2\}.$$

$$\underset{n\to \infty }{\lim }E\underset{\mathbf{t}\in \mathcal{T}}{\sup }|\text{textsc}sure(\mathbf{t})-{\ell }_n(\mathbf{t})|=0.$$

$${\mathbf{\delta }}_{\rho }^{\hat{\mathbf{t}}}\text{mbox}asin\text{eqref}eq:delt{a}^{t}with\hat{\mathbf{t}}=\underset{\mathbf{t}\in \mathcal{T}}{\mathrm{argmin}}\text{textsc}sure(\mathbf{t}).$$

$$\hat{\mathcal{T}}=\prod _{i=1}^n[X_{i1}-M{\sigma }_1,X_{i1}+M{\sigma }_1]\times [X_{i2}-M{\sigma }_2,X_{i2}+M{\sigma }_2]$$

$$\frac{\lambda }{\lambda +{\sigma }_1^2}{X}_{i1}+\frac{{\sigma }_1^2}{\lambda +{\sigma }_1^2}h(X_{i2})$$

$$I\left(\sum _{i=1}^n{\hat{\theta }}_i{G}_i/{\hat{s}}_i\ge c\right)$$

$${\delta }_i^{\star }({\mathbf{X}}_{\cdot 1})=\frac{{\sum }_{j=1}^n{\theta }_{j1}\varphi \{(X_{i1}-{\theta }_{j1})/{\sigma }_1\}/{\sigma }_1}{{\sum }_{j=1}^n\varphi \{(X_{i1}-{\theta }_{j1})/{\sigma }_1\}/{\sigma }_1},$$

$${\delta }_i^{\star }({\mathbf{X}}_{\cdot 1},{\mathbf{X}}_{\cdot 2})=\frac{{\sum }_{j=1}^n{\theta }_{j1}\exp \{-({\mathbf{X}}_{i\cdot }-{\mathbf{\theta }}_{i\cdot }{)}^{\top }{\Sigma }^{-1}({\mathbf{X}}_{i\cdot }-{\mathbf{\theta }}_{i\cdot })/2\}}{{\sum }_{j=1}^n\exp \{-({\mathbf{X}}_{i\cdot }-{\mathbf{\theta }}_{i\cdot }{)}^{\top }{\Sigma }^{-1}({\mathbf{X}}_{i\cdot }-{\mathbf{\theta }}_{i\cdot })/2\}},$$

$$\Sigma =\left(\begin{array}{cc}{\sigma }_1^2& r{\sigma }_1{\sigma }_2\\ r{\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right).$$

$$R_n(\mathbf{\theta },\mathbf{\delta })=$$

$$\frac{2}{n}\sum _{i=1}^n\int \{{\theta }_{i1}-{\delta }_i^{\star }(x_1,x_2)\}\{{\delta }_i^{\star }(x_1,x_2)-{\delta }_i(x_1,x_2)\}{p}_i^0(x_1,x_2)d{x}_{1}d{x}_2+$$

$$\frac{1}{n}\sum _{i=1}^n\int \{{\delta }_i^{\star }(x_1,x_2)-{\delta }_i(x_1,x_2){\}}^2{p}_i^0(x_1,x_2)d{x}_{1}d{x}_2,$$

$$R_n(\mathbf{\theta },\mathbf{\delta })=$$

$$\ge$$

$$R_n(\mathbf{\theta },{\mathbf{\delta }}_{\rho }^{\star })=$$

$$\frac{2}{n}\sum _{i=1}^n\int \frac{{\sum }_j({\theta }_{j1}-x_1)p_j^0(x_1,x_2)}{\rho +{\sum }_j{p}_j^0(x_1,x_2)}({\theta }_{i1}-x_1)p_i^0(x_1,x_2)d{x}_{1}d{x}_2+$$

$$\frac{1}{n}\sum _{i=1}^n\int {\left\{\frac{{\sum }_j({\theta }_{j1}-x_1)p_j^0(x_1,x_2)}{\rho +{\sum }_j{p}_j^0(x_1,x_2)}\right\}}^2{p}_i^0(x_1,x_2)d{x}_{1}d{x}_2.$$

$$\frac{2}{n}\sum _{i=1}^n\int \frac{{\sum }_j({\theta }_{j1}-x_1)p_j^0(x_1,x_2)}{\rho +{\sum }_j{p}_j^0(x_1,x_2)}({\theta }_{i1}-x_1)p_i^0(x_1,x_2)d{x}_{1}d{x}_2$$

Full text

1 Introduction

Methods for the integrative analysis of multiple datasets are becoming increasingly important. This is especially true in genetics and genomics, where petabytes of public data are readily available for integrative analysis (Richardson et al., 2016; Ritchie et al., 2015). For example, Pickrell et al. (2016) analyzed summary statistics from genome-wide association studies of 42 human traits and found that multiple traits were influenced by several hundred common genetic variants. In a cross-species example, Shpigler et al. (2017) combined results from a honey bee gene expression study with a database of autism-associated genetic variants and found evidence for evolutionary conservation of genes associated with both honey bee sociality and human autism spectrum disorder. Comparing and contrasting existing data, or combining them with new data, can lead to novel insights that would have been difficult or impossible to uncover with a single dataset alone (Tseng et al., 2015).

Integrative analysis strategies can take many forms, and one particularly common implementation is to leverage side information from one or several auxiliary studies for the purpose of improving the analysis of some primary dataset of interest. Examples abound in the multiple testing literature, where methods such as $p$-value weighting and false discovery rate regression incorporate auxiliary information to improve the power to detect true signals in a primary dataset (Genovese et al., 2006; Ramdas et al., 2017). In the genomic risk prediction literature, Hu et al. (2017) and Zhao (2017) showed that summary statistics from previously conducted genome-wide association studies can be used to improve the performance of polygenic risk scores.

Growing interest in these ideas gives rise to an important statistical question: what is the best way to leverage side information? This paper studies this question in a simple but nontrivial problem: the simultaneous estimation of a vector of normal means. The classical version of this problem considers a sequence of independent $X_{i1}\sim N(\theta_{i1},\sigma_{1}^{2})$ for $i=1,\ldots,n$ with known $\sigma_{1}^{2}$, where the goal is to estimate the $\theta_{i1}$ (Johnstone, 2017). The integrative version, studied here, investigates how a auxiliary sequence of Gaussian random variables can be used to improve estimation of the means $\theta_{i1}$ of the primary Gaussian sequence.

This classical Gaussian sequence model is simplistic, but studying data integration in this setting is nevertheless instructive. First, the model is still important for many applications (Cai, 2012; Johnstone, 2017). Second, more accurate estimation of the mean vector has immediate implications for high-dimensional classification in genomics (Greenshtein and Park, 2009), which will be demonstrated in Section 6. Finally, this simple problem can reveal general statistical phenomena that arise in integrative data analysis. More complicated variations of the Gaussian sequence model have been studied, for example involving unknown variances that differ across different indices $i$; see Section 2.2. Extensions of the present work to these more realistic settings are important directions for future work.

Section 2 formalizes this integrative estimation task as a compound decision problem and summarizes previous related work. The optimal way to leverage side information is derived in Section 3, which presents an oracle integrative decision rule that achieves the best risk within a certain class of estimators. This section also introduces a regularized version of the oracle rule that has the same asymptotic risk. A data-driven estimate of this regularized oracle rule is introduced in Section 4, and is shown to asymptotically achieve the optimal risk. Its good performance is illustrated in simulations in Section 5 and in two genomic risk prediction problems in Section 6. A discussion is presented in Section 7 and additional simulations and proofs can be found in the Appendix.

2 Normal means problem with side information

2.1 Problem statement

As in the classical Gaussian sequence problem, consider a sequence of independent $X_{i1}\sim N(\theta_{i1},\sigma_{1}^{2})$ for $i=1,\ldots,n$, with $\sigma_{1}^{2}$ known. The side information problem studied in this paper further supposes that a second sequence of independent $X_{i2}\sim N(\theta_{i2},\sigma_{2}^{2}),i=1,\ldots,n$ is available, with $\sigma_{2}^{2}$ known. The goal is to estimate the $\theta_{i1}$, just as in the classical problem, but here both the $X_{i1}$ and the $X_{i2}$ can be used for estimation. In this sense, the $X_{i1}$ play the role of a primary dataset, and the $X_{i2}$ act as side information from an auxiliary dataset. This paper assumes that the $X_{i1}$ are independent of the $X_{i2}$ for each $i$, though extensions to dependent $X_{i1}$ and $X_{i2}$ are discussed in Section 7.

This formulation is motivated by applications in integrative genomics. The indices $i$ represent different genomic features, such as different genes, and the $X_{i1}$ and $X_{i2}$ represent different measurements on feature $i$ from different studies. For example, in the genomics classification problem described in Section 6, each $X_{i1}$ estimates a classifier parameter $\theta_{i1}$ corresponding to the $i$th gene from a primary study of interest, and each $X_{i2}$ is the $Z$-score for the $i$th gene reported by an auxiliary study of a related phenotype. The goal is to improve classification accuracy in the primary study by leveraging both $X_{i1}$ and $X_{i2}$ to better estimate the $\theta_{i1}$.

In the above example, the $X_{i1}$ and $X_{i2}$ are paired for each $i$, as both correspond to the same genomic feature. The informativeness of this pairing is crucial for the good performance of data integration. For example, because the phenotypes considered by the two studies in Section 6 are related, genes with significant $Z$-scores in the auxiliary study are also likely to be important features for classification in the primary study, so combining the studies is likely to be fruitful. In contrast, if the phenotypes were unrelated, $X_{i2}$ would likely not be informative about $\theta_{i1}$. The challenge is to develop an estimation procedure that can make optimal use of $X_{i2}$, incorporating them when appropriate and discarding them otherwise. This is addressed by the method proposed in this paper.

To more formally state the problem, define $\mathbf{X}_{\cdot d}=(X_{1d},\ldots,X_{nd})$, $\boldsymbol{\theta}_{\cdot d}=(\theta_{1d},\ldots,\theta_{nd})$ for $d=1,2$, and $\boldsymbol{\theta}=(\boldsymbol{\theta}_{\cdot 1},\boldsymbol{\theta}_{\cdot 2})$. Then the normal means problem with side information is to find a decision rule $\boldsymbol{\delta}(\mathbf{X}_{\cdot 1},\mathbf{X}_{\cdot 2})=\{\delta_{1}(\mathbf{X}_{\cdot 1},\mathbf{X}_{\cdot 2}),\ldots,\delta_{n}(\mathbf{X}_{\cdot 1},\mathbf{X}_{\cdot 2})\}:\mathbb{R}^{2n}\rightarrow\mathbb{R}^{n}$ that minimizes the risk function

[TABLE]

over some class of decision rules. An important class, namely the class of separable estimators, will be considered in this paper and is discussed in Section 3. This paper adopts the frequentist framework where the $\boldsymbol{\theta}_{\cdot d}$ are fixed nonrandom constants. The auxiliary data are thus statistically independent of the primary data, and it is interesting that they can still provide useful information for estimating $\boldsymbol{\theta}_{\cdot 1}$.

To illustrate the complexities of this problem, first suppose that it were known that $\theta_{i2}=\theta_{i1}$ for all $i=1,\ldots,n$ and that $\sigma_{1}=\sigma_{2}$. The best way to integrate the auxiliary dataset would clearly be to apply existing optimal estimation methods for a single Gaussian sequence to the sequence of averaged observations $(X_{i1}+X_{i2})/2$. Next consider a slightly more complicated setting: $\theta_{i2}=\theta_{i1}$ for all but one $i$, but the $i$ for which $\theta_{i2}\neq\theta_{i1}$ is unknown. The auxiliary sequence is clearly still informative for estimating the $\boldsymbol{\theta}_{\cdot 1}$, but how it should be used is no longer obvious. Finally, consider an even more complicated scenario where $\theta_{i2}=h(\theta_{i1})+e_{i}$ for some unknown function $h(t)$, where the $e_{i}$ are unknown perturbations that exhibit no patterns with respect to $\theta_{i1}$. If the magnitudes of the $e_{i}$ are small relative to the $\theta_{i1}$, $\mathbf{X}_{\cdot 2}$ should still be useful when estimating $\boldsymbol{\theta}_{\cdot 1}$, but it is even less clear how to optimally integrate it into the estimation procedure. This paper provides one approach.

2.2 Previous work

The classical normal means estimation problem without side information, which aims to minimize the risk function (1) using decision rules that can depend only on $\mathbf{X}_{\cdot 1}$ and not $\mathbf{X}_{\cdot 2}$, has inspired an enormous literature (Johnstone, 2017). Stein (1956) found that the maximum likelihood estimator $\delta_{i}(\mathbf{X}_{\cdot 1})=X_{i1}$ is inadmissible, and since then research has focused on finding alternative estimators with better risk properties. Several different but intimately related perspectives on this problem have been developed.

The shrinkage perspective is exemplified by the James-Stein estimator (James and Stein, 1961; Stigler, 1990), which estimates $\theta_{i1}$ by scaling $X_{i1}$ towards zero. The empirical Bayes perspective (Robbins, 1964) treats the $\theta_{i1}$ as random draws from a prior distribution, uses the $X_{i1}$ to estimate any unknown parameters in the prior, then estimates each $\theta_{i1}$ by its posterior mean conditional on $X_{i1}$. Efron and Morris (1973) showed that the James-Stein estimator is an empirical Bayes estimator assuming a normal prior for the $\theta_{i1}$. The compound decision perspective (Robbins, 1951; Zhang, 1997) treats the $\theta_{i1}$ as nonrandom constants and directly derives the decision rule that minimizes the risk. Under certain conditions, the optimal solution from this perspective is closely related to nonparametric empirical Bayes estimators (Brown and Greenshtein, 2009; Jiang et al., 2009; Zhang, 2003).

More complicated versions of the classical normal means problem have also been intensely studied. For example, specialized methods have been developed for estimating sparse normal means, where most of the $\theta_{i1}$ are assumed to equal zero (Castillo et al., 2012; Donoho and Johnstone, 1994, 1995; Martin et al., 2014). Heteroscedastic normal sequences, where the $X_{i1}$ can have different variances for different indices $i$, have also been considered, both when the variances are known (Fu et al., 2019; Tan, 2016; Weinstein et al., 2018; Xie et al., 2012; Zhang and Bhattacharya, 2017) and when they are unknown but estimates are available (Feng and Dicker, 2018; Gu and Koenker, 2017; Jing et al., 2016).

So far, however, most work on the normal means problem and its variants has considered only a single sequence of observations $X_{i1}$, and it appears that the side information problem (1) has not yet been widely studied. Jiang et al. (2010), Cohen et al. (2013), Tan (2016), and Kou and Yang (2017) proposed methods that can integrate $X_{i2}$, but these essentially require knowledge of the nature of the relationship between $\theta_{i1}$ and $X_{i2}$, and may not work well when this relationship is misspecified. Banerjee et al. (2018) studied the side information problem, but only for sparse $\boldsymbol{\theta}_{\cdot 1}$. Very recently Saha and Guntuboyina (2017) and Koudstaal and Yao (2018) considered two or more Gaussian sequences, but minimized the risk of estimating the means of all of the sequences, rather than the means of just one of them as in (1). In contrast to existing work, this paper studies the optimal use of $\mathbf{X}_{\cdot 1}$ and $\mathbf{X}_{\cdot 2}$ for estimating possibly non-sparse $\boldsymbol{\theta}_{\cdot 1}$.

3 Oracle integrative separable rules

Without any restrictions, the optimal decision rule is simply $\delta_{i}(\mathbf{X}_{\cdot 1},\mathbf{X}_{\cdot 2})=\theta_{i1}$, which is not useful because the performance of this rule cannot realistically be achieved using the observed data alone. Instead, this paper only considers rules in the class

[TABLE]

where $f$ is some fixed real-valued function that is applied to each pair $(X_{i1},X_{i2})$ in order to estimate $\theta_{i1}$. In other words, the estimate of $\theta_{i1}$ is calculated by applying $f(x_{1},x_{2})$ to only the $i$th pair of observations $(X_{i1},X_{i2})$, and $f(x_{1},x_{2})$ cannot vary with $i$.

Rules in $\mathcal{S}$, called “separable” rules, are appealing because of their simplicity and have been extensively studied (Brown and Greenshtein, 2009; Cai, 2012; Robbins, 1951; Zhang, 2003). The maximum likelihood estimator $\boldsymbol{\delta}(\mathbf{X}_{\cdot 1},\mathbf{X}_{\cdot 2})=X_{i1}$ belongs to $\mathcal{S}$, and the James-Stein estimator approximates the optimal separable rule that is linear in $X_{i1}$ (Jiang et al., 2009). The minimum risk among all separable estimators has been shown to be asymptotically equivalent, in a certain sense, to the minimum achievable risk over the larger class of permutation invariant estimators (Greenshtein and Ritov, 2009).

The following proposition describes the oracle optimal integrative rule in $\mathcal{S}$ for estimating $\boldsymbol{\theta}_{\cdot 1}$, which will perform no worse than any separable rule that relies only on $\mathbf{X}_{\cdot 1}$. It is a direct consequence of the fundamental theorem of compound decision problems (Robbins, 1951; Jiang et al., 2009). Let $\phi(x)$ denote the standard normal density and define

[TABLE]

so that the density of $(X_{i1},X_{i2})$ can be abbreviated by $p_{i}^{0}(x_{1},x_{2})$. As mentioned in the problem statement in Section 2.1, this paper assumes that the $X_{i1}$ and $X_{i2}$ are independent, but the following result is easily extended to settings where $X_{i1}$ and $X_{i2}$ are correlated; see Section 7.

Proposition 1

Define the decision rule $\boldsymbol{\delta}^{\star}=(\delta^{\star}_{1},\ldots,\delta^{\star}_{n})$ where $\delta^{\star}_{i}(\mathbf{X}_{\cdot 1},\mathbf{X}_{\cdot 2})=f^{\star}(X_{i1},X_{i2})$ and

[TABLE]

Then $R_{n}(\boldsymbol{\theta},\boldsymbol{\delta})\geq R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star})$ for any $\boldsymbol{\delta}\in\mathcal{S}$ (2) for $R_{n}(\boldsymbol{\theta},\boldsymbol{\delta})$ in (1).

The oracle rule $\boldsymbol{\delta}^{\star}$ also has a useful interpretation as a Bayes rule. If the $\theta_{i1}$ are viewed as independent draws from the discrete prior distribution

[TABLE]

then the posterior expectation $E(\theta_{i1}\mid X_{i1},X_{i2})$ of $\theta_{i1}$ is exactly equal to (4). This is an example of the close connection between compound decision problems and nonparametric empirical Bayes procedures. The dependence between $\theta_{i1}$ and $\theta_{i2}$ under $G_{n}$ quantifies the amount of information that can be borrowed from $X_{i2}$.

While appealing, this Bayesian interpretation is not necessary for Proposition 1, which holds for fixed and constant $\boldsymbol{\theta}_{\cdot 1}$ and $\boldsymbol{\theta}_{\cdot 2}$. Interestingly, under this frequentist setting Proposition 1 shows that $\mathbf{X}_{\cdot 2}$ can improve the estimation of $\boldsymbol{\theta}_{\cdot 1}$ even though $\mathbf{X}_{\cdot 1}$ and $\mathbf{X}_{\cdot 2}$ are statistically independent, as long the sequences $\boldsymbol{\theta}_{\cdot 1}$ and $\boldsymbol{\theta}_{\cdot 2}$ are related in some sense. There need not be an obvious functional relationship between the two mean vectors.

The above view of side information is slightly different from that of existing frameworks. Previous methods (Jiang et al., 2010; Kou and Yang, 2017; Tan, 2016) posit some functional relationship, typically linear, between $\theta_{i1}$ and the observed $X_{i2}$, rather than between $\theta_{i1}$ and the true mean $\theta_{i2}$. For example, Kou and Yang (2017) assume that $\theta_{i1}=h(X_{i2})+e_{i}$ for some error term $e_{i}$, where $h(x)$ must be known up to a finite-dimensional parameter. These methods treat the $X_{i2}$ as fixed, while the proposed framework acknowledges that the $X_{i2}$ are random variables. The difference between existing work and the present setting is akin to the difference between classical regression methods and those that take into account covariate measurement error.

Figure 1 illustrates the oracle rule $\boldsymbol{\delta}^{\star}$ (4) and compares it to the best separable estimator that does not use $\mathbf{X}_{\cdot 2}$, which is the posterior expectation of $\theta_{i1}$ under the prior $G_{n}$ conditional only on $X_{i1}$ (Zhang, 2003). In both panels, $\theta_{i1}$ was generated by drawing $n=20$ values from a standard normal distribution. In the left panel, $\boldsymbol{\theta}_{\cdot 2}=\boldsymbol{\theta}_{\cdot 1}$, so $\mathbf{X}_{\cdot 2}$ was highly informative for $\boldsymbol{\theta}_{\cdot 1}$. Thus $f^{\star}(x_{1},3)$ gave the best estimates of $\theta_{i1}$ for large $X_{i1}$ and $f^{\star}(x_{1},-3)$ was most accurate for small $X_{i1}$. In the right panel of Figure 1, the $\theta_{i2}$ were generated from an independent standard normal so that $\mathbf{X}_{\cdot 2}$ was completely uninformative. In this non-informative setting, $\boldsymbol{\delta}^{\star}$ may not have the same performance as the optimal non-integrative separable rule for any given set of $\mathbf{X}_{\cdot 1}$ and $\mathbf{X}_{\cdot 2}$, but in expectation Proposition 1 guarantees that it will have equal or lower risk.

The oracle separable integrative rule $\boldsymbol{\delta}^{\star}$ described in (4) cannot be implemented in practice because it requires knowing the true $(\theta_{i1},\theta_{i2})$ up to permutation of the indices. Section 4 will introduce a data-driven rule that targets the performance of $\boldsymbol{\delta}^{\star}$, though for technical reasons it will be more convenient to target a regularized version of the oracle rule. This will be denoted by $\boldsymbol{\delta}^{\star}_{\rho}=(\delta^{\star}_{\rho 1},\ldots,\delta^{\star}_{\rho n})$, with $\delta^{\star}_{\rho i}(\mathbf{X}_{\cdot 1},\mathbf{X}_{\cdot 2})=f^{\star}_{\rho}(X_{i1},X_{i2})$ for

[TABLE]

and $\rho$ a small positive constant that prevents the denominator from being too close to zero. Under some assumptions, $\boldsymbol{\delta}^{\star}_{\rho}$ will have the same asymptotic risk as the oracle $\boldsymbol{\delta}^{\star}$.

Assumption 1

There exist positive constants $C$ and $\eta$ such that $|\theta_{id}|\leq Cn^{1/4-\eta}$ for $i=1,\ldots,n$ and $d=1,2$.

Theorem 1

Under Assumption 1, $\lim_{n\rightarrow\infty}\{R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star}_{\rho})-R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star})\}=0$.

Assumption 1 determines how quickly the magnitudes of $\theta_{id}$ can grow. To put this rate into perspective, if the $\theta_{id}$ were random draws from a normal distribution, then $\max_{i}|\theta_{id}|$ would be $O(\log^{1/2}n)$ almost surely. Related assumptions, which essentially restrict how variable the $\theta_{id}$ can be, have been made in previous work on normal means estimation without side information. For example, Xie et al. (2012) require $\lim n^{-1}\sum_{i}\theta_{i1}^{2}<\infty$, and Jiang et al. (2009) and Zhang (2009) control the rate of the $p$-th weak moment of the distribution function $n^{-1}\sum_{i}I(\theta_{i1}\leq t_{1})$.

4 Data-driven separable estimator

4.1 Existing nonparametric empirical Bayes approach

By Proposition 1 and Theorem 1, the regularized oracle $\boldsymbol{\delta}^{\star}_{\rho}$ (6) is asymptotically optimal within the class of separable estimators (2), but it cannot be implemented in practice. It therefore remains to develop a fully data-driven estimator for the $\theta_{i1}$. Two classes of approaches already exist. They have been termed $f$- and $g$-modeling (Efron, 2014, 2019) and are based on nonparametric empirical Bayes principles that pretend that the $(\theta_{i1},\theta_{i2})$ are random variables with prior distribution $G_{n}(t_{1},t_{2})$ (5).

In $f$-modeling, the oracle estimator (4) would be re-expressed as

[TABLE]

where $\tilde{p}^{\prime}(x_{1},x_{2})=\partial\tilde{p}/\partial x_{1}$ and

[TABLE]

with $p(x_{1},x_{2};t_{1},t_{2})$ from (3). If the $(\theta_{i1},\theta_{i2})$ were truly random, $\tilde{p}(x_{1},x_{2})$ could be interpreted as the marginal density of $(X_{i1},X_{i2})$, and $\tilde{p}(x_{1},x_{2})$ and $\tilde{p}^{\prime}(x_{1},x_{2})$ could be estimated nonparametrically using kernel density estimators. In $g$-modeling, the oracle estimator would be re-expressed as

[TABLE]

and if the $(\theta_{i1},\theta_{i2})$ were truly random, a nonparametric estimate of $G_{n}(t_{1},t_{2})$ could be obtained by maximizing the marginal log-likelihood (Kiefer and Wolfowitz, 1956)

[TABLE]

Both $f$- and $g$-modeling have been used in normal means problems without side information, where they are asymptotically optimal even in the frequentist framework where the $\theta_{i1}$ and $\theta_{i2}$ are nonrandom (Brown and Greenshtein, 2009; Feng and Dicker, 2018; Fu et al., 2019; Jiang et al., 2009; Koenker, 2014; Koenker and Mizera, 2014; Saha and Guntuboyina, 2017; Zhang, 2009). However, neither approach directly estimates the oracle decision rule, with $f$-modeling proceeding through the intermediate quantity $\tilde{p}(x_{1},x_{2})$ and $g$-modeling proceeding through $G_{n}(t_{1},t_{2})$.

4.2 Proposed direct risk minimization approach

This paper explores a more direct approach to estimating the oracle integrative separable classifier. Motivated by the regularized oracle (6), consider separable rules of the form $\boldsymbol{\delta}^{t}_{\rho}=(\delta^{t}_{\rho 1},\ldots,\delta^{t}_{\rho n})$ with

[TABLE]

for a given $\mathbf{t}=(t_{11},\ldots,t_{n1},t_{12},\ldots,t_{n2})$. By Theorem 1, the optimal $\mathbf{t}$ equals $(\theta_{11},\ldots,\theta_{n1},\theta_{12},\ldots,\theta_{n2})$, but the $\theta_{id}$ are not known. The challenge is to choose a $\mathbf{t}$ in a data-driven fashion that can still asymptotically achieve the optimal risk.

Choosing $\mathbf{t}$ to minimize the risk (1) of $\boldsymbol{\delta}^{\mathbf{t}}_{\rho}$ (7) should give an estimator with good performance. However, calculating the risk requires knowing the true $\theta_{id}$. On the other hand, Stein’s Lemma (Stein, 1981) can be used to give an unbiased estimate of the true risk as a function only of $\mathbf{t}$:

[TABLE]

The following theorem shows that $\textsc{sure}(\mathbf{t})$ is also a good approximation to the actual loss

[TABLE]

uniformly over the set

[TABLE]

Theorem 2

Under Assumption 1, if $0<\rho\leq 1$, then

[TABLE]

The tuning parameter $\mathbf{t}$ can now be chosen by minimizing this estimated risk, as a proxy for minimizing the unknown true risk. The proposed estimator is therefore defined to be

[TABLE]

This strategy of direct risk minimization is common in the compound decision literature (Jing et al., 2016; Kou and Yang, 2017; Tan, 2016; Weinstein et al., 2018; Xie et al., 2012, 2016; Zhang and Bhattacharya, 2017), but has not yet been used to approximate an optimal separable rule like (6). The following theorem shows that (11) can asymptotically achieve the same the performance as the optimal separable decision rule.

Theorem 3

Under the same conditions as Theorem 2, $\lim_{n\rightarrow\infty}\{E\ell_{n}(\hat{\mathbf{t}})-R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star})\}\leq 0$, where $E\ell_{n}(\hat{\mathbf{t}})$ is the risk of the proposed estimator $\boldsymbol{\delta}^{\hat{\mathbf{t}}}_{\rho}$ (11).

4.3 Implementation

The proposed estimator has been implemented in the R package cole, available at github.com/sdzhao/cole. In practice, the exact value of $\rho$ appears to make little difference, and $\rho=0$ works well in most cases. When the range of the $X_{id}$ is very large or the variances $\sigma_{d}^{2}$ are very small, problem may arise when calculating $\textsc{sure}(\mathbf{t})$ due to numerical precision, in which case setting $\rho=10^{-12}$ seems to be sufficient. Throughout this paper, the proposed method was implemented with $\rho=0$.

Because the value of $Cn^{1/4-\eta}$ that defines the feasible set $\mathcal{T}$ (10) is not known, in practice the minimization in (11) can be performed over

[TABLE]

for some sufficiently large positive constant $M$, so that $\hat{\mathcal{T}}$ contains $(\boldsymbol{\theta}_{\cdot 1},\boldsymbol{\theta}_{\cdot 2})$ with probability $\Phi(-M)^{n}$, where $\Phi$ is the cumulative distribution function of a standard normal. By default, cole uses $M=5$, so that $\hat{\mathcal{T}}$ contains $(\boldsymbol{\theta}_{\cdot 1},\boldsymbol{\theta}_{\cdot 2})$ with probability 0.99 when $n=10,000$. Optimizing $\textsc{sure}(\mathbf{t})$ over $\hat{\mathcal{T}}$ is sensible because it is known from Theorem 1 that $E\{\textsc{sure}(\mathbf{t})\}$ achieves a global minimum at $t_{jd}=\theta_{jd}$. This method works well, but bridging the gap between the theoretical procedure and its practical implementation is an important direction for future work.

Minimizing $\textsc{sure}(\mathbf{t})$ is difficult because it is a nonconvex function. The implementation in cole performs a simple coordinate descent. At initialization, $t_{id}$ is set to $X_{id}$, and at each iteration one $t_{id}$ is updated by optimizing over $K$ equally spaced candidates in $[X_{id}-M\sigma_{d},X_{id}+M\sigma_{d}]$. By default, cole uses $K=10$, and all analyses in this paper use $K=10$ unless otherwise stated. The coordinates of $\mathbf{t}$ are updated in the order $t_{11},t_{21},\ldots,t_{n1},t_{12},\ldots,t_{n2}$, and convergence is reached when all of the coordinates have been cycled through once without changing the value of $\textsc{sure}(\mathbf{t})$ by more than a small $\epsilon$, which cole sets to $10^{-5}$ by default.

5 Simulations

5.1 Normal means problem without side information

The direct risk minimization approach proposed in this paper for estimating optimal separable decision rules appears to be novel in the compound decision literature. This section thus first illustrates how this idea performs in the classical normal means problem without side information. The optimal separable estimator and its corresponding unbiased risk estimate will look like (4) and (8), respectively, but with the density $p(x_{1},x_{2};t_{1},t_{2})$ replaced by $\phi\{(x_{1}-t_{1})/\sigma_{1}\}/\sigma_{1}$, where $\phi(x)$ is the standard normal density. Similar to (11), a data-driven estimator of the oracle rule can be obtained by minimizing the risk estimate over $\mathbf{t}_{1}$ using the coordinate descent algorithm described in Section 4.3; this is available in the cole package. Analogs of Theorems 1–3 can also be proved.

The direct estimator was compared to the $g$-modeling procedure of Jiang et al. (2009), which is also asymptotically risk-optimal. One independent sequence $X_{i1},i=1,\ldots,1,000$ was generated from N$(\theta_{i1},1)$, with the goal of estimating $\boldsymbol{\theta}_{\cdot 1}$ using only $\mathbf{X}_{\cdot 1}$. The $\theta_{i1}$ equaled either 0 or $\mu$ and the number of nonzero $\theta_{i1}$ equaled either 5, 50, or 500. Table 1 displays the average total squared errors over 100 replications. Results for the estimator of Jiang et al. (2009) were taken directly from their Table 1, while the proposed estimator was implemented using a coordinate descent algorithm that optimized over $K=50$ candidates for each $t_{1j}$. The results show that both estimators had almost identical performance.

5.2 Settings for normal means problem with side information

The primary data $\mathbf{X}_{\cdot 1}=(X_{11},\ldots,X_{n1})$ were generated in four different ways, for three dense and one sparse configuration of their means $\boldsymbol{\theta}_{\cdot 1}$. To generate dense $\boldsymbol{\theta}_{\cdot 1}$, values of $\theta_{i1}$ were independently drawn from either a N$(0,1)$, a Unif$(-2,2)$, or an Exp$(1)$ distribution. To generate the sparse configuration, 10% of the coordinates of $\boldsymbol{\theta}_{\cdot 1}$ were set to 1.5 and the rest were set to 0. The observed primary data were generated as $X_{i1}=\theta_{i1}+\epsilon_{i1}$, where the $\epsilon_{i1}$ were independently drawn from standard normal random variables. The $\theta_{i1}$ were fixed across all replications.

For each of these four settings, the auxiliary data $\mathbf{X}_{\cdot 2}=(X_{12},\ldots,X_{n2})$ were generated in three different ways, to model different degrees of informativeness of $\boldsymbol{\theta}_{\cdot 2}$. First define $e_{i}$ to be independent draws from a Unif$(-4,4)$. To generate strongly, weakly, and non-informative side information, $\theta_{i2}$ was set to be either $2\theta_{i1}^{2}$, $\theta_{i1}^{2}+e_{i}$, or $e_{i}$, respectively. The observed auxiliary data were generated as $X_{i2}=\theta_{i2}+\epsilon_{i2}$, where the $\epsilon_{i}$ were again independently drawn from standard normal random variables. The $\theta_{i2}$ were fixed across all replications.

The proposed integrative normal means estimator (11) was compared to two existing approaches that can incorporate side information. One was the procedure of Banerjee et al. (2018). The other was estimator (1) of Kou and Yang (2017), defined as

[TABLE]

for some function $h(x)$ known up to a finite number of parameters. These unknown parameters, as well as $\lambda$, are chosen by minimizing an unbiased estimate of the risk of this estimator. This estimator is motivated by the regression model $\theta_{i1}=h(X_{i2})+e_{i}$ for some error terms $e_{i}$. However, it can be difficult to choose the correct regression function $h(x)$. For example, in some of the present simulation settings, the true relationship between the primary and auxiliary data is $\theta_{i2}=2\theta_{i1}^{2}+e_{i}$, which is difficult to translate into a regression model of $\theta_{i1}$ on $X_{i2}$. When implementing the method of Kou and Yang (2017), these simulations used both the nonlinear model $\theta_{i1}=\beta_{0}+\beta_{1}|X_{i2}|^{1/2}+e_{i}$ and the linear model $\theta_{i1}=\beta_{0}+\beta_{1}X_{i2}+e_{i}$.

Finally, two additional estimators for $\boldsymbol{\theta}_{\cdot 1}$ were also implemented to provide performance baselines. The first was the oracle (4), which attains the lowest possible risk of any separable decision rule that incorporates side information. The second was the $g$-modeling method of Jiang et al. (2009), which can asymptotically achieve the optimal risk of any separable rule that does not use side information.

5.3 Results for normal means problem with side information

Figure 2 illustrates the average losses, over 200 simulations, achieved by the competing methods for N$(0,1)$, Unif$(-2,2)$, Exp$(1)$, or sparse $\boldsymbol{\theta}_{\cdot 1}$ and non-informative, weakly informative, or strongly informative $\boldsymbol{\theta}_{\cdot 2}$. Comparing the performances of the oracle rule (4) and the method of Jiang et al. (2009) shows that including auxiliary data does not degrade estimation accuracy asymptotically when $\boldsymbol{\theta}_{\cdot 2}$ is non-informative, and can greatly improve it when $\boldsymbol{\theta}_{\cdot 2}$ is informative.

The performance of the proposed data-driven estimator $\boldsymbol{\delta}^{\hat{\mathbf{t}}}_{\rho}$ (11) indeed appeared to converge to the oracle performance as the number of observations $n$ increased, consistent with Theorem 3. Unlike the oracle, however, incorporating non-informative $\mathbf{X}_{\cdot 2}$ in $\boldsymbol{\delta}^{\hat{\mathbf{t}}}_{\rho}$ resulted in worse performance compared to the other methods when $n$ was small. This is expected, as non-informative $X_{i2}$ add extra noise without decreasing bias, and the data-driven method requires enough samples to learn that the $X_{i2}$ are not useful. In contrast, $\boldsymbol{\delta}^{\hat{\mathbf{t}}}_{\rho}$ regained its competitiveness for larger $n$, and when the auxiliary $\mathbf{X}_{\cdot 2}$ were at least weakly informative, it frequently achieved the lowest risk among all methods. These results suggest that incorporating $\mathbf{X}_{\cdot 2}$ using the proposed method is highly effective when $\mathbf{X}_{\cdot 2}$ is informative, and does not do too much harm when it is not.

The proposed $\boldsymbol{\delta}^{\hat{\mathbf{t}}}_{\rho}$ was sometimes outperformed by the two different implementations of the procedure of Kou and Yang (2017), for example when the $\theta_{i1}$ were generated from Exp$(1)$. This may be because this setting was particularly difficult for the proposed method. Out of the four configurations of $\boldsymbol{\theta}_{\cdot 1}$, the maximum value of $|\theta_{i1}|$ was largest under the Exp$(1)$ configuration, and Assumption 1 makes it clear that restricting this maximum value is important for the good performance of $\boldsymbol{\delta}^{\hat{\mathbf{t}}}_{\rho}$. On the other hand, when $n=1,000$ $\boldsymbol{\delta}^{\hat{\mathbf{t}}}_{\rho}$ had essentially the same risk as the methods of Kou and Yang (2017), and for other configurations of $\theta_{i1}$, $\boldsymbol{\delta}^{\hat{\mathbf{t}}}_{\rho}$ could perform significantly better.

Finally, the proposed rule performed extremely well with sparse $\boldsymbol{\theta}_{\cdot 1}$, even though it was not designed for this scenario. When the auxiliary data were strongly informative, it achieved the lowest risks among all implemented methods when $n\geq 200$. It would be interesting to explore extensions of the proposed procedure to estimate sparse normal means.

6 Data analysis

High-dimensional classification is an important problem in genomics. Shi et al. (2010) studied the effectiveness of using gene expression microarray data to develop classification rules for various phenotypes. This section focuses on classification of two of these phenotypes: estrogen receptor status and treatment response status in breast cancer patients. The training and validation datasets they used are publicly available from the Gene Expression Omnibus (Edgar et al., 2002) under accession number GSE20194.

Integrating auxiliary data may help improve classification accuracy. Wang et al. (2005) developed a gene expression signature for distant metastasis-free survival in estrogen receptor-positive and -negative breast cancer patients. It may be possible to leverage data from Wang et al. (2005), publicly available under accession number GSE2034, to more accurately classify the two outcomes from Shi et al. (2010). However, it is not clear how to best integrate these auxiliary data.

The normal means estimation problem using side information, studied in this paper, provides one approach. Greenshtein and Park (2009) showed that minimizing the squared error risk in the normal means problem is closely connected to minimizing the misclassification rate in high dimensional classification. Let $\bar{G}_{i}^{Y}$ denote the average expression level of gene $i$ in across all training subjects in class $Y=0,1$ and $\hat{s}^{Y}_{i}$ denote its estimated standard deviation. Greenshtein and Park (2009) considered classifying an observed gene expression vector $(G_{1},\ldots,G_{n})$ using

[TABLE]

for some cutoff $c$, where $\hat{s}_{i}=\{(\hat{s}^{1}_{i})^{2}/n_{1}+(\hat{s}^{0}_{i})^{2}/n_{0}\}^{1/2}$ and $\hat{\theta}_{i}$ is an estimate of the expected value of $Z_{i}=(\bar{G}^{1}_{i}-\bar{G}^{0}_{i})/\hat{s}_{i}$. They showed that using the $f$-modeling procedure of Brown and Greenshtein (2009) to obtain $\hat{\theta}_{i}$ can lead to more accurate classification compared to simply using $\hat{\theta}_{i}=Z_{i}$.

Combined with ideas in this paper, this framework leads to a natural integrative classifier. Let $X_{i1}$ equal $Z_{i}$ calculated for either estrogen receptor status or treatment response status from the Shi et al. (2010) study, and let $X_{i2}$ be the differential expression $Z$-score of the $i$th gene with respect to either estrogen receptor status or distant metastasis-free survival from the Wang et al. (2005) study. Integrating $X_{i2}$ into the estimate $\hat{\theta}_{i1}$ should lead to more accurate classification.

This integrative classification was implemented using the proposed rule $\boldsymbol{\delta}^{\hat{\mathbf{t}}}_{\rho}$ (11), the method of Kou and Yang (2017) using a model linear in $X_{i2}$, and the procedure of Banerjee et al. (2018) for sparse normal means. These were compared to five classifiers that do not make use of auxiliary information: 1) the method of Greenshtein and Park (2009) but implemented using the $g$-modeling procedure of Jiang et al. (2009), 2) the naive Bayes classifier, 3) logistic lasso using the R package glmnet (Friedman et al., 2010), 4) random forest using the R package ranger (Wright and Ziegler, 2017), and 5) the regularized optimal affine discriminant analysis of Fan et al. (2012) using the R package TULIP (Pan et al., 2019). Tuning parameters for lasso and the method of Fan et al. (2012) were chosen using 10-fold cross-validation while random forest was run using default parameters.

The integrative, naive Bayes, and Greenshtein and Park (2009) classifiers all assume that the $X_{id}$ are independent across $i$. For these procedures, screening was thus first performed to ensure that the magnitude of the correlation between every pair of genes in the training data was small, similar to what was done in Dicker and Zhao (2016). Specifically, genes were sorted from most to least significantly associated with the outcome in the training data, with $p$-values calculated using the R package limma (Smyth, 2005). Starting from the most significant gene, any other gene with correlation greater than 0.2 in magnitude was removed from the dataset. No screening was performed for lasso, random forest, or the method of Fan et al. (2012).

Misclassification rates for estrogen receptor and treatment response status were assessed using the same training and testing datasets used in Shi et al. (2010), and classification was also repeated after swapping the roles of the training and testing data. The averages of the two resulting misclassification rates for the different methods are displayed in Figure 3.

The results suggest that integrative classification can be a useful strategy. Intuitively, the survival results from Wang et al. (2005) should be most informative for predicting treatment response, while the ER status data from Wang et al. (2005) should be most useful for predicting ER status. Indeed, the proposed integrative classifier using survival $Z$-scores to predict treatment response gave the lowest misclassification rate among all methods. The proposed method integrating ER status $Z$-scores to predict ER status performed better than every method except random forest and lasso.

7 Discussion

This paper assumes that the primary data $\mathbf{X}_{\cdot 1}$ and the auxiliary data $\mathbf{X}_{\cdot 2}$ are statistically independent. However, in some practical settings $X_{i1}$ and $X_{i2}$ may be correlated for each $i$, for example if $\mathbf{X}_{\cdot 1}$ and $\mathbf{X}_{\cdot 2}$ arise from case-control studies with shared controls (Zaykin and Kozbur, 2010). The ideas proposed in this paper can be naturally extended to this correlated setting. Assuming $(X_{i1},X_{i2})$ were bivariate normal with a known correlation, the oracle integrative rule would be similar to (4) and is given in (13) in the Appendix. An asymptotically risk-optimal data-driven estimator could then be constructed by minimizing an unbiased risk estimate derived using Stein’s lemma.

As pointed out by one referee, this setting is especially interesting because when $X_{i1}$ and $X_{i2}$ are correlated, the $\mathbf{X}_{\cdot 2}$ provides useful information for estimating $\boldsymbol{\theta}_{\cdot 1}$ even when $\boldsymbol{\theta}_{\cdot 2}$ and $\boldsymbol{\theta}_{\cdot 1}$ are completely unrelated. This is not true when $X_{i1}$ and $X_{i2}$ are independent. This is verified by Figure 4 in the Appendix, where the oracle integrative rule performed better than the oracle non-integrative rule when $|cor(X_{i1},X_{i2})|=0.9$ even though $\boldsymbol{\theta}_{\cdot 2}$ was generated to be non-informative for $\boldsymbol{\theta}_{\cdot 1}$. Thus rules such as (4) can therefore take full advantage of information about $\boldsymbol{\theta}_{\cdot 1}$ contained the auxiliary $\mathbf{X}_{\cdot 2}$, whether that information comes in the form of informative $\theta_{i2}$, correlated $X_{i2}$, or both.

This paper considered only a single sequence of auxiliary data, but it is straightforward to extend the proposed procedure to multiple auxiliary sequences. However, this would result in theoretical and computational difficulties, for given $D-1$ auxiliary datasets, Assumption 1 would require $|\theta_{id}|\leq n^{1/(2D)-\eta}$ for $d=1,\ldots,D$, and the proposed procedure would require optimizing over $Dn$ parameters. It would be of great interest to study whether there exists a convex surrogate of the unbiased risk estimate (8). An alternative approach might be to use parametric or semiparametric methods, such as those proposed by Kou and Yang (2017), but to endow them with data-driven model selection capabilities.

It would be interesting to extend data integration ideas to other variants of the classical normal means problem, such as heteroscedastic sequences, sparse sequences, and non-normal observed data. It would also be interesting to consider broader applications of the compound decision framework beyond the simultaneous estimation of a mean vector, such as the integrative high-dimensional classification problem in Section 6.

Though this paper studied the highly stylized problem of normal means estimation with side information, its results reveal several general principles of integrative analysis. First, auxiliary data can be useful even if they are statistically independent of, and have no clearly expressible functional relationship with, primary data. The two datasets need only be related in the sense discussed in Section 3. Second, in principle, integrating auxiliary data can only help and not harm the primary analysis. This is because it is possible to learn from the data themselves the degree to which the auxiliary data are informative, and thus the degree to which they should influence inference in the primary data. Third, nonparametric methods, such as the proposed (11), can asymptotically achieve ideal performance.

Acknowledgments

The author thanks Professors Jiaying Gu, Gourab Mukherjee, and Roger Koenker for helpful discussions. This research was supported in part by NSF grant DMS-1613005.

Appendix A Simulations with correlated $X_{i1}$ and $X_{i2}$

When $X_{i1}$ and $X_{i2}$ are correlated, $\mathbf{X}_{\cdot 2}$ can provide useful information for estimating $\boldsymbol{\theta}_{\cdot 1}$ even if $\boldsymbol{\theta}_{\cdot 2}$ and $\boldsymbol{\theta}_{\cdot 1}$ are completely unrelated. This section illustrates this phenomenon in simulations that compare the oracle non-integrative rule

[TABLE]

where $\phi(x)$ is the standard normal density, to the oracle integrative rule when $(X_{i1},X_{i2})$ is bivariate normal with correlation $r$:

[TABLE]

with $\mathbf{X}_{i\cdot}=(X_{i1},X_{i2})^{\top}$, $\boldsymbol{\theta}_{i\cdot}=(\theta_{i1},\theta_{i2})^{\top}$, and

[TABLE]

The $\theta_{i1}$ and non-informative $\theta_{i2}$ were generated as in Section 5, and the $X_{id}=\theta_{id}+\epsilon_{id}$ for $d=1,2$, where $(\epsilon_{i1},\epsilon_{i2})$ are bivariate normal with mean zero and variance $\Sigma$. Figure 4 reports the average losses of these two oracle estimators over 200 simulations for $r=-0.9,0,0.9$ and shows that the integrative estimator significantly outperformed the non-integrative one when $r\neq 0$.

Appendix B Proof of Proposition 1

For any separable rule $\boldsymbol{\delta}=(\delta_{1},\ldots,\delta_{n})\in\mathcal{S}$ (2),

[TABLE]

where$p_{i}^{0}(x_{1},x_{2})$ equals the density of $(X_{i1},X_{i2})$ as defined in (3). Since $\boldsymbol{\delta}$ and $\boldsymbol{\delta}^{\star}$ are separable, the functions $\delta_{i}(x_{1},x_{2})$ and $\delta^{\star}_{i}(x_{1},x_{2})$ do not depend on the index $i$. Thus for the $\delta^{\star}_{i}$ given in (4), the middle term above equals zero, so

[TABLE]

for any $\boldsymbol{\delta}\in\mathcal{S}$. Therefore $\boldsymbol{\delta}^{\star}$ achieves the minimum risk within the class of separable rules.

Appendix C Proof of Theorem 1

Denoting the density of $(X_{i1},X_{i2})$ by $p_{i}^{0}(x_{1},x_{2})$ as in (3), the risk of the regularized oracle rule $\boldsymbol{\delta}^{\star}_{\rho}$ (6) can be written as

[TABLE]

The second term above obeys

[TABLE]

since $\rho>0$. Therefore

[TABLE]

Similarly, the risk of the unregularized oracle rule $\boldsymbol{\delta}^{\star}$ (4) is

[TABLE]

Since

[TABLE]

it follows that $R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star}_{\rho})-R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star})$ is upper bounded by

[TABLE]

To simplify notation, for any set $\mathcal{C}\subseteq\mathbb{R}^{2}$, define

[TABLE]

so that $R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star}_{\rho})-R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star})\leq\Delta_{n}(\mathbb{R}^{2})$. It will be shown below that $\Delta_{n}(\mathbb{R}^{2})$ tends to zero asymptotically. Since the regularized oracle rule $\boldsymbol{\delta}^{\star}_{\rho}$ is separable (2), $0\leq R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star}_{\rho})-R_{n}(\boldsymbol{\theta},\boldsymbol{\delta}^{\star})$ as well, by Proposition 1. This will conclude the proof.

To upper-bound $\Delta_{n}(\mathbb{R}^{2})$, define the set

[TABLE]

so that, roughly, $\sum_{i}p_{i}^{0}(x_{1},x_{2})$ is large on $\mathcal{L}_{n}$. Then $\Delta_{n}(\mathbb{R}^{2})=\Delta_{n}(\mathcal{L}_{n})+\Delta_{n}(\mathcal{L}_{n}^{c})$, and

[TABLE]

where the second inequality follows from Jensen’s inequality. It remains to show that $\Delta(\mathcal{L}_{n}^{c})$ goes to zero. To this end, further define the sets

[TABLE]

It will be shown that $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{1})$, $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{2})$, and $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{3})$ converge to zero.

First,

[TABLE]

Defining the function $g(x)=\exp\{x/(2\sigma_{1}^{2})\}$, Jensen’s inequality gives

[TABLE]

which is always at least 1. Combining this with another application of Jensen’s inequality implies

[TABLE]

Since $|\theta_{j1}|\leq Cn^{1/4-\eta}$ by Assumption 1, $\theta_{i1}-x_{1}\geq 0$ on $\mathcal{A}_{1}$, so

[TABLE]

Now for each value of $x_{2}$, define the function

[TABLE]

so that

[TABLE]

This implies that

[TABLE]

Since $\theta_{j1}-x_{1}>0$ for all $j$ on $\mathcal{A}_{1}$. Therefore, $u_{x_{2}}(x_{1})$ is increasing on $\mathcal{A}_{1}$ and $u_{x_{2}}(-Cn^{1/4-\eta})>0$. Then

[TABLE]

where the last line follows from the Cauchy-Schwarz inequality. Integrating by parts gives

[TABLE]

so

[TABLE]

Finally, first integrate the remaining integral over $x_{2}\leq-Cn^{1/4-\eta}$. Since $|\theta_{j2}|\leq Cn^{1/4-\eta}$ by Assumption 1, it follows that $\theta_{j2}-x_{2}\geq-Cn^{1/4-\eta}-x_{2}\geq 0$ for all $j=1,\ldots,n$ and thus that $(\theta_{j2}-x_{2})^{2}\geq(-Cn^{1/4-\eta}-x_{2})^{2}$. Therefore

[TABLE]

which implies that

[TABLE]

A similar calculation shows that on $Cn^{1/4-\eta}<x_{2}$,

[TABLE]

as well. Finally, on $-Cn^{1/4-\eta}<x_{2}\leq Cn^{1/4-\eta}$,

[TABLE]

Putting these results together shows that $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{1})$ is at most

[TABLE]

A very similar calculation can be used to show that $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{3})\rightarrow 0$ as well.

It remains to bound $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{2})$. Denote $u(x_{1},x_{2})=\sum_{j}p_{j}^{0}(x_{1},x_{2})$ and use (15) see that $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{2})$ is at most

[TABLE]

Now for any constant $K$,

[TABLE]

This implies that the function $u^{1/2}\log(K/u)$ increases for $u\leq Ke^{-2}$ and decreases thereafter. Now consider the function $2u^{1/2}$, which equals $u^{1/2}\log(K/u)$ at $u=Ke^{-2}$ and increases thereafter. Therefore $u^{1/2}\log(K/u)$ is upper-bounded by the nondecreasing function $u^{1/2}\{2\vee\log(K/u)\}$. Therefore $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{2})$ is at most

[TABLE]

To finish bounding $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{2})$, integrate over $x_{2}\leq-Cn^{1/4-\eta}$, $-Cn^{1/4-\eta}<x_{2}\leq Cn^{1/4-\eta}$, and $Cn^{1/4-\eta}<x_{2}$, as was done when bounding $\Delta_{n}(\mathcal{L}_{n}^{c}\cap\mathcal{A}_{1})$. Using Assumption 1, it can be shown that

[TABLE]

which means that

[TABLE]

Appendix D Proof of Theorem 2

D.1 Outline

For the parameter vector $\mathbf{t}=(t_{11},\ldots,t_{n1},t_{12},\ldots,t_{n2})$, define the function

[TABLE]

Then it is straightforward to show that

[TABLE]

By Jensen’s inequality,

[TABLE]

Next, the variable $\sum_{i}(X_{i1}-\theta_{i1})^{2}/\sigma_{1}^{2}$ has a $\chi^{2}_{n}$ distribution, which has mean $n$ and variance $2n$. Therefore

[TABLE]

The remainder of the proof uses empirical process techniques to show that

[TABLE]

D.2 Symmetrization

The terms $f(X_{i1},X_{i2},\theta_{i1};\mathbf{t})$ are independent but not identically distributed. Nevertheless, the proof of the symmetrization Lemma 2.3.1 of van der Vaart and Wellner (1996) still holds and implies that

[TABLE]

where the $\varepsilon_{i}$ are Rademacher random variables independent of the $X_{i1}$ and $X_{i2}$. Using Stein (1981), it can be shown that $Ef(X_{i1},X_{i2},\theta_{i1};\mathbf{t})=0$. Therefore it only remains to show that

[TABLE]

D.3 Truncation

This section first constructs a function $F_{n}(x_{1},\mu)$ that satisfies $|f(x_{1},x_{2},\mu;\mathbf{t})|\leq F_{n}(x_{1},\mu)$ as well as $E\{F_{n}(X_{i1},\mu)\}<\infty$, for $f$ defined in (16). By the triangle inequality,

[TABLE]

Using Jensen’s inequality as in (15), it follows that

[TABLE]

Let $y=p(x_{1},x_{2};t_{1},t_{2})$ and $C=(2\pi\sigma_{1}\sigma_{2})^{-1}n$ such that this upper bound can be written as

[TABLE]

where $0\leq y\leq C$ by the definition of $p(x_{1},x_{2};t_{1},t_{2})$ in (3). Since

[TABLE]

and the function $\{\rho/(\rho+y)\}\log(C/y)-1$ is monotone decreasing, $u(y)$ attains a unique maximum at a $y^{\star}$ that satisfies $u^{\prime}(y^{\star})=0$. This implies that

[TABLE]

for all $0\leq y\leq C$. The exact value of $y^{\star}$ is difficult to determine exactly, by assumption $0<\rho\leq 1$, so for every $n\geq 2\pi\sigma_{1}\sigma_{2}e^{1/e}$, $\rho+\log C\geq 0$, $\log\log C\geq-1$, and

[TABLE]

Since $u^{\prime}(y^{\star})=0$, it follows that $y^{\star}\leq\log C$ and $u(y)\leq(2\sigma_{1}^{2}\log C)/\rho$, or in other words,

[TABLE]

Furthermore, by Jensen’s inequality

[TABLE]

Therefore $|f(x_{1},x_{2},\mu;\mathbf{t})|\leq F_{n}(x_{1},\mu)$ for the function

[TABLE]

The next step to truncate $f$ (16) in expression (17), which will be useful when applying a maximal inequality later in the proof. Define

[TABLE]

By the triangle inequality and the construction of the function $F_{n}$ (19),

[TABLE]

For each $i$, $E\{F_{n}(X_{i1},\theta_{i1})I(|X_{i1}-\theta_{i1}|>C_{n})\}$ equals

[TABLE]

By Mill’s inequality,

[TABLE]

and furthermore,

[TABLE]

This implies that $n^{-1}\sum_{i}E\{F_{n}(X_{i1},\theta_{i1})I(|X_{i1}-\theta_{i1}|>C_{n})=O(C_{n}^{-1})\rightarrow 0$, and it remains to show that

[TABLE]

D.4 Approximation and covering number

For $Cn^{1/4-\eta}$ as in Assumption 1 and some integer $K_{n}$ that can grow with $n$, define

[TABLE]

Next define functions

[TABLE]

indexed by $\mathbf{u}\in\mathcal{U}(K_{n})$ and $\mathbf{w}\in\mathcal{W}(K_{n})$. This section constructs discrete subsets $\mathbb{U}(K_{n})\subset\mathcal{U}(K_{n})$ and $\mathbb{W}(K_{n})\subset\mathcal{W}(K_{n})$ such that for all $\mathbf{t}\in\mathcal{T}$ (10), there exists some $(\mathbf{u},\mathbf{w})\in\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$ such that conditional on the $(X_{i1},X_{i2})$, for any $\epsilon>0$,

[TABLE]

for $n$ sufficiently large, with $C_{n}$ defined in (20).

The $\epsilon$-covering number $\mathcal{N}(\mathcal{T})$ of the set of functions $f(x_{1},x_{2},\mu;\mathbf{t})$ indexed by $\mathbf{t}\in\mathcal{T}$ is defined to be the cardinality of the smallest set $\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$ such that (24) holds. This section establishes an upper bound for this covering number, which will be a function of the $(X_{i1},X_{i2})$ and $K_{n}$. The proof makes use of the following lemmas.

Lemma 1

For any $\gamma_{n}\leq e^{-e}$ and any integer $q>0$, define

[TABLE]

Then for each $\mathbf{t}\in\mathcal{T}$ (10) there exist $K_{n}$ support points $\boldsymbol{\mu}\in\mathcal{U}(K_{n})$ and $\boldsymbol{\omega}\in\mathcal{W}(K_{n})$, with $\mathcal{U}(K_{n})$ and $\mathcal{W}(K_{n})$ defined in (22), such that for each $(X_{i1},X_{i2})$,

[TABLE]

Proof. By Taylor expansion and because $L_{n}!\geq L_{n}^{L_{n}}e^{-L_{n}}$,

[TABLE]

for $d=1,2$. By assumption, $\mathbf{t}\in\mathcal{T}$, so

[TABLE]

Then choosing $L_{n}$ as in (25) guarantees that

[TABLE]

where the last inequality follows because since $\gamma_{n}<e^{-e}$ by assumption, $\log\gamma_{n}<0$ while $\log\log\gamma_{n}^{-1}\geq 1$. Therefore

[TABLE]

which also implies that

[TABLE]

This imply that

[TABLE]

Next, by binomial expansion,

[TABLE]

Lemma A.1 of Ghosal and Van Der Vaart (2001) can be used to show that for $K_{n}$ as in (25), there exist $\boldsymbol{\mu}\in\mathcal{U}(K_{n})$ and $\boldsymbol{\omega}\in\mathcal{W}(K_{n})$ (22) such that $\sum_{j=1}^{n}t_{j1}^{m}t_{j2}^{m^{\prime}}n=\sum_{k=1}^{K_{n}}\mu_{k1}^{m}\mu_{k2}^{m^{\prime}}\omega_{k}$ for every $0\leq m\leq 2\ell+q$ and $0\leq m^{\prime}\leq 2\ell^{\prime}$. Therefore

[TABLE]

Finally, by the same reasoning that led to (26),

[TABLE]

Combining (26), (27), and (28) concludes the proof. $\square$

Next, for any $\gamma_{n}>0$, define the discrete sets

[TABLE]

The set $\mathbb{U}(K_{n})$ can be constructed by specifying a grid of equally-spaced points on the interval $[-Cn^{1/4-\eta},Cn^{1/4-\eta}]$, where any two neighboring points are separated by at most a distance $\gamma_{n}$, for each of the $2K_{n}$ dimensions of the vector $\boldsymbol{\mu}\in\mathcal{U}(K_{n})$. Therefore the cardinality of the smallest possible $\mathbb{U}(K_{n})$ satisfies

[TABLE]

Arguments from Jiang et al. (2009) and Zhang (2009) imply that the cardinality of the smallest possible $\mathbb{W}(K_{n})$ satifies

[TABLE]

Lemma 2

For any $(\boldsymbol{\mu},\boldsymbol{\omega})\in\mathcal{U}(K_{n})\times\mathcal{W}(K_{n})$ (22) and any integer $q>0$, there exists $(\mathbf{u},\mathbf{w})\in\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$ (29) and a constant $D_{q}$ such that for each $(X_{i1},X_{i2})$,

[TABLE]

Proof. For any $(\boldsymbol{\mu},\boldsymbol{\omega})\in\mathcal{U}(K_{n})\times\mathcal{W}(K_{n})$, by construction of $\mathbb{W}(K_{n})$ (29) there exists a $\mathbf{w}\in\mathbb{W}(K_{n})$ such that

[TABLE]

and for any $\mathbf{u}\in\mathcal{U}(K_{n})$,

[TABLE]

Define the function

[TABLE]

for $q\geq 0$. Now if $\mathbf{u}$ is also in $\mathbb{U}(K_{n})$ (29), then by construction

[TABLE]

Since

[TABLE]

the first derivative of $\{h_{q}^{\prime}(x)\}^{2}$ is given by

[TABLE]

Therefore $\{h_{q}^{\prime}(x)\}^{2}$ has at most seven stationary points, is increasing for $x$ smaller than the smallest negative stationary point, and is decreasing for $x$ larger than the largest positive stationary point. Therefore the maximum value of $|h_{q}^{\prime}(x)|$ can only occur at one of the seven stationary points, at which $|h_{q}^{\prime}(x)|$ is finite. Therefore there is some constant $D_{q}$ such that

[TABLE]

Combining (32) and (33) concludes the proof. $\square$

For $K_{n}$ (25) large enough, the result (24), that $f$ (16) can be in a sense approximated by $g$ (23), can now be established using Lemmas 1 and 2. For each $(X_{i1},X_{i2})$,

[TABLE]

The terms (34)–(37) can be bounded using (18) and Lemmas 1 and 2. For $L_{n}$ as defined in (25) and

[TABLE]

there exist $(\mathbf{u},\mathbf{w})\in\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$ (29) such that for any $\gamma_{n}\leq e^{-e}$,

[TABLE]

so by Lemmas 1 and 2,

[TABLE]

Similarly,

[TABLE]

Next, using Jensen’s inequality and the bound (18),

[TABLE]

By similar reasoning as in the bounding of (34), it follows that

[TABLE]

and that

[TABLE]

Finally, choose

[TABLE]

so that as long as $n\log^{2}n>e^{e}$, $\gamma_{n}\leq n^{-e}$ and Lemmas 1 and 2 can be applied to obtain the bounds on (34)–(37). These bounds then imply that for $C_{n}$ as in (20), $M_{n}$ and $L_{n}$ as in (25), $K_{n}$ as in (38), and $\mathbb{U}(K_{n})$ and $\mathbb{W}(K_{n})$ constructed as in (29), for any $\mathbf{t}\in\mathcal{T}$ (10) there exists some $(\mathbf{u},\mathbf{w})\in\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$ such that conditional on the $(X_{i1},X_{i2})$,

[TABLE]

With $\gamma_{n}$ as in (39), for any $\epsilon>0$ there exists a sufficiently large $n$ such that the approximation result (24) holds. The cardinality bounds (30) and (31) imply that the $\epsilon$-covering number $\mathcal{N}(\mathcal{T})$ of the $f(x_{1},x_{2},\mu;\mathbf{t})$ indexed by $\mathbf{t}\in\mathcal{T}$ obeys

[TABLE]

where $C$ and $\eta$ are defined in Assumption 1 and $C_{1}$ and $C_{2}$ are positive constants. Since $K_{n}=4L_{n}^{2}+1$ from (38) and $L_{n}=(1+M_{n}^{2}e/2)\log(1/\gamma_{n})>1$ from (25), $K_{n}\leq 5L_{n}^{2}$ and there is a constant $C_{2}$ such that

[TABLE]

D.5 Maximal inequality

This section concludes the proof of Theorem 2 by showing (21). First, the supremum over the infinite set $\mathcal{T}$ can be converted to a maximum over the discrete set $\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$ constructed in (29). For any $\mathbf{t}\in\mathcal{T}$, using the approximation result (24) there exists a $(\mathbf{u},\mathbf{w})$ inside a suitably constructed set $\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$ (29) such that for any $\epsilon>0$, there exists a sufficiently large $n$ such that

[TABLE]

where $\varepsilon_{i}$ are Rademacher variables and $C_{n}=(2\sigma_{1}^{2}\log\log n)^{1/2}$ from (20). Since $\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$ is discrete, it follows that

[TABLE]

To bound the expectation of the maximum over $\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$, define the convex function $\psi(x)=\exp(x^{2})-1$. Using Lemma 2.2.2 of van der Vaart and Wellner (1996), it can be shown that there is a constant $A$ such that

[TABLE]

where $E_{\varepsilon}$ is the expectation over the $\varepsilon_{i}$ conditional on the $(X_{i1},X_{i2})$ and $\|X\|_{\psi}$ is the Orlicz norm of the random variable $X$ with respect to the function $\psi(x)$. It is straightforward to show that $|g(x_{1},x_{2},\mu;\mathbf{u},\mathbf{w})|\leq F_{n}(x_{1},\mu)$ for the function $F_{n}$ defined in (19), so by Hoeffding’s inequality and Lemma 2.2.7 of van der Vaart and Wellner (1996),

[TABLE]

for any $(\mathbf{u},\mathbf{w})\in\mathbb{U}(K_{n})\times\mathbb{W}(K_{n})$. Therefore, using Jensen’s inequality,

[TABLE]

To show (21), it remains to show that $E\log\mathcal{N}(\mathcal{T})=o(n/\log^{2}n)$, where the expectation is now over the $(X_{i1},X_{i2})$. From (40) and the inequality $\log(x)\leq mx^{1/m}$ for any $x>0$ and integer $m>0$, for $n$ sufficiently large there is some constant $A_{2}$ such that

[TABLE]

By the definition of $M_{n}$ (25) and the bounds on $\theta_{id}$ from Assumption 1,

[TABLE]

Therefore letting $Z_{i}$, $i=1,\ldots,2n$ denote independent standard normal random variables,

[TABLE]

Now define the convex function $\psi(x)=x^{p/(4+10/m)}$ for $p\geq 4+10/m$. Then by Lemma 2.2.2 of van der Vaart and Wellner (1996), there is some constant $A_{3}$ such that

[TABLE]

where $\|X\|_{\psi}$ is the Orlicz norm of $X$ with respect to the function $\psi(x)$. But

[TABLE]

which is a constant that does not depend on $n$, so choosing $p$ such that $p\geq 5$ and $1/p\leq 1/4-\eta$ gives

[TABLE]

Then

[TABLE]

Choosing $m$ large enough such that $4\eta-32/(4m)+10\eta/m>0$ implies that $E\log\mathcal{N}(\mathcal{T})=o(n/log^{2}n)$, concluding the proof. $\square$

Appendix E Proof of Theorem 3

Let $\boldsymbol{\theta}=(\theta_{11},\ldots,\theta_{n1},\theta_{12},\ldots,\theta_{n2})$. Then when $\mathbf{t}=\boldsymbol{\theta}$, the estimator $\boldsymbol{\delta}^{\mathbf{t}}_{\rho}$ (7) equals the oracle regularized estimator $\boldsymbol{\delta}^{\star}_{\rho}$. Therefore with $\textsc{sure}(\mathbf{t})$ defined as in (8),

[TABLE]

where the second line follows because $E\textsc{sure}(\hat{\mathbf{t}})\leq E\textsc{sure}(\theta)$, since $\hat{\mathbf{t}}$ (11) is defined to minimize $\textsc{sure}(\mathbf{t})$, and the third line follows by Theorem 1. Finally, since $\hat{\mathbf{t}}\in\mathcal{T}$ by construction, by Theorem 2,

[TABLE]