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.
100,771 characters · 25 sections · 67 citation commands
On the use of U-statistics for linear dyadic interaction models
\bibpunct{(}{)}{,}{a}{,}
Dyadic regression analysis is a common practice in several applications for network models. It is used, for instance, in the estimation of gravity models for international trade flows since its establishment by tinbergen1962shaping. As defined by graham2020dyadic, a dyadic dataset corresponds to a situation where the outcome of interest reflects a pairwise interaction among the sample units. Therefore, it is natural that datasets on trade flows are characterized by such a dyadic structure, as the value of imports and exports are determined by both the importer and the exporter countries. Other examples of applications of dyadic settings are, for instance, the estimation of models of migration, equity, international financial flows, and information flows jackson2013diffusion.
Even though dyadic regressions are widely used in empirical applications, the (asymptotic) properties of estimation methods only began to be studied recently in the literature. One key feature that all studies mentioned above contain is the presence of two-way unit-specific effects, one for each individual in the dyadic interaction, and the fact that the outcome variable (and, in most cases, explanatory variables) is double indexed. For linear models with the two-way unobserved heterogeneity and the idiosyncratic error term entering additively in the specification it is possible to estimate the model consistently and (asymptotically) unbiasedly using the two-way fixed effects estimator (as long as the model is correcly specified; see juodis2020shock). .
However, many economic models are non-linear, more specifically, in the context of network datasets, network formation models have the structure of discrete choice models (where the outcome of interest is binary), and outcomes of interest that are bounded below at zero (such as trade flows) can be approximated by a gravity equation in its multiplicative form. Naturally, one approach be to estimate such models with a probit/logit, and a poisson-pseudo maximum likelihood estimator silva2006log, respectively. The challenge is that, while it is desirable to treat the unit-specific effects as parameters to be estimated (i.e., fixed effects, such that the conditional distribution of the unobserved heterogeneity given the covariates is left unrestricted), as fernandez2016individual show, even if both dimensions of the (pseudo-)panel dataset grow with the sample size, these estimators suffer from the incidental parameter problem neyman1948consistent in the presence of two-way fixed effects. This problem occurs since, in non-linear models, the estimates of the coefficients of the covariates depend on the estimates of the fixed-effects, and the latter converges at a slower rate than the first, resulting in an asymptotic bias in the estimates (and, therefore, invalid inference).
To address the incidental parameter problem in estimates for non-linear models with two-way fixed effects, such as logit, probit and poisson-pseudo maximum likelihood estimators, fernandez2016individual proposed analytical and jackknife bias correction methods. However, for network formation models (discrete choice models), others propose a conditional maximum-likelihood approach under the logistic specification, such as in charbonneau2017multiple, jochmans2018semiparametric and graham2019network. The conditioning sets in this approach translates to a model where pairwise differences of the outcomes (and covariates) are taken such that the two-way fixed effects are differenced out from the objective function, eliminating the incidental parameter problem.
The advantage of the conditional maximum-likelihood approach as opposed to the bias-correction methods is that it accomodates sparse networks in the case of a network formation model jochmans2018semiparametric. However, when taking such differences in the model, in general, the summands of the influence functions will not be independent anymore, showing some dependence on the unit level and translating to the fact that usual law of large numbers and central limit theorems do not straightforwardly apply.
To overcome such obstacle, U-statistics tools are generally applied such that one can show the asymptotic properties of those estimators. Although the U-statistics properties are well-known for single-index variables (see serfling2009approximation, and van2000asymptotic), those of double-indexed variables are generally not treated (up to my knowledge) in textbooks, with few exceptions related to applications for dyadic contexts, such as graham2019network.
The main purpose of this paper is to illustrate, step-by-step in a comprehensive way, how to obtain the asymptotic properties of pairwise differences estimators (such that the fixed effects are cancelled out) for models of dyadic interactions using tools from the literature on U-statistics. More specifically, we show how to accomodate such tools to double-indexed variables, as the outcome is and the covariates are indexed by both individuals in the interaction. Even though, as mentioned earlier, the classical two-way fixed effects estimators delivers consistent and (asymptotically) unbiased estimates in the linear model (such that the pairwise differencing approach is not needed), for simplicity, we consider a linear two-way fixed effects, but the arguments can be generalized (with additional regularity conditions) to non-linear models or models with multiplicative individual heterogeneity and errors structure (as showed in jochmans2017two).
A U-statistic for single-indexed variables is formally defined as an unbiased estimator that is of the form of an average of a function (kernel) of i.i.d. random variables. The main idea to determine the asymptotic properties of this estimator is to define a projection of the U-statistic, the so-called H\'{a}jek projection, that is asymptotically equivalent to the U-statistic itself. This projection consists on the average of the conditional expected value of the U-statistics, where each summand is the expected value of the U-statistic conditional on each index. Thus, by conditional independence arguments the H\'{a}jek projection becomes a simple average of i.i.d. random variables, to which laws of large numbers and central limit theorems can be applied. This concept will become more clear in the following sections.
A key result provided in this paper, is that, for a directed dyadic structure, there can be different ways of defining the H\'{a}jek projection, which will lead to distinct, but equivalent, consistent estimators for the asymptotic variance of the proposed estimator. More specifically, we provide two possible projections depending on which random variables one conditions the expected value of the U-statistics. Central to both possibilities is the fact that the summands of the influence function have a conditional independence structure once conditioned on dyad level attributes and the individual heterogeneity of both individuals forming a dyad (a fundamental difference with respect to single-index contexts). This is intuitive to dyadic settings, where the dependence across the dyads arises only through the individual fixed effects and the possible correlated observations for the same individual in different pairwise interactions and thus, outcomes and covariates.
The organization of this paper is as follows: in Section 2 we define the linear model of directed dyadic interactions, in Section 3 we propose a pairwise differences estimator, in Section 4 we explain how some tools of U-statistics can be employed in this context and extended to dyadic settings, in Section 5 we discuss the asymptotic properties of the estimator and possible consistent estimates of its asymptotic variance and in Section 6 we demonstrate a Monte Carlo simulation exercise to investigate the finite sample properties of the estimator.\\ \\ Notation
Random variables are denoted by capital letters, specific realizations thereof by lower case, and their support by blackboard bold capital letters. That is, $Y$, $y$, and $\mathbb{Y}$ respectively denote a generic draw of, a specific value of, and the support of $Y$.
Calligraphic letters denote sets. For instance, denote by $\mathcal{N} = \{1,2,...N\}$ the set of indices for $N$ individuals (or nodes). Denote by $\mathcal{C}(\mathcal{N},4)$ the multiset containing all sets of combinations of four individuals from the sampled $N$ observations. Moreover, denote $\rvert \mathcal{C}(\mathcal{N},4) \rvert = {N \choose 4}$ the number of obtained combinations, and denote by $\mathcal{C}$ an unordered set formed by a given combination, say $\mathcal{C} = \{i,j,k,l\}$.
Set $\mathcal{P}(\mathcal{C},4)$ to be the multiset containing all sets of permutations containing four elements of a given combination $\mathcal{C}$. Also, let $\rvert \mathcal{P}(\mathcal{C},4) \rvert = 4!$ to be the number of possible permutations, and $\pi$ to be the ordered set formed by a given permutation, where $\pi_1$, $\pi_2$, $\pi_3$ and $\pi_4$ denotes its first, second, third and fourth elements. For instance, given a permutation $\pi = \{k,l,j,i\}$, we have that $\pi_1 = k$, $\pi_2 = l$, $\pi_3 = j$ and $\pi_4 = i$.
We consider a linear model of dyadic interactions between $N$ agents, in which we assume that all variables of all pairwise interactions are observed (therefore, there is no sample selection). Let $(y_{ij}, x_{ij})$ denote the realizations of the random vector of outcomes and covariates $(Y_{ij}, X_{ij})$ for the dyad $(i,j)$, i.e., related to the interaction between agents $i$ and $j$. Importantly, $Y_{ij}$ is an outcome variable generated by the interaction of the individuals and it is continuous in this setting. We allow for directed interactions, such that $(Y_{ij}, X_{ij})$ need not be equal to $(Y_{ji}, X_{ji})$, and we do not include self links. Therefore, for a set $\mathcal{N} = \{1,2,... N\}$ of $N$ agents, we have $N(N-1)$ observed dyads. Following the notation of graham2020dyadic, we denote that the first subscript on $Y_{ij}$ or $X_{ij}$ to be the ego, or sending agent, and the second to be the alter, or receiving agent.
Consider the following linear model of dyadic direct interactions taking into account two-way fixed effects:
For simplicity, we consider for now only one regressor $X_{ij}$ (which can easily be relaxed to a vector). We assume that an agent-level attribute $A_i$ (which also can be relaxed to be a vector), and an attribute $B_j$ are observed, such that $X_{ij} = f(A_i, B_j)$ is a constructed dyad-level attribute. On the other hand, the sequences of individual-level heterogeneity $\{\theta_i\}_{i=1}^N$ and $\{\xi_i\}_{i=1}^N$ (for the ego and for the alter, respectively) are unobserved and we treat them as fixed-effects. In particular, there are no restrictions on correlations between $\theta_i$, $\xi_j$ and $X_{ij}$. In other words, the joint distribution between the observed and unobserved agent-level characteristics, $\{A_i, B_i, \theta_i, \xi_i \}$ is left unrestricted, such that the model is semiparametric. Finally, $U_{ij}$ is an idiosyncratic component that is also not necessarily equal to $U_{ji}$.
Taking for instance the classical gravity model for international trade flows given by anderson2003gravity, and usually applied in the empirical literature, the variable $Y_{ij}$ would refer to the log of the value of exports from country $i$ to country $j$, $X_{ij}$ refers to characteristics of the dyad $(i,j)$, for instance, the distance between the two countries, and $\theta_i$ and $\xi_j$ refers to the so-called unobserved multilateral resistance terms. The latter terms refer, for example, to unmodeled export orientation of an economy, undervalued currencies and consumption taste.
We impose the following assumptions on this model:
Given the presence of the two-way fixed effects, namely $\theta_i$ and $\xi_j$, we have that conditional independence between the outputs of different dyads given the sequences of $A_i$ and $B_j$ (or, given the covariates $X_{ij}$) is unlikely to hold. Even conditioning on the sequence of covariates, outcomes that share the same ego or alter indices are likely to not be independent. For instance, the outcomes $Y_{12}$ and $Y_{34}$ are independent of each other, but the outcomes $Y_{12}$ and $Y_{13}$ are likely to be dependent, even after conditioning on $X_{12}$ and $X_{13}$. As pointed out by graham2020dyadic, in the international trade example, this translates to the fact that exports from Japan to Korea will likely covary with exports from Japan to the United States, even after controlling for covariates, due to the Japan exporter effect. graham2020dyadic denotes these patterns as dyadic dependence.
However, after conditioning also on the fixed-effects, that is, conditional on $\{ X_{12}, X_{13}, \theta_1, \xi_2, \xi_3 \}$, or, equivalently from Assumption (ref), conditional on $\{A_1, B_2, B_3, \theta_1, \xi_2, \xi_3\}$, the outcomes $Y_{12}$ and $Y_{13}$ are independent. This result follows from Assumption (ref). This conditional independence structure will be essential for the asymptotic properties of the estimator proposed in the following Section, mainly because this structure is well-suited for applying the tools of U-statistics. shalizi2016cid denotes such models of dyadic interactions with such independency structure as conditionally independent dyad models (CID).
Even though the model given by Equation (ref) and under Assumptions (ref) and (ref) could be consistently and (asymptotically) unbiasedly estimated with a two-way fixed-effects estimator, we propose an estimator that differences out the fixed effects through pairwise differences. This estimator builds up on differencing arguments for a similar model, however, non-linear, introduced by charbonneau2017multiple. She considers a model of network formation, where the outcome variable $Y_{ij}$ is binary, indicating whether an individual $i$ forms a directed link with individual $j$.
As mentioned before, fernandez2016individual shows that maximum likelihood estimators for nonlinear models with two-way fixed effects, such as probit/logit, suffer from the incidental parameter problem even if both dimensions of the (pseudo-)panel dataset tend to infinity. This is due to the fact that the dimensions of the vectors of nuisance parameters (how the fixed effects are treated in both charbonneau2017multiple and in this paper) grows with the number of observations. At this point, it is important to notice that datasets of dyadic interactions can be seen as a pseudo panel data, where both dimensions of the panel tend to infinity as the number of individuals grow. fernandez2016individual proposes analytical bias corrections to reduce the incidental parameter (asymptotic) bias, which was implemented by dzemski2019empirical to a network formation context. However, as explained by jochmans2018semiparametric, the problem with the bias correction approach is that for sparse networks the individual-specific parameters (the fixed effects) may not be consistently estimable or may be estimable only at a very slow rate.
The approach proposed by charbonneau2017multiple becomes very attractive for sparse networks, since, through a conditional maximum likelihood approach for logistic models, it delivers an estimator that differences out the fixed-effects. The estimator essentially is based on a set of conditions that translates to a transformation of the dependent and covariates, where pairwise differences are taken, such that the fixed-effects in the model are cancelled out. Even though a classic logit estimation can be used to obtain the estimates of the coefficients of the covariates, inference does not follow the textbook usual procedures, since agent-level dependencies arise when taking such pairwise differences. The asymptotic properties of this estimator are studied by jochmans2018semiparametric and are obtained by employing tools of U-statistics. He shows that this estimator is consistent, asymptotically unbiased and the estimated variances deliver correct sizes for the t-test.
The estimator that we introduce in this Section is based on a similar pairwise differences methodology for transforming the dependent variable and covariates to difference out fixed-effects as presented by charbonneau2017multiple, however, for a linear model. Our purpose when introducing this estimator is to provide a better understanding, through a simpler and linear model, on how to apply the tools of U-statistics to derive the asymptotic properties of estimators based on such pairwise differences for dyadic data.
First, we define the following notation for the random variable obtained by taking the specified pairwise differences among different dyads' outputs:
and analogously for $\Tilde{X}_{ijkl}$ and $\Tilde{U}_{ijkl}$.
If we substitute the expressions for each of the outcomes $Y_{ij}$, $Y_{ik}$, $Y_{lj}$ and $Y_{lk}$ given by the model in Equation (ref) to the expression for $\Tilde{Y}_{ijkl}$ in Equation (ref), we obtain:
where the fixed effects are differenced out. The equation above is simply a linear regression with the transformed variables obtained by taking the pairwise differences between the dyads $(i,j)$, $(i,k)$ and $(l,j)$, $(l,k)$.
This form of differencing out the fixed effects depends heavily on the fact that the individual-specific heterogeneity parameters (i.e., the fixed effects themselves) enter the model additively. For more general specifications, this transformation fails to difference out the fixed effects. However, other studies, such as chen2021nonlinear and jochmans2017two, study cases with interactive fixed-effects. The former proposes an analytical bias correction estimator and the latter provides also an argument for differencing out the individual-specific parameters.
Inspired by the same methodology of charbonneau2017multiple that is further studied by jochmans2018semiparametric, we can then estimate $\beta_1$ with an ordinary least squares estimator by taking into account the transformed variables $\Tilde{Y}_{ijkl}$ and $\Tilde{X}_{ijkl}$. Notice that the model given by Equation (ref), where the fixed effects are differenced out, holds for all combinations of quadruples of indices from the set $\mathcal{N} = \{1,\dots,N\}$ and its permutations. Therefore, we can write the pairwise differences OLS estimator as:
where, in the second line, we use the fact that summing over all possible permutations of quadruples is equivalent to summing over all possible combinations of quadruples and then all permutations of such combinations. Therefore, say we look at a specific combination given by $ \mathcal{C}$, then the multiset denoted by $\mathcal{P}(\mathcal{C},4)$ corresponds to all permutations of those indices. Then, given a permutation, $\pi = \{i,j,k,l\}$, we have that $\pi_1 = i$, $\pi_2 = j$, $\pi_3 = k$ and $\pi_4 = l$, such that $\pi_1$ refers to the index occuping the first position in the permutation set, and analogously for $\pi_2$, $\pi_3$ and $\pi_4$.
To obtain the properties of the estimator $\hat{\beta}_{1,PD}$, it is useful, as shown in the regular textbook case for OLS estimators to rewrite the previous expression in terms of its influence function:
In order to derive the asymptotic properties of this estimator, it is necessary to first derive the asymptotic properties of the last term in the equation above, namely: $$\left[ \frac{1}{{N \choose 4}} \sum_{\mathcal{C} \in \mathcal{C}(\mathcal{N},4)} \frac{1}{4!} \sum_{\pi \in \mathcal{P}(\mathcal{C},4)} \Tilde{X}_{\pi_1\pi_2\pi_3\pi_4} \Tilde{U}_{\pi_1\pi_2\pi_3\pi_4} \right]. $$ Notice that the transformed error terms $\Tilde{U}_{\pi_1\pi_2\pi_3\pi_4} \equiv (U_{\pi_1\pi_2} - U_{\pi_1\pi_3}) - (U_{\pi_4\pi_2} - U_{\pi_4\pi_3})$ are not independent over the dataset obtained when applying the transformation over all possible combinations and its permutations of quadruples, since the same dyads will appear in different terms, leading to a correlation amongst the terms. Therefore, the traditional application of LLNs and CLTs does not hold straightforwardly.
From now on we will denote the last term in Equation (ref) by:
Even though this term resembles a U-statistic, it is not strictly speaking. However, we can adapt tools used in the literature of U-statistics to obtain the properties of the term $U_N$. Namely, we employ a Hoeffding decomposition hoeffding1948central to obtain the variance of the term $U_N$, and we also propose two possibilities of H\'{a}jek projections of this term, which we prove both to be asymptotically equivalent to $U_N$. The reason why obtaining such projections is that the terms on it are i.i.d. such that it is possible to apply laws of large numbers and central limit theorems to obtain the asymptotic properties of $U_N$, and, thus of the proposed estimator.
In this paper we consider asymptotics under one single network growing, i.e., we consider that $N$ (the number of individuals in a network) tends to infinity when obtaining the asymptotic properties of the proposed estimator.
According to serfling2009approximation, the U-statistic is a generalization of the sample mean, i.e., a generalization of the notion of forming an average. The formal definition of the U-statistic is the following:
We can see that the term $U_N$, as defined in Equation (ref), contains elements of a U-statistic, resembling one at a first glace. However, it is not formally one given the definition above.
The shared properties to a U-statistic are related to having a similar dependence structure, such that it consists of a sum over all combinations of quadruples of individuals, evaluated at some given function, analogous to a fourth-order U-process. We can define the symmetric kernel for a given combination $\mathcal{C} = \{i,j,k,l\}$ in our case to be:
which is essentially the score of our estimator (being the reason why we denote by $s$, and not $h$). Note once again that the indices $k_1, k_2, k_3, k_4$ denote the elements of the permutations of a given combination of individuals ${i,j,k,l}$. Then, we can also see that another shared property with the U-statistic is that the kernel is permutation invariant and that the arguments of it, namely, the random variables $X_{ij}$ and $U_{ij}$, are identically distributed from Assumptions (ref) and (ref).
Another important property that the term $U_N$ has in common to a U-statistic is that, if we define a parametric function $\theta$ to be:
then, we also have in our context that $U_N$ is an estimator of $\theta$, and it is also unbiased, since,
where the second equality follows from linearity of expectations.
In spite of these similarities, the statistic $U_N$ is not an U-statistic as conventionally defined, since its kernel includes random variables at both the individual (since $X_{ij} = f(A_i, B_j)$) and dyad level ($U_{ij}$). Therefore, single-index U-statistics as the one defined in the definition above are not well-suited, hence, the tools need to be slightly modified to accomodate the dyadic structure.
Even more crucial is the fact that the observations $\{X_{ij}\}_{i=1,j\neq i}^N$ are not independent, due to the common individual characteristics $A_i$ or $B_j$. However, the fact that $U_{ij}$ is i.i.d. and that it is independent of $X_{ij}$ allows us to employ tools of U-statistics, such as the Hoeffding decomposition to obtain the variance of $U_N$, and the possibility to define a H\'{a}jek projection to obtain the asymptotic properties of this term (since we will demonstrate that $U_N$ and the projections are asymptotically equivalent). Importantly, to apply LLNs and CLTs to the H\'{a}jek projection, we will exploit the conditional independence structure of the projection. The conditional independence arguments extend straightforwardly to CID models, making them well suited for the use of U-statistic tools.
To derive the variance, we first use some arguments provided by serfling2009approximation, that are also employed by, for instance, graham2017econometric. First, we define:
Notice that independently of which pairs of combinations of quadruples we look at from the sampled individuals, the covariance between the two kernels evaluated at such combinations will only depend on the number of common individuals that the combinations share, namely, $q$. This follows from the fact that from Assumptions (ref) and (ref), $U_{ij}$, $A_i$ and $B_j$ are i.i.d. and that the kernel (score) $s$ is symmetric on its arguments. By working out further the expression for the covariance $\Delta_q$, one can see that the nonzero terms in the expression are mainly driven by the covariance between the idiosyncratic errors. This is due to: (i) $\{U_{ij}\}_{i=1,j \neq i}^N$ being independent of $\{X_{ij}\}_{i=1,j \neq i}^N$, and to (ii) the idiosyncratic errors being independent of each other, while, for instance, $X_{ij}$ is correlated with $X_{ik}$ due to the common individual factor $A_i$. Therefore, if there is no common dyad in the expressions of the kernels for both combinations, the covariance between them will be zero, since $U_{ij}$ is i.i.d. (importantly, for instance, $U_{ij}$ and $U_{ik}$ are independent). This argument will become clearer in the Appendix A.
Due to the dyadic structure and the fact that $U_{ij}$ is i.i.d., $\text{Cov} (s_{ijkl}, s_{mnop}) = 0$ whenever the quadruples share zero or only one individual in common. Therefore, $\Delta_0 = \Delta_1 = 0$ indicates that $U_N$ exhibits degeneracy of order one. As long as the combinations have two or more indices in common, since the kernels sums over all permutations of the combinations, the same idiosyncratic error (with the same indices $i$ and $j$) appears in both terms $s_{ijkl}$ and $s_{mnop}$, leading to a non-zero covariance.
Assuming further that:
Since the covariances $\Delta_q$ are constant across pairs of combinations sharing $q$ individuals in common, we can obtain the variance of $U_N$ through the Hoeffding decomposition hoeffding1948central, as the following Lemma states:
Proof. Provided in Appendix B.1.
Given the result provided by Lemma (ref), we can rescale the statistic $U_N$ by the factor $\sqrt{N(N-1)}$, and by taking into account that $\Delta_0 = \Delta_1 = 0$, and denoting: $$ \bar{s}_{ij} = \mathbb{E}[s_{ijkl} \rvert A_i, B_j, U_{ij}] \quad \text{and} \quad \bar{s}_{ji} = \mathbb{E}[s_{ijkl} \rvert A_j, B_i, U_{ji}], $$ $$ \delta_2 = \mathbb{E} [\bar{s}_{ij}^2] = \mathbb{E} [\bar{s}_{ji}^2], $$ we arrive at the following result:
Proof. Provided in Appendix B.2.
Where the terms of order $\mathcal{O}(1)$ in the above expression relates to the term $\Delta_2$, $\mathcal{O}\left(\frac{1}{N}\right)$ relates to the term $\Delta_3$ and $\mathcal{O}\left(\frac{1}{N^2}\right)$ relates to the term $\Delta_4$. Furthermore, given our simplified model, it is possible to further pin down the expression for $\Delta_2$. This result can be found in Appendix C.
As explained by serfling2009approximation, the appealing feature of a U-statistic as given by Definition (ref), is its simple structure as a sum of identically distributed random variables. But, even in the simpler context of a single-index U-statistic, if the kernel $h$ has a dimension $m>1$, then the summands in the statistic $U_N$ are not all independent, as the sample sampled observations are taken into account in different combinations. Therefore, it is not possible to directly employ LLNs and CLTs for sums of independent random variables, as it is customary done. However, serfling2009approximation and other textbooks on U-statistics show that it is possible to obtain a projection to which the U-statistic can be approximated to. The advantage is that such projection is a sum of i.i.d. random variables, to which classical limit theory can be applied.
In the following we will present the formal definition of this projection, the H\'{a}jek projection, and explain how this concept can be applied in our context. We also highlight that there are considerable differences between our approach and the classical textbook projection. Again, the main difference is that, while the standard definitions account for single-index variables, in our case of a dyadic setting the random variables forming the U-statistics have double-indices. Moreover, the pairwise differences structure in the kernel are formed by random variables reflecting the dyadic interactions originated by four individuals.
Therefore, to obtain the H\'{a}jek projection, instead of conditioning on a single-indexed random variable alone as it is done in textbooks, we condition on both the individual and dyad-level random variables given by the dyad indices. We will show that by doing so, we will still obtain a projection where the summands are conditionally independent, even if the sequence $\{X_{ij}\}_{i=1,j\neq i}^N$ is not formed by independent variables, since the idiosyncratic errors $U_{ij}$ are i.i.d. and independent of the former sequence. This relies on the previously mentioned arguments of conditional independence of CID models.
Besides, in the general textbook case or single-index variables, it is stated that the projection has no purpose in the case where $\Delta_1 = 0$, however, we will see that in our case, due to the dyadic structure, the projection proves to be useful even when $\Delta_1 = 0$ holds.
The most important result in this section is that we can propose two different forms of H\'{a}jek projections, depending whether we condition on all random variables generated by the combination $\{i,j\}$ of a dyad, or if we condition on the random variables generated by the permutation $\{i,j\}$, where the ordering of the indices matter.
According to serfling2009approximation, and following the same notation and framework as in Definition (ref), we have the following definition for a H\'{a}jek projection:
It is important to notice that, in the definition of $\hat{U}_n$ above, when taking the expectation of $U_n$ conditioning on each different $W_i$ for each summand, we are left with a sum of i.i.d. random variables, since $W_i$ are i.i.d. themselves.
We denote the first proposed H\'{a}jek projection by $\hat{U}_{N,1}$. In our context, we already derived before that $\theta = \mathbbm{E}[s_{ijkl}] = 0$ (see Equation (ref)), therefore, we only need to focus now on deriving the first term of the similar projection proposed in Definition (ref). In addition, we are working in a context of dyads, therefore, the sum is over the expected value of the statistic $U_N$ conditional on each of the dyad characteristics, namely, for a given dyad $\{i',j'\}$ we condition on $\{A_i', B_j', U_{i'j'}\}$. Therefore, we sum over all the possible dyads $N(N-1)$. Notice that the order of the indices in the dyad matter, since we have a directed network.
The main idea behind this projection is that the double sum $\sum_{i' = 1}^N \sum_{j' \neq i'}$ fixes the two indices of a dyad, and refers to the individual-level $\{A_{i'}, B_{j'}\}$ and the dyad level characteristics $U_{i'j'}$, which we condition the statistic $U_N$ on. In this case, the order of the indices $(i',j')$ matter to determine on which random variables we condition on.
For each summand of the double sum $\sum_{i' = 1}^N \sum_{j' \neq i'}$ we take the expectation of the statistic $U_N$ conditional on the variables described above. The statistic is essentially an average of the scores $s_{ijkl}$ evaluated at all possible combinations of quadruples $\{i,j,k,l\}$ from the set $\mathcal{N}$. From Assumption (ref), and more precisely, the fact that $U_{ij}$ is independent from the sequence $\{X_{ij}\}_{i=1, j \neq i}^N$, leads to the fact that the only non-zero summands are the ones where the combination $\mathcal{C}$ contains the elements $i'$ and $j'$, and any other two remaining elements. Since the kernel contains all permutations of the combination, inevitably the term $U_{i'j'}$ will appear in the expression for the kernel $s_{ijkl}$ in this case (where $\{i',j'\} \subset \{i,j,k,l\}$, with, for instance, $i=i'$ and $j=j'$), leading to a non-zero conditional expectation.
We can further boil down the expression of the projection $\hat{U}_{N,1}$ by first noting that, as shown in Appendix C, for a given combination $\{i',j',k,l\}$ for any value of $k$ and $l$, the conditional expected value of the kernel evaluated at such combination is of the form: $$ \mathbbm{E}[s_{i'j'kl} \rvert A_{i'}, B_{j'}, U_{i'j'}] = 8 [(X_{i'j'} - \mathbbm{E}[X_{i'j'} \rvert A_{i'}] - \mathbbm{E}[X_{i'j'} \rvert B_{j'}] + \mathbbm{E}[X_{i'j'}])U_{i'j'}]. $$ Moreover, there will be ${N-2 \choose 2}$ possible combinations of four elements of the set $\mathcal{N}$ containing the individuals $i'$ and $j'$. Then, we can rewrite the projection as:
In order to show in the following sections that the statistic $U_N$ and the projection $\hat{U}_{N,1}$ are asymptotically equivalent, we first need to derive the variance of the projection.
Proof. Proof provided in Appendix B.3.
Before defining the second possibility for the H\'{a}jek projection, notice first that, conditioning on individual and dyad-level random variables related to a directed dyad $(i',j')$ is different than conditioning on all individual and dyad-level random variables related to a combination $\{i',j'\}$. More specifically, the first comprises of the elements $A_{i'}$, $B_{j'}$ and $U_{i'j'}$, while the second comprises of $A_{i'}$, $A_{j'}$, $B_{i'}$, $B_{j'}$, $U_{i'j'}$ and $U_{j'i'}$.
Therefore, in this second proposed projection, instead of summing over all possible directed dyads, we sum over all possible combinations of indices $i'$ and $j'$, which amounts to $\frac{N(N-1)}{2}$ combinations. We therefore condition on all characteristics of these both indices:
Again, the double sum, $\sum_{i'=1}^N \sum_{j > i}$, fixes two indices $i'$ and $j'$ of the possible tetrads, and it runs over the conditioning terms. Then, we take the expectation of the statistic $U_N$ conditional on the terms described above. The structure of the second H\'{a}jek projection is essentially the same as the first, apart from which terms we condition on. Therefore, again, we will have that for all combinations of quadruples for which we take the average of the conditional expectation of the score function (kernel), only ${N-2 \choose 2}$ combinations will lead to non-zero expectated values. Those combinations refer again to the ones containing the elements $i'$ and $j'$.
The difference with respect to the previous projection, that is induced by the extra conditioning terms, boils down to the terms in the score function $s_{ijkl}$ (where $\{i',j'\} \subset \{i,j,k,l\}$, with, for instance, $i=i'$ and $j=j'$) that will be non-zero, since now the permutations that both the terms $U_{i'j'}$ and $U_{j'i'}$ will have non-zero terms.
Once again, we can further boil down the expression of the projection $\hat{U}_{N,2}$ by first noting that, as shown in Appendix C, for a given combination $\{i',j',k,l\}$ for any value of $k$ and $l$, the conditional expected value of the kernel evaluated at such combination is of the form:
Such that we can further simplify the expression for the projection:
Again, we proceed by deriving the variance of this second projection, which should be equivalent to the variance of the first proposed projection.
Proof. Provided in Appendix B.4.
The main idea of defining a H\'{a}jek projection is to obtain a statistic that is approximatelly and asymptotically close enough to the U-statistic to which central limit theorems and laws of large numbers can be applied.
serfling2009approximation provides readily applicable results for the asymptotic equivalence of the U-statistic given by Definition (ref) and the H\'{a}jek projection given by Definition (ref), and, consequently, to the asymptotic properties of the U-statistics, since in this case the projection is an average of i.i.d. random variables. However, in our case, as the statistic $U_N$ is not formally an U-statistic, such results cannot be immediately used.
To derive the asymptotic equivalence between $U_N$ and the two proposed projections, $\hat{U}_{N,1}$ and $\hat{U}_{N,2}$, we follow closely the arguments in graham2017econometric.
Following up this Remark, we have the following result:
Proof. Provided in Appendix B.5.
Hence, even though the statistic $U_N$ is not properly defined as an U-statistic, we still have the result that, under the assumptions needed for the results above, its limit distribution coincides with that of the proposed H\'{a}jek projections. This property is key to define the asymptotic properties of the pairwise differences estimator in the next section.
Considering the rewritten estimator, defined before:
We note that to obtain the asymptotic properties, it is key to obtain the convergence of the Hessian:
While in most of the studies on dyadic regressions and U-statistics tools associated with it the convergence of this Hessian is assumed (for instance, in graham2019network), we instead proof such convergence result. Observe that, even though this term also resembles a U-statistic, or, at least a term to which U-statitics tools can be applied to, it is not the case. This follows from not having a term such as $U_{ij}$ which is i.i.d. in the dyad level in such statistic, which would guarantee the conditional independence of summands. Therefore, the same tools applied to the statistic $U_N$ cannot be carried over here. Instead, our approach relies on deriving the variance of such term, and proving the convergence in probability through the Chebyshev's inequality.
Proof. Provided in Appendix B.6.
Given the result of the Proposition above, we can rewrite, by rescaling the expression of the rewritten estimator by $\sqrt{N(N-1)}$:
which follows by the continuous mapping theorem. Therefore, the asymptotic sampling properties of $\sqrt{N(N-1)} (\hat{\beta}_{1,PD} - \beta_1)$ will be driven by the behaviour of $\sqrt{N(N-1)} U_N$.
From Theorem (ref), we have that the statistic $U_N$ is asymptotically equivalent to the projections $\hat{U}_{N,1}$ and $\hat{U}_{N,2}$. Therefore, the asymptotic properties of those carry over to the asymptotic properties of $U_N$. Notice that, from Equation (ref) and its analogous for the second proposed projection, we have that the summands of the projections are uncorrelated, but not necessarily independently distributed. The dependence structure remains since the same individual characteristics $A_i$ and $B_j$ for a given $i$ and $j$ are still present in different summands, as for instance, we can have terms such as $\mathbb{E} [X_{ij'} \rvert B_{j'}]$ in one summand and $\mathbb{E} [X_{ik'} \rvert B_{k'}]$ in another summand.
However, as pointed out by graham2017econometric and jochmans2018semiparametric in their contexts, by law of iterated expectations, we can rewrite:
Such that the summands of the projections are conditionally independent when conditional on all indices $\{A_i\}_{i=1}^N$ and $\{B_j\}_{j=1}^N$, which is a characteristic of CID models, and that carries over to this context. Given this conditional independence of the random variables, we can assert:
Proof. Available in the next versions of this paper.
Following Proposition (ref) and Lemma (ref), we can establish first the consistency of the estimator, which is provided in the following theorem.
Proof. Provided in Appendix B.7.
From the same Proposition and Lemma, the asymptotic normality and the associated asymptotic variance of the estimator can be established. Also note that the estimator is asymptotically unbiased according to the following theorem.
From Theorem (ref), it trivially follows that: $$ \hat{\beta}_{1,PD} \overset{a}{\sim} N \left( \beta_1, \frac{1}{N(N-1)} \Gamma^{-1} 144 \delta_2 \Gamma^{-1}\right) $$ $$ \hat{\beta}_{1,PD} \overset{a}{\sim} N \left( \beta_1, \frac{1}{N(N-1)} \Gamma^{-1} 72 \Delta_2 \Gamma^{-1}\right) $$ Therefore, the asymptotic variance of the estimator $\hat{\beta}_{1,PD}$ can be estimated as: $$ \widehat{\text{AVar}}(\hat{\beta}_{1,PD}) = \frac{1}{N(N-1)} \hat{\Gamma}^{-1} 144 \hat{\delta}_2 \hat{\Gamma}^{-1} $$ $$ \widehat{\text{AVar}}(\hat{\beta}_{1,PD}) = \frac{1}{N(N-1)} \hat{\Gamma}^{-1} 72 \hat{\Delta}_2 \hat{\Gamma}^{-1}, $$ where: $$ \hat{\Gamma} = \frac{1}{{N \choose 4}} \sum_{\mathcal{C} \in \mathcal{C}(\mathcal{N},4)} \frac{1}{4!} \sum_{\pi \in \mathcal{P}(\mathcal{C},4)} \Tilde{X}_{\pi_1\pi_2\pi_3\pi_4} \Tilde{X}_{\pi_1\pi_2\pi_3\pi_4}',$$ where we then have that the asymptotic variance can be estimated using either a consistent estimator for $\delta_2$ or a consistent estimator for $\Delta_2$. In the following subsections we will propose consistent estimators for both.
As mentioned before, the definition of $\delta_2$ is:
where: $$\bar{s}_{ij} = \mathbbm{E}[s_{ijkl} \rvert A_{i}, B_{j}, U_{ij}]$$ Importantly, the elements which we condition on, namely, $A_{i}, B_{j}, U_{ij}$ have indices $i$ and $j$ that necessarily are in the combinations $\{i,j,k,l\}$ for any other elements $k$ and $l$. This reflects the fact that the term $\delta_2$ originates from the expression of the variance of the statistic $U_N$, considering the components in such variance that has two elements in common. To obtain a consistent estimator $\hat{\delta}_2$, we also need a consistent estimator $\hat{\bar{s}}_{ij}$.
graham2017econometric suggests that, for an undirected network the consistent estimators are: $$ \hat{\Delta}_{2,G} = \frac{1}{n} \sum_{i<j} \hat{\bar{s}}_{ij} \hat{\bar{s}}_{ij}', $$ $$ \hat{\bar{s}}_{ij,G} = \frac{1}{n - 2(N-1) + 1} \sum_{k<l, \{i,j\} \cap \{k,l\} = \emptyset} s_{ijkl}, $$ where $n = \frac{N(N-1)}{2}$ is the number of undirected dyads, therefore the expression for $\hat{\Delta}_{2,G}$ considers the average over all undirected dyads. The sum $\sum_{k<l, \{i,j\} \cap \{k,l\}}$ explicity means that, given two fixed indices $i$ and $j$ for the first dyad, we take the sum over all possible remaining different indiced $k$ and $l$, such that $k<l$, since in the context of graham2017econometric we have an undirected network, and therefore only the different combinations $\{k,l\}$ matters, but not its different permutations. Moreover, notice that $n - 2(N-1) + 1$ coincides with the ${N-2 \choose 2}$ tetrads that contain a fixed $i$ and $j$. Therefore, the expression for $\hat{\bar{s}}_{ij,G}$ averages over all the kernels of the combinations that contain $i$ and $j$.
As in our case we are looking at a directed network, some adjustments seem to be necessary. Especially, notice that, for a directed network we have that not necessarily $\bar{s}_{ij} = \bar{s}_{ji}$, since: $$\mathbbm{E}[s_{ijkl} \rvert A_i, B_j, U_{ij}] \neq \mathbbm{E}[s_{ijkl} \rvert A_j, B_i, U_{ji}]$$ That means that in the expression for the consistent estimator of $\delta_2$ we should average over all possible directed dyads:
One possibility is to work out further the expression for $\bar{s}_{ji}$, such that it does not simply boil down to, when estimated, the average over the kernels.
To be more precise, we can see that, when taking the expectation over the kernel $s_{ijkl}$ conditioning on the characteristics of a single dyad $\{i,j\}$, only some of its permutations (that are inside the kernel, and namely the ones that contain the idiosyncratic error term $U_{ij}$) will have a conditional expectation different than zero:
which is different than:
Therefore, if we take the sample analogue of those expressions applied to a given combination $\{i,j,k,l\}$ that contain the fixed elements $i, j$ and any elements $k, l$:
In the expressions above we plugged in the estimates of the idiosyncratic error terms, obtained from the estimated coefficient $\hat{\beta}_{1,PD}$, such that: $$ \hat{\Tilde{U}}_{ijkl} = \Tilde{Y}_{ijkl} - \hat{\beta}_{1,PD} \Tilde{X}_{ijkl}, $$ for any indices $i,j,k,l$.
Then, for these proposed consistent estimators we would have that $\hat{\bar{s}}_{ij} \neq \hat{\bar{s}}_{ji}$.
In this case, we have that the previous definition of $\Delta_2$ is:
As the H\'{a}jek projection in this case was obtained by summing all combinations (and not permutations) of indices $i$ and $j$, we have that the consistent estimator of $\Delta_2$ should average over all these possible combinations:
Moreover, remembering that $\bar{s}_{ij,2}$ is the kernel conditioning on all characteristics of $i$ and $j$, we have that its estimator, $\hat{\bar{s}}_{ij,2}$ is given by:
where, again, in the expression above we plugged in the estimates of the idiosyncratic error terms, obtained from the estimated coefficient $\hat{\beta}_{1,PD}$, such that.
With both consistent estimates of the covariances, it is then possible to conduct valid inference. Moreover, in the next Section we investigate the finite sample performance of both analytical estimates.
In this section, we explore the finite sample properties of the estimator $\hat{\beta}_{1,PD}$ through a Monte Carlo simulation exercise. We also aim to evaluate the finite sample properties of the estimator of the asymptotic variance of $\hat{\beta}_{1,PD}$, and the associated t-tests using both the consistent estimator $\hat{\delta}_2$, based on the first obtained Hájek projection, and the consistent estimator $\hat{\Delta}_2$, based on the second. In a nutshell, we find that: (i) the estimated slope parameters are unbiased in general, even when the fixed effects are correlated with the covariates; (ii) the estimated asymptotic variances using either estimators are very close to each other, which was expected; and (iii) the size of the t-tests are correct, indicating a valid inference procedure.
For simplifying purposes, for now, we consider the case of a single regressor $X_{ij}$ in the different proposed designs. In general, I follow closely the DGP specifications proposed by jochmans2018semiparametric and charbonneau2017multiple, who also consider a directed network model. Note, however, that in their cases, they consider a binary outcome variable, while we consider a continuous dependent variable.
The DGP in general follows:
$$ Y_{ij} = \beta_1 X_{ij} + \theta_i + \xi_j + U_{ij} $$
In all different designs, we take $\beta_1 =0$. The idiosyncratic error terms $U_{ij}$ are independently drawn from a standard normal distribution, $U_{ij} \sim N(0,1)$. In our case of a directed network, we specifically have that $U_{ij} \neq U_{ji}$, therefore for a simulation considering $N$ nodes we draw from the standard normal $N(N-1)$ idiosyncratic errors. The fixed effects $\theta_i$ and $\xi_j$ are also drawn from standard normal distributions.
The difference among the designs relies on how the regressor $X_{ij}$ is drawn.
Here we follow essentially the same DGP as proposed by jochmans2018semiparametric. We generate the single regressor as:
$$X_{ij} = - \rvert A_i - B_j \rvert $$
where $ A_i = V_i - \frac{1}{2}$, for $V_i \sim \text{Beta}(2,2)$, and the same for $B_j$. The covariate thus is generated in such a way that is dependent across both senders and receivers in the dyadic relation. The difference to the DGP proposed by jochmans2018semiparametric relies on the fact that we consider the individual effect of the alter, $A_i$, to be different and drawn independently from that of the ego, $B_j$, while jochmans2018semiparametric considers $B_j = A_j$.
We introduce a correlation between the regressor $X_{ij}$ and the fixed effects $\theta_i$ and $\xi_j$, such that:
$$X_{ij} = - \rvert A_i - B_j \rvert + \theta_i + \xi_j $$
where $ A_i = V_i - \frac{1}{2}$, for $V_i \sim \text{Beta}(2,2)$. Also, $B_j$ is drawn independently from $A_i$, such that $ B_j = V_j - \frac{1}{2}$, for $V_j \sim \text{Beta}(2,2)$. Note that the manner in which we introduce a correlation between the regressor and the fixed effects is similar to that of charbonneau2017multiple.
We now consider a binary regressor that is uncorrelated with the fixed effects. We generate the regressor according to:
$$X_{ij} = \mathbbm{1} \{ A_i - B_j > 0 \}$$
where $A_i$ and $B_j$ are drawn according to Designs 1 and 2.
We again consider a binary regressor, however, now it is correlated with the fixed effects, such that:
$$X_{ij} = \mathbbm{1} \{ A_i - B_j + \theta_i + \xi_j > 0 \}$$
where, again, $A_i$ and $B_j$ are drawn according to Designs 1 and 2.
We propose several settings of Monte Carlo simulations. For each design, we run simulations for $S \in \{1000, 5000, 10000\}$, where $S$ refers to the number of simulations, and for $N \in \{10, 20, 30, 50\}$.
In the tables below we show the results for the estimator $\hat{\beta}_{1,PD}$ in terms of biasedness, as well as its variance across the simulations. We also present the results for the average of the estimated asymptotic variance considering both the estimations taking into account $\hat{\delta}_2$, according to Equation (ref), and taking into account $\hat{\Delta}_2$, according to Equation (ref).
From the tables above, we point out two results: (i) the estimator $\hat{\beta}_1$ seems to be unbiased, and (ii) the mean of the estimated variaces is basically on spot when compared to the variance of $\hat{\beta}_1$ across the simulations and across the different designs. More specifically, while for all designs (except for Design 3) there is still some bias in the simulations for $N=10$ and $S=1000$, the bias essentially vanishes as we consider larger numbers of nodes $N$, or a larger number of simulations $S$.
Another feature that was already expected is that the average of the estimated asymptotic variances are very close when comparing to whether the variance was estimated using $\hat{\delta}_2$ or $\hat{\Delta}_2$. Moreover, we notice that in general those averages are almost spot on with the variances of the estimated $\hat{\beta}_1$ across simulations. The only exception are the simulations with $N=10$ for designs 1 and 2, however, as soon as $N$ is increased the results are again essentially the same. This indicates that the variance estimator captures well the small-sample variability in the point estimator, and that inference using such estimators is valid.
We next explore if normality might be a good approximation to the finite sample distribution of the proposed estimator $\hat{\beta}_1$. We present below the histograms and the QQ-plots of the estimated values for Designs 2 and 4, which are considered to be the most relevant, since it allows for correlations between the covariates and the fixed-effects. However, the plots for the other designs can be found in Appendix D.
From the histograms above there is evidence that the estimator of $\beta_1$ is normally distributed as the size of $N$ increases for all the number of simulations $S$. More specifically, for smaller values of $N$, we can see that the range of the histogram is wider than the one from a normal distribution. This is corroborated by the QQ-plots, that show that for any number of simulations $S$, for lower values of $N$, the distribution seems to have fatter tails than a normal distribution, but as $N$ increases it seems to be distributed as a normal.
We next examine the size of the t-tests where the test statistic use the asymptotic variance estimator proposed before. We test the null hypothesis that the coefficient $\beta_1$ is equal to its true value, $\beta_1 = 0$. The tables below shows the fractions of samples for which the null hypothesis is rejected at the $5\%$ statistical significance level.\\ \\
When we look at the size of the t-test for the different variance estimators, we see that, as expected from the previous findings, the sizes for the estimators are close to 0.05, however the estimates using $\hat{\Delta}_2$ are somewhat closer than those using $\hat{\delta}_2$.
In this paper we showed how one can adapt U-statistics tools to show the asymptotic properties of linear dyadic models for network data. More especifically, we proposed a linear model with two-way fixed effects that enter additively in the specification. While the usual two-way fixed effects estimator is consistent and asymptotically unbiased for this particular model, we propose an estimator that relies on pairwise differences, that completely eliminate the fixed effects from the objective (and influence) function(s). This choice of estimator was done with the purpose of demonstrating step-by-step in a simple model how one can adapt tools from U-statistics to this particular dyadic setting (with a pairwise differences estimator) to obtain an analytical form and an estimator for the asymptotic variance.
These specific tools are needed because the pairwise differencing approach introduces a dependence structure in the summands of the influence function of the estimator. A similar set of tools are also used in non-linear models that employ a similar estimation method, in particular in charbonneau2017multiple and in graham2017econometric. For non-linear models, differencing out the fixed effects is desirable to eliminate the incidental parameter problem, which would lead to asymptotically biased estimates of the coefficients of the covariates.
With a Monte Carlo exercise, we showed that the obtained estimates for the slope coefficients are unbiased in finite samples, and the estimated asymptotic variance delivers the correct size for the t-test. However, the model assumed in this paper can still be relaxed to allow for a richer dependence structure in the network in future research. For instance, we also could allow for dependencies across outcomes $Y_{ij}$ and $Y_{ji}$ by relaxing Assumption (ref) such that the idiosyncratic terms $U_{ij}$ and $U_{ji}$ are allowed to covary. In practice, this would have implications for how the tools of U-statistics are employed in this dyadic framework. Namely, the Hoeffding decomposition for the variance of the U-statistic and the H\'{a}jek projection would have to be modified to allow for these dependencies.
Another avenue of interest is that we could allow for, in Assumption (ref), the individual-level observed characteristics $A_i$ and $B_i$ to covary. That would allow, for instance, that exports from Japan to Korea might covary with those from Korea to Thailand. In our derivations, this would have implications for the probability limit of the Hessian of the proposed estimator.
Finally, the main computational challenge is the estimate of $\Delta_2$. It could also be of pratical use, in future research, to explore possible bootstrap procedures to obtain inference, such as in graham2019network and in menzel2018bootstrap.