A selective review of statistical methods using calibration information from similar studies

Abstract

In the era of big data, divide-and-conquer, parallel, and distributed inference methods have become increasingly popular. How to effectively use the calibration information from each machine in parallel computation has become a challenging task for statisticians and computer scientists. Many newly developed methods have roots in traditional statistical approaches that make use of calibration information. In this paper, we first review some classical statistical methods for using calibration information, including simple meta-analysis methods, parametric likelihood, empirical likelihood, and the generalized method of moments. We further investigate how these methods incorporate summarized or auxiliary information from previous studies, related studies, or populations. We find that the methods based on summarized data usually have little or nearly no efficiency loss compared with the corresponding methods based on all-individual data. Finally, we review some recently developed big data analysis methods including communication-efficient distributed approaches, renewal estimation, and incremental inference as examples of the latest developments in methods using calibration information.

Keywords:

• Calibration information
• empirical likelihood
• estimating equations
• generalized method of moments
• meta-analysis

1. Introduction

Statistical inference with big data can be extremely challenging owing to the high volume and large variety of observed quantities. Currently, one of the most popular approaches to this problem in statistics and computer science is the divide-and-conquer paradigm. The basic idea of this method is to break down a problem recursively into two or more sub-problems of the same or related type, such that each sub-problem becomes simple enough to be solved easily. The solution to the original problem is the optimal combination of the solutions to the sub-problems. A closely related statistical method is called parallel and distributed inference. In essence, large amounts of observed data are stored in different machines in a distributed manner. The computation is often relatively inexpensive in each machine. Then, communication is essential to enable assembly of the available results from all machines. Many related references can be found in, for example, Jordan et al. (2019). Although many new statistical methods have been developed for big data analysis, most of them have roots in traditional statistical methods of combining auxiliary information.

Combining information from similar studies has been and will continue to be an extremely important strategy in statistical inference. The most popular example of such methods is meta-analysis, in which the published results of multiple similar scientific studies are pooled to produce an enhanced estimate without using the raw individual data from each study. We refer to Borenstein et al. (2009) for a comprehensive introduction to meta-analysis. For various reasons such as privacy or capacity of computer storage, in massive data inference, only summarized data rather than the original individual data may be available. This poses a very challenging problem: how to conduct efficient updated inference by making full use of the summarized data? In recent years, many methods of combining information have been developed in economic studies, machine learning, and distributed statistical inference. The goal of this paper is to selectively review a few popular methods that are able to integrate information in different disciplines.

Utilizing external summary data or auxiliary information to obtain more accurate inference is an old and effective method in survey sampling. Owing to restrictions such as cost effectiveness or convenience, the variable of interest Y may be available for only a small portion of individuals. However, the explanatory variable X associated with Y may readily be available for all individuals. Cochran (1977) presented a comprehensive discussion on regression-type estimators making use of the summarized information from X. Chen and Qin (1993), Chen et al. (2002), and Wu and Sitter (2001) used empirical likelihood (EL; Owen, 1988) to incorporate such information in finite populations.

With advances in technology, many summarized statistical results have become available in public domains. For example, many aggregated demographic and socioeconomic status data are provided in the US census reports. The Surveillance, Epidemiology, and End Results (SEER) programme of the National Cancer Institute provides population-based cancer survival statistics such as covariate-specific survival probabilities. Imbens and Lancaster (1994) combined micro and macro data in economic studies through the generalized method of moments (GMM). Chaudhuri et al. (2008) showed that inclusion of population-level information could reduce bias and increase the efficiency of the parameter estimates in a generalized linear model setup. Wu and Thompson (2020) published an excellent monograph on combining auxiliary information in survey sampling.

In this paper, we consider two situations. In the first, the summarized information from different studies was derived using the same statistical model. Second, the summarized information was derived using statistical models that were similar but not exactly the same. In general, combining information in the former case is easier. The latter case is more complex, as one has to take into consideration the heterogeneity among different studies.

The rest of this paper is organized as follows. In Section 2, we briefly review two simple and popular meta-analysis methods for combining similar results. In Section 3, we review Owen's (1988) EL method and Qin and Lawless's (1994) over-identified parameter problem as examples of general tools for synthesizing information from summarized data. In particular, we present a new way of deriving the lower information bound for the over-identified parameter problem. Section 4 discusses enhanced inference by utilizing auxiliary information. Section 5 presents results on more flexible meta-analyses where information on different covariates are available in similar studies. Calibration of information from previous studies is described in Section 6. We discuss methods of using disease prevalence information for more efficient estimation in case–control studies in Section 7. The popular communication-efficient distributed statistical inference method used in machine learning is discussed in Section 8. Renewal estimation and incremental inference are briefly presented in Section 9. Finally, some further discussion is presented in Section 10.

2. Two simple information-combining methods

2.1. Convex combination

Suppose that ${\hat{\theta }}_1$ and ${\hat{\theta }}_2$ are two asymptotically unbiased estimators for θ from two independent studies, and that they satisfy ${\hat{\theta }}_i\sim N(\theta ,{\sigma }_i^2),i=1,2$. The most straightforward way of combining ${\hat{\theta }}_1$ and ${\hat{\theta }}_2$ is a convex combination, $\hat{\theta }=\alpha {\hat{\theta }}_1+(1-\alpha ){\hat{\theta }}_2,0<\alpha <1.$ The asymptotic variance of $\hat{\theta }$ is ${\sigma }^2={\alpha }^2{\sigma }_1^2+(1-\alpha {)}^2{\sigma }_2^2,$ which takes its minimum at $\alpha ={\sigma }_2^2/({\sigma }_1^2+{\sigma }_2^2)$. This suggests combining ${\hat{\theta }}_1$ and ${\hat{\theta }}_2$ by $\hat{\theta }=\frac{{\sigma }_2^2}{{\sigma }_1^2+{\sigma }_2^2}{\hat{\theta }}_1+\frac{{\sigma }_1^2}{{\sigma }_1^2+{\sigma }_2^2}{\hat{\theta }}_2=\frac{{\hat{\theta }}_1/{\sigma }_1^2+{\hat{\theta }}_2/{\sigma }_2^2}{1/{\sigma }_1^2+1/{\sigma }_2^2},$ an inverse-variance weighting estimator. In general, ${\sigma }_1^2$ and ${\sigma }_2^2$ are unknown; we may replace them by their estimators ${\hat{\sigma }}_1^2$ and ${\hat{\sigma }}_2^2$, respectively, which leads to $\hat{\theta }=\frac{{\hat{\sigma }}_2^2}{{\hat{\sigma }}_1^2+{\hat{\sigma }}_2^2}{\hat{\theta }}_1+\frac{{\hat{\sigma }}_1^2}{{\hat{\sigma }}_1^2+{\hat{\sigma }}_2^2}{\hat{\theta }}_2=\frac{{\hat{\theta }}_1/{\hat{\sigma }}_1^2+{\hat{\theta }}_2/{\hat{\sigma }}_2^2}{1/{\hat{\sigma }}_1^2+1/{\hat{\sigma }}_2^2}.$ As an alternative method, we may use the maximum likelihood method to argue that this is the best estimator. We can treat ${\hat{\theta }}_i$ as a direct observation from ${\hat{\theta }}_i|\theta \sim N(\theta ,{\sigma }_i^2)$, i = 1, 2. Then, the log-likelihood (regarding ${\sigma }_1^2$ and ${\sigma }_2^2$ as known constants) is $-({\hat{\theta }}_1-\theta {)}^2/(2{\sigma }_1^2)-({\hat{\theta }}_2-\theta {)}^2/\left(2{\sigma }_2^2\right).$ Maximizing this likelihood with respect to θ or setting the score function to be zero, we end up with the same inverse-variance weighting estimator.

2.2. Random-effect meta-analysis

Dersimonian and Laird (1986) proposed a moment-based estimation method using a random-effect model for meta-analysis. Let ${\hat{\theta }}_i$ be an estimator of ${\theta }_i$ from the i-th study, $i=1,2,\dots ,K$. For example, ${\hat{\theta }}_i$ could be the estimated mean response from the i-th study. When the sample size $n_i$ in the i-th study is reasonably large, we may assume that $\begin{array}{rl}& {\hat{\theta }}_i|{\theta }_i\sim N({\theta }_i,w_i^{-1}),{\theta }_i\sim N(\theta ,{\tau }^2),\\ & i=1,2,\dots ,K,\end{array}$ where the $w_i^{-1}$s are treated as known. Although the normal models hold to be true approximately, we assume that they are all true for ease of theoretical development. The goal here is to better estimate θ by combining the results from all the studies.

Unconditionally, we have ${\hat{\theta }}_i\sim N(\theta ,w_i^{-1}+{\tau }^2).$ Consider the following inverse-variance weighting estimator for θ: $\hat{\theta }=\frac{\sum _{i=1}^K{\hat{\theta }}_i{w}_i}{\sum _{i=1}^K{w}_i}$ with variance $\mathbb{V}\mathrm{a}\mathrm{r}\left(\hat{\theta }\right)=\sum _{i=1}^K{w}_i^2(w_i^{-1}+{\tau }^2)/{\left(\sum _{i=1}^K{w}_i\right)}^2.$ Define $\begin{array}{rl}Q& =\sum _{i=1}^K{w}_i({\hat{\theta }}_i-\hat{\theta }{)}^2=\sum _{i=1}^K{w}_i({\hat{\theta }}_i-\theta {)}^2\\ & -(\hat{\theta }-\theta {)}^2\sum _{i=1}^K{w}_i.\end{array}$ We can easily check that $\mathbb{E}\left(Q\right)=(K-1)+{\tau }^2\left(\sum _{i=1}^K{w}_i-\sum _{i=1}^K{w}_i^2/\sum _{j=1}^K{w}_j\right),$ which implies that a natural estimator of ${\tau }^2$ is ${\hat{\tau }}^2=\frac{Q-(K-1)}{\sum _{i=1}^K{w}_i-\sum _{i=1}^K{w}_i^2/\sum _{j=1}^K{w}_j}.$ For small sample sizes, there is no guarantee that this estimator is non-negative; one may replace it by $max({\hat{\tau }}^2,0).$

Alternatively, we may estimate τ using the likelihood approach. The joint likelihood based on the ${\hat{\theta }}_i$s is $\ell (\theta ,\tau )=-\frac{1}{2}\sum _{i=1}^K\frac{({\hat{\theta }}_i-\theta {)}^2}{{\tau }^2+w_i^{-1}}-\frac{1}{2}\sum _{i=1}^K\log ({\tau }^2+w_i^{-1}).$ Maximizing ℓ with respect to θ and ${\tau }^2$ gives their maximum likelihood estimators (MLEs).

Lin and Zeng (2010) compared the relative efficiency of using summary statistics versus individual-level data in meta-analysis. They found that in general there was no information loss when using the summarized information compared with inference based on the original individual data when available.

3. Empirical likelihood and general estimating equations

In this section we briefly review Owen's (1988) EL and Qin and Lawless' (1994) estimating equations approaches, as those methods represent general tools for assembly of information from different sources. The maximum likelihood method for regular parametric models is among the most popular methods in statistical inference, as it has many nice properties. However, model mis-specification is a major concern, as a mis-specified model may lead to biased results. For the case when the underlying distribution is multinomial, Hartely and Rao (1968) proposed a mean constrained estimator for the population total in survey sampling problems. To mimic the parametric likelihood but discard parametric model assumptions, Owen (1988) and Owen (1990) proposed the EL method, which is a natural generalization of the multinomial likelihood when the number of categories is equal to the sample size. The EL approach can be thought of as a bootstrap that does not resample, or as a likelihood without parametric assumptions (Owen, 2001).

3.1. Definition of empirical likelihood

Suppose that $X_1,\dots ,X_n$ are n independent and identically distributed observations from X, with cumulative distribution F. For convenience, we assume there are no ties, i.e., any two observations are unequal to each other. The techniques developed below can be easily adapted to handle ties. Let $dF\left(X_i\right),i=1,2,\dots ,n,$ be the jumps of $F\left(x\right)$ at the observed data points. The nonparametric likelihood is $L\left(F\right)=\prod _{i=1}^n{p}_i.$ It is clear that if any $p_i=0$, then $L\left(F\right)=0$, and if $\sum _{i=1}^n{p}_i<1$, then $L\left(F\right)<L\left(F_{\ast }\right)$, where $F_{\ast }\left(x\right)=\sum _{i=1}^n{p}_{i}I(X_i\le x)/\sum _{i=1}^n{p}_i$. According to the likelihood principle (that parameters with larger likelihoods are preferable), one need only consider the distribution functions $F\left(x\right)$ with $p_i>0$ and $\sum _{i=1}^n{p}_i=1$.

