# Adaptive and Efficient Learning with Blockwise Missing and Semi-Supervised
Data

Yiming Li  Department of Biostatistics, Columbia University, US  Xuehan Yang
Department of Biostatistics, Columbia University, US  Ying Wei  Department of
Biostatistics, Columbia University, US  Molei Liu111Corresponding email:
<EMAIL_ADDRESS>Department of Biostatistics, Columbia University, US

###### Abstract

Data fusion is an important way to realize powerful and generalizable analyses
across multiple sources. However, different capability of data collection
across the sources has become a prominent issue in practice. This could result
in the blockwise missingness (BM) of covariates troublesome for integration.
Meanwhile, the high cost of obtaining gold-standard labels can cause the
missingness of response on a large proportion of samples, known as the semi-
supervised (SS) problem. In this paper, we consider a challenging scenario
confronting both the BM and SS issues, and propose a novel Data-adaptive
projecting Estimation approach for data FUsion in the SEmi-supervised setting
(DEFUSE). Starting with a complete-data-only estimator, it involves two
successive projection steps to reduce its variance without incurring bias.
Compared to existing approaches, DEFUSE achieves a two-fold improvement.
First, it leverages the BM labeled sample more efficiently through a novel
data-adaptive projection approach robust to model misspecification on the
missing covariates, leading to better variance reduction. Second, our method
further incorporates the large unlabeled sample to enhance the estimation
efficiency through imputation and projection. Compared to the previous SS
setting with complete covariates, our work reveals a more essential role of
the unlabeled sample in the BM setting. These advantages are justified in
asymptotic and simulation studies. We also apply DEFUSE for the risk modeling
and inference of heart diseases with the MIMIC-III electronic medical record
(EMR) data.

Keywords: Missing Data; Data fusion; Control variate; Calibration; Linear
allocation; Statistical Efficiency.

##  1 Introduction

###  1.1 Background

Data fusion enables more powerful and comprehensive analysis by integrating
complementary views provided by different data sources, and has been more and
more frequently adopted with the increasing capacity in unifying and
synthesizing data models across different cohorts or institutions. As an
important example, in the field of health and medicine, there has been growing
interests and efforts in linking electronic medical records (EMRs) with
biobank data to tackle complex health issues (Bycroft et al.,, 2018; Castro et
al.,, 2022, e.g.). EMRs include detailed longitudinal clinical observations of
patients, offering a rich source of health information over time. By combining
these records with the genetics, genomics and other multi-omics data in
biobanks, this type of data fusion enables a much deeper understanding of
disease prognosis at the molecular level, potentially transforming patient
care and treatment strategies (Zengini et al.,, 2018; Verma et al.,, 2023;
Wang et al.,, 2024, e.g.). In addition to combine different types of data,
data fusion has also been scaled up to merge data from multiple institutions.
The UK Biobank (Sudlow et al.,, 2015) and the All-of-Us Research Program (The
All of Us Research Program Investigators,, 2019) in the U.S. are prime
examples.

To robustly and efficiently learn from multi-source fused data sets, several
main statistical challenges need to be addressed. One is blockwise missing
(BM) often occurring when certain variables are collected or defined
differently across the sources. Consequently, one or multiple entire blocks of
covariates are missing in the merged dataset. Dealing with the paucity of
accurate outcomes in EMR is another major problem. This is due to the
intensive human efforts or time costs demanded for collecting accurate
outcomes like the manual chart reviewing labels created by experts for certain
diseases or a long-term endpoint taking years to follow on. In this situation,
one is often interested in semi-supervised (SS) learning that combines some
small labeled subset with the large unlabeled sample collected from the same
population without observations of the true outcome. In the aforementioned EMR
applications, it is often to have the BM and SS problems occurring together,
e.g., our real example in Section 5. This new setup has a more complicated
data structure to be formally described in Section 1.2 and motivates our
methodological research in this paper.

###  1.2 Problem Setup

Let $Y$ be the outcome of interest and
$\boldsymbol{X}=(X_{1},X_{2},...,X_{p})^{\top}$ be a vector of $p$-dimensional
covariates. As illustrated on the left panel of Figure 1, there are three data
structures of the observations from the same population: (i)
$N_{\mathcal{L}\mathcal{C}}$ observations of $(Y,\boldsymbol{X})$ that are
labeled and complete, denoted as $\mathcal{L}\mathcal{C}$; (ii)
$N_{\mathcal{L}\mathcal{M}}$ labeled and covariates-missing observations of
$(Y,\boldsymbol{X}_{\mathcal{P}_{m^{c}}})$ denoted as $\mathcal{L}\mathcal{M}$
where $\mathcal{P}_{m}\subset\\{1,2,\ldots,p\\}$ is the index set of the
blockwise missing covariates and
$\mathcal{P}_{m^{c}}=\\{1,2,...,p\\}\setminus\mathcal{P}_{m}$; and (iii)
$N_{\mathcal{U}\mathcal{C}}$ unlabeled sample with complete observation of the
covariates denoted as $\mathcal{U}\mathcal{C}$. Note that we consider a
missing completely at random (MCAR) regime with the $\mathcal{L}\mathcal{C}$,
$\mathcal{L}\mathcal{M}$ and $\mathcal{U}\mathcal{C}$ samples having a
homogeneous underlying distribution. In the SS setting, we assume that
$N_{\mathcal{U}\mathcal{C}}\gg
N_{\mathcal{L}\mathcal{C}}+N_{\mathcal{L}\mathcal{M}}$. Denote by
$\rho={N_{\mathcal{L}\mathcal{M}}}/{N_{\mathcal{L}\mathcal{C}}}$. In our
application, $N_{\mathcal{L}\mathcal{C}}$ and $N_{\mathcal{L}\mathcal{M}}$
have similar magnitudes while our method and theory also accommodate the
$\rho\rightarrow 0$ and $\rho\rightarrow\infty$ scenarios. Let
$N=N_{\mathcal{L}\mathcal{C}}+N_{\mathcal{U}\mathcal{C}}+N_{\mathcal{L}\mathcal{M}}$,
and the observed samples as

$\boldsymbol{D}_{i}=\\{I(\mathcal{T}_{i}\neq\mathcal{U}\mathcal{C})Y_{i},I(\mathcal{T}_{i}\neq\mathcal{L}\mathcal{M})\boldsymbol{X}_{i,\mathcal{P}_{m}},\boldsymbol{X}_{i,\mathcal{P}_{m^{c}}}\\},\quad
i=1,2,\ldots,N.$

where
$\mathcal{T}_{i}\in\\{\mathcal{L}\mathcal{C},\mathcal{U}\mathcal{C},\mathcal{L}\mathcal{M}\\}$
indicates the data type of the observation $i$. Also, we are interested the
scenario with multiple BM structures:
$\mathcal{T}_{i}\in\\{\mathcal{L}\mathcal{C},\mathcal{U}\mathcal{C},\mathcal{L}\mathcal{M}_{1},\ldots,\mathcal{L}\mathcal{M}_{R}\\}$
with each $\mathcal{L}\mathcal{M}_{r}$ representing a subset of data with
observations of $(Y,\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}})$ and missing
covariates on the index set $\mathcal{P}_{m}^{r}$ satisfying
$\mathcal{P}_{m}^{r}\neq\mathcal{P}_{m}^{\ell}$ for any pair of $r\neq\ell$ in
$\\{1,\ldots,R\\}$; see the right panel of Figure 1.

Figure 1:  Data structures of the BM and SS setting considered in this paper.
The red block ($\mathcal{L}\mathcal{C}$) represents the labeled and complete
subset of data where both $Y$ and $\boldsymbol{X}$ are observed completely;
the yellow block ($\mathcal{L}\mathcal{M}$ or $\mathcal{L}\mathcal{M}_{r}$)
represent the labeled sample with $\boldsymbol{X}_{\mathcal{P}_{m}}$ or
$\boldsymbol{X}_{\mathcal{P}_{m}^{r}}$ missing; and the green block
($\mathcal{U}\mathcal{C}$) represents the subset of unlabeled data with a
complete observation of $\boldsymbol{X}$. The sample size of
$\mathcal{U}\mathcal{C}$ is usually much larger than the sample sizes of
$\mathcal{L}\mathcal{C}$ and $\mathcal{L}\mathcal{M}$. We consider both the
single BM scenario on the left penal and the multiple BM scenario on the
right.

Our primary interest lies in estimating and inferring the generalized linear
model (GLM):

$\mathrm{E}(Y|\boldsymbol{X})=g(\boldsymbol{\gamma}^{\top}\boldsymbol{X}),$
(1)

where $\boldsymbol{\gamma}=(\gamma_{1},\gamma_{2},...,\gamma_{p})^{\top}$ is a
vector of unknown coefficients and $g(\cdot)$ is a known link function. For
example, $g(a)=a$ corresponds to the Gaussian linear model and
$g(a)=e^{a}/(1+e^{a})$ for the logistic model. It is important to note that
our framework allows the model (1) to be misspecified, and we actually define
the model coefficients of our interests as $\bar{\boldsymbol{\gamma}}$, the
solution to the population-level moment equation:

$\mathbf{U}(\boldsymbol{\gamma})={\rm
E}\left[\boldsymbol{X}\left\\{Y-g\left(\boldsymbol{\gamma}^{\top}\boldsymbol{X}\right)\right\\}\right]=\mathbf{0}.$
(2)

The model parameter
$\boldsymbol{\gamma}=(\gamma_{1},\ldots,\gamma_{p})^{\top}$ can be identified
and estimated using the standard regression on the $\mathcal{L}\mathcal{C}$
sample but $\mathcal{L}\mathcal{M}$ and $\mathcal{U}\mathcal{C}$ cannot be
directly incorporated into this procedure, which results in loss of effective
samples. In this paper, our aim is to leverage both the
$\mathcal{L}\mathcal{M}$ and $\mathcal{U}\mathcal{C}$ samples to enhance the
estimation efficiency of $\boldsymbol{\gamma}$ without incurring bias.

###  1.3 Related Literature

With the increasing ability in collecting data with common models across
multiple institutions, integrative regression with blockwise missing
covariates has been an important problem frequently studied in recent
literature. One of the most well known and commonly used approach to address
missing data is Multiple Imputation with Chained Equations (MICE) White et
al., (2011). However, as a practical, general, and principal strategy, MICE is
sometimes subjected to high computation costs as well as potential bias and
inefficiency caused by the misspecification or excessive error in imputation,
e.g., see our simulation results.

Recently, Yu et al., (2020) proposed an approach named DISCOM for integrative
linear regression under BM. It approximates the covariance of
$(\boldsymbol{X},Y)$ through the optimal linear combination across different
data structures and incorporates it with an extended sparse regression to
derive an integrative estimator. Xue and Qu, (2021) proposed a multiple
blockwise imputation (MBI) approach that aggregates linear estimating
equations with multi-source imputed data sets constructed according to their
specific BM structures. They justified the efficiency gain of their proposal
over the complete-data-only estimator. Xue et al., (2021) extended this
framework work for high-dimensional inference with BM data sets. Inspired by
the idea of debiasing, Song et al., (2024) developed an approach for a similar
problem without imputing the missing blocks but showing better efficiency. Jin
and Rothenhäusler, (2023) proposed a modular regression approach accommodating
BM covariates under the conditional independence assumption across the blocks
of covariates given the response.

All works reviewed above primarily focused on the linear model, which is
subtly different from the GLM with a nonlinear link as the latter typically
requires the $\mathcal{L}\mathcal{C}$ sample to ensure identifiability of the
outcome model while the linear model does not. For integrative GLM regression
with BM data, Li et al., 2023b  proposed a penalized estimating equation
approach with multiple imputation (PEE–MI) that incorporates the correlation
of multiple imputed observations into the objective function and realizes
effective variable selection. Kundu and Chatterjee, (2023) and Zhao et al.,
(2023) proposed to incorporate the $\mathcal{L}\mathcal{M}$ sample through the
generalized methods of moments (GMM) that actually uses the score functions of
a reduced GLM (i.e., the parametric GLM for
$Y\sim\boldsymbol{X}_{\mathcal{P}_{m^{c}}}$) on the $\mathcal{L}\mathcal{M}$
sample for variance reduction. Such a “reduced model” strategy has been used
frequently in other work like Song et al., (2024), as well as for fusing
clinical and observational data to enhance the efficiency of causal inference
(Yang and Ding,, 2019; Hatt et al.,, 2022; Shi et al.,, 2023, e.g.).

As will be detailed in our asymptotic and numerical studies, the above-
introduced reduced model or imputation strategies does not characterize the
true influence of the $\mathcal{L}\mathcal{M}$ data and, thus, not achieving
semiparametric efficiency in the BM regression setting, which dates back to
the seminal work by Robins et al., (1994) and has been explored (Chen et al.,,
2008; Li and Luedtke,, 2023; Li et al., 2023a, ; Qiu et al.,, 2023, e.g.) in
various settings. In particular, Li and Luedtke, (2023) proposed derived the
canonical gradient for a general data fusion problem and justified its
efficient. Nevertheless, such semiparametric efficiency property is
established under an ideal scenario with a correct conditional model for
$\boldsymbol{X}_{\mathcal{P}_{m}}$, which is stringent in practice. Their
method could be less robust and effective in a more realistic case when the
imputation for $\boldsymbol{X}_{\mathcal{P}_{m}}$ is not consistent with the
truth. This issue has been recognized and addressed in recent literature
mainly focusing on the SS inference problem (Gronsbell et al.,, 2022; Schmutz
et al.,, 2022; Angelopoulos et al.,, 2023; Miao et al.,, 2023; Gan and Liang,,
2023; Wu and Liu,, 2023).

Our work is also closely related to recent progress in semi-supervised
learning (SSL) that aims at using the large unlabeled sample to enhance the
statistical efficiency of the analysis with relatively small labeled data
sets. On this track, Kawakita and Kanamori, (2013) and Song et al., (2023)
proposed to use some density ratio model matching the labeled and unlabeled
samples to reweight the regression on the labeled sample, which, counter-
intuitively, results in SSL estimators with smaller variances than their
supervised learning (SL) counterparts. Chakrabortty and Cai, (2018) and
Gronsbell et al., (2022) proposed to incorporate the unlabeled sample by
imputing its unobserved outcome $Y$ with an imputation model constructed with
the labeled data. Extending this imputation framework, Cai and Guo, (2020),
Zhang and Bradic, (2022) and Deng et al., (2023) studied the SSL problem with
high-dimensional $\boldsymbol{X}$, and Wu and Liu, (2023) addressed the SSL of
graphic models.

Though having different constructions, these recent SSL approaches achieve a
very similar adaptive property that they can produce SSL estimators more
efficient than SL when the outcome model for $\mathrm{E}(Y|\boldsymbol{X})$ is
misspecified and as efficient as SL when it is correct. Azriel et al., (2022)
discussed this property and realized it for linear models in a more flexible
and effective way. Notably, our SSL setting with BM designs has not been
studied in this track of SSL literature. As will be shown later, this problem
is essentially distinguished from existing work in SSL as the efficiency gain
by using the $\mathcal{U}\mathcal{C}$ sample is attainable even when the
outcome model (1) is correctly specified in our case. In addition, we notice
that Xue et al., (2021) and Song et al., (2024) also considered using large
$\mathcal{U}\mathcal{C}$ samples to assist learning of high-dimensional linear
models under the BM setting. Unlike us, they only used the
$\mathcal{U}\mathcal{C}$ sample to combat and reduce the excessive
regularization errors in sparse regression and debiased inference.
Nevertheless, for the inference of a low-dimensional projection of the model
parameters, the asymptotic variance still remains unchanged compared to the SL
scenario after incorporating the $\mathcal{U}\mathcal{C}$ sample in their
methods.

###  1.4 Our contribution

Regarding the limitation of existing work discussed in the previous section,
we propose a Data-adaptive projecting Estimation approach for data FUsion in
the SEmi-supervised setting (DEFUSE). Our method contains two major novel
steps, it initiates from the simple $\mathcal{L}\mathcal{C}$-only estimator,
and first incorporates the $\mathcal{L}\mathcal{M}$ sample to reduce its
asymptotic variance through a data-adaptive projection procedure for the score
functions. Then, it utilizes the large $\mathcal{U}\mathcal{C}$ sample in
another projecting step, to further improve the efficiency of the data-fused
estimator resulting from the previous step. In the following, we shall
summarize the main novelty and contribution of our work that have been
justified and demonstrated through the asymptotic, simulation, and real-world
studies.

  1. 1.

Our approach in combining the $\mathcal{L}\mathcal{M}$ and
$\mathcal{L}\mathcal{C}$ samples is shown to produce an asymptotically
unbiased and more efficient estimator compared to recent developments for BM
regression (Xue and Qu,, 2021; Zhao et al.,, 2023; Song et al.,, 2024, e.g.).
The newly proposed framework is inspired by the semiparametric theory and
allows the flexible use of general ML algorithms to impute the missing
variables. The efficiency gain of our method over existing methods in various
settings like model misspecification are justified from both the theoretical
and numerical perspectives.

  2. 2.

We develop a novel data-driven control function approach to realize more
adaptive and efficient estimation under potential misspecification or poor
quality of the imputation models, especially when there exists nonlinear
relationship not captured by them. Compared to recent work addressing
problematic nuisance models in a more simple SS scenario (Angelopoulos et
al.,, 2023; Miao et al.,, 2023, e.g.), our strategy can be generalized to
their settings and realize more flexible and efficient calibration.

  3. 3.

We propose a linear allocation approach to extend our framework in addressing
multiple BM structures, in which the classic semiparametric theory does not
apply.

  4. 4.

We establish a novel strategy to bring in large unlabeled samples under the SS
setting and leverage them to enhance the efficiency of the BM estimation. This
enables a more efficient use of the $\mathcal{U}\mathcal{C}$ sample leading to
an estimator of strictly smaller variance than that obtained only using
$\mathcal{L}\mathcal{C}$ and $\mathcal{L}\mathcal{M}$, which has not been
readily achieved in previous work on a similar setting (Xue et al.,, 2021;
Song et al.,, 2024). Our development is also shown to be more efficient than
existing approaches for the standard SS setting without BM (Gronsbell et al.,,
2022, e.g.), in both theoretical and numerical studies.

