Inference for high-dimensional exchangeable arrays

Inference for high-dimensional exchangeable arrays

\address[H. D. Chiang]{ Department of Economics, University of Wisconsin-Madison\\ William H. Sewell Social Science Building, 1180 Observatory Drive, Madison, WI 53706, USA.} \email{[email removed]}

\address[K. Kato]{ Department of Statistics and Data Science, Cornell University \\ 1194 Comstock Hall, Ithaca, NY 14853, USA.} \email{[email removed]}

\address[Y. Sasaki]{Department of Economics, Vanderbilt University\\ VU Station B \#351819, 2301 Vanderbilt Place, Nashville, TN 37235-1819, USA.} \email{[email removed]}

Introduction

Many recent statistical problems involve non-independent observations indexed by multiple interlocking sets of entities. Examples include dyadic/polyadic networks, bipartite networks, and multiway clustering. When the sets of entities that form each of these indices are different, as is the case with market-product data and book-reader data, a natural stochastic framework is separate exchangeability MacKinnonNielsenWebb2019. Separately exchangeable arrays include row-column exchangeable models Mccullagh2000, additive cross random effect models Owen2007,OwenEckles2012, and multiway clustering CGM2011. Meanwhile, when all indices belong to a common set of entities, as is the case with friendship network data, the underlying structure is well-captured by joint exchangeability BickelChen2009. Joint exchangeability covers nonparametric random graph models of BickelChen2009 for dyadic networks, which contain widely used models in the statistical network analysis literature such as stochastic block models.

Analysis of these types of data requires accounting for the underlying complex dependence structures induced by these exchangeability notions. Thus, developing valid inference methods for exchangeable arrays is challenging. The literature has witnessed some research on statistical inference that focuses on exchangeable arrays with low or fixed dimensions. For modern statistical learning methods, it is crucial to allow the dimension of data to increase with sample size. However, the existing literature has been silent about statistical inference for such high-dimensional exchangeable arrays.

This paper is concerned with the problem of inference for separately or jointly exchangeable high-dimensional arrays. We develop new high-dimensional central limit theorems (CLTs) over the rectangles for the sample mean under both exchangeability notions. Building on the high-dimensional CLTs, we propose new multiplier bootstrap methods tailored to separate and jointly exchangeable arrays and derive their nonasymptotic error bounds. Such nonasymptotic results can be translated into asymptotic results that hold uniformly over a large set of distributions, which is crucial in many high-dimensional statistical applications.

To derive these theoretical results, we develop several new technical tools, which are of independent interest and would be useful for other analyses of exchangeable arrays. Specifically, we develop novel Hoeffding-type decompositions for both separately and jointly exchangeable arrays and establish novel maximal inequalities for Hoeffding-type projections in both cases. Such maximal inequalities lead to sharp rates for degenerate components in Hoeffding-type decompositions in both cases and play a crucial role in establishing the high-dimensional CLTs and the validity of the bootstrap methods. The proofs of these technical results are highly nontrivial. For example, the proof of the symmetrization inequality for exchangeable arrays involves a careful induction argument (see Lemma (ref) in the Appendix) combined with a repeated conditioning argument. {\color{black}Furthermore, the proof of the maximal inequality for jointly exchangeable arrays involves a delicate conditioning argument combined with the decoupling inequalities for $U$-statistics with index-dependent kernels delaPenaGine1999.}

We illustrate applications of the bootstrap methods to a couple of concrete statistical problems. Specifically, 1) we develop a method to construct simultaneous or uniform confidence bands for density functions with jointly exchangeable dyadic arrays, and 2) we develop a method to choose a penalty level for $\ell_1$-penalized regression (Lasso) and establish error bounds for the Lasso with separately exchangeable arrays. These applications are also new in the literature.

We conduct extensive simulation studies, which demonstrate precise uniform coverage across various designs and under both notions of exchangeability, thereby supporting our theoretical results. Finally, we apply our bootstrap method to international trade network data to draw uniform confidence bands for trade flow volumes in 1990, 1995, 2000, and 2005. The results indicate that there have been increasing numbers of bilateral trading pairs with high flow volumes as time progresses.

Relation to the literature

There is now a large literature on high-dimensional CLTs and bootstraps with the “$p \gg n$” regime; see CCK2013AoS,CCK2014AoS,CCK2015PTRF,CCK2016SPA,CCK2017AoP, DengZhang2017, CCKK2019, kuchibotla2020, and fang2020 for the independent case, Chen2018, ChenKato2019b,ChenKato2019 for $U$-statistics and processes, and ZhangWu2017, ZhangCheng2018, CCK2019RES, Koike2019 for time series dependence. However, none of the above references covers extensions to exchangeable arrays. The present paper builds on and contributes to this literature by developing high-dimensional CLTs and bootstrap methods for exchangeable arrays.

Early applications of exchangeable arrays in statistics include arnold1979linear, bowman1995saturated, and Andrews2005, to name a few. For reviews, see, e.g. goldenberg2010survey, orbanz2014bayesian, and kuchibhotla2020exchangeability. Analysis of exchangeable random graphs has been an active research area in the recent statistics literature; see, e.g., DiaconisJanson2008, BickelChenLevina2011, lloyd2012random, choi2014co, caron2017sparse, choi2017co, zhang2017estimating, crane2018edge. Limit theorems for jointly exchangeable arrays (in the fixed dimensional case) date back to Silverman1976 and EaglesonWeber1978. FafchampsGubert2007 and CGM2011 derive standard error formulas for jointly exchangeable dyadic arrays and separately exchangeable arrays, respectively; see also CameronMiller2014, CM2015, AronowSamiiAssenova2015, and Tabord-Meehan2019 for further development. Menzel2017 studies inference for separately exchangeable arrays, covering both degenerate and non-degenerate cases. DDG2019 develop functional limit theorems for Donsker classes under separate and joint exchangeability. To the best of our knowledge, however, no existing work in this literature permits high-dimensional inference. We note that DDG2019 develop symmetrization inequalities different from ours. Specifically, symmetrization inequalities developed in DDG2019 are applied to the whole empirical process and do not lead to correct orders for degenerate components in Hoeffding-type decompositions (indeed, DDG2019 do not derive Hoeffding-type decompositions), thereby not powerful enough to derive our results; see Remarks (ref) and (ref) in the Appendix for details.

