Sparse network asymptotics for logistic regression

Sparse network asymptotics for logistic regression

Let $i=1,\ldots,N$ index a random sample of consumers and $j=1,\ldots,M$ a random sample of products. For each consumer-product pair $ij$ we observe $Y_{ij}=1$ if consumer $i$ purchases product $j$ and $Y_{ij}=0$ otherwise. Let $W_{i}$ be a vector of observed consumer attributes, $X_{j}$ a vector of product attributes and $n=M+N$ the total number of sampled consumers and products. The conditional probability that consumer $i$ buys product $j$ is given by

where $Z_{ij}\overset{def}{\equiv}z\left(W_{i},X_{j}\right)$ is a vector of known functions of $W_{i}$ and $X_{j}$, $\alpha_{0,n}$ an “intercept” parameter (which may vary with $n$), $\beta_{0}$ a vector of fixed “slope” parameters, and $e\left(\cdot\right)$ a known increasing function mapping the real line into the unit interval. Below I will emphasize the logit case with $e\left(v\right)=\exp\left(v\right)/\left[1+\exp\left(v\right)\right]$; this case is convenient and dominates empirical work, but nothing which follows hinges essentially upon it.

I am interested in settings where both the number of consumers, $N$, and the number of products, $M$, are very large. To motivate this focus, consider a large book retailer. Such a retailer may service many customers and also stock many books. Let $x$ be the attribute vector associated with a newly released book, $\hat{\theta}_{n}=\left(\hat{\alpha}_{n},\hat{\beta}'\right)'$ estimates of the parameters in (ref), constructed from some training sample, and $\hat{e}_{i}\left(x\right)=e\left(\hat{\alpha}_{n}+z\left(W_{i},x\right)'\hat{\beta}\right)$ the predicted probability that agent $i$ purchases a book of type $X_{j}=x$. With this knowledge the retailer might use

to predict total unit sales for the new book. This prediction could be useful for making wholesale purchase decisions. Other objects of interest include various average partial effects Chamberlain_HBE84,Wooldridge_IIEM05.

In this paper I present a method of estimating the coefficient vector $\theta_{0,n}=\left(\alpha_{0,n},\beta_{0}'\right)'$ as well as one for conducting inference on it. I also explore estimation and inference for aggregate effects, like (ref), as well as for average effects. The econometric framework outlined below is designed to accommodate two peculiarities of the setting described above that appear to be important in practice and also consequential for inference.

First, consider predicting whether randomly sampled consumer $i$ purchases book $j$, say The Clue in the Crossword Cipher, the forty-fourth novel in the celebrated Nancy Drew mystery series. Knowledge of the frequency with which other consumers $k=1,\ldots,i-1,i+1,\ldots,N$ purchase book $j$ will generally alter the econometrician's prediction of whether $i$ also purchases book $j$. That is, for any $k\neq i$, \[ \Pr\left(\left.Y_{ij}=1\right|Y_{kj}=1\right)>\Pr\left(Y_{ij}=1\right) \] or that $Y_{i_{1}j_{1}}$and $Y_{i_{2}j_{2}}$ will covary whenever the two transactions correspond to a common book (such that $j_{1}=j_{2}$).

Similarly, if the econometrician knew that consumer $i$ was a frequent book buyer, she might conclude that this consumer is also more likely to purchase some other book (relative to the average consumer). That is $Y_{i_{1}j_{1}}$and $Y_{i_{2}j_{2}}$ will also covary whenever the transactions correspond to a common buyer (such that $i_{1}=i_{2}$).

Importantly, dependence across $Y_{i_{1}j_{1}}$ and $Y_{i_{2}j_{2}}$ when $\left\{ i_{1},j_{j}\right\} $ and $\left\{ i_{2},j_{2}\right\} $ share a common buyer or book index may hold even conditional on observed consumer, $W_{i}$, and product attributes, $X_{j}$. Some consumers may have latent attributes (i.e., not contained in $W_{i}$) which induce them to buy many books and some books may be especially popular (for reasons not captured adequately by $X_{j}$). It might be, for example, that \[ \Pr\left(\left.Y_{ij}=1\right|Y_{kj}=1,W_{i},W_{k},X_{j}\right)>\Pr\left(\left.Y_{ij}=1\right|W_{i},X_{j}\right). \]

The structured form of dependence across the elements of $\left[Y_{ij}\right]_{1\leq i\leq N,1\leq j\leq M}$ described above is a feature of separately exchangeable random arrays Aldous_JMA81,Hoover_WP79. The inferential implications of such dependence, in the context of subgraph counts, were first considered by Holland_Leinhardt_SM76 almost fifty years ago. Bickel_et_al_AS11 make an especially important recent contribution in this area. In the context of regression models, the inferential implications of dyadic dependence have been considered by, among others, Fafchamp_Gubert_JDE07, Cameron_Miller_WP14, Aronow_et_al_PA17, Graham_Book_DR_Chap2020, Davezies_et_al_AS20 and Menzel_arXiv17 (see Graham_HBE2020 for a review and references). Dyadic dependence will generate distinct issues here.

The second peculiarity explored here is suggested by the observation that, even when presented with the opportunity to purchase many books, the typical customer will only purchase a few. Similarly, a retailer only sells a few copies of most titles in a given year. Put differently personal libraries are generally small and the market share of most books is, for all practical purposes, infinitesimally small. These observations have implications for what types of asymptotic approximations are likely to be useful in practice. In this paper I consider sequences where both $N$ and $M$ grow at the same rate such that, recalling that $n=M+N$, \[ M/n\rightarrow\phi\in\left(0,1\right) \] as $n\rightarrow\infty$.

Let $\rho_{0,n}\overset{def}{\equiv}\mathbb{E}\left[e\left(\alpha_{0,n}+Z_{ij}'\beta_{0}\right)\right]$ be probability that a randomly sampled consumer purchases a randomly sampled book. The average number of books purchased by the average consumer is then

In network parlance $\lambda_{0,n}^{c}$ corresponds to average consumer degree. If $\rho_{0,n}$ is bounded away from zero, then $\lambda_{0,n}^{c}\rightarrow\infty$ as $N,M\rightarrow\infty$. This implies that the number of actual book purchases and the number of possible book purchases should be of equal order. This not true in practice.\footnote{There are tens of millions of print titles available on, for example, Amazon, even consumers who buys hundreds of books in a year are completing only very small fraction of all possible purchases.} To develop a distribution theory which is concordant with the empirical regularity that consumers only purchase a small number of books (and, similarly, that retailers only sell a small number of copies of any given title) I let $\alpha_{0,n}\rightarrow-\infty$ at a rate which ensures that $\lambda_{0,n}^{c}$ converges to a non-zero and bounded constant $0<\lambda_{0}^{c}<\infty$ as $N,M\rightarrow\infty$. In language of networks I consider bi-partite graphs which are sparse.

The asymptotic analysis in this paper is, to my knowledge, novel, but it does connect with two important areas of prior research by others. The first is the literature on subgraph counts and dyadic regression cited above. However, with the partial exception of Bickel et al.'s Bickel_et_al_AS11 analysis of acyclic subgraph counts, this work has been, starting with Holland_Leinhardt_SM76, limited to to dense networks.\footnote{Graham_EM17 and Jochmans_JBES18 also considered regression in the context of graphs which are sparse in the limit. Both these papers utilize conditional likelihood type ideas; this has the effect of “conditioning away” some of the dependence which is central to the analysis below.} The second connection is to the literature on “rare events” analysis King_Zeng_PA2001; an area of special concern in political science and epidemology, but also increasingly relevant in economics (especially in the era of “Big Data”). An interesting feature of Theorem (ref) below is that it contains Wang's HaiYing_ICML2020 recent result for logistic regression with rare events and iid data as a special case.

The formal analysis of this paper is confined to bipartite networks, but adapting it to directed and/or undirected networks would be straightforward. Several consumer demand settings might be appropriately modeled with the methods described in this paper. Examples include (i) the listening behavior of streaming music service customers and (ii) the purchase behavior of big box store customers. A limitation vis-a-vis these applications, is that the basic set-up explored here is not useful for understanding complementary and substitution patterns across products Lewbel_Nesheim_Cemmap2019. Other possible applications include modeling plant locations in an industry where firms typically operate multiple plants (here $i=1,\ldots,N$ would index firms and $j=1,\ldots,M$ locations; see KPMG_CAReport2016). Other many-to-many matching problems that have drawn economists' interest include (i) bank-firm lending relationships Marotta_et_al_PlosOne2015, (ii) the matching of venture capital with start-ups Bengtsson_Hsu_JBV2015, and (iii) supply chain settings with strong bi-partite structure (e.g., automakers and parts supplies as in Fox_QE2018). When $i=1,\ldots,N$ and $j=1,\ldots,M$ index the same units with $M=N$, applications include the modeling of “rare events” in international relations data, such as interstate wars King_Zeng_PA2001. More generally, the methods developed in this paper, with minimal adaptation, can be used for link prediction in any sparse network setting: bi-partite, directed or undirected.\footnote{The methods outlined here are not appropriate for use in one-to-one matching settings.} There are numerous applications of link prediction in the other social sciences, the bench sciences as well as in industry.

In what follows random variables are denoted by capital Roman letters, specific realizations by lower case Roman letters and their support by blackboard bold Roman letters. That is $Y$, $y$ and $\mathbb{Y}$ respectively denote a generic random draw of, a specific value of, and the support of, $Y$. A “0” subscript on a parameter denotes its population value and may be omitted when doing so causes no confusion. In what follows I use graph, network and purchase graph to refer to $\mathbf{Y}\overset{def}{\equiv}\left[Y_{ij}\right]_{1\leq i\leq N,1\leq j\leq M}$. All graph theory terms and notation used below are standard Chartand_Zhang_GT2012.

Population and sampling assumptions

Let $i\in\mathbb{N}$ index consumers in an infinite population of interest. Associated with each consumer is the vector of observed attributes $W_{i}\in\mathbb{W}=\left\{ w_{1},\ldots,w_{J}\right\} .$ Let $j\in\mathbb{M}$ index products in a second infinite population of interest. The model is a two population one Graham_Imbens_Ridder_JBES18. Associated with each product is the vector of characteristics $X_{i}\in\mathbb{X}=\left\{ x_{1},\ldots,x_{K}\right\} $. The finite support assumption on $\mathbb{W}$ and $\mathbb{X}$ is not essential, but simplifies the discussion of exchangeability below.

Let $\sigma_{w}:\mathbb{N}\rightarrow\mathbb{N}$ be a permutation of a finite number of consumer indices which satisfies the restriction

Restriction (ref) implies that $\sigma_{w}$ only permutes indices across observationally identical consumers (i.e., with the same values of $W$). Let $\sigma_{x}:\mathbb{M}\rightarrow\mathbb{M}$ be an analogously constrained permutation of a finite number of product indices. Adapting the terminology of Crane_Towsner_JSL18, I assume that the purchase graph is $W$-$X$-exchangeable

Here $\overset{D}{=}$ denotes equality of distribution. One way to think about (ref) is as a requirement that any probability law for $\left[Y_{ij}\right]_{i\in\mathbb{N},j\in\mathbb{M}}$ should attach equal probability to all purchase graphs which are isomorphic as vertex-colored graphs. Here I associate $W_{i}$ and $X_{j}$ with the color of the corresponding consumer and product vertices in the overall purchase graph. Virtually all single-population micro-econometric models assume that agents are exchangeable, restriction (ref) extends this idea to the two-population setting considered here. Our probability law for the model should not change if we re-label observationally identical units.

Graphon

It is well-known that exchangeability implies restrictions on the structure of dependence across observations in the cross-section setting deFinetti_AN1931. Aldous_JMA81, Hoover_WP79 and Crane_Towsner_JSL18 showed that exchangeable random arrays also exhibit a special dependence structure. Let $\mu$, $\left\{ \left(W_{i},A_{i}\right)\right\} _{i\geq1}$, $\left\{ \left(X_{j},B_{j}\right)\right\} _{j\geq1}$ and $\left\{ V_{ij}\right\} _{i\geq1,j\geq1}$ be sequences of i.i.d. random variables, additionally independent of one another, and consider the purchase graph $\left[Y_{ij}^{*}\right]_{i\in\mathbb{N},j\in\mathbb{M}}$, generated according to

with $h:\left[0,1\right]\times\mathbb{W}\times\mathbb{X}\times\left[0,1\right]^{2}\rightarrow\left\{ 0,1\right\} $ a measurable function, henceforth referred to as a graphon (we can normalize $\mu$, $A_{i}$, $B_{j}$ and $V_{ij}$ to have support on the unit interval, uniformly distributed, without loss of generality).

The results of Crane_Towsner_JSL18, which extend the earlier work of Aldous_JMA81 and Hoover_WP79, show that, for any $W$-$X$-exchangeable random array $\left[Y_{ij}\right]_{i\in\mathbb{N},j\in\mathbb{M}}$, there exists another array $\left[Y_{ij}^{*}\right]_{i\in\mathbb{N},j\in\mathbb{M}}$, generated according to (ref), such that the two arrays have the same distribution. An implication of this result is that we may use (ref) as a nonparametric data generating process for $\left[Y_{ij}\right]_{i\in\mathbb{N},j\in\mathbb{M}}$.

Inspection of (ref) indicates that exchangeability implies a particular pattern of dependence across the elements of $\left[Y_{ij}\right]_{i\in\mathbb{N},j\in\mathbb{M}}$. In particular $Y_{i_{1}j_{1}}$ and $Y_{i_{2}j_{2}}$ may covary whenever $i_{1}=i_{2}$ or $j_{1}=j_{2}$; this covariance may be present even conditional on consumer and product attributes. This is, of course, precisely the dependence structure discussed in the introduction.

Sampling process

Let $\mathbf{Y}=\left[Y_{ij}\right]_{1\leq i\leq N,1\leq j\leq M}$ be the observed $N\times M$ matrix of consumer purchase decisions. Let $\mathbf{W}$ and $\mathbf{X}$ be the associated matrices of consumer and product regressors. I assume that $\mathbf{Y}$ is the adjacency matrix associated with the subgraph induced by a random sample of consumers and products from a $W$-$X$-exchangeable infinite population graph. Let $G_{\infty,\infty}$ denote this population network. Associated with this network is some graphon (ref). Let $\mathcal{V}_{c}$ and $\mathcal{V}_{p}$ denote the set of consumers and products randomly sampled by the econometrician from $G_{\infty,\infty}$. We have $\mathbf{Y}$ equal to the adjacency matrix of the network:

An implication of (ref), (ref) and (ref) is that we may proceed `as if' the adjacency matrix in hand was generated according to \[ Y_{ij}=h\left(\mu,W_{i},X_{j},A_{i},B_{j},V_{ij}\right) \] for $i=1,\ldots,N$ and $j=1,\ldots,M$. The marginal probability of the event, random consumer $i$, purchases random product $j$, is thus

Let $\left\{ G_{N,M}\right\} $ be a sequence of networks indexed by, respectively, the cardinality of the sampled consumer and product index sets, $N=\left|\mathcal{V}_{c}\right|$ and $M=\left|\mathcal{V}_{p}\right|$. The average number of products purchased per consumer, or average consumer degree,

will diverge as $M\rightarrow\infty$ when $\rho_{0}>0$. Likewise the average number of times a given product is purchased, or average product degree,

will also diverge as $N\rightarrow\infty$. A consequence of this divergence is that the number of possible purchases, and the number of actual purchases, will be of equal order. In practice, however, only a small fraction of all possible purchases are made. To capture this feature of the real world in our asymptotic approximations requires a slightly more elaborate thought experiment; which I outline next.

Instead of considering a sequence of graphs sampled from a fixed population, I consider a sequence of graphs sampled from a corresponding sequence of populations. The sequence of networks $\left\{ G_{N,M}\right\} $ is one where both $N$ and $M$ grow at the same rate such that, recalling that $n=M+N$, \[ M/n\rightarrow\phi\in\left(0,1\right) \] as $n\rightarrow\infty$. For each $N,M$ the graphon describing the infinite population sampled from is

This sequence of graphons/populations $\left\{ h_{N,M}\right\} $ has the property that network density \[ \rho_{0,N,M}=\mathbb{E}_{N,M}\left[h_{N,M}\left(\mu,W_{i},X_{j},A_{i},B_{j},V_{ij}\right)\right] \] may approach zero as $n\rightarrow\infty$. Under this setup the order of $\lambda_{0,N,M}^{c}=M\rho_{0,N,M}$ and $\lambda_{0,N,M}^{p}=N\rho_{0,N,M}$ will depend upon the speed with which $\rho_{0,N,M}$ approaches zero as $n\rightarrow\infty$. Here I use the notation $\mathbb{E}_{N,M}\left[\cdot\right]$ to emphasize that the probability law used to compute expectations may vary with the sample size.

As in other exercises in alternative asymptotics, indexing the population data generating process by the sample size is not meant to capture a literal feature of how the data are generated, rather it is done so that the limiting properties of the model share important features -- in this case sparseness -- with the actual finite network in hand. In other settings such an approach has led to more useful asymptotic approximations, a premise I maintain here Staiger_Stock_Em1997.

Composite likelihood estimator

The estimation target is the regression function of $Y_{ij}$ given $X_{i}$ and $W_{j}$. This is a predictive function and may, or may not, have structural economic meaning as well (see Graham_HBE2020). I assume that this regression function takes the parametric form

where $Z_{ij}\overset{def}{\equiv}z\left(W_{i},X_{j}\right)$ is a finite vector of known functions of $W_{i}$ and $X_{j}$. It would be interesting to extend what follows to semiparametric regression models, but this is not done here.

Assumption (ref) formalizes the population and sampling set-up of the previous section.

To allow the probability of making a purchase decline with $n$, let $\alpha_{0,n}=\ln\left(\alpha_{0}/n\right).$ This gives, after some manipulation,

and hence an expression for average consumer degree (ref) of, recalling that $M/n\approx\phi$,

which converges to a bounded constant as $n\rightarrow\infty$ as long as $\mathbb{E}\left[\exp\left(Z_{ij}'\beta_{0}\right)\right]<\infty$. By allowing $\alpha_{0,n}\rightarrow-\infty$ as $n\rightarrow\infty$ we ensure that, in the limit, the bipartite graph $\mathbf{Y}=\left[Y_{ij}\right]$ is sparse. Similar devices are used by Owen_JMLR2007 and HaiYing_ICML2020 to model “rare events” in cross-sectional binary outcome data.

The compact support assumption on $Z_{ij}$ is not essential, but simplifies the proofs. Here I focus on estimation of $\theta_{n}=$$\left(\alpha_{n},\beta'\right)'$ with $\alpha_{n}=\ln\left(\alpha/N\right)$. Observe that $\theta=\left(\alpha,\beta'\right)'$ does not vary with $n$, while $\theta_{n}$ does. Define $\theta_{0,n}=\left(\alpha_{0,n},\beta_{0}'\right)'$ with $\alpha_{0,n}=\ln\left(\alpha_{0}/n\right)$. Let $\bar{\alpha}=\sup\mathbb{A}$, the parameter space for $\theta_{n}$, $\Theta_{n}=\left(-\infty,\ln\left(\bar{\alpha}\right)\right]\times\mathbb{B}$, is the one induced by $\Theta=\mathbb{A}\times\mathbb{B}$ and the mapping from $\theta$ to $\theta_{n}$.

To estimate $\theta_{0,n}$we maximize the composite log-likelihood function \[ \hat{\theta}_{n}=\arg\underset{\theta\in\Theta_{n}}{\max}\thinspace L_{n}\left(\theta\right) \] with $L_{n}\left(\theta\right)\overset{def}{\equiv}\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M}l_{ij}\left(\theta\right)$ and \[ l_{ij}\left(\theta\right)=\left(2Y_{ij}-1\right)R_{ij}'\theta-\ln\left[1+\exp\left[\left(2Y_{ij}-1\right)R_{ij}'\theta\right]\right] \] the logit kernel function and $R_{ij}=\left(1,Z_{ij}'\right)'$. Observe that $\hat{\theta}_{n}$ is simply the coefficient vector associated with a logistic regression of $Y_{ij}$ onto a constant and $Z_{ij}$ using all $NM$ dyads in the network. Although, by virtue of Assumption (ref) above, $L_{n}\left(\theta\right)$ correctly represents the marginal (conditional) probability of $Y_{ij}$ for each element of $\mathbf{Y}$, it does not accurately reflect the dependence structure across these elements; hence the term “composite likelihood”. See Lindsay_CM88 an introduction to estimation by composite likelihood and Graham_HBE2020 for discussion in the contexts of network model estimation.

Sparse network asymptotics

If $\alpha_{0,n}$ equals a fixed constant, then $\rho_{0,n}$ - network density -- will also be fixed such that the network will be dense in the limit. The limit distribution of $\hat{\theta}_{n}$ under such “dense network asymptotics” was derived by Graham_HBE2020. More general results for dyadic M-estimators under dense network asymptotics, including results on the bootstrap, can be found in Menzel_arXiv17 and Davezies_et_al_AS20. None of these results apply here. To derive a result that does apply, begin with the mean value expansion \[ \sqrt{n}\left(\hat{\theta}_{n}-\theta_{0,n}\right)=\left[nH_{n}\left(\bar{\theta}_{n}\right)\right]^{+}\times n^{3/2}S_{n}\left(\theta_{0,n}\right). \] where

with $s_{ij}\left(\theta\right)=\frac{\partial l_{ij}\left(\theta\right)}{\partial\theta}=\left(Y_{ij}-e_{ij}\left(\theta\right)\right)R_{ij}$ and $e_{ij}\left(\theta\right)=e\left(\alpha+Z_{ij}'\beta\right)$, corresponds to the score vector of the composite likelihood and

the associated Hessian matrix. Here $\bar{\theta}_{n}$ is a mean value between $\theta_{0,n}$ and $\hat{\theta}_{n}$ which may vary from row to row.

Lemma (ref), stated and proved in Appendix (ref), shows that, after re-scaling by $n$, that $nH_{n}\left(\theta\right)$ converges uniformly to

An intuition for why $H_{n}\left(\theta\right)$ needs to be rescaled to ensure convergence is that, under sparse network asymptotics, information accrues at a slower rate: the effective sample size is not $NM=\left(n^{2}\right)$, but rather $O\left(n\right)$. I return to this point briefly at the end of the paper.

Assumption (ref) is a standard identification condition (see, for example, Amemiya_AE85). This assumption, in conjunction with Lemma (ref), gives the linear approximation \[ \sqrt{n}\left(\hat{\theta}_{n}-\theta_{n}\right)=\Gamma_{0}^{-1}\times n^{3/2}S_{n}\left(\theta_{0,n}\right)+o_{p}\left(1\right). \] To derive the limit distribution of $\sqrt{n}\left(\hat{\theta}_{n}-\theta_{n}\right)$ I show that the distribution $n^{3/2}S_{n}\left(\theta_{0,n}\right)$ is well-approximated by a Gaussian random variable. The main tool used is a martingale CLT for triangular arrays. That the variance stabilizing rate for $S_{n}\left(\theta_{0,n}\right)$ is $n^{3/2}$, like the need to rescale the Hessian, is non-standard. The need to “blow up” $S_{n}\left(\theta_{0,n}\right)$ at a faster than $\sqrt{n}$ rate is a consequence of the fact that the summands in $S_{n}\left(\theta_{0,n}\right)$ are $O\left(n^{-1}\right)$ since $\alpha_{0,n}\rightarrow-\infty$ as $n\rightarrow\infty$.

A detailed proof of Theorem (ref), stated below, is provided in Appendix (ref). Here I outline the main arguments. Begin with the following three part decomposition of the score vector

where

with $\bar{s}_{ij}\left(\theta\right)=\bar{s}\left(W_{i},X_{j},A_{i},B_{j};\theta\right)$ with $\bar{s}\left(w,x,a,b;\theta\right)=\mathbb{E}\left[\left.s_{ij}\left(\theta\right)\right|W_{i}=w,X_{j}=x,A_{i}=a,B_{j}=b\right]$ and

with $\bar{s}_{1}^{c}\left(w,a;\theta\right)=\mathbb{E}\left[\bar{s}\left(w,X_{j},a,B_{j};\theta\right)\right]$ and $\bar{s}_{1}^{p}\left(x,b;\theta\right)=\mathbb{E}\left[\bar{s}\left(W_{i},x,A_{i},b;\theta\right)\right]$.

Decomposition (ref) also features in Graham_Book_DR_Chap2020 and Menzel_arXiv17.\footnote{It is also implicit in the elegant proof in Bickel_et_al_AS11.} It can be derived by first projecting $S_{n}\left(\theta\right)$ on to $\mathbf{A}=\left[A_{i}\right]_{1\leq i\leq N}$, $\mathbf{W}=\left[W_{i}\right]_{1\leq i\leq N}$, $\mathbf{B}=\left[B_{j}\right]_{1\leq j\leq M}$, and $\mathbf{X}=\left[X_{i}\right]_{1\leq j\leq N}$ as follows:

Next observe that $\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M}\bar{s}_{ij}\left(\theta\right)$ is a two sample U-Statistic, albeit one defined partially in terms of the latent variables $A_{i}$ and $B_{j}$. Equation (ref) corresponds to the Hajek Projection of this U-statistic onto (separately) $\left\{ \left(W_{i}',A_{i}\right)\right\} _{i=1}^{N}$ and $\left\{ \left(X_{j}',B_{j}\right)\right\} _{j=1}^{M}$. Equation (ref) is the usual Hajek Projection error term.

Define $\phi_{n}=M/n$, $\bar{s}_{1ni}^{c}\overset{def}{\equiv}\bar{s}_{1i}^{c}\left(\theta_{0,n}\right),$ $\bar{s}_{1nj}^{p}\overset{def}{\equiv}\bar{s}_{1j}^{p}\left(\theta_{0,n}\right)$ and also $\bar{s}_{nij}\overset{def}{\equiv}\bar{s}_{ij}\left(\theta_{0,n}\right)$. Similarly let $S_{n}=S_{n}\left(\theta_{0,n}\right)$ and so on. Applying the variance operator to $S_{n}$ yields:

where

In the dense case $\Sigma_{1n}^{c}$, $\Sigma_{1n}^{p},$ $\Sigma_{2n}$ and $\Sigma_{3n}$ are all constant in $n$; hence the asymptotic properties of $S_{n}$ coincide with those of $U_{1n}$. Since $U_{1n}$ is a sum of independent random variables a standard argument gives

as long as $\Sigma_{1}^{c}$ and/or $\Sigma_{1}^{p}$ are non-zero.

Under the sparse network asymptotics considered here the order of $\Sigma_{1n}^{c}$, $\Sigma_{1n}^{p},$ $\Sigma_{2n}$ and $\Sigma_{3n}$ varies with $n$. This affects the order of the four variance terms in (ref) and, consequently, which components of $S_{n}$ contribute to its asymptotic properties. In Appendix (ref) I show the order of the four terms in (ref) are, respectively,

Since $\Sigma_{1}^{c}$ and $\Sigma_{1}^{p}$ are both $O\left(\rho_{n}^{2}\right)=O\left(n^{-2}\right)$ we can multiply them by $n^{2}$ to stabilize them. Define $\tilde{\Sigma}_{1}^{c}$ to be the limit of $n^{2}\Sigma_{1n}$ and $\tilde{\Sigma}_{1}^{p}$ to be the limit of $n^{2}\Sigma_{1n}^{p}$. Similarly we can define $\tilde{\Sigma}_{3}$ to be the limit of $n\Sigma_{3n}$, all as $n\rightarrow\infty$. Normalizing (ref) by $n^{3/2}$ therefore gives

where I also use the fact that $\Sigma_{2n}=O\left(n^{-2}\right)$.

Under sparse network asymptotics both $U_{1n}$ and $V_{n}$ matter. In Appendix (ref) I show that $U_{1n}+V_{n}$ is a martingale difference sequence (MDS) to which a martingale CLT can be applied; Theorem (ref) then follows.

Theorem (ref) indicates that under sparse network asymptotics there are additional sources of sampling variation in $\sqrt{n}\left(\hat{\theta}_{n}-\theta_{n}\right)$ relative to those which appear in the dense case. Not incorporating these into inference procedures will lead to tests with incorrect size and/or confidence intervals with incorrect coverage. A further advantage of considering sparse network asymptotics is that Theorem (ref) remains valid even under degeneracy of the graphon, $h_{N,M}\left(\mu,W_{i},X_{j},A_{i},B_{j},V_{ij}\right)$. For example, if the graphon is constant in $A_{i}$ and $B_{j}$ such that $Y_{ij}$and $Y_{ik}$ do not covary conditional on covariates (and likewise for $Y_{ji}$ and $Y_{ki}$), then $\tilde{\Sigma}_{1}^{c}=\tilde{\Sigma}_{1}^{p}=0$, but Theorem (ref) nevertheless remains valid. In contrast, under dense network asymptotics, degeneracy -- as elegantly shown by Menzel_arXiv17 -- generates additional complications. In that case the variance of $U_{1n}$ is identically equal to zero, while that of $U_{2n}$ and $V_{n}$ are of equal order. In some cases, the behavior of $U_{2n}$ may even induce a non-Gaussian limit distribution (see vanderVaart_ASBook00). In the sparse network cases, $U_{2n}$ is always negligible relative to $V_{n}$. Furthermore $V_{n}$ is well approximated -- after suitable scaling -- by a Gaussian distribution.

Extensions and discussion

In this section I connect Theorem (ref) to prior work on rare events logistic analysis, sketch some results about the estimation of aggregate and average effects, and, finally, close with a few ideas about possible areas of additional research.

Rare events with iid data

King_Zeng_PA2001 discuss, with a focus on finite sample bias, the behavior of logistic regression under “rare events” with iid data. Evidently binary choice analyses where the marginal frequency of positive events is quite small are common in empirical work.\footnote{The King_Zeng_PA2001 has close to four thousand citations on Google Scholar.} The properties of logistic regression under sequences where the number of “events” becomes small (i.e., “rare”) relative to the sample size as it grows were recently characterized by HaiYing_ICML2020 (see also Owen_JMLR2007). The main result in HaiYing_ICML2020 coincides with a special case Theorem (ref) above. To see this observe that if the graphon is constant in $A_{i}$ and $B_{j}$, then $\bar{s}_{nij}$ will be identically equal to zero for all $1\leq i\leq N$ and $1\leq j\leq M$. In this scenario there is no “dyadic dependence” (after conditioning on $W_{i}$ and $X_{j}$) and $\tilde{\Sigma}_{1}^{c}=\tilde{\Sigma}_{1}^{p}=0$. Inspection of the calculations in Appendix (ref) also reveals that, in this case, we further have an information-matrix type equality of $\tilde{\Sigma}_{3}=\Gamma_{0}$. Under these conditions Theorem 1 specializes to \[ \sqrt{n}\left(\hat{\theta}-\theta_{n}\right)\overset{D}{\rightarrow}\mathcal{N}\left(0,\Gamma_{0}^{-1}\right), \] as $n\rightarrow\infty$. This is precisely, up to some small differences in notation, the result given in Theorem 1 of HaiYing_ICML2020.\footnote{HaiYing_ICML2020 scales by the square root of the number of events or “ones” in the dataset. This is, of course, of the same order as $n$ as defined here. This difference leads to a minor difference in our two expressions for $\Gamma_{0}$. Making these adjustments the results coincide.}

In his analysis HaiYing_ICML2020 emphasizes that information accumulates more slowly under “rare event asymptotics”. In the present setting this is reflected in the need to rescale the Hessian matrix by $n$ to achieve convergence (see Lemma (ref) in Appendix (ref)). In the network setting dyadic dependence additionally slows down the rate of convergence Graham_Niu_Powell_WP2019. If a researcher is working with a sparse network and concerned about dyadic dependence, then she should base inference on Theorem (ref). If the graphon is degenerate or, more strongly, the elements of $\left[Y_{ij}\right]_{1\leq i\leq N,1\leq j\leq M}$ are, in fact, iid, then her inferences will remain valid (since Theorem (ref) specializes to the “rare events” result of HaiYing_ICML2020 in that case).

Aggregate effects

Define $e_{ni}\left(x\right)=e\left(\alpha_{0,n}+z\left(W_{i},x\right)'\beta_{0}\right)$ and $\hat{e}_{ni}\left(x\right)=e\left(\hat{\alpha}_{n}+z\left(W_{i},x\right)'\hat{\beta}\right)$ and consider an estimate of total unit sales for a product with attribute vector $X_{j}=x$ of

\[ \hat{\gamma}_{n}\left(x\right)=\sum_{i=1}^{N}\hat{e}_{ni}\left(x\right). \] Under dense network asymptotics this statistic would diverge as $n\rightarrow\infty$. Under sparse network asymptotics the sum $\sum_{i=1}^{N}e_{ni}\left(x\right)$ behaves like an average because its summands are $O\left(N^{-1}\right)$. Consequently, $\hat{\gamma}_{n}\left(x\right)$ has a well-defined probability limit of

This probability limit reflects the boundedness of average product degree in sparse networks (i.e., in expectation, total sales of a product are finite, even asymptotically). This result also holds within a sub-population of products with characteristics $X_{j}=x$. Parameter (ref) corresponds to a conditional version of $\lambda_{0,n}^{p}$, the average product degree parameter defined in Section (ref) above.

To derive the rate of convergence and limit distribution of $\hat{\gamma}_{n}\left(x\right)$ I proceed in the usual way. A mean value expansion and Theorem (ref) together yield

By the conditional mean zero property of the score function, the two terms in (ref) are asymptotically uncorrelated. The variance of the first term in (ref) is \[ \mathbb{V}\left(\sqrt{n}\sum_{i=1}^{N}\left\{ e_{ni}\left(x\right)-\gamma_{0}\left(x\right)\right\} \right)=O\left(n^{2}\left(1-\phi_{n}\right)\rho_{n}^{2}\right)=O\left(1\right), \] whereas the Jacobian in the second term has the approximation

Defining $\Phi_{0}\left(x\right)\overset{def}{\equiv}$$\alpha_{0}\left(1-\phi\right)\left(

'\right)$, $\Lambda_{0}\left(x\right)=\underset{n\rightarrow\infty}{\lim}n^{2}\left(1-\phi\right)\mathbb{V}\left(e_{ni}\left(x\right)-\gamma_{0}\left(x\right)\right)$, and $\Omega_{0}=\Gamma_{0}^{-1}\left[\frac{\tilde{\Sigma}_{1}^{c}}{1-\phi}+\frac{\tilde{\Sigma}_{1}^{p}}{\phi}+\frac{\tilde{\Sigma}_{3}}{\phi\left(1-\phi\right)}\right]\Gamma_{0}^{-1}$ suggest that \[ \sqrt{n}\left(\hat{\gamma}_{n}\left(x\right)-\gamma_{0}\left(x\right)\right)\rightarrow N\left(0,\Lambda_{0}\left(x\right)+\Phi_{0}\left(x\right)\Omega_{0}\Phi_{0}\left(x\right)'\right) \] as $n\rightarrow\infty$. Aggregate effects are estimable with the same degree of precision as the logit coefficients themselves.

Average partial effects

Next consider estimating the average marginal effect of unit increases in the elements of $Z_{ij}$ on making a purchase:

Recall that $e_{nij}=e\left(\alpha_{0,n}+Z_{ij}\beta_{0}\right)$ and $\hat{e}_{nij}=e\left(\hat{\alpha}_{n}+Z_{ij}'\hat{\beta}\right)$. Interest in average partial effects of this type is widespread in modern micro-econometric empirical research Blundel_Powell_WC03,Wooldridge_IIEM05. Since (ref) is an average of summands, each of which is $O\left(n^{-1}\right)$, we might expect some variance reduction relative to the aggregate case just discussed. In a certain sense, this conjecture appears to be correct.

Define $\gamma_{0,n}=\mathbb{E}_{N,M}\left[e_{nij}\left(1-e_{nij}\right)Z_{ij}\right]=O\left(\rho_{n}\right)$; a mean-value expansion and some re-scaling yields

As above, the conditional mean zero property of the score function ensures that the two terms in (ref) are asymptotically uncorrelated. We rescale the estimate and parameter using $\rho_{n}$ since $\gamma_{0,n}\rightarrow0$ as $n\rightarrow\infty$ Bickel_et_al_AS11. The need to rescale the Jacobian in (ref) stems from the observation that

Next observe that first term in (ref) is a two-sample U-Statistics. A Hoeffding variance decomposition gives \[ \mathbb{V}\left(\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M}e_{nij}\left(1-e_{nij}\right)Z_{ij}\right)=\frac{\Lambda_{1n}^{c}}{N}+\frac{\Lambda_{1n}^{p}}{M}+\frac{1}{NM}\left[\Lambda_{2n}-\Lambda_{1n}^{c}-\Lambda_{1n}^{p}\right] \] with

Inspection indicates that $\Lambda_{1n}^{c}=O\left(\rho_{n}^{2}\right)$, $\Lambda_{1n}^{p}=O\left(\rho_{n}^{2}\right)$ and $\Lambda_{2n}^{c}=O\left(\rho_{n}^{2}\right)$ and hence that \[ n^{3}\mathbb{V}\left(\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M}e_{nij}\left(1-e_{nij}\right)Z_{ij}\right)=\frac{\tilde{\Lambda}_{1}^{c}}{1-\phi}+\frac{\tilde{\Lambda}_{1}^{p}}{\phi}+O\left(n^{-1}\right). \]

Putting these calculations together suggests that

\[ \rho_{n}n^{3/2}\left(\frac{\hat{\gamma}_{n}-\gamma_{0,n}}{\rho_{n}}\right)\overset{D}{\rightarrow}\mathcal{N}\left(\frac{\tilde{\Lambda}_{1}^{c}}{1-\phi}+\frac{\tilde{\Lambda}_{1}^{p}}{\phi}+\Phi_{0}\Omega_{0}\Phi_{0}'\right) \] with $\Phi_{0}\overset{def}{\equiv}\alpha_{0}\mathbb{E}\left[\exp\left(Z_{ij}'\beta_{0}\right)\left[

\right]\right]$.

If we set $T=NM=O\left(n^{2}\right)$, then we have that $T^{1/4}\left(\hat{\theta}_{n}-\theta_{0,n}\right)$ has a Gaussian limit distribution. The rate of convergence of $\hat{\theta}_{n}$ toward $\theta_{0,n}$ is slow. For average partial effects we need to rescale in order ensure a meaningful probability limit. Let $\gamma_{0,n}^{*}=\gamma_{0,n}/\rho_{n}$ and similarly for $\hat{\gamma}_{n}$; the result above implies that $T^{1/4}\left(\hat{\gamma}_{n}^{*}-\gamma_{0,n}^{*}\right)$ is also Gaussian. In this sense the rates-of-convergence for the logit coefficients and their APEs coincide. However, if we think in terms of the resulting implied approximation to the finite sample distribution of the two parameter estimates, we have $\mathbb{V}\left(\hat{\theta}_{n}\right)=O\left(T^{-1/2}\right)=O\left(n^{-1}\right)$, but $\mathbb{V}\left(\hat{\gamma}_{n}\right)=O\left(T^{-3/2}\right)=O\left(n^{-3}\right)$. In this sense inference on APEs appears to be more precise.

Areas for additional research

For empirical researchers the main implication of this paper is to use an estimate for the variance of $S_{N}$ that includes all components -- even ones that are negligible under certain asymptotic sequences -- when constructing standard errors. This is not a new idea. In the context of U-statistics it goes back to Hoeffding_AMS48. It is implicit in Holland_Leinhardt_SM76 in their work on subgraph counts; see also the recent work on dyadic regression by Fafchamp_Gubert_JDE07, Cameron_Miller_WP14 and Aronow_et_al_PA17, as well as that on density weighted average derivatives by Cattaneo_et_al_ET14. However, the small amount of extant formal limit theory for dyadic regression (cited earlier) suggests different approaches to variance estimation. This paper has outlined an asymptotic framework that provides formal justification for one of the leading “practical” approaches to inference in the presence of dyadic dependence. Graham_Book_DR_Chap2020 discusses variance estimation for dyadic regression in detail, advocating a variant of the estimate proposed by Fafchamp_Gubert_JDE07, Cameron_Miller_WP14 and Aronow_et_al_PA17. Theorem (ref) provides a formal justification for this recommendation.

Many outstanding questions remain. Can the above framework be generalized to a generic dyadic M-estimation setting? What is the “general” notion of “sparseness” needed for this? Extensions to semiparametric regression models are also of interest. In recent work, Menzel_arXiv17 and Davezies_et_al_AS20 propose bootstrap procedures for dyadic regression. Are these procedures also valid under sparse network asymptotics and, if not, how might they be adapted to be so? The aggregate and average effect examples sketched above suggest that the systematic exploration of policy analysis questions -- considered under dense network asymptotics by Graham_HBE2020 -- would be interesting. Finally, although it seems likely that -- in the absence of imposing more structure -- that the composite maximum likelihood estimator is efficient, this is currently only a conjecture.