##  2 Method

###  2.1 Outline and Preliminary Framework

Without the loss of generality, we aim at estimating and inferring a linear
projection of the model coefficients $\boldsymbol{\gamma}$ denoted by
$\boldsymbol{c}^{\top}\boldsymbol{\gamma}$ where $\boldsymbol{c}$ is an
arbitrary vector satisfying $\|\boldsymbol{c}\|_{2}=1$. In this framework, one
can take $\boldsymbol{c}$ as the $j$-th unit vector in $\mathbb{R}^{p}$ to
infer each GLM effect $\gamma_{j}$ and can also take $\boldsymbol{c}$ as the
observed covariates of some new coming subject for the purpose of prediction.
Let $\widehat{{\rm E}}_{\mathcal{T}}[f(\boldsymbol{D})]$ represent the
empirical mean of some measurable function $f(\cdot)$ on the data set
$\mathcal{T}$ that can be taken as one or the union of more than one belonging
to
$\\{\mathcal{L}\mathcal{C},\mathcal{L}\mathcal{M},\mathcal{U}\mathcal{C}\\}$,
e.g.,

$\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[f]=\frac{1}{N_{\mathcal{L}\mathcal{C}}}\sum_{\mathcal{T}_{i}=\mathcal{L}\mathcal{C}}f(\boldsymbol{D_{i}}),\quad\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}\cup\mathcal{L}\mathcal{M}}[f]=\frac{1}{N_{\mathcal{L}\mathcal{C}}+N_{\mathcal{L}\mathcal{M}}}\sum_{\mathcal{T}_{i}\in\\{\mathcal{L}\mathcal{C},\mathcal{L}\mathcal{M}\\}}f(\boldsymbol{D_{i}}).$

Similarly, we use $\widehat{{\rm Var}}_{\mathcal{T}}[f]$ to represent the
empirical variance of $f(\boldsymbol{D})$ on data set $\mathcal{T}$, e.g.,
$\widehat{{\rm
Var}}_{\mathcal{L}\mathcal{C}}[f]=N_{\mathcal{L}\mathcal{C}}^{-1}\sum_{\mathcal{T}_{i}=\mathcal{L}\mathcal{C}}\\{f(\boldsymbol{D_{i}})-\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[f]\\}^{2}$. We start from introducing the DEFUSE
method for the scenario with a single $\mathcal{L}\mathcal{M}$ data set; see
the left panel of Figure 1. It includes three steps outlined in Figure 2 and
also in Algorithm 1 presented later.

Figure 2:  Outline of the proposed DEFUSE Method. See also Algorithm 1.

First, we solve the empirical version of the estimating equation (2) for
$\boldsymbol{\gamma}$ based on the $\mathcal{L}\mathcal{C}$ sample:

$\tilde{\mathbf{U}}_{\mathcal{L}\mathcal{C}}(\boldsymbol{\gamma})=\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[\boldsymbol{X}\\{Y-g(\boldsymbol{\gamma}^{\top}\boldsymbol{X})\\}]=\mathbf{0}.$
(3)

Denote the solution as $\tilde{\boldsymbol{\gamma}}$. It is also the maximum
likelihood estimator (MLE) for the Gaussian linear and logistic models when
$g(a)=a$ and $g(a)=e^{a}/(1+e^{a})$ respectively. Using the standard
asymptotic theory for M-estimation (Van der Vaart,, 2000), we have the
expansion

$\tilde{\boldsymbol{{\gamma}}}-\boldsymbol{\bar{\gamma}}\approx\frac{1}{N_{\mathcal{L}\mathcal{C}}}\sum_{\mathcal{T}_{i}={\mathcal{L}\mathcal{C}}}\bar{\boldsymbol{H}}^{-1}\boldsymbol{X}_{i}\\{Y_{i}-g(\boldsymbol{\bar{\gamma}}^{\top}\boldsymbol{X}_{i})\\}:=\frac{1}{N_{\mathcal{L}\mathcal{C}}}\sum_{\mathcal{T}_{i}={\mathcal{L}\mathcal{C}}}\bar{\boldsymbol{S}}_{i};\quad\bar{\boldsymbol{H}}=\mathbb{E}[\dot{g}(\boldsymbol{\bar{\gamma}}^{\top}\boldsymbol{X})\boldsymbol{X}\boldsymbol{X}^{\top}]$

which will be more strictly given by Lemma 1. Consequently, the asymptotic
variance of $\boldsymbol{c}^{\top}\tilde{\boldsymbol{{\gamma}}}$ is
proportional to that of the score function
$\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}_{i}$. Therefore, our basic idea of
incorporating the $\mathcal{L}\mathcal{M}$ sample is to leverage its
observations of $(Y_{i},\boldsymbol{X}_{i,\mathcal{P}_{m^{c}}})$ to form some
function explaining the variance of
$\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}_{i}$ as much as possible, and use
it as a control variate to adjust
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{{\gamma}}}$ for variance reduction.

In specific, we let $\tilde{\boldsymbol{H}}=\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[g(\tilde{\boldsymbol{{\gamma}}}^{\top}\boldsymbol{X})\boldsymbol{X}\boldsymbol{X}^{\top}]$
and
$\tilde{\boldsymbol{S}}_{i}=\tilde{\boldsymbol{H}}^{-1}\boldsymbol{X}_{i}\\{Y_{i}-g(\tilde{\boldsymbol{{\gamma}}}^{\top}\boldsymbol{X}_{i})\\}$.
Motivated by semiparametric theory (Robins et al.,, 1994, e.g.), the ideal
control function in this case is
$\varphi^{*}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)=\frac{\rho}{1+\rho}\mathbb{E}_{\boldsymbol{X}_{\mathcal{P}_{m}}}[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y]$
(with empirically form
$\frac{\rho}{1+\rho}\mathbb{E}_{\boldsymbol{X}_{\mathcal{P}_{m}}}[\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y]$),
with the corresponding estimator incorporating $\mathcal{L}\mathcal{M}$ formed
as

$\boldsymbol{c}^{\top}\tilde{\boldsymbol{{\gamma}}}-\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[\varphi^{*}]+\widehat{{\rm
E}}_{\mathcal{L}\mathcal{M}}[\varphi^{*}].$ (4)

Here, the augmentation term $-\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[\varphi^{*}]+\widehat{{\rm
E}}_{\mathcal{L}\mathcal{M}}[\varphi^{*}]$ can be viewed as a control variate
of mean zero and negatively correlated with the empirical mean of the score
$\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}$ on $\mathcal{L}\mathcal{C}$. Among
all the control variate in this form,
$\varphi^{*}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)$ provides the optimal
choice to achieve the semiparametric efficiency. However, a critical concern
in this classic framework is that computing
$\varphi^{*}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)$ requires the knowledge
of the conditional distribution ${\rm
P}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)$.
It is typically unknown in practice and needs to be replaced with some data
generation model
${\mu}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})$
where $\mu$ is pre-specified and $\tilde{\boldsymbol{\eta}}$ is an empirical
estimate of the unknown model parameters. Correspondingly, our imputation of
$\varphi^{*}$ is formed as

$\tilde{\varphi}_{1}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})=\frac{\rho}{1+\rho}\int\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\\{\boldsymbol{x}_{\mathcal{P}_{m}},\boldsymbol{X}_{\mathcal{P}_{m^{c}}}\\},Y){\mu}(\boldsymbol{x}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})\,d\boldsymbol{x}_{\mathcal{P}_{m}},$
(5)

and to be used as a control function for the variance reduction of
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$. In the following, we
remark and discuss on the constructing strategy and computation of
$\tilde{{\mu}}_{\boldsymbol{X}_{\mathcal{P}_{m}}}:={\mu}(\boldsymbol{x}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})$
and
$\tilde{\varphi}_{1,i}:=\tilde{\varphi}_{1}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})$
in practice.

######  Remark 1.

The data imputation model $\tilde{{\mu}}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$
can be either a simple parametric model with finite-dimensional
$\tilde{\boldsymbol{\eta}}$ or more sophisticated ones with diverging
parameter space. For example, for continuous
$\boldsymbol{X}_{\mathcal{P}_{m}}$, one can either (i) assume that
$\boldsymbol{X}_{\mathcal{P}_{m}}$ follows a (multivariate) Gaussian linear
model against $\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y$ and learn it through
the a standard MLE; or (ii) assume
$\boldsymbol{X}_{\mathcal{P}_{m}}=m(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)+\boldsymbol{\varepsilon}$
with $\boldsymbol{\varepsilon}\sim{\rm N}(0,\boldsymbol{\Sigma})$ independent
from $\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y$, then use some ML algorithm
(e.g., random forest or neural network) to learn the mean model
$m(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)$ and the standard method of moment
on the fitted residuals to estimate $\boldsymbol{\Sigma}$; or (iii) implement
generative ML like conditional generative adversarial network (GAN) (Mirza and
Osindero,, 2014) to learn ${\mu}$ without any model assumptions. Regarding
these choices, our asymptotic analysis generally accommodate all the three
types, and we use both (i) and (ii) in our simulation and real-world studies.

When using complex ML algorithms like random forests to learn
$\tilde{{\mu}}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$, we adopt cross-fitting to
avoid the over-fitting bias as frequently done in related contexts
(Chernozhukov et al.,, 2016; Liu et al.,, 2023, e.g.). In specific, we
randomly split the $\mathcal{L}\mathcal{C}$ sample into $K$ folds with an
equal size, denoted as
$\\{\mathcal{L}\mathcal{C}^{(1)},\mathcal{L}\mathcal{C}^{(2)},...,\mathcal{L}\mathcal{C}^{(K)}\\}$.
For each $k\in[K]$, we implement the ML method on
$\mathcal{L}\mathcal{C}^{(-k)}=\cup_{j\neq k}\mathcal{L}\mathcal{C}^{(j)}$ to
obtain an estimator $\tilde{\boldsymbol{\eta}}^{(-k)}$. Then for each $i$ from
$\mathcal{L}\mathcal{C}$, we use
$\tilde{\varphi}_{1}(\boldsymbol{X}_{i,\mathcal{P}_{m^{c}}},Y_{i};\tilde{\boldsymbol{\eta}}^{(-k)})$
as its score imputation since $\tilde{\boldsymbol{\eta}}^{(-k)}$ is
independent with $\boldsymbol{D}_{i}$. For each $i$ from
$\mathcal{L}\mathcal{M}$, we use the average imputation
$K^{-1}\sum_{k=1}^{K}\tilde{\varphi}_{1}(\boldsymbol{X}_{i,\mathcal{P}_{m^{c}}},Y_{i};\tilde{\boldsymbol{\eta}}^{(-k)})$.
One can find more details about this cross-fitting procedure in Appendix. As
shown in Section 3, this strategy is important to eliminate the over-fitting
bias of arbitrary ML methods. For the ease of presentation, we moderately
abuse the notation in the remaining paper by dropping the superscript of $k$
related to cross-fitting.

With $\tilde{{\mu}}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$, one can use either
numerical methods like Gaussian quadrature or Monte-Carlo (MC) sampling to
compute $\tilde{\varphi}_{1,i}$. When ${\boldsymbol{X}_{\mathcal{P}_{m}}}$ is
multiple or high dimensional, the computational burden and accuracy through
either way seems to be a challenging issue. However, inspecting the form of
the score function
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\boldsymbol{X},Y)$, we notice
that it depends on the entire $\boldsymbol{X}_{\mathcal{P}_{m}}$ only through
at most two-dimensional linear representations
$\tilde{\boldsymbol{\gamma}}_{\mathcal{P}_{m}}^{\top}\boldsymbol{X}_{\mathcal{P}_{m}}$
and
$(\boldsymbol{c}^{\top}\bar{\boldsymbol{H}}^{-1})_{\mathcal{P}_{m}}\boldsymbol{X}_{\mathcal{P}_{m}}$.
Thus, both Gaussian quadrature and MC can achieve good computational
performance thanks to this low dimensionality (Caflisch,, 1998; Kovvali,,
2022). For MC, we could also use the same set of MC samples on all subjects in
certain scenarios (e.g., Gaussian noise under homoscedasticity) to save
computational costs. In our case, we choose the MC method to compute
$\tilde{\varphi}_{1}$ as it facilitates our data-adaption procedure to be
introduced next.

###  2.2 Data-adaptive Calibration

The numerical approximation of (5) offers a convenient approach for estimating
the conditional expectation function $\varphi^{*}$. However, the nuisance
model $\tilde{{\mu}}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$, whether parametric
or nonparametric, can deviate from the truth ${\rm
P}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)$
in practice due to model misspecification or excessive ML errors. This
deviation typically does not cause excessive bias (as shown in our asymptotic
analysis) but could result in inefficient estimation in (4). Similar
observations have been made in recent missing data literature (Angelopoulos et
al.,, 2023; Wu and Liu,, 2023, e.g.). To address this issue, we propose a
novel data-adaptive calibration approach to further reduce the variance in the
presence of model misspecifications. The key idea is to introduce a broader
class of conditional model for $\boldsymbol{X}_{\mathcal{P}_{m^{c}}}$ with
additional unknown parameters $\boldsymbol{\delta}$, denoted as
$\mu^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\boldsymbol{\delta},\tilde{\boldsymbol{\eta}})$.
We can then calibrate through its free parameter $\boldsymbol{\delta}$ to
produce a more efficient estimator than (4) with wrong
$\tilde{{\mu}}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$.

In specific, let
$\\{\boldsymbol{X}_{i,\mathcal{P}_{m}}^{[q]}:q=1,2,\ldots,Q\\}$ be $Q$
independent random draw from the preliminary
$\mu(\boldsymbol{X}_{\mathcal{P}_{m}}\mid\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})$
for the $i$th subject in $\mathcal{L}\mathcal{C}\cup\mathcal{L}\mathcal{M}$.
The standard MC approximation of $\tilde{\varphi}_{1,i}$ defined in (5) is
given by
$Q^{-1}\sum_{q=1}^{Q}\frac{\rho}{1+\rho}\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\\{\boldsymbol{X}_{i,\mathcal{P}_{m^{c}}},\boldsymbol{X}_{i,\mathcal{P}_{m}}^{[q]}\\},Y_{i}).$
In Section 3, we show that a size $Q$ moderately increasing with the sample
size is sufficient to guarantee the asymptotic normality and efficiency of our
estimator. Inspired by importance weighting, we re-construct the control
function with our newly introduced
$\mu^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\boldsymbol{\delta},\tilde{\boldsymbol{\eta}})$
as follows

$\frac{\rho}{1+\rho}\int\frac{\mu^{\prime}(\boldsymbol{x}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\boldsymbol{\delta},\tilde{\boldsymbol{\eta}})}{\mu(\boldsymbol{x}_{\mathcal{P}_{m}}\mid\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})}\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\\{\boldsymbol{x}_{\mathcal{P}_{m}},\boldsymbol{X}_{\mathcal{P}_{m^{c}}}\\},Y)\left[\mu(\boldsymbol{x}_{\mathcal{P}_{m}}\mid\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})\,d\boldsymbol{x}_{\mathcal{P}_{m}}\right].$

Its MC approximation is given by

$\tilde{\varphi}_{1}^{\prime}(\boldsymbol{X}_{i,\mathcal{P}_{m^{c}}},Y_{i};\boldsymbol{\delta},\tilde{\boldsymbol{\eta}})=\frac{1}{Q}\sum_{q=1}^{Q}\frac{\rho}{1+\rho}w(\boldsymbol{X}_{i}^{[q]},Y_{i};\boldsymbol{\delta})\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\boldsymbol{X}_{i}^{[q]},Y_{i}),$

where
$\boldsymbol{X}_{i}^{[q]}=\\{\boldsymbol{X}_{i,\mathcal{P}_{m^{c}}},\boldsymbol{X}_{i,\mathcal{P}_{m}}^{[q]}\\}$
and $w(\cdot)$ is known as the importance weighting (IW) function such that

$w(\boldsymbol{X}_{i}^{[q]},Y;\boldsymbol{\delta})=\frac{\mu^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m}}^{[q]}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\boldsymbol{\delta},\tilde{\boldsymbol{\eta}})}{\mu(\boldsymbol{X}_{\mathcal{P}_{m}}^{[q]}\mid\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})}.$

Following the definition of $w$, we can construct $\mu^{\prime}$ as the
product of $\mu$ and an IW function $w$ that is only parameterized by
$\boldsymbol{\delta}$. It is easy to see that when $w\equiv 1$,
$\mu^{\prime}\equiv\mu$. Hence, $\mu^{\prime}$ encompasses a broader class of
conditional models of $\boldsymbol{X}_{P_{m}}$. We then replace $\varphi^{*}$
in (5) by $\tilde{\varphi}_{1}^{\prime}$, and calibrate it by optimizing
$\boldsymbol{\delta}$ as follows.

$\tilde{\boldsymbol{\delta}}=\arg\min_{\boldsymbol{\delta}}\widehat{{\rm
Var}}_{\mathcal{L}\mathcal{C}}\left[\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\boldsymbol{X},Y)-\tilde{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\boldsymbol{\delta},\tilde{\boldsymbol{\eta}})\right]+\frac{1}{\rho}\widehat{{\rm
Var}}_{\mathcal{L}\mathcal{C}}\left[\tilde{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\boldsymbol{\delta},\tilde{\boldsymbol{\eta}})\right].$
(6)

By design, $\tilde{\boldsymbol{\delta}}$ maximize the correlation between
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\boldsymbol{X},Y)$ and
$\tilde{\varphi}_{1}$ while minimizing the variance of $\tilde{\varphi}_{1}$
itself. Consequently, the new control variable $\tilde{\varphi}_{1}$ with
optimized $\tilde{\boldsymbol{\delta}}$ can further reduce the estimation
variance when $\mu$ is misspecified. We define the corresponding estimator in
a single BM scenario as

