Minimax Risk and Uniform Convergence Rates for Nonparametric Dyadic Regression

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

22,984 characters · 4 sections · 26 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

\singlespacing

abstract\singlespacing Let $i=1,\ldots,N$ index a simple random sample of units drawn from some large population. For each unit we observe the vector of regressors $X_{i}$ and, for each of the $N\left(N-1\right)$ ordered pairs of units, an outcome $Y_{ij}$. The outcomes $Y_{ij}$ and $Y_{kl}$ are independent if their indices are disjoint, but dependent otherwise (i.e., “dyadically dependent”). Let $W_{ij}=\left(X_{i}',X_{j}'\right)'$; using the sampled data we seek to construct a nonparametric estimate of the mean regression function $g\left(W_{ij}\right)\overset{def}{\equiv}\mathbb{E}\left[\left.Y_{ij}\right|X_{i},X_{j}\right].$ We present two sets of results. First, we calculate lower bounds on the minimax risk for estimating the regression function at (i) a point and (ii) under the infinity norm. Second, we calculate (i) pointwise and (ii) uniform convergence rates for the dyadic analog of the familiar Nadaraya-Watson (NW) kernel regression estimator. We show that the NW kernel regression estimator achieves the optimal rates suggested by our risk bounds when an appropriate bandwidth sequence is chosen. This optimal rate differs from the one available under iid data: the effective sample size is smaller and $d_W=\mathrm{dim}(W_{ij})$ influences the rate differently.

{\bf JEL codes:} C14

{\bf Keywords:} Networks, Exchangeable Random Graphs, Dyadic Regression, Kernel Regression, Minimax Risk, Uniform Convergence.

\setcounter{page}{1}

\onehalfspacing

Introduction

Let $i=1,\ldots,N$ index a simple random sample of units drawn from some large population. For each unit we observe the vector of regressors $X_{i}$ and, for each of the $N\left(N-1\right)$ ordered pairs of units, or directed dyads, we observe the “dyadic” outcome $Y_{ij}$ (e.g., total exports from country $i$ to country $j$). The outcomes $Y_{ij}$ and $Y_{kl}$ are independent if their indices are disjoint, but dependent otherwise (e.g., exports from Japan to Korea may covary with those from Japan to Vietnam).

Let $W_{ij}=\left(X_{i}',X_{j}'\right)'$; using the sampled data we seek to construct a nonparametric estimate of the mean regression function

equation[equation omitted — 135 chars of source]

We present two sets of results. First, we calculate lower bounds on the minimax risk for estimating the regression function at (i) a point and (ii) under the infinity norm. Second, we calculate (i) pointwise and (ii) uniform convergence rates for the dyadic analog of the familiar Nadaraya-Watson (NW) kernel regression estimator. We show that the NW kernel regression estimator achieves the optimal rates suggested by our risk bounds when an appropriate bandwidth sequence is chosen.

Analogous results are widely available in the i.i.d. setting. For nonparametric regression risk bounds see, for example, stone1980,stone1982 and ibragimov1982, ibragimov1984. tsybakov2008 provides a masterful synthesis of these results, from which we draw in formulating our own proofs.

Uniform convergence of kernel averages with i.i.d. data, as well as stationary strong mixing data, have been studied by, for example, newey1994 and hansen2008 respectively. The latter paper includes additional references to the extensive literature in this area. Our uniform convergence proofs build upon those of hansen2008. Nonparametric density estimation with dyadic data was first considered by graham2019; chiang2019 present uniform convergence results for dyadic density estimators.\footnote{It is possible that the methods of inference presented in chiang2019 could be adapted to our setting.}

Our results provide insight in the structure of dyadic nonparametric estimation problems. Our minimax risk bounds suggest that, $N$, the number of units, not $n \overset{def}{\equiv} N \times (N-1)$, the number of dyadic outcomes, is the relevant “sample size” for dyadic estimation problems. This is consistent with the long standing intuition among empirical researchers that dyadic dependence makes inference less precise (see aronow_samii_assenova_2017 and the references cited therein), as well as with a small, but growing, number of more formal rates-of-convergence results graham2020network.