If we maximize the log-likelihood (1) $\ell \left(F\right)=\sum _{i=1}^n\log p_i$(1) subject to the constraints (2) $\sum _{i=1}^n{p}_i=1,p_i\ge 0,$(2) then we obtain $p_i=1/n,i=1,2,\dots ,n.$ Therefore, the maximum EL estimator of F is $F_n\left(x\right)=\sum _{i=1}^n{p}_{i}I(X_i\le x)=n^{-1}\sum _{i=1}^{n}I(X_i\le x).$ This is why the empirical distribution is called the nonparametric MLE of $F\left(x\right)$.

Suppose we are interested in constructing a confidence interval for $\mu =\mathbb{E}\left(X\right)=\int x\mathrm{d}F\left(x\right)$, the mean of X. Since we have discretized F at each of the observed data points, the integral becomes $\mu =\sum _{i=1}^n{p}_i{X}_i$. Next, we maximize the nonparametric log-likelihood subject to an extra constraint: (3) $\sum _{i=1}^n{p}_i(X_i-\mu )=0.$(3) Maximizing the log-likelihood (1(1) $\ell \left(F\right)=\sum _{i=1}^n\log p_i$(1) ) subject to constraints (2(2) $\sum _{i=1}^n{p}_i=1,p_i\ge 0,$(2) ) and (3(3) $\sum _{i=1}^n{p}_i(X_i-\mu )=0.$(3) ), the Lagrange multiplier method gives the profile log-likelihood of μ, (4) ${\ell }_n\left(\mu \right)=-\sum _{i=1}^n\log \{1+{\lambda }^{\mathrm{\top }}(X_i-\mu \left)\right\}-n\log n,$(4) where λ is the solution to $\sum _{i=1}^n(X_i-\mu )/\{1+{\lambda }^{\mathrm{\top }}(X_i-\mu \left)\right\}=0.$ We can treat ${\ell }_n\left(\mu \right)$ as a parametric likelihood of μ. Based on this likelihood, the maximum EL estimator of μ is $\hat{\mu }=\overline{X}=n^{-1}\sum _{i=1}^n{X}_i$, which is exactly the sample mean. We define the likelihood ratio function as $R_n\left(\mu \right)=2\left\{\underset{\mu }{max}{\ell }_n\right(\mu )-{\ell }_n(\mu \left)\right\}=2\left\{{\ell }_n\right(\overline{X})-{\ell }_n(\mu \left)\right\}.$ Under the regularity conditions specified in Owen (1988) and Owen (1990), as n goes to infinity, $R_n\left({\mu }_0\right)$ converges to the ${\chi }^2$ distribution with p degrees of freedom, where p is the dimension of μ, and ${\mu }_0$ is the true value of μ.

3.2. General estimating equations

The original EL was mainly used to make inference for linear functionals of the underlying population distribution such as the population mean (Owen, 1988, 1990). Qin and Lawless (1994) applied this method to general estimating models, which greatly broadened its applications. Specifically, suppose the population of interest satisfies a general estimating equation (5) $\mathbb{E}\left\{g\right(X,\theta \left)\right\}=0,$(5) for a $r\times 1$ vector-valued function g and some θ, which is a $p\times 1$ parameter to be estimated. We assume $r\ge p$ as otherwise the true parameter value of θ would be undefined.

For general estimating equations with r>p or over-identified models, Hansen (1982) proposed the celebrated GMM, which has become one of the most popular methods in the econometric community. In essence, the GMM minimizes ${\left\{\sum _{i=1}^{n}g\right(X_i,\theta \left)\right\}}^{\mathrm{\top }}{\mathrm{\Sigma }}^{-1}\left\{\sum _{i=1}^{n}g\right(X_i,\theta \left)\right\}$ with respect to θ, where Σ is the variance matrix of the estimating equation $g(X,\theta )$. If Σ is unknown, we may replace it by the sample variance $\hat{\mathrm{\Sigma }}=\frac{1}{n}\sum _{i=1}^{n}g(X_i,\tilde{\theta })g^{\mathrm{\top }}(X_i,\tilde{\theta }),$ where $\tilde{\theta }$ is an initial and consistent estimate of θ.

Instead of GMM, Qin and Lawless (1994) used the EL to make inferences for parameters defined by a general estimating equation. For discretized $F\left(x\right)$ satisfying (2(2) $\sum _{i=1}^n{p}_i=1,p_i\ge 0,$(2) ), Equation (5(5) $\mathbb{E}\left\{g\right(X,\theta \left)\right\}=0,$(5) ) becomes (6) $\sum _{i=1}^n{p}_{i}g(X_i,\theta )=0.$(6) Maximizing the log-likelihood (1(1) $\ell \left(F\right)=\sum _{i=1}^n\log p_i$(1) ) subject to (2(2) $\sum _{i=1}^n{p}_i=1,p_i\ge 0,$(2) ) and (6(6) $\sum _{i=1}^n{p}_{i}g(X_i,\theta )=0.$(6) ), we have the following profile log-likelihood of θ (up to a constant): ${\ell }_n\left(\theta \right)=-\sum _{i=1}^n\log \{1+{\lambda }^{\mathrm{\top }}g(X_i,\theta \left)\right\},$ where λ is the Lagrange multiplier determined by $\sum _{i=1}^{n}g(X_i,\theta )/\{1+{\lambda }^{\mathrm{\top }}g(X_i,\theta \left)\right\}=0.$ We then estimate θ by the maximizer $\hat{\theta }=\arg \underset{\theta }{max}{\ell }_n\left(\theta \right)$, whose limiting distribution is established in the following theorem. Hereafter, we use ${\mathrm{\nabla }}_{\theta }$ to denote the differentiation operator with respect to θ.

Theorem 3.1

Qin & Lawless, 1994