$\hat{\beta}_{1}=\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}-\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[\tilde{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\delta}},\tilde{\boldsymbol{\eta}})]+\widehat{{\rm
E}}_{\mathcal{L}\mathcal{M}}[\tilde{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\delta}},\tilde{\boldsymbol{\eta}})],$
(7)

and denote it as DEFUSE1.

Just like $\mu$, the IW model $w$ can also be constructed using general
learners (cross-fitting needed for complex ML) as long as the constant $1$
belongs to their feasible spaces or sets. Without loss of generality, we
impose this constraint as $w(\boldsymbol{x}^{[q]},y;\boldsymbol{0})\equiv 1$.
As will be shown later, this can ensure that when
$\mu(\boldsymbol{X}_{\mathcal{P}_{m}}\mid\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})$
is consistent with the true conditional distribution of
$\boldsymbol{X}_{\mathcal{P}_{m}}$, we have $\mu^{\prime}\asymp\mu$ and the
resulting $\tilde{\gamma}_{1}$ to be asymptotically equivalent with the
semiparametric efficient estimator (4). Similarly, it is also tentative to
make $w\equiv 0$ feasible in our construction since it can ensure
$\hat{\beta}_{1}$ to have no larger variance than the
$\mathcal{L}\mathcal{C}$-only estimator
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$, i.e., avoiding negative
effects of incorporating $\mathcal{L}\mathcal{M}$. For the choice of $w$, one
can see Assumption 5 and Lemma 2 for more rigorous details.

Notably, even when $\mu^{\prime}$ has negative values or an integral not equal
to $1$, $\hat{\beta}_{1}$ can still be asymptotic unbiased and of effectively
reduced variance. In general, (6) can be a non-convex problem with a high
computational cost to optimize due to the MC sampling. To handle it, one may
use stochastic gradient descendant (SGD) that generates a few MC samples for
each observation (i.e., with small $Q$) and makes gradient update accordingly
in multiple iterations. We shall remark on practical choices on $w$ and
$\boldsymbol{\delta}$ subject to our computational concerns.

######  Remark 2.

To avoid the high computational costs and potential over-fitting, we recommend
setting $w$ to have a simpler form than $\mu$. A simple but usually effective
choice is to use exponential tilting as described in (Johns,, 1988), where
$w(\boldsymbol{x},y;\boldsymbol{\delta})=\exp\\{\boldsymbol{\delta}^{\top}\xi(\boldsymbol{x},y)\\}$
and $\xi(\boldsymbol{x},y)$ is a $d$-dimensional basis with
$\boldsymbol{\delta}\in\mathbb{R}^{d}$. However, this makes (6) non-convex and
hence challenging to solve. Therefore, we propose further simplifying it to
$w(\boldsymbol{x},y;\boldsymbol{\delta})=1+\boldsymbol{\delta}^{\top}\xi(\boldsymbol{x},y)$,
which makes (6) a convex problem with an explicit-form solution. This form
closely resembles the exponential one when
$\boldsymbol{\delta}^{\top}\xi(\boldsymbol{x},y)$ is small. Additionally, when
$\xi(\boldsymbol{x},y)$ include a constant term of $1$ for intercept, both
$w\equiv 0$ and $w\equiv 1$ are feasible solutions.

###  2.3 Semi-supervised Estimation

After obtaining DEFUSE1, we further incorporate the $\mathcal{U}\mathcal{C}$
sample to enhance its efficiency. This also involves a data-adaptive
projection step that shares a similar spirit as in DEFUSE1. For the
$\mathcal{L}\mathcal{C}$-sample terms composing $\hat{\beta}_{1}$ in (7), an
interesting observation is that the score function of
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$ is orthogonal to
$\boldsymbol{X}$ when the outcome model (1) is correct. This orthogonality
leads to limited efficiency gains in the existing SSL methods for standard
GLMs, such as (Chakrabortty and Cai,, 2018), because their estimates are
score-function based. However, the projection of
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$ on
$(Y,\boldsymbol{X}_{\mathcal{P}_{m^{c}}})$,
$\tilde{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\delta}},\tilde{\boldsymbol{\eta}})$
is typically not orthogonal to $\boldsymbol{X}$. Thus, unlike those existing
SSL methods, $\boldsymbol{X}$ can partially explain the variance of
$\hat{\beta}_{1}$ through the projection, even when the outcome model is
correct. This promises the benefit of leveraging the large
$\mathcal{U}\mathcal{C}$ sample in our estimates.

Inspired by this, we propose to model the conditional distribution $Y$ given
$\boldsymbol{X}$, denoted as
$\nu(Y|\boldsymbol{X};\tilde{\boldsymbol{\theta}})$, which can be either
directly specified as the outcome model (1) with some extra nuisance
parameters (e.g., variance of Gaussian noise), or estimated using ML as
discussed in Remark 1. Similar as in Section 2.2, we introduce an IW function

$h(\boldsymbol{X},Y;\boldsymbol{\zeta})=\frac{\nu^{\prime}(Y|\boldsymbol{X};\boldsymbol{\zeta},\tilde{\boldsymbol{\theta}})}{\nu(Y|\boldsymbol{X};\tilde{\boldsymbol{\theta}})},$

where
$\nu^{\prime}(Y|\boldsymbol{X};\boldsymbol{\zeta},\tilde{\boldsymbol{\theta}})$
stands for the calibrated conditional model for $Y$ to improve the efficiency
when $\nu(Y|\boldsymbol{X};\tilde{\boldsymbol{\theta}})$ is misspecified. The
optimal $\boldsymbol{\zeta}$ can be obtained by

$\tilde{\boldsymbol{\zeta}}=\arg\min_{\boldsymbol{\zeta}}\widehat{{\rm
Var}}_{\mathcal{L}\mathcal{C}}\left[\tilde{S}_{1}(\boldsymbol{X},Y)-\tilde{\varphi}_{2}^{\prime}(\boldsymbol{X};\boldsymbol{\zeta},\tilde{\boldsymbol{\theta}})\right],$
(8)

where
$\tilde{S}_{1}(\boldsymbol{X},Y)=\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\boldsymbol{X},Y)-\tilde{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\delta}},\tilde{\boldsymbol{\eta}})$
is the empirical score function of DEFUSE1 in the $\mathcal{L}\mathcal{C}$
sample. The term
$\tilde{\varphi}_{2}^{\prime}(\boldsymbol{X};\boldsymbol{\zeta},\tilde{\boldsymbol{\theta}})$
is constructed as
$Q^{-1}\sum_{q=1}^{Q}h(\boldsymbol{X},Y^{[q]};\boldsymbol{\zeta})\tilde{S}_{1}(\boldsymbol{X},Y^{[q]}),$
where $\\{Y^{[q]}:q=1,2,\ldots,Q\\}$ are i.i.d. random drawn from the
preliminary model $\nu(Y|\boldsymbol{X};\tilde{\boldsymbol{\theta}})$.

Comparing to (6), the objective function in (8) does not include the variance
term of $\tilde{\varphi}_{2}^{\prime}$. That is because the variance of
$\tilde{\varphi}_{2}^{\prime}$ on $\mathcal{U}\mathcal{C}$ is negligible in
the SS setting with $N_{\mathcal{U}\mathcal{C}}\gg
N_{\mathcal{L}\mathcal{C}}$. Finally, the DEFUSE estimator in the single BM
scenario is constructed as

$\hat{\beta}_{2}=\hat{\beta}_{1}-\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[\tilde{\varphi}_{2}^{\prime}(\boldsymbol{X};\boldsymbol{\zeta},\tilde{\boldsymbol{\theta}})]+\widehat{{\rm
E}}_{\mathcal{U}\mathcal{C}}[\tilde{\varphi}_{2}^{\prime}(\boldsymbol{X};\boldsymbol{\zeta},\tilde{\boldsymbol{\theta}})].$
(9)

Again, we impose $h(\boldsymbol{X},Y;\boldsymbol{0})=1$ so that
$\hat{\beta}_{2}$ can achieve semiparametric efficiency in this SS step when
the imputation distribution
$\nu(\cdot|\boldsymbol{X};\tilde{\boldsymbol{\theta}})$ is consistent with the
truth. Also, letting constant $0$ be in the feasible space of $h$ can ensure
that $\tilde{\varphi}_{2}^{\prime}=0$ is feasible and $\hat{\beta}_{2}$ has a
smaller or equal variance than $\hat{\beta}_{1}$. For specification of $h$, we
suggest using similar strategies as proposed in Remark 2. We end our
exposition for the single BM and SS setting with a remark on our constructing
strategy.

######  Remark 3.

Our framework starts from leveraging the $\mathcal{L}\mathcal{M}$ sample by
projecting $\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$ to
$(Y,\boldsymbol{X}_{\mathcal{P}_{m^{c}}})$ and then incorporates
$\mathcal{U}\mathcal{C}$ to further enhance the resulting DEFUSE1 estimator.
Reversing the order of these two projection steps is likely to result in a
less efficient estimator because the score of
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$ tends to be orthogonal to
$\boldsymbol{X}$ and directly using $\mathcal{U}\mathcal{C}$ on it may not be
helpful when the outcome model (1) is strictly or approximately correct. Also,
we notice that $\varphi^{*}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)$ is
orthogonal to $\boldsymbol{X}_{\mathcal{P}_{m^{c}}}$ when the models for $Y$
and $\boldsymbol{X}_{\mathcal{P}_{m}}$ are both correct and its empirical
version
$\tilde{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\delta}},\tilde{\boldsymbol{\eta}})$
tends to be poorly explained by $\boldsymbol{X}_{\mathcal{P}_{m^{c}}}$ on the
$\mathcal{L}\mathcal{M}$ data in our finite sample studies. Therefore, we do
not elaborate the variance contributed by the $\mathcal{L}\mathcal{M}$ sample
in our SS procedure to save computation costs.

###  2.4 DEFUSE in the Multiple BM Scenario

In this section, we extend the proposed DEFUSE method to the scenario of
$R\geq 2$ data sets $\\{\mathcal{L}\mathcal{M}_{r}:r=1,\ldots,R\\}$ with
different BM sets denoted as $\\{\mathcal{P}_{m}^{r}:r=1,\ldots,R\\}$ and
sample sizes $\\{N_{\mathcal{L}\mathcal{M}_{r}}:1,\ldots,R\\}$; see the right
panel of Figure 1. Same as the single BM scenario, we start with the standard
$\mathcal{L}\mathcal{C}$-only estimator
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$ and leverage
$\mathcal{L}\mathcal{M}_{r}$’s to improve its efficiency. To this end, we
first implement data-adaptive calibration proposed in Section 2.2 with each
$\mathcal{L}\mathcal{M}_{r}$ separately, to obtain the control function
$\tilde{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)$
analog to $\tilde{\varphi}_{1}^{\prime}$ in the single BM case. Here we drop
the notations of the nuisance parameters in $\tilde{\varphi}_{1,r}^{\prime}$
for simplification.

Denote by
$\tilde{\boldsymbol{\varphi}}_{1}^{\bullet}(\boldsymbol{X},Y)=\\{\tilde{\varphi}_{1,1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{1}},Y),\ldots,\tilde{\varphi}_{1,R}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{R}},Y)\\}^{\top}$.
Then we assemble the multiple control functions to form the DEFUSE1 estimator
as

$\hat{\beta}_{1}^{\bullet}=\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}+\sum_{r=1}^{R}\tilde{a}_{r}\left\\{\widehat{{\rm
E}}_{\mathcal{L}\mathcal{M}_{r}}\left[\tilde{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)\right]-\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}\left[\tilde{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)\right]\right\\},$
(10)

where the ensemble weights $\tilde{a}_{r}$’s in $\tilde{\boldsymbol{a}}$ are
solved from the quadratic programming:

$\min_{\boldsymbol{a}=(a_{1},\ldots,a_{R})^{\top}}~{}\widehat{{\rm
Var}}_{\mathcal{L}\mathcal{C}}\left[\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\boldsymbol{X},Y)-\boldsymbol{a}^{\top}\tilde{\boldsymbol{\varphi}}_{1}^{\bullet}(\boldsymbol{X},Y)\right]+\sum_{r=1}^{R}\frac{a_{r}^{2}}{\rho_{r}}\widehat{{\rm
Var}}_{\mathcal{L}\mathcal{M}_{r}}\left[\tilde{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)\right],$
(11)

and
$\rho_{r}=\min\\{N_{\mathcal{L}\mathcal{M}_{r}}/N_{\mathcal{L}\mathcal{C}},C\\}$.
The goal of (11) is to find the optimal allocation $\tilde{\boldsymbol{a}}$
minimizing the asymptotic variance of DEFUSE1 as formed in (10) that augments
the preliminary $\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$ with a
linear combination of the variance reduction term coming from every
$\mathcal{L}\mathcal{M}_{r}\cup\mathcal{L}\mathcal{C}$. Correspondingly, the
objective function in (11) is set as the empirical variance of such an
estimator with weights $\tilde{\boldsymbol{a}}$; see our asymptotic analysis
for details.

One can view (11) as a ridge regression with each coefficient $a_{r}$
penalized according to the sample size $N_{\mathcal{L}\mathcal{M}_{r}}$ and
variance of $\tilde{\varphi}_{1,r}^{\prime}$. For certain BM patterns,
$\tilde{\boldsymbol{\varphi}}_{1}^{\bullet}(\boldsymbol{X},Y)$ can have highly
correlated or collinear elements. In this case, the “ridge penalty” on
$\tilde{\boldsymbol{a}}$ can actually ensure the non-singularity and stability
of (11). Thus, we assign an upper bound $C$ to each $\rho_{r}$ in our
definition so that $\rho_{r}^{-1}$ cannot be too small to make this penalty
ineffective even when $N_{\mathcal{L}\mathcal{M}_{r}}\gg
N_{\mathcal{L}\mathcal{C}}$.

Similarly to the single BM scenario, we incorporate the
$\mathcal{U}\mathcal{C}$ sample to further reduce the variance of
$\hat{\beta}_{1}^{\bullet}$ using the same idea in Section 2.3. Note that the
empirical score function of $\hat{\beta}_{1}^{\bullet}$ on
$\mathcal{L}\mathcal{C}$ can be written as
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{S}}(\boldsymbol{X},Y)-\tilde{\boldsymbol{a}}^{\top}\tilde{\boldsymbol{\varphi}}_{1}^{\bullet}(\boldsymbol{X},Y)$.
Naturally, we use it to replace the score $\tilde{S}_{1}(\boldsymbol{X},Y)$ in
(8) and solve it to obtain corresponding calibration parameters. Then we use
this solution to impute the score and augment $\hat{\beta}_{1}^{\bullet}$ in
the same form as (9). This procedure can be regarded as a straightforward
extension of our proposal in Section 2.3.

We summarize the main steps of the DEFUSE approach in Algorithm 1. We notice
that some recent work like Kundu and Chatterjee, (2023) highlighted practical
needs to protect the privacy of individual-level data when fusing the data
sets
$\mathcal{L}\mathcal{C},\mathcal{L}\mathcal{M}_{1},\ldots,\mathcal{L}\mathcal{M}_{R}$
coming from different data institutions. This has been achieved by
transferring summary data only (e.g., mean vectors, model coefficients, and
covariance matrices) across the sites, which is often referred as the
DataSHIELD framework or constraint Wolfson et al., (2010). Our method can
naturally accommodate this situation as well. In specific, note that Steps 1,
2.1, and 3.1 in Algorithm 1 only rely on the $\mathcal{L}\mathcal{C}$ sample
while Steps 2.2 and 3.2 require summary data (i.e., the empirical mean of the
imputed score functions) from each $\mathcal{L}\mathcal{M}_{r}$ or
$\mathcal{U}\mathcal{C}$. In practice, $\mathcal{L}\mathcal{C}$ and
$\mathcal{U}\mathcal{C}$ are likely to be stored in the same site. Thus, under
the data-sharing constraint, the implementation of DEFUSE requires usually one
and at most two rounds of communication between $\mathcal{L}\mathcal{C}$ and
other sites, which is communication efficient and user-friendly.

Algorithm 1 DEFUSE

Samples $\mathcal{L}\mathcal{C}$,
$\\{\mathcal{L}\mathcal{M}_{r}:r=1,\ldots,R\\}$ ($R\geq 1$), and
$\mathcal{U}\mathcal{C}$; a loading vector $\boldsymbol{c}$.

[Step 1] Solving Equation (3) on $\mathcal{L}\mathcal{C}$ data to obtain
$\tilde{\boldsymbol{\gamma}}$ and compute $\tilde{\boldsymbol{S}}$.

[Step 2.1] For each $r\in\\{1,\ldots,R\\}$: learn the conditional distribution
$\tilde{\mu}_{\boldsymbol{X}_{\mathcal{P}_{m}^{r}}}$; then solve for the data-
adaptive calibration function via (6) on $\mathcal{L}\mathcal{C}$ (with
$\mathcal{P}_{m}$ set as each $\mathcal{P}_{m}^{r}$).

[Step 2.2] When $R=1$, construct the DEFUSE1 estimator $\hat{\beta}_{1}$ using
(7); when $R\geq 2$, obtain the optimal allocation $\tilde{\boldsymbol{a}}$ by
(11) and use it to construct $\hat{\beta}_{1}^{\bullet}$ as in (10).

[Step.3.1] Learn the conditional distribution
$\nu(Y|\boldsymbol{X};\tilde{\boldsymbol{\theta}})$ on
$\mathcal{L}\mathcal{C}$. Obtain the data-adaptive calibration function via
(8) on $\mathcal{L}\mathcal{C}$.

[Step.3.2] Construct the final DEFUSE estimator by (9) or in a similar way
when $R\geq 2$ as discussed in Section 2.4.

###  2.5 Extension on Auxiliary Covariates