More surprisingly, we find that the relevant dimension of our estimation problem is just $d_X = \mathrm{dim}(X_i)$, not $d_W = 2d_X$. We provide two intuitions for this fact. The first, described below, stems from the thought experiment underlying our minimax risk bound calculations. The second, arises from the fact that the H\'ajek projection of the NW estimator has a “partial-mean-like" structure. As is well known, averaging over the marginal distribution of some regressors, while holding the remaining ones fixed, improves rates-of-convergence newey1994, linton1995.

graham2020network surveys empirical studies in economics utilizing dyadic data. Interest in, as well as the availability of, such data are growing in economics, other academic fields, and in enterprise settings. This paper provides an initial set of results for nonparametric regression with dyadic data. These results are, of course, of direct interest. They should, as has been true with their i.i.d. predecessors, also be useful for proving consistency of two-step semiparametric M-estimators under dyadic dependence (see chiang2019 for some results on double machine learning with dyadic data).

Lower Bounds on the Minimax Risk

Let $i=1,\ldots,N$ index a simple random sample of units drawn from some large population. The econometrician observes the vector of regressors, $X_{i}$, for each sampled unit as well as the scalar outcome, $Y_{ij}$, for each directed pair of sampled units (i.e., each directed dyad). Let $\mathbf{Z}_N = (X_1, \ldots, X_N, Y_{ij}, 1 \leq i\neq j \leq N)$ be the observable data when $N$ units are sampled. The regression function of interest is (ref) above. The goal is to construct a nonparametric estimate of $g:\mathbb{R}^{d_{W}}\rightarrow\mathbb{R}$ where $d_{W}=2d_{x}$.

We assume that $Y_{ij}$ is generated according to the following conditionally independent dyad (CID) model graham2020network.

equation[equation omitted — 81 chars of source]

Random sampling ensures that $\left(X_{i},U_{i}\right)$ is independent and identically distributed for $i=1,\ldots,N$. We further assume that $\left\{ \left(V_{ij},V_{ji}\right)\right\} _{1\leq i<j\leq N}$ are i.i.d. and indepenent of $\mathbf{X}=\left(X_{1},\ldots,X_{N}\right)'$ and $\mathbf{U}=\left(U_{1},\ldots,U_{N}\right)$. Here $h$ is an unknown function, often called the graphon. This set-up, which can also be derived as an implication of more primitive exchangeability assumptions, has the following implications (see graham2020network,graham2020sparse for additional discussion):

enumerate• The $Y_{ij}$ are relatively exchangeable given the $W_{ij}$. Namely, the conditional distribution of $\mathbf{Y}$ is invariant across permutations of the indices $\sigma : \mathbb{N} \rightarrow \mathbb{N}$ satisfying the restriction $[W_{\sigma(i)\sigma(j)}] \overset{d}{=} [W_{ij}]$: \[ [Y_{ij}] \overset{d}{=} [Y_{\sigma(i)\sigma(j)}]. \] • $Y_{ij}$ and $Y_{kl}$ are independent if their indices are disjoint. • $Y_{ij}$ and $Y_{kl}$ are dependent (unconditionally or conditionally given $X_1,\ldots,X_N$) if they share at least one index in common.

The statistical problem is to estimate the regression function $g$ when the only prior restriction on it is that it belongs to the H\"older class of functions.