Denote $g=g(X,{\theta }_0)$ and ${\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g={\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g(X,{\theta }_0)$. Suppose that (1) $\mathbb{E}\left(g{g}^{\mathrm{\top }}\right)$ is positive definite, (2) ${\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g(X,\theta )$ is continuous in a neighbourhood of ${\theta }_0$, (3) $\Vert {\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g(X,\theta )\Vert$ and $\Vert g(X,\theta ){\Vert }^3$ can be bounded by some integrable function $G\left(X\right)$ in this neighbourhood, and (4) $\mathbb{E}\left({\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g\right)$ is of full rank. Then, as $n\to \mathrm{\infty }$, $\sqrt{n}(\hat{\theta }-\theta )\stackrel{d}{⟶}N(0,V),$ where $\stackrel{d}{⟶}$ means ‘convergence in distribution’ and (7) $V={\left\{\mathbb{E}\left({\mathrm{\nabla }}_{\theta }{g}^{\mathrm{\top }}\right)\right(\mathbb{E}g{g}^{\mathrm{\top }}{)}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g\right)\}}^{-1}.$(7)

3.3. Calculation of the information bound

Assuming that the parameter of interest satisfies the general estimating equation $\mathbb{E}\left\{g\right(X,\theta \left)\right\}=0$, we next consider how well we can estimate θ based on this model, and whether the maximum EL estimator is optimal. To answer these questions, we consider an ideal situation, where the probability function X has a parametric form $f(x,\theta )$, which is known up to θ. We define $\begin{array}{rl}h(x,\eta ,\theta )& =\exp \left\{{\eta }^{\mathrm{\top }}g\right(x,\theta \left)\right\}f(x,\theta )÷\\ & \int \exp \left\{{\eta }^{\mathrm{\top }}g\right(t,\theta \left)\right\}f(t,\theta )\mathrm{d}t,\end{array}$ implicitly assuming that $\int \exp \left\{{\eta }^{\mathrm{\top }}g\right(t,\theta \left)\right\}f(t,\theta )\mathrm{d}t<\mathrm{\infty }.$ Clearly, $h(x,\eta ,\theta )$ is an enlarged parametric model of $f(x,\theta )$ as it reduces to $f(x,\theta )$ when $\eta =0$. As the parametric form $f(x,\theta )$ is unknown in practice, we anticipate that any estimator based on the moment constraints $\mathbb{E}\left\{g\right(X,\theta \left)\right\}=0$ should have a variance that is no less than that of the MLE derived from the enlarged model. We show that even if the form of $f(x,\theta )$ is available, the MLE of θ based on $h(x,\eta ,\theta )$ has the same asymptotic variance as the maximum EL estimator.

With the parametric model h, we can estimate θ by maximizing $L(\theta ,\eta )=\prod _{i=1}^{n}h(X_i,\eta ,\theta )$ with respect to $(\theta ,\eta )$. We denote the resulting MLE by $(\tilde{\eta },\tilde{\theta })$. We show in Section 3.4 that under some regularity conditions on h (see, e.g., Theorems 14 and 23 of van de Vaart (2000)), as $n\to \mathrm{\infty }$, (8) $\sqrt{n}(\tilde{\theta }-\theta )\stackrel{d}{⟶}N(0,V),$(8) where V is defined in (7(7) $V={\left\{\mathbb{E}\left({\mathrm{\nabla }}_{\theta }{g}^{\mathrm{\top }}\right)\right(\mathbb{E}g{g}^{\mathrm{\top }}{)}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g\right)\}}^{-1}.$(7) ). In general, the parametric form $f(x,\theta )$ is unknown; hence, we expect that the best estimator of θ should have an asymptotic variance at least as large as V. As the maximum EL estimator of θ of Qin and Lawless (1994) has asymptotic variance V, we conclude that it achieves the lower information bound.

Remark 3.1

If $g(x,\theta )$ is an unbounded function of x for each θ, we may construct a new density $h(x,\theta ,\eta )=\frac{\psi \left\{{\eta }^{\mathrm{\top }}g\right(x,\theta \left)\right\}f(x,\theta )}{\int \psi \left\{{\eta }^{\mathrm{\top }}g\right(t,\theta \left)\right\}f(t,\theta )\text{d}t},$ where $\psi \left(x\right)=2(1+{\mathrm{e}}^{-2x}{)}^{-1}$ with $\psi \left(0\right)={\psi }^{\prime }\left(0\right)=1$. Clearly, ψ is bounded. We may go through the same derivations to get the same conclusion.

Remark 3.2

Back and Brown (1992) established a similar result by constructing an exponential family. In particular, they defined $h(x,\theta )=\exp \left\{{\xi }^{\mathrm{\top }}\right(\theta \left)g\right(x,{\theta }_0)-a(\theta \left)\right\}{f}_0\left(x\right),$ where $f_0\left(x\right)=f(x,{\theta }_0)$ and $\xi \left(\theta \right)$ is determined implicitly by the following conditions: $\xi \left({\theta }_0\right)=0$, $a\left({\theta }_0\right)=0$, $\int \exp \left\{{\xi }^{\mathrm{\top }}\right(\theta \left)g\right(x,{\theta }_0)-a(\theta \left)\right\}{f}_0\left(x\right)\mathrm{d}x=1,$ and $\int g(x,\theta )\exp \left\{{\xi }^{\mathrm{\top }}\right(\theta \left)g\right(x,{\theta }_0)-a(\theta \left)\right\}{f}_0\left(x\right)\mathrm{d}x=0.$ In Back & Brown's approach, $\xi \left(\theta \right)$ is determined implicitly by the above constraint equation, whereas in our new approach, η is an independent parameter.

3.4. A sketched proof of (8)

The log-likelihood based on the enlarged model is $\ell (\theta ,\eta )=\sum _{i=1}^n\log \left\{h\right(X_i,\eta ,\theta \left)\right\}$, where $\begin{array}{rl}\log \left\{h\right(x,\eta ,\theta \left)\right\}& ={\eta }^{\mathrm{\top }}g(x,\theta )+\log f(x,\theta )\\ & -\log [\int \exp \{{\eta }^{\mathrm{\top }}g(t,\theta )\left\}f\right(t,\theta \left)\mathrm{d}t\right].\end{array}$ If $\log \left\{h\right(x,\eta ,\theta \left)\right\}$ satisfies the conditions of Theorem 14 of van de Vaart (2000) on $m_{\theta }\left(x\right)$, then $(\tilde{\theta },\tilde{\eta })$ is consistent with $({\theta }_0,0)$.

Result (8(8) $\sqrt{n}(\tilde{\theta }-\theta )\stackrel{d}{⟶}N(0,V),$(8) ) follows from Theorem 23 of van de Vaart (2000). With tedious algebra, we find that $\begin{array}{rl}& {\mathrm{\nabla }}_{\theta }\log \left\{h\right(x,{\theta }_0,0\left)\right\}=0,\\ & {\mathrm{\nabla }}_{\eta }\log \left\{h\right(x,{\theta }_0,0\left)\right\}=g(x,{\theta }_0),\\ & \mathbb{E}{\mathrm{\nabla }}_{\theta {\theta }^{\mathrm{\top }}}\log \left\{h\right(X,{\theta }_0,0\left)\right\}=0,\\ & \mathbb{E}{\mathrm{\nabla }}_{\eta {\eta }^{\mathrm{\top }}}\log \left\{h\right(X,{\theta }_0,0\left)\right\}=-\mathbb{E}\left(g{g}^{\mathrm{\top }}\right),\\ & \mathbb{E}{\mathrm{\nabla }}_{\eta {\theta }^{\mathrm{\top }}}\log \left\{h\right(X,{\theta }_0,0\left)\right\}=\mathbb{E}{\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g(X,{\theta }_0)\\ & -\mathbb{E}\left\{{\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g\right(X,{\theta }_0)+g(X,{\theta }_0){\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}\log f(X,{\theta }_0\left)\right\}.\end{array}$ Under some mild assumptions, such as that $\int g(x,\theta )\times f(x,\theta )\mathrm{d}x=0$ holds for θ in a neighbourhood of ${\theta }_0$, differentiating both sides with respect to θ leads to $\mathbb{E}\left\{{\mathrm{\nabla }}_{\theta }g\right(X,\theta \left)\right\}+\mathbb{E}\left\{{\mathrm{\nabla }}_{\theta }g\right(X,\theta ){\mathrm{\nabla }}_{\theta }\log f(X,\theta \left)\right\}=0,$ which means $\mathbb{E}{\mathrm{\nabla }}_{\eta {\theta }^{\mathrm{\top }}}\log \left\{h\right(X,{\theta }_0,0\left)\right\}=\mathbb{E}\left\{{\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g\right(X,{\theta }_0\left)\right\}.$ As $(\tilde{\theta },\tilde{\eta })$ is consistent with $({\theta }_0,0)$, by Theorem 5.23 of van de Vaart (2000), we have (9) $\begin{array}{rl}& \sqrt{n}\left(\begin{array}{c}\tilde{\theta }-{\theta }_0\\ \tilde{\eta }-0\end{array}\right)={\left(\begin{array}{cc}0& \mathbb{E}\left({\mathrm{\nabla }}_{\theta }{g}^{\mathrm{\top }}\right)\\ \mathbb{E}\left({\mathrm{\nabla }}_{{\theta }^{\mathrm{\top }}}g\right)& -\mathbb{E}\left(g{g}^{\mathrm{\top }}\right)\end{array}\right)}^{-1}\\ & \times \left(\begin{array}{c}0\\ n^{-1/2}\sum _{i=1}^{n}g(X_i,{\theta }_0)\end{array}\right)+o_p\left(1\right).\end{array}$(9) This, together with the fact that $n^{-1/2}\sum _{i=1}^{n}g(X_i,{\theta }_0)\stackrel{d}{⟶}N(0,\mathbb{E}(g{g}^{\mathrm{\top }}\left)\right)$ as n goes to infinity, implies (8(8) $\sqrt{n}(\tilde{\theta }-\theta )\stackrel{d}{⟶}N(0,V),$(8) ).

3.5. Empirical entropy family

Again we assume that the available information is given by the estimating equation $\mathbb{E}\left\{g\right(X,\theta \left)\right\}=0$. The enlarged parametric model $h(x,\eta ,\theta )$ satisfies $\int h(x,\eta ,\theta )g(x,\theta )\mathrm{d}x=0$ only if $\eta =0$. Naturally, one may require $\eta =\eta \left(\theta \right)$ to satisfy $\int g(x,\theta )\exp \left\{{\eta }^{\mathrm{\top }}g\right(x,\theta \left)\right\}f(x,\theta )\mathrm{d}x=0.$ It is often too restrictive to assume a known underlying parametric model $f(x,\theta )$ in the construction of the enlarged parametric model $h(x,\eta ,\theta )$. We may replace the cumulative distribution function $F(x,\theta )={\int }_{-\mathrm{\infty }}^{x}f(t,\theta )\mathrm{d}t$ by the empirical distribution $F_n\left(x\right)=n^{-1}\sum _{i=1}^{n}I(X_i\le x)$. In this situation, $\eta =\eta \left(\theta \right)$ is the solution to $\sum _{i=1}^{n}g(x_i,\theta )\exp \left\{{\eta }^{\mathrm{\top }}g\right(x_i,\theta \left)\right\}=0.$

Let $H(x,\eta ,\theta )={\int }_{-\mathrm{\infty }}^{x}h(t,\eta ,\theta )\mathrm{d}t$. For fixed parameter values $(\eta ,\theta )$, the jump of H at $x=X_i$ is $\begin{array}{rl}\mathrm{d}H(X_i,\eta ,\theta )& =\exp \left\{{\eta }^{\mathrm{\top }}\right(\theta \left)g\right(X_i,\theta \left)\right\}÷\\ & [\sum _{j=1}^n\exp \{{\eta }^{\mathrm{\top }}\left(\theta \right)g(X_j,\theta )\left\}\right],\end{array}$ and the likelihood becomes $\prod _{i=1}^n\mathrm{d}H(X_i,\eta ,\theta )=\prod _{i=1}^n\frac{\exp \left\{{\eta }^{\mathrm{\top }}\right(\theta \left)g\right(X_i,\theta \left)\right\}}{\sum _{j=1}^n\exp \left\{{\eta }^{\mathrm{\top }}\right(\theta \left)g\right(X_j,\theta \left)\right\}}.$ In fact, this is equivalent to the EL $\prod _{i=1}^n{p}_i$, where the $p_i$s minimize the Kullback–Leibler divergence (up to a constant) or minus the exponential titling likelihood $\sum _{i=1}^n{p}_i\log \left(p_i\right)$ subject to the constraints $\sum _{i=1}^n{p}_i=1$, $p_i\ge 0$, and $\sum _{i=1}^n{p}_{i}g(X_i,\theta )=0.$ See Susanne (2007) for more details. We call this the empirical entropy family induced by the estimating equation $\mathbb{E}\left\{g\right(X,\theta \left)\right\}=0$.

4. Enhancing efficiency using auxiliary information

In this section, we discuss methods of incorporating auxiliary information to enhance estimation efficiency. This aspect was also investigated by Qin (2000). We assume a parametric model $f(y|x,\text{β})$ for the conditional density function of Y given X and leave the marginal distribution $G\left(x\right)$ of X unspecified. We wish to make inferences for β when some auxiliary information is summarized through an estimating equation $\mathbb{E}\left\{\varphi \right(X,\text{β}\left)\right\}=0.$ For example, if we know the mean μ of Y, then we can construct an estimating equation $\mathbb{E}(Y-\mu )=0.$ We can take $\begin{array}{rl}\varphi (X,\text{β})& =\int (y-\mu )f(y|X,\text{β})\mathrm{d}y\\ & =\int yf(y|X,\text{β})\mathrm{d}y-\mu .\end{array}$ Furthermore, we allow that the response Y may have missing values. Let D be the non-missingness indicator, which takes the value 1 if Y is available, and 0 otherwise. We assume a missing-at-random model $\begin{array}{rl}& \mathrm{p}\mathrm{r}(D=1|Y=y,X=x)=\mathrm{p}\mathrm{r}(D=1|X=x)\\ & =\pi \left(x\right),\end{array}$ where $\pi \left(x\right)$ depends only on x. We denote the observed data by $(d_i,d_i{y}_i,x_i)$ ($i=1,2,\dots ,n$) and $p_i=dG\left(x_i\right)$. The likelihood of $(\text{β},G)$ is $\begin{array}{rl}L& =\prod _{i=1}^n\left\{\pi \right(x_i\left)f\right(y_i|x_i,\text{β}\left)\mathrm{d}G\right(x_i){\}}^{d_i}\\ & \times \left[\right\{1-\pi \left(x_i\right)\left\}dG\right(x_i){]}^{1-d_i}\\ & =\prod _{j=1}^n\left\{\pi \right(x_j\left){\}}^{d_j}\right\{1-\pi \left(x_j\right){\}}^{1-d_j}\\ & \cdot \prod _{i=1}^n\left\{f\right(y_i|x_i,\text{β}){\}}^{d_i}\cdot p_i.\end{array}$ We can maximize this likelihood subject to the constraints $\sum _{i=1}^n{p}_i=1,p_i\ge 0,\sum _{i=1}^n{p}_i\varphi (x_i,\text{β})=0.$ As $\prod _{j=1}^n\left\{\pi \right(x_j\left){\}}^{d_j}\right\{1-\pi \left(x_j\right){\}}^{1-d_j}$ is not a function of β, the profile hybrid empirical log-likelihood (up to a constant) is (10) $\begin{array}{r}\ell \left(\text{β}\right)=\sum _{i=1}^n[d_i\log f(y_i|x_i,\text{β})-\log \{1+{\lambda }^{\mathrm{\top }}\varphi (x_i,\text{β})\left\}\right],\end{array}$(10) where λ is the Lagrange multiplier determined by (11) $\sum _{i=1}^n\frac{\varphi (x_i,\text{β})}{1+{\lambda }^{\mathrm{\top }}\varphi (x_i,\text{β})}=0.$(11) For the special case where data are missing completely at random, i.e., $\pi \left(x\right)$ is a constant function of x, Qin (2000) established the following theorem.

Theorem 4.1

Let ${\text{β}}_0$ be the true parameter value, let $\widehat{\text{β}}$ be the maximum hybrid EL estimator, i.e., the maximizer of (10(10) $\begin{array}{r}\ell \left(\text{β}\right)=\sum _{i=1}^n[d_i\log f(y_i|x_i,\text{β})-\log \{1+{\lambda }^{\mathrm{\top }}\varphi (x_i,\text{β})\left\}\right],\end{array}$(10) ), and let $\hat{\lambda }$ be the corresponding Lagrange multiplier. Denote $\varphi =\varphi (X,{\text{β}}_0)$, ${\mathrm{\nabla }}_{\text{β}}\varphi ={\mathrm{\nabla }}_{\text{β}}\varphi (X,{\text{β}}_0),$ and $\begin{array}{rl}J& =-\mathbb{E}\{d_i{\mathrm{\nabla }}_{\text{β}{\text{β}}^{\mathrm{\top }}}\log f(y_i|x_i,{\text{β}}_0\left)\right\}\\ & =\mathbb{V}\mathrm{a}\mathrm{r}\{d_i{\mathrm{\nabla }}_{\text{β}}\log f(y_i|x_i,{\text{β}}_0\left)\right\}.\end{array}$ Under some regularity conditions, when n goes to infinity, we have $\sqrt{n}\left(\right(\widehat{\text{β}}-{\text{β}}_0{)}^{\mathrm{\top }},{\hat{\lambda }}^{\mathrm{\top }}{)}^{\mathrm{\top }}\stackrel{d}{⟶}N(0,\mathrm{\Sigma }),$ where $\mathrm{\Sigma }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathrm{\Sigma }}_{11},{\mathrm{\Sigma }}_{22})$ with (12) $\begin{array}{rl}{\mathrm{\Sigma }}_{11}& =\{J+\mathbb{E}({\mathrm{\nabla }}_{\text{β}}{\varphi }^{\mathrm{\top }}\left)\right(\mathbb{E}\varphi {\varphi }^{\mathrm{\top }}{)}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}\varphi \right){\}}^{-1},\\ {\mathrm{\Sigma }}_{22}& =\left\{\mathbb{E}\right({\mathrm{\nabla }}_{\text{β}}{\varphi }^{\mathrm{\top }}\left)J^{-1}\mathbb{E}\right({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}\varphi )+\mathbb{E}(\varphi {\varphi }^{\mathrm{\top }}){\}}^{-1}.\end{array}$(12)

Remark 4.1

Imbens and Lancaster (1994) studied the same problem using GMM. In particular, they directly combined the conditional score estimating equation ${\mathrm{\nabla }}_{\text{β}}\log f(y|x,\text{β})$ and $\varphi (x,\text{β})$. Even though the first-order large-sample results are the same, the hybrid EL based approach is more appealing as it respects the parametric conditional likelihood and replaces only the marginal likelihood with the EL. See Qin (2000) for numerical comparisons of results of the two methods.

5. Combining summary information: a more flexible method for meta-analysis

Developing systematic methods for combining published information is one of the main goals of meta-analysis, which has become increasingly popular since little extra cost is needed. The main restriction in meta-analysis is that all studies must include the same variables in their analyses. The only difference allowed is in the sample sizes. Thus, studies must be discarded if they contain different variables from those in other studies.

Summarized information is often available from publications such as census reports and results of national health studies. For reasons including confidentiality, it is typically not possible to gain access to the original data, only the summarized reports. Suppose we are interested in conducting a new study that may contain some new variables of interest that are not available in the summarized information, for example, a genetic study involving newly discovered biomarkers or genes. Below we discuss a more flexible method that could be used to combine published information and individual study data for enhanced inference in such cases. Chatterjee et al. (2016) discussed a related problem on the utilization of auxiliary information. As Han and Lawless (2016) pointed out, however, their methodology and theoretical results had already been developed by Imbens and Lancaster (1994) and Qin (2000) in the absence of selection bias in sampling.

We consider two cases. (I) The sample size for the summarized information is much larger than that of the new study. (II) Sample sizes from the two data sources are comparable. In Case I, we can treat the summarized information as known, i.e., the variation in the summarized data is negligible compared with the variation in the new study. In Case II, we have to take the variation in the summarized information into consideration as it is comparable to the variation in the new study. We focus on Case I in this section and study Case II in Section 6.

5.1. Setup and solution

Suppose that the summarized results were obtained from statistical analyses of response Y and covariate variables X (although the original data are not available), and that the new study includes an extra covariate Z in addition to $(Y,X)$. We are interested in fitting a parametric model $f(y|x,z,\text{β})$ for the conditional density function of Y given X and Z. Let $(y_1^{\ast },x_1^{\ast }),\dots .,(y_N^{\ast },x_N^{\ast })$ be the historic data even though they are unavailable. The published information can be summarized in two ways:

1. $\overline{h}=N^{-1}\sum _{i=1}^{N}h(y_i^{\ast },x_i^{\ast })$ is known; and

2. ${\gamma }^{\ast }$ is the solution of an estimating equation $\sum _{i=1}^{N}h(y_i^{\ast },x_i^{\ast },\gamma )=0,$ where the function $h(y,x,\gamma )$ is known up to γ.

Let $(y_1,x_1,z_1),\dots ,(y_n,x_n,z_n)$ be observed data from the new study. The basic assumption is that $(y_i,x_i),i=1,2,\dots ,n$, and $(y_i^{\ast },x_i^{\ast })$ have the same distribution. To utilize the summarized information, we can define estimating functions $\begin{array}{c}g=(g_1,g_3),g_1(y,x,z)={\mathrm{\nabla }}_{\text{β}}\log f(y|x,z,\text{β}),\\ g_3(y,x)=h(y,x)-\overline{h}\end{array}$ in Scenario (I), and $\begin{array}{rl}g& =(g_1,g_3),\\ g_1(y,x,z)& ={\mathrm{\nabla }}_{\text{β}}\log f(y|x,z,\text{β}),\\ g_3(y,x)& =h(y,x,{\gamma }^{\ast })\end{array}$ in Scenario (II). We consider only the situation where $n/N\to 0$. In other words, the variation in the auxiliary information is negligible.

The EL approach amounts to maximizing $\sum _{i=1}^n\log p_i$ subject to the constraint $\sum _{i=1}^n{p}_{i}g(y_i,x_i,z_i,\text{β})=0,p_i\ge 0,\sum _{i=1}^n{p}_i=1.$ According to Qin and Lawless (1994), the asymptotic variance of the maximum EL estimator $\widehat{\text{β}}$ based on estimating equation g is $\left[\mathbb{E}\right({\mathrm{\nabla }}_{\text{β}}{g}^{\mathrm{\top }}\left)\right\{\mathbb{E}\left(g{g}^{\mathrm{\top }}\right){\}}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}g\right){]}^{-1},$ where ${\mathrm{\nabla }}_{\text{β}}g=\mathrm{\partial }g(y,x,z,\text{β})/\mathrm{\partial }\text{β}{|}_{\text{β}={\text{β}}_0}$, $g=g(y,x,z,{\text{β}}_0)$, and ${\text{β}}_0$ is the true value of β. We denote $\begin{array}{rl}A& =\mathbb{E}\left(g{g}^{\mathrm{\top }}\right)=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{12}^{\mathrm{\top }}& A_{22}\end{array}\right),\\ A_{22.1}& =A_{22}-A_{12}^{\mathrm{\top }}{A}_{11}^{-1}{A}_{12}.\end{array}$ Equivalently, the asymptotic variance can be written as $\begin{array}{rl}& \left[\mathbb{E}\right({\mathrm{\nabla }}_{\text{β}}{g}_1^{\mathrm{\top }}\left)A_{11}^{-1}\mathbb{E}\right({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}{g}_1)+\mathbb{E}({\mathrm{\nabla }}_{\text{β}}{g}_1^{\mathrm{\top }})A_{11}^{-1}{A}_{12}\\ & \times A_{22.1}^{-1}{A}_{21}{A}_{11}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}{g}_1\right){]}^{-1},\end{array}$ or $(J+A_{12}{A}_{22.1}^{-1}{A}_{21}{)}^{-1},$ where $A_{11}=J$ is Fisher's information matrix.