We shall end the method section with brief description and discussion on the
extension of DEFUSE to accommodate the situation with auxiliary features
frequently encountered in many real-world applications. Consider the scenario
with $\boldsymbol{X}=\\{\boldsymbol{A},\boldsymbol{W}\\}$ consisting of some
primary risk factors $\boldsymbol{A}$ and auxiliary or surrogate features
$\boldsymbol{W}$. Instead of the outcome model against the whole
$\boldsymbol{X}$, now the focus is on $Y\sim\boldsymbol{A}$ with the
population equation ${\rm
E}\left[\boldsymbol{A}\left\\{Y-g\left(\boldsymbol{\gamma}^{\top}\boldsymbol{A}\right)\right\\}\right]=\mathbf{0}$.
Meanwhile, the missing structure remains the same as described in Section 1.2,
i.e., one or multiple subsets of covariates $\boldsymbol{X}_{\mathcal{P}_{m}}$
partially missing on the labeled sample and $Y$ being unlabeled on
$\mathcal{U}\mathcal{C}$ with the complete $\boldsymbol{X}$. Here, covariates
$\boldsymbol{A}$ can be either included by, intersected with, or disjoint with
$\boldsymbol{X}_{\mathcal{P}_{m}}$.

This is a common setup in real-world studies. For example, on EMR linked
biobank data, the interest lies in the relationship between certain genetic
variants $\boldsymbol{A}$ and disease status $Y$. $\boldsymbol{W}$ is taken as
EMR proxies (e.g., relevant codified features) for $Y$ and acts as auxiliary
features not included in the genetic risk model. In clinical or
epidemiological studies, $\boldsymbol{A}$ is the treatment or key risk
factors, $Y$ is a long-term outcome, and $\boldsymbol{W}$ can be some early
endpoints or surrogate outcomes (e.g., tumor response rates) occurring and
observed post the baseline with only $\boldsymbol{A}$. In this case, one is
interested in $Y\sim\boldsymbol{A}$ to make clinical decision at the baseline
without $\boldsymbol{W}$. In both examples, the auxiliary $\boldsymbol{W}$
could still be informative to $Y$ and, thus, included to boost the statistical
efficiency.

Our framework can be naturally and easily extended to accommodate the setting
introduced above. In specific, the only change to make is to replace the Step
1 in Algorithm 1 with the one solving for $\tilde{\boldsymbol{{\gamma}}}$ from
$\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[\boldsymbol{A}\\{Y-g(\boldsymbol{\gamma}^{\top}\boldsymbol{A})\\}]=\mathbf{0}$
and deriving the corresponding score function
$\tilde{\boldsymbol{S}}_{i}=\tilde{\boldsymbol{H}}^{-1}\boldsymbol{A}_{i}\\{Y_{i}-g(\tilde{\boldsymbol{{\gamma}}}^{\top}\boldsymbol{A}_{i})\\}$
where $\tilde{\boldsymbol{H}}=\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[g(\tilde{\boldsymbol{{\gamma}}}^{\top}\boldsymbol{A})\boldsymbol{A}\boldsymbol{A}^{\top}]$.
All the conditional modeling, imputation, and calibration procedures in the
following Steps 2 and 3 in Algorithm 1 can be conducted in exactly the same
forms as presented in previous sections. In this case, $\boldsymbol{W}$ could
often be of high-dimensionality as considered in Zhou et al., (2024). Our
theoretical framework to be established in the next section can also be
applied to this setting.

##  3 Asymptotic Analysis

###  3.1 Notations and Assumptions

In this section, we conduct theoretical analysis on the asymptotic normality
and efficiency gain of DEFUSE. Justification of all the results presented in
this section can be found in Appendix. We begin with notations and main
assumptions. Let ${\rm P}(\boldsymbol{Z})$ denote the true probability density
or mass function of random variables $\boldsymbol{Z}$ and ${\rm
P}(\boldsymbol{Z}_{1}|\boldsymbol{Z}_{2})$ be the conditional probability
function, ${\rm Cov}(\cdot)$ be the population covariance operator, and ${\rm
Diag}(u_{k})_{k=1,\ldots,d}$ be the $d\times d$ diagonal matrix with its
$k$-th diagonal entry being $u_{k}$. For simplicity of presentation, denote by
$n=N_{\mathcal{L}\mathcal{C}}$. For two sequences $a_{n}$ and $b_{n}$ with
index $n\rightarrow\infty$, we say $a_{n}<\infty$ when there exists some
constant $M>0$ such that $a_{n}<M$, $a_{n}=O(b_{n})$ when
$\lim_{n\rightarrow\infty}|a_{n}/b_{n}|<\infty$, $a_{n}=o(b_{n})$ for
$\lim_{n\rightarrow\infty}a_{n}/b_{n}=0$, and $a_{n}=o_{\sf P}(b_{n})$ or
$a_{n}=O_{\sf P}(b_{n})$ for $a_{n}=o(b_{n})$ or $a_{n}=O(b_{n})$ with the
probability approaching $1$. We say $f(x)$ is a Lipschitz function if there
exists a constant $A>0$ such that for any $x_{1}$ and $x_{2}$ in the domain of
$f$, $|f(x_{1})-f(x_{2})|\leq A|x_{1}-x_{2}|$.

We first introduce the assumptions for analyzing DEFUSE in the single BM
scenario. Assume that $\rho=N_{\mathcal{L}\mathcal{M}}/n$ converges to some
$\bar{\rho}\in[0,\infty]$,
$n+N_{\mathcal{L}\mathcal{M}}=o(N_{\mathcal{U}\mathcal{C}})$, and the number
of MC sample $Q>n^{\epsilon}$ with an arbitrary constant $\epsilon>0$ that can
be set small to save computational costs without having impacts on the
asymptotic properties of DEFUSE.

######  Assumption 1.

The derivative of the link function $\dot{g}(u)$ is Lipschitz;
${\boldsymbol{\gamma}}$ belongs to a compact space including
$\bar{\boldsymbol{\gamma}}$ as its interior; and the Hessian
$\bar{\boldsymbol{H}}$ has all its eigenvalues staying away from $0$ and
$\infty$. The joint distribution of $(Y,\boldsymbol{X})$ has a continuously
differentiable density for the continuous components and satisfies that
$\mathrm{E}\left[Y^{4}+\|\boldsymbol{X}\|_{2}^{4}+\\{g(\boldsymbol{X}^{\top}\bar{\boldsymbol{\gamma}})\\}^{2}+\\{\dot{g}(\boldsymbol{X}^{\top}\bar{\boldsymbol{\gamma}})\\}^{2}\right]<\infty.$

######  Assumption 2.

There exists some population level parameters $\bar{\boldsymbol{\eta}}$ and
$\bar{\boldsymbol{\theta}}$ such that

$\displaystyle{\rm E}_{\bar{\mu}_{\mathcal{P}_{m}}\times{\rm
P}_{Y,\mathcal{P}_{m^{c}}}}\left[\Delta_{\mu}(\tilde{\boldsymbol{\eta}},\bar{\boldsymbol{\eta}})\right]^{2}=o_{p}(1);\quad{\rm
E}_{\bar{\nu}_{Y}\times{\rm
P}_{\boldsymbol{X}}}[\Delta_{\nu}(\tilde{\boldsymbol{\theta}},\bar{\boldsymbol{\theta}})]^{2}=o_{p}(1);$
$\displaystyle{\rm E}_{\bar{\mu}_{\mathcal{P}_{m}}\times{\rm
P}_{Y,\mathcal{P}_{m^{c}}}}\left[Y^{4}+\|\boldsymbol{X}\|_{2}^{4}+\\{g(\boldsymbol{X}^{\top}\bar{\boldsymbol{\gamma}})\\}^{2}+\\{\dot{g}(\boldsymbol{X}^{\top}\bar{\boldsymbol{\gamma}})\\}^{2}\right]<\infty;$
$\displaystyle{\rm E}_{\bar{\nu}_{Y}\times{\rm
P}_{\boldsymbol{X}}}\left[Y^{4}+\|\boldsymbol{X}\|_{2}^{4}+\\{g(\boldsymbol{X}^{\top}\bar{\boldsymbol{\gamma}})\\}^{2}+\\{\dot{g}(\boldsymbol{X}^{\top}\bar{\boldsymbol{\gamma}})\\}^{2}\right]<\infty,$

where we denote by

$\displaystyle\Delta_{\mu}(\tilde{\boldsymbol{\eta}},\bar{\boldsymbol{\eta}})=\left|1-\frac{{\mu}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\tilde{\boldsymbol{\eta}})}{{\mu}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\bar{\boldsymbol{\eta}})}\right|;\quad\Delta_{\nu}(\tilde{\boldsymbol{\theta}},\bar{\boldsymbol{\theta}})=\left|1-\frac{\nu(Y|\boldsymbol{X};\tilde{\boldsymbol{\theta}})}{\nu(Y|\boldsymbol{X};\bar{\boldsymbol{\theta}})}\right|;$
$\displaystyle\bar{\mu}_{\mathcal{P}_{m}}\times{\rm
P}_{Y,\mathcal{P}_{m^{c}}}={\mu}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\bar{\boldsymbol{\eta}})\times{\rm
P}(Y,\boldsymbol{X}_{\mathcal{P}_{m^{c}}});\quad\bar{\nu}_{Y}\times{\rm
P}_{\boldsymbol{X}}=\nu(Y|\boldsymbol{X};\bar{\boldsymbol{\theta}})\times{\rm
P}(\boldsymbol{X}).$

######  Assumption 3.

There exists some population level parameters $\bar{\boldsymbol{\delta}}$ and
$\bar{\boldsymbol{\zeta}}$ such that

$\displaystyle{\rm E}_{\bar{\mu}_{\mathcal{P}_{m}}\times{\rm
P}_{Y,\mathcal{P}_{m^{c}}}}[\Delta_{w}(\tilde{\boldsymbol{\delta}},\bar{\boldsymbol{\delta}})]^{2}=o_{p}(1);\quad{\rm
E}_{\bar{\nu}_{Y}\times{\rm
P}_{\boldsymbol{X}}}[\Delta_{h}(\tilde{\boldsymbol{\zeta}},\bar{\boldsymbol{\zeta}})]^{2}=o_{p}(1),$

where
$\Delta_{w}(\tilde{\boldsymbol{\delta}},\bar{\boldsymbol{\delta}})=|w(\boldsymbol{X},Y;\tilde{\boldsymbol{\delta}})-w(\boldsymbol{X},Y;\bar{\boldsymbol{\delta}})|$
and
$\Delta_{h}(\tilde{\boldsymbol{\zeta}},\bar{\boldsymbol{\zeta}})=|h(\boldsymbol{X},Y;\tilde{\boldsymbol{\zeta}})-h(\boldsymbol{X},Y;\bar{\boldsymbol{\zeta}})|$.

######  Remark 4.

Assumption 1 include mild and common regularity conditions for M-estimation
(Van der Vaart,, 2000). Assumptions 2 and 3 convey a message that to establish
the asymptotic normality of DEFUSE, all the nuisance estimators used for
imputation and calibration need to be consistent with some limiting values.
Compared to the double machine learning (DML) framework (Chernozhukov et al.,,
2018), we consider a more stringent MCAR regime, under which the estimation
forms (7) and (9) automatically correcting the excessive bias of the nuisance
estimators (with cross-fitting for general ML). Consequently, our assumption
on the quality of the nuisance estimators is substantially weaker than DML
from two aspects. First, we only require them to converge to some population
limits but not necessarily the true distributions required by DML. Second, we
assume an $o_{p}(1)$-convergence rate essentially slower than the
$o_{p}(n^{-1/4})$ rate assumed in DML.

######  Remark 5.

Assumption 2 on the consistency of the preliminary imputation models is
reasonable and justifiable across broad settings with
$\bar{\mu}_{\mathcal{P}_{m}}$ and $\bar{\nu}_{Y}$ constructed using either
classic regression or modern ML tools. One can refer to comprehensive
literature for the consistency of various learning approaches under
potentially misspecification, e.g., parametric GLMs (Tian et al.,, 2007,
e.g.); kernel methods (Zhang et al.,, 2023, e.g.); high-dimensional sparse
regression (Negahban et al.,, 2012, e.g.); and random forest (Athey et al.,,
2019, e.g.).

For the calibration estimators $\tilde{\boldsymbol{\delta}}$ and
$\tilde{\boldsymbol{\zeta}}$, we present Proposition 1 to justify their
consistency in Assumption 3 under a parametric form satisfying mild regularity
conditions.

######  Proposition 1.

Under the specification
$w(\boldsymbol{x},y;\boldsymbol{\delta})=1+\boldsymbol{\delta}^{\top}\xi(\boldsymbol{x},y)$
and
$h(\boldsymbol{x},y;\boldsymbol{\zeta})=1+\boldsymbol{\zeta}^{\top}\phi(\boldsymbol{x},y)$
with some prespecified basis functions $\xi$ and $\phi$, Assumptions 1 and 2,
as well as the regularity Assumption A1 in Appendix, our Assumption 3 holds.

Based on the population-level parameters introduced in Assumptions 2 and 3, we
further define the population control and score functions as

$\begin{split}&\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)=\frac{\bar{\rho}}{1+\bar{\rho}}{\rm
E}_{\bar{\mu}_{\mathcal{P}_{m}}(\cdot|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\bar{\boldsymbol{\eta}})}[w(\boldsymbol{X},Y;\bar{\boldsymbol{\delta}})\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}\mid\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y];\\\
&\bar{S}_{1}(\boldsymbol{X},Y)=\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}-\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y);\quad\bar{\varphi}_{2}^{\prime}(\boldsymbol{X})={\rm
E}_{\nu(\cdot|\boldsymbol{X};\bar{\boldsymbol{\theta}})}h(\boldsymbol{X},Y;\bar{\boldsymbol{\zeta}})\bar{S}_{1}(\boldsymbol{X},Y),\end{split}$
(12)

which will be used later to characterize the asymptotic properties of DEFUSE.
Finally, we introduce an additional Assumption 4 used for asymptotic analysis
of the multiple BM scenario in Section 2.4.

######  Assumption 4.

In the multiple BM setting, it is satisfied: (i) Our data and construction for
all BM components in $\\{\mathcal{L}\mathcal{M}_{r}:r=1,\ldots,R\\}$ satisfy
Assumptions 1–3 for $\mathcal{L}\mathcal{M}$ in the single BM scenario; (ii)
Each
$\rho_{r}=\min\\{N_{\mathcal{L}\mathcal{M}_{r}}/N_{\mathcal{L}\mathcal{C}},C\\}$
converges to some $\bar{\rho}_{r}$ for $r=1,\ldots,R$; (iii) Given that (i)
and (ii) hold, define
$\bar{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{1}},Y)$ as
the limiting function for each $\mathcal{L}\mathcal{M}_{r}$ analog to
$\bar{\varphi}_{1}^{\prime}$ in (12), and denote by
$\bar{\boldsymbol{\varphi}}_{1}^{\bullet}(\boldsymbol{X},Y)=\\{\bar{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{1}},Y)\\}_{r=1,\ldots,R}^{\top}$.
Then matrix

${\rm
Cov}\left[\bar{\boldsymbol{\varphi}}_{1}^{\bullet}(\boldsymbol{X},Y)\right]+{\rm
Diag}\left\\{\frac{1}{\bar{\rho}_{r}}{\rm
Var}\left[\bar{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{1}},Y)\right]\right\\}_{r=1,\ldots,R}$

has all its eigenvalues staying away from $0$.

Assumption 4 is introduced to guarantee that (11) can produce a consistent
estimation of the unique population-level optimal allocation

$\bar{\boldsymbol{a}}=(\bar{a}_{1},\ldots,\bar{a}_{R})^{\top}=\arg\min_{\boldsymbol{a}}~{}{\rm
Var}\left[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)-\boldsymbol{a}^{\top}\bar{\boldsymbol{\varphi}}_{1}^{\bullet}(\boldsymbol{X},Y)\right]+\sum_{r=1}^{R}\frac{a_{r}^{2}}{\bar{\rho}_{r}}{\rm
Var}\left[\bar{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)\right].$

This assumption holds when every data set $\mathcal{L}\mathcal{M}_{r}$ is
informative enough to
$\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)$ and the variance
of $\bar{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{1}},Y)$
is not too small.

###  3.2 Asymptotic Normality

In this section, we present the asymptotic properties of the proposed
estimators, with their efficiency analyzed and articulated in the next
section. We start from the standard M-estimation result of the
$\mathcal{L}\mathcal{C}$-only estimator $\tilde{\boldsymbol{\gamma}}$ in Lemma
1 that is already well-established (Van der Vaart,, 2000, e.g.).

######  Lemma 1.

Under Assumptions 1–3,
$\tilde{\boldsymbol{\gamma}}\overset{p}{\to}\boldsymbol{\bar{\gamma}}$ and
$n^{1/2}(\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}-\boldsymbol{c}^{\top}\boldsymbol{\bar{\gamma}})=n^{1/2}\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}]+o_{\sf
P}(1),$ weakly converges to the normal distribution ${\rm N}\\{0,{\rm
Var}(\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}_{i})\\}$.

Next, we establish the consistency and asymptotic normality for the DEFUSE1
and DEFUSE estimators in the single BM scenario.

######  Theorem 1.

Under Assumptions 1–3, we have both
$\hat{\beta}_{1},\hat{\beta}_{2}\overset{p}{\to}\boldsymbol{c}^{\top}\boldsymbol{\bar{\gamma}}$
and

$\displaystyle\sqrt{n}(\hat{\beta}_{1}-\boldsymbol{c}^{\top}\boldsymbol{\bar{\gamma}})$
$\displaystyle={\sqrt{n}}\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}\left[\bar{S}_{1}(\boldsymbol{X},Y)\right]+{\sqrt{n}}\widehat{{\rm
E}}_{\mathcal{L}\mathcal{M}}\left[\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)\right]+o_{\sf
P}(1)$
$\displaystyle\sqrt{n}(\hat{\beta}_{2}-\boldsymbol{c}^{\top}\boldsymbol{\bar{\gamma}})$
$\displaystyle={\sqrt{n}}\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}\left[\bar{S}_{1}(\boldsymbol{X},Y)-\bar{\varphi}_{2}^{\prime}(\boldsymbol{X})\right]+{\sqrt{n}}\widehat{{\rm
E}}_{\mathcal{L}\mathcal{M}}\left[\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)\right]+o_{\sf
P}(1),$

