Multiway empirical likelihood

Multiway empirical likelihood

\address{Department of Economics, University of Wisconsin-Madison, 1180 Observatory Drive Madison, WI 53706-1393, USA.} \email{[email removed]} \address{Graduate School of Economics, Hitotsubashi University, 2-1 Naka, Kunitachi, Tokyo 186-8601, Japan.} \email{[email removed]} \address{Department of Economics, London School of Economics, Houghton Street, London, WC2A 2AE, UK.} \email{[email removed]}

\allowdisplaybreaks

Introduction

Many important econometric problems feature multiway data, in which observations are indexed by multiple sets of entities, often arranged as rows and columns. Examples include longitudinal data liang1986longitudinal, classical random effect models searle2009variance, row-column exchangeable models Mccullagh2000, nonnested multilevel data miglioretti2007marginal, bipartite networks choi2014co, and multiway clustering CGM2011, to list a few. Observations in such datasets, when correspond to a same set of entity, can exhibit strong dependence that does not diminish as certain distance measure increases, which invalidates conventional asymptotic theory. Although Eicker-White type multiway cluster robust standard errors have been developed for statistical inference on multiway data, and are frequently used in empirical research,\footnote{For example, according to Google Scholar, as of 13th of October, 2021, CGM2011 receives over $3,100$ citations, while thompson2011simple has over $1,400$ citations.} (i) their derivations are largely case-by-case, (ii) the resulting inference may not be reliable in finite samples, especially when one of the index dimensions contains only a moderate number of units, (iii) in situations with weak or no cluster dependence, the resulting inference often demonstrates significantly less precision than in the cases with strong dependence, and (iv) they often underestimate the variance in finite samples and lead to distortions in the size or coverage properties.

This paper develops a general framework to conduct inference on statistical models for various multiway data. In particular, inspired by the idea of jackknife empirical likelihood of JingYuanZhang2009, we propose a novel multiway empirical likelihood (MEL). Unlike the conventional leave-one-out operation for jackknifing, one leaves all the observations in a column or a row out at a time to construct the leave-out pseudo values. The resulting MEL function is computationally attractive and is shown to be asymptotically pivotal when the linear terms of the Hoeffding type decomposition of the statistical object of interest dominate the quadratic terms (called the non-degenerate case). In multiway data and models, however, degeneracy often occurs so that the quadratic terms in the Hoeffding type decomposition emerge in the first-order, and the MEL statistic loses its asymptotic pivotalness. This phenomenon can be understood as an analogy of emergence of the Efron-Stein bias for the jackknife variance estimator in the multiway context even though in the original setup of EfronStein1981, the bias is of second-order. To recover asymptotic pivotalness, we modify the baseline MEL statistic by incorporating leave one column and row out adjustments, which may be considered as an extension of the leave-two-out bias correction idea in hinkley1978improving and EfronStein1981 to our two-way setup. Under mild regularity conditions, this modified MEL statistic converges to a chi-square distribution regardless of the degeneracy.

To further motivate our modified MEL approach, we investigate higher-order properties of the modified MEL statistic in a simplified setup. In this analysis, we find that the second-order term of the asymptotic expansion for the modified MEL statistic is closer to zero than that of the Wald statistic with the Eicker-White type robust standard error, which is always negative. Although the results are derived under a restrictive setting, our higher-order analysis illustrates an advantage of our modified MEL inference and also provides some explanation on the oversize phenomenon of the Wald test based on Eicker-White cluster robust standard errors that has been well-documented in the literature CGM2011,thompson2011simple,mackinnon2021wild.

We illustrate wide-applicability of the modified MEL method by various econometric applications, including sparse bipartite network formation models BickelChenLevina2011, Graham2020logit, and generalized estimating equations liang1986longitudinal, xie2003asymptotics, balan2005asymptotic. Throughout these different contexts, the proposed methodology can be applied without modification. We also generalize the modified MEL method to a three-way index setting. Finally, we conduct simulation studies over various settings for classical random effect models and bipartite stochastic block models. The results suggest that the finite sample performance of the modified MEL significantly dominates the Eicker-White multiway cluster robust standard error. The difference is especially profound when cluster dependence is weak or absent. Furthermore, in contrast to the Eicker-White procedure, the modified MEL delivers reliable coverage probabilities even when one of the index dimensions contains only a moderate number of observations.