Methodologically, this paper is also related to the recent literature on high-dimensional $U$-statistics, such as Chen2018, ChenKato2019,ChenKato2019b, among others. Under suitable assumptions, the data of our interest can be written as $U$-statistic-like latent structure (in distribution) via the Aldous-Hoover-Kallenberg representation Aldous1981,Hoover1979,Kallenberg2006, i.e. the data can be written as a kernel function of some latent independent random variables. However, unlike in $U$-statistics, neither the kernel nor the latent independent random variables is known to us. In addition, we need to cope with the existence of extra higher-order shocks in the latent structure. Both aspects present extra challenges.

Regarding our bootstraps, Mccullagh2000 shows that no resampling scheme for the raw data is consistent for variance of a sample mean under separate exchangeability. A Pigeonhole bootstrap is subsequently proposed by Owen2007 and its different variants are further investigated in OwenEckles2012, Menzel2017 and DDG2019. Whether the pigeonhole bootstrap works for increasing or high-dimensional test statistics remains unknown to us. We therefore develop a novel bootstrap method in this paper which we argue works for high-dimensional data.

Notations and organization

Let $\mathbb{N}$ denote the set of positive integers. We use $\left\| \cdot \right\|, \left\| \cdot \right\|_{0}, \left\| \cdot \right\|_{1}$, and $\left\| \cdot \right\|_{\infty}$ to denote the Euclidean, $\ell_{0}$, $\ell_1$, and $\ell^{\infty}$-norms for vectors, respectively (precisely, $\left\| \cdot \right\|_{0}$ is not a norm but a seminorm). For two real vectors $\bm{a}= (a_{1},\dots,a_{p})^{T}$ and $\bm{b} = (b_{1},\dots,b_{p})^{T}$, the notation $\bm{a} \le \bm{b}$ means that $a_{j} \le b_{j}$ for all $1 \le j \le p$. Let $\operatorname{supp} (\bm{a})$ denote the support of $\bm{a} = (a,\dots,a_p)^{T}$, i.e., $\operatorname{supp}(\bm{a}) = \{ j : a_j \ne 0\}$. We denote by $\odot$ the Hadamard (element-wise) product, i.e., for ${\bm{i}} = (i_1,\dots,i_K)$ and $\bm{j} = (j_1,\dots,j_K)$, ${\bm{i}} \odot \bm{j}= (i_1 j_1,\dots,i_K j_K)$. For any $a,b \in \mathbb{R}$, let $a \vee b = \max \{ a,b \}$. For $0 < \beta < \infty$, let $\psi_{\beta}$ be the function on $[0,\infty)$ defined by $\psi_{\beta} (x) = e^{x^{\beta}}-1$. Let $\| \cdot \|_{\psi_{\beta}}$ denote the associated Orlicz norm, i.e., $\| \xi \|_{\psi_\beta}=\inf \{ C>0: \mathbb{E}[ \psi_{\beta}( | \xi | /C)] \leq 1\}$ for a real-valued random variable $\xi$. “Constants” refer to nonstochastic and finite positive numbers.

The rest of the paper is organized as follows . In Section (ref), we develop a high-dimensionl CLT (over the rectangles) and a bootstrap method for separately exchangeable arrays. In Section (ref), we develop analogous results to jointly exchangeable arrays. We illustrate {\color{black}two applications} in Section (ref), present simulation results in Section (ref), and demonstrate an empirical application in Section (ref). We defer all the technical proofs to the Appendix.

Separately exchangeable arrays

In this section, we consider separately exchangeable arrays that correspond to multiway clustered data. Pick any $K \in \mathbb{N}$. With ${\bm{i}} = (i_{1},\dots, i_{K}) \in \mathbb{N}^{K}$, we consider a $K$-array $( {\bm{X}}_{{\bm{i}}})_{{\bm{i}} \in \mathbb{N}^{K}}$ consisting of random vectors in $\mathbb{R}^{p}$ {\color{black}with $p \ge 2$}. We denote by $X_{{\bm{i}}}^{j}$ the $j$-th coordinate of ${\bm{X}}_{{\bm{i}}}$: ${\bm{X}}_{{\bm{i}}} = (X_{{\bm{i}}}^{1},\dots,X_{{\bm{i}}}^{p})^{T}$. We say that the array $({\bm{X}}_{{\bm{i}}})_{{\bm{i}} \in \mathbb{N}^K}$ is separately exchangeable if the following condition is satisfied Kallenberg2006.

See Appendix (ref) in the supplementary material for more details, discussions, and examples. From the Aldous-Hoover-Kallenberg representation Kallenberg2006, any separately exchangeable array $({\bm{X}}_{{\bm{i}}})_{{\bm{i}} \in \mathbb{N}^{K}}$ is generated by the structure \[ {\bm{X}}_{{\bm{i}}} = \mathfrak{f} ((U_{{\bm{i}} \odot {\bm{e}}})_{{\bm{e}} \in \{ 0,1 \}^{K}}), \ {\bm{i}} \in \mathbb{N}^{K}, \quad \{ U_{{\bm{i}} \odot {\bm{e}}} : {\bm{i}} \in \mathbb{N}^{K}, {\bm{e}} \in \{ 0,1 \}^{K} \} \stackrel{i.i.d.}{\sim} U[0,1] \] for some Borel measurable map $\mathfrak{f}:[0,1]^{2^{K}} \to \mathbb{R}^{p}$.

The latent variable $U_{\bm{0}}$ appears commonly in all ${\bm{X}}_{{\bm{i}}}$'s. In the present paper, as in Andrews2005 and Menzel2017, we consider inference conditional on $U_{\bm{0}}$ and treat it as fixed. In the rest of Section (ref), we will assume (without further mentioning) that the array $({\bm{X}}_{{\bm{i}}})_{{\bm{i}} \in \mathbb{N}^K}$ has mean zero (conditional on $U_{\bm{0}}$) and is generated by the structure

where $\mathfrak{g}$ is now a map from $[0,1]^{2^{K}-1}$ into $\mathbb{R}^{p}$.