In the above approach, the estimating equation $g_3=h(y,x)-\overline{h}$ does not involve the parameter β. However, there are ways to achieve higher efficiency. For example, we define $g_2(x,z,\text{β})=\psi (x,z,\text{β})$ with $\begin{array}{rl}\psi (x,z,\text{β})& =\mathbb{E}\left\{h\right(Y,X)|X=x,Z=z\}-\overline{h}\\ & =\int h(y,x)f(y|x,z,\text{β})\mathrm{d}y-\overline{h}.\end{array}$ Then, $\mathbb{E}\left\{g_2\right(x,z,\text{β}\left)\right\}=0.$ If we combine the empirical log-likelihood based on the estimating equation $g_2$ and the log-likelihood $\sum _{i=1}^n\log f(y_i|x_i,z_i,\text{β})$ as in the previous section (see Equation (12(12) $\begin{array}{rl}{\mathrm{\Sigma }}_{11}& =\{J+\mathbb{E}({\mathrm{\nabla }}_{\text{β}}{\varphi }^{\mathrm{\top }}\left)\right(\mathbb{E}\varphi {\varphi }^{\mathrm{\top }}{)}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}\varphi \right){\}}^{-1},\\ {\mathrm{\Sigma }}_{22}& =\left\{\mathbb{E}\right({\mathrm{\nabla }}_{\text{β}}{\varphi }^{\mathrm{\top }}\left)J^{-1}\mathbb{E}\right({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}\varphi )+\mathbb{E}(\varphi {\varphi }^{\mathrm{\top }}){\}}^{-1}.\end{array}$(12) )), then the asymptotic variance of the resulting MLE $\widehat{\text{β}}$ is given by $\{J+\mathbb{E}({\mathrm{\nabla }}_{\text{β}}{\psi }^{\mathrm{\top }}\left)\right(\mathbb{E}\psi {\psi }^{\mathrm{\top }}{)}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}\psi \right){\}}^{-1}.$ In general, this approach can achieve better efficiency.

5.2. A comparison

Given two pairs of estimation functions, $\{g_1,g_3\}$ and $\{g_1,g_2\}$, we may wonder combining which pair leads to a better estimator if we directly compare their asymptotic variance formulae. Alternatively, we may enquire whether we should combine all three constraints $g=(g_1,g_2,g_3)$ together. Write $g_{12}=g_{21}^{\mathrm{\top }}=(g_1,g_2)$, $a=\mathbb{E}\left\{h^{\mathrm{\top }}\right(y,x){\mathrm{\nabla }}_{\text{β}}\log f(y|x,z,\text{β}\left)\right\}$, and $\begin{array}{rl}\mathbb{E}\left(g{g}^{\mathrm{\top }}\right)& =\left(\begin{array}{ccc}J& 0& a\\ 0& \mathbb{E}\left(\psi {\psi }^{\mathrm{\top }}\right)& \mathbb{E}\left(\psi {\psi }^{\mathrm{\top }}\right)\\ a^{\mathrm{\top }}& \mathbb{E}\left(\psi {\psi }^{\mathrm{\top }}\right)& \mathbb{E}\left(h{h}^{\mathrm{\top }}\right)\end{array}\right)\\ & =\left(\begin{array}{cc}{B}_{11}& B_{12}\\ B_{12}^{\mathrm{\top }}& B_{22}\end{array}\right),\\ B_{11}& =\left(\begin{array}{cc}J& 0\\ 0& \mathbb{E}\left(\psi {\psi }^{\mathrm{\top }}\right)\end{array}\right),B_{12}=\left(\begin{array}{c}a\\ \mathbb{E}\left(\psi {\psi }^{\mathrm{\top }}\right)\end{array}\right).\end{array}$ Using results from Qin and Lawless (1994) and $\begin{array}{rl}{\left(\begin{array}{cc}{B}_{11}& B_{12}\\ B_{21}& B_{22}\end{array}\right)}^{-1}& =\left(\begin{array}{cc}I& -B_{11}^{-1}{B}_{12}\\ 0& I\end{array}\right)\\ & \times \left(\begin{array}{cc}{B}_{11}^{-1}& 0\\ 0& B_{22.1}^{-1}\end{array}\right)\left(\begin{array}{cc}I& 0\\ -B_{21}{B}_{11}^{-1}& I\end{array}\right)\end{array}$ with $B_{22.1}=B_{11}-B_{12}^{\mathrm{\top }}{B}_{11}^{-1}{B}_{12}$, we find that the asymptotic variance of $\widehat{\text{β}}$ obtained by combining the three estimating equations and $\sum _{i=1}^n\log f(y_i|x_i,z_i,\text{β})$ is $\begin{array}{rl}& [J+\mathbb{E}({\mathrm{\nabla }}_{\text{β}}{\psi }^{\mathrm{\top }}\left)\right\{\mathbb{E}\left(\psi {\psi }^{\mathrm{\top }}\right){\}}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}\psi \right)\\ & +\mathbb{E}\left({\mathrm{\nabla }}_{\text{β}}{g}_{21}\right)B_{11}^{-1}{B}_{12}{B}_{22.1}^{-1}{B}_{21}{B}_{11}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}{g}_{12}\right){]}^{-1}.\end{array}$ It can be shown that $\mathbb{E}\left({\mathrm{\nabla }}_{\text{β}}g\right)=(-J,\mathbb{E}({\mathrm{\nabla }}_{\text{β}}\psi ),0)$ and $\mathbb{E}\left({\mathrm{\nabla }}_{\text{β}}{g}_{12}\right)=(-J,a).$ Immediately, we have $\begin{array}{rl}\mathbb{E}\left({\mathrm{\nabla }}_{\text{β}}{g}_{12}\right)B_{11}^{-1}{B}_{12}& =(-J,a)\left(\begin{array}{cc}{J}^{-1}& 0\\ 0& \left\{\mathbb{E}\right(\psi {\psi }^{\mathrm{\top }}){\}}^{-1}\end{array}\right)\\ & \times \left(\begin{array}{c}a\\ \mathbb{E}\left(\psi {\psi }^{\mathrm{\top }}\right)\end{array}\right)=0,\end{array}$ which implies that the asymptotic variance in the case where $g_1$, $g_2$, and $g_3$ are combined is the same as that in the case where $g_1$ and $g_2$ only are combined. This indicates that taking $g_3$ into account leads to no efficiency gain in the estimation of β.

The method of combining $g_2$ and the parametric likelihood $\prod _{i=1}^{n}f(y_i|x_i,z_i,\text{β})$ is better than that of combining $g_1$, $g_3$, and the parametric likelihood. To see this, recall that the asymptotical variances for the MLEs of β with the two methods are $V_1=\{J+\mathbb{E}({\mathrm{\nabla }}_{\text{β}}{\psi }^{\mathrm{\top }}\left)\right(\mathbb{E}\psi {\psi }^{\mathrm{\top }}{)}^{-1}\mathbb{E}\left({\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}\psi \right){\}}^{-1}$ and $V_2=(J+A_{12}{A}_{22.1}^{-1}{A}_{21}{)}^{-1}.$ It suffices to show that $V_2-V_1\ge 0$, namely, $V_2-V_1$ is non-negative definite.

5.3. Proof of $V_2-V_1\ge 0$

For convenience, we assume that $\mathbb{E}\left(h\right)=0$. As $\mathbb{E}\left({\mathrm{\nabla }}_{\text{β}}{\psi }^{\mathrm{\top }}\right)=A_{12}$ and $\psi =\mathbb{E}(h|X,Z)$, it suffices to show that (13) $\begin{array}{rl}& A_{22.1}-\mathbb{E}\left(\psi {\psi }^{\mathrm{\top }}\right)=(A_{22}-A_{21}{A}_{11}^{-1}{A}_{12})\\ & -\mathbb{E}\left[\right\{\mathbb{E}(h|X,Z){\}}^{\otimes 2}]\ge 0.\end{array}$(13) Let ${\mathbb{E}}_{\ast }$ and ${\mathbb{V}\mathrm{a}\mathrm{r}}_{\ast }$ denote $\mathbb{E}(\cdot |X,Z)$ and $\mathbb{V}\mathrm{a}\mathrm{r}(\cdot |X,Z)$, respectively. As $\begin{array}{rl}\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)& =\mathbb{E}\left\{{\left(\begin{array}{c}{g}_1\\ h\end{array}\right)}^{\otimes 2}\right\}\\ & =\mathbb{E}\left\{{\mathbb{V}\mathrm{a}\mathrm{r}}_{\ast }\left(\begin{array}{c}{g}_1\\ h\end{array}\right)\right\}+\mathbb{V}\mathrm{a}\mathrm{r}\left\{{\mathbb{E}}_{\ast }\left(\begin{array}{c}{g}_1\\ h\end{array}\right)\right\}\end{array}$ and ${\mathbb{E}}_{\ast }\left(g_1\right)=0$, it follows that $\begin{array}{rl}\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)& \ge \mathbb{V}\mathrm{a}\mathrm{r}\left\{{\mathbb{E}}_{\ast }\left(\begin{array}{c}{g}_1\\ h\end{array}\right)\right\}\\ & =\mathbb{E}\left(\begin{array}{cc}0& 0\\ 0& {\mathbb{E}}_{\ast }\left(h\right){\mathbb{E}}_{\ast }\left(h^{\mathrm{\top }}\right)\end{array}\right).\end{array}$ Multiplying both sides by $(-A_{21}{A}_{11}^{-1},I)$ from the left and by $(-A_{21}{A}_{11}^{-1},I{)}^{\mathrm{\top }}$ from the right, we arrive at $A_{22}-A_{21}{A}_{11}^{-1}{A}_{12}\ge \mathbb{E}\left\{{\mathbb{E}}_{\ast }\right(h\left){\mathbb{E}}_{\ast }\right(h^{\mathrm{\top }}\left)\right\},$ that is, inequality (13(13) $\begin{array}{rl}& A_{22.1}-\mathbb{E}\left(\psi {\psi }^{\mathrm{\top }}\right)=(A_{22}-A_{21}{A}_{11}^{-1}{A}_{12})\\ & -\mathbb{E}\left[\right\{\mathbb{E}(h|X,Z){\}}^{\otimes 2}]\ge 0.\end{array}$(13) ) holds, which implies $V_2-V_1\ge 0$.

