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.
101,251 characters · 19 sections · 75 citation commands
Statistical inference in large multi-way networks
\affil[$\dagger$]{{CREST, ENSAE, Institut Polytechnique de Paris, France}} \affil[$\S$]{{ESSEC Business School, France}}
{\bf JEL codes} C13 C31 C33 C55
Network data, which are becoming available at increasingly granular levels, are receiving a great deal of attention in economics, see Grah_dePa_Book_2020. In this paper, we consider “multi-way” networks that involve interactions between entities of different nature such as importing and exporting countries, buyers and suppliers, teachers and schools, doctors and patients. The strength of the interactions in such networks is commonly measured at disaggregated levels, e.g., industries or products for trade data, consultations and medical procedures for health data, patents and citations for innovation data, etc. Multi-way network data, sometimes referred to as “polyadic”, are indexed by multidimensional indices that represent the relevant dimensions in each case, for instance exporter, importer, product, and time in the trade example.
To model connections in multi-way networks and control for unobserved heterogeneity along various dimensions, recent applied research has gradually considered models with richer structures of fixed effects, involving higher-dimensional interactions. For instance, three-way gravity models, with exporter-year, importer-year, and exporter-importer fixed effects, are common in the modern trade literature.\footnote{Recent studies recognize the economic importance of the sector or product levels, potentially leading to even richer structures of fixed effects, breinlich2024trade and Delb_Dina_WD_2020.} Yet as the structure of fixed effects becomes more complex, maximum likelihood estimators may be plagued by incidental parameter problems, see fernandez2016individual and weidner2021bias. Specifically, as the sample size grows, so too does the number of nuisance parameters representing the fixed effects, possibly creating bias in standard maximum likelihood estimation of the parameters of interest, as first described by Neym_Scot_Ecma_1948.
In this paper, we propose a novel estimator that does not suffer from the incidental parameter problem and handles multi-way gravity models. In the spirit of graham2017econometric, we regard data through the lens of graph theory. The key insight is that certain configurations of outcomes within subgroups of observations, which we call polyads, have identical sufficient statistics for the fixed effects, making their relative likelihood independent of the fixed effects. Our framework differs from graham2017econometric in two important dimensions. First, while Graham models undirected graphs, we use multipartite graphs to model multi-way networks. Second, Graham's network formation model considers only the extensive margin, i.e., the probability that potential links are realized. By contrast, we study the strength of connections in weighted networks, thus modeling both the intensive and extensive margins.
This study is connected to the strand of the gravity literature starting with silva2006log.\footnote{Their seminal paper shows that traditional log-linear OLS estimation suffers from bias under heteroskedasticity, particularly when many flows are zero, which greatly motivated the adoption of Poisson models for gravity. The properties of Poisson pseudo-maximum likelihood estimators have first been investigated by gourieroux1984pseudo in the absence of fixed effects, with one-way fixed effect panel data applications being pioneered by HHG. Recently, chen2024logs argue that log-like transformations can also distort the interpretation of coefficients as percentage effects, since they depend on the units of the outcome.} Our method compares favorably with recent econometric studies along several dimensions. First, it accommodates multi-way models with an arbitrary number of node groups, in contrast to debiasing methods such as fernandez2016individual, jochmans2017two, and weidner2021bias, which are restricted to two- and three-way structures. Second, unlike PPML estimation, our estimator has no incidental parameter problem by construction. Even relative to bias-corrected PPML procedures, our approach remains advantageous: the corrections may themselves be biased in finite samples, as weidner2021bias argue, or computationally prohibitive in large datasets, as in zylkin2024bootstrap. Third, compared with approaches such as charbonneau2012multiple, our estimator better exploits the available variability while remaining computationally feasible, and the convexity of our loss function delivers strong numerical performance relative to more general semiparametric method-of-moments procedures (e.g., jochmans2017two, yang2023three).
Our study is also connected to the literature on discrete choice models for panel and network data (Rasc_1960, Ande_JRSS_1970,Cham_Hand_1984, Hono_Kyri_Ecma_2000, Magn_Ecma_2004). Presenting the conditional likelihood methods used in these settings, the recent review of Dano_Hono_Weid_WP_2025 highlights how identification strategies relate to difference-in-differences approaches. Specifically, they provide the differencing vectors that are valid to identify the parameters of interest.\footnote{Muris_Pakel_WP_2025 follow this approach to study the formation of triadic networks. They introduce a hexad logit estimator that extends the tetrad logit estimator of graham2017econometric.} We proceed the same way for count data and multi-way networks. The polyad estimator can be thought of as a nonlinear version of a difference-in-differences estimator. Contrary to the above cited literature, the polyad method handles count data and recovers both the existence and intensity of relationships. In the special case where the count data is right-censored at one, i.e., where the intensive margin is ignored, the polyad estimator maximizes the same conditional likelihood as graham2017econometric and the studies mentioned in this paragraph .
For researchers working with sparse networks --- i.e., where most potential connections are not realized ---, our approach offers a distinct computational advantage: polyads can be constructed by looping over pairs of edges with strictly positive counts (i.e., realized connections), allowing the estimation procedure to scale with the number of observed relationships rather than the number of potential relationships. This is a significant step forward, since the usage of tetrad-based methods has been limited by its computational cost: the available methods for computing tetrad-based statistics, which work only for two- or three-way models, either (i) require looping over all pairs of edges, including unrealized connections (e.g., graham2017econometric, Muris_Pakel_WP_2025), or (ii) rely on matrix multiplications that remain slow for sparse networks (e.g., jochmans2017two). In particular, as exemplified by our experiments, our computational implementation enables the use of the polyads method on large administrative datasets. See Section (ref) for a more detailed discussion.
We establish consistency and asymptotic normality under mild assumptions, extending the current theoretical framework in two main directions: relaxing sampling assumptions and better exploiting convexity. First, unlike graham2017econometric,graham2020dyadic, we do not require the existence of a distribution to sample the nodes, thus treating the fixed effects as true parameters and revealing fundamental geometrical properties that a graph should satisfy for consistency and asymptotic normality (see Assumptions (ref) and (ref)). Second, we obtain consistency without boundedness assumptions by modifying classical results from newey1994large under the light of convex analysis tools from rockafellar-1970a and asymptotic statistics results from Andersen1982. We also avoid the existence of a limiting risk for asymptotic normality exploring results from victor, which are based on MR1026303, MR1186263.
The practical limitations of our method merit clear statement. First, we do not consider interdependencies between observations beyond those captured by fixed effects and observed covariates. Second, our approach is designed specifically for count data, rather than continuous weights. Third, the method becomes computationally inefficient in dense networks where most potential relationships are realized. Within these constraints, however, our estimator provides a powerful tool for inference in multi-way networks with high-dimensional fixed effects structures.
The remainder of the paper proceeds as follows. Section (ref) introduces the Poisson model and the structure of fixed effects. Section (ref) introduces the polyad estimator. Section (ref) establishes its theoretical properties, demonstrating consistency and asymptotic normality. Section (ref) develops the computational implementation, emphasizing how the algorithm efficiently exploits sparsity. Section (ref) documents the finite-sample properties of the estimator using artificial data and healthcare claims data.
We consider a count variable $Y_{i_1i_2\cdots i_D}\in\mathbb{N}$ that is indexed by $(i_1, \dots, i_D)\in \mathcal{I} = [n_1] \times \cdots \times [n_D]$. Using the $D$-dimensional index ${\bf i} = (i_1, \dots, i_D) \in \mathcal{I}$, we represent the variable as $Y_{\bf i}$.
Throughout the paper, we think of $Y=(Y_{\bf i})_{\bf i}\in \mathbb{N}^\mathcal{I}$ as a random $D$-partite graph, with the sets $[n_1],[n_2],\dots[n_D]$ representing the nodes of each category, the multidimensional index ${\bf i}$ representing a potential (hyper-)edge of a the graph, and $Y_{{\bf i}}$ being the number of connections along edge ${\bf i}$, see the concrete examples below. We denote by $E$ the set of positive edges, i.e., the set of $D$-dimensional indices ${\bf i}$ such that $Y_{{\bf i}}>0$. The graph is sparse when the data contains many zeros, a case where our method delivers especially good results.
We assume that the dependent variable $Y_{{\bf i}}$ depends on a set of $p$ explanatory variables $X_{\bf i} \in \mathbb{R}^p$ and a set of fixed effects. A level of fixed effect is represented by a proper subset $g$ of $[D]$. By abuse of notation, we set $g({\bf i}) = (i_d : d \in g)$ and the fixed effects for level $g$ are denoted as $\theta_{g({\bf i})}^g \in \mathbb{R}$. The structure of the fixed effects in the model is represented by a collection $\mathcal{G}$ of fixed effects levels. Of particular interest to us is the structure $\mathcal{G}^{{\rm max}}$ consisting of the $D$ subsets of $[D]$ of cardinal $D-1$; in this particular case, each level of fixed effect absorbs the variations of $Y_{\bf i}$ in all but one dimension of ${\bf i}$. The set of all fixed effects is denoted by $\theta^\mathcal{G} = \{\theta^g_\rho : g\in \mathcal{G}, \rho\in g({\bf i}) \}$.
The log-likelihood of $(\beta, \theta^{\cal G})$ in the model (ref) at the observed graph $y$ is
involves a potentially high number of fixed effects. The MLE estimator of the parameter of interest $\beta_\star$ has been shown to be asymptotically biased for $D\geq 3$, see weidner2021bias as well as the experiments presented in Section (ref).
The following examples show how our framework encompasses gravity models and other classical econometric models.
To avoid the incidental parameter problem mentioned above, we first condition the likelihood on a sufficient statistics for the fixed effects given by the degrees of the nodes, see Subsection (ref). This approach has been followed in two-way contexts by charbonneau2012multiple, graham2017econometric and jochmans2018semiparametric. We construct a loss function based on a set of `directions' that generate all the variability in the data for given degrees.
Consider a level of fixed effect $g\in\mathcal{G}$. For $\rho\in g(\mathcal{I})$, the fixed effect $\theta_{\rho}^g$ enters the log-likelihood (ref) only through the quantity \[ \delta^g_\rho(y) = \sum_{{\bf i} \in \mathcal{I} : g({\bf i}) = \rho } y_{\bf i}, \] which we call the degree of $\rho$ relative to the fixed effect level $g$ in the graph $y$. This quantity generalizes the notion of degree used in graham2017econometric, capturing the total number of connections between edges ${\bf i}$ for which $g({\bf i}) = \rho$.\footnote{Consider for instance Example (ref), where the index $i_1,i_2,i_3$ are denoted $i,j,t$ as in many gravity models. The degree of (4,5) relative to the fixed effect level $u_{ij}$ is $\delta^{1,2}_{4,5}(y)=\sum_{t} y_{45t}$.} The set of degrees for the fixed effects level $g \in \mathcal{G}$ is denoted by $\delta^g(y) = \left[ \delta^g_\rho(y) \right]_{\rho \in g(\mathcal{I})}$. The set of all degrees is denoted by $\delta(y)=\left[ \delta^g(y) \right]_{g\in \mathcal{G}}$.
As announced above, we condition the likelyhood on the set of degrees $\delta(Y)$:
Rewriting the conditional likelihood as
shows that it does not depend on the fixed effects $\theta^\mathcal{G}$ regardless of their structure $\mathcal{G}$. In the next subsection, we characterize the support of the distribution of the $Y$ conditional on all degrees $\delta(Y)$.
An edge of a polyad $\xi$ is an index ${\bf i}=(i_d)_d\in{\cal I}$ such that $i_d\in\{j_d, j_d^\prime\}$ for all $d\in[D]$. We denote by ${\cal E}(\xi)$ the set of all edges of $\xi$. Any polyad has $|{\cal E}|= 2^D$ edges. In other words, a polyad $\xi$ induces a subgraph of $Y$ made of $2^D$ edges with weights $(y_{{\bf i}})_{{\bf i}\in {\cal E}(\xi)}$.
Polyads as introduced in Definition (ref) are generalizations to the $D$-dimensional framework of tetrads from charbonneau2012multiple, graham2017econometric, jochmans2018semiparametric defined for $D=2$. The total number of polyads, i.e. the size of $\Xi$, is $\prod_{d=1}^D n_d(n_d-1)$. In the upper left corner of Figure (ref) we exemplify a polyad on $D=2$,
In particular, $s_\xi({\bf i}) =0$ when ${\bf i}$ is not an edge of $\xi$ and $s_\xi({\bf i})\in\{\pm1\}$ otherwise. The main purpose of the sign function $s_\xi$ is that it gives the signs of the diff-in-diff property stated below. The sign function is represented in the upper right corner of Figure (ref), each edge has a sign associated to it, notice that the signs sum zero for all axis and that the sign of ${\bf j}$ is always 1.
We now define a class of transformations indexed by polyads. These transformations act on $D$-partite graphs with integer edge weights, i.e., on the set $\mathbb{Z}^\mathcal{I}$. Recall that we see the dependent variable $Y \in \mathbb{N}^\mathcal{I}$ as a graph with nonnegative edge weights. This discrepancy plays an important role in the analysis developed below.
For any polyad $\xi$, the transformation $T_\xi$ alters only the subgraph of $y$ induced by the polyad $\xi$. In other words, $y'_{\bf i}=y_{\bf i}$ for all ${\bf i}\notin {\cal E}(\xi)$. The proposition below states the polyad transformations preserve degrees and, conversely, that they allow to generate all graphs sharing the same degrees as a given graph.\footnote{The converse result is stated in Proposition (ref) only for $\mathcal{G}=\mathcal{G}^{{\rm max}}$. In the appendix, we consider any fixed effect structure $\mathcal{G}$.}
According to Proposition (ref), the conditioning set in the likelihood (ref) can be written as \[ \left\{\,Y \,|\, \delta(Y)= \delta(y)\,\right\}= \left\{\, Y= T_{\xi_1}^{r_1} \circ \cdots \circ T_{\xi_m}^{r_m} (y) : m\geq 0, (\xi_1,\dots,\xi_m)\in\Xi^m, (r_1,\dots,r_m)\in\mathbb{Z}^m \,\right\}, \] which can be represented only at a prohibitively high computational cost.
Rather than attempting to exploit the data variations within this whole set, we propose to exploit variations in the directions induced by each polyad separately, i.e., to restrict attention to sequences of polyads of length $m=1$. Because the count variable $Y_{\bf i}$ takes nonnegative values, the transformations of the graph $T_\xi^r(Y)$ that do not belong to $\mathbb{N}^\mathcal{I}$ are irrelevant because they occur with zero probability. We thus further restrict the conditioning set. For any polyad $\xi\in\Xi$, we introduce the nonnegative integers $m_\xi(y)$ and $M_\xi(y)$ given by \[ m_\xi(y) = \bigwedge_{{\bf i} : s_\xi({\bf i}) = 1} y_{\bf i} \ \ \ \text{ and }\ \ \ M_\xi(y) = \bigwedge_{{\bf i} : s_\xi({\bf i}) = -1} y_{\bf i}. \] The range of integers $r$ such that $T_\xi^r(y)\in\mathbb{N}^\mathcal{I}$ is $\{-m_\xi(y), \dots, M_\xi(y)\}$. Accordingly, we define the orbit $\mathcal{O}_\xi(y)$ of a graph $y\in\mathbb{N}^\mathcal{I}$ with respect to $\xi$ as
All graphs $y$ in the orbit $\mathcal{O}_\xi(Y)$ of the observed graph have the same degrees as $Y$, $\delta(y)=\delta(Y)$, and coincide with $Y$ except in the subgraph induced by $\xi$, $Y_{\bf i}=y_{\bf i}$ for ${\bf i}\notin\mathcal{E}(\xi)$.
We introduce the loss function at the level of the polyad $\xi$:
A polyad $\xi$ is not informative if the orbit $\mathcal{O}_\xi(Y)$ is a singleton, i.e., if both $m_\xi(Y)$ and $M_\xi(Y)$ are zero. For non-informative polyads, we have $\ell_\xi(y|X,\beta) =0$ for all $\beta \in \mathbb{R}^p$. We can thus restrict attention to informative (or “active”) polyads $\xi$ for which the corresponding orbit $\mathcal{O}_\xi(Y)$ has room for potential variation in the data, i.e., contains at least two distinct elements. Formally, active polyads satisfy: $|\mathcal{O}_\xi(Y)| = m_\xi(Y)+M_\xi(Y)+1\geq 2$. To simplify notations, we denote the transformed graph $T^r_\xi(y)$ by $y^r$, hence $y^r_{\bf i}=y_{\bf i}+r s_\xi({\bf i})$. The same computation as in (ref) yields
As already explained, all graphs in the orbit $\mathcal{O}_\xi(Y)$ coincide with $Y$ except in the subgraph induced by the polyad $\xi$. Formally, for ${\bf i}\notin\mathcal{E}(\xi)$, we have $s_\xi({\bf i})=0$ and hence $y_{\bf i}^r=y_{\bf i}$ for all $r$ between $-m_\xi$ and $M_\xi$. In the above sums over edges ${\bf i}$, we can thus restrict attention to edges ${\bf i}\in\mathcal{E}(\xi)$. From (ref), the loss function associated with the polyad $\xi$ is
where the generalized “difference-in-differences” (DiD) operator $\widetilde{X}_\xi$ is defined as
In Example (ref), consider the polyad \[ \xi=
. \] In this two-way example, we get the usual DiD formula: \[ \widetilde{X}_\xi= X_{22}- X_{21} -(X_{12} - X_{11})= X_{22}- X_{21} -X_{12} + X_{11}. \] In Example (ref), consider the polyad \[ \xi=
. \] In this three-way example, we get (the opposite of) a “triple difference” formula:
Finally, we call $\Xi_a$ the set of informative (or active) polyads and form the loss function taking $y = Y$:
Two immediate observations will play an important role in our analysis. First, because the LogSumExp function is convex,\footnote{Recall that the LogSumExp function $\mathbb{R}_{+*}^S\rightarrow\mathbb{R}$\ is given by $\mbox{LSE}(U_0,\cdots,U_{S-1};S) = \ln\left(\sum_{s=0}^{S-1} e^{U_s}\right)$.} the function $\ell_\xi(y | X, \beta)$ is convex in $\beta$ for any polyad $\xi$ and hence the loss function is convex. Second, given two active polyads $\xi$ and $\xi'$, the terms $\ell_\xi(Y|X,\beta)$ and $\ell_{\xi'}(Y|X,\beta)$ in the above some are not independent if the two polyads share at least one edge, i.e., if $|\mathcal{E}(\xi)\cap\mathcal{E}(\xi')|\geq 1$.
We are now in a position to define our loss function and the associated estimator.
The above polyad estimator can be interpreted in relation to a logit classification problem. Specifically, given a polyad $\xi$ and the observed orbit $\mathcal{O}_\xi(Y)$, the minimal rank $m_\xi(Y)$ is distributed according to a conditional Logit model, see McFa_Chap_1974. To see this, consider the transformed graph $\underline{y}=T^{-m_\xi(Y)}(Y)$, which can be thought of as the “minimal” graph in the orbit of the observed graph $y$. In the example of Figure (ref), the graph $\underline{y}$ is the second graph among the six graphs shown on the bottom line (that is for $r=-m_\xi(y)=-1$). The observed orbit $\mathcal{O}_\xi(Y)$ can thus be represented as \[ \mathcal{O}_\xi(Y) := \left\{ \, T_\xi^m(\underline{y}): 0\leq m \leq |\mathcal{O}_\xi(Y)|-1 \,\right\}. \] If the polyad $\xi$ is active, the size of the orbit $|\mathcal{O}_\xi(Y)|= m_\xi(Y)+ M_\xi(Y)+1$ is greater than 2. Changing the indices $m=r+m_\xi(Y)$ in (ref) yields the conditional logit structure for the distribution of $m_\xi(Y)$
The conditional logit distribution of $m_\xi(Y)$ yields the following result.
As mentioned above, the loss function $L(y | X, \beta)$ defined by (ref) is convex in $\beta$. Thanks to Lemma (ref), we can now compute its Hessian \[ \nabla^2_\beta \widehat{L}_\Xi(y | X, \beta) = \sum_{\xi \in \Xi_a} \mathbb{V}_\beta\left[ m_\xi(Y) | Y \in \mathcal{O}_\xi(y) \right] \widetilde{X}_\xi\widetilde{X}_\xi^\top \] and discuss strict convexity.
In other words, $\widehat{L}_\Xi(y|X,\beta)$ is strictly convex in $\beta$ if the DiD-features vary in the data.
In this section we establish consistency and asymptotic normality of the polyads estimator introduced in Section (ref). The concept of an active polyad is central to both our theoretical results and the implementation of the method (see Section (ref)). Recall that a polyad $\xi$ is active when its orbit satisfies $|\mathcal{O}_\xi(Y)| > 1$, meaning it exhibits variation from which the parameters can be identified. This observation motivates normalizing the loss function by the expected number of active polyads and suggests deriving limiting behavior as the number of active polyads increases.
Let $\widehat{N}_a = |\Xi_a|$ denote the observed number of active polyads and let $N_a$ denote its conditional expectation given $X$ under Assumption (ref). Formally,
As mentioned in Remark (ref), $N_a$ plays the role of the average of number of data, we therefore normalize the loss function by $N_a$, defining
In machine learning $\beta\to Q_\Xi(\beta)$ is refered as the risk function.
Recall that $n = \prod_{d=1}^D n_d$ is the size of the graph indexed by $\mathcal{I} = [n_1] \times [n_2] \times \dots \times [n_D]$ and that if $\Xi$ is the set of polyads given by the indices $\mathcal{I}$, then $|\Xi|$ has order $n^2$. Through this section, one is given a sequence $\left(\Xi^{(1)}, \Xi^{(2)}, \dots\right)$ of families of polyads, defined over growing sets of indices. Let $n_d^{(k)}$ be the size of the $d$-th dimension associated with $\Xi^{(k)}$. We impose no restrictions on how each dimension $n_d^{(k)}$ grows, requiring only that $n^{(k)} = \prod_{d=1}^D n_d^{(k)} \to \infty$ --- or, equivalently, $|\Xi^{(k)}|\to\infty$. Notably, our results accommodate short panels where one dimension remains bounded while others diverge. For three-way models, this includes settings where the PPML estimator suffers from the incidental parameter problem weidner2021bias. We also do not assume the existence of a limiting function $Q_\infty$, thus avoiding restrictive assumptions on the asymptotic behavior of fixed effects and covariates. In particular, we do not impose that fixed effects or covariates are drawn from any probability distribution. To circumvent a convoluted notation we use $(k)$ as index in place of $\Xi^{(k)}$, for instance, $\widehat{\beta}_{\Xi^{(k)}}$ is written as $\widehat{\beta}^{(k)}$. Also, notice that associated with each $(k)$ we also have covariates $X^{(k)}$, fixed effects $\theta^{(k)}$ and random variables $Y^{(k)}$.
Finally, unlike graham2017econometric,jochmans2018semiparametric, we do not require the components of $\left[n_d^{(k)}\right]$ to arise from random sampling. This permits a more general interpretation of the data generating process. Consider a two-way model of doctor-patient interactions. Under the framework of graham2017econometric,jochmans2018semiparametric, one assumes the existence of a large population graph containing all doctors and patients, from which patients and doctors are randomly sampled, inducing distributions on both fixed effects and covariates. Our approach is more constructive: for each $(k)$, the distribution of fixed effects and covariates may differ entirely. In the doctor-patient example, our framework accommodates data collection that expands geographically. For instance, initially observing patients and doctors in Paris only, then adding doctors from Marseilles, then adding also patients from Marseilles, and so forth. Imposing a random sampling structure would limit this possibility, as any finite sample could contain doctors from Marseilles with positive probability.
The polyads estimator (ref) is defined as the minimizer of a loss function that is itself the sum of possibly non-independent terms. Our consistency result proceeds in two steps. First, we show that the normalized empirical loss $\widehat{Q}^{(k)}$ converges to its expectation $Q^{(k)}$ as $k\to\infty$. Second, we establish regularity of $Q^{(k)}$ around its minimum $\beta_\star$. Note that at no point do we require the existence of a limit risk function $Q_\infty$.
Assumption (ref) controls the approximation between the empirical and expected versions of $\widehat{Q}^{(k)}$. Observe that $\ell_\xi$ and $\ell_{\xi'}$ are independent as long as $\xi$ and $\xi'$ share no edges. Moreover, a polyad $\xi$ contributes to the loss only when it is active. Thus, Assumption (ref) controls the amount of dependence across the terms of $\widehat{Q}^{(k)}$, enabling a law of large numbers for the normalized losses.
Assumption (ref) is not restrictive: as $n^{(k)}$ grows, the number of positive entries $Y_{\bf i}^{(k)}$ with disjoint indices should also increase, and therefore, the number of pairs of active polyads sharing no edge should outnumber those sharing at least one edge. Notice that if one assumes a distribution on the fixed effects and on the covariates this assumption is automatically true, as in this case the probability that $\xi$ and $\xi'$ are both active is always upper and lower bounded by a constant, thus it is just a matter of comparing the number of pairs of polyads -- which has order $(n^{(k)})^4$ -- and the number of pairs of polyads sharing an edge -- which has order $(n^{(k)})^3$. This last approach is used by graham2017econometric and jochmans2018semiparametric.
Next, Assumption (ref) ensures enough regularity of $Q^{(k)}$ around the minimum $\beta_\star$ so that it can be identified.
We are now in position to state the following consistency result:
Before introducing the assumptions and main result on the asymptotic normality of the polyads estimator we first investigate the error $\beta^{(k)} - \beta_\star$. A Taylor approximation yields, for some $\overline{\beta}^{(k)}$ between $\beta_\star$ and $\widehat{\beta}^{(k)}$, \[ \nabla \widehat{Q}^{(k)}\left(\beta_\star\right) - \nabla \widehat{Q}^{(k)}\left(\widehat{\beta}^{(k)}\right) = \nabla^2 \widehat{Q}^{(k)}\left(\overline{\beta}^{(k)}\right) \left(\beta_\star - \widehat{\beta}^{(k)}\right), \] which implies
If $\widehat{\beta}^{(k)}$ is close to $\beta_\star$, $\widehat{Q}^{(k)}$ is a good approximation of $Q^{(k)}$ and $\nabla^2 Q^{(k)}(\beta_\star) \to \Gamma$ for some invertible $\Gamma$ we may approximate \[ \widehat{\beta}^{(k)} - \beta_\star \approx - \Gamma^{-1} \nabla \widehat{Q}^{(k)}\left(\beta_\star\right).\] The last approximation shows that to control the fluctuations of $\widehat{\beta}^{(k)} - \beta_\star$ we need to control the fluctuations of $\nabla \widehat{Q}^{(k)}\left(\beta_\star\right)$, which is not a sum of independent random variables. Following graham2017econometric, jochmans2018semiparametric, we use Hájek projections to approximate $\nabla \widehat{Q}^{(k)}\left(\beta_\star\right)$ by a sum of independent random variables for which classical central limit theorems hold. The following assumption is useful to control the error of this approximation:
Assumption (ref) bears resemblance to Assumption (ref). Assumption (ref) asks for the active polyads to be distributed among the edges of $Y_{\bf i}$ in a way that no edge ${\bf i}$ has a substantial amount of active polyads $\xi$ such that ${\bf i} \in \mathcal{E}(\xi)$. Meanwhile and up to the weights $w_\xi$, which are zero when $\xi$ is not active, Assumption (ref) looks at the cases where two polyads share one edge and requires that, up to some weights, in most of these cases only one edge is being shared.
Assumption (ref) below is a technical assumption to control the convergence of the Hessian:
In order to present the main result we introduce the quantity $\Sigma^{(k)}$, which is the covariance of the Hájek projection of $\nabla \widehat{Q}^{(k)}\left(\beta_\star\right)$ and is given by
where \[ \bar{s}_{{\bf i}}^{(k)} = \sum_{\xi:{\bf i}\in{\cal E}(\xi)} {\mathbb E}_{\beta^*}\left[\nabla \ell_\xi\left(Y^{(k)} | X^{(k)}, \beta_\star\right) \left| X^{(k)}, Y^{(k)}_{\bf i}\right.\right]. \]
Under Assumption (ref), $\Sigma^{(k)}$ will play the role of the asymptotic variance of the polyads estimator if it does not vanish to zero as $k\to+\infty$ and if a third moment condition holds. We gather these last conditions in Assumption (ref) below:
Assumption (ref) is useful to bound the third moment of the Hájek projection of $\nabla \widehat{Q}^{(k)}\left(\beta_\star\right)$ in a central limit theorem and is not restrictive. It is satisfied as long as the quantities $\bar{s}_{\bf i}$ (and their variances) are not concentrated on a limited number of edges ${\bf i}$ as $k\to\infty$.
In Theorem (ref), the asymptotic normality is written in terms of $\Sigma^{(k)}$, which cannot be directly evaluated from the sample. In Section (ref) we discuss two alternatives to approximate $\Sigma^{(k)}$ from the sample; in Section (ref) we provide empirical validation that the confidence intervals obtained by each approach are accurate.
Lemma (ref) provides the tools to solve the optimization problem (ref) defining the polyads estimator $\widehat{\beta}_\Xi$. Since each term in the loss function is convex and we have expressions for both their gradients and their Hessians, we can efficiently solve the minimization problem using Newton's method. However, computing the gradient and Hessian involves summing over all polyads $\xi \in \Xi$, which naively requires looping over $\prod_{d=1}^D n_d (n_d - 1)$ terms, an operation that quickly becomes computationally expensive. The main goal of this section is to reduce this complexity by avoiding unnecessary iterations, which is done characterizing the set of active polyads $\Xi_a$. Besides that, we also discuss the implementation of two approximations of the variance. These approximations are essential to construct confidence intervals. The methods presented in this section are of special interest when $Y$ is sparse. In particular, the computational complexity of our methods outperforms the PPML alternative correia2019ppmlhdfe when the size of $E = \{ {\bf i} : Y_{\bf i} > 0 \}$ is of order smaller than $\sqrt{n}$.
A key tool that we explore to obtain efficient computational implementations of the polyads estimator is their invariance to permutations. A permutation of a polyad defined by $({\bf i}, {\bf i}')$ is obtained by flipping some indices between ${\bf i}$ and ${\bf i}'$. For instance, \[ \xi'=
is a permutation of \xi=
on indices d\in\{2,3\}. \] We say that a permutation is odd when an odd number of indices are flipped and even when an even number of indices are flipped. Notice that each polyad $\xi$ has a total of $2^D$ unique permutations, including itself. Lemma (ref) below collects the main tools that will be necessary in this section:
A direct conclusion of Lemma (ref) is that if $\xi \in \Xi_a$, then exists a permutation $\xi'$ of $\xi$ such that $i_d < i_d'$ for all $d = 2, \dots, D$ and $m_{\xi'}(y)>0$. If $M_{\xi'} = 0$ then this permutation is unique, but if $M_{\xi'} > 0$ then there are two such permutations, one with $i_1 < i_1'$ and another with $i_1 > i_1'$. Based on this observation we define the following set: \[ \Xi_a^\star = \left\{ \xi = ({\bf i}, {\bf i}') : m_\xi(y) > 0 \text{ and }i_d < i_d' \,\forall d=2,\dots,D \text{ and } (M_\xi(y) = 0 \text{ or } i_1 < i_1') \right\}. \]
The set $\Xi_a^\star$ contains exactly one permutation of each active polyad $\xi \in \Xi_a$ and, by part (iii) of Lemma (ref),
Thus, to solve (ref) it suffices to look at all polyads in $\Xi_a^\star$. The definition of $\Xi_a^\star$ also leads to an efficient method to construct it. Notice that $m_\xi(y)$ it is positive if and only if $y_{\bf i} > 0$ for all ${\bf i}$ with $s_\xi({\bf i}) = 1$. In particular, we need at least $y_{\bf i}$ to be positive to have $m_\xi(y)$ positive. This suggests looping over pairs ${\bf i},{\bf i}' \in E$. The procedure to do it differs slightly depending on the parity of $D$.
First, take $D=2$ and let $i_1 \neq i_1'$ be given. To have $({\bf i}, {\bf i}') \in \Xi_a^\star$ we need to find $i_2 \neq i_2'$ such that $y_{i_1 i_2} \wedge y_{i_1' i_2'} > 0$. Now let $D=3$ and $i_1 \neq i_1'$ be given. We search for $i_2 \neq i_2'$ and $i_3 \neq i_3'$ satisfying $y_{i_1 i_2 i_3} \wedge y_{i_1 i_2' i_3'} \wedge y_{i_1' i_2 i_3'} \wedge y_{i_1' i_2' i_3} > 0$, in particular, $y_{i_1 i_2 i_3} \wedge y_{i_1 i_2' i_3'} > 0$. More generally, given $i_1$ we let \[E_{i_1} = \{ (j_2, \dots, j_D) : y_{i_1j_2\dots j_D} > 0 \},\] it holds that if $({\bf i}, {\bf i}') \in \Xi_a^\star$, then (i) $(i_2, \dots, i_D) \in E_{i_1}$ and $(i_2',\dots, i_D') \in E_{i_1'}$ when $D$ is even; or (ii) $(i_2, \dots, i_D), (i_2', \dots, i_D') \in E_{i_1}$ when $D$ is odd. Thus, one only needs to loop over the pairs $i_1 \neq i_1'$ and, for each of these pairs, over $((i_2, \dots, i_D), (i_2', \dots, i_D')) \in E_{i_1} \times E_{i_1'}$ (if $D$ is even) or $((i_2, \dots, i_D), (i_2', \dots, i_D')) \in E_{i_1} \times E_{i_1}$ (if $D$ is odd). Then, one simply verifies the remaining conditions for $({\bf i}, {\bf i}') \in \Xi_a^\star$. This procedure is summarized in Algorithm (ref).
The next theorem establishes the computational complexity of constructing $\Xi_a^\star$ using Algorithm (ref).
Once the set of polyads $\Xi_a^\star$ is computed, we minimize the loss (ref) using Newton's method. Lemma (ref) provides closed-form expressions for the gradient and Hessian of each $\ell_\xi$, so a Newton step can be computed exactly. When the loss is strictly convex (see Lemma (ref)), Newton's method converges from any initial value $\beta^0$. Algorithm (ref) displays one update step from $\beta^t$ to $\beta^{t+1}$ for $t\geq0$.
The function EvaluateMoments in Algorithm (ref) must return the expectation and variance of $m_\xi(Y)$ conditioned on $Y \in \mathcal{O}_\xi(y)$ when $Y$ has law parametrized by $\beta = \beta^t$. Notice that from (ref), for all $\beta$,
where \[ v(m,\xi;\beta) = m \,\beta^\top \widetilde{X}_\xi \;-\; \sum_{{\bf i} \in \mathcal{E}(\xi)} \ln \bigl(\underline{y}_{\bf i}^{\,m}! \bigr), \] and $\underline{y}_{\bf i}^{\,m} = y_{\bf i}+(m-m_\xi(y))s_\xi({\bf i})$ . Directly evaluating $v(m,\xi;\beta)$ for each $m$ would require computing log factorials, we avoid this computation by noticing that \[ v(m,\xi;\beta) - v(m-1,\xi;\beta) = \beta^\top \widetilde{X}_\xi + \sum_{{\bf i} : s_\xi({\bf i}) = -1} \ln\!\bigl(y_{\bf i} - (m-1-m_\xi(y))\bigr) - \sum_{{\bf i} : s_\xi({\bf i}) = 1} \ln\!\bigl(y_{\bf i} + (m-m_\xi(y))\bigr), \] so the values $v(m,\xi;\beta) - v(0, \xi; \beta)$ can be computed sequentially by cumulative summation without a log factorial.
The procedure EvaluateMoments, given in Algorithm (ref), computes the moments $\mu,\sigma^2$ required in Algorithm (ref) using approximately $2^D |\mathcal{O}_\xi(y)|$ operations. Since this cost scales linearly with the orbit size, evaluating all $Y \in \mathcal{O}_\xi(y)$ may become prohibitive when the orbit is large. In practice, whenever $|\mathcal{O}_\xi(y)|$ exceeds a predefined threshold $L$, we approximate the conditional distribution of $m_\xi(Y)$ by restricting the computation to the truncated set \[ m \in \Bigl[ m_\xi(y) - \bigl( L/2 \wedge m_\xi(y) \bigr), \dots, m_\xi(y) + \bigl(L/2 \wedge M_\xi(y)\bigr) \Bigr]. \] This truncation has negligible numerical effects, because the distribution of $m_\xi(Y)$ is concentrated around $m_\xi(y)$, and extreme values contribute essentially nothing to the expectation or variance.
We now discuss two approaches for evaluating the covariance matrix of the polyads estimator. To shorten notation let $\nabla \ell_\xi\left(Y | X, \widehat{\beta}_\Xi\right)$ be denoted by $\nabla \widehat{\ell}_\xi$ and define \[ \widehat{\Gamma} = \nabla^2 L_\Xi\left(Y | X, \widehat{\beta}_\Xi\right). \] In Theorem (ref), the asymptotic normality is written in terms of $\Sigma^{(k)}$, which can not be directly evaluated from the sample. In practice, (ref) suggest to approximate it by $\widehat{\Sigma}$ given by
We also implement and empirically verify the performance of another variance estimator. Equation (ref) suggests approximating the covariance of $\left( \nabla^2 \widehat{Q}^{(k)}\left(\overline{\beta}^{(k)}\right) \right) \left( \widehat{\beta}^{(k)} - \beta_\star\right)$ by the expectation of \[ \left( \nabla \widehat{Q}^{(k)}\left(\beta_\star\right)\right) \left( \nabla \widehat{Q}^{(k)}\left(\beta_\star\right)\right)^\top = \left( \widehat{N}_a^{(k)} \right)^{-2}\sum_{\xi, \xi' \in \Xi^{(k)}} \left( \nabla \ell_\xi\left(Y^{(k)} | X^{(k)}, \beta_\star\right) \right) \left( \nabla \ell_{\xi'}\left(Y^{(k)} | X^{(k)}, \beta_\star\right) \right)^\top. \] Notice that if $\xi$ and $\xi'$ share no edges, the expectation of their corresponding term is zero since it is the product of independent quantities with zero mean. This suggests approximating $\Sigma^{(k)}$ by
The difference between $\Omega$ and $\Omega'$ is subtle. The next lemma illuminates this difference and provides computationally tractable expressions for $\Omega$ and $\Omega'$.
Indeed, (ref) and (ref) make explicit the two main differences between $\widehat{\Omega}$ and $\widehat{\Omega}'$. First, $\widehat{\Omega}$ contains duplicates of certain pairs of polyads. A closer inspection of the proof reveals that these duplicates arise precisely on pairs that share strictly more than one edge. This clarifies the role of Assumption (ref) in Theorem (ref): for the projection strategy to be valid, the covariance of the Hájek projection of $\nabla \widehat{Q}^{(k)}(\beta_\star)$ --- which is approximately $\widehat{N}_a^{-2}\mathbb{E}\,\widehat{\Omega}$ --- and the true covariance --- approximately $\widehat{N}_a^{-2}\mathbb{E}\,\widehat{\Omega}'$ --- must converge to each other. Second, computing $\widehat{\Omega}$ is less costly than computing $\widehat{\Omega}'$. For each $\mathbf{i} \in \mathcal{I}_a$, evaluating $\widehat{\Omega}$ requires only a single pass over each $\mathbf{i}'$ such that $(\mathbf{i}, \mathbf{i}') \in\Xi_a$. In contrast, computing $\widehat{\Omega}'$ requires an extra loop over ${\bf i}''$ such that $(\mathbf{i}, \mathbf{i}'') \in\Xi_a$.
We now use equations (ref) and (ref) to obtain an algorithm for computing $\widehat{\Sigma}$ and $\widehat{\Sigma}'$. We need to be able to loop over all ${\bf i} \in \mathcal{I}_a$ and, given ${\bf i}$, to efficiently loop over all ${\bf i}'$ such that $({\bf i}, {\bf i}') \in \Xi_a$. Recall that $\Xi_a^\star$ contains exactly one permutation of each active polyad $\xi \in \Xi_a$, in fact, $\widehat{N}_a = |\Xi_a| = 2^D|\Xi_a^\star|$. One can loop over each $\xi \in \Xi_a^\star$ and compute all permutations of $\xi$. By updating a dictionary containing for each key ${\bf i}$ the corresponding set of ${\bf i}'$s we can easily construct the data structure needed to evaluate the covariances. In practice we also store a pointer to the original $\xi$ so that we can profit from the already evaluated $\widehat{X}_\xi$ and $\{Y_{\bf i} : {\bf i} \in \mathcal{E}(\xi)\}$. The following result gives the computational complexity of computing each variance alternative.
As a consequence of this section's discussion, we can provide a complete computational cost analysis of our method:
In practice, Newton's method converges in less than $10$ iterations, yielding $O(|E|^2)$ operations. This quantity is to be compared with the fast implementation of PPML from correia2019ppmlhdfe, which is $O(n)$. Our analysis suggests that our method is faster when $|E| \ll \sqrt{n}$ and competitive when $|E|$ is of order $\sqrt{n}$. Although being the standard practice when reporting the complexity of algorithms, the big-O notation hides a constant that matters to practitioners. In the next section, we empirically verify that, as predicted by our cost analysis, our method outperforms PPML in terms of computational time when $|E|$ is smaller than $\sqrt{n}$ and remains competitive as $|E|$ grows. Indeed, in our computational setup (see Section (ref)) the running time of our method is shorter than that of PPML as long as $|E| \leq 15 \sqrt{n}$. This is, for example, the case of a bipartite network with $n_1=n_2$ such that the average degree of a node is at most $15$.
We provide experiments comparing our method with PPML and with the analytical debias proposed by zylkin2024bootstrap. We consider artificial and real data. With artificial data we investigate the impact of the incidental parameter problem while knowing the correct value of $\beta_\star$. With real data we display evidence of the incidental parameter bias and show how it may lead towards wrong conclusions in inference. We also discuss the computational time and the effect of sparsity for all methods.
Our polyads estimator is implemented as described in Section (ref), and confidence intervals use the covariance approximation $\widehat{\Sigma}'$ from Section (ref). All reported running times correspond to executions of $25$ seeds in parallel on a server equipped with a 32vCPU Intel Xeon Gold 6444Y and 3 TB of RAM.
We conduct computational experiments using a three-way data-generating process inspired by weidner2021bias; see Example (ref). The dimensions are $n_1 = n_2$ (varied) and $n_3 = 5$ (fixed). The fixed effects $u_{ij}, w_{it}, v_{jt}$ are i.i.d.\ $\mathcal{N}(0,1/16)$, and $\beta_\star = 1$. The covariates $X_{ijt}$ are correlated with the fixed effects and, along the third axis, with their own past values: \[ X_{ijt} =
\] The mean of the ${\bf i}=(i,j,t)$ edge's weight satisfies \[ \mathbb{E}(Y_{ijt}) = \exp\!\left(c + \beta_\star X_{ijt} + u_{ij} + w_{it} + v_{jt}\right), \] where the constant $c$ can be selected to control the density $|E|/n$ of the graph. We generate both Poisson data (satisfying Assumption (ref) with intensity $\lambda_{\bf i} = \lambda_{ijt}$) and non-Poisson data. To obtain non-Poisson outcomes, we generate $Y_{{\bf i}}$ as Negative Binomial via a Gamma--Poisson mixture: the rate of the Gamma controls the variance, while its shape is scaled to match the desired mean $\lambda_{{\bf i}}$. Setting the rate to $\infty$ recovers the Poisson model; for experiments with overdispersion, we take the rate equal to $0.1$. We executed $600$ replications of each configuration.
We compare three estimators: PPML, PPML debiased, and our polyads estimator. For PPML we use the fast implementation of correia2019ppmlhdfe; for PPML debiased we use the analytical correction of weidner2021bias via their Stata package. Our first experiment varies the graph density with \[ |E| \in \{0.02n,\ 0.03n,\ 0.04n,\ 0.05n,\ 0.1n\}, \quad\text{and}\quad n_1 = n_2 \in \{50,100\}. \] Each configuration is replicated $600$ times. Figure (ref) summarizes the results. The first row shows the distributions of the normalized errors $\sqrt{n}\,(\widehat{\beta}_\Xi - \beta_\star)$ at densities 2%, 5%, and 10%. The second row reports, from left to right, the empirical coverage of the 95% confidence interval (which should be close to 95%), the convergence rate, and the running time. We observe:
We next consider a second experiment aimed at evaluating performance in the sparse regime. Here we let $n$ vary from $2\times 10^{5}$ to $32\times 10^{5}$ and choose the constant $c$ so that \[ |E| \approx 4\sqrt{n}, \] thus forcing the density $\frac{|E|}{n}$ to shrink towards zero as the sample size grows. This setting is substantially sparser than the fixed-density designs considered above. Figure (ref) reports the results for the Poisson case and Figure (ref) reports the results for the Negative Binomial case. The qualitative conclusions of both cases are similar and consistent with those found in the low-density examples from Figure (ref), but become even more pronounced under sparsity:
Across bias, coverage, computation time, and convergence, the results in this sparse regime reinforce the findings from the low-density experiment: PPML exhibits severe incidental parameter bias and essentially zero inferential validity; the debiased PPML estimator improves upon PPML but continues to suffer from miscoverage under sparsity; the polyads estimator remains accurate, fast, and statistically reliable, even when sparsity increases with sample size.
Comparing Figures (ref) and (ref) we notice that the Negative Binomial case closely match the Poisson one: PPML remains biased, while both the debiased PPML and our polyads estimator remove the bias and achieve better coverage. Thus, although our theory assumes a Poisson model, these results indicate that the polyads estimator is empirically robust to overdispersion and model misspecification.
We study health insurance claims data that cover the universe of physician consultations in metropolitan France over the years 2016 to 2018.\footnote{More details on data are provided in Appendix (ref).} In May 2017, the fees charged by general practitioners (GPs) belonging to the regulated sector\footnote{The majority of French GPs are subject to fee regulation. In January 2016, only 8.6% of GPs can charge fees in excess of the regulated level.} have increased by 8.7%. Using a difference-in-differences approach at the doctor level, we find that the stronger financial incentives have caused physician activity (as measured by number of visits) to rise by approximately 10% (see Appendix (ref) for details about the considered control groups and specifications, as well as Table (ref)).
Given the strong policy concern about care accessibility, it is important to understand the impact of the reform on potential patients. Have certain categories of patients (patients with chronic diseases, women, elderly patients, patients living in medically undeserved areas, low-income patients, etc.) been disproportionately affected by the reform? Has the reform caused patients to seek care further away from their home? To answer these questions, we need to perform the analysis at the doctor-patient level rather than at the doctor level. Using a three-way model and controlling for dyads fixed effects allows to estimate how the reform has changed the doctor-patient connections. To illustrate the method, we highlight in this subsection the role of two variables: gender and geographic distance (spatial accessibility).
The outcome $Y_{ijt}$ is the number of visits by patient $i$ to doctor $j$ on month $t$. We aggregate data at the city-sex level. The high number of potential patients in each city-sex group makes to the Poisson assumption plausible.\footnote{The aggregate number of consultations in each group is the sum of (possibly heterogeneous) Bernoulli distributions that represent the occurrence of a consultation for all potential patients in the group. This sum follows approximately a Poisson distribution under the conditions exposed in le1960approximation. The assumption that the individual occurrences of a consultation are independent across potential patients is relaxed in galambos1973general and serfling1978some. } The index $i$ (resp. $j$) stands thus for the set of patients (resp. doctors) in a given municipality with given gender. The treatment $T_{j}$ is a binary variable equal to 1 for sector 1 GPs, and to 0 for direct access specialists. The reform has been implemented from May 2017 onward, hence the definition of $\text{Post}_t$, a dummy variable equal to 1 after that date. The model includes three features that account for any possible policy-induced change in homophily preferences as regards gender and spatial dimensions: the interactions between $\text{Post}_t \times T_j$ (as in any difference-in-differences approach) and (i) a dummy variable equals to 1 when patients and doctors have the same sex, (ii) a dummy variable that is equal to 1 when the medical office chosen is located in the same municipality as the patient's home, and (iii) travel time (measured in minutes between the centroids of municipalities).
We thus consider the three-way Poisson model $Y_{ijt}\sim\mathcal{P}(\lambda_{ijt})$ with intensity given by \[\ln\lambda_{ijt}=\big(\beta_\texttt{d} d_{ij}+\beta_\texttt{sc} \mathds{1}\{\text{city}_i=\text{city}_j\}+\beta_\texttt{ss} \mathds{1}\{\text{sex}_i=\text{sex}_j\}\big) \times \text{Post}_t \times T_j +u_{ij}+v_{jt}+w_{it}.\]
We compare three estimators: the standard PPML estimator, the analytical correction of weidner2021bias, and our polyads estimator. Because the full dataset contains $n_1 = 69{,}265$ patient groups, $n_2 = 16{,}941$ doctors, and $n_3 = 34$ months --- corresponding to roughly $n \approx 40$ billion edges and $|E| = 56{,}034{,}015$ --- direct estimation on the full graph is computationally infeasible for all methods. We therefore adopt a subsampling strategy combined with meta-analysis.
For each configuration, we sample a proportion $s \in \{2\%, 3\%, 4\%\}$ of patient groups and the same proportion of doctors, while always retaining all $34$ months\footnote{Notice that we sample the groups of patients, the doctors and then get all edges through the sampled nodes, keeping all times. This sampling strategy is consistent with Assumption (ref) and also with the sampling assumptions in fernandez2016individual, weidner2021bias, graham2017econometric, jochmans2018semiparametric. One could also propose sampling directly the edges and not the nodes, but there is evidence that this procedure can lead to bias, see subsamplegraphs for details.}. Each subsample consists of independent random draws of patient and doctor groups; doctors in the treatment and control groups are sampled independently to ensure comparability. We repeat the procedure independently to obtain $100$ subsamples. The final estimates are obtained via a random-effects meta-analysis using the default implementation in the statsmodels Python library, which is based on the iterated method of paule1982consensus, a refinement of the classical method of dersimonian1986meta.
The results are the following. First, Table (ref) and Figure (ref) show that debiased PPML has many outliers, mostly on `small' samples (admittedly, not that small since $|E|\approx 20,000$ with the 2% sampling), which results in (very) large confidence intervals and point estimates. As expected, the variance decreases with the subsample size, regardless of the estimation method used.
Second, as reported on Table (ref), the running time of PPML and debiased PPML scales with $n$ while the polyads time scales with $|E|^2$.
Third, a closer look at the 4% subsample coefficients suggests that: (i) PPML estimates of $\beta_d$ yields confidence intervals that do not contain the point estimate of debiased PPML and of Polyads, which may be a sign of incidental parameter bias. This is to be compared with the artificial experiments where PPML yields small variance around a biased point, missing the point estimate. All methods point towards a positive coefficient, implying that after the reform some patients travel further away from home to get healthcare; (ii) as regards the preference coefficient for being treated in the same municipality $\beta_{sc}$, the Polyads estimator yields similar results to the PPML, while debiased PPML gives a different point estimate, although the latter confidence intervals are overall consistent with PPML and Polyads. All methods point towards a negative coefficient, meaning that after the reform patients are more likely to visit a doctor outside their own city ; (iii) as regards gender homophily, we see some difference up to the level of signs. The PPML would lead to estimate a negative impact on the reform on that homophily (this diagnosis holding true regardless of the subsample size considered here), while Polyads and debiased PPML both contain $0$ in their confidence intervals in the $4\%$ sample. Note also that, PPML yields a confidence interval that is close to the lower half of Polyads confidence interval, this may be a sign of the incidental parameter problem as PPML may be giving a narrow variance around a biased point estimate and leading to possibly wrong conclusions (in this case, that gender homophily between patients and doctors decreased).
A fast python-based implementation of the polyads method here presented is available at the Github repository lucasresenderc/polyads, available at \href{https://github.com/lucasresenderc/polyads}{https://github.com/lucasresenderc/polyads}.
We are grateful to Cl\'ement de Chaisemartin, Laurent Davezies, Xavier D'Haultf\oe uille, Francis Kramarz, and Yannick Guyonvarch for helpful comments. We thank the Agence Nationale de la Recherche for financial support (ANR-23-CE36-0014).
\numberwithin{equation}{section}
\setcounter{equation}{0} \setcounter{table}{0} \setcounter{figure}{0}
\numberwithin{figure}{section} \numberwithin{lemma}{section} \numberwithin{table}{section}