Suppose that we observe $\{ {\bm{X}}_{{\bm{i}}} : {\bm{i}} \in [\bm N] \}$ with ${\bm{N}} = (N_1,\dots,N_K)$ and $[\bm N] = \prod_{k=1}^{K} \{ 1,\dots,N_{k} \}$. We are interested in approximating the distribution of the sample mean \[ {\bm{S}}_{{\bm{N}}} =\frac{1}{\prod_{k=1}^{K}N_{k}}\sum_{{\bm{i}} \in [{\bm{N}}]} {\bm{X}}_{{\bm{i}}} \] in the high-dimensional setting where the dimension $p$ is allowed to entail $p \gg \min \{N_1,\dots,N_K\}$.

For later convenience, we fix some additional notations. Let $n = \min_{1 \le k \le K} N_{k}$ and $\overline{N} = \max_{1 \le k \le K}N_{k}$ denote the minimum and maximum cluster sizes, respectively. For $1 \le k \le K$, denote by ${\mathcal{E}}_{k} = \{ {\bm{e}}= (e_1,\dots,e_K) \in \{ 0,1 \}^{K}: \sum_{k=1}^{K} e_k = k \}$ the set of vectors in $\{ 0,1 \}^{K}$ whose support has cardinality $k$. Let ${\bm{e}}_{k} \in \mathbb{R}^{K}$ denote the vector such that the $k$-th coordinate of ${\bm{e}}_{k}$ is $1$ and the other coordinates are $0$. For a given ${\bm{e}} \in \{ 0,1 \}^{K}$, define \[ I_{{\bm{e}}} ([{\bm{N}}]) = \{ {\bm{i}} \odot {\bm{e}} : {\bm{i}} \in [{\bm{N}}] \} \subset \mathbb{N}_{0}^{K} \quad \text{with} \ \mathbb{N}_0 = \mathbb{N} \cup \{ 0 \}. \] The following decomposition of the sample mean ${\bm{S}}_{{\bm{N}}}$ will play a fundamental role in our analysis, which is reminiscent of the Hoeffding decomposition for $U$-statistics Lee1990, delaPenaGine1999.

The proof of this lemma can be found in Appendix (ref).

The leading term in the decomposition ((ref)) is \[ \sum_{{\bm{e}} \in \mathcal{E}_{1}} \frac{1}{\prod_{k' \in \operatorname{supp} ({\bm{e}})} N_{k'}} \sum_{{\bm{i}} \in I_{{\bm{e}}}([{\bm{N}}])} \hat{{\bm{X}}}_{{\bm{i}}} =\sum_{k=1}^{K}N_{k}^{-1} \sum_{i_{k}=1}^{N_{k}} \mathbb{E}[ {\bm{X}}_{{\bm{i}}} \mid U_{(0,\dots,0,i_{k},0,\dots,0)} ], \] which we call the H\'{a}jek projection of ${\bm{S}}_{{\bm{N}}}$. With this in mind, define ${\bm{W}}_{k,i_{k}} = \mathbb{E}[ {\bm{X}}_{{\bm{i}}} \mid U_{(0,\dots,0,i_{k},0,\dots,0)} ]$ for $k=1,\dots,K$.

High-dimensional CLT for separately exchangeable arrays

We first establish a high-dimensional CLT for ${\bm{S}}_{{\bm{N}}}$ over the class of rectangles, $\mathcal{R}= \{ \prod_{j=1}^{p} [a_{j},b_{j}] : - \infty \le a_{j} \le b_{j} \le \infty, \ 1 \le j \le p \}$. This high-dimensional CLT will be a building block for establishing the validity of the multiplier bootstrap (cf. Section (ref)).

We start with discussing regularity conditions. Denote by $\bm{1} = (1,\dots,1)$ the vector of ones. Let $D_{{\bm{N}}} \ge 1$ be a given constant that may depend on the cluster sizes ${\bm{N}}$ {\color{black}(and $p$; when we consider asymptotics we have in mind that $p$ is a function of ${\bm{N}}$ or $n$ so we omit the dependence of $D_{{\bm{N}}}$ on $p$)}, and let $\underline{\sigma} > 0$ be another given constant independent of the cluster sizes ${\bm{N}}$. We will assume either of the following moment conditions.

We will also assume the following condition.

Condition ((ref)) requires that each coordinate of ${\bm{X}}_{\bm 1}$ is sub-exponential. By Jensen's inequality, Condition ((ref)) implies that $\max_{1 \le j \le p;1 \le k \le K}\| W_{k,1}^j \|_{\psi_1} \le D_{{\bm{N}}}$. Condition ((ref)) is an alternative moment condition on ${\bm{X}}_{\bm{1}}$. Condition ((ref)) is satisfied for example under the following situation: Suppose that ${\bm{X}}_{{\bm{i}}}$ is given by ${\bm{X}}_{{\bm{i}}} = \varepsilon_{{\bm{i}}} \bm{Z}_{{\bm{i}}}$ where $\varepsilon_{{\bm{i}}}$ is a scalar “error” variable while $\bm{Z}$ is a vector of “covariates”. If each coordinate of $\bm{Z}_{{\bm{i}}}$ is bounded by a constant $\overline{D}$ and $\varepsilon_{{\bm{i}}}$ has finite $q$-th moment, then $\mathbb{E}[\| {\bm{X}}_{{\bm{i}}} \|_{\infty}^{q}] \le \overline{D}^{q} \mathbb{E}[|\varepsilon_{{\bm{i}}}|^{q}]$. {\color{black}Also Condition ((ref)) is satisfied if, in the discretized empirical process application (cf. Example (ref)), the function class possesses an envelope function with finite $q$-th moment.} Again, by Jensen's inequality, Condition ((ref)) implies that $\max_{1 \le k \le K}\mathbb{E}[\| {\bm{W}}_{k,1} \|_{\infty}^{q}] \le D_{{\bm{N}}}^{q}$. {\color{black}The restriction $q > 4$ is needed to guarantee that Condition ((ref)) appearing in Theorem (ref) to be non-void.}

Condition ((ref)) requires the maximum of third (respectively, fourth) moment across coordinates to be increasing at speed no faster than the first (respectively, second) power of $D_{{\bm{N}}}$. By Jensen's inequality, the first part of Condition ((ref)) is satisfied if $\max_{1 \le j \le p} \mathbb{E}[|X_{\bm{1}}^{j}|^{2+\kappa}] \le D_{{\bm{N}}}^{\kappa}$ for $\kappa=1,2$. The second part of Condition ((ref)) guarantees that the H\'ajek projection is nondegenerate.