6. Calibration of information from previous studies

We consider calibration of information using parametric likelihood, EL (Owen, 1988), and GMM (Hansen, 1982). When only summary information from previous studies is available, these three well-known methods can be used to calibrate such summary information and to make inferences about the unknown parameters of interest. We may wonder whether doing so results in efficiency loss compared with inferences based on the pooled data if they were all available. Zeng and Lin (2015) found that parametric-likelihood-based meta-analysis of summarized information retained first-order asymptotic efficiency compared with analysis based on individual data. We show here that EL and GMM also possess this property. This is extremely important, as individual data may involve privacy issues, whereas summarized information does not.

6.1. Efficiency comparison

Suppose that $(Y_{ij},X_{ij})$ ($j=1,2,\dots ,n_i;i=1,2,\dots ,K$) are independent observations from the same population. We consider two scenarios according to the model's assumption about the population.

1. The conditional probability function (i.e., the probability density/mass function of a continuous/discrete random variable) of Y given X has a parametric form $f(y|x,\text{β})$.

2. The population satisfies $\mathbb{E}\left\{g\right(Y,X,\text{β}\left)\right\}=0.$

Here, β is a finite-dimensional unknown parameter, and ${\text{β}}_{\ast }$ is its true value. Assume that data are available batch by batch, and that $n_i/n={\rho }_i\in (0,1)$, where $n=\sum _{i=1}^K{n}_i$. For the i-th batch ($i=1,2,\dots ,K$) of data:

1. under assumption (I), the parametric log-likelihood function of β is ${\ell }_{i,\mathrm{P}\mathrm{L}}\left(\text{β}\right)=\sum _{j=1}^{n_i}\log \left\{f\right(Y_{ij}|X_{ij},\text{β}\left)\right\};$

2. under assumption (II), we define an empirical log-likelihood function $\begin{array}{rl}{\ell }_{i,\mathrm{E}\mathrm{L}}\left(\text{β}\right)& =sup\{\sum _{j=1}^{n_i}\log (n_i{p}_j):p_j\ge 0,\\ & \sum _{j=1}^{n_i}{p}_j=1,\\ & \sum _{j=1}^{n_i}{p}_{j}g(Y_{ij},X_{ij};\text{β})=0\}\\ & =-\sum _{j=1}^{n_i}\log \{1+{\lambda }_i^{\mathrm{\top }}g(Y_{ij},X_{ij};\text{β}\left)\right\}\\ & -n_i\log \left(n_i\right),\end{array}$ where ${\lambda }_i$ satisfies $\sum _{j=1}^{n_i}\frac{g(Y_{ij},X_{ij};\text{β})}{1+{\lambda }_i^{\mathrm{\top }}g(Y_{ij},X_{ij};\text{β})}=0$;

3. under assumption (II), we define the objective function of the GMM method (GMM log-likelihood for short) as $\begin{array}{rl}{\ell }_{i,\mathrm{G}\mathrm{M}\mathrm{M}}\left(\text{β}\right)& =-{\left\{\sum _{j=1}^{n_i}g\right(Y_{ij},X_{ij};\text{β}\left)\right\}}^{\mathrm{\top }}\\ & \times {\mathrm{\Omega }}^{-1}\left\{\sum _{j=1}^{n_i}g\right(Y_{ij},X_{ij};\text{β}\left)\right\},\end{array}$ where $\mathrm{\Omega }=\mathbb{V}\mathrm{a}\mathrm{r}\left\{g\right(Y,X,{\text{β}}_{\ast }\left)\right\}$ and ${\text{β}}_{\ast }$ is the true value of β. In practice, ${\text{β}}_{\ast }$ is generally replaced by a consistent estimator of β in the expression for Ω. Using the true value ${\text{β}}_{\ast }$ of β does not affect the theoretical analysis presented in this section.

Let ${\ell }_i\left(\text{β}\right)={\ell }_{i,\mathrm{P}\mathrm{L}}\left(\text{β}\right)$, ${\ell }_{i,\mathrm{E}\mathrm{L}}\left(\text{β}\right)$, or ${\ell }_{i,\mathrm{G}\mathrm{M}\mathrm{M}}\left(\text{β}\right)$. Under certain regularity conditions, it can be verified that for $\text{β}={\text{β}}_{\ast }+O_p\left(n^{-1/2}\right)$, (14) $\begin{array}{rl}{\ell }_i\left(\text{β}\right)& =U_i^{\mathrm{\top }}\sqrt{n_i}(\text{β}-{\text{β}}_{\ast })-\frac{n_i}{2}(\text{β}-{\text{β}}_{\ast }{)}^{\mathrm{\top }}\\ & \times V(\text{β}-{\text{β}}_{\ast })+o_p\left(1\right).\end{array}$(14) In Case (a), $\begin{array}{rl}{U}_i& =n_i^{-\frac{1}{2}}\sum _{j=1}^{n_i}{\mathrm{\nabla }}_{\text{β}}\log \left\{f\right(Y_{ij}|X_{ij},{\text{β}}_{\ast }\left)\right\},\\ V& =\mathbb{V}\mathrm{a}\mathrm{r}[{\mathrm{\nabla }}_{\text{β}}\log \{f(Y|X,{\text{β}}_{\ast })\left\}\right].\end{array}$ In Case (b) $U_i=n_i^{-\frac{1}{2}}\sum _{j=1}^{n_i}g(Y_{ij},X_{ij};{\text{β}}_{\ast }),V=A_{12}{A}_{22}^{-1}{A}_{21},$ where $\begin{array}{rl}A& =\left(\begin{array}{cc}0& \mathbb{E}\left\{{\mathrm{\nabla }}_{\text{β}}{g}^{\mathrm{\top }}\right(Y,X;{\text{β}}_{\ast }\left)\right\}\\ \mathbb{E}\left\{{\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}g\right(Y,X;{\text{β}}_{\ast }\left)\right\}& \mathbb{E}\left\{g\right(Y,X;{\text{β}}_{\ast }\left)g\right(Y,X;{\text{β}}_{\ast }\left)\right\}\end{array}\right)\\ & \equiv \left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right).\end{array}$ In Case (c), $\begin{array}{rl}{U}_i& =-\left\{\mathbb{E}{\mathrm{\nabla }}_{\theta }{g}^{\mathrm{\top }}\right(Y,X,{\text{β}}_{\ast }\left)\right\}{\mathrm{\Omega }}^{-1}{n}_i^{-\frac{1}{2}}\sum _{j=1}^{n_i}g(Y_{ij},X_{ij},{\text{β}}_{\ast }),\\ V& =\left\{\mathbb{E}{\mathrm{\nabla }}_{\theta }{g}^{\mathrm{\top }}\right(Y,X,{\text{β}}_{\ast }\left)\right\}{\mathrm{\Omega }}^{-1}\left\{\mathbb{E}{\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}g\right(Y,X,{\text{β}}_{\ast }\left)\right\}.\end{array}$ We denote the MLE of β based on the i-th batch of data by ${\widehat{\text{β}}}_i=\arg max{\ell }_r\left(\text{β}\right).$ The above approximation implies that $\sqrt{n_i}({\widehat{\text{β}}}_i-{\text{β}}_{\ast })=V^{-1}{U}_i+o_p\left(1\right)\stackrel{d}{⟶}N(0,V^{-1}).$ When the K-th batch of individual data are available, we no longer have access to the individual data of the previous K−1 batches but only have summarized information $({\widehat{\text{β}}}_i,{\hat{\mathrm{\Sigma }}}_i),i=1,2,\dots ,K-1$, where ${\widehat{\text{β}}}_i$ is the MLE based on the i-th batch of data and ${\hat{\mathrm{\Sigma }}}_i=V^{-1}/n_i+o\left(n^{-1}\right)$, and we can define an augmented log-likelihood ${\ell }_A\left(\text{β}\right)={\ell }_K\left(\text{β}\right)-\frac{1}{2}\sum _{i=1}^{K-1}({\widehat{\text{β}}}_i-\text{β}{)}^{\mathrm{\top }}{\hat{\mathrm{\Sigma }}}_i^{-1}({\widehat{\text{β}}}_i-\text{β})$ and the corresponding MLE ${\widehat{\text{β}}}_A=\arg max{\ell }_A\left(\text{β}\right).$ For $\text{β}={\text{β}}_{\ast }+O_p\left(n^{-1/2}\right)$, using the approximation in (14(14) $\begin{array}{rl}{\ell }_i\left(\text{β}\right)& =U_i^{\mathrm{\top }}\sqrt{n_i}(\text{β}-{\text{β}}_{\ast })-\frac{n_i}{2}(\text{β}-{\text{β}}_{\ast }{)}^{\mathrm{\top }}\\ & \times V(\text{β}-{\text{β}}_{\ast })+o_p\left(1\right).\end{array}$(14) ), we have $\begin{array}{rl}{\ell }_A\left(\text{β}\right)& =U_K^{\mathrm{\top }}\sqrt{n_K}(\text{β}-{\text{β}}_{\ast })-\frac{n_K}{2}(\text{β}-{\text{β}}_{\ast })V(\text{β}-{\text{β}}_{\ast })\\ & -\frac{1}{2}\sum _{i=1}^{K-1}{n}_i(\text{β}-{\text{β}}_{\ast }{)}^{\mathrm{\top }}V(\text{β}-{\text{β}}_{\ast })\\ & +\sum _{i=1}^{K-1}{n}_i({\widehat{\text{β}}}_i-{\text{β}}_{\ast }{)}^{\mathrm{\top }}V(\text{β}-{\text{β}}_{\ast })+C+o_p(1)\\ & =n^{-1/2}\sum _{i=1}^K\sqrt{n_i}{U}_i^{\mathrm{\top }}\cdot \sqrt{n}(\text{β}-{\text{β}}_{\ast })\\ & -\frac{n}{2}(\text{β}-{\text{β}}_{\ast }{)}^{\mathrm{\top }}V(\text{β}-{\text{β}}_{\ast })+C+o_p(1),\end{array}$ where the constant C differs in different equations.