which weakly converge to ${\rm N}\\{0,{\rm
Var}[\bar{S}_{1}(\boldsymbol{X},Y)]+{\rm
Var}[\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)]/\bar{\rho}\\}$
and ${\rm N}\\{0,{\rm
Var}[\bar{S}_{1}(\boldsymbol{X},Y)-\bar{\varphi}_{2}^{\prime}(\boldsymbol{X})]+{\rm
Var}[\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)]/\bar{\rho}\\}$
respectively.

Theorem 1 states that our DEFUSE estimators are $n^{-1/2}$-consistent,
asymptotically unbiased and normal no matter the imputation distributions are
correct or not and even when the nuisance estimators have slow convergence
rates. The efficiency of $\hat{\beta}_{1}$ and $\hat{\beta}_{2}$ will be
analyzed and discussed based on Theorem 1 later. As an extension, we also
establish an asymptotic theorem for $\hat{\beta}_{1}^{\bullet}$ in the
multiple BM setting. Its downstream SS estimation can be justified in
basically the same way as the analysis on $\hat{\beta}_{2}$ in Theorem 1 so we
omit this part to avoid redundancy.

######  Theorem 2.

Under Assumption 4,
$\hat{\beta}_{1}^{\bullet}\overset{p}{\to}\boldsymbol{c}^{\top}\boldsymbol{\bar{\gamma}}$
and

$\sqrt{n}(\hat{\beta}_{1}^{\bullet}-\boldsymbol{c}^{\top}\boldsymbol{\bar{\gamma}})={\sqrt{n}}\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}\left[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)-\bar{\boldsymbol{a}}^{\top}\bar{\boldsymbol{\varphi}}_{1}^{\bullet}(\boldsymbol{X},Y)\right]+{\sqrt{n}}\sum_{r=1}^{R}\bar{a}_{r}\widehat{{\rm
E}}_{\mathcal{L}\mathcal{M}_{r}}\left[\bar{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)\right]+o_{\sf
P}(1),$

weakly weakly converges to the normal distribution of mean $0$ and variance

${\rm
Var}\left[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)-\bar{\boldsymbol{a}}^{\top}\bar{\boldsymbol{\varphi}}_{1}^{\bullet}(\boldsymbol{X},Y)\right]+\sum_{r=1}^{R}\frac{\bar{a}_{r}^{2}}{\bar{\rho}_{r}}{\rm
Var}\left[\bar{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)\right].$

###  3.3 Relative Efficiency

Now we shall demonstrate the relative efficiency of the proposed DEFUSE
estimators compared to a broad set of existing ones including: (a) the
$\mathcal{L}\mathcal{C}$-only estimator
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$; (b) the standard
semiparametric efficient estimator (Robins et al.,, 1994) based on (4) with
$\varphi^{*}$ simply imputed by $\tilde{\varphi}_{1}$; (c) BM estimators
constructed in different strategies (Xue and Qu,, 2021; Zhao et al.,, 2023),
as well as Xue et al., (2021); Song et al., (2024) also considering the SS
setting; (d) SS methods only incorporating the $\mathcal{U}\mathcal{C}$ sample
(Gronsbell et al.,, 2022, e.g.). In addition, we will remark on the connection
of our data-adaptive calibration strategy to some recent work for the SS
inference. We begin with the assumptions and lemma on the feasible set and
regularity of the calibration models.

######  Assumption 5.

The population-level parameters $\bar{\boldsymbol{\delta}}$ and
$\bar{\boldsymbol{\zeta}}$ defined in Assumption 2 are the unique minimizer of
the variance functions $V_{1}({\boldsymbol{\delta}})$ and
$V_{2}(\boldsymbol{\zeta})$ respectively, where

$\displaystyle V_{1}({\boldsymbol{\delta}})=$ $\displaystyle{\rm
Var}\left[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)-\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\boldsymbol{\delta},\bar{\boldsymbol{\eta}})\right]+{\bar{\rho}}^{-1}{\rm
Var}\left[\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\boldsymbol{\delta},\bar{\boldsymbol{\eta}})\right];$
$\displaystyle V_{2}(\boldsymbol{\zeta})=$ $\displaystyle{\rm
Var}\left[\bar{S}_{1}(\boldsymbol{X},Y)-\bar{\varphi}_{2}^{\prime}(\boldsymbol{X};\boldsymbol{\zeta},\bar{\boldsymbol{\theta}})\right]+{\bar{\rho}}^{-1}{\rm
Var}\left[\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\bar{\boldsymbol{\delta}},\bar{\boldsymbol{\eta}})\right].$

######  Assumption 6.

Constant $1$ belongs to the feasible sets of the calibration models and we can
set $w(\boldsymbol{x},y;\boldsymbol{0})\equiv 1$ and
$h(\boldsymbol{x},y;\boldsymbol{0})\equiv 1$ without loss of generality. If
some $f(\boldsymbol{x},y)$ belongs to the feasible set of
$w(\boldsymbol{x},y;\boldsymbol{\delta})$ or
$h(\boldsymbol{x},y;\boldsymbol{\zeta})$, then the function
$cf(\boldsymbol{x},y)$ also belongs to that feasible set for any
$c\in\mathbb{R}$.

######  Lemma 2.

Assumptions 5 and 6 will hold when all assumptions used in Proposition 1 hold
and the forms of the calibration models are taken as
$w(\boldsymbol{x},y;\boldsymbol{\delta})=1+\boldsymbol{\delta}^{\top}\xi(\boldsymbol{x},y)$
and
$h(\boldsymbol{x},y;\boldsymbol{\zeta})=1+\boldsymbol{\zeta}^{\top}\phi(\boldsymbol{x},y)$
with $\xi(\boldsymbol{x},y)$ and $\phi(\boldsymbol{x},y)$ including constant
$1$.

Though imposing linear specification on $\omega$ and $h$, Lemma 2 can be
generalized to other forms like
$w(\boldsymbol{x},y;\delta_{1},\boldsymbol{\delta}_{2})=\delta_{1}\exp\\{\boldsymbol{\delta}_{2}^{\top}\xi(\boldsymbol{x},y)\\}$
but this may encounter local minima issues in practice. The regularity
assumptions on $\xi$ and $\phi$ in Lemma 2, i.e., Assumption A1 is mild and
tends to hold under proper choices.

######  Remark 6.

We notice that some recent approaches for power enhancement under misspecified
nuisance models leverage similar properties as Assumptions 5 and 6. For
example, Gronsbell et al., (2022) combined the supervised and SS estimators
through their optimal convex allocation with minimum variance. We can show
that their strategy is essentially similar to setting our
$h(\boldsymbol{x},y;\boldsymbol{\zeta})=1+\zeta_{1}\in[0,1]$ with only one
parameter $\zeta_{1}$ to fit. This easily satisfies above-introduced
assumptions on $h(\cdot)$ (with $c\in[0,1]$ in Assumption 6). As a result, our
following corollaries about efficiency may apply in their case. Similar
constructions could be found also in Miao et al., (2023) and others.
Consequently, when taking $h(\boldsymbol{x},y;\boldsymbol{\zeta})$ as a nested
model of the most simple one $1+\zeta_{1}$, e.g.,
$h(\boldsymbol{x},y;\boldsymbol{\zeta})=1+\boldsymbol{x}^{\top}\boldsymbol{\zeta}$,
one could get more efficient estimators through our strategy. That means, our
work actually implies a more flexible and powerful framework for nuisance
model calibration in a broader set of missing data problems like SS inference.

Now we shall introduce several corollaries of Section 3.2 to characterize the
efficiency of the DEFUSE estimators. We define
$\mbox{RE}(\hat{b}_{1}|\hat{b}_{2})={\rm AVar}(\hat{b}_{2})/{\rm
AVar}(\hat{b}_{1})$ as the relative efficiency between two asymptotically
normal estimators $\hat{b}_{1}$ and $\hat{b}_{2}$, where ${\rm AVar}(\cdot)$
represents their asymptotic variances. First, we focus on the single BM
scenario and compare the asymptotic variances of the DEFUSE1 and
$\mathcal{L}\mathcal{C}$-only estimators in Corollary 1. We show that their
$\mbox{RE}(\hat{\beta}_{1}|\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}})$
is ensured to be larger than $1$ and increase along with the ratio
$\bar{\rho}=\lim N_{\mathcal{L}\mathcal{M}}/N_{\mathcal{L}\mathcal{C}}$. This
is a natural result as $w(\boldsymbol{x},y;\boldsymbol{\delta}_{0})\equiv 0$
will reduce our construction of DEFUSE1 to the $\mathcal{L}\mathcal{C}$-only
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$ and utilizing more
$\mathcal{L}\mathcal{M}$ samples can amplify our efficiency gain over
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$ not depending on
$\mathcal{L}\mathcal{M}$.

######  Corollary 1 (DEFUSE1 v.s. $\mathcal{L}\mathcal{C}$-only).

Under Assumptions 1–3, 5 and 6, there exists some $\boldsymbol{\delta}_{0}$
such that $w(\boldsymbol{x},y;\boldsymbol{\delta}_{0})\equiv 0$, and we have

$\mbox{RE}(\hat{\beta}_{1}|\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}})=\frac{V_{1}(\boldsymbol{\delta}_{0})}{\min_{\boldsymbol{\delta}}V_{1}({\boldsymbol{\delta}})}\geq
1.$

When $\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)$ and
$\bar{\varphi}_{1}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\boldsymbol{\delta},\bar{\boldsymbol{\eta}})$
have a non-zero correlation,
$\mbox{RE}(\hat{\beta}_{1}|\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}})>1$
and
$\mbox{RE}(\hat{\beta}_{1}|\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}})$
monotonically increases with $\bar{\rho}$ (the sample size ratio of
$\mathcal{L}\mathcal{M}$ to $\mathcal{L}\mathcal{C}$) with all other
parameters unchanged.

Then, still in the single BM scenario, we connect DEFUSE1 with the classic
semiparametric efficient estimator $\hat{\beta}_{\sf
semieff}:=\boldsymbol{c}^{\top}\tilde{\boldsymbol{{\gamma}}}-\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}[\tilde{\varphi}_{1}]+\widehat{{\rm
E}}_{\mathcal{L}\mathcal{M}}[\tilde{\varphi}_{1}]$ obtained by replacing the
truth $\varphi^{*}$ with its estimate $\tilde{\varphi}_{1}$ in (4). We show
that DEFUSE1 will be as efficient as $\hat{\beta}_{\sf semieff}$ in the ideal
case that $\tilde{\varphi}_{1}$ is consistent with $\varphi^{*}$. More
importantly, $\hat{\beta}_{1}$ is guaranteed to be no worse and typically more
efficient than $\hat{\beta}_{\sf semieff}$ when the conditional model for
$\boldsymbol{X}_{\mathcal{P}_{m}}$ is wrongly specified or erroneous and
$\tilde{\varphi}_{1}$ does not converge to the truth. Similar efficiency
properties can be justified for our step of leveraging
$\mathcal{U}\mathcal{C}$ to construct $\hat{\beta}_{2}$; see Appendix.

######  Corollary 2.

Under Assumptions 1–3, 5 and 6, and no matter
$\tilde{\mu}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$ is consistent with the truth
or not, we have

$\mbox{RE}(\hat{\beta}_{1}|\hat{\beta}_{\sf
semieff})=\frac{V_{1}(\boldsymbol{0})}{\min_{\boldsymbol{\delta}}V_{1}({\boldsymbol{\delta}})}\geq
1$

When
${\mu}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\bar{\boldsymbol{\eta}})={\rm
P}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)$,
we have $\mbox{RE}(\hat{\beta}_{1}|\hat{\beta}_{\sf semieff})=1$ and

${\rm AVar}(\hat{\beta}_{1})={\rm
Var}\left[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)\big{|}\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y\right]+\frac{{\rm
Var}\left\\{{\rm
E}\left[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)\big{|}\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y\right]\right\\}}{1+{\bar{\rho}}}.$

######  Remark 7.

By Corollary 2, DEFUSE1 reaches the semiparametric efficiency bound (Robins et
al.,, 1994) when $\tilde{\mu}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$ is
consistent with the truth. This desirable property was not readily achieved by
recent BM fusion approaches relying on other fusing strategies (Xue and Qu,,
2021; Zhao et al.,, 2023; Song et al.,, 2024, e.g.). This implies that with
correct nuisance models, ${\rm AVar}(\hat{\beta}_{1})$ is ensured to be
smaller those of the existing methods. Our numerical results in next sections
support this point. Nevertheless, our method requires the observation of
$\mathcal{L}\mathcal{C}$ samples while methods like Xue and Qu, (2021)
consider a more general case not necessarily having the
$\mathcal{L}\mathcal{C}$ data. Also, our method tends to have higher
computational costs due to MC sampling.

In Corollary 3, we quantify the efficiency gain of incorporating the
$\mathcal{U}\mathcal{C}$ sample to further enhance $\widehat{\beta}_{1}$.
Based on this, all above corollaries and remarks on $\widehat{\beta}_{1}$ will
naturally hold for our final DEFUSE estimator $\hat{\beta}_{2}$. We also
present the asymptotic variance of $\hat{\beta}_{2}$ in the ideal case that
both the conditional models $\tilde{\mu}$ and $\tilde{\nu}$ are consistent
with the true distributions.

######  Corollary 3.

Under Assumptions 1–3, 5 and 6, there exists some $\boldsymbol{\zeta}_{0}$
such that $h(\boldsymbol{x},y;\boldsymbol{\zeta}_{0})\equiv 0$ and

$\mbox{RE}(\hat{\beta}_{2}|\hat{\beta}_{1})=\frac{V_{2}(\boldsymbol{\zeta}_{0})}{\min_{\boldsymbol{\zeta}}V_{2}({\boldsymbol{\zeta}})}\geq
1,$

with a strict $>$ holding when the correlation of
$\bar{S}_{1}(\boldsymbol{X},Y)$ and
$\bar{\varphi}_{2}^{\prime}(\boldsymbol{X})$ is not zero. When both
${\mu}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y;\bar{\boldsymbol{\eta}})={\rm
P}(\boldsymbol{X}_{\mathcal{P}_{m}}|\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y)$
and $\nu(Y|\boldsymbol{X};\bar{\boldsymbol{\theta}})={\rm
P}(Y|\boldsymbol{X})$, we further have

${\rm AVar}(\hat{\beta}_{2})={\rm
Var}\left\\{\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)-{\rm
E}\left[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)\big{|}\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y\right]\big{|}\boldsymbol{X}\right\\}+\frac{{\rm
Var}\left\\{{\rm
E}\left[\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)\big{|}\boldsymbol{X}_{\mathcal{P}_{m^{c}}},Y\right]\right\\}}{1+{\bar{\rho}}}.$

In this case, $\mbox{RE}(\hat{\beta}_{2}|\hat{\beta}_{1})$ will monotonically
increase with $\bar{\rho}$ and we will have ${\rm
AVar}(\hat{\beta}_{2})<\sigma_{\sf ss}^{2}:={\rm
Var}\left\\{\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)\big{|}\boldsymbol{X}\right\\}$,
where $\sigma_{\sf ss}^{2}$ is the semiparametric efficiency bound for the
estimation of $\boldsymbol{c}^{\top}\bar{\boldsymbol{\gamma}}$ under the SS
setting with only $\mathcal{L}\mathcal{C}$ and $\mathcal{U}\mathcal{C}$.

######  Remark 8.

Recent work including Xue et al., (2021) and Song et al., (2024) also
considered the SS and BM setting and only leveraged the large
$\mathcal{U}\mathcal{C}$ sample to assist debiased inference of high-
dimensional linear models. The efficiency gain of the SS estimator
$\hat{\beta}_{2}$ over $\hat{\beta}_{1}$ stated in Corollary 3 has not been
readily achieved in earlier work on this track. For the standard SS setting,
methods proposed by Chakrabortty and Cai, (2018) and others have achieved the
semiparametric efficiency bound $\sigma_{\sf ss}^{2}$ when
$\nu(Y|\boldsymbol{X};\bar{\boldsymbol{\theta}})$ is correct. We show in
Corollary 3 that ${\rm AVar}(\hat{\beta}_{2})<\sigma_{\sf ss}^{2}$ when both
$\tilde{\mu}$ and $\tilde{\nu}$ are consistent with the truths. Meanwhile, we
demonstrate in our simulation and real-world studies that such an efficiency
gain of DEFUSE over existing SS approaches (Gronsbell et al.,, 2022; Miao et
al.,, 2023, e.g.) is generally significant, no matter the nuisance models are
correct or not.

At last, we introduce the efficiency property of $\hat{\beta}_{1}^{\bullet}$
in the multiple BM setting. In Corollary 4, we establish the optimality of the
linear allocation weights $\tilde{\boldsymbol{a}}$ for
$\mathcal{L}\mathcal{M}_{r}$’s, and show that $\hat{\beta}_{1}^{\bullet}$ is
more efficient than the $\mathcal{L}\mathcal{C}$-only
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}$. Although the classic
semiparametric theory is only applied to the single BM scenario (see Section 6
for more discussion), we still find in numerical studies that our way of
constructing the control functions
$\tilde{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)$
for each $\mathcal{L}\mathcal{M}_{r}$ as well as the linear allocation
strategy on $\mathcal{L}\mathcal{M}_{1},\ldots,\mathcal{L}\mathcal{M}_{R}$ can
help to improve the performance $\hat{\beta}_{1}^{\bullet}$ over existing
multiple BM methods not relying on the semiparametric construction; see
Section 4.2 for details.

######  Corollary 4.

Under Assumption 4 and that each $\bar{\rho}_{r}<C$,
$\hat{\beta}_{1}^{\bullet}$ achieves the smallest asymptotic variance among
the linear allocation estimators formed as