Let $\gamma = N(\bm{0},\Sigma)$ with $\Sigma = \sum_{k=1}^K (n/N_k) \Sigma_{W_k}$ and $\Sigma_{W_k} = \mathbb{E}[{\bm{W}}_{k,1}{\bm{W}}_{k,1}^T]$ for $k=1,\dots,K$.

Multiplier bootstrap for separately exchangeable arrays

Let $\{ \xi_{1,i_{1}} \}_{i_1=1}^{N_{1}}, \dots, \{ \xi_{K,i_{K}} \}_{i_K=1}^{N_K}$ be independent $N(0,1)$ random variables independent of the data. Ideally, we want to make use of the bootstrap statistic $ {\color{black}\sum_{k=1}^K} N_k^{-1}\sum_{i_k=1}^{N_k}\xi_{k,i_k}({\bm{W}}_{k,i_k}-{\bm{S}}_{{\bm{N}}}). $ However, this bootstrap is infeasible as ${\bm{W}}_{k,i_k} = \mathbb{E}[{\bm{X}}_{{\bm{i}}} \mid U_{(0,\dots,i_{k},\dots,0)}]$ are unknown to us. Estimation of ${\bm{W}}_{k,i_k}$ is nontrivial as $U_{(0,\dots,i_{k},\dots,0)}$ is a latent variable. We propose to estimate each ${\bm{W}}_{k,i_{k}}$ by \[ \overline{{\bm{X}}}_{k,i_{k}} = \frac{1}{\prod_{k' \ne k} N_{k'}} \sum_{i_{1},\dots,i_{k-1},i_{k+1},\dots,i_{K}} {\bm{X}}_{{\bm{i}}}, \ i_{k} = 1,\dots,N_{k}; k=1,\dots,K, \] i.e., the sample mean taken over all indices but $i_k$. Then, we apply the multiplier bootstrap to $\overline{{\bm{X}}}_{k,i_k}$ in place of ${\bm{W}}_{k,i_k}$ \[ {\bm{S}}_{{\bm{N}}}^{MB} = \sum_{k=1}^{K} N_{k}^{-1} \sum_{i_{k}=1}^{N_{k}} \xi_{k,i_{k}} (\overline{{\bm{X}}}_{k,i_k} - {\bm{S}}_{{\bm{N}}}). \] To the best of our knowledge, this multiplier bootstrap for separately exchangeable arrays is new in the literature. We will formally study the validity of this multiplier bootstrap for high-dimensional separately exchangeable arrays with $p \gg n$.

We are now in position to establish the validity of the proposed multiplier bootstrap for separately exchangeable arrays. Let $\mathbb{P}_{|{\bm{X}}_{[{\bm{N}}]}}$ denote the law conditional on the data ${\bm{X}}_{[{\bm{N}}]} = ({\bm{X}}_{{\bm{i}}})_{{\bm{i}}\in[{\bm{N}}]}$ and $\overline \sigma=\max_{1\le j \le p; 1 \le k \le K} \sqrt{\mathbb{E}[|W_{k,1}^j|^2]}$.

Jointly exchangeable arrays

In this section, we consider another class of exchangeable arrays, namely, jointly exchangeable arrays. The notations in the current section are independent from those in Section (ref) unless otherwise noted. Joint exchangeability induces a more complex dependence structure on arrays than separate exchangeability, but still we are able to develop analogous results to the preceding section for jointly exchangeable arrays as well. It should be noted, however, that we do require a different bootstrap and technical tools (cf. Appendix (ref)) to accommodate a specific dependence structure induced from joint exchangeability.

Pick any $K \in \mathbb{N}$. For a given positive integer $n \ge K$, let $I_{n,K} = \{ (i_1,\dots,i_K) : 1 \le i_1, \dots, i_K \le n \ \text{and $i_1,\dots,i_K$ are distinct} \}$. Also let $I_{\infty,K} = \bigcup_{n=K}^{\infty} I_{n,K}$. For any ${\bm{i}} = (i_1,\dots,i_K) \in \mathbb{N}^K$, let $\{ {\bm{i}} \}^{+}$ denote the set of distinct nonzero elements of $(i_1,\dots,i_K)$.

In this section, we consider a $K$-array $( {\bm{X}}_{{\bm{i}}})_{{\bm{i}} \in I_{\infty,K}}$ consisting of random vectors in $\mathbb{R}^{p}$ {\color{black}with $p \ge 2$}. We say that the array $({\bm{X}}_{{\bm{i}}})_{{\bm{i}} \in I_{\infty,K}}$ is jointly exchangeable if the following condition is satisfied Kallenberg2006.

See Appendix (ref) in the supplementary material for more details, discussions, and examples. From the Aldous-Hoover-Kallenberg representation Kallenberg2006, any jointly exchangeable array $({\bm{X}}_{{\bm{i}}})_{{\bm{i}} \in I_{\infty,K}}$ is generated by the structure \[ {\bm{X}}_{{\bm{i}}} = \mathfrak{f} ( (U_{\{{\bm{i}} \odot {\bm{e}}\}^{+}})_{{\bm{e}} \in \{ 0,1 \}^{K}}), \ {\bm{i}} \in I_{\infty,K}, \quad \{ U_{\{{\bm{i}} \odot {\bm{e}}\}^{+}} : {\bm{i}} \in I_{\infty,K}, {\bm{e}} \in \{ 0,1 \}^{K} \} \stackrel{i.i.d.}{\sim} U[0,1] \] for some Borel measurable map $\mathfrak{f}: [0,1]^{2^{K}} \to \mathbb{R}^{p}$. Here the coordinates of the vector $(U_{\{{\bm{i}} \odot {\bm{e}}\}^{+}})_{{\bm{e}} \in \{ 0,1 \}^{K}}$ are understood to be properly ordered, so that, e.g., when $K=2$, ${\bm{X}}_{(i_1,i_2)} = \mathfrak{f} (U_{\varnothing},U_{ i_1}, U_{ i_2}, U_{\{ i_1,i_2 \}})$ and ${\bm{X}}_{(i_2,i_1)} = \mathfrak{f} (U_{\varnothing}, U_{i_2}, U_{i_1}, U_{\{ i_1,i_2 \}})$ differ (although they have the identical distribution).