For comparison, based on the pooled data, in Case (a) we define the parametric log-likelihood as ${\ell }_{\mathrm{P}\mathrm{L}}\left(\text{β}\right)=\sum _{i=1}^K\sum _{j=1}^{n_i}\log \left\{f\right(Y_{ij}|X_{ij},\text{β}\left)\right\};$ in Case (b) we define the empirical log-likelihood function as $\begin{array}{rl}{\ell }_{\mathrm{E}\mathrm{L}}\left(\text{β}\right)& =sup\{\sum _{i=1}^K\sum _{j=1}^{n_i}\log (n{p}_{ij}):p_{ij}\ge 0,\\ & \sum _{i=1}^K\sum _{j=1}^{n_i}{p}_{ij}=1,\\ & \sum _{i=1}^K\sum _{j=1}^{n_i}{p}_{ij}g(Y_{ij},X_{ij};\text{β})=0\}\\ & =-\sum _{i=1}^K\sum _{j=1}^{n_i}\log \{1+{\lambda }^{\mathrm{\top }}g(Y_{ij},X_{ij};\text{β}\left)\right\}\\ & -\sum _{i=1}^K{n}_i\log \left(n_i\right),\end{array}$ where λ satisfies $\sum _{i=1}^K\sum _{j=1}^{n_i}\frac{g(Y_{ij},X_{ij};\text{β})}{1+{\lambda }^{\mathrm{\top }}g(Y_{ij},X_{ij};\text{β})}=0$; and in Case (c) we define the GMM log-likelihood as $\begin{array}{rl}{\ell }_{\mathrm{G}\mathrm{M}\mathrm{M}}\left(\text{β}\right)& =-{\left\{\sum _{i=1}^K\sum _{j=1}^{n_i}g\right(Y_{ij},X_{ij};\text{β}\left)\right\}}^{\mathrm{\top }}\\ & \times {\mathrm{\Omega }}^{-1}\left\{\sum _{i=1}^K\sum _{j=1}^{n_i}g\right(Y_{ij},X_{ij};\text{β}\left)\right\}.\end{array}$ Let the log-likelihood based on the pooled data be ${\ell }_{\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}}\left(\text{β}\right)={\ell }_{\mathrm{P}\mathrm{L}}\left(\text{β}\right)$, ${\ell }_{\mathrm{E}\mathrm{L}}\left(\text{β}\right)$, and ${\ell }_{\mathrm{G}\mathrm{M}\mathrm{M}}\left(\text{β}\right)$ in Cases (a), (b), and (c), respectively. Then, it can be shown that $\begin{array}{rl}{\ell }_{\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{d}}\left(\text{β}\right)& =n^{-1/2}\sum _{j=1}^K\sqrt{n_j}{U}_j^{\mathrm{\top }}\cdot \sqrt{n}(\text{β}-{\text{β}}_{\ast })\\ & -\frac{n}{2}(\text{β}-{\text{β}}_{\ast }{)}^{\mathrm{\top }}V(\text{β}-{\text{β}}_{\ast })+C+o_p(1),\end{array}$ for some constant C. Let ${\widehat{\text{β}}}_{\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{d}}=\arg max{\ell }_{\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{d}}\left(\text{β}\right)$. By comparing ${\ell }_{\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{d}}\left(\text{β}\right)$ and ${\ell }_A\left(\text{β}\right)$, we obtain ${\ell }_{\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{d}}\left(\text{β}\right)={\ell }_A\left(\text{β}\right)+C+o_p\left(1\right)$ and $\begin{array}{rl}\sqrt{n}({\widehat{\text{β}}}_A-{\text{β}}_{\ast })& =\sqrt{n}({\widehat{\text{β}}}_{\mathrm{p}\mathrm{o}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{d}}-{\text{β}}_{\ast })+o_p\left(1\right)\\ & =V^{-1}\cdot n^{-1/2}\sum _{j=1}^K\sqrt{n_j}{U}_j^{\mathrm{\top }}+o_p\left(1\right)\\ & \stackrel{d}{⟶}N(0,V^{-1}).\end{array}$ This indicates that compared with the methods, including parametric likelihood, EL, and GMM, based on all individual data, the calibration method based on the last batch of individual data and all summary results of the previous batches has no efficiency loss.

6.2. When nuisance parameters are present

For batch i, assume that the data $(Y_{ij},X_{ij})$ ($j=1,2,\dots ,n_i$) satisfy either $\mathrm{p}\mathrm{r}(Y_{ij}=y|X_{ij}=x)=f(y|x,\text{β},{\gamma }_i)$ or $\mathbb{E}\left\{g\right(Y,X,\text{β},{\gamma }_i\left)\right\}=0,$ where β is common but ${\gamma }_i$ is a batch-specific parameter. We define ${\ell }_r(\text{β},{\gamma }_r)$ in the same way as ${\ell }_r\left(\text{β}\right)$. Let $({\widehat{\text{β}}}_i,{\hat{\gamma }}_i)$ be the MLE of $(\text{β},{\gamma }_i)$ based on the i-th batch of data, and assume that approximately $\left(\right({\widehat{\text{β}}}_i-\text{β}{)}^{\mathrm{\top }},({\hat{\gamma }}_i-{\gamma }_i{)}^{\mathrm{\top }}{)}^{\mathrm{\top }}\sim N(0,{\hat{\mathrm{\Sigma }}}_i)$ with ${\hat{\mathrm{\Sigma }}}_i=({\hat{\mathrm{\Sigma }}}_{i,rs}{)}_{1\le r,s\le 2}$.

We have two ways of combining information from previous studies. If we use all the previous summary information, we can define (15) $\begin{array}{rl}& {\ell }_A^{\left(1\right)}(\text{β},{\gamma }_1,\dots ,{\gamma }_K)={\ell }_K(\text{β},{\gamma }_K)-\frac{1}{2}\sum _{i=1}^{K-1}\left(\right({\widehat{\text{β}}}_i-\text{β}{)}^{\mathrm{\top }},\\ & ({\hat{\gamma }}_i-{\gamma }_i{)}^{\mathrm{\top }}){\hat{\mathrm{\Sigma }}}_i^{-1}\left(\right({\widehat{\text{β}}}_i-\text{β}{)}^{\mathrm{\top }},({\hat{\gamma }}_i-{\gamma }_i{)}^{\mathrm{\top }}{)}^{\mathrm{\top }}.\end{array}$(15) As ${\widehat{\text{β}}}_i|{\hat{\gamma }}_i\sim N(\text{β},{\hat{\mathrm{\Sigma }}}_{i,11\cdot 2}),$ where ${\hat{\mathrm{\Sigma }}}_{i,11\cdot 2}={\hat{\mathrm{\Sigma }}}_{i,11}-{\hat{\mathrm{\Sigma }}}_{i,12}{\hat{\mathrm{\Sigma }}}_{i,22}^{-1}{\hat{\mathrm{\Sigma }}}_{i,21},$ using only this summary information, we can define $\begin{array}{rl}{\ell }_A^{\left(2\right)}(\text{β},{\gamma }_K)& ={\ell }_K(\text{β},{\gamma }_K)-\frac{1}{2}\sum _{i=1}^{K-1}({\widehat{\text{β}}}_i-\text{β}{)}^{\mathrm{\top }}\\ & \times {\hat{\mathrm{\Sigma }}}_{i,11\cdot 2}^{-1}({\widehat{\text{β}}}_i-\text{β}).\end{array}$ Below we show that the MLEs of β based on these two likelihoods are actually equal to each other. In other words, there is no efficiency loss when estimating β based on ${\ell }_A^{\left(2\right)}(\text{β},{\gamma }_K)$ instead of ${\ell }_A^{\left(1\right)}(\text{β},{\gamma }_1,\dots ,{\gamma }_K)$.