$\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}}+\sum_{r=1}^{R}a_{r}\left\\{\widehat{{\rm
E}}_{\mathcal{L}\mathcal{M}_{r}}\left[\tilde{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)\right]-\widehat{{\rm
E}}_{\mathcal{L}\mathcal{C}}\left[\tilde{\varphi}_{1,r}^{\prime}(\boldsymbol{X}_{\mathcal{P}_{m^{c}}^{r}},Y)\right]\right\\},\quad
a_{1},\ldots,a_{R}\in\mathbb{R}.$

Also,
$\mbox{RE}(\hat{\beta}_{1}^{\bullet}|\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}})\geq
1$, where $>$ holds when at least one
$\bar{{\varphi}}_{1,r}(\boldsymbol{X},Y)$ among $r=1,\ldots,R$ have non-zero
correlations with
$\boldsymbol{c}^{\top}\bar{\boldsymbol{S}}(\boldsymbol{X},Y)$.

##  4 Simulation Study

###  4.1 Settings and Benchmarks

In this section, we conduct numerical studies to assess the validity,
efficiency, and robustness of DEFUSE and compare it with a broad set of
existing methods. Our simulation studies include various designs and data
generation settings introduced as follows. In every scenario, we consider
$\boldsymbol{X}=(X_{1},X_{2},X_{3},X_{4},X_{5})^{\top}$ with
$(X_{1},X_{2},X_{3})$ fully observed (in $\mathcal{P}_{m^{c}}$) and
$(X_{4},X_{5})$ possibly missing (in $\mathcal{P}_{m}$).

#### (I) Binary $Y$ and mixture linear ${\rm
P}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$.

In this setting, we first generate binary $Y$ from the bernoulli distribution
with probability $0.5$, and then generate
$\boldsymbol{X}_{\mathcal{P}_{m^{c}}}$ from the following mixture Gaussian:

$\boldsymbol{X}_{\mathcal{P}_{m^{c}}}\mid
Y\sim\begin{cases}N\\{(0.4,-0.5,0.5),\mathbb{I}_{3}\\},&\mbox{ if }Y=1\\\
N\\{(0.6,-0.3,0.7),\mathbb{I}_{3}\\},&\mbox{ if }Y=0,\end{cases}$

where $\mathbb{I}_{3}$ represents the identity matrix. With
$\boldsymbol{X}_{\mathcal{P}_{m^{c}}}$ and $Y$, generate

$\boldsymbol{X}_{\mathcal{P}_{m}}=\begin{cases}\omega\cdot\begin{pmatrix}1&-2&3\\\
2&-2&-1\end{pmatrix}\boldsymbol{X}_{\mathcal{P}_{m^{c}}}+\epsilon&\mbox{ if
}Y=1\\\ \omega\cdot\begin{pmatrix}-1&3&-2\\\
-3&2&2\end{pmatrix}\boldsymbol{X}_{\mathcal{P}_{m^{c}}}+\epsilon&\mbox{ if
}Y=0\\\ \end{cases}$

where $\epsilon\sim{\rm N}(0,\mathbb{I}_{2})$ is an independent noise term.
The hyperparameter $\omega$ is set to range from $0.1$ to $1$ characterizing
the dependence of $\boldsymbol{X}_{\mathcal{P}_{m}}$ on the other variables.
The sample sizes of the three data sets are set to be
$N_{\mathcal{L}\mathcal{C}}=500$, $N_{\mathcal{U}\mathcal{C}}=30000$, and
$N_{\mathcal{L}\mathcal{M}}=\rho\cdot 500$ with the hyperparameter
$\rho\in\\{2,3,4,5\\}$.

#### (II) Binary $Y$ and nonlinear ${\rm
P}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$.

We generate $(Y,\boldsymbol{X}_{\mathcal{P}_{m^{c}}})$ in the same way as in
(I) but generate $\boldsymbol{X}_{\mathcal{P}_{m}}$ through a nonlinear model:
$\boldsymbol{X}_{\mathcal{P}_{m^{c}}}$ and $Y$, generate

$\boldsymbol{X}_{\mathcal{P}_{m}}=\begin{cases}\begin{pmatrix}0.3&-0.6&0.9\\\
0.6&-0.6&-0.3\end{pmatrix}(X_{1}^{2},X_{2}^{2},X_{3}^{2})^{\top}+\begin{pmatrix}-0.5\\\
0.5\end{pmatrix}+\epsilon&\mbox{ if }Y=1\\\ \begin{pmatrix}-0.3&0.9&-0.6\\\
-0.9&0.6&0.6\end{pmatrix}(X_{1}^{2},X_{2}^{2},X_{3}^{2})^{\top}+\epsilon&\mbox{
if }Y=0\\\ \end{cases}$

where $\epsilon\sim{\rm N}(0,\mathbb{I}_{2})$ is an independent noise. Set
$N_{\mathcal{L}\mathcal{C}}=500$, $N_{\mathcal{U}\mathcal{C}}=30000$, and
$N_{\mathcal{L}\mathcal{M}}=2500$.

#### (III) Continuous $Y$ and linear model of ${\rm
P}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$.

First, we generate $(X_{1},X_{2},X_{3})$ independently from the uniform
distribution on $[-1,1]$, and then generate the rest of $\boldsymbol{X}$ by

$X_{\ell}=(0.5,0.5,0.5)\boldsymbol{X}_{\mathcal{P}_{m^{c}}}+\epsilon_{\ell},\quad\epsilon_{\ell}\sim
N(0,1),\quad\mbox{for }\ell=4,5.$

The following linear model is used to simulate the response $Y$:

$\displaystyle Y=-X_{1}-X_{2}+X_{3}+\omega X_{4}+\omega
X_{5}+\epsilon_{Y},\quad\epsilon_{Y}\sim N(0,1).$

The hyperparameter $\omega$ is set to range from $0.1$ to $1$ characterizing
the association between $\boldsymbol{X}_{\mathcal{P}_{m}}$ and $Y$ conditional
on $\boldsymbol{X}_{\mathcal{P}_{m}^{c}}$. The sample sizes of the three data
sets are set to be
$N_{\mathcal{L}\mathcal{C}}=500,N_{\mathcal{U}\mathcal{C}}=30000$, and
$N_{\mathcal{L}\mathcal{M}}=\rho\cdot 500$ with the hyperparameter
$\rho\in\\{2,3,4,5\\}$.

#### (IV) Linear regression with multiple BM scenario.

First, we generate $(Y,\boldsymbol{X})$ in the same way as in (III). In this
case, set $R=2$, $\boldsymbol{X}_{\mathcal{P}_{m}}^{1}=X_{4}$, and
$\boldsymbol{X}_{\mathcal{P}_{m}}^{2}=X_{5}$. The sample sizes of the four
data sets are set to be
$N_{\mathcal{L}\mathcal{C}}=N_{\mathcal{L}\mathcal{M}_{1}}=N_{\mathcal{L}\mathcal{M}_{2}}=500$
and $N_{\mathcal{U}\mathcal{C}}=30000$.

In the settings (I) and (II), the outcome model is the logistic model for
$Y\sim\boldsymbol{X}$. We also include (III) and (IV) with linear models for
$Y$ to facilitate the comparison with a broader set of existing methods for BM
regression. In (I) and (III), the hyperparameters $\omega$ and $\rho$ are
introduced and varied to study how the efficiency gain of our estimators
changes according to the dependence between the missing and the observed
variables as well as to the sample size ratio between $\mathcal{L}\mathcal{M}$
and $\mathcal{L}\mathcal{C}$. In (II), the hyperparameters $\omega$ and $\rho$
are fixed, and the dependence of missing conditioning on the observed variable
is introduced through the quadratic function. In (IV), the hyperparameters
$\omega$ and $\rho$ are fixed, and the more complex missing structure is
brought into the simulation to appraise the robustness and efficiency of the
estimators. In each scenario, we evaluate and compare the accuracy and
efficiency of the estimators through their bias and mean square errors on each
coefficient or averaged over all coefficients.

#### Methods under comparison.

We summarize all methods under comparison with some implementation details.
$\mathcal{L}\mathcal{C}$-only: The standard estimator
$\tilde{\boldsymbol{\gamma}}$ obtained using $\mathcal{L}\mathcal{C}$ data
only. DEFUSE1: Our proposed estimator $\hat{\beta}_{1}$ based on
$\mathcal{L}\mathcal{C}$ and $\mathcal{L}\mathcal{M}$ data. DEFUSE: Our final
estimator $\hat{\beta}_{2}$. MICE: Multiple imputation by chained equations
with the default settings and tuning parameters of the R package _mice_ (van
Buuren and Groothuis-Oudshoorn,, 2011). MBI: The multiple blockwise imputation
method proposed by Xue and Qu, (2021) combining $\mathcal{L}\mathcal{C}$ and
$\mathcal{L}\mathcal{M}$ data. HTLGMM: The data fusion approach proposed by
Zhao et al., (2023) based on the generalized method of moments. SSB: The SS
and BM fusion method proposed by Song et al., (2024) using the
$\mathcal{L}\mathcal{C}$, $\mathcal{L}\mathcal{M}$ and
$\mathcal{U}\mathcal{C}$ samples. Benchmarks MBI, HTLGMM, and SSB are all
implemented under our low-dimensional regression settings without using any
sparse penalties. SemiEff: The classic semiparametric efficient estimator
$\hat{\beta}_{\sf semieff}$ constructed following Robins et al., (1994). SSL:
The SS estimator using $\mathcal{L}\mathcal{C}$ and $\mathcal{U}\mathcal{C}$
(Gronsbell et al.,, 2022, e.g.), with its imputation model for $Y$ constructed
through natural spline regression. We summarize all these methods in Table 1.

Method | Nonlinear Model | Single BM | Multiple BM | Leverage $\mathcal{U}\mathcal{C}$
---|---|---|---|---
DEFUSE1 | $\checkmark$ | $\checkmark$ | $\checkmark$ | –
DEFUSE | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\checkmark$
MICE | $\checkmark$ | $\checkmark$ | $\checkmark$ | –
MBI | – | $\checkmark$ | $\checkmark$ | –
SSL | $\checkmark$ | – | – | $\checkmark$
HTLGMM | $\checkmark$ | $\checkmark$ | – | –
SSB | – | $\checkmark$ | $\checkmark$ | $\checkmark$
$\mathcal{L}\mathcal{C}$-only | $\checkmark$ | – | – | –
SemiEff | $\checkmark$ | $\checkmark$ | – | –
Table 1:  A list of the methods under comparison in our simulation studies. We
also make a summary on their applicability to nonlinear (logistic) GLMs,
whether they use the single or multiple BM data sets, and whether they use the
unlabeled complete sample $\mathcal{U}\mathcal{C}$. All methods can be
implemented and evaluated in our linear model and single BM setting (III) and
some of them are only applicable to a subset of our simulation settings, e.g.,
MBI only applied to (III) and (IV) targeting the linear model.

For the conditional distribution of $\boldsymbol{X}_{\mathcal{P}_{m}}$ used in
multiple approaches, including ours, we separately consider using both the
parametric regression and the nonparametric kernel machine method to learn its
conditional mean function and the standard method of moment to estimate its
conditional covariance matrix as described in Remark 1. For ease of
presentation, we will mainly report the results obtained using the parametric
strategy and the results of kernel machine are provided in the Appendix. For
parametric regression, we fit the ordinary multivariate linear regression for
the missing covariates. In (I) and (II), we fit Gaussian linear models for
$\boldsymbol{X}_{\mathcal{P}_{m}}\sim\boldsymbol{X}_{\mathcal{P}_{m}^{c}}$
stratified by the binary $Y$. In (III) and (IV), we fit Gaussian linear
regression for
$\boldsymbol{X}_{\mathcal{P}_{m}}\sim(\boldsymbol{X}_{\mathcal{P}_{m}^{c}},Y)$.
In settings (I), (III), and (IV), our parametric model for
$\boldsymbol{X}_{\mathcal{P}_{m}^{c}}$ is correctly specified while it is
wrong under (II) due to the quadratic effects of
$\boldsymbol{X}_{\mathcal{P}_{m}^{c}}$. We use the parametric modeling
strategy for the conditional distribution of $Y$ used for imputation in our
method. For the data-adaptive calibration functions, we set
$w(\boldsymbol{x},y;\boldsymbol{\delta})=1+\boldsymbol{\delta}^{\top}(1,\boldsymbol{x}_{\mathcal{P}_{m}^{c}},y)$
and $h(\boldsymbol{x},y;\boldsymbol{\zeta})=1+\zeta_{1}$ in all settings.

###  4.2 Comparison with existing methods

In settings (I) and (III), we evaluate the accuracy and efficiency of the
included estimators with a fixed $\omega=0.3$ in (I) and $\omega=1$ in (III),
and the sample size ratio $\rho\in\\{2,3,4,5\\}$. The results are summarized
in Table 2 for setting (I) and Table 3 for (III). These two tables present the
relative efficiency (RE) to the $\mathcal{L}\mathcal{C}$-only estimator
defined as the inverse ratio of the average mean squared errors (MSE) over the
coefficients $\gamma_{1},\ldots,\gamma_{5}$. They also include the average
mean absolute bias overall coefficients.

All methods except MICE produce small bias and valid estimation. Compared to
recent BM data fusion methods, DEFUSE1 attains significantly higher efficiency
than HTLGMM in the logistic model setting (I) and higher efficiency than MBI,
SSB, and HTLGMM in the linear setting (III). For example, the RE of DEFUSE1 is
around $30\%$ larger than that of MBI in setting (III) with $\rho=2$. This
improvement becomes even larger when
$\rho=N_{\mathcal{L}\mathcal{M}}/N_{\mathcal{L}\mathcal{C}}$ grows. As another
example, DEFUSE1 is $40\%$ more efficient than HTLGMM in the logistic setting
(I) with $\rho=3$. As seen from Table 5, our method also attains better
performance over both MBI and SSB in setting (IV) with the linear outcome
model and multiple BM data structure. For instance, DEFUSE1 has around $50\%$
higher RE than MBI and SSB on $\gamma_{1}$ in setting (IV).

Meanwhile, DEFUSE1 shows close REs to the classic semiparametric efficient
estimator (SemiEff) in both settings (I) and (III) with a correct model
$\mu_{\boldsymbol{X}_{\mathcal{P}_{m}}}$ for
$\boldsymbol{X}_{\mathcal{P}_{m}}$. This finite-sample result is consistent
with our Corollary 2 and Remark 7 that when
$\tilde{\mu}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$ is consistent with the truth,
DEFUSE1 will be asymptotically equivalent to SemiEff. Moreover, through
incorporating the large $\mathcal{U}\mathcal{C}$ sample, DEFUSE achieves
around $10\%$ improved average efficiency over DEFUSE1 in all settings. This
agrees with our Corollary 3. Such improvements will be demonstrated and
discussed more comprehensively and carefully in Section 4.3. We also notice
that the standard SSL method has no efficiency gain over the
$\mathcal{L}\mathcal{C}$-only estimator in settings (I) and (III) because it
only incorporates $\mathcal{U}\mathcal{C}$, which has been shown to contribute
zero improvement when the outcome model for $Y$ is correctly specified
(Chakrabortty and Cai,, 2018). This finding also corresponds to Remark 3 that
reversing our Steps 2 and 3 in Algorithm 1 tends to result in less efficient
estimator.

In Table 4, we present the RE (to $\mathcal{L}\mathcal{C}$-only) on each
coefficient in the simulation setting (II) with a nonlinear model for
$\boldsymbol{X}_{\mathcal{P}_{m}}$ not captured by the parametric conditional
model. To demonstrate the importance of properly specifying the calibration
function $w(\boldsymbol{x},y;\boldsymbol{\delta})$, we take
$w(\boldsymbol{x},y;\boldsymbol{\delta})=1+\boldsymbol{\delta}^{\top}\xi(\boldsymbol{x},y)$
and include two choices on $\xi(\boldsymbol{x},y)$ for the construction of
DEFUSE1: (1) $\xi(\boldsymbol{x},y)=1$ and (2)
$\xi(\boldsymbol{x},y)=(1,\boldsymbol{x}_{\mathcal{P}_{m}^{c}},y)$. Choice
$\xi(\boldsymbol{x},y)=1$ only involves an “intercept” term for calibration
and is actually relevant to proposals in recent work like Miao et al., (2023).
Interestingly, DEFUSE1 with
$\xi(\boldsymbol{x},y)=(1,\boldsymbol{x}_{\mathcal{P}_{m}^{c}},y)$ outperforms
both SemiEff and DEFUSE1 with $\xi(\boldsymbol{x},y)=1$. For example, on
$\gamma_{3}$, DEFUSE1 with
$\xi(\boldsymbol{x},y)=(1,\boldsymbol{x}_{\mathcal{P}_{m}^{c}},y)$ has around
$50\%$ higher RE than both SemiEff and DEFUSE1 with $\xi(\boldsymbol{x},y)=1$.
This demonstrates that our delicate data-adaptive calibration approach can
effectively improve the estimation efficiency over the classic SemiEff when
the imputation model $\tilde{\mu}_{\boldsymbol{X}_{\mathcal{P}_{m}}}$ is
wrongly specified or of poor quality. To realize this advantage, it is
important to specify the calibration function
$w(\boldsymbol{x},y;\boldsymbol{\delta})$ in a proper way, which can be seen
by comparing DEFUSE1 under the two different choices on
$\xi(\boldsymbol{x},y)$ in this setting.