As in the separately exchangeable case, we consider inference conditional on $U_{\varnothing}$, and in what follows, we will assume that the array $({\bm{X}}_{{\bm{i}}})_{{\bm{i}} \in I_{\infty,K}}$ has mean zero (conditional on $U_{\varnothing}$) and is generated by the structure

where $\mathfrak{g}$ is now a map from $[0,1]^{2^{K}-1}$ into $\mathbb{R}^{p}$.

Suppose that we observe $\{ {\bm{X}}_{{\bm{i}}}:{\bm{i}}\in I_{n,K} \}$ with $n \ge K$ and are interested in distributional approximation of the polyadic sample mean

in the high-dimensional setting where the dimension $p$ is allowed to entail $p \gg n$.

As in Section (ref), define ${\mathcal{E}}_{k} = \{ {\bm{e}}= (e_1,\dots,e_K) \in \{ 0,1 \}^{K}: \sum_{k=1}^{K} e_k = k \}$ for $1 \le k \le K$. The analysis of the jointly exchangeable array relies on the following decomposition

It turns out that the first term on the right-hand side, which we call the the H\'{a}jek projection of ${\bm{S}}_{n}$, is a dominant term. Defining $h_k(u)=\mathbb{E}[{\bm{X}}_{(1,\dots,K)} \mid U_{k}=u]$ for $k=1,\dots,K$, we can simplify the H\'ajek projection into $n^{-1}\sum_{i=1}^n {\bm{W}}_j$ where ${\bm{W}}_j = \sum_{k=1}^K h_k(U_j)$.

High-dimensional CLT for jointly exchangeable arrays

We consider to approximate the distribution of $\sqrt {n} {\bm{S}}_{n}$ by a Gaussian distribution on the set of rectangles $\mathcal{R}$ as defined in Section (ref).

Let $D_n \ge 1$ be a given constant that may depend on $n$, and $\underline{\sigma} > 0$ be another given constant independent of $n$. We will assume either of the following moment conditions.

We will also assume the following condition.

The conditions required here are similar to those in the case of separate exchangeability in Section (ref). The main difference is that Condition ((ref)) is now imposed on ${\bm{W}}_1$.

Let $\gamma_{\Sigma} = N(\bm{0},\Sigma)$ with $\Sigma =\mathbb{E}\left[{\bm{W}}_{1}{\bm{W}}_{1}^{T}\right]$.

Multiplier bootstrap for jointly exchangeable arrays

Let $\{ \xi_j \}_{j=1}^{n}$ be independent $N(0,1)$ random variables independent of the data. Ideally, we want to make use of the multiplier bootstrap statistic $ n^{-1}\sum_{j=1}^n \xi_j( {\bm{W}}_{j}-K{\bm{S}}_n). $ This is infeasible, however, as the projections ${\bm{W}}_j$ are unknown. As an alternative, we replace each ${\bm{W}}_j$ by its estimate \[ \hat{{\bm{W}}}_j =\frac{(n-K)!}{(n-1)!}\sum_{k=1}^K\sum_{{\bm{i}}\in I_{n,K}: i_k=j} {\bm{X}}_{{\bm{i}}}, \] and apply the multiplier bootstrap to $\hat{{\bm{W}}}_{j}$, i.e.,

When $K=2$ (dyadic), this mulitplier bootstrap coincides with the multiplier bootstrap statistic considered in Section 3.2 of DDG2019. However, DDG2019 do not consider the extension to general $K$ arrays, and focus on the empirical process indexed by a Donsker class, which excludes the high-dimensional sample mean. We will study the validity of this multiplier bootstrap for jointly exchangeable arrays.

Let $\mathbb{P}_{|{\bm{X}}_{I_{n,K}}}$ denote the law conditional on the data $({\bm{X}}_{{\bm{i}}})_{{\bm{i}} \in I_{n,K}}$ and $\overline \sigma=\max_{1\le \ell \le p}\sqrt{\mathbb{E}[|W_{1}^\ell|^2]}$.

Applications

In this section, we illustrate a couple of applications of our bootstrap methods. Section (ref) is concerned with construction of confidence bands for densities of flows in dyadic data. Section (ref) is concerned with penalty choice for the Lasso and the performance of the corresponding estimate.

Confidence bands for densities of flows in dyadic data

Researchers are often interested in “the densities of migration across states, trade across nations, liabilities across banks, or minutes of telephone conversation among individuals” Graham2019kernel. Densities of these flow measures use dyadic data. We illustrate an application of our method in Section (ref) to constructing confidence bands for such density functions. We refer the reader to bickel1973, claeskens2003, cck2014density as references on confidence bands for density estimation with i.i.d. data.

Following Graham2019kernel, suppose that we observe the dyadic data $\{Y_{ij} : 1 \le i \ne j \le n \}$ that admits the structure

where $\mathfrak{g}$ is symmetric in the first two arguments and hence $Y_{ij}=Y_{ji}$. We are interested in inference on the density of $Y_{ij}$. However, in certain empirical applications, such as international trade head2014gravity, a proportion of the variable of interest is zero. Hence we assume that $Y_{ij}$ has a probability mass at zero, i.e. $Y_{ij}$ is such that $\mathbb{P}(Y_{ij}\ne 0)=a\in(0,1]$, and $Y_{ij}\sim f$ when $Y_{ij}\ne 0$, where $f$ is a density function on $\mathbb{R}$. Let $b(y)=af(y)$ denote the scaled density. We may estimate $f(\cdot) = b(\cdot)/a$ by $\hat f(\cdot) =\hat b(\cdot)/\hat a$, where $ \hat a={n\choose 2}^{-1}\sum_{1\le i<j\le n} \mathbbm{1}(Y_{ij}\ne 0)$ and $\hat b(y)={n\choose 2}^{-1}\sum_{1\le i<j\le n} K_h(y- Y_{ij})\mathbbm{1}(Y_{ij}\ne 0). $ Here $K:\mathbb{R}\to \mathbb{R}$ is a kernel function (a function that integrates to one), $K_h(\cdot):=h^{-1} K(\cdot/h)$, and $h=h_n\to 0$ is a bandwidth.