This paper also contributes to the literature of empirical likelihood (Owen1988; see Owen2001Book for an overview). After the seminal work by JingYuanZhang2009, jackknife empirical likelihood and its variants have been extended to various econometric and statistical problems, e.g., GongPengQi2010JMA, zhang2013empirical, zhong2014jackknife, among others. In particular, MatsushitaOtsu2020 proposed modified jackknife empirical likelihood to cope with the Efron-Stein bias and established its asymptotic pivotalness under both conventional and non-standard asymptotics. Under the conventional asymptotics, empirical likelihood inference has been studied and extended to various contexts; see e.g., bertail2006empirical, zhu2006empirical, HjortMcKeagueVan_Keilegom2009AoS, BravoEscancianoVan_Keilegom2020AoS, and a review by ChenVan_Keilegom2009Test, among many others.

This paper is organized as follows. Section (ref) presents our basic theoretical results on the MEL and its modification for inference on the means of two-way data. Sections (ref) and (ref) study the first and higher-order asymptotic properties, respectively. In Section (ref), we extend our MEL approach to a bipartite network model (Section (ref)), generalized estimating equations (Section (ref)), and three-way data (Section (ref)). Section (ref) illustrates the proposed method by two simulation examples. All proofs are contained in the Appendix.

Benchmark case: Two-way empirical likelihood for mean

First-order asymptotic theory

As a benchmark, for each $N,M\in \mathbbm N$, we first consider a two-way sample of $d$-dimensional random vectors $\{X_{ij}:i=1,\ldots,N,j=1,\ldots,M\}$ generated by

for $i=1,\ldots,N$ and $j=1,\ldots,M$, where $\{U_{i0},U_{0j},U_{ij}:i=1,\ldots,N,j=1,\ldots,M\}$ are i.i.d. unobservable latent shocks that can be normalized to $U[0,1]$,\footnote{Independence over $(U_{0j})_{j\in \mathbb{N}}$ is assumed for simplicity, and can be replaced by a martingale difference sequence type condition.} and $\tau=\tau_{(N,M)}$ is an unknown $\mathbb R^d$-valued Borel-measurable map that may vary with $(N,M)$. Hereafter all population objects, such as the distribution and moments of $X_{ij}$, depend on $(N,M)$ through $\tau(\cdot)$, but we suppress the dependence on $(N,M)$ for notational brevity. In this benchmark setup, we consider statistical inference on the mean vector $\theta=E[X_{11}]$ by using the MEL method.

A common set of sufficient conditions for the representation in ((ref)) is that $\{X_{ij}:i=1,\ldots,N,j=1,\ldots,M\}$ is dissociated and is embedded into an infinite two-way separately exchangeable array $(X_{ij})_{(i,j)\in\mathbb{N}^{2}}$. A set of random variables $\{X_{ij}: =1,...,N, j=1,...M\}$ is dissociated if for any two disjoint sets $A,B\subset \{1,...,N\}\times \{1,...,M\}$, $\{X_{ij}:(i,j)\in A\}$ and $\{X_{ij}:(i,j)\in B\}$ are independent. An infinite array $(X_{ij})_{(i,j)\in\mathbb{N}^{2}}$ is called separately exchangeable if for any two permutations of positive integers $\pi_{1},\pi_{2}:\mathbb{N}\to\mathbb{N}$ and any finite subset $A\subset\mathbb{N}^{2}$, it holds $(X_{ij})_{(i,j)\in A}\overset{d}{=}(X_{\pi_{1}(i)\pi_{2}(j)})_{(i,j)\in A}$. Under such conditions, ((ref)) follows from the celebrated Aldous-Hoover-Kallenberg representation for separately exchangeable arrays (e.g., Corollary 7.23 of Kallenberg2006). This representation is widely used in modern statistics, e.g., diaconis2008graph, BickelChenLevina2011, choi2014co, bhattacharyya2015subsampling, gao2015rate, caron2017sparse, choi2017co, zhang2017estimating, lauritzen2018random, veitch2019sampling, davezies2021, mackinnon2021wild, menzel2021bootstrap, and many more. See also the review by orbanz2014bayesian.

The representation in ((ref)) is useful for our theoretical development since it allows us to establish a Hoeffding type decomposition for the estimation error of the point estimator $\hat{\theta}=\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M}X_{ij}$, that is

where