Methods | ratio $\rho$ = 2 | ratio $\rho$ = 3 | ratio $\rho$ = 4 | ratio $\rho$ = 5
---|---|---|---|---
RE | Bias | RE | Bias | RE | Bias | RE | Bias
DEFUSE1 | 1.74 | 0.011 | 1.90 | 0.012 | 2.15 | 0.015 | 2.23 | 0.015
DEFUSE | 1.88 | 0.007 | 2.09 | 0.010 | 2.45 | 0.012 | 2.56 | 0.012
MICE | 0.93 | 0.190 | 1.00 | 0.214 | 1.18 | 0.229 | 1.30 | 0.240
SSL | 1.01 | 0.020 | 1.01 | 0.020 | 1.01 | 0.020 | 1.01 | 0.020
HTLGMM | 1.42 | 0.029 | 1.46 | 0.033 | 1.52 | 0.035 | 1.54 | 0.036
$\mathcal{L}\mathcal{C}$-only | 1 | 0.019 | 1 | 0.019 | 1 | 0.019 | 1 | 0.019
SemiEff | 1.76 | 0.010 | 1.92 | 0.010 | 2.18 | 0.012 | 2.26 | 0.012
Table 2:  Average absolute bias and relative efficiency (RE) to $\mathcal{L}\mathcal{C}$-only of the methods in simulation setting (I) with $\omega=0.3$ and $\rho\in\\{2,3,4,5\\}$. Methods | ratio $\rho$ = 2 | ratio $\rho$ = 3 | ratio $\rho$ = 4 | ratio $\rho$ = 5
---|---|---|---|---
RE | Bias | RE | Bias | RE | Bias | RE | Bias
DEFUSE1 | 1.32 | 0.006 | 1.38 | 0.007 | 1.39 | 0.007 | 1.42 | 0.007
DEFUSE | 1.41 | 0.004 | 1.51 | 0.004 | 1.56 | 0.004 | 1.62 | 0.004
MICE | 0.06 | 0.294 | 0.04 | 0.381 | 0.03 | 0.453 | 0.02 | 0.506
MBI | 1.11 | 0.006 | 1.15 | 0.007 | 1.14 | 0.007 | 1.16 | 0.007
SSL | 0.98 | 0.005 | 0.98 | 0.005 | 0.98 | 0.005 | 0.98 | 0.005
HTLGMM | 1.15 | 0.005 | 1.21 | 0.005 | 1.20 | 0.006 | 1.22 | 0.006
SSB | 0.95 | 0.007 | 1.03 | 0.009 | 1.07 | 0.006 | 1.11 | 0.006
$\mathcal{L}\mathcal{C}$-only | 1 | 0.005 | 1 | 0.005 | 1 | 0.005 | 1 | 0.005
SemiEff | 1.32 | 0.006 | 1.38 | 0.007 | 1.39 | 0.007 | 1.42 | 0.007
Table 3:  Average absolute bias and relative efficiency (RE) to $\mathcal{L}\mathcal{C}$-only in simulation setting (III) with $\omega=1$ and $\rho\in\\{2,3,4,5\\}$. Method | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$ | $\gamma_{4}$ | $\gamma_{5}$
---|---|---|---|---|---
DEFUSE1 $\xi(\boldsymbol{x},y)=1$ | 3.48 | 3.24 | 2.62 | 0.99 | 0.98
DEFUSE1 $\xi(\boldsymbol{x},y)=(1,\boldsymbol{x}_{\mathcal{P}_{m}^{c}},y)$ | 4.35 | 4.35 | 3.59 | 1.21 | 1.11
SemiEff | 3.53 | 3.13 | 2.35 | 0.98 | 0.99
SSL | 1.00 | 1.10 | 0.99 | 0.93 | 1.38
HTLGMM | 2.31 | 2.64 | 1.64 | 1.10 | 1.13
Table 4: Relative efficiency (RE) to $\mathcal{L}\mathcal{C}$-only on each coefficient in $\boldsymbol{\gamma}$ under Setting (II). Method | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$ | $\gamma_{4}$ | $\gamma_{5}$
---|---|---|---|---|---
DEFUSE1 | 1.60 | 1.78 | 2.39 | 1.49 | 1.67
DEFUSE | 1.69 | 1.83 | 2.55 | 1.51 | 1.68
MICE | 0.52 | 0.63 | 0.69 | 0.20 | 0.22
MBI | 1.06 | 1.09 | 1.06 | 0.95 | 1.01
SSB | 1.07 | 1.12 | 1.51 | 0.60 | 0.81
SSL | 1.01 | 1.00 | 1.00 | 0.99 | 1.02
Table 5: Relative efficiency (RE) to $\mathcal{L}\mathcal{C}$-only on each
coefficient in $\boldsymbol{\gamma}$ under Setting (IV).

###  4.3 Patterns of relative efficiency

In this section, we present more detailed results in settings (I) and (III) to
assess DEFUSE under different values of
$\rho=N_{\mathcal{L}\mathcal{M}}/N_{\mathcal{L}\mathcal{C}}$ and $\omega$
quantifying the dependence of $\boldsymbol{X}_{\mathcal{P}_{m}}$ on
$(Y,\boldsymbol{X}_{\mathcal{P}_{m}^{c}})$. We only present results on the
coefficients $\gamma_{1,3,5}$. The results on other coefficients are similar
and can be found in Appendix. We focus on the RE of DEFUSE1 to the
$\mathcal{L}\mathcal{C}$-only estimator as well as the RE of DEFUSE to
DEFUSE1.

In Figures 3 and 5, we present the REs of DEFUSE1 to
$\mathcal{L}\mathcal{C}$-only. A larger efficiency gain of DEFUSE1 over
$\mathcal{L}\mathcal{C}$-only can be observed with an increasing sample size
ratio $\rho=N_{\mathcal{L}\mathcal{M}}/N_{\mathcal{L}\mathcal{C}}$; see also
Tables 2 and 3. This finding agrees with our Corollary 1 that
$\mbox{RE}(\hat{\beta}_{1}|\boldsymbol{c}^{\top}\tilde{\boldsymbol{\gamma}})$
monotonically increases with $\rho$. On $\gamma_{5}$ corresponding to the
partially missing $X_{5}$, the efficiency gain of DEFUSE1 tends to be higher
with an increasing dependence parameter $\omega$. Intuitively, when
$\omega=0$, $X_{5}$ will be independent from
$(Y,\boldsymbol{X}_{\mathcal{P}_{m}^{c}})$ and the $\mathcal{L}\mathcal{M}$
sample with only $(Y,\boldsymbol{X}_{\mathcal{P}_{m}^{c}})$ observed cannot
provide any information on $X_{5}$. Consequently, incorporating
$\mathcal{L}\mathcal{M}$ could not help much on $\gamma_{5}$ whose score
function is largely dependent on $X_{5}$. On the other hand, under a large
$\omega$, $(Y,\boldsymbol{X}_{\mathcal{P}_{m}^{c}})$ observed on
$\mathcal{L}\mathcal{M}$ can be informative to the score function of
$\gamma_{5}$ and improve the estimation efficiency on it. Interestingly, for
$\gamma_{1}$ and $\gamma_{3}$ corresponding to covariates in
$\boldsymbol{X}_{\mathcal{P}_{m}^{c}}$ observed on $\mathcal{L}\mathcal{M}$,
increasing $\omega$ tends to cause a lower efficiency gain of our method,
which is opposite to the pattern of REs on $\gamma_{5}$. For example, in
setting (III) with $\rho=1$, the RE of DEFUSE1 on $\gamma_{3}$ decrease from
$4$ to $2$ when $\omega$ changes from $0.1$ to $0.8$; see Figure 5. This is
because when $\omega$ becomes larger, the score functions of $\gamma_{1,3}$
tend to be more dependent on $\boldsymbol{X}_{\mathcal{P}_{m}}$, which makes
the missingness of $\boldsymbol{X}_{\mathcal{P}_{m}}$ on
$\mathcal{L}\mathcal{M}$ more impactful.

In Figures 4 and 6, we explore the RE trends of DEFUSE to DEFUSE1 under
different sample size ratios $\rho$ and weights $\omega$. Similarly, we find
that a larger $\rho$ results in more efficiency gain of DEFUSE over DEFUSE1,
which agrees with Corollary 3 that
$\mbox{RE}(\hat{\beta}_{1}|\hat{\beta}_{2})$ will increase with $\rho$ when
the nuisance $\tilde{\mu}$ and $\tilde{\nu}$ are of good quality. Also, on
$\gamma_{5}$ corresponding to the BM covaraite $X_{5}$, the RE of DEFUSE to
DEFUSE1 tends to become larger when the dependence of $X_{5}$ on
$(Y,\boldsymbol{X}_{\mathcal{P}_{m}^{c}})$ characterized by $\omega$ gets
larger. Nevertheless, the overall improvements of DEFUSE over DEFUSE1 is
relatively small on $\gamma_{5}$. Interestingly, on $\gamma_{1}$ and
$\gamma_{3}$, we find that the RE first increases and then decreases with
$\omega$ being larger. Future theoretical analysis is needed to help
understanding and explaining this non-monotonic trend better.

Figure 3: Relative efficiency (RE) of DEFUSE1 to $\mathcal{L}\mathcal{C}$-only
on $\gamma_{1,3,5}$ in setting (I), with $\omega\in\\{0.1,0.2,\ldots,1\\}$ and
$\rho\in\\{2,3,4,5\\}$. The horizontal black line represents $\mbox{RE}=1$.
Figure 4: Relative efficiency (RE) of DEFUSE to DEFUSE1 on $\gamma_{1,3,5}$ in
setting (I), with $\omega\in\\{0.1,0.2,\ldots,1\\}$ and
$\rho\in\\{2,3,4,5\\}$. The horizontal black line represents $\mbox{RE}=1$.
Figure 5: Relative efficiency (RE) of DEFUSE1 to $\mathcal{L}\mathcal{C}$-only
on $\gamma_{1,3,5}$ in setting (III), with $\omega\in\\{0.1,0.2,\ldots,1\\}$
and $\rho\in\\{2,3,4,5\\}$. The horizontal black line represents
$\mbox{RE}=1$. Figure 6: Relative efficiency (RE) of DEFUSE to DEFUSE1 on
$\gamma_{1,3,5}$ in setting (III), with $\omega\in\\{0.1,0.2,\ldots,1\\}$ and
$\rho\in\\{2,3,4,5\\}$. The horizontal black line represents $\mbox{RE}=1$.

##  5 Heart Disease Risk Modeling with MIMIC-III

To illustrate the utility of DEFUSE, we apply it to MIMIC-III (Medical
Information Mart for Intensive Care III; Johnson et al., (2016)), a large
database comprising de-identified electronic health records (EMRs) of over
forty thousand patients associated with their stay in critical care units of
the Beth Israel Deaconess Medical Center between 2001 and 2012. The database
includes various types of information such as demographics, vital sign
measurements made at the bedside, laboratory test results, procedures,
medications, caregiver notes, imaging reports, and mortality in and out of the
hospital. As EMRs include large and diverse patient populations and cover
comprehensive phenotypes, they have been widely used by researchers to
identify underlying risk factors for certain disease conditions. Our goal is
to make statistical inference on the risk prediction model for heart diseases
(HD) against some potential risk factors related to the health condition of a
subject.

Our response of interest $Y\in\\{0,1\\}$ is the HD status identified from the
manual chart reviewing labels obtained in a previous study on MIMIC-III data
(Gehrmann et al.,, 2018); see their paper for the detailed definition of the
phenotype HD. We greatly thank the authors for making their gold standard
labels accessible. For $\boldsymbol{X}$, we include age, encounter of the
International Classification of Diseases (ICD) codes of low high-density
lipoproteins (HDL), ICD encounter of type II diabetes, and several HD-related
laboratory (lab) variables selected using KESER (Hong et al.,, 2021), an
online tool for extracting clinical concepts and features relevant to a given
phenotype. We exclude lab variables too closely linked with HDL or diabetes
such as A1C. For the convenience of downstream risk modeling, we also derive a
lab risk score (LRS) through weighting and combining the selected lab
variables with the main ICD code of HD as a surrogate outcome on the large
$\mathcal{U}\mathcal{C}$ sample. Since the LRS is derived with large sample,
it has an negligible estimation uncertainty and can be viewed as given. See
the list of selected lab variables and details about the LRS in Appendix.

Noticing that each lab variable is not observed on a fraction of patients, we
separately consider two settings of missingness: (i) a single BM task with one
BM covariate taken as the derived LRS; (ii) a multiple BM task with two BM lab
variables thyroid-stimulating hormone (TSH) and hypochromia (HPO), with joint
missingness on $\mathcal{L}\mathcal{M}_{1}$, TSH-only missingness on
$\mathcal{L}\mathcal{M}_{2}$, HPO-only missing on
$\mathcal{L}\mathcal{M}_{3}$. On task (i), there are
$N_{\mathcal{L}\mathcal{C}}=467$ labeled subject with complete observation of
$\boldsymbol{X}$, $N_{\mathcal{L}\mathcal{M}}=678$ labeled subjects with BM
covariates, and $N_{\mathcal{U}\mathcal{C}}=7037$ unlabeled subjects with
complete $\boldsymbol{X}$. Similar to the simulation, we compare our
approaches to four applicable benchmarks to this case, including
$\mathcal{L}\mathcal{C}$-only, SSL, Semieff, and HTLGMM. We do not include
MICE since it tends to result in excessive bias, as seen from the simulation
results. For (ii), there are $N_{\mathcal{L}\mathcal{C}}=486$ labeled subject
with complete observation of $\boldsymbol{X}$,
$N_{\mathcal{L}\mathcal{M}_{1}}=231$, $N_{\mathcal{L}\mathcal{M}_{2}}=175$,
and $N_{\mathcal{L}\mathcal{M}_{3}}=153$ for three types of missing, and
$N_{\mathcal{U}\mathcal{C}}=7014$. For the implementation of DEFUSE, we follow
the same strategy as our simulation with the conditional models constructed
through kernel machine and the data-adaptive calibration functions
$w(\boldsymbol{x},y;\boldsymbol{\delta})=1+\boldsymbol{\delta}^{\top}(1,\boldsymbol{x}_{\mathcal{P}_{m}^{c}},y)$
and $h(\boldsymbol{x},y;\boldsymbol{\zeta})=1+\zeta_{1}$.

Table 6 provides a summary of the estimated coefficients on the single BM task
(i), as well as the bootstrap standard errors and $p$-values. Table 7 reports
the relative efficiency (RE) of the data-fusing estimators to the standard
$\mathcal{L}\mathcal{C}$-only estimator. DEFUSE achieves the highest RE on all
the coefficients among the included methods. As an example, for the estimation
of the HDL coefficient, our method reduces more than $50\%$ variance compared
to $\mathcal{L}\mathcal{C}$-only and its RE to HTLGMM is around $1.4$.

Through the output of DEFUSE, one could conclude that given the age and lab
measures on other risk factors, a higher risk of HD is associated with having
low HDL at the level $0.1$ as well as type II diabetes at the level $0.05$.
Notably, with the methods $\mathcal{L}\mathcal{C}$-only and SSL, both the
effects of HDL and diabetes are shown to be insignificant even at the $0.1$
level. Also, diabetes turns out to be not significant when using HTLGMM. This
suggests that among all the methods under comparison, DEFUSE is the only one
powerful enough to discover both signals. Interestingly, it is well-documented
in existing biomedical literature on the associations between heart disease
and diabetes (Haffner,, 2000; Peters et al.,, 2014, e.g.) as well as low HDL
(Després et al.,, 2000; Rader and Hovingh,, 2014, e.g.), which supports the
findings of DEFUSE. Table 8 summarizes the RE to the
$\mathcal{L}\mathcal{C}$-only estimator in the multiple BM scenario introduced
above. DEFUSE attains more than $25\%$ efficiency improvement over
$\mathcal{L}\mathcal{C}$-only on all coefficients and more than $50\%$
variance reduction on that of age, which is again more significant than SSL.
One can find the summary table of the estimates and $p$-values in this
scenario in Appendix.

Method | Age | HDL | Diabetes | LRS
---|---|---|---|---
$\text{coe}_{\text{se}}$ |  $p$-value | $\text{coe}_{\text{se}}$ |  $p$-value | $\text{coe}_{\text{se}}$ |  $p$-value | $\text{coe}_{\text{se}}$ |  $p$-value
$\mathcal{L}\mathcal{C}$-only | $0.41_{0.120}$ |  $<$0.001 | $0.35_{0.219}$ | 0.120 | $0.15_{0.215}$ | 0.48 | $0.25_{0.113}$ | 0.028
DEFUSE1 | $0.46_{0.0823}$ |  $<$0.001 | $0.25_{0.148}$ | 0.086 | $0.34_{0.146}$ | 0.020 | $0.22_{0.112}$ | 0.048
DEFUSE | $0.45_{0.0821}$ |  $<$0.001 | $0.25_{0.147}$ | 0.087 | $0.33_{0.145}$ | 0.022 | $0.22_{0.111}$ | 0.051
SSL | $0.41_{0.140}$ | 0.003 | $0.36_{0.240}$ | 0.129 | $0.15_{0.238}$ | 0.532 | $0.25_{0.113}$ | 0.027
HTLGMM | $0.47_{0.102}$ |  $<$0.001 | $0.37_{0.192}$ | 0.053 | $0.20_{0.193}$ | 0.295 | $0.25_{0.117}$ | 0.032
Table 6: Fitted coefficient (coe), empirical standard error (se), and $p$-value for the coefficient of each covariate in the single BM scenario of our real-world study.  Method | Age | HDL | Diabetes | Lab Score
---|---|---|---|---
DEFUSE1 | 2.11 | 2.19 | 2.18 | 1.00
DEFUSE | 2.13 | 2.21 | 2.19 | 1.10
SSL | 1.00 | 1.01 | 1.00 | 1.01
HTLGMM | 1.89 | 1.60 | 1.47 | 1.00
Table 7: Relative efficiency (RE) to the standard $\mathcal{L}\mathcal{C}$-only estimator on each coefficient in the single BM scenario of our real-world study. Method | Age | HDL | Diabetes | TSH | Hypochromia
---|---|---|---|---|---
DEFUSE1 | 2.41 | 1.29 | 1.30 | 1.71 | 1.27
DEFUSE | 2.42 | 1.29 | 1.31 | 1.82 | 1.28
SSL | 1.00 | 0.99 | 1.06 | 0.94 | 1.13
Table 8: Relative efficiency (RE) to the standard
$\mathcal{L}\mathcal{C}$-only estimator on each coefficient in the multiple BM
scenario of our real-world study.