We consider to construct simultaneous confidence intervals (bands) for $f$ over the set of design points $y_1,\dots,y_p$, where $p=p_n \to \infty$ is allowed. Define

for $\ell=1,\dots,p$. Then, the multiplier bootstrap statistic is given by

where $\sum_{j \ne i} = \sum_{j \in \{1,\dots,n\} \setminus \{ i \}}$. For a given $\alpha\in (0,1)$, consider the $(1-\alpha)$-simultaneous confidence intervals defined by

where $\tilde\sigma_\ell^2 = n^{-1}\sum_{i=1}^n( \tilde W_i^\ell - 2\tilde S_n^\ell )^2 $, $\tilde\Lambda=\operatorname{diag}(\tilde\sigma_1^2,\dots,\tilde \sigma_p^2)$, $\tilde c(1-\alpha)$ is the conditional $(1-\alpha)$-quantile of $\|\sqrt{n}\tilde{\bm{S}}_n^{MB}\|_\infty$, and $\tilde c^N(1-\alpha)$ is the conditional $(1-\alpha)$-quantile of $\|\sqrt{n}\hat\Lambda^{-1/2}\tilde{\bm{S}}_n^{MB}\|_\infty$. The first method $\mathcal I(1-\alpha)$ is a constant-length confidence band, while the second method $\mathcal I^N(1-\alpha)$ is a variable-length confidence band based on Studentization.

The following proposition establishes asymptotic validity of the confidence bands. We will assume that there exists a conditional density of $Y_{ij}$ given $U_i$ and $Y_{ij} \ne 0$, denoted by $f_{Y_{12} \mid U_1, Y_{12} \ne 0}(y \mid u)$ (more formally, we assume that the conditional distribution of $Y_{ij}$ given $U_i$ is $\mathbb{P} (Y_{ij} \in dy \mid U_i) = \mathbb{P} (Y_{ij} = 0 \mid U_i) \delta_{0}(dy) + \mathbb{P} (Y_{ij} \ne 0 \mid U_i) f_{Y_{12} \mid U_{1},Y_{12} \ne 0} (y \mid U_i) dy$, where $\delta_{0}$ is the Dirac delta at $0$). Let $\overline f_h(y)=\int K_h(y-z)f(z)dz $ and $\overline f_h (y\mid u)=\int K_{h}(y-z) f_{Y_{12}\mid U_1, Y_{12}\ne 0}(z \mid u)dz$ denote the surrogate density and conditional density, respectively. Recall that a kernel $K$ is an $r$-th order kernel for some $r \ge 2$ if $\int y^{t} K(y) dy = 0$ for $t =1,\dots,r-1$ and $\int |y^r K(y)| dy < \infty$. Let $M$, $h_0$, $\sigma_0$, and $a\in(0,1]$ be given positive constants independent of $n$.

Some comments on the proposition are in order.

Penalty choice for Lasso under separate exchangeability

Consider a regression model

where $Y_{\bm{i}}$ is a scalar outcome variable, $\bm{Z}_{\bm{i}}\in \mathbb{R}^d$ is a $d$-dimensional vector of covariates, $f:\mathbb{R}^d \to \mathbb{R} $ is an unknown regression function of interest, and $\varepsilon_{\bm{i}}$ is an error term. We approximate $f$ by a linear combination of technical controls $\bm{X}_{\bm{i}}=P(\bm{Z}_{\bm{i}})$ for some transformation $P:\mathbb{R}^d \to \mathbb{R}^p$, i.e., $f(\bm{Z}_{\bm{i}})=\bm{X}_{\bm{i}}^T\beta_0 + r_{\bm{i}},\, \bm{i}\in [\bm{N}]$, where $r_{\bm{i}}$ is a bias term. The dimension $p$ can be much larger than the cluster sizes ${\bm{N}}$, but we assume that the vector $\beta_0\in \mathbb{R}^p$ is sparse in the sense that $\|\beta_0\|_{0}=s\ll n$ with $n = \min_{1 \le k \le K} N_k$. Suppose that the array $\big ( (Y_{\bm{i}},\bm{Z}_{\bm{i}}^T)^T \big)_{{\bm{i}} \in \mathbb{N}^{K}}$ is separately exchangeable and generated as \[ (Y_{\bm{i}},\bm{Z}_{\bm{i}}^T)^T= \mathfrak{g} ((U_{\bm{i} \odot \bm{e}})_{\bm{e} \in \{ 0,1 \}^{K} \setminus \{ \bm {0} \}}), \ {\bm{i}} \in \mathbb{N}^{K}, \quad \{ U_{\bm{i} \odot \bm{e}} : \bm{i} \in \mathbb{N}^{K}, \bm{e} \in \{ 0,1 \}^{K} \setminus \{ \bm {0} \} \} \stackrel{i.i.d.}{\sim} U[0,1], \] for some Borel measurable map $\mathfrak{g}: [0,1]^{2^{K}-1} \to \mathbb{R}^{1+d}$.

Arguably, one of the most popular estimation methods for such a high-dimensional regression problem is the Lasso tibshirani1996; we refer to vandegeer2011, giraud2015, wainwright2019 as standard references on high-dimensional statistics. Let $N=\prod_{k=1}^K N_k$ denote the total sample size. The Lasso estimate for $\beta_0$ is defined by \[ \hat \beta^\lambda=\operatorname*{arg\,min}_{\beta\in \mathbb{R}^p}\left\{ \frac{1}{N}\sum_{\bm{i}\in [\bm{N}]}(Y_{\bm{i}}-\bm{X}_{\bm{i}}^T\beta)^2+ \lambda\|\beta\|_1\right\}, \] where $\lambda>0$ is a penalty level. We estimate the vector $\bm{f} = (f_{{\bm{i}}})_{{\bm{i}} \in [{\bm{N}}]} = (f(\bm{Z}_{{\bm{i}}}))_{{\bm{i}} \in [{\bm{N}}]}$ by $\hat{\bm{f}}^{\lambda} = ({\bm{X}}_{{\bm{i}}}^{T}\hat{\beta}^{\lambda})_{{\bm{i}} \in [{\bm{N}}]}$. Let $\| \bm{t} \|_{N,2}^2=N^{-1}\sum_{\bm{i}\in [\bm{N}]} t_{\bm{i}}^2$ for $\bm{t} = (t_{\bm{i}})_{\bm{i} \in [ \bm{N} ]}$.