definition(H\"older Class) Given a vector $s = (s_1, \ldots, s_{d})$, define $|s| = s_1 + \cdots + s_{d}$ and \[D^s = \frac{\partial^{s_1 + \cdots + s_{d}}}{\partial^{s_1} w_1 \cdots \partial^{s_{d}} w_{d}}.\] Let $\beta$ and $L$ be two positive numbers. The H\"older class $\Sigma(\beta, L)$ on $\mathbb{R}^{d}$ is defined as the set of $l = \floor{\beta}$ times differentiable functions $g:\mathbb{R}^{d} \rightarrow \mathbb{R}$ whose partial derivative $D^s g$ satisfies \[ |D^s g(w) - D^s g(w')| \leq L||w - w'||_{\infty}^{\beta - l}, \quad \forall w, w'\in \mathbb{R}^{d} \] for all $s$ such that $|s| = \floor{\beta}$. $\floor{\beta}$ denotes the greatest integer strictly less than the real number $\beta$.

An estimator $\hat{g}_N$ is a function $w \mapsto \hat{g}_N(w) = \hat{g}_N(w, \mathbf{Z}_N)$ measurable with respect to $\mathbf{Z}$. Our first result establishes a lower bound on the minimax risk for estimating the regression function at a single point and under the infinity norm. We state this result under a Gaussian error assumption, which simplifies the proof.

theorem(Minimax Risk Lower Bound) Suppose that $\beta > 0$ and $L > 0$; $X_i$ is continuously distributed on $\mathbb{R}^{d_X}$ with density $f$ and $\sup_x f(x) \leq B_3 < \infty$; and $Y_{ij}$ is generated according to the following nonparametric regression model: \begin{align*} Y_{ij} & = g\left(W_{ij}\right) + e_{ij}, \quad i \neq j, \end{align*} with $e_{ij} = U_i + U_j + V_{ij}$, $U_i \mathrel{\stackrel{\makebox[0pt]{\mbox{\normalfont\tiny iid}}}{\sim}} \text{N}(0, 1)$, and $V_{ij}\mathrel{\stackrel{\makebox[0pt]{\mbox{\normalfont\tiny iid}}}{\sim}} \text{N}(0, 1)$, then \begin{enumerate}[label=(\roman*)] • For all $w \in \mathbb{R}^{d_W}$, \begin{align*} \liminf_{N\rightarrow \infty} \inf_{\hat{g}_N} \sup_{g \in \Sigma(\beta, L)} \mathbb{E}_{g} \left[N^{\frac{2\beta}{2\beta + d_X}}\left(\hat{g}_N(w) - g(w)\right)^2\right] \geq c_1, \end{align*} where $c_1 > 0$ depends only on $\beta$ and $L$. • \begin{align*} \liminf_{N\rightarrow \infty} \inf_{\hat{g}_N} \sup_{g \in \Sigma(\beta, L)} \mathbb{E}_{g} \left[\left(\frac{N}{\ln N}\right)^{\frac{2\beta}{2\beta + d_X}}\left|\left|\hat{g}_N - g\right|\right|_{\infty}^2\right] \geq c_2, \end{align*} where $c_2 > 0$ also depends only on $\beta$ and $L$. \end{enumerate}

Our proof follows the general recipe outlined in Chapter 2 of tsybakov2008. The lower bound at a point is based on Le Cam's method of two hypotheses. The lower bound under the infinity norm is based on Fano's method of multiple hypotheses.

The key, and novel, step in our proof involves constructing hypotheses close enough to one other in terms of Kullback-Leibler (KL) divergence while being at the same time different enough in terms of the target regression function.

An essential feature of our construction is additive separability of the regression functions. In the hypotheses we consider, $Y_{ij} = k(X_i) + k(X_j) + U_i + U_j + V_{ij}$. Next suppose we also observe $T_i \overset{def}{\equiv} k(X_i) + U_i$. Observe that $(X_i, T_i, i = 1,\ldots, N)$ is sufficient with respect to $(X_i, T_i, i = 1,\ldots, N, Y_{kl}, 1\leq k\neq l \leq N)$ for the parameter $k$.

It is well-known that the optimal rates of convergence for estimating $k$ using iid data $(X_i, T_i, i = 1,\ldots, N)$ are $N^{-\frac{\beta}{2\beta + d_X}}$ pointwise and $\left(\frac{N}{\ln N}\right)^{-\frac{\beta}{2\beta + d_X}}$ for the infinity norm. We expect the rates for estimating $g$ to be no faster than these. The proof of Theorem (ref) makes this intuition rigorous.

Relative to its iid counterpart, there are two distinctive features of Theorem (ref). First, the relevant sample size is not the number of observed dyadic outcomes $n = N \times (N-1)$, but instead the number of sampled units, $N$. Dependence across outcomes sharing indices in common is strong enough to slow down the feasible rate of convergence. Second, although the regression function has $d_W = 2d_X$ arguments, the relevant dimension reflected in the rate of convergence result is just $d_X$ (i.e., just half of what might naively be expected).

The form of our constructed hypotheses provides one intuition for this second finding: clearly the relevant dimension of the problem of estimating $k(x)$ is just $d_X$. Relatedly this finding is consistent with those of linton1995 in their analysis of additively separable, but otherwise nonparametric, regression functions (see also newey1994).

The pairwise structure of dyadic data results in apparent data abundance (sample $N$ agents, but observe $O(N^2)$ outcomes!). This abundance is both illusory, in the sense that the effective sample size is indeed just $N$, and real, in the sense that availability of the pairwise outcome data allows for an effective reduction in the dimensionality of the problem via partial mean like average (as in newey1994 and linton1995 in a different context).

Kernel Estimator of Dyadic Regression

In this section we study the properties of a specific nonparametric regression estimator. Namely, the dyadic analog of the well-known Nadaraya-Watson (NW) kernel regression estimator. While our results are specific to this estimator, they could, for example, be extended to apply to local linear regression hansen2008.

The dyadic NW kernel regression estimator is

equation[equation omitted — 134 chars of source]

where \[ K_{ij, N}(w) := \frac{1}{h_N^{d_W}}K\left(\frac{W_{ij} - w}{h_N}\right), \] $K$ is a fixed multivariate kernel function, and $h_N$ is a vanishing bandwidth sequence.

We first develop a sequence of results useful for bounding the variance of kernel objects of the form

align[align omitted — 97 chars of source]

and then apply these results to the NW regression estimator. We then bound the NW estimator's bias and combine the two sets of results to formulate a risk bound.

Variance Bound and Uniform Convergence

Here we are interested in bounding the deviation of $\hat{\Psi}_N(w)$ from its mean. We begin with a presentation of our maintained assumptions.

assumption[Model] The data generating process is as described in Section (ref) with \begin{enumerate}[label=(\roman*)] • $X_i$ continuously distributed with marginal density $f(x)$ s.t. $\sup_{x\in \mathbb{R}^{d_X}}f(x) \leq B_3 < \infty$; • $\sup_{x_1, x_2 \in \mathbb{R}^{d_X}} \mathbb{E}\left[|Y_{12}|^2\big|(X_1, X_2)=(x_1, x_2)\right] \cdot f(x_1) f(x_2) \leq B_4 < \infty$, \\$\sup_{x_1, x_2, x_3 \in \mathbb{R}^{d_X}} \mathbb{E}\left[|Y_{12}Y_{13}|\big|(X_1, X_2, X_3)=(x_1, x_2, x_3)\right] \cdot f(x_1) f(x_2) f(x_3) \leq B_5 < \infty$. \end{enumerate}

Condition (i) is a standard condition in the context of kernel estimation, while (ii) ensures that various second moments appearing in our variance calculations are finite.

assumption[Kernel, Part A] $\sup_{w \in \mathbb{R}^{d_W}} |K(w)| \leq K_{\text{max}} < \infty$, $\int_{w\in \mathbb{R}^{d_W}} |K(w)|\mathrm{d}w \leq B_1 < \infty$, and $\sup_{x\in \mathbb{R}^{d_X}}\int |K(x, x')|\mathrm{d}x' \leq B_2 < \infty$.

Assumption (ref) is satisfied by many widely-used multivariate kernel functions. Our first result holds under Assumptions (ref) and (ref).

theorem[Variance Bound] Under Assumptions (ref) and (ref), and the bandwidth condition $Nh_N^{d_X} \rightarrow \infty$ as $N\rightarrow \infty$, there exists a constant $M_0 < \infty$ such that for $N$ sufficiently large \[ \operatorname{Var}\left(\hat{\Psi}_N(w)\right) \leq \frac{M_0}{Nh_N^{d_X}} \] for all $w\in \mathbb{R}^{d_W}$.

A proof is available in the appendix. Mirroring our risk bound results, two features of Theorem (ref) merit comment. First, $N$ not $n = N \times (N-1)$ appears in the denominator. This is due to the effects of dependence across dyads sharing units in common. Second, the relevant dimension of the problem is $d_X$, not $d_W=2d_X$, this reflects the U-statistic like structure of kernel weighted averages and the partial mean like averaging this structure induces.

To establish uniform convergence, we need additional moment conditions on $Y_{ij}$ as well as some smoothness conditions on the kernel $K$. As in hansen2008, we require the kernel to either have bounded support and be Lipschitz or have bounded derivatives and an integrable tail. See hansen2008 for additional discussion about these conditions. As with Assumption (ref) above, most commonly used kernels satisfy these conditions.

assumption[Regularity Condition] \begin{enumerate}[label=(\roman*)] • For some $s>2$, $\mathbb{E}|Y_{12}|^s < \infty$ and $\sup_{x_1, x_2 \in \mathbb{R}^{d_X}} \mathbb{E}\left[|Y_{12}|^s\big|(X_1, X_2)=(x_1, x_2)\right] \cdot f(x_1, x_2) \leq B_{4,s} < \infty$; • For some $\Lambda_1 < \infty$ and $L < \infty$, either (a) or (b) holds \begin{itemize} • $K(w) = 0$ for $||w|| > L$, and $|K(w) - K(w')| \leq \Lambda_1||w-w'||$ for all $w, w'\in \mathbb{R}^{2d}$ • $K(w)$ is differentiable, $\left|\left|\frac{\partial}{\partial w} K(w)\right|\right| \leq \Lambda_1$, where $\left|\left|\frac{\partial}{\partial w} K(w)\right|\right| = \left|\left|\left(\frac{\partial}{\partial w_1} K(w), \ldots \frac{\partial}{\partial w_{2d}} K(w)\right)\right|\right|_{\infty}$, and for some $\nu > 1$, $\left|\left|\frac{\partial}{\partial w} K(w)\right|\right| \leq \Lambda_1 ||w||^{-\nu}$ for $||w|| > L$. \end{itemize} \end{enumerate}

Part (ii) coincides with Assumption 3 in hansen2008. This assumption implies that for all $||w_1 - w_2|| \leq \delta \leq L$,

align*[align* omitted — 53 chars of source]

where $K^*(u)$ satisfies Assumption (ref). If case (a) holds, then $K^*(u) = 2d\Lambda_1 \mathds{1}(||u||\leq 2L)$. If case (b) holds, then, $K^*(u) = 2d[\Lambda_1 \mathds{1}(||u||\leq 2L) + \left( ||u|| - L\right)^{-\nu} \mathds{1}(||u||> 2L)]$. In both cases $K^*$ is bounded and integrable and therefore satisfies Assumption (ref).

Define \[ a_N := \left(\frac{\ln N}{N h_N^{d_X}}\right)^{1/2}. \]

theorem[Weak Uniform Convergence] Under Assumptions (ref), (ref), (ref), and the bandwidth conditions $\max\left\{\min\left\{(a_N h_N^{2d_X})^{-\frac{1}{s-1}}, [N^2(\ln (\ln N))^2 \ln N]^{\frac{1}{s}}\right\}, a_N^{-\frac{1}{s-1}}\right\} \ll \min\left\{a_N^{-1}, \frac{N}{\ln N} h_N^{\frac{3}{2}d_X} \right\}$ and $\frac{N }{\ln N}h_N^{d_X} \rightarrow \infty$, we have for any $q > 0$, $c_N = N^q$, \[ \sup_{||w|| \leq c_N} \left|\hat{\Psi}_N(w) - \mathbb{E}\hat{\Psi}_N(w)\right| = O_P(a_N). \]

This theorem establishes uniform convergence of $\hat{\Psi}_N(w)$ to its mean in probability over an expanding set with radius growing at a polynomial rate.

In the proof, we decompose $\hat{\Psi}_N(w)$ into two parts \[ \hat{\Psi}_N(w) = \tilde{\Psi}_N(w) + R_N(w), \] in which $\tilde{\Psi}_N(w) = \frac{1}{N(N-1)}\sum_{1\leq i \neq j\leq N} Y_{ij}\cdot \mathds{1}\left(|Y_{ij}| < \tau_N\right) K_{ij, N} $ is a truncated version of $\hat{\Psi}_N(w)$ with a carefully chosen threshold parameter $\tau_N$ and $R_N(w)$ is a residual. The boundedness induced by this truncation is technically convenient as it facilitates the application of various concentration inequalities. To establish concentration of $\tilde{\Psi}_N$, we apply Bernstein inequality to its H\'{a}jek Projection (i.e., to the first-order terms in the Hoeffding decomposition) and apply arcones1993's concentration inequalities for degenerate U-statistics to the second-order terms in the Hoeffding decompositon. Both these bounds requires the truncation threshold to be small enough. To bound the magnitude of the residual $R_N$, we can either apply a triangular inequality to bound the sup of its first moment or use the Borel-Cantelli Lemma to bound its probability of being nonzero. Both these bounds requires the truncation threshold to be large.

A proper truncation threshold satisfying both requirements exists only if the bandwidth sequence satisfies the condition \[ \max\left\{\min\left\{(a_N h_N^{2d_X})^{-\frac{1}{s-1}}, [N^2(\ln (\ln N))^2 \ln N]^{\frac{1}{s}}\right\}, a_N^{-\frac{1}{s-1}}\right\} \ll \min\left\{a_N^{-1}, \frac{N}{\ln N} h_N^{\frac{3}{2}d_X} \right\}. \] The complicated form of this condition is technical in nature. When all (conditional) moments of $Y_{12}$ are bounded, such that $s = \infty$ (of Assumption (ref) above), this condition simplifies to $\frac{N}{\ln N} h_N^{\frac{3}{2} d_X} \gg 1$.

In order to state the weak uniform convergence result for the kernel regression estimator $\hat{g}_N$, we need additional smoothness assumptions on the kernel. As in other applications of kernel estimation, these assumptions are employed for bias reduction purpose.

assumption[Kernel, Part B] \begin{align*} \int_{\mathbb{R}^{d_W}} w_1^{l_1} \cdots w_{d_W}^{l_{d_W}} K(w) dw = \left\{ \begin{array}{ll} 1, & if l_1 = \cdots = l_{d_W} = 0\\ 0, & if (l_1, \ldots, l_{d_W})' \in \mathbb{Z}_+^{d_W} and l_1 + \cdots + l_{d_X} < \beta \end{array} \right. \end{align*}

We can now give a uniform convergence result for the NW regression estimator under dyadic dependence over a sequence of expanding sets.

theoremSuppose $f_W, g \in \Sigma(\beta, L)$ and $\delta_N = \inf_{||w||\leq C_N} f_W(w)>0$, $\delta_N^{-1} a_N^*\rightarrow 0$ where $a_N^* := \left(\frac{\ln N}{N h_N^{d_X}}\right)^{1/2} + h_N^{\beta}$. Under the Assumptions of Theorem, (ref) and Assumption (ref) \begin{align*} \sup_{||w||\leq C_N} |\hat{g}_N(w) - g(w)| & = O_p(\delta_N^{-1} a_N^*). \end{align*} The optimal convergence rate is \begin{align*} \sup_{||w||\leq C_N} |\hat{g}_N(w) - g(w)| & = O_p\left(\delta_N^{-1} \left(\frac{\ln N}{N }\right)^{\frac{\beta}{2\beta + d_X}}\right). \end{align*}

As in the iid case, the KW estimator achieves the optimal rate suggested by Theorem (ref) for a compact set with $C_N = C$. If we look at a sequence of expanding sets approaching the entire space $\mathbb{R}^{d_W}$, then there is an additional penalty term $\delta_N$ due to the presence of the denominator $f_W(w)$.