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.
80,065 characters · 19 sections · 33 citation commands
Specification Testing for Dyadic Regression Models
\noindentKeywords: dyadic data, specification testing, residual-marked empirical process, exchangeable arrays, multiplier bootstrap.\\ JEL Classification: C12, C14, C21, C31.
Dyadic data are increasingly common in empirical economics. International trade is observed between pairs of countries, credit and ownership relationships connect pairs of firms or financial institutions, and social and professional networks record interactions between pairs of individuals. A central feature of such data is that two observations may be dependent whenever their dyads share a node. Bilateral trade flows involving the same country, for example, may respond to common country-level shocks even when the trading partners differ. This shared-node dependence is the defining inferential challenge in dyadic regression; see, e.g., fafchamps2007formation, aronow2015cluster, and tabord2019inference.
Most existing methods for dyadic regression conduct inference under a maintained specification of the conditional mean. In practice, however, the functional form of that conditional mean is rarely known. Linear specifications remain attractive because they are transparent, easy to estimate, and convenient for counterfactual analysis. Their conclusions can nevertheless be misleading when economically relevant nonlinearities or interactions are omitted; see, for example, white1981consequences and, in a related dependent-data setting, su2010profile. This concern is particularly important in network applications, where the effect of one characteristic may depend on the characteristics of both nodes. A researcher therefore needs a way to assess the maintained regression specification before interpreting its coefficients.
This paper develops omnibus specification tests for linear conditional-mean models with undirected dyadic data. Under the null hypothesis, the conditional mean of the dyadic outcome is linear in the included regressors, whereas the alternative leaves it unrestricted. We characterize the null through residual moments indexed by lower orthants of the regressor distribution and construct Kolmogorov-Smirnov (KS) and Cram\'er-von Mises (CvM) statistics from the resulting residual-marked empirical process.
Extending residual-based specification tests to dyadic data is not a matter of replacing an independent-observation variance estimator with a dyadic-robust one. A dyadic empirical process contains two sources of sampling variation. The first is generated by latent node characteristics and is shared by all dyads incident to the same node. The second is specific to individual dyads. When the first component is present, the effective sample size is governed by the number of nodes. When it vanishes, the dyad-specific component may become leading, and the effective sample size can instead be governed by the number of dyads. The normalization, limiting covariance, and appropriate bootstrap can therefore change with the variance regime.
This distinction is difficult to handle for an omnibus test. For a fixed residual mark, a standard projection decomposition identifies the node-level component and a degenerate second-order remainder. An omnibus statistic, however, depends on a continuum of marks. Pointwise negligibility of the remainder does not establish that its supremum or integrated square is negligible. Furthermore, a bootstrap that correctly reproduces the node component need not reproduce the dyad-specific component when the former disappears.
First, we establish a uniform projection theorem for dyadic empirical processes indexed by a VC-type class. The theorem shows that the second-order component is negligible uniformly over the index set at the node-level rate. The proof is nontrivial because the remainder contains both dyad-specific variation and variation generated jointly by the two node characteristics. After separating these components, we control the former as a Rademacher process and the latter as a decoupled second-order Rademacher chaos. The result reduces the leading dyadic process uniformly to an ordinary empirical process over the latent nodes and provides the functional weak convergence required for both the KS and CvM statistics.
Second, we develop a raw first-order node-multiplier bootstrap and prove its validity under nondegenerate shared-node dependence. The relevant node projections are functions of latent variables and must be recovered from averages over the dyads incident to each node. We establish this recovery uniformly over the class of residual marks. We then control the additional effects of estimated residuals and estimated orthogonalization terms by embedding the random estimated class in a deterministic enlargement with controlled entropy.
Third, we identify a failure of the raw node-multiplier bootstrap under first-order degeneracy and propose a covariance correction. Because every dyad is incident to two nodes, it enters two node averages. The raw node bootstrap consequently counts the same-dyad variation twice. This duplication is negligible when node-level variation dominates, but it becomes first-order when the node projection vanishes. In the independent-dyad case, the raw bootstrap therefore has twice the correct asymptotic covariance.
An exact finite-sample covariance decomposition separates the same-dyad contribution from the covariance between distinct dyads sharing a node. This decomposition leads to a corrected Gaussian bootstrap that retains the shared-node component while removing the extra same-dyad contribution. We prove that the correction is asymptotically irrelevant under nondegenerate shared-node dependence and essential under independent dyads. The corrected procedure thus agrees with the raw node bootstrap in the usual dyadic regime while also recovering the faster independent-dyad limit when node-level variation disappears.
The theory permits partial degeneracy: the node-level covariance may vanish at some evaluation points without invalidating the functional approximation. It also covers the totally first-order degenerate benchmark of mutually independent dyads. We do not claim validity for every totally degenerate exchangeable array. In more general cases, the second-order component may become leading and may have a non-Gaussian limit. The paper makes this boundary explicit rather than imposing a Gaussian approximation where it is not justified.
Building on these results, we establish the asymptotic validity and power of the proposed tests. They are consistent against fixed violations of the conditional-mean restriction and have nontrivial power against local alternatives at the rate appropriate to each variance regime. The local-power analysis also accounts for re-estimation of the linear coefficient and identifies the component of a local departure that remains detectable by the specification process. Despite the nonstandard theory, implementation on a fixed evaluation grid requires only one OLS estimation and no bandwidth selection or repeated estimation across bootstrap draws.
The Monte Carlo experiments examine the finite-sample size and local power of the proposed tests across different sample sizes and strengths of shared-node dependence. We compare the raw node-multiplier and the corrected procedures with tests that incorrectly treat dyads as independent. The results illustrate the practical importance of accounting for dyadic dependence and show how the relative performance of the procedures changes as the node-level component becomes weaker. In particular, the simulations confirm the theoretical distinction between nondegenerate shared-node dependence and the independent-dyad benchmark.
We also apply the proposed tests to a law-firm network. The empirical analysis examines whether professional links among lawyers can be adequately represented by an additive linear conditional mean based on individual and pairwise characteristics. The results indicate that professional connections cannot be adequately described by adding the separate contributions of experience, demographic similarity, and organizational proximity. Instead, these characteristics appear to operate jointly. The application also demonstrates how the proposed procedures can be used as formal diagnostic tools before interpreting a linear dyadic regression and illustrates their ability to detect economically meaningful departures from the maintained specification.
Our paper is related to several strands of literature. General empirical-process and bootstrap theory for exchangeable arrays is developed by davezies2021, but our results are not a direct application of that theory. We establish an explicit uniform projection theorem that separately controls the dyad-specific empirical-process component and the node-pair Rademacher chaos, uniformly recover the unobserved node projections from incident-dyad averages, and derive conditional maximal and contraction inequalities needed to handle estimated residual marks. Other related contributions to inference for exchangeable arrays and dyadic data include the high-dimensional Gaussian approximations of chiang2023inference, the analysis of degeneracy under two-way clustering by menzel2021bootstrap, and nonparametric methods for dyadic regression and density estimation developed by graham2021minimax and graham2024kernel, respectively.
The paper also contributes to the literature on omnibus conditional-moment specification testing, including bierens1982consistent, bierens1990, stute1997, whang2000consistent, and escanciano2006. Kernel-based approaches include zheng1996 and fan1996consistent, while specification testing under spatial dependence is studied by su2017specification, gupta2024consistent, yang2024model, and lee2025heteroskedasticity. Unlike these settings, dyadic sampling generates competing node-level and dyad-level variance components, which determine both the functional limit and the validity of the bootstrap.
The remainder of the paper is organized as follows. Section (ref) introduces the sampling framework and the residual-marked process. Section (ref) develops the uniform projection theory and derives the functional limits. Section (ref) presents the raw and covariance-corrected procedures, establishes their validity, and studies local power. Section (ref) shows the finite-sample simulation results. Section (ref) further demonstrates the practical relevance of the proposed tests through an empirical application. The Appendix contains all proofs and auxiliary results.
The following notation is used throughout. For a probability measure \(Q\) and a measurable function \(f\), write \( Qf=\int f\,dQ\) and \(\|f\|_{Q,p}=(Q|f|^p)^{1/p}. \) For vectors, \(\|\cdot\|\) denotes the Euclidean norm; for matrices, it denotes the induced operator norm. A prime denotes transposition, and \(\lambda_{\min}(A)\) denotes the smallest eigenvalue of a symmetric matrix \(A\). We write \(\mathbf{1}\{\cdot\}\) for the indicator function and write \(a\preceq b\) for coordinatewise inequality, i.e., every coordinate of \(a\) is no larger than the corresponding coordinate of \(b\).
For an index set \(T\), \(\ell^\infty(T)\) denotes the Banach space of bounded real-valued functions on \(T\), equipped with the supremum norm \(\|z\|_\infty=\sup_{t\in T}|z(t)|\). Thus a process \(Z_n=\{Z_n(t):t\in T\}\) is viewed as an \(\ell^\infty(T)\)-valued random element, and \( Z_n\rightsquigarrow Z\) in \(\ell^\infty(T) \) denotes weak convergence of the entire process. We use \(\to_p\), \(o_p(1)\), and \(O_p(1)\) for convergence in probability, convergence to zero in probability, and boundedness in probability, respectively; \(a_n\asymp b_n\) means that \(a_n/b_n\) is bounded above and away from zero. Conditional on the observed dyadic sample, expectation and probability are denoted by \(E^*\) and \(P^*\); \(o_{p^*}(1)\) denotes convergence to zero in conditional probability, in outer probability when measurability requires it. Finally, \(BL_1\) is the set of real-valued functions bounded by one and Lipschitz with constant at most one on the relevant metric space.
Let \(i,j\in\{1,\ldots,n\}\) index nodes. We observe one undirected dyadic observation \(Z_{ij}=(Y_{ij},X_{ij}')'\) for each unordered pair \((i,j)\) with \(1\leq i<j\leq n\), where \(X_{ij}\in\mathbb{R}^d\) includes an intercept. Let \(\mathcal D_n=\{(i,j):1\leq i<j\leq n\}\) and \(N_n=|\mathcal D_n|=n(n-1)/2\). For notational convenience, whenever \(j<i\), set \(Z_{ij}=Z_{ji}\), and use the same convention for all dyad-level functions, marks, and latent dyad variables. For a measurable function \(f\), define \(\mathbb{P}_{n} f=N_n^{-1}\sum_{i<j}f(Z_{ij})\) and \(Pf=E[f(Z_{12})]\). For a function class \(\mathcal{F}\), the node-scaled centered dyadic empirical process is the random map
In the independent-dyad scenario introduced below, the nondegenerate normalization is instead \(\sqrt{N_n}\).
The dyadic dependence assumption allows two dyadic scores $s_{ij}$ and $s_{pq}$ to be dependent only when the two dyads share at least one endpoint, that is, when \[ \{i,j\}\cap\{p,q\}\neq \varnothing . \] The array is modeled as jointly exchangeable and dissociated. Joint exchangeability means that relabeling the nodes does not change the joint distribution. Dissociation means that collections of dyads involving disjoint sets of nodes are independent. To accommodate this feature, we describe the dependence structure using the Aldous-Hoover-Kallenberg (AHK, aldous1981; hoover1979; kallenberg1989representation) representation.
Under Assumption (ref), the dyadic observations $Z_{ij}$ are identically distributed but need not be independent. Define \(Q=E[X_{12}X_{12}']\) and suppose that \(Q\) is nonsingular. Whether or not the null is true, define the population linear projection coefficient \(\beta_*=Q^{-1}E[X_{12}Y_{12}]\) and the projection residual \(\varepsilon_{ij}=Y_{ij}-X_{ij}'\beta_*\). By construction, \(E[X_{12}\varepsilon_{12}]=0\).
We test
against the unrestricted alternative that \(E[Y_{12}\mid X_{12}]\) is not linear almost surely. Equivalently, \(H_0\) states that \(E[\varepsilon_{12}\mid X_{12}]=0\) almost surely.
Write \(X_{ij}=(1,W_{ij}')'\). Define \(M(x)=E[X_{12}\mathbf{1}\{W_{12}\preceq x\}]\) and
The population residual-marked moment is $ \Delta(x)=E[\varepsilon_{12}\mathbf{1}\{W_{12}\preceq x\}].$ The adjustment in (ref) yields a Neyman-orthogonal moment. Since \(E[X_{12}\varepsilon_{12}]=0\), the following equality holds regardless of whether the null hypothesis is satisfied:
Lemma (ref) shows that lower-orthant indicators generate an omnibus collection of instruments. No smoothing parameter is needed. This is an important difference from kernel-based quadratic-form tests such as zheng1996.
Define the OLS estimator \(\widehat\beta=Q_n^{-1}\mathbb{P}_{n}(XY)\), where \(Q_n=\mathbb{P}_{n}(XX')\). Let \(\widehat\varepsilon_{ij}=Y_{ij}-X_{ij}'\widehat\beta\) and \(M_n(x)=\mathbb{P}_{n}[X\mathbf{1}\{W\preceq x\}]\). The feasible residual-marked process is
Define the sample orthogonalized instrument
The OLS normal equations imply \(\mathbb{P}_{n}(X\widehat\varepsilon)=0\), and therefore the following equality is exact whenever \(Q_n\) is nonsingular:
Therefore, $\widehat{R}_n(x)$ is the sample analog of $\Delta(x)$ as given in (ref).
Bounded regressors are imposed to keep the empirical-process verification transparent. They can be replaced by suitable moment and weighted-entropy conditions. The moment condition is used to estimate the two covariance components in the corrected bootstrap.
The asymptotic argument has two layers. The first is probabilistic: a dyadic empirical process can be reduced uniformly to its first-order node projection. The second is statistical: the process based on estimated OLS residuals is related uniformly to a fixed population-indexed class. We treat these layers in turn.
For a square-integrable function \(f\) of a dyadic observation, let \(\mu_f=Pf=E[f(Z_{ij})]\) and define its first- and second-order projection terms by
The subscripts “1” and “2” indicate the orders of the corresponding projection terms. By symmetry of the dyadic sampling structure, the two first-order projections coincide: \( f_1(u)=E[f(Z_{ij})\mid U_j=u]-\mu_f. \) Moreover, by construction, \[ E[f_1(U_i)]=0,\qquad E[f_2(U_i,U_j,U_{ij})\mid U_i]=0,\qquad E[f_2(U_i,U_j,U_{ij})\mid U_j]=0. \]
Lemma (ref) in the Appendix establishes the exact decomposition
where \( \mathbb U_n f_2 = N_n^{-1}\sum_{i<j}f_2(U_i,U_j,U_{ij}). \) For each fixed \(f\), the degeneracy of \(f_2\) implies that the second term is asymptotically negligible. The following result strengthens this pointwise conclusion by establishing uniform negligibility over the function class.
The proof of Lemma (ref) is nontrivial. It first decomposes \(f_2\) into a component driven by the dyad-specific latent variable and a component driven by the two node variables. After symmetrization, the former yields an ordinary Rademacher process, whereas the latter yields a decoupled Rademacher chaos of order two and therefore requires a separate chaining argument based on hypercontractivity. VC-type entropy bounds provide uniform control of both components.
Combining this reduction with the i.i.d. empirical-process limit for \(\{f_1:f\in\mathcal{F}\}\) yields the following general result for the dyadic empirical process.
For the inference analysis in Section (ref), we introduce an ideal node-multiplier process indexed by the fixed population class \(\mathcal{F}\). This process is infeasible in the specification-testing application because the relevant marks \(r_x\) depend on unknown population quantities, but it provides the benchmark for the feasible bootstrap developed later.
Let \((\xi_i)_{i=1}^n\) be i.i.d. multipliers, independent of the dyadic sample, such that
where \(C_\xi<\infty\). Define the ideal node-multiplier process by
The centered incident-dyad average in braces estimates the first-order node projection \(f_1(U_i)\). The following lemma shows that the multiplier process therefore reproduces the functional limit in Lemma (ref).
Lemma (ref) establishes validity for the ideal process. Section (ref) shows that replacing the population marks \(r_x\) by their feasible estimated counterparts is asymptotically negligible. We next turn from this general process theory to the feasible residual-marked sample process, first deriving its functional limits across the relevant variance regimes before constructing its feasible bootstrap counterpart.
For \(x\in\mathcal{X}\), define $r_x(Z)=(Y-X'\beta_*) q_x(X)=\varepsilon q_x(X)$ and $\mathcal{R}=\{r_x:x\in\mathcal{X}\}$. We now specialize the general projection theory to the marked class \(\mathcal{R}\). Appendix (ref) verifies that \(\mathcal{R}\) is VC type with a square-integrable envelope and establishes the preliminary uniform laws for \(Q_n\), \(M_n\), and \(\widehat\beta\). The following lemma connects the feasible residual-marked process to the population-indexed dyadic empirical process.
Lemma (ref) reduces the asymptotic behavior of the feasible process to that of \(\{\mathbb{G}_{n} r_x:x\in\mathcal{X}\}\). By Lemma (ref), the behavior of this process at the \(\sqrt n\) rate is determined by the first-order node projections $\frac{2}{ n}\sum_{i=1}^n\psi_x(U_i)$, where
Define the node-level covariance kernel
and the dyad-level covariance kernel
The relative importance of these two covariance components determines the appropriate normalization and, later, the appropriate bootstrap procedure.
Scenario (i) allows \(\Omega(x,x)=0\) at some indices and therefore permits partial degeneracy, but it requires the node projection to be nontrivial somewhere on \(\mathcal{X}\). Scenario (ii) instead imposes total first-order degeneracy while retaining nontrivial dyad-level variation. These two scenarios are the leading cases considered in the literature, although typically not for the more general supremum-type process studied here; see, for example, mackinnon2021wild, chiang2023standard, and hounyo2024wild. Two regimes lead to different effective sample sizes and are treated separately below.
Lemma (ref) reduces the asymptotic behavior of the feasible process to that of the population-indexed dyadic empirical process \(\{\mathbb{G}_{n} r_x:x\in\mathcal{X}\}\). By Lemma (ref), its limiting behavior at the \(\sqrt n\) rate is determined by the first-order node projection of \(r_x\). Combining Lemmas (ref) and (ref) gives the node-scale functional limit.
Theorem (ref) does not require \(\Omega(x,x)\) to be strictly positive at every \(x\). If \(\Omega(x,x)=0\) at an isolated point or on a subset of \(\mathcal{X}\), then \(\mathbb G_R(x)=0\) at those indices, but the functional convergence remains valid. If \(\Omega(x,x)=0\) for every \(x\in\mathcal{X}\), the theorem remains correct but yields only the degenerate limit $\sqrt n\{\widehat R_n-\Delta\} \rightsquigarrow0$ in $\ell^\infty(\mathcal{X}).$ Scenario (ii) of Assumption (ref) identifies an important case in which a nondegenerate limit can instead be obtained at the faster \(\sqrt{N_n}\) rate.
Theorem (ref) complements Theorem (ref) by describing a variance regime in which the first-order node projection vanishes. When all dyads are independent, \( E[r_x(Z_{12})\mid U_1]=\Delta(x) \) almost surely, so that \(\psi_x(U_1)=0\) and consequently \(\Omega(x_1,x_2)=0\) for every \(x_1,x_2\in\mathcal{X}\). Theorem (ref) therefore remains valid, but yields only the degenerate conclusion \[ \sqrt n\{\widehat R_n-\Delta\}\rightsquigarrow0 \quad\text{in }\ell^\infty(\mathcal{X}). \] Theorem (ref) identifies the next nondegenerate order: the appropriate normalization is \(\sqrt{N_n}\) rather than \(\sqrt n\). The limiting covariance is then generated by the dyad-level marks themselves, as represented by \(\Gamma\) in (ref), rather than by their first-order node projections.
We conclude this section by recording the behavior of the feasible process under fixed alternatives. By Lemma (ref), failure of \(H_0\) implies that \(\Delta\) is not identically zero.
Theorem (ref) gives the full-index consistency result that underlies the tests: the continuum KS functional detects every fixed violation of \(H_0\), while the continuum CvM functional detects every violation satisfying \(\int\Delta(x)^2d\nu(x)>0\). On any fixed grid \(\mathcal{X}_G=\{x_1,\ldots,x_G\}\), with nonnegative weights \(w_g\) summing to one, the corresponding separation conditions are \(\max_{g\leq G}|\Delta(x_g)|>0\) for KS and \(\sum_{g=1}^G w_g\Delta(x_g)^2>0\) for CvM. Establishing consistency of the implemented grid tests additionally requires the grid to capture the departure and the bootstrap critical values to be valid. We construct those critical values next.
The preceding section establishes the functional limit of the feasible residual-marked process. We now turn that limit theory into an implementable testing procedure. We first construct the raw node-multiplier bootstrap and then a covariance-corrected Gaussian bootstrap. We apply both procedures to KS and CvM functionals and establish their respective validity regions, then characterize local asymptotic power. The full-index processes are useful for developing the functional theory, but both procedures are implemented by evaluating the sample and bootstrap processes on the same prespecified finite grid.
Theorem (ref) shows that the first-order law of the residual-marked process is generated by the latent node projections \(\psi_x(U_i)\) in (ref). Because \(\psi_x(U_i)\) is unobserved, these projections cannot be used directly. They can, however, be recovered from incident-dyad averages. Conditional on \(U_i\), averaging over all dyads containing node \(i\) integrates out the latent variation associated with the other endpoint and the dyad-specific shock. Centering the resulting averages across nodes then removes the unknown population moment.
To implement this idea, define the feasible residual mark
For each node \(i\), let
Here, the residual mark $\widehat{r}_{ij}(x)$ and the estimated first-order node projection $\widehat \psi_i(x)$ are the sample analogs of $r_x(Z_{ij})$ and $\psi_x(U_i)$, respectively. Since every dyad enters exactly two incident-dyad averages, \( \frac{1}{n}\sum_{i=1}^n\widehat{\bar r}_{x,i} = \mathbb{P}_{n}\widehat r(x) = \widehat R_n(x). \) Hence, \(n^{-1}\sum_{i=1}^n\widehat\psi_i(x)=0\) exactly.
Recall the multipliers specified in (ref). Replacing the infeasible population marks in (ref) with the estimated node projections \(\widehat\psi_i(x)\) gives the feasible raw node-multiplier process
The factor two reflects the two symmetric positions in which each node enters an undirected dyad. Although \(\widehat R_n^*(x)\) is defined for every \(x\in\mathcal{X}\) to support the functional limit theory, its implementation requires evaluation only at the selected grid points.
The procedure can be based on either of two complementary functionals. The KS statistic captures the largest localized departure from the null, whereas the CvM statistic aggregates departures over the index set. Their continuum versions are
and, for a fixed finite Borel measure \(\nu\) on \(\mathcal{X}\),
These continuum statistics provide a convenient full-index formulation of the asymptotic theory; they are not the statistics computed in the implementation. In practice, both bootstrap procedures evaluate the sample and bootstrap processes on the same prespecified finite grid \(\mathcal{X}_G\), replacing the supremum over \(\mathcal{X}\) by the maximum over \(\mathcal{X}_G\). For the CvM statistic, we use the discrete measure \(\nu_G=\sum_{g=1}^G w_g\delta_{x_g}\), where \(\delta_{x_g}\) denotes the probability measure placing unit mass at \(x_g\), \(w_g\geq0\), and \(\sum_{g=1}^G w_g=1\); equal weights \(w_g=1/G\) are the default choice.\footnote{There is no universally optimal \(G\); we use a sufficiently fine grid and assess robustness to further refinement. Alternatively, following stute1997, one may evaluate the process at the observed marks and use their empirical distribution for the CvM measure. Its randomness is asymptotically negligible under the maintained uniform-convergence conditions.} The implemented sample statistics are
These sample statistics are common to the raw and corrected procedures; only the bootstrap statistics and their critical values differ.
For the auxiliary full-index formulation of the raw procedure, define \[ T_n^{KS,\rm raw,*} = \sup_{x\in\mathcal{X}}|\sqrt{n}\widehat R_n^*(x)|,\qquad T_n^{CvM,\rm raw,*} = n\int_{\mathcal{X}}[\widehat R_n^*(x)]^2\,d\nu(x). \] The raw bootstrap statistics used in implementation are instead \[ T_{n,G}^{KS,\rm raw,*} = \max_{g\leq G}|\sqrt n\,\widehat R_n^*(x_g)|, \qquad T_{n,G}^{CvM,\rm raw,*} = n\sum_{g=1}^G w_g\{\widehat R_n^*(x_g)\}^2. \] For \(S\in\{KS,CvM\}\), let \(c_{n,1-\alpha}^{S,\rm raw,*}\) and \(c_{n,G,1-\alpha}^{S,\rm raw,*}\) denote the conditional \((1-\alpha)\)-quantiles of \(T_n^{S,\rm raw,*}\) and \(T_{n,G}^{S,\rm raw,*}\), respectively. The full-index comparison rejects \(H_0\) whenever \( T_n^S>c_{n,1-\alpha}^{S,\rm raw,*}, \) and is retained below as a theoretical consequence of functional weak convergence. The implemented raw test rejects whenever \( T_{n,G}^S>c_{n,G,1-\alpha}^{S,\rm raw,*}. \)
\paragraph{Algorithm 1: Raw node-multiplier bootstrap.} Fix a grid \(\mathcal{X}_G=\{x_g:g\leq G\}\subset\mathcal{X}\), a statistic \(S\in\{KS,CvM\}\), a nominal level \(\alpha\), and the number of bootstrap draws \(B\). Proceed as follows.
The ideal multiplier process in Lemma (ref) is indexed by the infeasible population marks \(r_x\). In contrast, \(\widehat R_n^*\) uses estimated residuals, \(Q_n\), and \(M_n(x)\). The next lemma shows that these plug-in operations are uniformly negligible. Its proof embeds the random estimated class in a deterministic VC-type enlargement and then combines a dyad-to-node contraction with a conditional maximal inequality.
Lemma (ref) transfers the conditional limit of the ideal multiplier process to its feasible counterpart. Together with Lemma (ref), it yields the following theorem.
When \(\mathbb G_R\) is nontrivial, this theorem gives a valid bootstrap for the \(\sqrt n\)-scaled sample process. Under independent dyads, \(\mathbb G_R\equiv0\), so convergence at the node scale does not justify critical values for the nondegenerate \(\sqrt{N_n}\)-scaled statistic.
The preceding node multiplier reproduces the node projection but counts each dyad once through each endpoint. This is asymptotically harmless when the node component has order one, but it doubles the leading same-dyad variance when all dyads are independent. We therefore introduce a correction that keeps the shared-node term and removes exactly one of the two copies of the same-dyad term.
The coefficient is obtained from an exact counting identity. For any centered vector mark \(\bm a_{ij}\), let \(V_0=E[\bm a_{12}\bm a_{12}']\) and \(V_1=E[\bm a_{12}\bm a_{13}']\). Dissociation implies
The first coefficient counts the \(N_n\) same-dyad terms; the second counts the \(n(n-1)(n-2)\) ordered pairs of distinct dyads sharing one node.
Let \(\mathcal{X}_G=\{x_1,\ldots,x_G\}\) be a fixed deterministic grid and write \(\widehat{\bm r}_{ij}=(\widehat r_{ij}(x_1),\ldots, \widehat r_{ij}(x_G))'\), \(\widehat{\bm R}_n= (\widehat R_n(x_1),\ldots,\widehat R_n(x_G))'\), and \(\widehat{\bm a}_{ij}=\widehat{\bm r}_{ij}-\widehat{\bm R}_n\). Define the same-dyad and shared-node covariance estimators
The covariance matrices of the raw node-multiplier vector and its finite-sample-corrected counterpart on \(\mathcal{X}_G\) are, respectively,
The corrected matrix (ref) mirrors the finite-sample covariance decomposition in (ref). Indeed, expanding the incident sums in (ref) gives \[ \frac4n\sum_{i=1}^n\widehat{\bm\psi}_i\widehat{\bm\psi}_i' =\widehat K_n^{\rm raw}, \qquad \widehat K_n^{\rm raw}-\widehat K_n^{\rm FS} =\frac2{n-1}\widehat V_{0,n}\succeq0, \] where \(\widehat{\bm \psi}_{i}=(\widehat \psi_{i}(x_1),\ldots, \widehat \psi_{i}(x_G))'\).
Sampling noise can make \(\widehat K_n^{\rm FS}\) slightly indefinite. Let \(\Pi_+(A)\) replace the negative eigenvalues of a symmetric matrix \(A\) by zero and set \(\widehat K_{n,+}^{\rm FS}=\Pi_+(\widehat K_n^{\rm FS})\). Because \(\widehat K_{n,+}^{\rm FS}\) estimates the covariance of \(\sqrt n\,\widehat{\bm R}_n\), let \(\widehat{\bm R}_{n,G}^{\rm corr,*} =(\widehat R_{n,G}^{\rm corr,*}(x_1),\ldots, \widehat R_{n,G}^{\rm corr,*}(x_G))'\) denote the unscaled corrected bootstrap process on the grid, drawn according to
Its bootstrap statistics are \[ T_{n,G}^{KS,\rm corr,*} = \max_{g\leq G} \left|\sqrt n\,\widehat R_{n,G}^{\rm corr,*}(x_g)\right|, \qquad T_{n,G}^{CvM,\rm corr,*} = n\sum_{g=1}^G w_g \{\widehat R_{n,G}^{\rm corr,*}(x_g)\}^2. \] Let \(c_{n,G,1-\alpha}^{S,\rm corr,*}\) denote the conditional \((1-\alpha)\)-quantile of \(T_{n,G}^{S,\rm corr,*}\). The corrected test compares this critical value with the same sample statistic \(T_{n,G}^S\) defined in (ref) and rejects \(H_0\) whenever \( T_{n,G}^S>c_{n,G,1-\alpha}^{S,\rm corr,*}. \) Under independent dyads, one may equivalently replace the node-scale factor \(\sqrt n\) by the dyad-scale factor \(\sqrt{N_n}\) in both the sample and bootstrap statistics. This common rescaling leaves every bootstrap comparison, critical-value decision, and Monte Carlo \(p\)-value unchanged.
\paragraph{Algorithm 2: Covariance-corrected Gaussian bootstrap.} Fix a grid \(\mathcal{X}_G=\{x_g:g\leq G\}\subset\mathcal{X}\), a statistic \(S\in\{KS,CvM\}\), a nominal level \(\alpha\), and the number of bootstrap draws \(B\). Proceed as follows.
Having stated each procedure together with its own statistic and rejection rule, we now compare their validity. Both implemented tests use a fixed deterministic grid. Under nondegenerate shared-node dependence, the raw and corrected grid tests are asymptotically equivalent. The corrected grid test also remains valid when the dyads are independent, whereas the raw grid test does not. For completeness, we additionally state the continuum validity of the raw procedure as a theoretical consequence of its full-index functional limit. The following regularity condition translates convergence of the sample and bootstrap laws into consistency of their critical values.
Assumption (ref) imposes two regime-specific requirements. First, continuity and strict increase of the limiting distribution at the relevant quantile ensure that convergence of the bootstrap law translates into convergence of the bootstrap critical value and asymptotically exact size. Second, the grid conditions ensure that the selected evaluation points capture nontrivial sampling variation, so that the limiting grid-based KS and CvM statistics are not degenerate. The continuum condition is used only for the auxiliary full-index raw result under scenario (i); under scenario (ii), where \(\mathbb G_R\equiv0\), only the grid conditions for \(\mathbb G_D\) are imposed. All implementation results rely on the grid conditions. The next theorem establishes asymptotic validity.
Theorem (ref) shows that the raw and corrected procedures are asymptotically equivalent under nondegenerate shared-node dependence, where the node component determines the limiting law. Under independent dyads, however, the raw node multiplier counts each dyad through both endpoints and therefore produces twice the correct limiting covariance. The corrected procedure removes this duplication and yields a valid grid test in both regimes. Although the appropriate rate changes from \(\sqrt n\) to \(\sqrt{N_n}\) under independence, this common rescaling of the sample and bootstrap statistics does not affect rejection decisions or bootstrap \(p\)-values.
Because the centered bootstrap critical values are \(O_p(1)\) under the corresponding normalization, Theorem (ref) and the bootstrap results above imply consistency against fixed alternatives. Fixed-alternative consistency does not, however, describe the ability of the tests to detect departures that shrink with the sample size. The relevant rate is \(s_n^{-1}\): \(n^{-1/2}\) under nondegenerate shared-node dependence and \(N_n^{-1/2}\) under independent dyads. Accordingly, write \[ s_n=\sqrt n\quad\text{in scenario (i)},\qquad s_n=\sqrt{N_n}\quad\text{in scenario (ii)}. \] Consider
where \(\Delta_0\) is square integrable.
Part of \(\Delta_0(X)\) may lie in the linear span of \(X\) and is therefore absorbed by re-estimation of the linear projection coefficient. Define \[ b_\Delta = Q^{-1}E[X_{12}\Delta_0(X_{12})], \qquad \widetilde\Delta_0(X) = \Delta_0(X)-X'b_\Delta. \] The component of the local departure that remains visible to the test is
Thus \(\mu_\Delta\) is the local departure after projection onto the orthogonalized instrument class.
Let \(c_{1-\alpha}^{KS}\) and \(c_{1-\alpha}^{CvM}\) denote the \((1-\alpha)\)-quantiles of \(\|\mathbb G_R\|_\infty\) and \(\int_{\mathcal{X}}\mathbb G_R(x)^2d\nu(x)\), respectively. Let \(c_{G,1-\alpha}^{KS,R}\) and \(c_{G,1-\alpha}^{CvM,R}\) denote the \((1-\alpha)\)-quantiles of \(\max_{g\leq G}|\mathbb G_R(x_g)|\) and \(\sum_{g=1}^G w_g\mathbb G_R(x_g)^2\), respectively. Define \(c_{G,1-\alpha}^{KS,D}\) and \(c_{G,1-\alpha}^{CvM,D}\) analogously using \(\mathbb G_D\).
Theorem (ref) shows that a departure of order \(s_n^{-1}\) enters the limiting sample process through the deterministic drift \(\mu_\Delta\), without changing the first-order covariance kernel of the centered process. Thus the shifted limit is \(\mathbb G_R+\mu_\Delta\) under shared-node dependence and \(\mathbb G_D+\mu_\Delta\) under independent dyads. Both bootstraps remain centered; the corrected bootstrap estimates the corresponding null Gaussian law in both scenarios, while the raw bootstrap does so only in scenario (i). The continuum raw results in part (a) describe the full-index theory, whereas all fixed-grid results correspond directly to implementation.
The limiting rejection probabilities of the raw and corrected tests under scenario (i), and of the corrected tests under scenario (ii), are at least \(\alpha\). For the fixed-grid tests, if the covariance matrix of the relevant Gaussian vector is positive definite, \[ (\mu_\Delta(x_1),\ldots,\mu_\Delta(x_G))'\neq0, \] and \(w_g>0\) for every \(g\) in the CvM case, then the corresponding limiting rejection probabilities are strictly greater than \(\alpha\). Under scenario (ii), when \(\mu_\Delta=0\), the limiting rejection probabilities of the raw grid tests are strictly below \(\alpha\). A nonzero drift raises rejection probability relative to this conservative null limit, but the limiting local power need not exceed \(\alpha\) when the drift is small. This variance inflation explains the loss of local power of the raw tests under independent dyads.
This section examines the finite-sample size and local power of the proposed specification tests. We compare three bootstrap procedures: the covariance-corrected Gaussian bootstrap, the raw node-multiplier bootstrap, and a naive dyad-level multiplier bootstrap that treats all dyads as independent. For each bootstrap procedure, we consider both the KS and CvM statistics.
For each node \(i=1,\ldots,n\) and regressor \(k=1,2\), generate mutually independent random variables \(U_{k,i}^W\sim U(-1,1)\), and independently generate \(U^\varepsilon_i\sim N(0,1)\). For each dyad \(i<j\), independently generate \(U^W_{k,ij}\sim U(-1,1)\), \(k=1,2\), and \(U^\varepsilon_{ij}\sim N(0,1)\), independently of all node-level variables. The regressors and regression disturbance are constructed as
where \(\omega\geq0\) controls the strength of shared-node dependence, and all node- and dyad-level components are mutually independent across \(k=1,2\). Because \(U^\varepsilon_i\), \(U^\varepsilon_j\), and \(U^\varepsilon_{ij}\) are independent of the regressors and have mean zero, the conditional mean restriction \(E(\varepsilon_{ij}\mid W_{ij})=0\) continues to hold.
We consider two data-generating processes:
DGP 1 therefore introduces an omitted quadratic term, whereas DGP 2 introduces an omitted interaction between the two regressors. The nonlinear components are orthogonal to the regressors included in the fitted linear model. In DGP 1, symmetry implies \(E(W_{1,ij})=E(W_{1,ij}^{3})=0\), so \(W_{1,ij}^{2}-E(W_{1,ij}^{2})\) is orthogonal to both the intercept and \(W_{1,ij}\). In DGP 2, independence and centering imply that \(W_{1,ij}W_{2,ij}\) is orthogonal to the intercept, \(W_{1,ij}\), and \(W_{2,ij}\). Consequently, these departures are not absorbed by re-estimation of the linear regression coefficients.
For the size experiments, we impose the null by setting \(\gamma_n=0\). We first vary the number of nodes over \[ n\in\{10,12,15,20,25,30,40,50,70,100\} \] while holding \(\omega=1\). We then fix \(n=50\) and vary the strength of dyadic dependence over \[ \omega\in\{0,0.2,0.4,\ldots,2\}. \] When \(\omega=0\), the complete dyadic observations are mutually independent. Positive values of \(\omega\) induce dependence between dyads sharing a node, with larger values representing stronger shared-node dependence.
For the power experiments, we consider the local alternatives
Thus, \(h=0\) corresponds to the null, while larger values of \(h\) represent increasingly pronounced local departures from linearity. We fix \(n=50\) and \(\omega=1\). For DGP 1, we use \(h\in\{0,0.15,\ldots,1.5\}\), while for DGP 2 we use \(h\in\{0,0.4,\ldots,4\}\). The different grids account for the different scales of the quadratic and interaction departures.
All results are based on \(10,000\) Monte Carlo replications and \(B=399\) bootstrap repetitions. The nominal significance level is \(5\%\). All procedures use the same prespecified evaluation grid. The raw procedure uses the node-multiplier bootstrap in Algorithm 1; the corrected procedure uses the bootstrap procedure in Algorithm 2; and the naive procedure attaches independent multipliers directly to the \(\binom n2\) dyads.
Figure (ref) reports rejection probabilities under the null. Panels (a) and (c) vary \(n\) under DGPs 1 and 2, respectively, whereas panels (b) and (d) vary the dependence parameter \(\omega\).
Several clear patterns emerge. First, the corrected KS test provides the most stable size control across the two DGPs and dependence regimes. Under DGP 1, its rejection probability remains close to the nominal \(5\%\) level over the full range of \(n\) and under weak or moderate dyadic dependence. Its rejection probability increases only modestly as \(\omega\) becomes large. Under DGP 2, the corrected KS test exhibits some overrejection, especially for small \(n\), but its size distortion remains substantially smaller than that of the naive tests.
The corrected CvM test generally overrejects more than the corrected KS test. Its rejection probability is approximately \(7.5\%\)-\(10\%\) in many designs and increases further under strong shared-node dependence. This pattern suggests that the integrated CvM functional may be more sensitive to finite-sample covariance-estimation error than the supremum-based KS functional.
Second, the raw node tests are severely conservative when dyadic dependence is absent or weak. At \(\omega=0\), their rejection probabilities are close to zero. This behavior agrees with the factor-of-two variance discrepancy established in (ref): when the dyads are independent, the raw node-multiplier bootstrap overestimates the sampling variance and consequently produces critical values that are too large. The distortion declines as shared-node dependence becomes stronger because the first-order node component becomes increasingly important.
An interesting exception occurs for the raw CvM test under DGP 2. When \(\omega\) is large, its rejection probability moves close to the nominal level and is sometimes more accurate than that of the corrected procedures. Thus, the raw CvM test performs well in this particular strongly dependent design. Its performance is not stable across regimes, however: it remains markedly conservative under independence and weak dependence. It therefore cannot be recommended when the strength of dyadic dependence is unknown.
Third, the naive dyad-level tests are reliable only near \(\omega=0\), where the dyads are genuinely independent. Their rejection probabilities rise rapidly with \(\omega\), and the distortion becomes especially severe for the CvM statistic and under DGP 2. The distortion also becomes more visible as \(n\) increases. Treating the \(\binom n2\) dyads as independent understates the sampling variation generated by shared nodes, leading to critical values that are too small. The resulting overrejection demonstrates that the large number of dyads cannot be interpreted as an equally large number of independent observations.
Figure (ref) reports rejection probabilities under the local alternatives \(\gamma_n=h/\sqrt n\). Because the values at \(h=0\) reproduce the finite-sample size of each procedure, power comparisons must be interpreted together with the size results in Figure (ref).
All six tests exhibit increasing rejection probabilities as \(h\) increases, confirming that the residual-marked process detects both the omitted quadratic term in DGP 1 and the omitted interaction in DGP 2. The corrected procedures have nontrivial power against departures of order \(n^{-1/2}\), in agreement with Theorem (ref).
Within each bootstrap procedure, the CvM statistic generally rejects more frequently than the corresponding KS statistic. In particular, the corrected CvM test has somewhat higher raw power than the corrected KS test in both DGPs. Part of this difference, however, reflects the larger null rejection probability of the corrected CvM test. The corrected KS test starts considerably closer to the nominal level and nevertheless develops power rapidly as \(h\) increases. It therefore provides a more favorable balance between size accuracy and power.
The raw node tests have substantially lower rejection probabilities for small and moderate values of \(h\). This should not be interpreted as evidence that they are intrinsically less sensitive to the alternatives. Rather, their low power largely reflects their severe underrejection under the null. As the departure becomes sufficiently large, their rejection probabilities eventually approach one, but they require a larger value of \(h\) than the corrected procedures.
The naive dyad tests display the highest unadjusted rejection probabilities in much of Figure (ref). These curves do not represent valid power gains because the same procedures already overreject strongly at \(h=0\). Their apparent advantage is therefore largely generated by underestimated critical values rather than by superior detection of nonlinear alternatives.
Overall, the corrected KS test delivers the best finite-sample performance across the designs considered here. It is valid at the independent-dyad boundary, remains relatively well sized as shared-node dependence strengthens, and retains substantial local power against both quadratic and interaction alternatives. The raw CvM test can provide particularly accurate size under strong dyadic dependence in DGP 2, but this advantage is confined to that regime and is accompanied by severe conservativeness when dependence is weak. When the dependence regime is not known in advance, the corrected KS test is consequently the most reliable default specification test.
A central and extensively debated question in the network literature is why individuals form professional and social connections. A prominent explanation is homophily: individuals with similar demographic, professional, or organizational characteristics tend to interact more frequently; see, among others, mcpherson2001birds, jackson2008social, and graham2017econometric. In professional networks, common office locations, practice areas, career status, and educational backgrounds may facilitate communication and cooperation. At the same time, the relationship between these characteristics and link formation need not be additive. For example, working in the same office may be particularly important for lawyers in the same practice area, while differences in age or seniority may have different implications across organizational groups. This creates a direct specification question: can the probability of a professional connection be adequately represented by an additive linear conditional mean, or are nonlinearities and interactions required?
We use the well-known Lazega law-firm network data collected from a corporate law firm in the northeastern United States between 1988 and 1991; see lazega2001collegial. The data contain information on \(n=71\) lawyers, including partners and associates, and are publicly available through the lazegalaw data set in the amen R package.\footnote{The data and variable documentation are available at \url{https://pdhoff.github.io/amen/reference/lazegalaw.html}. An alternative description of the original data is available at \url{https://www.stats.ox.ac.uk/ snijders/siena/Lazega_lawyers_data.htm}.} Although the source data are recorded in directed form, for the professional-connection outcome used here we verified that \(Y_{ij}=Y_{ji}\) for every pair. We therefore retain one binary observation for each unordered pair. The resulting undirected network contains \( N_n=\binom{71}{2}=2,485 \) unordered dyads.
Let \(Y_{ij}\) denote the observed professional connection between each pair \(\{(i,j):i<j\}\), where \(Y_{ij}=1\) indicates that the two lawyers are connected and \(Y_{ij}=0\) indicates that they are not connected. Because the outcome is binary, \[ E[Y_{ij}\mid X_{ij}] = \Pr(Y_{ij}=1\mid X_{ij}), \] so the conditional-mean model directly describes the probability of a professional connection.
The baseline dyadic regressors include the sum of the two lawyers' seniority levels, their absolute seniority difference, their absolute age difference, and indicators for whether they have the same office, practice area, professional status, gender, and law-school category. Thus, the regressors capture both the overall experience of a pair and several dimensions of professional, organizational, and demographic similarity.
We consider three increasingly flexible specifications. The first is the additive linear model \[ E[Y_{ij}\mid X_{ij}] = X_{ij}'\beta_0. \] The second augments the baseline model with the squared seniority-gap term, allowing the connection probability to vary nonlinearly with the seniority difference. The third specification additionally includes the interaction between shared office and shared practice: \[ E[Y_{ij}\mid X_{ij}] = X_{ij}'\beta_0+\text{seniority\_gap}_{ij}^2 \cdot\gamma_0+\text{shared\_office}_{ij}\times\text{shared\_practice}_{ij}\cdot\delta_0. \] This interaction allows the association between organizational proximity and connectivity to depend on whether the lawyers share a practice area.
For each specification, we compute the grid-based KS and CvM tests using \(9,999\) bootstrap repetitions. The raw and corrected procedures are evaluated on the same fixed grid. We report results from three resampling procedures. The covariance-corrected Gaussian bootstrap is our preferred procedure because it accommodates both nondegenerate shared-node dependence and the independent-dyad case. The raw node-multiplier bootstrap is valid under nondegenerate shared-node dependence, while the naive independent-dyad bootstrap is reported only as a benchmark because it ignores dependence between dyads sharing the same lawyer.
The results provide strong evidence against the additive linear specification. The corrected-bootstrap \(p\)-values are \(0.010\) for the KS statistic and \(0.002\) for the CvM statistic. Thus, both tests reject the null hypothesis that the conditional probability of a connection is linear and additive in the included dyadic characteristics at the \(5\%\) significance level. The rejection implies that the linear model leaves a systematic component of professional connectivity unexplained.
Adding the squared seniority-gap term does not resolve the misspecification. Under the quadratic model, the corrected KS and CvM \(p\)-values are \(0.004\) and \(0.002\), respectively. Both tests continue to reject strongly. Consequently, the failure of the linear model cannot be explained solely by simple curvature in the seniority difference. This finding suggests that an additive model remains too restrictive even after allowing this effect to be nonlinear.
The conclusion changes substantially when the shared-office-by-shared-practice interaction is included. For the quadratic-plus-interaction specification, the corrected \(p\)-values increase to \(0.733\) for KS and \(0.306\) for CvM. Neither statistic rejects at conventional significance levels. The raw node bootstrap gives the same qualitative conclusion, with \(p\)-values of \(0.974\) and \(0.727\). Thus, after allowing for this interaction, the residual-marked process contains no statistically detectable systematic departure from the proposed conditional mean.
Economically, the results indicate that professional connections cannot be adequately described by adding the separate contributions of experience, demographic similarity, and organizational proximity. Instead, these characteristics appear to operate jointly. The relevance of one source of similarity depends on other characteristics of the lawyer pair. For example, organizational proximity may be more strongly associated with connectivity for lawyers sharing a practice area or professional status. The specification tests therefore favor a model based on complementarities among pair characteristics over a purely additive homophily model.
Failure to reject the quadratic-plus-interaction model does not prove that it is the unique correct model. It means that, relative to the omnibus collection of lower-orthant moment restrictions considered by the test, the data provide no statistically significant evidence of remaining conditional-mean misspecification. In contrast, the linear and quadratic models are clearly rejected by the same restrictions.
Finally, the comparison of bootstrap methods demonstrates the empirical importance of accounting for dyadic dependence. For the richest specification, the naive KS test does not reject, whereas the naive CvM test produces a \(p\)-value of \(0.030\) and rejects at the \(5\%\) level. This rejection is not supported by either the corrected or raw node procedures. Because multiple dyads contain the same lawyer, treating all \(2,485\) dyads as independent overstates the effective amount of independent information and can produce misleading inference. The disagreement in the last row therefore provides a concrete illustration of why a dyadic-dependence-robust specification test is needed.
This paper develops omnibus specification tests for linear conditional-mean models with undirected dyadic data. The tests use an orthogonalized residual-marked empirical process indexed by lower orthants of the regressor distribution. We establish a uniform projection result showing that, under nondegenerate shared-node dependence, the process is asymptotically governed by its latent node projections. We also study the independent-dyad benchmark, where the node projection vanishes and the relevant convergence rate increases from \(\sqrt n\) to \(\sqrt{N_n}\).
The raw node-multiplier bootstrap is valid in the nondegenerate regime but counts the same-dyad variation twice under independence. An exact covariance decomposition motivates a corrected Gaussian bootstrap that retains the shared-node component while removing the additional same-dyad contribution. Both procedures are implemented on a common finite grid. The raw process also admits a full-index formulation used to establish the functional theory, whereas the corrected procedure is constructed directly on the grid. The resulting corrected grid tests are valid in both regimes, consistent against fixed alternatives, and have nontrivial power against rate-appropriate local alternatives.
The simulations show that the corrected KS test provides the most stable size control while retaining substantial local power. In the Lazega law-firm network, the tests reject additive linear and quadratic specifications but do not reject a richer model containing a squared seniority-gap term and a shared-office-by-shared-practice interaction. More general totally degenerate exchangeable arrays may contain a leading Gaussian-chaos component and therefore fall outside the scope of the proposed Gaussian procedures.