In what follows, we discuss the statistical performance of the Lasso estimate. Following bickel2009, we say that Condition RE$(s,c_0)$ holds (RE refers to “restricted eigenvalue”) if, for a given positive constant $c_0\ge 1$, the inequality \[ \kappa(s,c_0)=\min_{\substack{J\subset\{1,\dots,p\}\\mathbbm{1}\le |J|\le s}}\inf_{\substack{\theta\in \mathbb{R}^p,\,\theta\ne 0 \\\|\theta_{J^c}\|_1\le c_0\|\theta_{J}\|_1}} \frac{\sqrt{s N^{-1}\sum_{\bm{i}\in [\bm{N}]}(\theta^T\bm{X}_{\bm{i}})^2}}{\|\theta_{J}\|_1}>0 \] holds with $J^c= \{1,\dots,p\}\setminus J$. Here for $\theta = (\theta_1,\dots,\theta_p)^{T}$ and $J \subset \{1,\dots,p \}$, $\theta_{J} = (\theta_j)_{j \in J}$.

In addition, to guarantee fast rates for the Lasso, it is important to choose the penalty level $\lambda$ in such a way that $\lambda\ge 2c\| \bm{S}_{{\bm{N}}} \|_\infty$ with $\bm{S}_{\bm{N}}=N^{-1}\sum_{\bm{i}\in [\bm{N}]}\varepsilon_{{\bm{i}}} {\bm{X}}_{{\bm{i}}}$ for some $c > 1$ bickel2009,BC2013. To this end, we shall estimate the $(1-\eta)$-quantile of $2c \|\bm{S}_{{\bm{N}}}\|_{\infty}$ for some small $\eta>0$. We first estimate the error terms $\varepsilon_{{\bm{i}}}$ by pre-estimating $\beta_0$ by the preliminary Lasso estimate $\tilde{\beta} = \hat{\beta}^{\lambda_{0}}$ with penalty $ \lambda^0 =\tau_{n} (n^{-1}\log p)^{1/2}$ for some slowing growing sequence $\tau_{n} \to \infty$. In the following, we take $\tau_n = \log n$ for the sake of simplicity but other choices also work. We apply the multiplier bootstrap to $\tilde{\bm{S}}_{{\bm{N}}}=N^{-1}\sum_{{\bm{i}} \in [{\bm{N}}]}\tilde\varepsilon_{\bm{i}} {\bm{X}}_{{\bm{i}}}$ instead of $\bm{S}_{{\bm{N}}}$.