Observe that $\{L_{i0}:i=1,\ldots,N\}$ and $\{L_{0j}:j=1,\ldots,M\}$ are i.i.d., and $\{R_{ij}:i=1,\ldots,N,j=1,\ldots,M\}$ is also i.i.d. conditional on $\{U_{i0},U_{0j}:i=1,\ldots,N,j=1,\ldots,M\}$.

Throughout this paper, we regard $S(\theta)=\hat{\theta}-\theta$ as an estimating equation for $\theta$ and construct the MEL function to conduct inference on $\theta$. More precisely, let us define $n=N+M$ and introduce the leave-out pseudo value: \[ V_{l}(\theta)=nS(\theta)-(n-1)S_{l}(\theta), \] for $l=1,\ldots,n$, where $S_{l}(\theta)=\hat{\theta}^{(l)}-\theta$ and $\hat{\theta}^{(l)}$ is the leave one column or row out counterpart of $\hat{\theta}$ defined as { \[ \hat{\theta}^{(l)}=\left\{

\right. \]} Unlike the conventional leave-one-out operation for jackknifing, we leave all the observations that have a specific $i$ or $j$ out at a time. Thus the number of leave-out pseudo values is $n=N+M$ instead of $NM$, which indicates computational attractiveness of our MEL method.

By considering the leave-out pseudo value $V_{l}(\theta)$ as an estimating function for $\theta$, the MEL function for $\theta$ is constructed as \[ \ell(\theta)=-2\sup_{w_{1},\ldots,w_{n}}\sum_{l=1}^{n}\log(nw_{l})\qquad\text{s.t.}\quad w_{l}\ge0,\quad\sum_{l=1}^{n}w_{l}=1,\quad\sum_{l=1}^{n}w_{l}V_{l}(\theta)=0. \] { Although this constrained optimization involves $n$ variables $\{w_{1}\ldots,w_{n}\}$, we can apply the Lagrange multiplier method to derive its dual form (see Section 3.14 of Owen, 2001), that is

where $\lambda$ is a vector of the Lagrange multipliers for the constraints $\sum_{l=1}^{n}w_{l}V_{l}(\theta)=0$, and $\lambda'$ denotes the transpose of $\lambda$. In practice $\ell(\theta)$ can be computed by the dual form since the dimension of $\lambda$ is much lower than $n$.}

To study asymptotic properties of the MEL statistic $\ell(\theta)$, we impose the following assumptions. Let $\underline{n}=N\wedge M$, $\sigma_{L}^{2}=n\{Var(L_{10})/N+Var(L_{01})/M\}$ and { $\sigma_{R}^{2}=Var(R_{11})$} be the variance matrices of the components in the Hoeffding type decomposition in ((ref)), and $\lambda_{\min}(A)$ be the minimum eigenvalue of a matrix $A$. We say the sequence of data-generating processes is {\it non-degenerate} if $\lambda_{\min}(\underline{n}\sigma_{L}^{2})\to\infty$ and {\it nearly degenerate} if $\lambda_{\min}(\underline{n}\sigma_{L}^{2})=O(1)$.

Under these assumptions, the limiting distribution of the MEL statistic $\ell(\theta)$ is obtained as follows.

In order to conduct statistical inference based on the MEL statistic $\ell(\theta)$, we need to employ different critical values for the different cases. In particular, for the nearly degenerate case, we need to estimate $Var(\xi)$ and $\Omega$. Thus, it is desirable to modify the MEL statistic to have the same limiting distribution for both cases.

Motivated by the bias correction method in EfronStein1981, we develop a modified version of the MEL statistic as follows. For each $l=1,\ldots,N$ and $l_{1}=1,\ldots,M$, let $\hat{\theta}^{(l,l_{1})}=\frac{1}{(N-1)(M-1)}\sum_{i\ne l}\sum_{j\ne l_{1}}X_{ij}$ be the leave one-column and one-row out counterparts of $\hat{\theta}$, $S_{l,l_{1}}(\theta)=\hat{\theta}^{(l,l_{1})}-\theta$, and

where $\mathcal C(N,M)=\frac{(N-1)(M-1)n}{(NM)(n-2)}$. Note that the term $Q_{ll_{1}}$ is different from the leave-two-out counterpart employed in EfronStein1981 to correct the higher-order bias of the jackknife variance estimator since we delete a whole column and row of the data matrix $(X_{ij})$. Therefore, the total number of leave-out estimators required here is of order $O(n^2)$, similar to the usual leave-one-out procedures, in contrast to ${NM \choose 2}=O(n^4)$ of the conventional leave-two-out methods. The factor $\mathcal C(N,M)$ is a finite sample adjustment to make the coefficient of the leading term $(W_{ll_{1}}+R_{ll_{1}})$ of $Q_{ll_{1}}$ to be “$\frac{n}{NM}$” as in ((ref)) in the Appendix so that the leading term of the adjustment term $\frac{MN}{n^{2}}\sum_{l=1}^{N}\sum_{l_{1}=1}^{M}Q_{ll_{1}}Q_{ll_{1}}^{\prime}$ will be $\frac{1}{NM}\sum_{l=1}^{N}\sum_{l_{1}=1}^{M}(W_{ll_1}+R_{ll_1})^{2}$, which is unbiased for $Var(W_{ll_1}+R_{ll_1})$.

Based on $Q_{ll_{1}}$, the modified estimating function is defined as \[ V_{l}^{m}(\theta)=V_{l}(\hat{\theta})-\hat{\Gamma}\tilde{\Gamma}^{-1}\{V_{l}(\hat{\theta})-V_{l}(\theta)\}, \] for $l=1,\ldots,n$, where $\hat{\Gamma}$ and $\tilde{\Gamma}$ are so that

By using $V_{l}^{m}(\theta)$ as a moment function, the modified MEL statistic is defined as

and the asymptotic property of this statistic is obtained as follows.

This theorem shows that the modified MEL statistic $\ell^{m}(\theta)$ has the asymptotically pivotal distribution of $\chi_{d}^{2}$ for both asymptotic regimes. We emphasize that the modified MEL inference only requires the estimators, $\hat{\theta}$, $\hat{\theta}^{(l)}$, and $\hat{\theta}^{(l,l_{1})}$, and circumvents estimation of $Var(\xi)$ and $\Omega$ in Theorem (ref). Based on this theorem, the asymptotic $100(1-\alpha)\%$ modified MEL confidence set can be constructed as $\{\theta:\ell^{m}(\theta)\le\chi_{d,\alpha}^{2}\}$, where $\chi_{d,\alpha}^{2}$ is the $(1-\alpha)$-th quantile of the $\chi_{d}^{2}$ distribution.

Higher-order properties

In this subsection, we provide some theoretical justification for desirable accuracy of the modified MEL statistic by the asymptotic $\chi^{2}$ calibration based on the higher-order property { in a simplified setup:

where $\{\varepsilon_{ij}\}$ is a scalar sequence of i.i.d. random variables\footnote{ As such, the results in Theorem 3 are not directly comparable with the classical Edgeworth expansion results for U-statistics, such as those in helmers1991edgeworth, putter1998empirical, etc.} for $i=1,\ldots,N$ and $j=1,\ldots,M$. In this particular scenario, our simulation studies below present substantial disparities in performance between the modified MEL method and the conventional inference procedure relying on the Eicker-White variance estimator. Although this setup covers some specific examples, such as the random effect model with $\sigma^2=0$ in the simulation study below and the Erd\H{o}s-R\'{e}nyi model, it is arguably restrictive. In the context of network data analysis, a recent paper by Zhang and Xia (2022) thoroughly studied higher-order properties of the network moments via novel Edgeworth expansion techniques. As clarified in Zhang and Xia (2022), the network moment statistics are considered as noisy U-statistics that involve edgewise errors to observe the adjacency matrix, and their higher-order analysis is substantially different from the conventional noiseless case. It is beyond the scope of this paper to extend Zhang and Xia's (2022) higher-order analysis to the modified MEL and Wald statistics under our multiway setting.}

We compare the second-order terms of the distributions of the modified MEL statistic $\ell^{m}(\theta)$ and the Wald statistic or $t$-ratio based on the Eicker-White type cluster robust variance estimator, i.e., $T(\theta)=(\hat{\theta}-\theta)^{\prime}\hat{\Sigma}^{-1}(\hat{\theta}-\theta)$, where

The variance estimator $\hat{\Sigma}$ is a two-way version of the cluster robust variance estimator of liang1986longitudinal from miglioretti2007marginal, CGM2011, and thompson2011simple. Its asymptotic properties are subsequently investigated in davezies2021 and mackinnon2021wild. In terms of the first-order asymptotic property, both $\ell^{m}(\theta)$ and $T(\theta)$ converge in distribution to the $\chi_{d}^{2}$ distribution.

Denote $\Phi$ and $\phi$ be the standard normal cumulative distribution and density functions, respectively, and $\mathrm{i}=\sqrt{-1}$. Higher-order properties of $\ell^{m}(\theta)$ and $T(\theta)$ under the simplified setup in ((ref)) are presented as follows.

Several remarks follow. First, the asymptotic expansion for the (signed root of) Wald statistic $T(\theta)$ based on the Eicker-White type cluster robust variance estimator shows that its second-order term is of order $O(n^{-1})$ and takes a negative value. This result suggests undercoverage of the Wald-type confidence interval based on $T(\theta)$ in finite samples as illustrated in our simulation studies in Section (ref).

Second, the asymptotic expansion for the modified MEL statistic $\ell^{m}(\theta)$ shows that the second-order term is also of order $O(n^{-1})$ but is closer to zero. This result indicates desirable accuracy of the modified MEL statistic by the $\chi^{2}$ calibration for the degenerate case.

Finally, even without the factor $\mathcal{C}(N,M)$ in the adjustment term in ((ref)), the corresponding modified MEL (say, $\tilde{\ell}^{m}(\theta)$) yields a smaller second-order term than that of the Wald statistic as

However, the refinement of $\tilde{\ell}^{m}(\theta)$ in the coefficient of $t\phi(t)$ is smaller than that of $\ell^{m}(\theta)$.

Generalizations

Bipartite network

In this section, we extend our (modified) MEL inference method for slope parameters in the logistic regression model for sparse bipartite network models investigated by Graham2020logit. While the asymptotic properties of the maximum composite likelihood estimator under sparse network asymptotics have been studied in the literature, no inference method has been proposed for this estimator.

Let $\mu$, $\{(W_{i0},A_{i0}):i=1,\ldots,N\}$, $\{(W_{0j},A_{0j}):j=1,\ldots,M\}$, and $\{V_{ij}:i=1,\ldots,N,j=1,\ldots,M\}$ be i.i.d. sequences, where $W_{i0}$ and $W_{0j}$ are observed attributes with supports $\mathcal{W}_{1}$ and $\mathcal{W}_{2}$, respectively, and $(\mu,A_{i0},A_{0j},V_{ij})$ are unobserved shocks. Suppose the random bipartite graph $\{Y_{ij}\in\{0,1\}:i=1,\ldots,N,j=1,\ldots,M\}$ is generated according to

where $h_{N,M}:[0,1]\times\mathcal{W}_{1}\times\mathcal{W}_{2}\times[0,1]^{3}\to\{0,1\}$ is a graphon unknown to the researcher. Suppose the researcher observes $\{Y_{ij},Z_{ij}:i=1,\ldots,N,j=1,\ldots,M\}$, where $Z_{ij}=z(W_{i0},W_{0j})$ is a vector of known transformations of $W_{i0}$ and $W_{0j}$. Consider the logistic network formation model

where $\Lambda(u)=\exp(u)/(1+\exp(u))$, and $\alpha_{n}$ is an intercept, which may vary with $n=N+M$. The asymptotics is understood as $n\to\infty$. Suppose $M/n\to\phi\in(0,1)$. Let $\rho_{n}=E[\Lambda(\alpha_{n}+Z_{11}^{\prime}\beta)]$ be the marginal link formation probability, and $\lambda_{n}^{1}=M\rho_{n}$ and $\lambda_{n}^{2}=N\rho_{n}$ be the average degrees of the first and second cluster dimensions, respectively. Suppose $\alpha_{n}=\log(\eta/n)$ for some constant $\eta$. Note that under such setting, $\alpha_{n}\to-\infty$ and $\lambda_{n}^{1}\to\lambda^{1}\in(0,\infty)$ as $n\to\infty$, that is, both average degrees stay finite in the limit.

We estimate the model by the composite maximum likelihood

where $\mathcal{L}_{ij}(\alpha,\beta)=Y_{ij}\log\Lambda(\alpha_{n}+Z_{ij}^{\prime}\beta)+(1-Y_{ij})\log(1-\Lambda(\alpha_{n}+Z_{ij}^{\prime}\beta))$. Suppose the object of interest is a $d$-dimensional subvector $\theta$ of $(\alpha,\beta^{\prime})^{\prime}$. The MEL function $\ell(\theta)$ for $\theta$ is obtained as in ((ref)) by setting $S(\theta)=\hat{\theta}-\theta$ and $S^{(l)}(\theta)=\hat{\theta}^{(l)}-\theta$, where $\hat{\theta}$ is the corresponding subvector of the estimator in ((ref)) and $\hat{\theta}^{(l)}$ is the corresponding subvector of \[ (\hat{\alpha}^{(l)},\hat{\beta}^{(l)})=\left\{

\right. \] Similarly, the modified MEL function $\ell^{m}(\theta)$ can be obtained as in ((ref)) by setting $S_{l,l_{1}}(\theta)=\hat{\theta}^{(l,l_{1})}-\theta$, where $\hat{\theta}^{(l,l_{1})}$ is the corresponding subvector of $\arg\max_{\alpha,\beta}\sum_{i\neq l}\sum_{j\neq l_{1}}\mathcal{L}_{ij}(\alpha,\beta)$.

The asymptotic property of $\ell^{m}(\theta)$ is obtained under the following conditions. Note that we do not need to impose a high-level condition that corresponds to Assumption (ref)(iii), as it holds automatically under the current setting Graham2020logit.

Similar comments to Theorem (ref) apply. As shown in Graham2020logit, the asymptotic variance of the maximum composite likelihood estimator $\hat{\theta}$ involves several terms due to non-negligible contributions from the higher-order terms in the Hoeffding type decomposition ((ref)). In contrast, our modified MEL approach only requires $\hat{\theta}$, $\hat{\theta}^{(l)}$'s, and $\hat{\theta}^{(l,l_{1})}$'s, and circumvents estimation of such variance components.

Generalized estimating equations under cluster dependence

In Section (ref), we consider inference on a logistic regression model for bipartite network data, where the modified MEL function is constructed based on the composite maximum likelihood estimator. More generally, our MEL method can be applied to conduct inference on parameters defined via generalized estimating equations (GEEs) for longitudinal data liang1986longitudinal. The existing literature on the GEE mostly focuses on the case where the cluster size is fixed and there is no dependence across clusters. A notable exception is xie2003asymptotics who investigated the asymptotic properties of the GEE estimators under the asymptotic regime of $N\to \infty$ and $M$ being either fixed or diverges to infinity at some appropriate rates while maintaining independence across clusters. Thus it is an interesting open question whether we can conduct valid inference for parameters under both growing cluster sizes and dependence across clusters.

To fix the idea, consider a generalized linear model based on the density $f(Y_{ij}|Z_{ij},\theta,\phi)=\exp[\{Y_{ij}u(Z_{ij}^{\prime}\theta)-a(u(Z_{ij}^{\prime}\theta))+b(Y_{ij})\}/\phi]$ for $i=1,\ldots,N$ and $j=1,\ldots,M$, where $u$, $a$, and $b$ are known functions and $\phi$ is a known constant. To conduct inference on $\theta$ when $M\to\infty$ and $(Y_{ij},Z_{ij})$ is embedded into a separately exchangeable array, we employ the estimating equations using the independent working correlation matrix \[ \sum_{i=1}^{N}\sum_{j=1}^{M}u^{(1)}(Z_{ij}^{\prime}\theta)Z_{ij}\{Y_{ij}-a^{(1)}(u(Z_{ij}^{\prime}\theta))\}=0, \] where $u^{(1)}$ and $a^{(1)}$ are the derivatives of $u$ and $a$, respectively. Letting $X_{ij}(\theta)=u^{(1)}(Z_{ij}^{\prime}\theta)Z_{ij}\{Y_{ij}-a^{(1)}(u(Z_{ij}^{\prime}\theta))\}$, the modified MEL function $\ell^{m}(\theta)$ is defined as in ((ref)) by setting $S(\theta)=\frac{1}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M}X_{ij}(\theta)$,

Then as far as the assumptions for Theorem (ref) are satisfied for $X_{ij}(\theta)$, we obtain $\ell^{m}(\theta)\overset{d}{\to}\chi_{\dim (\theta)}^{2}$ at the true value of $\theta$. Also the modified MEL statistic for the composite null hypothesis $H_{0}:r(\theta)=0$ can be obtained by $\min_{\theta:r(\theta)=0}\ell^{m}(\theta)$, which converges to $\chi_{\dim (r(\theta))}^{2}$ by adapting the argument in qin1994empirical.

Multiway MEL

Let us now extend the (modified) MEL approach to the three way case $\{X_{ijt}:i=1,\ldots,N,j=1,\ldots,M,t=1,\ldots,T\}$. Similar modification works for $K$-way for any fixed $K\in\mathbb{N}$. Consider the sample mean $\hat{\theta}=(NMT)^{-1}\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{t=1}^{T}X_{ijt}$ for the population mean $\theta=E[X_{111}]$. Note that following the Aldous-Hoover representation as in ((ref)) and the Hoeffding type decomposition ChiangKatoSasaki2020 for general $K$-way mean, we have

where

Suppose it holds that \[ \sum_{i=1}^{N}\sum_{j=1}^{M}A_{ij0}=\sum_{i=1}^{N}\sum_{j=1}^{M}R_{ij0}\{1+o_{p}(1)\}, \] where $R_{ij0}=g_{12}(U_{ij0})$, $g_{12}:[0,1]\to\mathbb{R}^{\dim(\theta)}$ is an unknown Borel-measurable mapping. Similarly, suppose that $\sum_{i=1}^{N}\sum_{t=1}^{T}A_{i0t}=\sum_{i=1}^{N}\sum_{t=1}^{T}\{R_{i0t}+o_{p}(1)\}$ with $R_{i0t}=g_{13}(U_{i0t})$, $\sum_{j=1}^{M}\sum_{t=1}^{T}A_{0jt}=\sum_{j=1}^{M}\sum_{t=1}^{T}\{R_{0jt}+o_{p}(1)\}$ with $R_{0jt}=g_{23}(U_{0jt})$, and $\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{t=1}^{T}A_{ijt}=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{t=1}^{T}\{R_{ijt}+o_{p}(1)\}$ with $R_{ijt}=g_{123}(U_{ijt})$ for some unknown Borel-measurable mappings $g_{13}$, $g_{23}$, and $g_{123}$. This is analogous to Assumption (ref)(iii). Only in this subsection, let $n=N+M+T$. Then a central limit theorem implies that as $\min\{N,M,T\}\to\infty$, we have $\sqrt{n}(\hat{\theta}-\theta)\overset{d}{\to}N(0,\sigma_{*}^{2})$, where \[ \sigma_{*}^{2}=n\left\{ \frac{Var(L_{100})}{N}+\frac{Var(L_{010})}{M}+\frac{Var(L_{001})}{T}+\frac{Var(R_{110})}{NM}+\frac{Var(R_{101})}{NT}+\frac{Var(R_{011})}{MT}+\frac{Var(R_{111})}{NMT}\right\}. \]

Now, define $S(\theta)=\hat{\theta}-\theta$, $S_{l}(\theta)=\hat{\theta}^{(l)}-\theta$, and the leave-one-index-out estimators \[ \hat{\theta}^{(l)}=\left\{

\right. \] Further, define $V_{l}(\theta)=nS(\theta)-(n-1)S^{(l)}(\theta)$, $V_{l}^{m}(\theta)=V_{l}(\hat{\theta})-\hat{\Gamma}\tilde{\Gamma}^{-1}\{V_{l}(\hat{\theta})-V_{l}(\theta)\},$ for $l=1,\ldots,N$, where $\hat{\Gamma}$ and $\tilde{\Gamma}$ are so that

where

and the leave-two-index-out and leave-three-index-out estimators are defined by

Under some regularity conditions, it can be shown similarly as in the proof of Theorem (ref) that

We can then obtain $\ell^{m}(\theta)$ by following the same definition as in ((ref)) with corresponding components replaced by those defined in this section. It can be shown that under regularity conditions, the modified MEL statistic has the pivotal asymptotic distribution \[ \ell^{m}(\theta)\stackrel{d}{\to}\chi_{\dim(\theta)}^{2}. \] In fact, this proposed procedure is much less computationally intensive in comparison with the corresponding jackknife procedures for $U$-statistics since the number of all leave-out estimators in the proposed procedure is $O(n^3)$, the same as the order of all different leave-one-out estimators for i.i.d. data. On the other hand, the conventional leave-three-out estimators with sample size $NMT$ consist of ${NMT \choose 3}=O(n^9)$ possibilities.

Simulation

This section conducts a simulation study to evaluate the finite sample properties of the proposed MEL inference methods. In particular, we consider a random effect model (Section (ref)) and bipartite stochastic block model (Section (ref)). We shall focus on simple means based on the reasoning as in Owen2007 that one can expect a method that gives the correct variance for a mean to be reliable for more complicated statistics such as smooth functions of means and estimating equation parameters.

Random effect model

We first consider the random effect model studied in Owen2007 and searle2009variance: \[ X_{ij}=\theta+a_{i}+b_{j}+\varepsilon_{ij}, \] where $\theta=1$ and $(a_{i},b_{j},\varepsilon_{ij})$ are mutually independent random variables with $a_{i},b_{j}\sim N(0,\sigma^{2})$ and $\varepsilon_{ij}\sim N(0,1)$. The estimator considered here is the sample mean. We vary $\sigma^{2}\in\{1,0.1,0\}$ to examine the performance under non-degenerate, nearly degenerate, and degenerate cases, respectively. We set $N=50$ and $M\in\{5,10,15,20,30,50\}$.

We compare seven methods of constructing confidence intervals: (i) multiway empirical likelihood (MEL), (ii) modified MEL (mMEL), (iii) Wald confidence interval with modified multiway variance estimator from Remark (ref) (mMW), (iv) Wald with Eicker-White type multiway cluster robust variance estimator (EWW), (v) bootstrap with model selection (MBS), (vi) conservative bootstrap (MBC), and (vii) Wald with i.i.d. variance estimator (IID). The methods (i)-(iii) are our developments, (iv) is a conventional method, (v)-(vi) are proposed by menzel2021bootstrap, and (vii) is asymptotically invalid (except for the degenerate case) but included for comparison. The nominal coverage is set as $0.95$.

Table 1 reports empirical coverages of the methods (i)-(vii) based on $5,000$ Monte Carlo replications. Our findings are summarized as follows. First, IID does not work at all except for the degenerate case (i.e., $\sigma^{2}=0$). Since IID is asymptotically invalid, its size distortion remains even for $N,M=50$. Second, EWW exhibits severe under-coverages when $M$ is small, as predicted by the higher-order analysis in Section (ref). Third, mMEL, which is asymptotically valid for all cases, outperforms in almost all cases. Even if $M$ is small, mMEL performs well. Fourth, MEL works well for non-degenerate case (i.e., $\sigma^{2}=1$) but over-covers for nearly degenerate and degenerate cases. This result is expected from Theorem (ref). As shown in Theorem (ref), mMEL recovers asymptotic pivotalness for all cases and our simulation result clearly illustrates this point. Fifth, for Wald-type confidence intervals, the proposed mMW works better than the conventional EWW but slightly under-covers. Finally, both MBS and MBC show significant over-coverages. Overall we recommend mMEL, which exhibits accurate coverages and is robust for all cases.

Bipartite stochastic block model

We next consider a stochastic block model, which is an adapted version of bhattacharyya2015subsampling for bipartite graphs with two distinctive community dimensions. First, each $i$ is randomly assigned to a membership $a\in\{1,2\}$ of the first community dimension with probabilities $\pi_{1}=(0.7,0.3)^{\prime}$ and each $j$ is randomly assigned to a membership $b\in\{1,2\}$ of the second community dimension with probabilities $\pi_{2}=(0.2,0.8)^{\prime}$. Then consider the following edge formation probabilities \[ F_{ab}=\Pr(X_{ij}=1|i\in A_a,j\in B_b)=s_{\theta}S_{ab},\text{ for }a\in\{1,2\}\text{ and }b\in\{1,2\}, \] where the blocks $A_a$ and $B_b$ satisfy $A_1\cup A_2=\{1,...,N\}$, $A_1\cap A_2=\emptyset$, $B_1\cup B_2=\{1,...,M\}$, $B_1\cap B_2=\emptyset$, $S_{ab}$'s are elements of $S=

$, and $s_{\theta}$ is chosen to satisfy $\theta=\pi_{1}^{\prime}F\pi_{2}\in\{0.5,0.1,0.05\}$.

Similar to the last subsection, we consider the seven confidence intervals (i)-(vii), as introduced in the previous subsection, for $\theta$ with the nominal coverage $0.95$ for the cases of $N=50$ and $M\in\{5,10,15,20,30,50\}$. Table 2 reports empirical coverages of the methods (i)-(vii) based on $5,000$ Monte Carlo replications. The results are qualitatively similar to the ones in the last subsection. IID has size distortions, EWW exhibits under-coverages when $M$ is small, both MBS and MBC over-covering across dense and sparse cases, and mMEL outperforms the rest for almost all cases. It is worthy of pointing out that in this setup, MEL shows over-coverages for all cases, including the relatively dense case (i.e., $\theta=0.5$). Therefore, we recommend to use mMEL for this simulation study, too.