##  6 Discussion

In this paper, we develop DEFUSE, a novel approach for robust and efficient
data-fusion in the presence of blockwise missing covariates and large
unlabeled samples. Its validity and relative efficiency compared to a
comprehensive set of existing methods are justified through theoretical,
numerical, and real-world studies. The principle of our construction strategy
is to leverage the general and powerful theory on semiparametric efficient
estimation while incorporating the key data-adaptive calibration and optimal
linear allocation procedures to maintain robustness and effectiveness to a
more realistic scenario that the nuisance distributional models are imperfect
or problematic. We shall point out the limitation and future direction of our
work.

Our linear allocation strategy for the multiple BM scenario does not achieve
semiparametric efficiency. We notice that the general canonical gradient
method proposed by Li and Luedtke, (2023) could be leveraged to further
improve $\hat{\beta}_{1}^{\bullet}$ in this direction, even with non-nested
sets of observed covariates in
$\mathcal{L}\mathcal{M}_{1},\ldots,\mathcal{L}\mathcal{M}_{R}$. Nevertheless,
it is still an open question about how to incorporate our data-adaptive
calibration approach and the SS setting into this more complicated framework.
In the simulation study, we adopt the parametric and kernel methods to
construct the conditional models $\tilde{\mu}$ and $\tilde{\nu}$. As was
pointed out, DEFUSE allows the use of more complex tools like generative
adversarial models to estimate the unknown conditional distributions, which
could achieve better performance for relatively high-dimensional missing
covariates. However, care must be taken to avoid over-fitting and overly high
computational costs.

Also, our current method requires $\mathcal{L}\mathcal{C}$ sample that may not
be accessible in certain scenarios (Song et al.,, 2024, e.g.). It is not hard
to see that (2) can actually be identified and constructed with the
$\mathcal{U}\mathcal{C}$ and
$\mathcal{L}\mathcal{M}_{1},\ldots,\mathcal{L}\mathcal{M}_{R}$ samples as long
as the union set of the observed covariates in $\mathcal{L}\mathcal{M}_{1}$’s
is complete. It is interesting to generalize our method to accommodate this
setting without observation of $\mathcal{L}\mathcal{C}$. One possible way is
to first utilize a part of $\mathcal{L}\mathcal{M}$ data sets to derive a
preliminary estimator that is asymptotically normal but not efficient, and
then implement similar steps to Algorithm 1 to reduce its variance. Other more
complicated missing structures, e.g., single or multiple missing blocks in the
unlabeled sample may also warrant future research.

We notice a comprehensive set of recent literature in addressing covariate
shift of $\boldsymbol{X}$ (Liu et al.,, 2023; Qiu et al.,, 2023, e.g.), label
or model shift of $Y$ (Li et al., 2023a, ; Zhao et al.,, 2023, e.g.), as well
as high-dimensionality (Xue et al.,, 2021; Song et al.,, 2024, e.g.) in data
fusion under various types of missingness. DEFUSE can also be generalized to
address these practical issues. For high-dimensional sparse regression model
of $Y$, we can use debiased Lasso (Zhang and Zhang,, 2014; Van de Geer et
al.,, 2014; Javanmard and Montanari,, 2014, e.g.) in Step 1 of Algorithm 1 to
derive the preliminary estimator
$\boldsymbol{c}^{\top}\tilde{\boldsymbol{{\gamma}}}$ that is
$n^{-1/2}$-consistent and asymptotically normal. Then the remaining Steps 2
and 3 can be implemented in a similar way as Algorithm 1. Since our main
theorems only require the $o_{\sf P}(1)$ convergence of the nuisance models,
construction procedures of our nuisance parameters could naturally incorporate
sparse regularization to accommodate high-dimensionality without introducing
excessive bias. For potential covariate shift across the
$\mathcal{L}\mathcal{C}$, $\mathcal{L}\mathcal{M}$ and
$\mathcal{U}\mathcal{C}$ samples, one can introduce density ratio (propensity
score) models for adjustment (Qiu et al.,, 2023; Liu et al.,, 2023, e.g.)
based on our current MCAR construction. This can result in doubly robust forms
requiring stronger assumption on the correctness and convergence of the
nuisance models. For the conditional model shift of $Y$, adaptive data fusion
under our setting can be challenging, especially for statistical inference. We
notice some related recent literature like Han et al., (2021) and Li et al.,
2023a , and leave this open problem to future research.

## References

  * Angelopoulos et al., (2023) Angelopoulos, A. N., Duchi, J. C., and Zrnic, T. (2023).  Ppi++: Efficient prediction-powered inference.  arXiv preprint arXiv:2311.01453.
  * Athey et al., (2019) Athey, S., Tibshirani, J., and Wager, S. (2019).  Generalized random forests.  The Annals of Statistics, 47(2):1148–1178.
  * Azriel et al., (2022) Azriel, D., Brown, L. D., Sklar, M., Berk, R., Buja, A., and Zhao, L. (2022).  Semi-supervised linear regression.  Journal of the American Statistical Association, 117(540):2238–2251.
  * Bycroft et al., (2018) Bycroft, C., Freeman, C., Petkova, D., Band, G., Elliott, L. T., Sharp, K., Motyer, A., Vukcevic, D., Delaneau, O., O’Connell, J., et al. (2018).  The uk biobank resource with deep phenotyping and genomic data.  Nature, 562(7726):203–209.
  * Caflisch, (1998) Caflisch, R. E. (1998).  Monte carlo and quasi-monte carlo methods.  Acta numerica, 7:1–49.
  * Cai and Guo, (2020) Cai, T. T. and Guo, Z. (2020).  Semisupervised inference for explained variance in high dimensional linear regression and its applications.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(2):391–419.
  * Castro et al., (2022) Castro, V. M., Gainer, V., Wattanasin, N., Benoit, B., Cagan, A., Ghosh, B., Goryachev, S., Metta, R., Park, H., Wang, D., et al. (2022).  The mass general brigham biobank portal: an i2b2-based data repository linking disparate and high-dimensional patient data to support multimodal analytics.  Journal of the American Medical Informatics Association, 29(4):643–651.
  * Chakrabortty and Cai, (2018) Chakrabortty, A. and Cai, T. (2018).  Efficient and adaptive linear regression in semi-supervised settings.  The Annals of Statistics, 46(4):1541–1572.
  * Chen et al., (2008) Chen, X., Hong, H., and Tarozzi, A. (2008).  Semiparametric efficiency in gmm models with auxiliary data.
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal, 21(1):C1–C68.
  * Chernozhukov et al., (2016) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., and Newey, W. K. (2016).  Double machine learning for treatment and causal parameters.  Technical report, cemmap working paper.
  * Deng et al., (2023) Deng, S., Ning, Y., Zhao, J., and Zhang, H. (2023).  Optimal and safe estimation for high-dimensional semi-supervised learning.  Journal of the American Statistical Association, pages 1–12.
  * Després et al., (2000) Després, J.-P., Lemieux, I., Dagenais, G.-R., Cantin, B., and Lamarche, B. (2000).  Hdl-cholesterol as a marker of coronary heart disease risk: the quebec cardiovascular study.  Atherosclerosis, 153(2):263–272.
  * Gan and Liang, (2023) Gan, F. and Liang, W. (2023).  Prediction de-correlated inference.  arXiv preprint arXiv:2312.06478.
  * Gehrmann et al., (2018) Gehrmann, S., Dernoncourt, F., Li, Y., Carlson, E. T., Wu, J. T., Welt, J., Foote Jr, J., Moseley, E. T., Grant, D. W., Tyler, P. D., et al. (2018).  Comparing deep learning and concept extraction based methods for patient phenotyping from clinical narratives.  PloS one, 13(2):e0192360.
  * Gronsbell et al., (2022) Gronsbell, J., Liu, M., Tian, L., and Cai, T. (2022).  Efficient evaluation of prediction rules in semi-supervised settings under stratified sampling.  Journal of the Royal Statistical Society: Series B: Statistical Methodology.
  * Haffner, (2000) Haffner, S. M. (2000).  Coronary heart disease in patients with diabetes.
  * Han et al., (2021) Han, L., Hou, J., Cho, K., Duan, R., and Cai, T. (2021).  Federated adaptive causal estimation (face) of target treatment effects.  arXiv preprint arXiv:2112.09313.
  * Hatt et al., (2022) Hatt, T., Berrevoets, J., Curth, A., Feuerriegel, S., and van der Schaar, M. (2022).  Combining observational and randomized data for estimating heterogeneous treatment effects.  arXiv preprint arXiv:2202.12891.
  * Hong et al., (2021) Hong, C., Rush, E., Liu, M., Zhou, D., Sun, J., Sonabend, A., Castro, V. M., Schubert, P., Panickan, V. A., Cai, T., et al. (2021).  Clinical knowledge extraction via sparse embedding regression (keser) with multi-center large scale electronic health record data.  NPJ digital medicine, 4(1):1–11.
  * Javanmard and Montanari, (2014) Javanmard, A. and Montanari, A. (2014).  Confidence intervals and hypothesis testing for high-dimensional regression.  The Journal of Machine Learning Research, 15(1):2869–2909.
  * Jin and Rothenhäusler, (2023) Jin, Y. and Rothenhäusler, D. (2023).  Modular regression: Improving linear models by incorporating auxiliary data.  Journal of Machine Learning Research, 24(351):1–52.
  * Johns, (1988) Johns, M. V. (1988).  Importance sampling for bootstrap confidence intervals.  Journal of the American Statistical Association, 83(403):709–714.
  * Johnson et al., (2016) Johnson, A. E., Pollard, T. J., Shen, L., Lehman, L.-w. H., Feng, M., Ghassemi, M., Moody, B., Szolovits, P., Anthony Celi, L., and Mark, R. G. (2016).  Mimic-iii, a freely accessible critical care database.  Scientific data, 3(1):1–9.
  * Kawakita and Kanamori, (2013) Kawakita, M. and Kanamori, T. (2013).  Semi-supervised learning with density-ratio estimation.  Machine learning, 91(2):189–209.
  * Kovvali, (2022) Kovvali, N. (2022).  Theory and applications of Gaussian quadrature methods.  Springer Nature.
  * Kundu and Chatterjee, (2023) Kundu, P. and Chatterjee, N. (2023).  Logistic regression analysis of two-phase studies using generalized method of moments.  Biometrics, 79(1):241–252.
  * (28) Li, S., Gilbert, P. B., and Luedtke, A. (2023a).  Data fusion using weakly aligned sources.  arXiv preprint arXiv:2308.14836.
  * Li and Luedtke, (2023) Li, S. and Luedtke, A. (2023).  Efficient estimation under data fusion.  Biometrika, 110(4):1041–1054.
  * (30) Li, Y., Yang, H., Yu, H., Huang, H., and Shen, Y. (2023b).  Penalized estimating equations for generalized linear models with multiple imputation.  The Annals of Applied Statistics, 17(3):2345–2363.
  * Liu et al., (2023) Liu, M., Zhang, Y., Liao, K. P., and Cai, T. (2023).  Augmented transfer regression learning with semi-non-parametric nuisance models.  Journal of Machine Learning Research, 24(293):1–50.
  * The All of Us Research Program Investigators, (2019) The All of Us Research Program Investigators (2019).  The “all of us” research program.  New England Journal of Medicine, 381(7):668–676.
  * Miao et al., (2023) Miao, J., Miao, X., Wu, Y., Zhao, J., and Lu, Q. (2023).  Assumption-lean and data-adaptive post-prediction inference.  arXiv preprint arXiv:2311.14220.
  * Mirza and Osindero, (2014) Mirza, M. and Osindero, S. (2014).  Conditional generative adversarial nets.  arXiv preprint arXiv:1411.1784.
  * Negahban et al., (2012) Negahban, S. N., Ravikumar, P., Wainwright, M. J., Yu, B., et al. (2012).  A unified framework for high-dimensional analysis of $m$-estimators with decomposable regularizers.  Statistical Science, 27(4):538–557.
  * Peters et al., (2014) Peters, S. A., Huxley, R. R., and Woodward, M. (2014).  Diabetes as risk factor for incident coronary heart disease in women compared with men: a systematic review and meta-analysis of 64 cohorts including 858,507 individuals and 28,203 coronary events.  Diabetologia, 57:1542–1551.
  * Qiu et al., (2023) Qiu, H., Tchetgen, E. T., and Dobriban, E. (2023).  Efficient and multiply robust risk estimation under general forms of dataset shift.  arXiv preprint arXiv:2306.16406.
  * Rader and Hovingh, (2014) Rader, D. J. and Hovingh, G. K. (2014).  Hdl and cardiovascular disease.  The Lancet, 384(9943):618–625.
  * Robins et al., (1994) Robins, J. M., Rotnitzky, A., and Zhao, L. P. (1994).  Estimation of regression coefficients when some regressors are not always observed.  Journal of the American statistical Association, 89(427):846–866.
  * Schmutz et al., (2022) Schmutz, H., Humbert, O., and Mattei, P.-A. (2022).  Don’t fear the unlabelled: safe semi-supervised learning via debiasing.  In The Eleventh International Conference on Learning Representations.
  * Shi et al., (2023) Shi, X., Pan, Z., and Miao, W. (2023).  Data integration in causal inference.  Wiley Interdisciplinary Reviews: Computational Statistics, 15(1):e1581.
  * Song et al., (2023) Song, S., Lin, Y., and Zhou, Y. (2023).  A general m-estimation theory in semi-supervised framework.  Journal of the American Statistical Association, pages 1–11.
  * Song et al., (2024) Song, S., Lin, Y., and Zhou, Y. (2024).  Semi-supervised inference for block-wise missing data without imputation.  Journal of Machine Learning Research, 25(99):1–36.
  * Sudlow et al., (2015) Sudlow, C., Gallacher, J., Allen, N., Beral, V., Burton, P., Danesh, J., Downey, P., Elliott, P., Green, J., Landray, M., et al. (2015).  Uk biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age.  PLoS medicine, 12(3):e1001779.
  * Tian et al., (2007) Tian, L., Cai, T., Goetghebeur, E., and Wei, L. (2007).  Model evaluation based on the sampling distribution of estimated absolute prediction error.  Biometrika, 94(2):297–311.
  * van Buuren and Groothuis-Oudshoorn, (2011) van Buuren, S. and Groothuis-Oudshoorn, K. (2011).  mice: Multivariate imputation by chained equations in r.  Journal of Statistical Software, 45(3):1–67.
  * Van de Geer et al., (2014) Van de Geer, S., Bühlmann, P., Ritov, Y., Dezeure, R., et al. (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  The Annals of Statistics, 42(3):1166–1202.
  * Van der Vaart, (2000) Van der Vaart, A. W. (2000).  Asymptotic statistics, volume 3.  Cambridge university press.
  * Verma et al., (2023) Verma, A., Huffman, J. E., Rodriguez, A., Conery, M., Liu, M., Ho, Y.-L., Kim, Y., Heise, D. A., Guare, L., Panickan, V. A., et al. (2023).  Diversity and scale: genetic architecture of 2,068 traits in the va million veteran program.  medRxiv.
  * Wang et al., (2024) Wang, X., Liu, M., Nogues, I.-E., Chen, T., Xiong, X., Bonzel, C.-L., Zhang, H., Hong, C., Xia, Y., Dahal, K., et al. (2024).  Heterogeneous associations between interleukin-6 receptor variants and phenotypes across ancestries and implications for therapy.  Scientific Reports, 14(1):8021.
  * White et al., (2011) White, I. R., Royston, P., and Wood, A. M. (2011).  Multiple imputation using chained equations: issues and guidance for practice.  Statistics in medicine, 30(4):377–399.
  * Wolfson et al., (2010) Wolfson, M., Wallace, S. E., Masca, N., Rowe, G., Sheehan, N. A., Ferretti, V., LaFlamme, P., Tobin, M. D., Macleod, J., Little, J., et al. (2010).  DataSHIELD: resolving a conflict in contemporary bioscience performing a pooled analysis of individual-level data without sharing the data.  International journal of epidemiology, 39(5):1372–1382.
  * Wu and Liu, (2023) Wu, D. and Liu, M. (2023).  Robust and efficient semi-supervised learning for ising model.  arXiv preprint arXiv:2311.15031.
  * Xue et al., (2021) Xue, F., Ma, R., and Li, H. (2021).  Statistical inference for high-dimensional linear regression with blockwise missing data.  arXiv preprint arXiv:2106.03344.
  * Xue and Qu, (2021) Xue, F. and Qu, A. (2021).  Integrating multisource block-wise missing data in model selection.  Journal of the American Statistical Association, 116(536):1914–1927.
  * Yang and Ding, (2019) Yang, S. and Ding, P. (2019).  Combining multiple observational data sources to estimate causal effects.  Journal of the American Statistical Association.
  * Yu et al., (2020) Yu, G., Li, Q., Shen, D., and Liu, Y. (2020).  Optimal sparse linear prediction for block-missing multi-modality data without imputation.  Journal of the American Statistical Association, 115(531):1406–1419.
  * Zengini et al., (2018) Zengini, E., Hatzikotoulas, K., Tachmazidou, I., Steinberg, J., Hartwig, F. P., Southam, L., Hackinger, S., Boer, C. G., Styrkarsdottir, U., Gilly, A., et al. (2018).  Genome-wide analyses using uk biobank data provide insights into the genetic architecture of osteoarthritis.  Nature genetics, 50(4):549–558.
  * Zhang and Zhang, (2014) Zhang, C.-H. and Zhang, S. S. (2014).  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):217–242.
  * Zhang et al., (2023) Zhang, H., Li, Y., Lu, W., and Lin, Q. (2023).  On the optimality of misspecified kernel ridge regression.  In International Conference on Machine Learning, pages 41331–41353. PMLR.
  * Zhang and Bradic, (2022) Zhang, Y. and Bradic, J. (2022).  High-dimensional semi-supervised learning: in search of optimal inference of the mean.  Biometrika, 109(2):387–403.