The H\'{a}jek projection to $\bm{S}_{{\bm{N}}}$ is given by $\sum_{k=1}^{K} N_{k}^{-1}\sum_{k=1}^{N_{k}}\bm{V}_{k,i_{k}}$, where $\bm{V}_{k,i_{k}}$ is given by ${\bm{V}}_{k,i_k}=\mathbb{E}[ \varepsilon_{(1,\dots,1,i_k,1,\dots,1)}{\bm{X}}_{(1,\dots,1,i_k,1,\dots,1)}\mid U_{(0,\dots,0,i_k,0,\dots,0)}]$. We estimate ${\bm{V}}_{k,i_k}$ by \\ $\tilde{\bm{V}}_{k,i_k}=\big(\prod_{k'\ne k} N_{k'} \big)^{-1}\sum_{i_1,\dots,i_{k-1},i_{k+1},\dots,i_K} \tilde\varepsilon_{\bm{i}} {\bm{X}}_{{\bm{i}}}$. Let $\{ \xi_{1,i_{1}} \}_{i_1=1}^{N_{1}}, \dots, \{ \xi_{K,i_{K}} \}_{i_K=1}^{N_K}$ be i.i.d. $N(0,1)$ variables independent of the data, and consider

We propose to choose $\lambda$ as $\lambda=\lambda(\eta)=2c\Lambda_{\bm{N}}^\xi(1-\eta)$, where $\Lambda_{\bm{N}}^\xi(1-\eta)$ denotes the conditional $(1-\eta)$-quantile of $\Lambda_{\bm{N}}^\xi$. We allow $\eta$ to decrease with $n$, i.e, $\eta = \eta_n \to 0$.

The following proposition establishes the asymptotic validity of our choice of $\lambda$ (as $n \to \infty$) under separate exchangeability. In what follows, we understand that $s,p,{\bm{N}},\eta$ are functions of $n$ while other parameters such as $c, q, \underline{\kappa}$ are independent of $n$.

The proof of Proposition (ref) does not follow directly from the results of Section (ref), as we have to take care of the estimation error of the preliminary Lasso estimate $\tilde{\beta}$, which requires extra work.

Condition (iii) in the preceding proposition is a high-level condition on the sample gram matrix. The following proposition provides primitive sufficient conditions for Condition (iii) to hold for the case of $K=2$.

Under Condition (i) of Proposition (ref), $B_{{\bm{N}}} \le \overline{N}^{1/q}D_{{\bm{N}}}$, so that $sB_{{\bm{N}}}^2 \log^4 (p\overline{N}) = o(n)$ reduces to $s\overline{N}^{1/q}D_{{\bm{N}}}\log^{4}(p\overline{N}) = o(n)$, which is implied by Condition (iv) of Proposition (ref).

Simulation studies

In this section, we present simulation studies to evaluate the finite sample performance of the proposed multiplier bootstrap methods.

We first describe the simulation design for separately exchangeable arrays. With $\Sigma_{\bm{Z}}$ denoting the $p \times p$ covariance matrix consisting of elements of the form $4^{-\left\vert r-c \right\vert}$ in its $(r,c)$-th position, separately exchangeable data with $K=2$ indices are generated according to $ {\bm{X}}_{{\bm{i}}} = \frac{1}{4} \left( \bm{Z}_{(i_1,0)} + \bm{Z}_{(0,i_2)} \right) + \frac{1}{2} \bm{Z}_{(i_1,i_2)}, $ where $ \bm{Z}_{{\bm{i}} \odot {\bm{e}}} \sim B N(\bm{0}, \Sigma_{\bm{Z}}) + (1-B)N(\bm{0}, 2\Sigma_{\bm{Z}}) $ and $ B \sim \text{Bernoulli}(0.5) $ independently for ${\bm{i}} \in \{(i_1,i_2) \in \mathbb{N}^2: 1 \le i_1 \le N_1, 1 \le i_2 \le N_2\}$ and ${\bm{e}} \in \{0,1\}^2$. For this data generating design, we run 2,500 Monte Carlo iterations to compute the uniform coverage frequencies of $\mathbb{E}[{\bm{X}}_{{\bm{i}}}]$ for the nominal probabilities of 90% and 95% using our proposed multiplier bootstrap for separately exchangeable arrays with 2,500 bootstrap iterations.

We next describe the simulation design for jointly exchangeable arrays. We shall focus on the the most common case in practice, the dyadic data, i.e. $K=2$. With $\Sigma_{\bm{Z}}$ denoting the $p \times p$ covariance matrix consisting of elements of the form $4^{-\left\vert r-c \right\vert}$ in its $(r,c)$-th position, dyadic samples are generated symmetrically in $i$ and $j$ according to $ {\bm{X}}_{i,j} = \frac{1}{4} \left( \bm{Z}_{(i,0)} + \bm{Z}_{(j,0)} \right) + \frac{1}{2} \bm{Z}_{(i,j)}, $ where $ \bm{Z}_{{\bm{i}} \odot {\bm{e}}} \sim B N(\bm{0}, \Sigma_{\bm{Z}}) + (1-B) N(\bm{0}, 2\Sigma_{\bm{Z}}) $ and $ B \sim \text{Bernoulli}(0.5) $ independently for ${\bm{i}} \in \{(i,j) \in \mathbb{N}^2: 1 \le i < j \le n \}$ and ${\bm{e}} \in \{1\} \times \{0,1\}$. We run 2,500 Monte Carlo iterations to compute the uniform coverage frequencies of ${\bm{S}}_{n}$ for the nominal probabilities of 90% and 95% using our proposed multiplier bootstrap with 2,500 bootstrap iterations.

Table (ref) summarizes simulation results under the separate exchangeability. The columns consist of the dimension $p$ of ${\bm{X}}$ and the two-way sample size $(N_1,N_2)$. The displayed numbers indicate the simulated uniform coverage frequencies for the nominal probabilities of 90% and 95%. For each dimension $p \in \{25,50,100\}$, sample sizes vary as $(N_1,N_2) \in \{(25,25), (50,50),(100,100)\}$. Table (ref) summarizes simulation results under the joint exchangeability. The columns consist of the dimension $p$ of ${\bm{X}}$, and the dyadic sample size $N$. The displayed numbers indicate the simulated uniform coverage frequencies for the nominal probabilities of 90% and 95%. For each dimension $p \in \{25,50,100\}$, sample sizes vary as $n \in \{50, 100, 200\}$.

Observe that, for each simulation design and for each nominal probability, the uniform coverage frequencies approach the nominal probability as the sample size increases. These results support the theoretical property of our multiplier bootstrap method. We ran many other sets of simulations with various designs and sample sizes not presented here, but this observed pattern to support our theory remains invariant across all the different sets of simulations -- see Appendix (ref). In Appendix (ref), we further experiment with the separate exchangeability with $K=3$ indices.

Real data analysis

In this section, we present an empirical application of the method proposed in Section (ref) to constructing uniform confidence bands for the density functions of bilateral trade volumes in the international trade, with a similar motivation to that stated in Graham2019kernel,GrahamNiuPowell2020kernel. Recall that our method extends those by Graham2019kernel in that we can draw uniform confidence bands as opposed to point-wise confidence intervals. From this analysis, we can learn about the evolution of the distributions of international trade volumes over time.

We employ the international trade data used in head2014gravity, that come from the Direction of Trade Statistics (DoTS). This data set contains information about bilateral trade flows among 208 economies for 59 years from 1948 to 2006. In this analysis, we will focus on the relatively recent years, 1990, 1995, 2000 and 2005. Our measure of the bilateral trade volume $Y_{ij}$ is defined as the logarithm of the sum of the trade flow from economy $i$ to economy $j$ and the trade flow from economy $j$ to economy $i$. We perform simulation studies on confidence bands for densities in Appendix (ref), confirm that the method works as desired, and thus use the same software code here to draw confidence bands of the probability density function of $Y_{ij}$. Since there is a probability mass at zero in the international trade volumes, what we estimate is precisely the Lebesgue-Radon-Nikodym derivative of the continuous part of the distribution, rather than the probability density function. Specifically, we use $\hat b(y)$ defined in Section (ref) for estimation, and confidence bands are constructed by setting $\hat a = 1$. That said, we shall call it a density for conciseness.

Figure (ref) illustrates estimates and confidence bands of the density functions of $Y_{ij}$ in each of the years 1990, 1995, 2000 and 2005. Each panel of the figure displays the kernel density estimates in a solid curve and the 95% uniform confidence bands in a gray shade. In addition, we also display the proportion of zero bilateral trade volumes to the left of the kernel density plots so we can get an idea of the complementary proportion that consists the density of the continuously distributed part of the distribution. Although we treat $Y_{ij}$ as the logarithm of the bilateral trade volumes in estimation and inference, we use the original scale on the horizontal axis for ease of reading the graphs.

Observe that the proportion of the zero trade volume is decreasing over time, and the density function is accordingly moving upward over time. Despite this pattern of the changes over time, the shapes of the density functions are rather similar across time in the middle of the distribution. This observation entails a high level of confidence given the reasonably tight confidence bands. On the other hand, notice that the right tail of the distribution becomes fatter as time progresses, implying that there is an increasing number of bilateral trading pairs with very large trade volumes.

Summary

In this paper, we have developed methods and theories for inference about high-dimensional parameters with separately/jointly exchangeable arrays. Building on the high-dimensional CLTs over the rectangles, we have proposed bootstrap methods and established their finite sample validity for both notions of exchangeability. Simulation studies support the theoretical properties of the methods. We have illustrated a couple of applications of the bootstrap methods. First, extending Graham2019kernel, we have applied our method to construction of uniform confidence bands for density functions of dyadic data. Second, we have demonstrated an application of our method to penalty choice for $\ell_1$-penalized regression under the separate exchangeability. As such, the results in the present paper pave the way for a variety of applications to analyses of separately and jointly exchangeable arrays.