To see this, it suffices to show that (16) $\underset{{\gamma }_1,\dots ,{\gamma }_{K-1}}{sup}{\ell }_A^{\left(1\right)}(\text{β},{\gamma }_1,\dots ,{\gamma }_K)={\ell }_A^{\left(2\right)}(\text{β},{\gamma }_K).$(16) We denote the inverse matrix of ${\mathrm{\Sigma }}_i$ by ${\mathrm{\Sigma }}_i^{-1}=({\mathrm{\Sigma }}_i^{rs}{)}_{1\le r,s\le 2},$ where $\begin{array}{rl}{\mathrm{\Sigma }}_i^{11}& ={\mathrm{\Sigma }}_{i,11\cdot 2}^{-1},{\mathrm{\Sigma }}_i^{21}=-{\mathrm{\Sigma }}_{i,22}^{-1}{\mathrm{\Sigma }}_{i,21}{\mathrm{\Sigma }}_{i,11\cdot 2}^{-1},\\ {\mathrm{\Sigma }}_i^{12}& =-{\mathrm{\Sigma }}_{i,11\cdot 2}^{-1}{\mathrm{\Sigma }}_{i,12}{\mathrm{\Sigma }}_{i,22}^{-1},\\ {\mathrm{\Sigma }}_i^{22}& ={\mathrm{\Sigma }}_{i,22}^{-1}+{\mathrm{\Sigma }}_{i,22}^{-1}{\mathrm{\Sigma }}_{i,21}{\mathrm{\Sigma }}_{i,11\cdot 2}^{-1}{\mathrm{\Sigma }}_{i,12}{\mathrm{\Sigma }}_{i,22}^{-1}.\end{array}$ It can be seen that $\begin{array}{rl}& {\ell }_A^{\left(1\right)}(\text{β},{\gamma }_1,\dots ,{\gamma }_K)\\ & ={\ell }_K(\text{β},{\gamma }_K)-\frac{1}{2}\sum _{i=1}^{K-1}({\widehat{\text{β}}}_i-\text{β}{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_i^{11}({\widehat{\text{β}}}_i-\text{β})\\ & +\sum _{i=1}^{K-1}({\widehat{\text{β}}}_i-\text{β}{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_i^{12}({\gamma }_i-{\hat{\gamma }}_i)\\ & -\frac{1}{2}\sum _{i=1}^{K-1}({\gamma }_i-{\hat{\gamma }}_i{)}^{\mathrm{\top }}{\mathrm{\Sigma }}_i^{22}({\gamma }_i-{\hat{\gamma }}_i).\end{array}$ Setting $\mathrm{\partial }{\ell }_A^{\left(1\right)}(\text{β},{\gamma }_1,\dots ,{\gamma }_K)/\mathrm{\partial }{\gamma }_i=0$ ($1\le i\le K-1$) gives $({\gamma }_i-{\hat{\gamma }}_i)=\left({\mathrm{\Sigma }}_i^{22}{)}^{-1}{\mathrm{\Sigma }}_i^{21}\right({\widehat{\text{β}}}_i-\text{β}).$ Putting this back into ${\ell }_A^{\left(1\right)}(\text{β},{\gamma }_1,\dots ,{\gamma }_K)$ gives $\begin{array}{rl}& \underset{{\gamma }_1,\dots ,{\gamma }_{K-1}}{sup}{\ell }_A^{\left(1\right)}(\text{β},{\gamma }_1,\dots ,{\gamma }_K)\\ & ={\ell }_K(\text{β},{\gamma }_K)-\frac{1}{2}\sum _{i=1}^{K-1}({\widehat{\text{β}}}_i-\text{β}{)}^{\mathrm{\top }}\\ & \times \{{\mathrm{\Sigma }}_i^{11}-{\mathrm{\Sigma }}_i^{12}({\mathrm{\Sigma }}_i^{22}{)}^{-1}{\mathrm{\Sigma }}_i^{21}\left\}\right({\widehat{\text{β}}}_i-\text{β})+C\\ & ={\ell }_K(\text{β},{\gamma }_K)-\frac{1}{2}\sum _{i=1}^{K-1}({\widehat{\text{β}}}_i-\text{β}{)}^{\mathrm{\top }}\\ & \times {\mathrm{\Sigma }}_{i,11\cdot 2}^{-1}({\widehat{\text{β}}}_i-\text{β})+C,\end{array}$ where we used the definition of ${\mathrm{\Sigma }}_{i,11\cdot 2}$ in the last equation. We arrive at Equation (16(16) $\underset{{\gamma }_1,\dots ,{\gamma }_{K-1}}{sup}{\ell }_A^{\left(1\right)}(\text{β},{\gamma }_1,\dots ,{\gamma }_K)={\ell }_A^{\left(2\right)}(\text{β},{\gamma }_K).$(16) ) after comparing this with the definition of ${\ell }_A^{\left(2\right)}(\text{β},{\gamma }_K)$.

7. Using covariate-specific disease prevalent information

As discussed in the previous section, summarized statistics from previous studies can sometimes be utilized to enhance the estimation efficiency in a current study. This is especially important in the big data era, when many types of information can be found through the internet. More specifically, suppose the prevalence of a disease is known at various levels of a known risk factor X. In this section, we combine this type of information in a case–control biased sampling setup.

7.1. Induced estimating equations under case–control sampling

Case–control sampling is among the most popular methods in cancer epidemiological studies. This is mainly because it is the most convenient, economic, and effective method. In the study of rare diseases in particular, one has to collect large samples in order to get a reasonable number of cases by using prospective sampling, which may not be practical. Using case–control sampling, a pre-specified number of cases ($n_1$) and controls ($n_0$) are collected retrospectively from case and control populations, respectively. Typically, this can be accomplished by sampling cases from hospitals and controls from the general disease-free population.

For a given risk factor X, let $F_i\left(x\right)=\mathrm{p}\mathrm{r}(X\le x|D=i)$ for i = 0, 1. Given X in a range $(a,b]$, the disease prevalence is $\mathrm{p}\mathrm{r}(D=1|a<X\le b)=\varphi (a,b),$ where $\varphi (a,b)$ is known. Using Bayes' formula, we have $\begin{array}{rl}\varphi (a,b)& =\frac{\pi {\int }_a^b\mathrm{d}{F}_1\left(x\right)}{\mathrm{p}\mathrm{r}(a<X\le b)},\\ 1-\varphi (a,b)& =\frac{(1-\pi ){\int }_a^b\mathrm{d}{F}_0\left(x\right)}{\mathrm{p}\mathrm{r}(a<X\le b)}\end{array}$ with $\pi =\mathrm{p}\mathrm{r}(D=1).$ It follows that ${\int }_a^b\mathrm{d}{F}_1\left(x\right)=\frac{1-\pi }{\pi }\frac{\varphi (a,b)}{1-\varphi (a,b)}{\int }_a^b\mathrm{d}{F}_0\left(x\right),$ or $\begin{array}{rl}& {\mathbb{E}}_1\left[I\right(a<X\le b\left)\right]\\ & =\frac{1-\pi }{\pi }\frac{\varphi (a,b)}{1-\varphi (a,b)}{\mathbb{E}}_0\left[I\right(a<X\le b\left)\right],\end{array}$ where ${\mathbb{E}}_0$ and ${\mathbb{E}}_1$ denote the expectation operators with respect to $F_0$ and $F_1$, respectively.

We assume that given covariates X and Y, the underlying disease model is given by the conventional logistic regression (17) $\begin{array}{r}\mathrm{p}\mathrm{r}(D=1|x,y)=\frac{\exp ({\alpha }^{\ast }+x\text{β}+y\gamma +yx\xi )}{1+\exp ({\alpha }^{\ast }+x\text{β}+y\gamma +yx\xi )}.\end{array}$(17) Let $\alpha ={\alpha }^{\ast }-\eta$ with $\eta =\log \left\{\right(1-\pi \left)/\pi \right\}$. It can be shown (see Qin, 2017) that this is equivalent to the exponential tilting model $\begin{array}{rl}{f}_1(x,y)& =f(x,y|D=1)\\ & =\exp (\alpha +x\text{β}+y\gamma +yx\xi )f_0(x,y),\end{array}$ where $f_0(x,y)=f(x,y|D=0)$. As a consequence, $\begin{array}{rl}& {\mathbb{E}}_0\left\{I\right(a<X\le b)\exp (\eta +\alpha +\text{β}X+\gamma Y+\xi XY\left)\right\}\\ & =\frac{1-\pi }{\pi }\frac{\varphi (a,b)}{1-\varphi (a,b)}{\mathbb{E}}_0\left[I\right(a<X\le b\left)\right]\end{array}$ or (18) $\begin{array}{rl}& {\mathbb{E}}_0\left\{I\right(a<X\le b)\exp (\alpha +\text{β}X+\gamma Y+\xi XY)\\ & -\frac{\varphi (a,b)}{1-\varphi (a,b)}I(a<X\le b)\}=0.\end{array}$(18) We denote $g_0(X,Y)=\exp (\eta +\alpha +\text{β}{X}_i+\gamma Y_i+\xi X_i{Y}_i)-1$ and the summarized auxiliary information equations as $\begin{array}{rl}{g}_i(X,Y)& =I(a_{i-1}<X\le a_i)\exp (\alpha +\text{β}X+\gamma Y\\ & +\xi XY)-\frac{\varphi (a_{i-1},a_i)}{1-\varphi (a_{i-1},a_i)}\\ & \times I(a_{i-1}<X\le a_i)\end{array}$ with $i=1,2,\dots ,I$. Then ${\mathbb{E}}_0\left\{g\right(X,Y\left)\right\}=0$, where $g(X,Y)=\left(g_0\right(X,Y),g_1(X,Y),\dots ,g_I(X,Y){)}^{\mathrm{\top }}.$

7.2. Empirical likelihood approach

The log-likelihood is (19) $\begin{array}{r}\ell =\sum _{i=1}^n{d}_i(\eta +\alpha +\text{β}{x}_i+\gamma y_i+\xi x_i{y}_i)+\sum _{i=1}^n\log \left(p_i\right),\end{array}$(19) where $p_i=d{F}_0\left(x_i\right),i=1,2,\dots ,n$, and the constraints are $p_i\ge 0,\sum _{i=1}^n{p}_i=1,\sum _{i=1}^n{p}_{i}g(x_i,y_i)=0.$ The profile log-likelihood is $\begin{array}{rl}\ell & =\sum _{i=1}^n{d}_i(\eta +\alpha +\text{β}{x}_i+\gamma y_i+\xi x_i{y}_i)\\ & -\sum _{i=1}^n\log \{1+{\lambda }^{\mathrm{\top }}g(x_i,y_i\left)\right\},\end{array}$ where the Lagrange multiplier λ is determined by $\sum _{i=1}^n\frac{g(x_i,y_i)}{1+{\lambda }^{\mathrm{\top }}g(x_i,y_i)}=0.$ Finally, the underlying parameters can be obtained by maximizing ℓ.

If the overall disease prevalence probability $\pi =\mathrm{p}\mathrm{r}(D=1)$ is known, then $\eta =\log \left\{\right(1-\pi \left)/\pi \right\}$ is known. On the other hand, if it is unknown but $I\ge 1$, then π is identifiable. If I>1, then we have an over-identified equation problem. This can be treated as a generalization of the EL method for estimating functions (Qin & Lawless, 1994) for biased sampling problems. Qin et al. (2015) considered the case where η is unknown and $I\ge 1$.

Let $\omega =(\eta ,\alpha ,\text{β},\gamma ,\xi ,\lambda )$ and let $\hat{\omega }$ be its maximum EL estimator. As the first estimating function $g_0$ corrects biased sampling in a case–control study, the remaining estimating functions $g_1,\dots ,g_I$ are used for improving efficiency. When n goes to infinity, it can be shown that the limit of λ is a $(I+1)$-dimensional vector where the first component is $\underset{n\to \mathrm{\infty }}{lim}\left(n_1/n\right)$ and the remainder are all zero. Qin et al. (2015) showed that if $\rho =n_1/n_0$ remains constant as $n\to \mathrm{\infty }$ and $\rho \in (0,1)$, then under suitable regularity conditions $\sqrt{n}(\hat{\omega }-{\omega }_0)$ is asymptotically normally distributed with mean zero. Moreover, the estimation of the logistic regression parameters $(\text{β},\gamma ,\xi )$ improves as the number I of estimating functions increases. This means that a richer set of auxiliary information leads to better estimators. In practice, however, this consideration must be balanced with the numerical difficulty of solving a larger number of equations.

Notably, auxiliary information is informative for estimating β and ξ but not for estimating γ. This can be observed through the following equations: $\begin{array}{rl}& \int I(a<x<b)\exp (\alpha +\text{β}x+\gamma y+\xi xy)\mathrm{d}{F}_0(x,y)\\ & =\int I(a<x<b)\exp \{\alpha +\text{β}x+s+\xi x(s/\gamma \left)\right\}\\ & \times \mathrm{d}{F}_0(x,s/\gamma ).\end{array}$ As the underlying distribution $F_0(x,y)$ is unspecified, we can treat $F_0(x,s/\gamma )$ as a new underlying distribution $F_0^{\ast }(x,s)$. With $F_0^{\ast }$ profiled out, the auxiliary information equation does not involve γ if $\xi =0$. Hence, even if $\xi \ne 0$, the information for γ is minimal as γ and ξ cannot be separated.

7.3. Generalizations

The simulation results of Qin et al. (2015) indicate that when covariate-specific auxiliary information is employed, the estimator of the coefficient β of X has the maximum variance reduction, whereas the variance reductions for other coefficients are small. If the auxiliary information $\mathrm{p}\mathrm{r}(D=1|b_{j-1}<Y\le b_j)={\psi }_j,j=1,2,\dots ,J$ is also available, we can combine them through estimating equations $\begin{array}{rl}{g}_i(X,Y)& =I(a_{i-1}<X\le a_i){\mathrm{e}}^{\alpha +\text{β}X+\gamma Y+\xi XY)}\\ & -\frac{\varphi (a_{i-1},a_i)}{1-\varphi (a_{i-1},a_i)}I(a_{i-1}<X\le a_i),\\ h_j(X,Y)& =I(b_{j-1}<Y\le b_j){\mathrm{e}}^{\alpha +\text{β}X+\gamma Y+\xi XY}\\ & -\frac{\psi (b_{j-1},b_j)}{1-\psi (b_{j-1},b_j)}I(b_{j-1}<Y\le b_j).\end{array}$ It would be more informative if the auxiliary information $\mathrm{p}\mathrm{r}(D=1|a<X<b,c<Y<d)$ is available.

7.4. More on the use of auxiliary information

Under a logistic regression model, the case and control densities are linked by the exponential tilting model (20) $\begin{array}{rl}& \mathrm{p}\mathrm{r}(x,y|D=1)=\mathrm{p}\mathrm{r}(x,y|D=0)\\ & \times \exp (\alpha +x\text{β}+y\gamma +\xi xy).\end{array}$(20) Suppose that for the general population $\mathbb{E}\left(X\right)={\mu }_1$, $\mathbb{E}\left(Y\right)={\mu }_2$, and $\mathbb{E}\left(XY\right)={\mu }_3$ are all known, and $\pi =\mathrm{p}\mathrm{r}(D=1)$ is known or can be estimated using external data. Under the exponential tilting model (20(20) $\begin{array}{rl}& \mathrm{p}\mathrm{r}(x,y|D=1)=\mathrm{p}\mathrm{r}(x,y|D=0)\\ & \times \exp (\alpha +x\text{β}+y\gamma +\xi xy).\end{array}$(20) ), the density $f(x,y)$ in the general population and the density $\mathrm{p}\mathrm{r}(x,y|D=0)$ in the control population are linked by $\begin{array}{rl}\mathrm{p}\mathrm{r}(x,y)& =\{\pi {\mathrm{e}}^{\alpha +x\text{β}+y\gamma +\xi xy}+(1-\pi \left)\right\}\\ & \times \mathrm{p}\mathrm{r}(x,y|D=0).\end{array}$ As a consequence $\mathbb{E}\left(X\right)={\mathbb{E}}_0\left[X\right\{\pi {\mathrm{e}}^{\alpha +X\text{β}+Y\gamma +\xi XY}+(1-\pi )\left\}\right]={\mu }_1,$ where ${\mathbb{E}}_0$ is an expectation with respect to $\mathrm{p}\mathrm{r}(x,y|D=0)$. Let $h(x,y)=(x-{\mu }_1,y-{\mu }_2,xy-{\mu }_3)$ with known ${\mu }_1$,${\mu }_2$, and ${\mu }_3$. The log-likelihood under case–control data is still (19(19) $\begin{array}{r}\ell =\sum _{i=1}^n{d}_i(\eta +\alpha +\text{β}{x}_i+\gamma y_i+\xi x_i{y}_i)+\sum _{i=1}^n\log \left(p_i\right),\end{array}$(19) ), where the $p_i$s satisfy the following constraints: $\begin{array}{rl}& \sum _{i=1}^n{p}_i=1,p_i\ge 0,\sum _{i=1}^n{p}_i{\mathrm{e}}^{\alpha +x_i\text{β}+y_i\gamma +x_i{y}_i\xi }=1,\\ & \sum _{i=1}^n{p}_{i}h(x_i,y_i)\{\pi {\mathrm{e}}^{\alpha +x_i\text{β}+y_i\gamma +x_i{y}_i\xi }+(1-\pi \left)\right\}=0.\end{array}$ More generally, any information in the general population such as $\mathbb{E}\left[\psi \right(Y,X\left)\right]=0$ can be converted to an equation for the control population, ${\mathbb{E}}_0\left[\right\{\pi {\mathrm{e}}^{\alpha +X\text{β}+Y\gamma +\xi XY}+(1-\pi )\left\}\psi \right(Y,X\left)\right]=0.$ Therefore, the results developed by Qin et al. (2015) can be applied. The results of Chatterjee et al. (2016) for case–control data can be considered as a special case of Qin et al. (2015).

8. Communication-efficient distributed inference

In the era of big data, it is commonplace for data analyses to run on hundreds or thousands of machines, with the data distributed across those machines and no longer available in a single central location. Recently, parallel and distributed inference has become popular in the statistical literature in both frequentist and Bayesian settings. In essence, the data-parallel procedures are intended to break the overall dataset into subsets that are processed independently. To the extent that communication-avoiding procedures have been discussed explicitly, the focus has been on one-shot or embarrassingly parallel approaches that use only one round of communication in which estimators or posterior samples are first obtained in parallel on local machines, then communicated to a centre node, and finally combined to form a global estimator or approximation to the posterior distribution (Lee et al., 2017; Neiswanger et al., 2015; Wang & Dunson, 2015; Zhang et al., 2013). In the frequentist setting, most one-shot approaches rely on averaging (Zhang et al., 2013), where the global estimator is the average of the local estimators. Lee et al. (2017) extend this idea to high-dimensional sparse linear regression by combining local debiased Lasso estimates (van de Geer et al., 2014). Recent work by Duchi et al. (2015) shows that under certain conditions, these averaging estimators can attain the information-theoretic complexity lower bound for linear regression, and at least $O\left(dk\right)$ bits must be communicated in order to attain the minimax rate of parameter estimation, where d is the dimension of the parameter and k is the number of machines. This result holds even in the sparse setting (Braverman et al., 2016).

The method of Jordan et al. (2019) proceeds as follows. Suppose the big data consists of N observations and there are k machines. For the convenience of presentation, we assume that each machine has n observations, i.e., N = nk. Denote the full-data likelihood by ${\ell }_N\left(\theta \right)=\frac{1}{k}\sum _{j=1}^k{\ell }_j\left(\theta \right),$ where ${\ell }_j\left(\theta \right)$ is the log-likelihood based on the data from the j-th machine. For θ near its target value $\overline{\theta }$, $\begin{array}{rl}{\ell }_N\left(\theta \right)& ={\ell }_N\left(\overline{\theta }\right)+{\mathrm{\nabla }}_{\theta }{\ell }_N\left(\theta \right){|}_{\theta =\overline{\theta }}(\theta -\overline{\theta })+R_N\left(\theta \right),\\ {\ell }_1\left(\theta \right)& ={\ell }_1\left(\overline{\theta }\right)+{\mathrm{\nabla }}_{\theta }{\ell }_1\left(\theta \right){|}_{\theta =\overline{\theta }}(\theta -\overline{\theta })+R_1\left(\theta \right),\end{array}$ where $R_N\left(\theta \right)$ and $R_1\left(\theta \right)$ are remainders. Observing that $R_N\approx R_1$, define a surrogate log-likelihood $\begin{array}{rl}\overline{\ell }\left(\theta \right)& ={\ell }_N\left(\overline{\theta }\right)+(\theta -\overline{\theta }{)}^{\mathrm{\top }}{\mathrm{\nabla }}_{\theta }{\ell }_N(\theta ){|}_{\theta =\overline{\theta }}\\ & +\left\{{\ell }_1\right(\theta )-{\ell }_1(\overline{\theta })-(\theta -\overline{\theta }{)}^{\mathrm{\top }}{\mathrm{\nabla }}_{\theta }{\ell }_1\left(\theta \right){|}_{\theta =\overline{\theta }}\}.\end{array}$ Ignoring the constant terms, the surrogate log-likelihood is $\overline{\ell }\left(\theta \right)={\ell }_1\left(\theta \right)+{\theta }^{\mathrm{\top }}\left\{{\mathrm{\nabla }}_{\theta }{\ell }_N\right(\theta ){|}_{\theta =\overline{\theta }}-{\mathrm{\nabla }}_{\theta }{\ell }_1(\theta \left){|}_{\theta =\overline{\theta }}\right\}.$ The score equation based on the surrogate likelihood is $\begin{array}{rl}{\mathrm{\nabla }}_{\theta }\overline{\ell }\left(\theta \right)& ={\mathrm{\nabla }}_{\theta }{\ell }_1\left(\theta \right)+\left\{{\mathrm{\nabla }}_{\theta }{\ell }_N\right(\theta ){|}_{\theta =\overline{\theta }}-{\mathrm{\nabla }}_{\theta }{\ell }_1(\theta \left){|}_{\theta =\overline{\theta }}\right\}\\ & =0.\end{array}$ Let $\hat{\theta }$ be the solution. Expanding it at ${\theta }_0$ and using the fact that $\begin{array}{rl}& N^{-1}\left\{{\mathrm{\nabla }}_{\theta {\theta }^{\mathrm{\top }}}{\ell }_1\right({\theta }_0)-{\mathrm{\nabla }}_{\theta {\theta }^{\mathrm{\top }}}{\ell }_N({\theta }_0\left)\right\}\\ & \to 0\text{in probability},\end{array}$ we can easily show that if $\overline{\theta }-{\theta }_0=O_p\left(N^{-1/2}\right)$ then $(\hat{\theta }-{\theta }_0)={\left\{{\mathrm{\nabla }}_{\theta {\theta }^{\mathrm{\top }}}{\ell }_N\right({\theta }_0\left)\right\}}^{-1}{\mathrm{\nabla }}_{\theta }{\ell }_N\left({\theta }_0\right)+o_p\left(N^{-1/2}\right).$ If we let $\overline{\theta }$ be the MLE based on ${\ell }_1\left(\theta \right)$, the surrogate log-likelihood can be simplified to $\overline{\ell }\left(\theta \right)={\ell }_1\left(\theta \right)+{\theta }^{\mathrm{\top }}{\mathrm{\nabla }}_{\theta }{\ell }_N\left(\overline{\theta }\right),$ because ${\mathrm{\nabla }}_{\theta }{\ell }_1\left(\overline{\theta }\right)=0$.

If the dimension of θ is high, one may add a penalty function in the surrogate log-likelihood and estimate θ by $\tilde{\theta }=\arg \underset{\theta \in \mathrm{\Theta }}{max}\left\{\overline{\ell }\right(\theta )-\lambda \Vert \theta {\Vert }_1\},$ where $\Vert \theta {\Vert }_1$ is the $L_1$-norm of θ. Similarly, Bayesian inference can be adapted to the surrogate likelihood as well.

Duan et al. (2020) proposed distributed algorithms that account for heterogeneous distributions by allowing site-specific nuisance parameters. The proposed methods extend the surrogate likelihood approach (Jordan et al., 2019; Wang et al., 2017) to the heterogeneous setting by applying a novel density ratio tilting method to the efficient score function. Asymptotically, the approach described in Section 6.2 on nuisance parameters is equivalent to that of Duan et al. (2020).

9. Renewal estimation and incremental inference

Let $U(D_1,\text{β})={\mathrm{\nabla }}_{\text{β}}M(D_1,\text{β})$ be a score function of β based on some objective function $M(D_1,\text{β})$ from the first batch of data, where M can be either the log-likelihood $M(D_1,\text{β})=\sum _{i=1}^{n_1}\log f(y_{1i}|x_{1i},\text{β})$ or a pseudo log-likelihood.

Let ${\widehat{\text{β}}}_1$ be the solution to $U(D_1,\text{β})=0,$ when only the first batch of data $D_1$ is available. Let $D_2$ denote the second batch of data. If both of them are available, we let ${\widehat{\text{β}}}_2$ be the solution to the pooled score equation, $U(D_1,\text{β})+U(D_2,\text{β})=0.$ Clearly, ${\widehat{\text{β}}}_2$ is the most efficient estimator of β when $D_1$ and $D_2$ are both available.

If $D_2$ is available but $D_1$ is not, with only some summary information ${\widehat{\text{β}}}_1$ and ${\hat{\mathrm{\Sigma }}}_1$ in its place, how can we utilize the summary information efficiently? It is not feasible to estimate β by directly solving $U\left(\text{β}\right)\equiv U(D_1,\text{β})+U(D_2,\text{β})=0,$ which involves the individual data of the unavailable $D_1$. Luo (2020) considers expanding $U(D_1,\text{β})$ at $\text{β}={\widehat{\text{β}}}_1$, i.e., $\begin{array}{rl}U(D_1,\text{β})& =U(D_1,{\widehat{\text{β}}}_1)+(\text{β}-{\widehat{\text{β}}}_1{)}^{\mathrm{\top }}{\mathrm{\nabla }}_{\text{β}}U(D_1,{\widehat{\text{β}}}_1)\\ & +O(\Vert \text{β}-{\widehat{\text{β}}}_1{\Vert }^2)\end{array}$ for β close to ${\widehat{\text{β}}}_1$. As $U(D_1,{\widehat{\text{β}}}_1)=0$, it follows that $\begin{array}{rl}U\left(\text{β}\right)& =U(D_2,\text{β})+(\text{β}-{\widehat{\text{β}}}_1{)}^{\mathrm{\top }}{\mathrm{\nabla }}_{\text{β}}U(D_1,{\widehat{\text{β}}}_1)\\ & +O(\Vert \text{β}-{\widehat{\text{β}}}_1{\Vert }^2).\end{array}$ Luo (2020) proposes obtaining an updated estimator of β by solving (21) $(\text{β}-{\widehat{\text{β}}}_1{)}^{\mathrm{\top }}{\mathrm{\nabla }}_{\text{β}}U(D_1,{\widehat{\text{β}}}_1)+U(D_2,\text{β})=0.$(21) Alternatively, we may understand this renewal estimation strategy in the manner of Zhang et al. (2020), who propose estimating β by maximizing (22) $\begin{array}{r}\sum _{i=1}^{n_2}\log f(y_{2i}|x_{2i},\text{β})-\frac{1}{2}{n}_1({\widehat{\text{β}}}_1-\text{β}{)}^{\mathrm{\top }}\mathrm{\Sigma }({\widehat{\text{β}}}_1-\text{β}),\end{array}$(22) where $\mathrm{\Sigma }=\mathbb{E}\{{\mathrm{\nabla }}_{\text{β}}\log f(Y|X,\text{β}){\mathrm{\nabla }}_{{\text{β}}^{\mathrm{\top }}}\log f(Y|X,\text{β}\left)\right\}$ is the Fisher information. If both batches are available, the score for β is $\begin{array}{rl}S\left(\text{β}\right)& =\sum _{i=1}^{n_1}{\mathrm{\nabla }}_{\text{β}}\log f(y_{1i}|x_{1i},\text{β})\\ & +\sum _{i=1}^{n_2}{\mathrm{\nabla }}_{\text{β}}\log f(y_{2i}|x_{2i},\text{β}).\end{array}$ After recording ${\widehat{\text{β}}}_1$ and Σ, we no longer have the raw data $\left\{\right(y_{1i},x_{1i}),i=1,2,\dots ,n_1\}$. As $\begin{array}{rl}{\widehat{\text{β}}}_1-\text{β}& =-n_1^{-1}{\mathrm{\Sigma }}^{-1}\sum _{i=1}^{n_1}{\mathrm{\nabla }}_{\text{β}}\log f(y_{1i}|x_{1i},\text{β})\\ & +o_p\left(n_1^{-1/2}\right),\end{array}$ differentiating (22(22) $\begin{array}{r}\sum _{i=1}^{n_2}\log f(y_{2i}|x_{2i},\text{β})-\frac{1}{2}{n}_1({\widehat{\text{β}}}_1-\text{β}{)}^{\mathrm{\top }}\mathrm{\Sigma }({\widehat{\text{β}}}_1-\text{β}),\end{array}$(22) ) with respect to β gives $\begin{array}{rl}& \sum _{i=1}^{n_2}{\mathrm{\nabla }}_{\text{β}}\log f(y_{2i}|x_{2i},\text{β})-n_1\mathrm{\Sigma }({\widehat{\text{β}}}_1-\text{β})\\ & =\sum _{i=1}^{n_1}{\mathrm{\nabla }}_{\text{β}}\log f(y_{1i}|x_{1i},\text{β})\\ & +\sum _{i=1}^{n_2}{\mathrm{\nabla }}_{\text{β}}\log f(y_{2i}|x_{2i},\text{β})+o_p\left(n^{1/2}\right).\end{array}$ Here, we have assumed that $n_1=O\left(n_2\right)=O\left(n\right)$. This indicates that estimating β by maximizing (22(22) $\begin{array}{r}\sum _{i=1}^{n_2}\log f(y_{2i}|x_{2i},\text{β})-\frac{1}{2}{n}_1({\widehat{\text{β}}}_1-\text{β}{)}^{\mathrm{\top }}\mathrm{\Sigma }({\widehat{\text{β}}}_1-\text{β}),\end{array}$(22) ) results in no efficiency loss asymptotically compared with the MLE based on all individual data, where the latter is infeasible in the current situation.

10. Concluding remarks

Rapid growth in hardware technology has made data collection much easier and more effective. In many applications, data often arrive in streams and chunks, which leads to batch-by-batch data or streaming data. For example, web sites served by widely distributed web servers may need to coordinate many distributed clickstream analyses, e.g., to track heavily accessed web pages as part of their real-time performance monitoring. Other examples include financial applications, network monitoring, security, telecommunications data management, manufacturing, and sensor networks (Babcock et al., 2002; Nguyen et al., 2021). The continuous arrival of such data in multiple, rapid, time-varying, possibly unpredictable and unbounded streams not only yields many fundamentally new research problems but provides contains various forms of auxiliary information.

Assembling information from different data sources has become indispensable in big data and artificial intelligence research. Statistical tools play an essential part in updating information. In this paper, we have presented a selective review of several traditional statistical methods, including meta-analysis, calibration information methods in survey sampling, and EL together with over-identified estimating equations and GMM. We have also briefly reviewed some recently developed statistical methods, including communication-efficient distributed statistical inference and renewal estimation and incremental inference, which can be regarded as the latest developments of calibration information methods in the era of big data. Although these methods were developed in different fields and in different statistical frameworks, in principle, they are asymptotically equivalent to well-known methods developed for meta-analysis. These methods result in almost no or little information loss compared with the case when full data are available.

Finally, we apologize to people whose work has inadvertently have been left out of our reference list.

Acknowledgments

The authors thank the editor and two referees for constructive comments and suggestions that led to significant improvements in this paper.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This research was supported by the National Natural Science Foundation of China [grant numbers 71931004, 12171157, and 32030063], the 111 Project [grant number B14019], the Development Fund for Shanghai Talents and the Natural Sciences and Engineering Research Council of Canada (grant number RGPIN-2020-04964).