Submitted to the Annals of Statistics
arXiv: arXiv:1701.04889

EFFICIENT AND ADAPTIVE LINEAR REGRESSION IN
SEMI-SUPERVISED SETTINGS

arXiv:1701.04889v2 [stat.ME] 19 Aug 2017

By Abhishek Chakrabortty ∗,† and Tianxi Cai†
University of Pennsylvania and Harvard University
We consider the linear regression problem under semi-supervised
settings wherein the available data typically consists of: (i) a small or
moderate sized ‘labeled’ data, and (ii) a much larger sized ‘unlabeled’
data. Such data arises naturally from settings where the outcome, unlike the covariates, is expensive to obtain, a frequent scenario in modern studies involving large databases like electronic medical records
(EMR). Supervised estimators like the ordinary least squares (OLS)
estimator utilize only the labeled data. It is often of interest to investigate if and when the unlabeled data can be exploited to improve
estimation of the regression parameter in the adopted linear model.
In this paper, we propose a class of ‘Efficient and Adaptive SemiSupervised Estimators’ (EASE) to improve estimation efficiency. The
EASE are two-step estimators adaptive to model mis-specification,
leading to improved (optimal in some cases) efficiency under model
mis-specification, and equal (optimal) efficiency under a linear model.
This adaptive property, often unaddressed in the existing literature,
is crucial for advocating ‘safe’ use of the unlabeled data. The construction of EASE primarily involves a flexible ‘semi-non-parametric’
imputation, including a smoothing step that works well even when
the number of covariates is not small; and a follow up ‘refitting’ step
along with a cross-validation (CV) strategy both of which have useful practical as well as theoretical implications towards addressing
two important issues: under-smoothing and over-fitting. We establish
asymptotic results including consistency, asymptotic normality and
the adaptive properties of EASE. We also provide influence function
expansions and a ‘double’ CV strategy for inference. The results are
further validated through extensive simulations, followed by application to an EMR study on auto-immunity.

1. Introduction. In recent years, semi-supervised learning (SSL) has
emerged as an exciting new area of research in statistics and machine learning. A detailed discussion on SSL including its practical relevance, the pri∗

Corresponding author; previously at Harvard University during the time of this work.
This research was partially supported by National Institutes of Health grants R01
HL089778, U54 HG007963, and U01 HL121518.
MSC 2010 subject classifications: 62F35, 62J05, 62F12, 62G08.
Keywords and phrases: Semi-supervised Linear Regression, Semi-parametric Inference,
Model Mis-specification, Adaptive Estimation, Semi-non-parametric Imputation.
†

1

2

A. CHAKRABORTTY AND T. CAI

mary question of interest in SSL, and the existing relevant literature can
be found in Chapelle, Schölkopf and Zien (2006) and Zhu (2008). A typical semi-supervised (SS) setting is characterized by two types of available
data: (i) a small or moderate sized ‘labeled’ data, L, containing observations for both an outcome Y and a set of covariates X of interest, and (ii)
an ‘unlabeled’ data, U, of much larger size but having observations only
for the covariates X. By virtue of its large size, U essentially gives us the
distribution of X, denoted henceforth by PX . Such a setting arises naturally
whenever the covariates are easily available so that unlabeled data is plentiful, but the outcome is costly or difficult to obtain, thereby limiting the size
of L. This scenario is directly relevant to a variety of practical problems,
especially in the modern ‘big data’ era, with massive unlabeled datasets (often electronically recorded) becoming increasingly available and tractable.
A few familiar examples include machine learning problems like text mining,
web page classification, speech recognition, natural language processing etc.
Among biomedical applications, a particularly interesting problem where
SSL can be of great use is the statistical analysis of electronic medical records
(EMR) data. Endowed with a wealth of de-identified clinical and phenotype
data for large patient cohorts, EMR linked with bio-repositories are increasingly gaining popularity as rich resources of data for discovery research (Kohane, 2011). Such large scale datasets obtained in a cost-effective and timely
manner are of great importance in modern medical research for addressing
important questions such as the biological role of genetic variants in disease
susceptibility and progression (Kohane, 2011). However, one major bottleneck impeding EMR driven research is the difficulty in obtaining validated
phenotype information (Liao et al., 2010) since they are labor intensive or
expensive to obtain. Thus, gold standard labels and genomic measurements
are typically available only for a small subset nested within a large cohort. In
contrast, digitally recorded data on the clinical variables are often available
on all subjects, highlighting the necessity and utility of developing robust
SSL methods that can leverage such rich source of auxiliary information to
improve phenotype definition and estimation precision.
SSL primarily distinguishes from standard supervised methods by making
use of U, an information that is ignored by the latter. The ultimate question
of interest in SSL is to investigate if and when the information on PX in U
can be exploited to improve the efficiency over a given supervised approach.
In recent years, several graph based non-parametric SSL approaches have
been proposed (Zhu, 2005; Belkin, Niyogi and Sindhwani, 2006) for regression or classification. These approaches essentially target non-parametric SS
estimation of E(Y |X) and therefore, for provable improvement guarantees,

SEMI-SUPERVISED LINEAR REGRESSION

3

must rely implicitly or explicitly on assumptions relating PX to PY |X (the
conditional distribution of Y given X), as duly noted and characterized more
formally in Lafferty and Wasserman (2007). For non-parametric classification problems, the theoretical underpinnings of SSL including its scope and
the consequences of using U have been also studied earlier by Castelli and
Cover (1995, 1996). More parametric SS approaches, still aimed mostly at
prediction, have also been studied for classification, including the ‘generative
model’ approach (Nigam et al., 2000; Nigam, 2001) which is based on modeling the joint distribution of (Y, X) as an identifiable mixture of parametric
models, thereby implicitly relating PY |X and PX . However, these approaches
depend strongly on the validity of the assumed mixture model, violation of
which can actually degrade their performance compared to the supervised
approach (Cozman and Cohen, 2001; Cozman, Cohen and Cirelo, 2003).
However SS estimation problems, especially from a semi-parametric point
of view, has been somewhat less studied in SSL. Such problems are generally
aimed at estimating some (finite-dimensional) parameter θ0 ≡ θ0 (P), where
P = (PY |X , PX ), and the key to the potential usefulness of U in improving
estimation of θ0 lies in understanding when θ0 (P) relates to PX . For simple
parameters like θ0 (P) = E(Y ), unless E(Y |X) is a constant, θ0 clearly depends on PX and hence, improved SS estimation is possible compared to the
supervised estimator Y L , the sample mean of Y based on L. The situation is
however more subtle for other choices of θ0 , especially those where θ0 is the
target parameter corresponding to an underlying parametric working model
for PY |X . This includes the least squares parameter, as studied in this paper, targeted by a working linear model for E(Y |X). Such models are often
adopted due to their appealing simplicity and interpretability.
In general, for such cases, if the adopted working model for PY |X is correct and θ0 is not related to PX , then one cannot possibly gain through SSL
by using the knowledge of PX (Zhang and Oles, 2000; Seeger, 2002). On
the other hand, under model mis-specification, θ0 may inherently depend on
PX , and thus imply the potential utility of U in improving the estimation.
However, inappropriate use of U may lead to degradation of the estimation
precision. This therefore signifies the need for robust and efficient SS estimators that are adaptive to model mis-specification, so that they are as efficient
as the supervised estimator under the correct model and more efficient under model mis-specification. To the best of our knowledge, work done along
these lines is relatively scarce in the SSL literature, one notable exception
being the recent work of Kawakita and Kanamori (2013), where they use a
very different approach based on density ratio estimation, building on the
more restrictive approach of Sokolovska, Cappé and Yvon (2008). However,

4

A. CHAKRABORTTY AND T. CAI

as we observe in our simulation studies, the extent of the efficiency gain actually achieved by these approaches can be quite incremental, at least in finite
samples. Further, the seemingly unclear choice of the ideal (nuisance) model
to be used for density ratio estimation can also have a significant impact on
the performance, both finite sample and asymptotic, of these estimators.
We propose here a class of Efficient and Adaptive Semi-Supervised Estimators (EASE) in the context of linear regression problems. We essentially
adopt a semi-parametric perspective wherein the adopted linear ‘working’
model can be potentially mis-specified, and the goal is to obtain efficient and
adaptive SS estimators of the regression parameter through robust usage of
U. The EASE are two-step estimators with a simple and scalable construction based on a first step of ‘semi-non-parametric’ (SNP) imputation which
includes a smoothing step and a follow-up ‘refitting’ step. In the second step,
we regress the imputed outcomes against the covariates using the unlabeled
data to obtain our SNP imputation based SS estimator, and then further
combine it optimally with the supervised estimator to obtain the final EASE.
Dimension reduction methods are also employed in the smoothing step to
accommodate higher dimensional X, if necessary. Further, we extensively
adopt cross-validation (CV) techniques in the imputation, leading to some
useful theoretical properties (apart from practical benefits) typically not observed for smoothing based two-step estimators. We demonstrate that EASE
is guaranteed to be efficient and adaptive in the sense discussed above, and
also achieves semi-parametric optimality whenever the SNP imputation is
‘sufficient’ or the linear model holds. We also provide data adaptive methods
to optimally select the directions for smoothing when dimension reduction
is desired, and tools for inference with EASE.
The rest of this paper is organized as follows. In Section 2, we formulate
the SS linear regression problem. In Section 3, we construct a family of SS
estimators based on SNP imputation and establish all their properties, and
further propose the EASE as a refinement of these estimators. For all our
proposed estimators, we also address their associated inference procedures
based on ‘double’ CV methods. In Section 4, we discuss a kernel smoothing
based implementation of the SNP imputation and establish all its properties. In Section 5, we discuss SS dimension reduction techniques, useful for
implementing the SNP imputation. Simulation results and an application to
an EMR study are shown in Section 6, followed by concluding discussions
in Section 7. Proofs of all theoretical results and associated technical materials, and further numerical results and discussions are distributed in the
Appendix and the Supplementary Material [Chakrabortty and Cai (2017)].

SEMI-SUPERVISED LINEAR REGRESSION

5

2. Problem Set-up.
Data Representation. Let Y ∈ R denote the outcome random variable and
X ∈ Rp denote the covariate vector, where p is fixed, and let Z = (Y, X0 )0 .
Then the entire data available for analysis can be represented as S = (L∪U),
where L = {Zi ≡ (Yi , X0i )0 : i = 1, . . . , n} consists of n independent and
identically distributed (i.i.d.) observations from the joint distribution PZ of
Z, U = {Xi : i = n + 1, . . . , n + N } consists of N i.i.d. observations from PX ,
and L ⊥
⊥ U. Throughout, for notational convenience, we use the subscript ‘j’
to denote the unlabeled observations, and re-index without loss of generality
(w.l.o.g.) the N observations in U as: U = {Xj : j = n + 1, . . . , n + N }.
Assumption 2.1 (Basic Assumptions). (a) We assume that Z has finite
2nd moments and Σ ≡ Var(X) is positive definite, denoted as Σ  0. We
also assume, for simplicity, that X has a compact support X ⊆ Rp .
(b) We assume N  n i.e. n/N → 0 as n, N → ∞, and L and U arise
from the same underlying distribution, i.e. Z ∼ PZ for all subjects in S.
→
−→
−
−
Notations. Let Γ = E( X X 0 )  0, where ∀ v ∈ Rp , →
v = (1, v0 )0 ∈ R(p+1) .
Let L2 (PX ) denote the space of all R-valued measurable functions of X
having finite L2 norm with respect to (w.r.t.) PX , and for any g(.) ∈ L2 (PX ),
→
−→
−
let Σ(g)  0 denote the (p+1)×(p+1) matrix Γ−1 E[ X X 0 {Y −g(X)}2 ]Γ−1 .
Lastly, let k · k denote the L2 vector norm, and for any integer a ≥ 1, let
Ia denote the identity matrix of order a, and Na [µ, Ω] denote the a-variate
Gaussian distribution with mean µ ∈ Ra and covariance matrix Ωa×a  0.
Remark 2.1. Assumption 2.1 (b) enlists some fundamental characteristics of SS settings. Indeed, the condition of L and U being equally distributed has usually been an integral part of the definition of SS settings
(Chapelle, Schölkopf and Zien, 2006; Kawakita and Kanamori, 2013). Interpreted in missing data terminology, it entails that Y in U are ‘missing
completely at random’ (MCAR), with the missingness/labeling being typically by design. Interestingly, the crucial assumption of MCAR, although
commonly required, has often stayed implicit in the SSL literature (Lafferty
and Wasserman, 2007). It is important to note that while the SS set-up can
be viewed as a missing data problem, it is quite different from standard
ones, since with n/N → 0 i.e. |U|  |L|, the proportion of Y observed in S
tends to 0 in SSL. Hence, the ‘positivity assumption’ typical in missing data
theory, requiring this proportion to be bounded away from 0, is violated
here. It is also worth noting that owing to such violations, the analysis of
SS settings under more general missingness mechanisms such as ‘missing at

6

A. CHAKRABORTTY AND T. CAI

random’ (MAR) is considerably more complicated and to our knowledge, the
literature for SS estimation problems under such settings is virtually nonexistent. Furthermore, for such problems, the traditional goal in SSL, that
of improving upon a ‘supervised’ estimator, can become unclear without
MCAR, unless an appropriately weighted version of the supervised estimator is considered. Given these subtleties and the traditional assumptions
(often implicit) in SSL, the MCAR condition is assumed for most of this
paper, although a brief discussion on possible extensions of our proposed SS
estimators to MAR settings is provided in the Supplementary Material.
2.1. The Target Parameter and Its Supervised Estimator. We consider
the linear regression working model given by:
→
−
(2.1)
Y = X 0 θ + , with E( | X) = 0,
where, θ ∈ R(p+1) is an unknown regression parameter. Accounting for the
potential mis-specification of the working model (2.1), we define the target
parameter of interest as a model free parameter, as follows:
Definition 2.1. The target parameter θ 0 for linear regression may be
→
−
→
−
defined as the solution to the normal equations: E{ X(Y − X 0 θ)} = 0 in
→
−
θ ∈ R(p+1) , or equivalently, θ 0 = argmin E(Y − X 0 θ)2 .
θ∈R(p+1)

→
−
Existence and uniqueness of θ 0 in 2.1 is clear. Further, X 0 θ 0 is the L2
projection of E(Y |X) ∈ L2 (PX ) onto the subspace of all linear functions of X
and hence, is the best linear predictor of Y given X. The linear model (2.1) is
correct (else, mis-specified ) if and only if E(Y |X) lies in this space (in which
→
−
case, E(Y |X) = X 0 θ 0 ). When the model is correct, θ 0 depends only on PY |X ,
not on PX . Hence, improved estimation of θ 0 through SSL is impossible in
this case unless further assumptions relating θ 0 to PX are made. On the
other hand, under model mis-specification, the normal equations defining
θ 0 inherently depend on PX , thereby implying the potential utility of SSL
in improving the estimation of θ 0 in this case.
b the solution
The usual supervised estimator of θ 0 is the OLS estimator θ,
Pn →
−
→
−0
−1
in θ to the equation: n
i=1 X i (Yi − X i θ) = 0, the normal equations
based on L. Under Assumption 2.1 (a), it is well known that as n → ∞,
(2.2)

1

1

b − θ 0 ) = n− 2
n 2 (θ

n
X

 1
d
ψ 0 (Zi ) + Op n− 2
→ N(p+1) [0, Σ(gθ0 )],

i=1

→
−
→
−
→
−
where ψ 0 (Z) = Γ−1 { X(Y − X 0 θ 0 )} and gθ (X) = X 0 θ ∀ θ ∈ R(p+1) .

SEMI-SUPERVISED LINEAR REGRESSION

7

Our primary goal is to obtain efficient SS estimators of θ 0 using the entire
b It is worth noting that the
data S and compare their efficiencies to that of θ.
estimation efficiency of θ 0 also relates to the predictive performance of the
fitted linear model since its out-of-sample prediction error is directly related
to the mean squared error (w.r.t. the Σ metric) of the parameter estimate.
3. A Family of Imputation Based Semi-Supervised Estimators.
If Y in U were actually observed, then one would simply fit the working
model to the entire data in S for estimating θ 0 . Our general approach is
precisely motivated by this intuition. We first attempt to impute the missing
Y in U based on suitable training of L in step (I). Then in step (II), we fit the
linear model (2.1) to U with the imputed outcomes. Clearly, the imputation
is critical. Inaccurate imputation would lead to biased estimate of θ 0 , while
inadequate imputation would result in loss of efficiency. We next consider SS
estimators constructed under two imputation strategies for step (I) including
a fully non-parametric imputation based on kernel smoothing (KS), and a
semi-non-parametric (SNP) imputation that involves a smoothing step and
a follow up ‘refitting’ step. Although the construction of the final EASE is
based on the SNP imputation strategy, it is helpful to begin with a discussion
of the first strategy in order to appropriately motivate and elucidate the
discussion on EASE and the SNP imputation strategy.
3.1. A Simple SS Estimator via Fully Non-Parametric Imputation. We
present here an estimator based on a fully non-parametric imputation involving KS when p is small. For simplicity, we shall assume here that X is continuous with a density f (·). Let m(x) = E(Y | X = x) and l(x) = m(x)f (x).
Consider the local constant KS estimator of m(x),
1 Pn
b
l(x)
i=1 {Kh (Xi , x)}Yi
nhp
(3.1)
m(x)
b
=
=
,
P
n
1
fb(x)
i=1 Kh (Xi , x)
nhp
where Kh (u, v) = K{(u−v)/h} with K : Rp → R being some suitable kernel
function and h = h(n) > 0 being the bandwidth. With m(·)
b
as defined in
(3.1), we now fit (2.1) to the imputed unlabeled data: [{m(X
b j ), X0j }0 : j =
bnp of θ 0 as the solution in θ to:
n + 1, ..., n + N ] and obtain a SS estimator θ
(3.2)

n+N
−
→
−
1 X →
X j {m(X
b j ) − X 0j θ} = 0.
N
j=n+1

Here and throughout in our constructions of SS estimators, L with either
the true or the imputed Y is not included in the final fitting step mostly

8

A. CHAKRABORTTY AND T. CAI

due to technical convenience in the asymptotic analysis of our estimators,
and also due to the fact that the contribution of L, included in any form, in
the final fitting step is asymptotically negligible since n/N → 0.
bnp , we require uniform (in L∞ norm)
In order to study the properties of θ
convergence of m(·)
b
to m(·), a problem that has been extensively studied
in the non-parametric statistics literature (Newey, 1994; Andrews, 1995;
Masry, 1996; Hansen, 2008) under fairly general settings and assumptions.
In particular, we would assume the following regularity conditions to hold:
Assumption 3.1. (i) K(·) is a symmetric q th order kernel for some integer q ≥ 2. (ii) K(·) is bounded, Lipschitz continuous and has a bounded
support K ⊆ Rp . (iii) E(|Y |s ) < ∞ for some s > 2. E(|Y |s | X = x)f (x) and
f (x) are bounded on X . (iv) f (x) is bounded away from 0 on X . (v) m(·)
and f (·) are q times continuously differentiable with bounded q th derivatives
on some open set X0 ⊇ X . (vi) For any δ > 0, let Aδ ⊆ Rp denote the set
{(x − X)/δ : x ∈ X }. Then, for small enough δ, Aδ ⊇ K almost surely (a.s.).
Conditions (i)-(v) are fairly standard in the literature. In (v), the set X0 is
needed mostly to make the notion of differentiability well-defined, with both
m(·) and f (·) understood to have been analytically extended over (X0 \X ).
Condition (vi) implicitly controls the tail behaviour of X, requiring that
perturbations of X in the form of (X + δφ) with φ ∈ K (bounded) and δ
bnp .
small enough, belong to X a.s. [PX ]. We now present our result on θ
1

1

Theorem 3.1. Suppose n 2 hq → 0 and (log n)/(n 2 hp ) → 0 as n → ∞,
1
1
1
and let rn = n 2 hq + (log n)/(n 2 hp ) + (n/N ) 2 . Then, under Assumption 3.1,
(3.3)

n


X
1
d
bnp − θ 0 = n− 12
n2 θ
ψ eff (Zi ) + Op (rn ) → N(p+1) [0, Σ(m)],
i=1

→
−
where ψ eff (Z) = Γ−1 [ X{Y − m(X)}].
Remark 3.1. Theorem 3.1 establishes the efficient and adaptive nabnp . The asymptotic variance Σ(m) of θ
bnp satisfies Σ(g) − Σ(m) 
ture of θ
2
0 ∀ g(·) ∈ L (X) and the inequality is strict unless g(·) = m(·) a.s. [PX ].
bnp is asymptotically optimal among the class of all regular and
Hence, θ
asymptotically linear (RAL) estimators of θ 0 with influence function (IF) of
→
−
bnp is more
the form: Γ−1 [ X{Y − g(X)}] with g(·) ∈ L2 (PX ). In particular, θ
b
efficient than θ whenever (2.1) is mis-specified, and equally efficient when
(2.1) is correct i.e. m(·) = gθ0 (·). Further, it can also be shown that ψ eff (Z)

SEMI-SUPERVISED LINEAR REGRESSION

9

is the ‘efficient’ IF for estimating θ 0 under the semi-parametric model MX ≡
{(PY |X , PX ) : PX is known, PY |X is unrestricted upto Assumption 2.1 (a)}.
bnp also globally achieves the semi-parametric efficiency bound under
Thus, θ
MX . Lastly, note that at any parametric sub-model in MX that corresponds
b also achieves optimality, thus showing that under
to (2.1) being correct, θ
b if the linear model is correct.
MX , it is not possible to improve upon θ
Remark 3.2. The asymptotic results in Theorem 3.1 require a kernel of
order q > p and h smaller in order than the ‘optimal’ bandwidth order hopt =
O(n−1/(2q+p) ). This under-smoothing requirement, often encountered in twostep estimators involving a first-step smoothing (Newey, Hsieh and Robins,
1998), generally results in sub-optimal performance of m(.).
b
The optimal
under-smoothed bandwidth order for Theorem 3.1 is given by: O(n−1/(q+p) ).
3.2. SS Estimators Based on Semi-Non-Parametric (SNP) Imputation.
The simple and intuitive imputation strategy in Section 3.1 based on a fully
non-parametric p-dimensional KS is however often undesirable in practice
owing to the curse of dimensionality. In order to accommodate larger p, we
now propose a more flexible SNP imputation method involving a dimension
reduction, if needed, followed by a non-parametric calibration. An additional
‘refitting’ step is proposed to reduce the impact of bias from non-parametric
estimation and possibly inadequate imputation due to dimension reduction.
We also introduce some flexibility in terms of the smoothing methods, apart
from KS, that can be used for the non-parametric calibration.
Let r ≤ p be a fixed positive integer and let Pr = [p1 , .., pr ]p×r be any
rank r transformation matrix. Let XPr = P0r X. Given (r, Pr ), we may now
consider approximating the regression function E(Y |X) by smoothing Y
over the r dimensional XPr instead of the original X ∈ Rp . In general, Pr
can be user-defined and data dependent. A few reasonable choices of Pr
are discussed in Section 5. If Pr depends only on the distribution of X, it
may be assumed to be known given the SS setting considered. If Pr also
depends on the distribution of Y , then it needs to be estimated from L and
the smoothing needs to be performed using the estimated Pr .
For approximating E(Y |X), we may consider any reasonable smoothing
technique T . Some examples of T include KS, kernel machine regression and
smoothing splines. Let m(x; Pr ) denote the ‘target function’ for smoothing
Y over XPr using T . For notational simplicity, the dependence of m(x; Pr )
and other quantities on T is suppressed throughout. For T := KS, the appropriate target function is given by: m(x; Pr ) = mPr (P0r x), where mPr (z) ≡
E(Y | XPr = z). For basis function expansion based methods, m(x; Pr ) will
typically correspond to the L2 projection of m(x) ≡ E(Y | X = x) ∈ L2 (PX )

10

A. CHAKRABORTTY AND T. CAI

onto the functional space spanned by the basis functions associated with
T . The results in this section apply to any choice of T that satisfies the
required conditions. In Section 4, we provide more specific results for the
implementation of our methods using T := KS.
Note that we do not assume m(x; Pr ) = m(x) anywhere, and hence the
name ‘semi-non-parametric’. Clearly, with Pr = Ip and T := KS, it reduces
to a fully non-parametric approach. We next describe the two sub-steps
involved in step (I) of the SNP imputation: (Ia) smoothing, and (Ib) refitting.
b r and
(Ia) Smoothing Step. With Pr and m(x; Pr ) as defined above, let P
b
m(x;
b
Pr ) respectively denote their estimators based on L. In order to address
potential overfitting issues in the subsequent steps, we further consider generalized versions of these estimators based on K-fold CV for a given fixed integer K ≥ 1. For any K ≥ 2, let {Lk }K
k=1 denote a random partition of L into
−
K disjoint subsets of equal sizes, nK = n/K, with index sets {Ik }K
k=1 . Let Lk
−
denote the set excluding Lk with size nK = n − nK and respective index set
b r,k and m
b r,k ) denote the corresponding estimators based on
Ik− . Let P
b k (x; P
−
Lk . Further, for notational consistency, we define for K = 1, Lk = L−
k = L;
−
−
b
b
b
b r ).
b k (x; Pr,k ) = m(x;
b
P
Ik = Ik = {1, ..., n}; nK = nK = n; Pr,k = Pr and m
(Ib) Refitting Step. In this step, we fit the linear model to L using X as predictors and the estimated m(X; Pr ) as an offset. To motivate this, we recall
that the fully non-parametric imputation given in Section 3.1 consistently
estimates E(Y |X), the L2 projection onto a space that always contains the
→
−
working model space, i.e. the linear span of X. This need not be true for
the SNP imputation, since we do not assume m(X; Pr ) = m(X) necessarily.
The refitting step essentially ‘adjusts’ for this so that the final imputation,
combining the predictions from these two steps, targets a space that contains the working model space. In particular, for T := KS with r < p, this
step is critical to remove potential bias due to inadequate imputation.
Interestingly, it turns out that the refitting step should always be performed, even when m(X; Pr ) = m(X). It plays a crucial role in reducing
the bias of the resulting SS estimator due to the inherent bias from nonparametric curve estimation. In particular, for T := KS with any r ≤ p, it
ensures that a bandwidth of the optimal order can be used, thereby eliminating the under-smoothing issue as encountered in Section 3.1. The target
parameter for the refitting step is simply the regression coefficient obtained
from regressing the residual Y −m(X; Pr ) on X and may be defined as: η Pr ,
→
−
→
−
the solution in η ∈ R(p+1) to the equation: E[ X{Y − m(X; Pr ) − X 0 η}] = 0.

SEMI-SUPERVISED LINEAR REGRESSION

11

b (Pr ,K) , the solution in η to the equation:
For any K ≥ 1, we estimate η Pr as η
(3.4)

n

−1

K X
X
→
−
−0
b r,k ) − →
X i {Yi − m
b k (Xi ; P
X i η} = 0.
k=1 i∈Ik

For Xi ∈ Lk , the estimate of m(Xi ; Pr ) to be used as an offset is obtained
b r,k ) that is based on data in L− . For K ≥ 2, with L− ⊥
from m
b k (· ; P
k
k ⊥ Lk , the
residuals are thus estimated in a cross-validated manner. For K = 1 however,
b r ) is estimated using the entire L which can lead to considerable
m(·
b ;P
underestimation of the true residuals owing to over-fitting and consequently,
substantial finite sample bias in the resulting SS estimator of θ 0 . This bias
can be effectively reduced by using the CV approach with K ≥ 2. We next
estimate the target function for the SNP imputation given by:
−
µ(x; P ) = m(x; P ) + →
x 0η
(3.5)
as:
r

(3.6)

br,K ) = K−1
µ
b(x; P

r

K
X

Pr

−
b r,k ) + →
b (Pr ,K) ,
m
b k (x; P
x 0η

k=1

b r,k }K . For notational simplicity, we suppress throughout
br,K = {P
where P
k=1
K
the inherent dependence of µ
b(· ; ·) itself on K and {L−
k }k=1 . Note that
similar to m(X; Pr ), we also do not assume µ(X; Pr ) = m(X). Apart from
the geometric motivation for the refitting step and its technical role in bias
→
−
reduction, it also generally ensures the condition: E[ X{Y − µ(X; Pr }] = 0,
regardless of the true underlying m(X). This condition is a key requirement
for the asymptotic expansions, in Theorem 3.2, of our resulting SS estimabr,K ), we now construct our final SS estimator as follows.
tors. Using µ
b(· ; P
SS Estimator from SNP Imputation. In step (II), we fit the linear model to
br,K ), X0 }0 : j = n + 1, ..., n + N ]
the SNP imputed unlabeled data: [{b
µ(Xj ; P
j
b
and obtain a SS estimator θ (Pr ,K) of θ 0 given by:
n+N
−
−0
1 X →
b
br,K ) − →
(3.7) θ (Pr ,K) is the solution in θ to
X j {b
µ(Xj ; P
X j θ} = 0.
N
j=n+1

For convenience of further discussion, let us define: ∀ k ∈ {1, . . . , K},
(3.8)
(3.9)

b r,k ) = m
b r,k ) − m(x; Pr ) ∀ x ∈ X , and
b k (x; Pr , P
∆
b k (x; P
−b
−
b k (x) = →
b r,k ) − EX {→
b r,k )} ∀ x ∈ X ,
b k (x; Pr , P
G
x∆
X ∆k (X; Pr , P

b k (·)
where EX (·) denotes expectation w.r.t. X ∈ U. The dependence of G
b
on (Pr , Pr,k ) and PX is suppressed here for notational simplicity. We now
b(P ,K) .
present our main result summarizing the properties of θ
r

12

A. CHAKRABORTTY AND T. CAI

Theorem 3.2. Suppose that T satisfies: (i) supx∈X |m(x; Pr )| < ∞ and
b r ) − m(x; Pr )| = Op (cn ) for some cn = o(1). With G
b k (.)
(ii) supx∈X |m(x;
b
P
− 21 PK P
b
as in (3.9), define Gn,K = n
k=1
i∈Ik Gk (Xi ). Then, for any K ≥ 1,
(3.10)

n

1
2



n

X
− 12
b
θ (Pr ,K) − θ 0 = n
ψ(Zi ; Pr ) − Γ−1 Gn,K + Op (c∗n,K ),
i=1

1
1
→
−
where ψ(Z; Pr ) = Γ [ X{Y − µ(X; Pr )}] and c∗n,K = cn− + n− 2 + (n/N ) 2
K
= o(1). Further, for any fixed K ≥ 2, Gn,K = Op (cn− ), so that

−1

K

(3.11)





1
b(P ,K) − θ 0 = n− 12
n2 θ
r

n
X

ψ(Zi ; Pr ) + Op (cn− + c∗n,K ),

i=1

K

which converges in distribution to N(p+1) [0, Σ{µ(· ; Pr )}].
Remark 3.3. If the imputation is ‘sufficient’ so that µ(x; Pr ) = m(x),
b(P ,K) , for any K ≥ 2, enjoys the same set of optimality properties as
then θ
r
bnp (while requiring less stringent assumptions
those noted in Remark 3.1 for θ
about K(·) and h, if KS is used). If µ(x; Pr ) 6= m(x), then it is however unb(P ,K) is always more efficient than θ.
b This will be addressed
clear whether θ
r
in Section 3.3 where we develop the final EASE.
Remark 3.4. Apart from the fairly mild condition (i), Theorem 3.2 only
b r ) w.r.t. m(· ; Pr ) for establishing the
requires uniform consistency of m(·
b ;P
1
b(P ,K) for any K ≥ 2.
n 2 -consistency and asymptotic normality (CAN) of θ
r
The uniform consistency typically holds for a wide range of smoothing methods T under fairly general conditions. For T := KS in particular, we provide
explicit results in Section 4 under mild regularity conditions that allow the
use of any kernel order and the associated optimal bandwidth order. This is
a notable relaxation from the stringent requirements for Theorem 3.1 that
necessitate under-smoothing and the use of higher order kernels.
b(P ,1) has not yet been established.
Remark 3.5. The CAN property of θ
r
The term Gn,K in (3.10) behaves quite differently when K = 1, compared to
K ≥ 2 when it has a nice structure due to the inherent ‘cross-fitting’ involved,
and can be controlled easily, and quite generally, under mild conditions as
noted in Remark 3.4. For K = 1 however, Gn,K is simply a centered empirical
process devoid of any such structure and in general, controlling it requires
stronger conditions and the use of empirical process theory (see for instance
b(P ,1)
Van der Vaart (2000) for relevant results). We derive the properties of θ
r
for the case of T := KS in Theorem 4.2 using a different approach however,
specialized for KS estimators, in order to control Gn,1 .

13

SEMI-SUPERVISED LINEAR REGRESSION

3.3. Efficient and Adaptive Semi-Supervised Estimators (EASE). To ensure adaptivity even when µ(x; Pr ) 6= m(x), we now define the final EASE
b and θ
b(P ,K) . Specifically, for any fixed
as an optimal linear combination of θ
r
b(P ,K) (∆) = θ
b + ∆(θ
b(P ,K) − θ)
b is a CAN
(p + 1) × (p + 1) matrix ∆, θ
r
r
b and θ
b(P ,K) are, and an optimal ∆ can be
estimator of θ 0 whenever θ
r
selected easily to minimize the asymptotic variance of the combined estimator. For simplicity, we focus here on ∆ being a diagonal matrix with
bE
b
b
∆ = diag(δ1 , ..., δp+1 ). Then the EASE is defined as θ
(Pr ,K) ≡ θ (Pr ,K) (∆)
b being any consistent estimator (see Section 3.4 for details) of the
with ∆
minimizer ∆ = diag(δ 1 , ..., δ p+1 ), where ∀ 1 ≤ l ≤ (p + 1),
n
o
Cov ψ 0[l] (Z), ψ [l] (Z; Pr ) − ψ 0[l] (Z)
n
o
(3.12)
δ l = − lim
,
↓0
Var ψ [l] (Z; Pr ) − ψ 0[l] (Z) + 
and for any vector a, a[l] denotes its lth component. Note that in (3.12),
the  and the limit outside are included to formally account for the case:
ψ 0[l] (Z) = ψ [l] (Z, Pr ) a.s. [PZ ], when we define δ l = 0 for identifiability.
bE
b(P ,K) (∆) are asymptotiIt is straightforward to show that θ
and θ
(Pr ,K)

r

bE
cally equivalent, so that θ
(Pr ,K) is a RAL estimator of θ 0 satisfying:
n
 E

X
1
d
− 12
b
−
θ
=
n
ψ(Zi ; Pr , ∆) + op (1) → N(p+1) [0, ΣPr (∆)],
n2 θ
0
(Pr ,K)
i=1

as n → ∞, where ψ(Z; Pr , ∆) = ψ 0 (Z) + ∆{ψ(Z; Pr ) − ψ 0 (Z)} and
ΣPr (∆) = Var{ψ(Z; Pr , ∆)}. Note that when either the linear model holds
or the SNP imputation is sufficient, then ψ(Z; Pr , ∆) = ψ eff (Z), so that
bE
θ
is asymptotically optimal in the sense of Remark 3.1. Further, when
(Pr ,K)

bE
neither cases hold, θ
(Pr ,K) is no longer optimal, but is still efficient and
b Lastly, if the imputation is certain to be sufficient
adaptive compared to θ.
bE
b(P ,K) .
(for example, if r = p and T := KS), we may simply define θ
=θ
(Pr ,K)

r

Remark 3.6. It can be shown that under MX , defined in Remark 3.1,
the class of all possible IFs achievable by RAL estimators of θ 0 is given by:
IF θ0 ,MX = {ψ g (Z) ≡ ψ eff (Z)+g(X) : E{g(X)} = 0, g[j] (·) ∈ L2 (PX ) ∀ j}.
b θ
b(P ,K) and θ
bE
The IFs achieved by θ,
are clearly members of this class.
r

(Pr ,K)

The SNP imputation, for various choices of the imputation function µ(· ; Pr ),
b(P ,K) , θ
bE
therefore equips us with a family of RAL estimator pairs {θ
}
r

(Pr ,K)

14

A. CHAKRABORTTY AND T. CAI

bE
for estimating θ 0 . The IF of θ
(Pr ,K) is further guaranteed to dominate that
b and when µ(· ; Pr ) = m(·), it also dominates all other IFs ∈ IFθ ,M .
of θ,
0

X

3.4. Inference for EASE and the SNP Imputation Based SS Estimators.
b(P ,K)
We now provide procedures for making inference about θ 0 based on θ
r
E
b
and θ
obtained using K ≥ 2. We also employ a ‘double’ CV to over(Pr ,K)

come bias in variance estimation due to over-fitting. A key step involved in
the variance estimation is to obtain reasonable estimates of {µ(Xi ; Pr )}ni=1 .
b (Pr ,K) in (3.4) was constructed via CV, the corresponding estiAlthough η
br,K ) in (3.6), of µ(x; Pr ) is likely to be over-fitted for Xi ∈ L.
mate, µ
b(x; P
To construct bias corrected estimates of µ(Xi ; Pr ), we first obtain K separate doubly cross-validated estimates of η Pr , {b
η k(Pr ,K) : k = 1, ..., K}, with
P
b k(Pr ,K) , for each k, being the solution in η to k0 6=k Sk0 (η) = 0, where
η
Sk0 (η) =

P

→
−

i∈Ik0 X i {Yi

−0
b r,k0 ) − →
−m
b k0 (Xi ; P
X i η}

∀ k 0 ∈ {1, . . . , K}.

For each k and k 0 6= k, Sk0 (η) is constructed such that {Zi : i ∈ Ik0 } used for
b r,k0 ) that is based on L−0 ⊥
0
b k(Pr ,K) is independent of m
obtaining η
b k0 (· ; P
k ⊥ Lk .
Then, for each Xi ∈ Lk and k ∈ {1, . . . , K}, we may estimate µ(Xi ; Pr ) as:
−0 k
b r,k ) + →
br,K ) = m
b (Pr ,K) .
µ
bk (Xi ; P
b k (Xi ; P
X iη
b k(Pr ,K) to reduce over-fitting bias
We exclude Sk (η) in the construction of η
br,K )} which we now use for estimating the IFs.
in the residuals {Yi − µ
bk (Xi ; P
For each Zi ∈ Lk and k ∈ {1, .., K}, we estimate ψ 0 (Zi ) and ψ(Zi ; Pr ),
b and θ
b(P ,K) , respectively as:
the corresponding IFs of θ
r
−
→
− b
−
b (Zi ) = Γ
b (Zi ; Pr ) = Γ
b −1 {→
b −1 [→
br,K )}],
ψ
X i (Yi − X 0i θ)}
and ψ
X i {Yi −b
µk (Xi ; P
0
k
b denotes any consistent estimator of Γ from L and/or U (for example,
where Γ
Pn+N →
− →
−0
− →
−0
−1
b = Γn ≡ n−1 Pn →
b
Γ
i=1 X i X i based on L, or Γ = ΓN ≡ N
j=n+1 X j X j
based on U). Then, Σ{µ(· ; Pr )} in (3.11) may be consistently estimated as:
b
Σ{µ(·
; Pr )} = n−1

K X
X

b (Zi ; Pr )ψ
b 0 (Zi ; Pr ).
ψ
k
k

k=1 i∈Ik

To estimate the combination matrix ∆ in (3.12) and the asymptotic variance, ΣPr (∆), of EASE consistently, let us define, ∀ 1 ≤ l ≤ (p + 1),
P
P
b
b
b
σ
bl,12 = − n−1 K
k=1
i∈Ik ψ 0[l] (Zi ){ψ k[l] (Zi ; Pr ) − ψ 0[l] (Zi )},
P
P
2
b
b
σ
bl,22 = n−1 K
k=1
i∈Ik {ψ k[l] (Zi ; Pr ) − ψ 0[l] (Zi )} ,

15

SEMI-SUPERVISED LINEAR REGRESSION
1

and δbl = σ
bl,12 /(b
σl,22 + n ) for some sequence n → 0 with n 2 n → ∞. Then,
b = diag(δb1 , ..., δbp+1 ) and
we estimate ∆ and ΣPr (∆) respectively as: ∆
b P (∆)
b = n−1
Σ
r

K X
X

b (Zi ; Pr , ∆)
b 0 (Zi ; Pr , ∆),
b ψ
b
ψ
k
k

k=1 i∈Ik

b (Z; Pr , ∆)
b (Z) + ∆{
b (Z; Pr ) − ψ
b (Z)} ∀ k ∈ {1, . . . , K}.
b =ψ
b ψ
where ψ
k
0
k
0
Normal confidence intervals (CIs) for the parameters of interest can also be
constructed accordingly based on these variance estimates.
4. Implementation Based on KS. We next detail the specific implementation of the SNP imputation based on KS estimators. With T := KS,
the target function for the smoothing is given by: m(x; Pr ) = mPr (P0r x) ≡
E(Y | XPr = P0r x). For simplicity, we assume that XPr is continuous with a
density fPr (·) and support XPr ≡ {P0r x : x ∈ X } ⊆ Rr . Let us now consider
the following class of local constant KS estimators for m(x; Pr ):
1 P
b0
b0
− r
i∈Ik− {Kh (Pr,k Xi , Pr,k x)}Yi
nK h
b
(4.1) m
b k (x; Pr,k ) =
∀ 1 ≤ k ≤ K,
1 P
b0
b0
i∈I − Kh (Pr,k Xi , Pr,k x)
n− hr
K

k

where Kh (·) and h are as in Section 3.1 with K(·) now being a suitable
kernel on Rr . In the light of Theorem 3.2, we focus primarily on establishing
b r) ≡ m
b r,1 ) in (4.1) with K = 1,
the uniform consistency of m(x;
b
P
b 1 (x; P
b r . For establishing the
accounting for the additional estimation error from P
desired result, we shall assume the following regularity conditions to hold:
Assumption 4.1. (i) K(·) is a symmetric kernel of order q ≥ 2 with
finite q th moments. (ii) K(·) is bounded, integrable and is either Lipschitz
continuous with a compact support or, has a bounded derivative ∇K(·)
which satisfies: k∇K(z)k ≤ Λkzk−ρ ∀ z ∈ Rr with kzk > L, where Λ > 0,
L > 0 and ρ > 1 are some fixed constants, and k.k denotes the standard
L2 vector norm. (iii) XPr ⊆ Rr is compact. E(|Y |s ) < ∞ for some s > 2.
E(|Y |s | XPr = z)fPr (z) and fPr (z) are bounded on XPr . (iv) fPr (z) is
bounded away from 0 on XPr . (v) mPr (z) and fPr (z) are both q times
continuously differentiable with bounded q th derivatives on some open set
X0,Pr ⊇ XPr . Additional Conditions (required only when Pr needs to be
estimated): (vi) K(·) has a bounded and integrable derivative ∇K(·). (vii)
∇K(·) satisfies: k∇K(z1 ) − ∇K(z2 )k ≤ kz1 − z2 k φ(z1 ) ∀ z1 , z2 ∈ Rr such
that kz1 − z2 k ≤ L∗ , for some fixed constant L∗ > 0, and some bounded
and integrable function φ : Rr → R+ . (viii) ∇K(·) is Lipschitz continuous

16

A. CHAKRABORTTY AND T. CAI

on Rr . (ix) E(X | XPr = z) and E(XY | XPr = z) are both continuously
differentiable with bounded first derivatives on X0,Pr ⊇ XPr .
Assumption 4.1, mostly adopted from Hansen (2008), imposes some mild
smoothness and moment conditions most of which are fairly standard, except
perhaps the conditions on K(·) in (vi)-(viii) all of which are however satisfied
by the Gaussian kernel among others. We now propose the following result.
b r − Pr ) = Op (αn ) for some αn = o(1) with
Theorem 4.1. Suppose (P
αn = 0 identically if Pr is known. Let q be the order of the kernel K(.) in
(4.1) for some integer q ≥ 2. Define:

1
1

log n 2
log n 2
2 −(r+2)
an,1 = αn
+ αn h
+ αn , an,2 =
+ hq
nhr+2
nhr
and assume that each of the terms involved in an,1 = o(1) and an,2 = o(1).
b r ), based on (4.1), satisfies:
Then, under Assumption 4.1, m(x;
b
P
(4.2)

b r ) − m(x; Pr )| = Op (an,1 + an,2 ).
supx∈X |m(x;
b
P

b r)
Remark 4.1. Theorem 4.1 establishes the L∞ error rate of m(x;
b
P
under mild regularity conditions and restrictions on h. Among its various
b r ) at the
implications, the rate also ensures uniform consistency of m(x;
b
P
−1/(2q+r)
optimal bandwidth order: hopt = O(n
) for any kernel order q ≥ 2
and any r ≤ p, as long as αn = o(n−(r+2)/(4q+2r) ) which always includes:
1
αn = O(n− 2 ) and αn = 0. These two cases are particularly relevant in prac1
tice as Pr being finite dimensional, n 2 -consistent estimators of Pr should
typically exist. For both cases, using hopt results in an,1 to be of lower order
(for q > 2) or the same order (for q = 2) compared to that of the main term
an,2 , so that the usual optimal rate prevails as the overall error rate.
b(P ,K) for K = 1. We now address the CAN property of
Properties of θ
r
b(P ,K) for K = 1 under the KS framework. Based on (3.10) and Remark
θ
r
3.5, the only step required for this is to effectively control the term Gn,K in
(3.10). The following result is in this regard. It involves Lemmas A.2-A.3 as
the main technical tools which may themselves be of independent interest.
b r)
Theorem 4.2. Let K = 1, T := KS, Gn,K be as in (3.10), and m(x;
b
P
be the KS estimator based on (4.1). Let αn , an,1 and an,2 be as in Theorem
b r − Pr ) = Op (αn ). Assume that a∗ and a∗ are o(1), where
4.1 with (P
n,1
n,2
a∗n,1 = αn +

αn
1
2

n h(r+1)

1

1

1

+n 2 αn2 h−2 +n 2 a2n,1 +n 2 an,1 an,2

and

1

a∗n,2 = n 2 a2n,2 .

SEMI-SUPERVISED LINEAR REGRESSION

17

Then, under Assumption 4.1, Gn,K = Op (a∗n,1 + a∗n,2 ) = op (1). Further, let
c∗n,K be as in Theorem 3.2 with cn = (an,1 + an,2 ). Then, using (3.10),
(4.3)

n


X
1
b(P ,K) − θ 0 = n− 12
n2 θ
ψ(Zi , Pr ) + Op (c∗n,K + dn ),
r
i=1

1
d
b(P ,K) − θ 0 ) →
N(p+1) [0, Σ{µ(· ; Pr )}].
where dn = a∗n,1 + a∗n,2 . Hence, n 2 (θ
r

Remark 4.2. Note that the term a∗n,2 always requires q > r/2 in order
to converge to 0, thus showing the contrasting behavior of the case K = 1
compared to K ≥ 2 where no such higher order kernel restriction is required.
1
Nevertheless, when αn = O(n− 2 ) or αn = 0, the optimal bandwidth order:
hopt = O(n−1/(2q+r) ) can indeed be still used as long as q > r/2 is satisfied.
Despite these facts and all the theoretical guarantee in Theorem 4.2, emb(P ,1) can be substantially
pirical evidence however seems to suggest that θ
r
biased in finite samples, in part due to over-fitting.
Remark 4.3. Technical benefits of refitting and CV: Suppose that
Pr = Ip , so that the SNP imputation with T := KS is indeed sufficient.
Further, assume that all of Theorems 3.1-4.2 hold, so that the estimators
bnp , θ
b(P ,1) , and θ
b(P ,K) (K ≥ 2) are comparable and all asymptotically
θ
r
r
optimal. However, their constructions are quite different which can signifibnp is based on KS only, and
cantly affect their finite sample performances. θ
requires stringent under-smoothing and a kernel of order q > p (Remark
b(P ,1) is based on KS and refitting (although the KS itself is certain to
3.2); θ
r
be sufficient), and requires no under-smoothing but needs a (weaker) kernel
b(P ,K) (K ≥ 2) additionorder condition (q > p/2) (Remark 4.2); while θ
r
ally involves CV, and requires no under-smoothing or higher order kernel
conditions (Remark 3.4). This highlights the critical role played by refitting
and CV, apart from their primary roles in the SNP imputation, in removing any under-smoothing and/or higher order kernel restrictions when T :=
KS, and this continues to hold for any other (r, Pr ) as well. In particular,
it shows, rather surprisingly, that refitting should be performed in order to
avoid under-smoothing even if the smoothing is known to be sufficient.
Remark 4.4. As mentioned in Section 3.2, T := KS along with possible
dimension reduction is just one reasonable choice of T for implementing the
SNP imputation, all technical requirements for which have been thoroughly
established in Section 4. In general, other smoothing methods, as long as
the requirements are satisfied, can also be equally used as choices of T . One
such choice could be kernel machine (KM) regression (with possibly no use

18

A. CHAKRABORTTY AND T. CAI

of dimension reduction, as KM uses penalization to effectively regularize the
target even with Pr = Ip ). We leave its implementation details to the reader
as they are readily available in a multitude of references, and also skip any
theoretical treatment, considering the primary goal and scope of this paper.
However, detailed numerical results are presented in Section 6 for this choice
of T as well to illustrate the wider applicability of our proposed methods.
5. Choices of Pr : Dimension Reduction Techniques. We next discuss choosing and estimating the matrix Pr (r < p) to be used for dimension
reduction, if required, in the SNP imputation, and which can play an important role in the sufficiency of the imputation. Simple choices of Pr include r
leading principal component directions of X or any r canonical directions of
X. Note that under the SS setting, Pr is effectively known if it only involves
the distribution of X, as is true for these choices. We now focus primarily on
the case where Pr also depends on the distribution of Y and hence, is unknown. Such a choice of Pr is often desirable to ensure that the imputation is
as ‘sufficient’ as possible for predicting Y . Several reasonable choices of such
Pr and their estimation are possible based on sufficient dimension reduction
(s.d.r.) methods like Sliced Inverse Regression (SIR) (Li, 1991), Principal
Hessian Directions (PHD) (Li, 1992; Cook, 1998), Sliced Average Variance
Estimation (SAVE) (Cook and Weisberg, 1991; Cook and Lee, 1999) etc.
In particular, we focus here on SIR where the choice of Pr is given by: P0r
1
= Σ− 2 Pr , with Pr being the r leading eigenvectors of M = Var{E(X | Y )},
1
where X = Σ− 2 (X − µ), with µ = E(X), denotes the standardized version
of X. It is well known (Li, 1991) that these directions lead to an optimal
(in some appropriate sense) r-dimensional linear transformation of X that
can be predicted by Y . Apart from these general optimality, they also have
deeper implications in the context of s.d.r. We refer the reader to Li (1991)
and other relevant references in the s.d.r. literature for further details.
For estimating P0r , we consider the SIR algorithm of Li (1991) and further
− b
K
propose a SS modification to it. With K and {L−
k , Ik , Pr,k }k=1 as before,
−
b k ) denote the estimates of (µ, Σ) based on L and define X(k) =
let (b
µk , Σ
k
− 12
b (X − µ
b k ). Then, the original SIR algorithm estimates P0r based on L−
Σ
k
k
as follows: (i) Divide the range of {Yi }i∈I − into H slices {I1 , .., IH }, where
k

H may depend on n−
bh,k denote the proportion of
K . For 1 ≤ h ≤ H, let p
b
{Yi }i∈I − in slice Ih ; (ii) For each Ih , let Mh,k denote the sample average of
k
PH
(k)
b
b h,k M
b0
the set: {Xi ∈ L−
bh,k M
h=1 p
h,k
k : Yi ∈ Ih }; (iii) Estimate M as: Mk =
1

b0 = Σ
b−2 P
b r,k , where P
b r,k denotes the r leading eigenvectors
and P0r as: P
k
r,k

19

SEMI-SUPERVISED LINEAR REGRESSION

b k . However, the SIR algorithm often tends to give unstable estimates
of M
of P0r , especially for the directions corresponding to the smaller eigenvalues
of M. To improve the efficiency in estimating P0r , we now propose a semisupervised SIR (SS-SIR) algorithm as follows.
∗

− b
K
b ) denote the esSS-SIR Algorithm. Given {L−
µ∗k , Σ
k
k , Ik , Pr,k }k=1 , let (b
1
∗−
−
2
(k∗)
b
b k ). Then
timates of (µ, Σ) based on Lk ∪ U and define X
= Σk (X − µ
the SS-SIR proceeds as follows. Step (i) stays the same as in SIR. In step
∗ =Y
(ii), for each k, and each j ∈ {n + 1, ..., n + N }, we impute Yj as Yj,k
bij,k ,
(k∗)
(k∗)
b ∗ be the samwhere bij,k = argmin − kX
− X k2 . For each Ih , let M
i
(k∗)
{Xi

i∈Ik

ple average of the set:

in step (iii), we estimate M
∗− 12

j
h,k
(k∗)
−
∗
∈ Lk : Yi ∈ Ih } ∪ {Xj ∈ U : Yj,k
P
0
H
b∗
b∗ M
b∗ =
bh,k M
as: M
h=1 p
h,k h,k and
k

∈ Ih }. Then
then, P0r as:

b ∗.
b
b 0∗ = Σ
b ∗ denotes the r leading eigenvectors of M
b ∗ , where P
P
P
k
r,k
k
r,k
r,k
The SS-SIR algorithm aims to improve the estimation of P0r by making
use of U in step (ii) through a nearest neighbour approximation for the un−
observed Y in U using L−
k . With nK large enough and m(·) smooth enough,
the imputed and the true underlying Y should belong to the same slice with
a high probability. Thus, the set of X’s belonging to a particular slice is
now ‘enriched’ and consequently, improved estimation of M and P0r is expected. The proposed method based on a nearest neighbor approximation is
also highly scalable and while other smoothing based approximations may
be used, they can be computationally intensive. The SS-SIR algorithm is
1
fairly robust to the choice of H, and H = O(n 2 log n) seems to give fairly
satisfactory performance. The slices may be chosen to have equal width or
1
equal number of observations. For SIR, n 2 -consistency of the estimates are
well established (Li, 1991; Duan and Li, 1991; Zhu and Ng, 1995) for various
formulations under fairly general settings (without any model based assumptions). The theoretical properties of SS-SIR, although not derived here, are
expected to follow similarly. Our simulation results (not shown here) further
suggest that SS-SIR significantly outperforms SIR, leading to substantially
improved estimation of θ 0 from the proposed methods.
6. Numerical Studies.
6.1. Simulation Studies. We conducted extensive simulation studies to
examine the finite sample performance of our proposed point and interval estimation procedures as well as to compare with existing methods. Throughout we let n = 500, N = 10000, and considered p = 2, 10 and 20. For our
CV based methods, we let K = 5. The true values of the target parameter

20

A. CHAKRABORTTY AND T. CAI

θ 0 were estimated via monte carlo with a large sample size of 50, 000. For
each configuration, the results were summarized based on 500 replications.
Results for p = 2 are summarized in the Supplementary Material, and the
discussions below focus primarily on p = 10 and 20.
We generated X ∼ Np [0, Ip ] and restricted X to [−5, 5]p to ensure its
boundedness. Given X = x, we generated Y ∼ N1 [m(x), 1], where we considered four different choices of m(x) :
(i) Linear : m(x) = x0 bp ;
(ii) Non-linear one component (NL1C): m(x) = (x0 bp ) + (x0 bp )2 ;
(iii) Non-linear two component (NL2C): m(x) = (x0 bp )(1 + x0 δ p ); and
(iv) Non-linear three component (NL3C): m(x) = (x0 bp )(1 + x0 δ p ) + (x0 ω p )2 ;
(1)

where, for each setting, we considered bp = bp
(2)

≡ (10p/2 , 00p/2 )0 and bp =

bp ≡ 1p , and set δ p = (00p/2 , 10p/2 )0 and ω p = (1, 0, 1, 0, . . . , 1, 0)0p×1 , where
for any a, 1a = (1, . . . , 1)0a×1 and 0a = (0, . . . , 0)0a×1 . Through appropriate choices of bp , δ p and ω p , as applicable, these models can incorporate
commonly encountered linear, quadratic and interaction effects.
For each setting, we used two choices of the smoothing method: (a)
T := KS2,P2 denoting KS with 2-dimensional smoothing over P0r X ≡ P02 X,
where P2 was estimated via SIR with H = 100 slices of equal width, followb r,k )}K were obtained via KS using a Gaussian kernel;
ing which {m
b k (x; P
k=1
(b) T := KM where we let Pr = Ip and then estimated {m
b k (x; Ip )}K
k=1 using
kernel machine (KM) regression based on a radial basis function (RBF) kernel. Throughout, h for KS, and all tuning parameters for KM were selected
via least squares CV. For (a), with X ∼ Np [0, Ip ], results from Li (1991) imply that the SNP imputation with r = 2 is sufficient for models (i)-(iii), and
insufficient for model (iv). For comparison, we also implemented two other
SS estimators: the density ratio based “DRESS” estimator of Kawakita and
Kanamori (2013) and the estimator of Sokolovska, Cappé and Yvon (2008)
called “MSSL” by Kawakita and Kanamori (2013). The density ratio estimation for the DRESS estimator was implemented using either (i) linear bases
3
{1, (X[j] )pj=1 } (DRESS1 ); or (ii) cubic bases {1, (Xd[j] )p,
j=1,d=1 } (DRESS3 ).
First, we compare the various estimators with respect to their efficiencies
based on empirical mean squared error. In Table 1, we present the efficiencies
of the proposed SNP and EASE estimators as well as other SS estimators
relative to the OLS. As expected, under model mis-specification, our estimators are substantially more efficient than the OLS with the relative efficiency
(RE) as high as near 5 fold when p = 10 and 3 fold when p = 20, for the
non-linear models. The efficiency gain is generally lower for p = 20 than for
p = 10, likely a consequence of overfitting of the non-parametric estimators

SEMI-SUPERVISED LINEAR REGRESSION

21

involved in the SNP imputation for larger p. Comparing EASE to SNP, the
EASE generally perform better for both linear and non-linear settings, as
expected. Comparing the two smoothers, it appears that T := KM generally
attains higher efficiency compared to that of T := KS2,P2 . This is in part
due to the high variability in the SIR direction estimation which impacts
the performance of the resulting SS estimator in finite samples. Interestingly,
none of the existing SS estimators perform well with REs ranging only from
about 0.9 to 1.1 across all settings.
(a) p = 10
Setting
(I)

(II)

Models
Linear
NL1C
NL2C
NL3C
Linear
NL1C
NL2C
NL3C

OLS
(Ref.)
1
1
1
1
1
1
1
1

SNP EASE
(T := KS2,P2 )
0.895
0.983
4.481
4.424
2.683
2.700
2.772
2.795
0.841
0.989
4.511
4.585
3.596
3.634
3.280
3.301

SNP EASE
(T := KM)
0.772
0.985
4.501
5.543
4.268
5.055
4.481
5.560
0.657
0.993
4.416
5.471
4.405
5.497
4.636
5.566

Other SS Estimators
DRESS1 DRESS3 MSSL
0.982
0.927
0.982
1.136
1.110
1.135
1.120
1.016
1.119
1.102
1.025
1.103
0.981
0.924
0.981
1.132
1.030
1.130
1.127
1.042
1.128
1.110
1.079
1.109

(b) p = 20
Setting
(I)

(II)

Models
Linear
NL1C
NL2C
NL3C
Linear
NL1C
NL2C
NL3C

OLS
(Ref.)
1
1
1
1
1
1
1
1

SNP EASE
(T := KS2,P2 )
0.673
0.986
2.256
2.288
1.414
1.388
1.539
1.531
0.519
0.991
2.290
2.346
1.899
1.917
1.937
1.949

SNP EASE
(T := KM)
0.740
0.981
2.680
3.630
2.661
3.544
2.605
3.510
0.609
0.989
2.669
3.660
2.766
3.963
2.682
3.702

Other SS Estimators
DRESS1 DRESS3 MSSL
0.956
0.866
0.956
1.035
0.920
1.035
1.032
0.922
1.033
1.049
0.931
1.051
0.958
0.872
0.958
1.032
0.908
1.031
1.036
0.917
1.036
1.046
0.958
1.046

b(P ,K) (SNP) and θ
bE
Table 1: Efficiencies of θ
(Pr ,K) (EASE) using T := KS2,P2
r
b (OLS)
or T := KM, as well as DRESS1 , DRESS3 and MSSL, relative to θ
with respect to the empirical mean squared error (MSE) under models (i),
(1)
(2)
(ii), (iii) and (iv) each with: (I) bp = bp or, (II) bp = bp .
We next examine the performance of the proposed inference procedures. In
Table 2(a) and (b), we present the bias, empirical standard error (ESE), the
average of the estimated standard error (ASE) and the coverage probability
(CovP) of the 95% CIs for each component of θ 0 when p = 10 under the
linear and NL2C models. In general, the EASE with both the KS and the KM
smoothers have negligible biases although the KM based estimator appears
to have slightly lower biases. The ASEs are close to the ESEs and the CovPs
are close to the nominal level, suggesting that the variance estimators work

22

A. CHAKRABORTTY AND T. CAI

well in practice with K = 5.
(a) OLS and EASE for the linear model.
Parameter
α0 = 0
β01 = 1
β02 = 1
β03 = 1
β04 = 1
β05 = 1
β06 = 0
β07 = 0
β08 = 0
β09 = 0
β010 = 0

OLS
Bias
-0.001
0.002
-0.001
-0.001
-0.002
-0.004
-0.000
0.003
-0.001
-0.002
0.003

b
(θ)
ESE
0.043
0.044
0.044
0.046
0.045
0.048
0.045
0.046
0.045
0.047
0.045

EASE
Bias
-0.001
-0.003
-0.005
-0.005
-0.006
-0.008
-0.001
0.003
-0.001
-0.002
0.003

E

b(P ,K) ; T := KS2,P )
(θ
2
r
ESE
ASE CovP
0.043 0.044
0.95
0.045 0.044
0.94
0.044 0.044
0.94
0.046 0.044
0.95
0.045 0.044
0.94
0.049 0.044
0.92
0.045 0.044
0.94
0.046 0.044
0.93
0.045 0.044
0.95
0.048 0.044
0.94
0.045 0.044
0.94

bE
EASE (θ
(Pr ,K) ; T := KM)
Bias
ESE
ASE CovP
0.000 0.043 0.044
0.96
0.004 0.047 0.044
0.93
0.000 0.045 0.044
0.95
-0.004 0.045 0.044
0.94
0.001 0.047 0.044
0.94
-0.001 0.046 0.044
0.95
0.001 0.045 0.044
0.95
0.001 0.043 0.044
0.96
-0.000 0.048 0.044
0.94
0.000 0.045 0.045
0.95
-0.002 0.045 0.045
0.94

(b) OLS and EASE for the NL2C model.
Parameter
α0 = 0
β01 = 1
β02 = 1
β03 = 1
β04 = 1
β05 = 1
β06 = 0
β07 = 0
β08 = 0
β09 = 0
β010 = 0

OLS
Bias
-0.015
0.000
-0.004
-0.015
-0.001
0.013
-0.010
0.006
-0.008
0.002
-0.008

b
(θ)
ESE
0.239
0.260
0.269
0.249
0.267
0.260
0.281
0.277
0.277
0.279
0.272

EASE
Bias
-0.016
0.015
0.017
0.018
0.016
0.019
0.008
0.002
-0.004
0.003
0.002

E

b(P ,K) ; T := KS2,P )
(θ
2
r
ESE
ASE CovP
0.146 0.136
0.93
0.159 0.160
0.96
0.173 0.158
0.93
0.156 0.158
0.95
0.164 0.159
0.94
0.164 0.158
0.94
0.164 0.155
0.94
0.166 0.155
0.93
0.167 0.156
0.94
0.160 0.157
0.95
0.160 0.155
0.95

bE
EASE (θ
(Pr ,K) ; T := KM)
Bias
ESE
ASE CovP
0.013 0.105 0.096
0.93
-0.004 0.124 0.112
0.93
0.010 0.127 0.113
0.93
-0.000 0.118 0.113
0.95
0.007 0.124 0.113
0.93
0.002 0.120 0.113
0.94
0.005 0.119 0.112
0.94
0.011 0.116 0.111
0.95
-0.001 0.120 0.112
0.95
0.007 0.118 0.113
0.95
0.004 0.130 0.111
0.91

(c) All other SS estimators for the models in (a) and (b) above.
DRESS1
Bias
ESE
-0.001 0.043
-0.002 0.044
0.000 0.045
0.006 0.045
0.003 0.045
-0.004 0.047
-0.001 0.045
-0.000 0.048
-0.004 0.043
-0.001 0.048
-0.003 0.047

Linear Model
DRESS3
Bias
ESE
-0.001 0.044
-0.001 0.046
0.001 0.046
0.006 0.047
0.003 0.046
-0.004 0.049
-0.001 0.046
-0.001 0.050
-0.003 0.044
-0.001 0.049
-0.003 0.049

MSSL
Bias
ESE
-0.001 0.043
-0.002 0.044
0.000 0.045
0.006 0.045
0.003 0.045
-0.004 0.047
-0.001 0.045
-0.000 0.048
-0.004 0.043
-0.001 0.048
-0.003 0.047

DRESS1
Bias
ESE
-0.004 0.223
-0.014 0.266
0.005 0.257
-0.013 0.256
-0.005 0.262
0.002 0.250
-0.017 0.239
-0.022 0.260
-0.011 0.241
-0.020 0.256
-0.020 0.252

NL2C Model
DRESS3
Bias
ESE
-0.003 0.226
-0.009 0.279
0.006 0.266
-0.019 0.281
-0.007 0.274
-0.007 0.266
-0.009 0.247
-0.019 0.270
-0.013 0.261
-0.019 0.259
-0.022 0.269

MSSL
Bias
ESE
-0.004 0.223
-0.014 0.266
0.006 0.257
-0.011 0.256
-0.005 0.262
0.002 0.252
-0.017 0.239
-0.022 0.260
-0.010 0.241
-0.020 0.256
-0.020 0.252

Table 2: Coordinate-wise bias, ESE, ASE and CovP of EASE, obtained using
T := KS2,P2 or T := KM, for estimating θ 0 under the linear and NL2C
(1)
models with p = 10 and bp = bp . Shown also are the corresponding bias
and ESE of the OLS, as well as the DRESS1 , DRESS3 and MSSL estimators.

SEMI-SUPERVISED LINEAR REGRESSION

23

As shown in Table 2(c), the other SS estimators tend to have slightly
larger biases and substantially larger standard errors (SEs) compared to our
estimators under the NL2C model.
6.2. Application to EMR Data. We applied our proposed SS procedures
to an EMR study of rheumatoid arthritis (RA), a systemic auto-immune
disease (AD), conducted at the Partners HealthCare (Liao et al., 2010).
The study cohort consists of 3854 RA patients with blood samples stored.
The outcome of interest is the (logarithm of) anti-CCP (antibodies to cyclic
citrullinated polypeptide), a biomarker that is often used to determine subtypes of RA. Due to cost constraints, anti-CCP was measured only for a
random subset of n = 355 patients, thereby leading to a SS set-up. To investigate the validity of the MCAR assumption, we report in Table II in
the Supplementary Material summary measures of the distributions in the
labeled and unlabeled data for each of the predictors, as well as p-values
from various tests for assessing equality of those distributions. The results
suggest that the MCAR assumption is appropriate in this study.
We relate the log anti-CCP level to a set of p = 24 clinical variables X
related to ADs, including age, gender, race; total counts of codified and/or
narrative mentions extracted from physicians’ notes via natural language
processing (NLP) for various RA related conditions including RA, Lupus,
Polymyalgiarheumatica (PmR), Spondyloarthritis (SpA), as well as various
RA medications; indicators of seropositivity and radiological evidence of
erosion; mentions of rheumatoid factor (RF), as well as anti-CCP positivity
from prior medical history. Since the tests for RF and anti-CCP were not
always ordered, missing indicators for these variables were also included. All
count variables were transformed as: x → log(1 + x) to increase stability of
the model fitting. All predictors were normalized to have unit variance.
We obtained the OLS, the EASE using both T := KS2,P2 and T := KM
in the smoothing step, as well as the DRESS1 estimator for comparison.
For EASE, we again used K = 5 and for the KS2,P2 smoother, P2 was
obtained using SIR with H = 80 slices of equal width. In Table 3, we present
the coordinate-wise estimates of the regression parameters along with their
estimated SEs and the corresponding p-values based on these estimates.
Overall, the point estimators from all methods are quite close to each other.
Our proposed EASE, with both KS and KM smoothers, is substantially
more efficient than the OLS across all coordinates with efficiency ranging
from about 1.4 to 2.4. The DRESS1 estimator improved estimation for a
few coordinates but the efficiency remains comparable to the OLS for most
coordinates. This again suggests the advantage of our proposed estimators

24

A. CHAKRABORTTY AND T. CAI

compared to both OLS and other SS estimators.
Predictors
Age
Gender
Race
Lupus
PmR
RA
SpA
Other ADs
Erosion
Seropositivity
Anti-CCPprior
Anti-CCPmiss
RF
RFmiss
Azathioprine
Enbrel
Gold salts
Humira
Infliximab
Leflunomide
Methotrexate
Plaquenil
Sulfasalazine
Other meds.

Est
.105
-.032
-.041
.038
-.075
.015
-.137
-.022
.076
.056
.572
.527
.128
.085
-.080
.138
.138
-.051
.003
-.027
-.021
-.043
-.114
-.042

b
OLS (θ)
SE
.076
.059
.065
.066
.044
.089
.102
.078
.070
.062
.136
.123
.081
.085
.071
.070
.050
.068
.069
.069
.073
.069
.074
.074

Pval
.168
.589
.534
.563
.088
.862
.177
.775
.278
.370
.000
.000
.113
.316
.263
.048
.006
.453
.968
.697
.775
.540
.125
.570

EASE (KS2,P2 )
Est
SE
Pval
.106 .064 .099
-.028 .050 .570
-.042 .055 .452
.048 .052 .359
-.076 .031 .013
.012 .076 .879
-.133 .072 .063
-.024 .058 .679
.078 .059 .184
.054 .053 .310
.557 .110 .000
.520 .097 .000
.125 .066 .059
.084 .070 .233
-.074 .056 .185
.133 .058 .021
.136 .043 .002
-.049 .057 .391
.008 .057 .887
-.023 .058 .693
-.024 .061 .699
-.038 .057 .503
-.116 .063 .064
-.052 .060 .385

RE
1.40
1.41
1.40
1.59
2.07
1.37
2.02
1.79
1.44
1.37
1.54
1.61
1.49
1.46
1.62
1.48
1.37
1.43
1.50
1.40
1.42
1.47
1.39
1.52

Est
.094
-.027
-.044
.021
-.074
.008
-.128
-.018
.085
.041
.567
.508
.149
.137
-.075
.136
.147
-.057
.000
-.031
-.025
-.044
-.105
-.052

DRESS1
SE
.073
.058
.067
.063
.031
.080
.075
.067
.069
.061
.123
.115
.079
.080
.062
.073
.050
.067
.067
.071
.073
.070
.072
.071

Pval
.199
.638
.511
.731
.016
.923
.089
.792
.221
.496
.000
.000
.059
.088
.225
.064
.003
.389
.994
.660
.728
.532
.145
.466

RE
1.09
1.04
.95
1.11
2.10
1.23
1.82
1.35
1.03
1.05
1.23
1.15
1.05
1.12
1.33
.91
1.01
1.03
1.07
.93
1.01
.98
1.06
1.10

EASE (KM)
Est
SE
Pval
.104 .064 .103
-.031 .049 .524
-.040 .055 .462
.037 .051 .464
-.075 .030 .014
.016 .075 .832
-.136 .066 .038
-.022 .056 .692
.076 .058 .189
.055 .052 .296
.568 .107 .000
.523 .096 .000
.127 .066 .054
.084 .070 .231
-.079 .053 .132
.137 .057 .017
.137 .042 .001
-.051 .056 .360
.003 .055 .959
-.026 .057 .644
-.022 .060 .720
-.042 .057 .460
-.113 .061 .065
-.042 .059 .473

RE
1.42
1.44
1.41
1.70
2.04
1.30
2.37
1.93
1.47
1.41
1.60
1.64
1.49
1.48
1.83
1.49
1.40
1.49
1.57
1.45
1.46
1.48
1.45
1.59

Table 3: Estimates (Est) of the regression coefficients based on OLS, EASE
obtained using either T := KS2,P2 or T := KM, as well as DRESS1 , along
with their estimated standard errors (SE) and the corresponding p-values
(Pval.) for testing the null effect of each predictor. Shown also are the relative
efficiencies (RE) of all the estimators compared to the OLS.
We also estimated the prediction errors (PEs) for each of the fitted linear
models based on the aforementioned estimation methods via CV. To remove
potential randomness in the CV partitions, we averaged over 10 replications
of leave-5-out CV estimates. The PE was about 1.28 for EASE with both
smoothers, 1.29 for OLS and 1.30 for DRESS1 . For prediction purposes, we
may also directly employ non-parametric estimates of the conditional mean
rather than the fitted linear models. The PE in fact is slightly larger when we
b r ) or µ
br,K ). The PE was 1.34 for KS and 1.33 for KM based
use m(x;
b
P
b(x; P
b
br,K ). This
on m(x;
b
Pr ), and 1.30 for KS and 1.28 for KM based on µ
b(x; P
confirms that while the linear model may be mis-specified, it may often be
preferable to non-parametric models in practice as it may achieve simplicity
without substantial loss in prediction performance.
7. Discussion. We have developed in this paper an efficient and adaptive estimation strategy for the SS linear regression problem. The adaptive
property possessed by the proposed EASE is crucial for advocating ‘safe’ use
of the unlabeled data and is often unaddressed in the existing literature. In

SEMI-SUPERVISED LINEAR REGRESSION

25

general, the magnitude of the efficiency gain with EASE depends on the inherent degree of non-linearity in E(Y |X) and the extent of sufficiency of the
underlying SNP imputation. In particular, if the imputation is sufficient or
bE
the working linear model is correct, θ
(Pr ,K) is further optimal among a wide
class of estimators. We obtained theoretical results along with IF expansions
b(P ,K) and θ
bE
for θ
(Pr ,K) substantiating all our claims and also validated them
r
based on numerical studies. The double CV method further facilitates accurate inference, overcoming potential over-fitting issues in finite samples due
to smoothing. An R code for implementing EASE is available upon request.
The proposed SNP imputation, the key component of EASE, apart from
being flexible and scalable, enjoys several useful properties. The refitting
b(P ,K) , and for
step and CV play a crucial role in reducing the bias of θ
r
T := KS in particular, eradicate any under-smoothing or higher order kernel requirements: two undesirable, yet often inevitable, conditions required
1
for n 2 -consistency of two-step estimators based on a first step of smoothing.
b(P ,1) compared
Theorem 4.2, apart from showing the distinct behaviour of θ
r
b(P ,K) for K ≥ 2, also highlights the key role of CV in completely removto θ
r
ing kernel order restrictions, apart from addressing over-fitting issues. The
error rates in the results of Theorems 4.1-4.2 are quite sharp and account
b r . The regularity conditions required are also
for any estimation error from P
fairly mild and standard in the literature. The continuity assumption on X
in Sections 3.1 and 4 is mostly for the convenience of proofs, and the results
continue to hold for more general X. Lastly, while we have focussed here on
linear regression for simplicity, our methods can indeed be easily adapted to
other regression problems such as logistic regression for binary outcomes.
When the goal is solely that of prediction, one obviously does not have to
employ linear regression models, and models that incorporate non-linear efb r ) or µ
br,K ),
fects can be helpful. For such settings, the estimators m(x;
b
P
b(x; P
obtained as by-products of our SNP imputation, can themselves serve as potentially useful non-linear predictors. These SNP estimators may substantially outperform naive non-parametric estimators such as a p-dimensional
KS estimator, as demonstrated in Table III of the Supplementary Material
for the models considered in our simulation studies. In practice, when the
covariates are substantially correlated and the dimension of p is not small
as in the EMR example, it is unclear whether non-linear models necessarily
provide better prediction performance than the linear models. Under such
settings, the linear model also has a clear advantage due to its simplicity. Furthermore, while prediction is a vitally important goal of predictive
modeling, association analysis under interpretable models is key to clinical
studies for discovery research and efficient estimation of the corresponding

26

A. CHAKRABORTTY AND T. CAI

model parameters remains an important task.
b(P ,K) . While (3.11)
We end with a comment on the choice of K ≥ 2 in θ
r
holds for any K ≥ 2, the error term in (3.11) depends on K through cn−
1

K

cn− are
and more precisely, through e
cn− = K 2 cn− . Since K is fixed, cn− and e
K
K
K
K
asymptotically equivalent. But for a given n, cn− is expected to decrease with
K
K, while e
cn− is likely to increase. It is however desirable that both are small
K
since cn− inherently controls the efficiency of the SNP imputation, while e
cn−
K
K
b
directly controls the bias of θ (P ,K) . Hence, a reasonable choice of K ≥ 2 may
r

c2n− ) for some λ ≥ 0. Since the (first order)
be based on minimizing: (c2n− +λe
K

K

b(P ,K) is independent of K, this is equivalent to a
asymptotic variance of θ
r
b(P ,K) with λ denoting
penalized minimization of the asymptotic MSE of θ
r
the weightage of the (lower order) bias relative to the (first order) variance.
In general, the optimal K should be inversely related to λ. Conversely, choice
of any K may be viewed to have an associated regularization effect (through
λ) resulting in a ‘variance-bias trade-off’ with smaller K leading to lower bias
at the cost of some efficiency, and higher K leading to improved efficiency
in lieu of some bias. In practice, we find that K = 5 works well, and K = 10
tends to give slightly smaller MSE at the cost of increased bias.
APPENDIX A
A.1. Preliminaries. The following Lemmas A.1-A.3 would be useful
in the proofs of the main theorems. The proofs of these lemmas, as well as
Theorems 3.1, 4.1 and 4.2, can be found in the Supplementary Material.
Lemma A.1. Let Z ∈ Rl be any random vector and g(Z) ∈ Rd be any
measurable function of Z, where l and d are fixed. Let Sn = {Zi }ni=1 ⊥
⊥ Sm =
{Zj }m
be
two
random
samples
of
n
and
m
i.i.d.
observations
of
Z
respecj=1
bn (·) be any estimator of g(·) based on Sn such that the random
tively. Let g
sequence: Tbn ≡ supz∈χ kb
gn (z)k is Op (1), where χ ⊆ Rl denotes the support
b n,m denote the (double) random sequence: m−1 P
bn (Zj ),
of Z. Let G
Zj ∈Sm g
b n,m ) = EZ {b
and let Gn denote the random sequence: ESm (G
gn (Z)}, where
EZ (·) denotes expectation w.r.t. Z ∈ Sm ⊥
⊥ Sn , and all expectations involved
are assumed to be finite almost surely (a.s.) [Sn ] ∀ n.
1
Then: (a) Gn,m − Gn = Op (m− 2 ), and (b) as long as g(.) has finite 2nd
P
1
moments, m−1 Zj ∈Sm g(Zj ) − EZ {g(Z)} = Op (m− 2 ).
Controlling Empirical Processes Indexed by KS Estimators. The next two
lemmas would be useful in the proof of Theorem 4.2. They may also be

27

SEMI-SUPERVISED LINEAR REGRESSION

of more general use in other applications that involve controlling empirical
processes indexed by kernel estimators - both linear and ratio-type estimators, where the smoothing is further allowed to be performed over a possibly
lower dimensional and estimated transformation of the original covariate X.
These allowances make the technical analyses of such processes considerably
more involved and nuanced. The results of these lemmas and the techniques
used in their proofs may therefore be of independent general interest.
Suppose Assumption 2.1 (a) holds, and consider the KS framework introP
(%)
duced in Section 4. Let ϕ
ePr (w) = (nhr )−1 ni=1 Kh (w, P0r Xi )Yi% , for % =
(0)
(1)
0, 1. Let feP (·) = ϕ
e (·), e
lP (·) = ϕ
e (·) and m
e P (·) = e
lP (·)/feP (·). Next,
r

let

(0)
ϕPr (·)

Pr

= fPr (·) and

r

(1)
ϕPr (·)

Pr

r

r

r

= lPr (·), where lPr (·) = mPr (·)fPr (·). For
(%)

(%)

0
each % ∈ {0, 1}, let
e(%) (x; Pr ) = ϕ
ePr (P0r x).
r ) = ϕPr (Pr x) and ϕ
Further, let fe(·) = ϕ
e(0) (·), e
l(·) = ϕ
e(1) (·) and m(·)
e =e
l(·)/fe(·).
Lastly, let Pn denote the empirical probability measure on Rp based on
{Xi }ni=1 , and for any measurable (and possibly vector-valued) function γ(·)
1 R
of X, where γ(·) can be random itself, let G∗n (γ) = n 2 γ(x)(Pn − PX )(dx),
1
the (centered) n 2 -scaled empirical process indexed by γ(·). Lemmas A.2-A.3
together, among other more general implications, establish explicit rates of
g (·)−g(·)}, for any linear or ratio-type kernel
convergence of the quantity G∗n {b
estimator gb(·) of the type discussed above and its corresponding target g(·).

ϕ(%) (x; P

Lemma A.2. Consider the set-up introduced above. For any fixed integer
d ≥ 1, let λ(·) be any Rd -valued measurable function of X that is bounded
1
1
1
(1)
a.s. [PX ]. Define: bn = n− 2 h−r + hq and an,2 = (log n) 2 (nhr )− 2 + hq .
1
(1)
Assume bn = o(1) for (A.1) and n 2 a2n,2 = o(1) for (A.2) below. Then,
under Assumption 4.1 (i)-(v), and ∀ % ∈ {0, 1},
(A.1)

G∗n [λ(·){ϕ
e(%) (· ; Pr ) − ϕ(%) (· ; Pr )}] = Op (b(1)
n ) = op (1), and

(A.2)

G∗n [λ(·){m(·
e ; Pr ) − m(· ; Pr )}] = Op (n 2 a2n,2 ) = op (1).

1

b r ) = (nhr )−1 Pn Kh (P
b 0 x, P
b 0 Xi )Y % ∀ % ∈ {0, 1}, where P
br
Let ϕ
b(%) (x; P
r
r
i=1
i
is as in Section 3.2 and all other notations are the same as in the set-up of
b r) = ϕ
b r ) and b
b r) = ϕ
b r ). Then:
Lemma A.2. Let fb(x; P
b(0) (x; P
l(x; P
b(1) (x; P
b r ) be as
Lemma A.3. Consider the set-up of Lemma A.2. Let ϕ
b(%) (x; P
b r − Pr ) = Op (αn ) for
above, and let λ(·) be as in Lemma A.2. Suppose (P
1
(2)
(2)
some αn = o(1). Assume bn = o(1), where bn = αn + n− 2 αn h−(r+1) +
1
n 2 αn2 (h−2 + n−1 h−(r+2) ). Then, under Assumption 4.1,
b r) − ϕ
(A.3) G∗n [λ(·){ϕ
b(%) (· ; P
e(%) (· ; Pr )}] = Op (b(2)
n ) = op (1) ∀ % ∈ {0, 1}.

28

A. CHAKRABORTTY AND T. CAI

A.2. Proof of Theorem 3.2. Let Γn =
n

T(1)
n =

1
n

Pn →
− →
−0
i=1 X i X i , and

K

−
− b
1 X→
1XX→
(2)
b r,k ).
X i {Yi − µ(Xi ; Pr )} , Tn,K =
X i ∆k (Xi ; Pr , P
n
n
i=1

k=1 i∈Ik

Then, using (3.4)-(3.8), it is straightforward to see that:
(A.4)
(A.5)

→
−
→
−
→
−
E[ X{Y − µ(X; Pr )}] ≡ E[ X{Y − m(X; Pr ) − X 0 η Pr }] = 0, and


(2)
b (Pr ,K) − η Pr = T(1)
Γn η
n − Tn,K .

Under (A.4), Assumptions 2.1 (a) and (i), it follows from Lemma A.1 (b)
1
(1)
that Tn = Op (n− 2 ). Next, due to assumption (ii) and boundedness of X,
(2)

kTn,K k ≤ n−1

K X
X

−
b r,k )|} = Op (c − ).
b k (x; Pr , P
supx∈X {k→
x k |∆
n
K

k=1 i∈Ik
1

Finally, under Assumption 2.1 (a), we have: Γn = Γ+Op (n− 2 ) using Lemma
1
−1
A.1 (b). Further, since Γn is invertible a.s., Γ−1
+ Op (n− 2 ). Using all
n =Γ
(1)
(2)
(1)
−1
these facts, we then have: (b
η (Pr ,K) − η Pr ) = Γ−1
n (Tn − Tn,K ) = Γ (Tn −
1

(2)

1

Tn,K ) + Op {n− 2 (n− 2 + cn− )}. Thus,
K

(A.6)

(2)

1

−1
(b
η (Pr ,K) − η Pr ) = Γ−1 (T(1)
+ n− 2 cn− ).
n − Tn,K ) + Op (n
K

Next, let us define:
ΓN = N

−1

n+N
X

n+N
X →
→
− →
−0
−
→
−
(1)
−1
X j X j , RN = N
X j {µ(Xj ; Pr ) − X 0j θ 0 },

j=n+1

b (K) = N −1
and R
N,n

j=n+1
n+N
X

→
−
br,K ) − µ(Xj ; Pr )}.
X j {b
µ(Xj ; P

j=n+1

Then, using (3.7), we have:
b(P ,K) − θ 0 ) = N −1
ΓN ( θ
r

n+N
X

→
−
−0
(1)
b (K) .
br,K ) − →
X j [b
µ(Xj ; P
X j θ 0 ] = RN + R
N,n

j=n+1

b (K) = ΓN (b
b (K) , where
Next, using (3.4)-(3.8), we have: R
η (Pr ,K) − η Pr ) + S
N,n
N,n
− b
b (K) = K−1 PK {N −1 Pn+N →
b
S
k=1
j=n+1 X j ∆k (Xj ; Pr , Pr,k )}.
N,n

29

SEMI-SUPERVISED LINEAR REGRESSION
(1)
b(P ,K) − θ 0 ) = ΓN (b
b (K) .
Hence, we have: ΓN (θ
η (Pr ,K) − η Pr ) + RN + S
r
N,n

Now, under assumptions (i)-(ii) and Assumption 2.1 (a), we have:
PK
→
−b
b
(I)
k=1 supx∈X k x ∆k (x; Pr , Pr,k )k = Op (1),
− 21
b (K) = K−1 PK S
b∗
so that using Lemma A.1 (a), S
), where
k=1 n,k + Op (N
N,n
→
−
∗
b = EX { X ∆
b r,k )} ∀ 1 ≤ k ≤ K;
b k (X; Pr , P
S
n,k

1
1
→
−
→
−
(1)
(II) RN = E[ X{µ(X; Pr ) − X 0 θ 0 }] + Op (N − 2 ) = Op (N − 2 )
→
−
→
−
from Lemma A.1 (b) and E[ X{µ(X; Pr ) − X 0 θ 0 }] = 0 due to (A.4) and 2.1;
1
−1
and lastly, (III) Γ−1
+ Op (N − 2 ). It then follows from (I)-(III) that
N =Γ

b(P ,K) − θ 0 = (b
θ
η (Pr ,K) − η Pr ) + K−1 Γ−1
r

(A.7)

K
X

b ∗ + Op (N − 21 ).
S
n,k

k=1

Using (A.6) and (3.9) in (A.7), we then have:
n
K n
o
X
X
1 X b
−1 1
b(P ,K) − θ 0 = 1
θ
G
(X
)
+ Op (bn,K ),
ψ(Z
;
P
)
−
Γ
i
i
r
k
r
n
K
nK
i=1

where bn,K = n−1 + n

i∈Ik

k=1

− 21

cn− + N

− 12

. It follows, as claimed in (3.10), that

K

1

1

b(P ,K) − θ 0 ) = n− 2
n 2 (θ
r

(A.8)

n
X

ψ(Zi ; Pr ) − Γ−1 Gn,K + Op (c∗n,K )

i=1

We next show that Gn,K = Op (cn− ) for any fixed K ≥ 2. To this end,
K
1 P
(n)
b k (Xi ), D
b r,k )| and C =
b k = supx∈X |∆
b k (x; Pr , P
let Tk = (nK )− 2 i∈Ik G
→
−
supx∈X k x k < ∞. For any subset A ⊆ L, let PA denote the joint distribution
of the observations in A, and let EA (·) denote expectation w.r.t. PA . By
1 P
(n)
definition, Gn,K = K− 2 K
T = Op (cn− ) if and only if given any  > 0,
 k=1 k
 K
∃ M > 0 such that P kGn,K k > M cn− ≤  ∀ n. Note that for any M > 0,
K





− 12

P kGn,K k > M cn− ≤ P K
K

≤

K
X

− 21

P K

k=1
p+1
K X
X

(A.9) =

k

k=1 l=1

M cn−

(n)

EL− PLk

!
(n)
kTk k

> M cn−
K

k=1

kTk k >
"

K
X

p+1
K X
X

!

K

≤

K

(n)

P |Tk[l] | >

k=1 l=1

(
(n)
|Tk[l] |

(

>

)#

M cn−
K

1

1

K 2 (p + 1) 2

L−
k

,

)

M cn−
K

1
2

1

K (p + 1) 2

30

A. CHAKRABORTTY AND T. CAI

where the steps follow from repeated use of Bonferroni’s inequality and other
1

(n)

standard arguments. Now, conditional on L−
⊥ Lk , with K ≥ 2), nK2 Tk is
k (⊥
→
− b
b r,k )}i∈I which,
a centered sum of the i.i.d. random vectors { X i ∆k (Xi ; Pr , P
k
bk < ∞
due to assumption (ii) and the compactness of X , are bounded by: C D
(n)
a.s. [PL− ] ∀ k, n. Hence, applying Hoeffding’s inequality to Tk[l] ∀ l, we have:
k


(
)
2 c2


M
M
c
−
−
nK
nK
(n)
−
L
≤
2
exp
−
(A.10) PLk |Tk[l] | > 1
1
k
b2 
 2(p + 1)KC 2 D
K 2 (p + 1) 2
k
a.s. [PL− ] ∀ n; for each k ∈ {1, ..., K} and ∀ 1 ≤ l ≤ (p + 1).
k

b k = Op (c − ), (c − /D
b k ) ≥ 0 is stochastically bounded away
Now, since D
nK
nK
from 0. Thus, ∀ k, and for any given  > 0, ∃ δ(k, ) > 0 (independent of n)
b k ) ≤ δ(k, )} ≤ ∗ ∀ n, where ∗ = /{4K(p+1)} > 0.
such that: PL− {(cn− /D
K
k
e
Let δ(K,
) = min{δ(k, ) : k = 1, ..., K} > 0 (as K is fixed). Let A(k, )
e
b k ) ≤ δ(K,
denote the event: {(c − /D
)}, and let Ac (k, ) be its complement.
nK

e
b k ) > δ(K,
). Thus, the
Then, PL− {A(k, )} ≤ ∗ , while on Ac (k, ), (cn− /D
K
k
bound in (A.10) is dominated by: 2 exp[−M 2 δe2 (K, )/{2(p + 1)KC 2 }] on
Ac (k, ), and trivially by 2 on A(k, ) ∀ k. Plugging the bound of (A.10) into
(A.9) and using all these facts, we then have:



p+1
K X


M 2 c2n−


X
K


P kGn,K k > M cn− ≤
EL− 2 exp −
K
k
b2 
 2(p + 1)KC 2 D
k
k=1 l=1




2
2
p+1
K X


M cn−
X
K
1 c
+ 1A(k,) 
=
EL− 2 exp −
k
b 2  A (k,)
 2(p + 1)KC 2 D
k
k=1 l=1
"
)
(
#
p+1
K
XX
M 2 δe2 (K, )
≤
2 exp −
PL− {Ac (k, )} + 2 PL− {A(k, )}
k
k
2(p + 1)KC 2
k=1 l=1
)
#
"
(
M 2 δe2 (K, )
+ ∗
≤ 2K(p + 1) exp −
2(p + 1)KC 2


(A.11) ≤ + =  (with some suitable choice M for M ),
2 2
where the last step follows from noting the definition of ∗ and choosing M
to be any M large enough such that 4 exp[− M 2 δe2 (K, )/{2(p + 1)KC 2 }] ≤
/{K(p + 1)}. Thus, (A.11) shows Gn,K = Op (cn− ) for any fixed K ≥ 2. This
K
further establishes (3.11) and all its associated implications. The proof of
Theorem 3.2 is now complete.

SEMI-SUPERVISED LINEAR REGRESSION

31

ACKNOWLEDGEMENTS
The authors would like to thank Dr. James Robins and Dr. Eric Tchetgen
Tchetgen for many helpful discussions throughout the progress of this paper,
as well as the editor Dr. Edward George, the anonymous associate editor
and the two referees for their useful comments and suggestions that helped
significantly in improving and revising the original version of this article.
SUPPLEMENTARY MATERIAL
Supplement to “Efficient and Adaptive Linear Regression in
Semi-Supervised Settings”
(doi: COMPLETED BY THE TYPESETTER; .pdf). The supplement includes: (i) Supplementary results for the simulation studies and the real data
analysis; (ii) Discussions on generalization of the proposed SS estimators to
MAR settings; (iii) Proof of Lemma A.1; (iv) Proof of Theorem 3.1; (v)
Proof of Theorem 4.1; and (vi) Proofs of Lemmas A.2-A.3 and Theorem 4.2.

SEMI-SUPERVISED LINEAR REGRESSION

i

SUPPLEMENT TO “EFFICIENT AND ADAPTIVE LINEAR
REGRESSION IN SEMI-SUPERVISED SETTINGS”
BY ABHISHEK CHAKRABORTTY AND TIANXI CAI
University of Pennsylvania and Harvard University
This supplementary document contains additional numerical results as well
as technical materials, including proofs etc., that could not be accommodated in the main article.
I. Numerical Studies: Supplementary Results.
I.1. Simulation Results for p = 2. For p = 2, we investigated three choices
of m(x) as follows:
(Linear): m(x) = x1 + x2 ;
(NL-Iλ(k) ): m(x) = x1 + x2 + λ(k) x1 x2 for λ(1) = 0.5 and λ(2) = 1; and
(NL-Qγ (k) ): m(x) = x1 + x2 + γ (k) (x21 + x22 ) for γ (1) = 0.3 and γ (2) = 1.
Since the dimension is low, we implemented EASE using the KS and KM
smoothers with P2 = I2 for both, i.e. without any dimension reduction.
For comparison, the other SS estimators were also obtained. In Table I,
we summarize the efficiencies of all the estimators relative to OLS, based
e the
on the empirical mean squared error (MSE), where for any estimator θ,
2
e 0 k averaged over the 500 replications.
empirical MSE is summarized as kθ−θ
Models
Linear
NL-Iλ(1)
NL-Iλ(2)
NL-Qγ (1)
NL-Qγ (2)

OLS
(Ref.)
1
1
1
1
1

SNP EASE
(T := KS)
0.897
0.995
1.229
1.243
2.261
2.261
2.241
2.215
4.096
4.144

SNP EASE
(T := KM)
0.920
0.988
1.338
1.355
2.301
2.267
2.500
2.550
4.612
4.641

Other SS Estimators
DRESS1 DRESS3 MSSL
0.993
0.963
0.993
1.072
1.039
1.072
1.217
1.181
1.216
1.187
2.063
1.187
1.352
3.217
1.352

Table I: Efficiencies of SNP and EASE, obtained using T := KS or KM, as
well as DRESS1 , DRESS3 and MSSL, relative to OLS with respect to the
empirical MSE under the various models considered with p = 2.
For this setting, all estimators have comparable efficiency under the linear
model, as expected. Under the non-linear models, the EASE estimators are
substantially more efficient than the OLS and also more efficient than the
other SS estimators. For the non-linear models with quadratic effects, the
DRESS3 is also substantially more efficient than the OLS while our EASE
estimator performs even better. For the non-linear models with interaction

ii

A. CHAKRABORTTY AND T. CAI

effects, the efficiency gain was very modest when employing existing SS
estimation procedures while it was quite substantial for EASE.
I.2. Supplementary Results for the Data Example. We present in Table
II some summary measures of the distributions in the labeled and unlabeled
data for each of the predictors in the data example, and also report p-values
for diagnostic tests aimed at detecting any possible differences in the labeled
and unlabeled data distributions for each of the predictors.
Predictors
Age
Gender
Race
Lupus
PmR
RA
SpA
Other ADs
Erosion
Seropositivity
Anti-CCPprior
Anti-CCPmiss
RF
RFmiss
Azathioprine
Enbrel
Gold salts
Humira
Infliximab
Leflunomide
Methotrexate
Plaquenil
Sulfasalazine
Other meds.
(Intercept)

Labeled Data
Mean
Sd
4.090 0.241
0.786 0.411
0.696 0.461
0.230 0.520
0.057 0.336
4.171 1.071
0.073 0.343
0.251 0.642
0.577 0.495
0.369 0.483
0.386 0.629
0.645 0.479
0.949 0.772
0.307 0.462
0.121 0.397
0.738 0.858
0.346 0.568
0.856 0.833
0.386 0.673
0.549 0.740
1.417 0.638
0.248 0.464
0.535 0.752
0.163 0.380
–
–

Unlabeled Data
Mean
Sd
4.070
0.272
0.799
0.401
0.673
0.469
0.251
0.600
0.078
0.382
4.084
1.079
0.066
0.313
0.271
0.690
0.567
0.494
0.395
0.489
0.405
0.718
0.610
0.488
0.897
0.846
0.324
0.468
0.137
0.419
0.722
0.822
0.336
0.555
0.917
0.835
0.400
0.672
0.555
0.743
1.389
0.669
0.273
0.496
0.554
0.734
0.189
0.403
–
–

P-values from Diagnostic Tests
T-test Wilcoxon Test PS Model
0.151
0.373
0.475
0.556
0.547
0.574
0.371
0.378
0.456
0.461
0.877
0.689
0.269
0.326
0.255
0.144
0.055
0.255
0.716
0.932
0.780
0.582
0.780
0.540
0.709
0.701
0.743
0.331
0.335
0.231
0.592
0.471
0.564
0.192
0.198
0.089
0.227
0.184
0.199
0.516
0.520
0.410
0.463
0.449
0.734
0.740
0.921
0.574
0.729
0.791
0.962
0.193
0.190
0.177
0.711
0.607
0.789
0.895
0.912
0.973
0.435
0.584
0.549
0.331
0.463
0.398
0.661
0.458
0.724
0.223
0.192
0.322
–
–
0.000

Table II: Comparison of the means and standard deviations (sd) from the
labeled and unlabeled data for each predictor in the data example. Shown
also are the p-values obtained from various diagnostic tests, testing for possible differences in the distributions of each of the predictors in the labeled
and unlabeled data, including a two-sample T-test (with possibly unequal
variances in the two populations), a Wilcoxon rank sum test, and a test
obtained by fitting a parametric logistic regression model for the propensity
score (PS) of missingness, with all the predictors included as covariates, and
then testing for the null effect of each of the predictors in the fitted model.

iii

SEMI-SUPERVISED LINEAR REGRESSION

I.3. Simulation Results on the Prediction Error. In Table III, we present
the out-of-sample mean squared prediction error for various SNP imputation
estimators under all the models considered in the simulation studies for
p = 2, 10 and 20. The results suggest that the SNP imputation estimators
for both KS and KM based smoothers perform substantially better than the
naive non-parametric estimator based on a p-dimensional kernel smoothing.
(a) p = 2
Models
Linear
NL-Iλ(1)
NL-Iλ(2)
NL-Qγ (1)
NL-Qγ (2)

T := KS
m
bT
µ
bT
0.36 0.35
0.34 0.33
0.29 0.28
0.29 0.29
0.20 0.20

T := KM
m
bT
µ
bT
0.35 0.35
0.32 0.32
0.27 0.27
0.27 0.27
0.17 0.18

(b) p = 10 and 20

Models
(I)

(II)

Linear
NL1C
NL2C
NL3C
Linear
NL1C
NL2C
NL3C

T := KS2,P2
m
bT
µ
bT
0.186 0.174
0.143 0.131
0.241 0.228
0.275 0.254
0.106 0.096
0.147 0.131
0.179 0.163
0.218 0.199

p = 10
T := KM
m
bT
µ
bT
0.214 0.204
0.147 0.154
0.150 0.156
0.138 0.146
0.133 0.126
0.135 0.143
0.135 0.142
0.137 0.144

KSp
0.260
0.290
0.554
0.458
0.313
0.658
0.487
0.527

T := KS2,P2
m
bT
µ
bT
0.143 0.099
0.337 0.332
0.557 0.566
0.543 0.552
0.095 0.054
0.341 0.336
0.423 0.422
0.424 0.424

p = 20
T := KM
m
bT
µ
bT
0.130 0.119
0.280 0.317
0.271 0.299
0.282 0.327
0.083 0.073
0.283 0.331
0.272 0.302
0.282 0.318

KSp
0.698
0.921
1.000
0.974
0.755
1.095
0.980
0.996

Table III: Mean squared prediction errors (PEs), relative to Var(Y ), of the
b r ) and µ
br,K ), with T :=
SNP imputation estimators m
b T = m(x;
b
P
bT = µ
b(x; P
KS2,P2 or T := KM, under the various models discussed in the simulation
studies. Shown also are the corresponding PEs for the fully non-parametric
KS estimator, KSp , for comparison, in the case of p = 10 and 20.
II. Generalization to the Missing at Random (MAR) Case. Our
SS estimation methods proposed so far assume that the underlying Y for
subjects in U are MCAR, a standard (and often implicit) assumption made
in SSL. In this section, we provide some discussions on possible generalizations of our SS methods to the MAR case. Such generalizations might be
desirable for settings where the availability of Y is not determined by design. To this end, let N̄ = N + n denote the sample size of the entire data
S = L ∪ U. Then, S = {Zi ≡ (Ri , Ri Yi , Xi ) : i = 1, . . . , N̄ } consists of N̄
i.i.d. realizations of Z = (R, RY, X), where R ∈ {0, 1} denotes the indicator

iv

A. CHAKRABORTTY AND T. CAI

of Y being observed. As opposed to the stronger MCAR setting with the
assumption R ⊥
⊥ (Y, X) and the probability law of S being determined by
the law of (Y, X), we now have: under the MAR setting, R ⊥
⊥ Y | X and the
probability law of S is determined by PZ , the law of Z. For notational ease,
we also let PN̄ denote the empirical measure for S, and for any function e(·)
P
of Z, possibly random andR vector-valued, we let PN̄ (e) = N̄ −1 N̄
i=1 e(Zi ),
and PZ (e) = EZ {e(Z)} = e(z)dPZ (z).
P
Under a SS set-up as above, we have: n = N̄
i=1 Ri is a random quantity
and n/N̄ → 0 in probability. It is important to note that πN̄ ≡ P(R = 1)
must depend on N̄ . Let πN̄ (X) = P(R = 1 | X) be the “propensity score”,
assumed to be strictly greater than 0 almost surely (a.s.) for any given N̄
1
−2
and let bN̄ = [E{πN̄
(X)}]− 2 . Then, under the above set-up, we assume that
E{πN̄ (X)} = πN̄ → 0, bN̄ → 0 and N̄ bN̄ → ∞ as N̄ → ∞. This decaying
sampling probability is the main factor that distinguishes SSL from standard
missing data problems and contributes to the complexity of devising and
analyzing SS estimators which, even for the MCAR setting, were seen to only
1
1
have a convergence rate of (N̄ bN̄ )− 2 , rather than N̄ − 2 , with bN̄ = (n/N̄ ).
For simplicity, we shall first assume that πN̄ (X) is known and next detail
how we may extend our proposed procedures to obtain SS estimators of θ 0 ,
→
−
→
−
the solution to: φ(θ) ≡ E{ X(Y − X 0 θ)} = 0, under the MAR setting.
To derive efficient SS estimators of θ 0 based on S under MAR, we first
b(P ,K) in Section 3.2 in fact remains
note that our proposed SNP estimator θ
r
valid even under the MAR setting, whenever the imputation is sufficient
i.e. µ(X; Pr ) equals the true conditional mean m(X) ≡ E(Y |X). For the
general case allowing for insufficient imputation, we need to modify our
SNP imputation to account for the MAR setting. To this end, we note that


→
−
→
−0
−
→
−0
R →
φ(θ) ≡ E{ X(Y − X θ)} = E
(1)
X(Y − X θ) ,
πN̄ (X)
→
−
→
−
(2)
= E[ X{m(X) − X 0 θ}],


→
−
→
−0
−
R →
(3)
= E[ X{m(X) − X θ}] + E
X{Y − m(X)} .
πN̄ (X)
More generally, for any µ(·) ∈ L2 (PX ) satisfying:


−
→
−
R →
X{Y − µ(X)} ≡ E[ X{m(X) − µ(X)}] = 0,
(4)
E
πN̄ (X)
it is easy to see that the following representation of φ(·) holds under MAR:


→
−
→
−0
−
R →
(5)
φ(θ) = E[ X{µ(X) − X θ}] + E
X{Y − µ(X)} .
πN̄ (X)

v

SEMI-SUPERVISED LINEAR REGRESSION

Let µ
b(·) be any estimator of µ(·) based on S. Motivated by (5), we may
bMAR ≡ θ
bMAR, µ(·) , as the solution to:
then devise a SS estimator of θ 0 , θ
(6)

b (θ) ≡ N̄ −1
φ
N̄

N̄ 
X
→
−
→
−
X i {b
µ(Xi ) − X 0i θ} +
i=1


−
Ri →
X i {Yi − µ
b(Xi )} .
πN̄ (Xi )

P →
− →
−0
Then, letting ΓN̄ = N̄ −1 N̄
i=1 X i X i , it is straightforward to show that:


bMAR − θ 0 = P (T) + P (S) − P (b
(7)
ΓN̄ θ
N̄
N̄
N̄ e), where
→
−
→
−
→
−
where T(Z) = X{µ(X) − X 0 θ 0 }, S(Z) = {R/πN̄ (X)} X{Y − µ(X)}, and


→
−
R
b
− 1 X{b
µ(X) − µ(X)}.
e(Z) ≡
πN̄ (X)
Convergence rates of the terms in (7) need careful analysis as the asymp1
totics here is non-standard, with the dominating rate being slower than N − 2 .
To this end, note that PN̄ (T) is a simple centered i.i.d. average of vari1
ables with bounded variance. Hence, TN̄ = O(N − 2 ) indeed. On the other
hand, PN̄ (S) has a slower convergence rate since the variance of S, VN̄ ,
diverges due to the πN̄ (·) ↓ 0 appearing in the denominator. Under mild
moment conditions, it can be shown that bN̄ VN̄ converges to a positive definite matrix V with kVk < ∞. Hence, using concentration inequalities, and
assuming N̄ bN̄ → ∞, it can be shown that the convergence rate of PN̄ (S) is
1
O{(N bN̄ )− 2 }. Further, using CLT for triangular arrays, it can be shown un1

d

der suitable conditions that (N bN̄ ) 2 PN̄ (S) → N(p+1) [0, V]. Lastly, to control
the term PN̄ (b
e), note that PZ (b
e) = 0. Therefore, GN̄ is a centered empirical
P

process indexed by µ
b(·) − µ(·). Hence, as long as EX [{b
µ(X) − µ(X)}2 ] → 0,
and µ
b(·) − µ(·) lies in a P−Donsker class with probability → 1, it can be
shown using results from empirical process theory (see Van der Vaart (2000)
1
for instance) that (N bN̄ ) 2 PN̄ (b
e) = op (1). Finally, note that ΓN̄  0 a.s., and
−1
−1
− 21
ΓN̄ = Γ +Op (N ). Hence, under suitable regularity conditions, we have:
1
bMAR − θ 0 )
(N̄ bN̄ ) 2 (θ

=

1

(N̄ bN̄ ) 2 Γ−1

N̄
1 X
S(Zi ) + op (1)
N̄ i=1

d

→ N(p+1) [0, Γ−1 VΓ−1 ], with V as defined above.
Having now provided an abstract sketch of the construction of the estimators and their properties, we next briefly discuss the choice of the

vi

A. CHAKRABORTTY AND T. CAI

‘imputation’ function µ(·), and its estimator µ
b(·) inherent in the construcbMAR ≡ θ
bMAR, µ(·) . With (r, Pr , P
b r ) as defined in Section 3.2 and
tion of θ
{K(·), h, Kh (·, ·)} as in Section 4, we may modify the SNP estimator in Section 3.2 under the MAR setting as follows: consider µ(X) ≡ µ(X; Pr ) =
→
−
m(X; Pr ) + X 0 η Pr , where m(X; Pr ) = E(Y |P0r X), and η Pr satisfies


−
→
−0
R →
E
X{Y − m(X; Pr ) − X η Pr } = 0.
πN̄ (X)
b r)
This will ensure that (4) holds. Then, we may estimate µ(X; Pr ) as µ
b(X; P
→
−
b r ) + X 0η
b Pr , where
= m(X;
b
P
PN̄
b r) =
m(x;
b
P

Ri
b0
b0
i=1 πN̄ (Xi ) Yi Kh (Pr Xi , Pr x)
,
PN̄
Ri
b0
b0
i=1 πN̄ (Xi ) Kh (Pr Xi , Pr x)

b Pr satisfies: N̄ −1
η

N̄
X
i=1

and

−
−0
Ri →
b r) − →
b Pr } = 0.
X i {Yi − m(X
b i; P
X iη
πN̄ (Xi )

Thus, to accommodate the MAR setting, one essentially needs to implement
appropriately weighted versions of both the smoothing and the refitting steps
in our original SNP imputation. Of course, while we have chosen the smoothing method T to be the weighted KS here for illustration, other reasonable
choices of T such as an appropriately weighted KM may also be used. Under
bMAR
MCAR, with πN̄ (X) ≡ πN̄ ≡ n/N̄ and bN̄ = (n/N̄ ), the estimator θ
indeed becomes (asymptotically) equivalent to the SNP estimators obtained
earlier in Section 3.2. Further, with various choices of µ(·), the SNP imputation strategy again equips us with a family of SS estimators of θ 0 under
the MAR setting, with µ(·) = m(·) leading to the optimal estimator.
bMAR is derived with a known π (·), for simplicity.
The above estimator θ
N̄
In practice, πN̄ (·) is typically unknown and a consistent estimator π
bN̄ (·)
b
may be constructed. Then, one may modify θ MAR by replacing πN̄ (·) with
π
bN̄ (·) in all the steps. The resulting estimator will have an expansion similar
to (7) but with extra error terms accounting for the variability in π
bN̄ (·),
which need to be properly controlled. The theoretical analysis will be more
involved since establishing the convergence rates and asymptotic expansion
for π
bN̄ (·) − πN̄ (·) is also non-standard due to bN̄ → 0. Lastly, the score
bMAR has the additional benefit of ‘double
equation (5) used to construct θ
robustness’, in the sense that even if π
bN̄ (·) is inconsistent for πN̄ (·), as long
bMAR, µ(·) is consistent for θ 0 . On the other
as µ
b(·) estimates the true m(·), θ
hand, as long as π
bN̄ (·) targets the true πN̄ (·), then for any choice of µ(·),

SEMI-SUPERVISED LINEAR REGRESSION

vii

bMAR, µ(·) is consistent for θ 0 . For the MCAR case, π
θ
bN̄ (·) ≡ π
bN̄ = n/N̄
is always consistent for πN̄ (·) and in fact this is exactly what allowed us
to achieve a family of SNP estimators, all consistent for θ 0 , for various
choices of the SNP imputation function µ(·). For a specific choice of µ(·), this
b Note however that under MAR,
also included the supervised estimator θ.
b
θ simply denotes the naive ‘complete case’ estimator and in general, is not
consistent for θ 0 unless πN̄ (·) is trivially a constant or m(·) is exactly linear.
III. Proof of Lemma A.1. Firstly, since d is fixed, it suffices to prove
(j)
b (j)
b
the result for any arbitrary scalar coordinate G
n,m ≡ Gn,m (say) and Gn ≡
b n,m and Gn respectively, for any j ∈ {1, . . . , d}. For any data
G n (say) of G
∗
S and S , we let PS and PS,S∗ denote the joint probability distributions of
the observations in S and (S, S∗ ) respectively, ES (·) denote the expectation
w.r.t PS , and PS|S∗ denote the conditional probability distribution of the
observations in S given S∗ .
1
To show that Gbn,m − G n = Op (m− 2 ), we first note that since Sn ⊥
⊥ Sm ,


n
o

1
1
PSn ,Sm |Gbn,m − G n | > m− 2 t = ESn PSm |Gbn,m − G n | > m− 2 t Sn ,
b n,m − Gn is a centered average of
for any t > 0. Now, conditional on Sn , G
m
{b
gn (Zj )}j=1 which are i.i.d. and bounded by Tbn < ∞ a.s. [PSn ] ∀ n. Hence,
applying Hoeffding’s inequality, we have for any n and m,
!


2 t2
1
2m
a.s. [PSn ].
(8) PSm |Gbn,m − G n | > m− 2 t Sn ≤ 2 exp −
4m2 Tbn2
Now, since Tbn ≥ 0 is Op (1), we have: for any given  > 0, ∃ δ() > 0 such
that: PSn {Tbn > δ()} ≤ /4 ∀ n. Let A() denote the event: {Tbn > δ()} and
let Ac () denote its complement. Then, using (8), we have: ∀ n and m,
(
!)


2 t2
1
2m
−
PSn ,Sm |Gbn,m − G n | > m 2 t ≤ ESn 2 exp −
4m2 Tbn2
(
!)
"
!
#

t2
t2
= ESn 2 exp −
= ESn 2 exp −
1Ac () + 1A()
2Tbn2
2Tbn2




t2
c
≤ 2 exp − 2
PSn {A ()} + 2 PSn {A()}
2δ ()


t2



≤ 2 exp − 2
+
≤ + =  (for some suitable choice of t),
2δ ()
2
2 2

viii

A. CHAKRABORTTY AND T. CAI

where the last step follows by choosing t ≡ t to be any large enough t
such that exp{−t2 /2δ 2 ()} ≤ /4. Such a choice of t clearly exists. This
establishes the first claim (a) in Lemma A.1. The second claim (b) in Lemma
A.1 is a trivial consequence of the Central Limit Theorem (CLT).
IV. Proof of Theorem 3.1. To show Theorem 3.1, we first note that
1
1
under Assumption 3.1 (i)-(v), and letting an = (log n) 2 (nhp )− 2 + hq , the
following holds:
(9)

supx∈X |m(x)
b
− m(x)| = Op (an ) = supx∈X |fb(x) − f (x)|.

(9) is a fairly standard result and we only provide a sketch of its proof as
follows. Under Assumption 3.1 (ii)-(iii), using Theorem 2 of Hansen (2008),
supx∈X |b
l(x) − EL {b
l(x)}| = Op (a∗n ) = supx∈X |fb(x) − EL {fb(x)}|, where a∗n =
1
1
(log n) 2 (nhp )− 2 . Next, using standard arguments based on Taylor series
expansions of l(·) and m(·) under their assumed smoothness, and noting
that K(·) is a q th order kernel having finite q th moments, we obtain:
supx∈X |EL {b
l(x)} − l(x)| = O(hq ) = supx∈X |EL {fb(x)} − f (x)|.
Combining these two results, and the definitions of m(.) and m(.)
b along with
Assumption 3.1 (iv), we have (9). Next, note that using (3.2), we have:
bnp − θ 0 ) = EU [N −1
ΓN (θ

n+N
X

1
→
−
→
−
X j {m(X
b j ) − X 0j θ 0 }] + Op (N − 2 )

j=n+1
1
→
−
= EX [ X{m(X)
b
− m(X)}] + Op (N − 2 ),

−
→
−
where the first step is due to Lemma A.1 (a) with supx∈X k→
x {m(x)−
b
x 0 θ 0 }k
→
−
→
−
0
b
− m(x)| + |m(x) − x θ 0 |}] = Op (1) due to (9) and
≤ supx∈X [k x k{|m(x)
→
−
the boundedness of X and m(·), while the last step uses: EX [ X{m(X) −
→
−0
X θ 0 }] = 0 which follows from the definitions of θ 0 and m(·). It then follows
1
−1
further, using Γ−1
+ Op (N − 2 ), that
N =Γ
 n 1
1
−
2
bnp − θ 0 ) = n 12 Γ−1 EX [→
n 2 (θ
X{m(X)
b
− m(X)}] + Op
.
N
P
Letting φn (X) = (nhp )−1 ni=1 K{(X − Xi )/h}{Yi − m(X)}, and expanding
the first term in the above equation, we now obtain:
(10)





 1
1
bnp − θ 0 = Γ−1 T(1) + T(2) + Op n 2 ,
n2 θ
n,1
n,1
N

SEMI-SUPERVISED LINEAR REGRESSION

ix

1
→
−
(1)
where Tn,1 = n 2 EX { Xφn (X)/f (X)} and

h→
i
1
−
(2)
Tn,1 = n 2 EX Xφn (X){fb(X)−1 − f (X)−1 }
1
→
−
= n 2 EX [ X {m(X)
b
− m(X)} {f (X) − fb(X)}/f (X)]
n
o
 1 
1
−
(11) ≤ n 2 supx∈X k→
x k |m(x)
b
− m(x)| fb(x)/f (x) − 1
= Op n 2 a2n ,
where the last step in (11) follows from (9), Assumption 3.1 (iv) and the
(1)
boundedness of X. For Tn,1 , we have:
(1)

1

Z

1

n
X

Tn,1 = n 2
(12)

1
→
−
x φn (x)dx = n− 2

X

= n2

n Z
X
i=1

n−1

Z
Ai,n

i=1

→
−
x h−p Kh (x − Xi ) {Yi − m(x)} dx

X

−−−−−−−→
(Xi + hψ i ) K (ψ i ) {Yi − m(Xi + hψ i )} dψ i ,

where ψ i = (x − Xi )/h and Ai,n = {ψ i ∈ Rp : (Xi + hψ i ) ∈ X }. Now, since
K(·) is zero outside the bounded set K, the ith integral in (12) only runs over
(Ai,n ∩ K). Further, since h = o(1), using Assumption 3.1 (vi), Ai,n ⊇ K a.s.
[PL ] or, (Ai,n ∩ K) = K a.s. [PL ] ∀ 1 ≤ i ≤ n with n large enough. Thus, for
large enough n, (12) can be written as:
(1)
Tn,1

(13)
(14)

=n

− 21

1

= n2

n Z
X
−−−−−−−→
(Xi + hψ i ) K (ψ i ) {Yi − m(Xi + hψ i )} dψ i a.s. [PL ]

i=1
n
X

K

h→
i
−
n−1 X i {Yi − m(Xi )} + Op (hq )

i=1
n
X

1

= n− 2

 1 
→
−
X i {Yi − m(Xi )} + Op n 2 hq ,

i=1

where (13), and hence (14), follows from standard arguments based on Taylor
series expansions of m(Xi + hψ i ) around m(Xi ) under the assumed smoothness of m(·), and using the fact that K(·) is a q th order kernel. Combining
1
1
(10), (11) and (14), and noting that under our assumptions, (n 2 a2n + n 2 hq )
1
1
= O{n 2 hq + (log n)(n 2 hp )−1 }, the result of Theorem 3.1 now follows.

x

A. CHAKRABORTTY AND T. CAI
1

1

V. Proof of Theorem 4.1. Let an,2 = (log n) 2 (nhr )− 2 + hq . Then, we
first note that
(%)

(%)

supw∈XPr |ϕ
ePr (w) − ϕPr (w)| = Op (an,2 ),

(15)

∀ % ∈ {0, 1}.

To see this, note that under Assumption 4.1 (ii)-(iii), Theorem 2 of Hansen
1
1
(2008) applies, and we have for dn = (log n) 2 (nhr )− 2 ,
(%)

(%)

supw∈XPr |ϕ
ePr (w) − EL {ϕ
ePr (w)}| = Op (dn )

∀ % ∈ {0, 1}.

Next, using standard arguments based on a q th order Taylor series expansion
(%)
of ϕPr (·) and noting that K(·) is a q th order kernel, we obtain:
(%)

(%)

ePr (w)} − ϕPr (w)| = O(hq )
supw∈XPr |EL {ϕ

∀ % ∈ {0, 1}.

Combining these two results gives (15). Further,
supx∈X |m(x;
e
Pr ) − m(x; Pr )| = supw∈XPr |m
e Pr (w) − mPr (w)|
(
)
e
lPr (w) − lPr (w)
|lPr (w)| |lPr (w)|
≤ supw∈XPr
+ supw∈XPr
−
fPr (w)
feP (w)
feP (w)
r

r

(16) = Op (an,2 ),
where the last step follows from repeated use of (15) and Assumption 4.1
b r )− ϕ
(iii)-(iv). Next, we aim to bound supx∈X |ϕ
b(%) (x; P
e(%) (x; Pr )| to account
b r . Using a first order Taylor series
for the potential estimation error of P
expansion of K(.) under Assumption 4.1 (vi), we have: ∀ % ∈ {0, 1},


n
X
x
−
X
1
i
0
0
b r − Pr )
b r) − ϕ
∇K (wi,x )(P
Yi%
ϕ
b (x; P
e (x; Pr ) =
nhr
h
i=1
o
o
n
n
0 c (1)
b r − Pr )0 M
c (2)
b
(17)
= trace (Pr − Pr ) Mn,%,x + trace (P
n,%,x ,
(%)

(%)

where


 0
n
0
1 X
0 Pr x − Pr X i
(x − Xi ) ∇K
Yi% and
nhr+1
h
i=1


 0
n
0
X
1
0
0 Pr x − Pr Xi
=
(x − Xi ) ∇K (wi,x ) − ∇K
Yi% ,
r+1
nh
h

c (1)
M
n,%,x =
c (2)
M
n,%,x

i=1

with wi,x ∈ Rr being ‘intermediate’ points satisfying: kwi,x −P0r (x−Xi )h−1 k
b 0 (x−Xi )h−1 −P0 (x−Xi )h−1 k ≤ Op (αn h−1 ). The last bound, based on
≤ kP
r
r

xi

SEMI-SUPERVISED LINEAR REGRESSION

b r − Pr ) = Op (αn ) and the compactness of X , is uniform in (i, x). For any
(P
matrix A = [aij ], let kAkmax denote the max-norm of A, and |A| denote the
matrix [|aij |]. Now, Assumption 4.1 (viii) implies: k∇K(w1 ) − ∇K(w2 )k ≤
Bkw1 −w2 k ∀ w1 , w2 ∈ Rr , for some constant B < ∞. Then using the above
c (2)
arguments, we note that ∀ % ∈ {0, 1}, k supx∈X |M
n,%,x |kmax is bounded by:
)
(
n
B X
P0r x − P0r Xi
supx∈X
|Yi% |
kx − Xi k wi,x −
nhr+1
h
i=1
(
)
n
b r − Pr )0 (x − Xi )
B X
(P
%
≤ supx∈X
|Yi |
kx − Xi k
nhr+1
h
i=1

≤

sup

n
o
b r − Pr )0 (x − X)k
kx − Xkk(P

x∈X ,X∈X

n
 α 
B X %
n
|Y
|
≤
O
.
p
i
nhr+2
hr+2
i=1

The first two steps above use the triangle inequality, the Lipschitz continuity
of ∇K(·) and the definition of wi,x , while the next two use the compactness
of X , the uniform bound obtained in the last paragraph, the Law of Large
b r − Pr ) = Op (αn ). Thus, we have:
Numbers (LLN), and that (P
 2 
o
n
(2)
0
b r − Pr ) M
c n,%,x = Op αn
(18) supx∈X trace (P
∀ % ∈ {0, 1}.
hr+2
c (1)
c (1)
c (1,1) c (1,2)
Now for bounding M
n,%,x , let us first write it as: Mn,%,x = Mn,%,x − Mn,%,x ,
P
(1,2)
c n,%,x
c (1,1)
=
where M
= (nhr+1 )−1 ni=1 x∇K 0 {P0r (x − Xi )/h}Yi% and M
n,%,x
P
%
n
0
0
r+1
−1
(nh )
i=1 Xi ∇K {Pr (x − Xi )/h}Yi ∀ % ∈ {0, 1}. Then, under Assumption 4.1 (iii), (vi) and (vii), using Theorem 2 of Hansen (2008) along
with the compactness of X , we have: for each s ∈ {1, 2} and % ∈ {0, 1},
(19)

supx∈X



c (1,s)
c (1,s)
M
n,%,x − EL Mn,%,x


max

≤ Op

log n
nhr+2

1
2

.

Now, ∀ % ∈ {0, 1}, let ν (%) (w) = E{Y % | XPr = w}fPr (w) and ξ (%) (w) =
E{XY % | XPr = w}fPr (w). Further, let {∇ν (%) (w)}r×1 and {∇ξ (%) (w)}p×r
denote their respective first order derivatives. Then, ∀ % ∈ {0, 1}, we have:


c (1,1)
supx∈X EL M
n,%,x
max
 0

Z
x
(%)
0 Pr x − w
dw
= supx∈X r+1 ν (w) ∇K
h
h
max
Z

(%)0
0
(20)
= O(1),
Pr x + hψ K(ψ)dψ
= supx∈X x ∇ν
max

xii

A. CHAKRABORTTY AND T. CAI

=
(21)

=



(1,2)
c n,%,x
supx∈X EL M
max
 0

Z
(%)
0 Pr x − w
−(r+1)
ξ (w) ∇K
dw
supx∈X h
h
Z

supx∈X D(x) ∇ξ (%) P0r x + hψ K(ψ)dψ

max

= O(1),

max

where, ∀ x ∈ X , D(x) denotes the p × p diagonal matrix: diag(x[1] , . . . , x[p] ).
In both (20) and (21), the first step follows from definition, the second from
standard arguments based on integration by parts (applied coordinate-wise)
and change of variable, while the last one is due to compactness of X and
a medley of the conditions in Assumption 4.1 namely, boundedness and
integrability of K(·) and ∇K(·), (iii) and (v) for (20) so that ∇ν (%) (·) is
bounded on XPr , and (ix) for (21). It now follows that for each % ∈ {0, 1},


c (1)
supx∈X EL M
(22)
= O(1).
n,%,x
max

1

1

Letting d∗n = (log n) 2 (nhr+2 )− 2 , we now have from (19) and (22):
o
n
b 0 − P0 )M
c (1)
(23) supx∈X trace (P
= Op (αn d∗n + αn ) ∀ % ∈ {0, 1}.
n,%,x
r
r
Applying (23) and (18) to (17) using the triangle inequality, we have ∀ %,
(
)
1
2
2
(log
n)
α
n
(%)
(%)
b r )− ϕ
+ αn
(24) supx∈X |ϕ
b (x; P
e (x; Pr )| = Op
.
1 + αn
hr+2
(nhr+2 ) 2
b r) = b
b r )/fb(x; P
b r) = ϕ
b r )/ϕ
b r ).
Finally, note that m(x;
b
P
l(x; P
b(1) (x; P
b(0) (x; P
Repeated use of (24), along with (16) and Assumption 4.1 (iii)-(iv), leads
to:
b r ) − m(x; Pr )
supx∈X m(x;
b
P
b r ) − m(x;
≤ supx∈X m(x;
b
P
e
Pr ) + supx∈X |m(x;
e
Pr ) − m(x; Pr )|
(
)
b
e
e
b r) − e
l(x; P
l(x; Pr )
l(x; Pr )
l(x; Pr )
≤ supx∈X
+
−
+ Op (an,2 )
b r)
b r)
fb(x; P
fe(x; Pr ) fb(x; P
(
)
1
αn2
(log n) 2
(25) ≤ Op
+ αn
+ Op (an,2 ) = Op (an,1 + an,2 ).
1 + αn
hr+2
(nhr+2 ) 2
The proof of Theorem 4.1 is now complete.

SEMI-SUPERVISED LINEAR REGRESSION

xiii

VI. Proofs of Lemmas A.2-A.3 and Theorem 4.2.
VI.1. Proof of Lemma A.2. First note that for each % ∈ {0, 1},
Z
n X
n
X
(n,%)
ϕ
e(%) (x; Pr )Pn (dx) = n−2
Hi1 ,i2
i1 =1 i2 =1
(n,%)

is a V-statistic, where Hi1 ,i2 = h−r λ(Xi1 )Yi%2 K{P0r (Xi1 − Xi2 )/h}. Using
the V-statistic projection result given in Lemma 8.4 of Newey and McFadden
(1994), it then follows that for each % ∈ {0, 1},
n
o
G∗n λ(·)[ϕ
e(%) (· ; Pr ) − EL {ϕ
e(%) (· ; Pr )}]
h
i
 1

1
1
(n,%)
(n,%)
= n− 2 Op E(kHi1 ,i1 k) + {E(kHi1 ,i2 k2 )} 2 = Op n− 2 h−r ,
(26)
% having finite
The last step follows from K(·) and λ(·)nbeing bounded and Y o
1
(%)
e? (· ; Pr )}] is a centered
2nd moments. Now, observe that n 2 G∗n λ(·)[EL {ϕ
sum of i.i.d. random vectors bounded by:
n
o
(%)
Dn,% = supx∈X kλ(x)k |EL {ϕ
e? (x; Pr )}| = O(hq ) ∀ % ∈ {0, 1},

˜ with population limit ξ(·), we use
where throughout, for any estimator ξ(·)
˜ − ξ(·).
the notation ξ˜? (·) to denote its centered version given by: ξ˜? (·) = ξ(·)
q
e? (x; Pr )}| =
Here, Dn,% = O(h ) since λ(·) is bounded and supx∈X |EL {ϕ
(%)
(%)
q
ePr (w)} − ϕPr (w)| = O(h ), as argued while proving (15).
supw∈XPr |EL {ϕ
Hence, ∃ a constant κ% > 0 such that hq /Dn,% ≥ κ% ∀ n. Then, using
Hoeffding’s Inequality, we have: ∀ n, given any  > 0 and any M = M ()
large enough,




d
n
o
X
M 2 h2q
M hq
(%)
∗
≤ 2d exp −
⇒
P Gn λ[l] (·)[EL {ϕ
e? (· ; Pr )}] > 1
2
2dDn,%
d2
l=1
!
h
n
o
i
M 2 κ2%
(%)
∗
q
e? (· ; Pr )}] > M h ≤ 2d exp −
≤ ⇒
P Gn λ(·)[EL {ϕ
2d
n
o
(%)
(27)
G∗n λ(·)[EL {ϕ
e? (· ; Pr )}] = Op (hq ) ∀ % ∈ {0, 1}.
Combining (26) and (27) using the linearity of G∗n (·), we then have (A.1).
Next, to show (A.2), let f (x; Pr ) = ϕ(0) (x; Pr ) and l(x; Pr ) = ϕ(1) (x; Pr ).
Then, we write
e (1) (·) − T
e (2) (·) − T
e (3) (·) + T
e (4) (·)}],
G∗n [λ(·){m
e ? (· ; Pr )}] = G∗n [λ(·){T
n,Pr
n,Pr
n,Pr
n,Pr

xiv

A. CHAKRABORTTY AND T. CAI

where
e
e (1) (x) = l? (x; Pr ) ,
T
n,Pr
f (x; Pr )
e
e
e (3) (x) = l? (x; Pr )f? (x; Pr ) , and
(28) T
n,Pr
fe(x; Pr )f (x; Pr )
(1)

e
e (2) (x) = f? (x; Pr )l(x; Pr ) ,
T
n,Pr
f (x; Pr )2
2
e
e (4) (x) = l(x; Pr )f? (x; Pr ) .
T
n,Pr
fe(x; Pr )f (x; Pr )2

(2)

Since λPr (x) ≡ λ(x)f (x; Pr )−1 and λPr (x) ≡ λ(x)l(x; Pr )f (x; Pr )−2 are
bounded a.s. [PX ] due to Assumption 4.1 (iii)-(iv) and the boundedness of
λ(·), using these as choices of ‘λ(·)’ in (A.1), we have:
(1)
e (1) (·)} = Op (b(1) ),
G∗n {λPr (·)e
l? (· ; Pr )} = G∗n {λ(·)T
n
n,Pr
(2)
e (2) (·)} = Op (b(1) ).
G∗n {λPr (·)fe? (· ; Pr )} = G∗n {λ(·)T
n
n,Pr

e (s) (x)k ≤ Op (a2 ) which follows
Further, for each s ∈ {3, 4}, supx∈X kT
n,2
n,Pr
from repeated use of (15) along with Assumption 4.1 (iii)-(iv). Consequently,
e (s) (·)} is bounded
with λ(·) bounded a.s. [PX ], for each s ∈ {3, 4}, G∗n {λ(·)T
n,Pr
1

by: Op (n 2 a2n,2 ). Combining all these results using the linearity of G∗n (·), we
1

(1)

1

e ? (· ; Pr )} = Op (bn + n 2 a2n,2 ) = Op (n 2 a2n,2 ), thus
finally obtain: G∗n {λ(·)m
leading to (A.2). The proof of the lemma is now complete.
VI.2. Proof of Lemma A.3. Throughout this proof, all additional notations introduced, if not explicitly defined, are understood to have been
adopted from the proof of Theorem 4.1 in Section V. Now, using (17),
0 c (2)
b r )−ϕ
b 0 −P0 )M
c (1)
b0
ϕ
b(%) (x; P
e(%) (x; Pr ) = trace{(P
n,%,x }+trace{(Pr −Pr )Mn,%,x },
r
r
c (1)
c (1,1) c (1,2)
and M
n,%,x = Mn,%,x − Mn,%,x , as defined in Section V. Thus,
n (1,1)
o
b
b(1,2) (·) + ζ
b(2) (·) ,
b r) − ϕ
G∗n [λ(·){ϕ
b(%) (· ; P
e(%) (· ; Pr )}] = G∗n ζ
(·)
−
ζ
n,%,λ
n,%,λ
n,%,λ
where ∀ (ω) ∈ {(1, 1), (1, 2), (2)}, % ∈ {0, 1}, and x ∈ X ,
n
o
b(ω) (x) = λ(x) trace (P
b 0 − P0 )M
c (ω)
(29)
ζ
n,%,x .
n,%,λ

r

r

Then, ∀ s ∈ {1, 2} and l ∈ {1, . . . , d}, each element of
Z
n X
n
X
(n,s)
−2
c (1,s)
Hl,% (i1 , i2 )
λ[l] (x)M
n,%,x Pn (dx) = n
i1 =1 i2 =1

is a V-statistic, where
(n,s)

Hl,% (i1 , i2 ) = h−(r+1) λ[l] (Xi1 )Yi%2 U(s) (i1 , i2 )∇K 0 {P0r (Xi1 − Xi2 )/h}

xv

SEMI-SUPERVISED LINEAR REGRESSION

with U(1) (i1 , i2 ) = Xi1 and U(2) (i1 , i2 ) = Xi2 . Hence, similar to the proof
of (26), using Lemma 8.4 of Newey and McFadden (1994) with X compact,
∇K(·) and λ(·) bounded, and Y % having finite 2nd moments, we have: for
each l ∈ {1, . . . , d}, s ∈ {1, 2} and % ∈ {0, 1},
 1
h
n
oi

− 2 −(r+1)
c (1,s) − EL λ[l] (·)M
c (1,s)
=
O
n
G∗n λ[l] (·)M
h
.
p
n,%,(·)
n,%,(·)
max

b r − Pr ) = Op (αn ) that for each s and %,
It then follows from (P
h (1,s)
n (1,s) oi


1
b
b
(30)
G∗ ζ
(·) − EL ζ
(·) = Op αn n− 2 h−(r+1) .
n

n,%,λ

n,%,λ

1
c (1,s) }] is
Next, for any given l, s and %, each element of n 2 G∗n [EL {λ[l] (·)M
n,%,(·)
a centered sum of i.i.d. random variables which are bounded by:
n
o
c (1,s)
supx∈X kλ(x)k |EL (M
= O(1),
n,%,x )|

max

where the order follows from (20), (21) and the boundedness of λ(·). Hence,
b r −Pr ) =
similar to the proof of (27), using Hoeffding’s inequality and that (P
Op (αn ), we have: ∀ l ∈ {1, . . . , d}, s ∈ {1, 2} and % ∈ {0, 1},
(31)

G∗n

h

n
oi
c (1,s)
EL λ[l] (·)M
n,%,(·)

max

= Op (1) ⇒

G∗n

 

(1,s)
b
EL ζ n,%,λ (·)
= Op (αn ).

For any matrix A, let us denote by A[a,b] the (a, b)th element of A. Now,
b(2) (.)} in (29), note that kG∗ {ζ
b(2) (·)}k is bounded by:
to control G∗ {ζ
n

n,%,λ

n

1

n 2 supx∈X kλ(x)k

XZ

n,%,λ



c (2)
b 0 − P0 )[b,a] M
(P
n,%,x
r
r

a,b
1
2

≤ n rp

b r − Pr
{kλ(x)k kx − Xk} P

sup

x∈X ,X∈X

(32)

[a,b]

max

(Pn + PX )(dx)

b %∗
Z
n




1
b %∗ ,
≤ Op n 2 αn Z
n

b r − Pr ) = Op (αn ) and the boundedness
where the last step follows from (P
R (%)
%∗
bn = Z
b n (x) (Pn + PX )(dx) with
of X and λ(·), and Z
b (%) (x) = n−1
Z
n

 0

n
X
|Yi% |
Pr (x − Xi )
∇K(wi,x ) − ∇K
.
hr+1
h
i=1

b r − Pr )0 (x − Xi )h−1 k ≤ Op (αn h−1 )
Now, kwi,x − P0r (x − Xi )h−1 k ≤ k(P
uniformly in (i, x), as noted while proving (18). Further, with L∗ , as defined

xvi

A. CHAKRABORTTY AND T. CAI

b r −Pr )0 (x−Xi )h−1 k ≤
in Assumption 4.1 (vii), let An denote the event: {k(P
∗
b
L ∀ x ∈ X , i = 1, .., n}. Then, with (Pr − Pr ) = Op (αn ), X compact and
1
αn h−1 = o(1) since n 2 αn2 h−2 = o(1) as assumed, it follows that P(An ) → 1.
Using these along with Assumption 4.1 (vii) and the function φ(.) defined
therein, we have: on An with P(An ) → 1,

 0
n
X
b r − Pr )0 (x − Xi )
|Yi% | (P
Pr (x − Xi )
φ
nhr+1
h
h
i=1

 0
n
%
X
|Y
|
P
(x
−
X
)
√
i
r
i
b r − Pr
≤ rp sup kx − Xk P
.
φ
nhr+2
h
max
x∈X ,X∈X

b (%) (x) ≤
Z
n

i=1





b %∗
e %∗
e %∗
Thus, Z
n ≤ Op αn Zn , where Zn =
e (%) (x)
Z
n

=n

−1

n
X

e (%) (x; Zi ),
Z
n

and

R

e (%)
Z
n (x)(Pn + PX )(dx),

e (%) (x; Z)
Z
n

i=1

|Y % |
= r+2 φ
h



P0r (x − X)
h


.

0

Let Z0 ≡ (Y 0 , X0 )0 ∼ PZ be generated independent of L, and define:
e (1) = n−1
U
n,%

n
X

e (%) (X0 ; Zi )},
EX0 {Z
n

e (2) = n−1
U
n,%

e (%) (Xi ; Z0 )},
EZ0 {Z
n

i=1

i=1

e (1,1) = E{Z
e (%) (X0 ; Z0 )},
U
n,%
n
Then, first note that:
Z

n
X

R

and

e (k) = E{Z
e (%) (X0 ; Z)k } for k = 1, 2.
V
n,%
n

e (%)
e (1)
Z
n (x)PX (dx) = Un,% . Further, since

e (%) (x)Pn (dx) = n−2
Z
n

n X
n
X

e (%) (Xi ; Zi )
Z
n
1
2

i1 =1 i2 =1

is a V-statistic, we have:
Z
e (%) (x)Pn (dx) = U
e (1) + U
e (2) − V
e (1) + Op {n−1 U
e (1,1) + n−1 (V
e (2) ) 21 }
Z
n
n,%
n,%
n,%
n,%
n,%
using Lemma 8.4 of Newey and McFadden (1994). Then, with all notations
as above, we have:


e (1,1) + n−1 (V
e (2) ) 21 ≤ Op n−1 h−(r+2) ,
(33)
n−1 U
and
n,%
n,%

SEMI-SUPERVISED LINEAR REGRESSION

xvii



Z
n
1 X %
w − P0r Xi
fPr (w)dw
φ
|Yi |
nhr+2
h
XPr
i=1
(
)
Z
n
BPr X
%
≤
|Yi |
φ(ψ i )dψ i ,
nh2
An
Xi
i=1
)
Z
(
n
X

BPr
%
−1
φ(ψ)dψ
n
|Yi | ≤ Op h−2 ,
≤ 2
h
Rr

e (1) =
U
n,%

(34)

i=1

where ψ i = h−1 (w − P0r Xi ) ∀ i, Anx = {ψ : (P0r x + hψ) ∈ XPr } ∀ x ∈ X ,
and BPr = supw∈XPr fPr (w) < ∞. The error rate in (33) follows since φ(·)
e (1)
is bounded and Y % has finite 2nd moments, while that of U
n,% follows from
Assumption 4.1 (iii), integrability of φ(·), and LLN applied to the sequence
e (2)
e (1)
{Yi% }ni=1 having finite 2nd moments. Now, note that U
n,% − Vn,% is a centered
0
n
e (%)
average of [EZ0 {Z
n (Xi ; Z )}]i=1 which are i.i.d. and bounded by:
 0

Z
Pr x − w
1
(%)
(%)
e
φ
mPr (w)fPr (w)dw,
sup EZ {Zn (x; Z)} = sup r+2
h
x∈X
x∈X h
XPr
(%)

where mPr (w) = E(|Y |% | XPr = w) ∀ % ∈ {0, 1} and w ∈ XPr . Using the
integrability of φ(·), we then have:
 0

(%) Z
CP r
Pr x − w
φ
dw
r+2
h
x∈X h
XPr

Z
(%) Z

CPr
φ(−ψ) dψ ≤ 2
φ(ψ)dψ = O h−2 ,
h
An
Rr
x

e (%) (x; Z)} ≤ sup
sup EZ {Z
n

x∈X

(%)

C
≤ sup P2r
x∈X h

(%)

(%)

where CPr = supw∈XPr mPr (w)fPr (w) < ∞ due to Assumption 4.1 (iii),
and Anx = {ψ : (P0r x + hψ) ∈ XPr }, as before. It then follows, similar to the
proof of (27), from a simple application of Hoeffding’s inequality that


e (2) − V
e (1) = Op n− 21 h−2 .
(35)
U
n,%
n,%
−2 + n−1 h−(r+2) ). Hence,
e %∗
Using (33)-(35), we finally have: Z
n = Op (h
Z




%∗
b
b (%) (x) (Pn + PX )(dx) ≤ Op αn Z
e %∗ = Op αn + αn
(36) Zn = Z
,
n
n
h2
nhr+2
!
1
1
n
o
n 2 αn2
n 2 αn2
∗ b(2)
(37) and
Gn ζ n,%,λ (·) ≤ Op
+
∀ % ∈ {0, 1},
h2
nhr+2

xviii

A. CHAKRABORTTY AND T. CAI

where the final bound in (37) follows from (32). The desired result in (A.3)
now follows by applying (30), (31) and (37) to (29) using the linearity of
G∗n (·). The proof of the lemma is now complete. (Note that conditions (i),
(iv) and (viii) in Assumption 4.1 were actually not used in this proof).
VI.3. Proof of Theorem 4.2. Finally, to establish the result of Theorem
−
4.2, let λ0 (x) = →
x which is measurable and bounded on X . Further, with
G∗n (·) as defined in Appendix A.1, note that Gn,K for K = 1 is given by:
b r ) − m(·
Gn,K = G∗n {λ0 (·)m
e ? (· ; Pr )} + G∗n [λ0 (·){m(·
b ;P
e ; Pr )}],

(38)

due to linearity of G∗n (·). Now, using Lemma A.2, we have:
1

G∗n {λ0 (·)m
e ? (· ; Pr )} = Op (n 2 a2n,2 ) = Op (a∗n,2 ).

(39)

b r ) − m(·
The second term G∗n [λ0 (·){m(·
b ;P
e ; Pr )}] in (38) can be written as:
(1)

(2)

(3)

(4)

b
b
b
b
G∗n [λ0 (·){T
n,Pr (·) − Tn,Pr (·) − Tn,Pr (·) + Tn,Pr (·)}]



1
1
2
∗
2
2
= Op b(2)
n + n an,1 + n an,1 an,2 = Op an,1 ,

(40)

where with slight abuse of notation,
b e
a−e
a
b (2) (x) = a(b − b) ,
b (1) (x) = b
,
T
T
n,Pr
n,Pr
b
b2
a−e
a)(eb − b) (b
a−e
a)(bb − eb)
b (3) (x) = (b
T
+
, and
n,Pr
eb bb
b eb
a(bb − eb)2
(e
a − a)(bb − eb) a(bb − eb)(eb − b)(b + eb)
b (4) (x) = e
T
+
,
−
n,Pr
bb b2
b2
(b eb)(b bb)
with (a, b) = {l(x; Pr ), f (x; Pr )}, (e
a, eb) = {e
l(x; Pr ), fe(x; Pr )} and (b
a, bb) =
b r ), fb(x; P
b r )}}.
{b
l(x; P
For (40), the starting expansion is due to a linearization similar to (28),
(1)
while the final rate is due to the following: note that λ0,Pr (·) ≡ b−1 λ0 (·) and
(2)

λ0,Pr (·) ≡ ab−2 λ0 (·) are both bounded a.s. [PX ] due to Assumption 4.1 (iii)(iv) and the boundedness of λ0 (·). Hence, using these as choices of ‘λ(·)’ in
(1)
b (1) (·)}] = Op (b(2)
Lemma A.3, we have: G∗n {(b
a−e
a)λ0,Pr (·)} = G∗n [λ0 (·){T
n )
n,Pr
(2)
(2)
(2)
∗
∗
b
e
b
and G {(b − b)λ
(·)} = G [λ0 (·){T
(·)}] = Op (bn ) respectively. Furn

0,Pr

n

n,Pr

b (s) (x)k ≤ Op (a2 + an,1 an,2 )
ther, note that for each s ∈ {3, 4}, supx∈X kT
n,1
n,Pr
which follows from repeated use of (15), (24) along with Assumption 4.1

SEMI-SUPERVISED LINEAR REGRESSION

xix

(iii)-(iv). Consequently, with λ0 (x) bounded a.s. [PX ], for each s ∈ {3, 4},
b (s) (·)}] is bounded by: Op (n 21 a2 + n 12 an,1 an,2 ). Combining all
G∗n [λ0 (·){T
n,1
n,Pr
these results using the linearity of G∗n (·) and noting that with a∗n,2 = o(1),
(2)

1

1

(bn + n 2 a2n,1 + n 2 an,1 an,2 ) = O(a∗n,1 ), (40) now follows and, along with (39)
and (38), implies: Gn,K = Op (a∗n,1 + a∗n,2 ) as claimed in Theorem 4.2. Lastly,
using this in (3.10), the expansion in (4.3) and its associated implications
follow. The proof of Theorem 4.2 is now complete.
References.
Andrews, D. W. K. (1995). Nonparametric Kernel Estimation for Semiparametric Models. Econometric Theory 11 560-586.
Belkin, M., Niyogi, P. and Sindhwani, V. (2006). Manifold Regularization : A Geometric Framework for Learning from Labeled and Unlabeled Examples. The Journal of
Machine Learning Research 7 2399-2434.
Castelli, V. and Cover, T. M. (1995). The Exponential Value of Labeled Samples.
Pattern Recognition Letters 16 105-111.
Castelli, V. and Cover, T. M. (1996). The Relative Value of Labeled and Unlabeled
Samples in Pattern Recognition with an Unknown Mixing Parameter. IEEE Transactions on Information Theory 42 2102-2117.
Chakrabortty, A. and Cai, T. (2017). Supplement to “Efficient and Adaptive Linear
Regression in Semi-Supervised Settings”. DOI: Completed by the typesetter.
Chapelle, O., Schölkopf, B. and Zien, A. (2006). Semi-Supervised Learning. MIT
Press, Cambridge, MA, USA.
Cook, R. D. (1998). Principal Hessian Directions Revisited (with Discussion). Journal of
the American Statistical Association 93 84-100.
Cook, R. D. and Lee, H. (1999). Dimension Reduction in Binary Response Regression.
Journal of the American Statistical Association 94 1187-1200.
Cook, R. D. and Weisberg, S. (1991). Discussion of “Sliced Inverse Regression” by
K.-C. Li. Journal of the American Statistical Association 86 328-332.
Cozman, F. G. and Cohen, I. (2001). Unlabeled Data Can Degrade Classification Performance of Generative Classifiers. Technical Report No. HPL-2001-234, HP Laboratories,
Palo Alto, CA, USA.
Cozman, F. G., Cohen, I. and Cirelo, M. C. (2003). Semi-Supervised Learning of
Mixture Models. In Proceedings of the Twentieth ICML 99-106.
Duan, N. and Li, K.-C. (1991). Sliced Regression: A Link-Free Regression Method. The
Annals of Statistics 19 505-530.
Hansen, B. E. (2008). Uniform Convergence Rates for Kernel Estimation with Dependent
Data. Econometric Theory 24 726-748.
Kawakita, M. and Kanamori, T. (2013). Semi-Supervised Learning with Density-Ratio
Estimation. Machine Learning 91 189-209.
Kohane, I. S. (2011). Using Electronic Health Records to Drive Discovery in Disease
Genomics. Nature Reviews Genetics 12 417-428.
Lafferty, J. D. and Wasserman, L. (2007). Statistical Analysis of Semi-Supervised
Regression. Advances in Neural Information Processing Systems 20 801-808.
Li, K.-C. (1991). Sliced Inverse Regression for Dimension Reduction. Journal of the American Statistical Association 86 316-327.

xx

A. CHAKRABORTTY AND T. CAI

Li, K.-C. (1992). On Principal Hessian Directions for Data Visualization and Dimension
Reduction: Another Application of Stein’s Lemma. Journal of the American Statistical
Association 87 1025-1039.
Liao, K. P., Cai, T., Gainer, V. et al. (2010). Electronic Medical Records for Discovery
Research in Rheumatoid Arthritis. Arthritis Care and Research 62 1120-1127.
Masry, E. (1996). Multivariate Local Polynomial Regression for Time Series: Uniform
Strong Consistency and Rates. Journal of Time Series Analysis 17 571-600.
Newey, W. K. (1994). Kernel Estimator of Partial Means and a Generalized Variance
Estimator. Econometric Theory 10 1-21.
Newey, W. K., Hsieh, F. and Robins, J. (1998). Undersmoothing and Bias Corrected
Functional Estimation. Technical Report No. 98-17, Dept. of Economics, MIT, USA.
Newey, W. K. and McFadden, D. (1994). Large Sample Estimation and Hypothesis
Testing. Handbook of Econometrics 4 2111–2245.
Nigam, K. P. (2001). Using Unlabeled Data to Improve Text Classification. PhD thesis,
Carnegie Mellon University, USA. CMU-CS-01-126.
Nigam, K., McCallum, A. K., Thrun, S. and Mitchell, T. (2000). Text Classification
from Labeled and Unlabeled Documents Using EM. Machine Learning 39 103-134.
Seeger, M. (2002). Learning with Labeled and Unlabeled Data. Technical Report No.
EPFL-REPORT-161327, University of Edinburgh, UK.
Sokolovska, N., Cappé, O. and Yvon, F. (2008). The Asymptotics of Semi-Supervised
Learning in Discriminative Probabilistic Models. In Proceedings of the Twenty Fifth
ICML 984-991.
Van der Vaart, A. W. (2000). Asymptotic Statistics 3. Cambridge University Press.
Zhang, T. and Oles, F. J. (2000). The Value of Unlabeled Data for Classification Problems. In Proceedings of the Seventeenth ICML 1191-1198.
Zhu, X. (2005). Semi-Supervised Learning through Graphs. PhD thesis, Carnegie Mellon
University, USA. CMU-LTI-05-192.
Zhu, X. (2008). Semi-Supervised Learning Literature Survey. Technical Report No. 1530,
Computer Sciences, University of Wisconsin-Madison, USA.
Zhu, L.-X. and Ng, K. W. (1995). Asymptotics of Sliced Inverse Regression. Statistica
Sinica 5 727-736.
Abhishek Chakrabortty
Department of Statistics
University of Pennsylvania
3730 Walnut Street
Jon M. Huntsman Hall, 4th Floor
Philadelphia, PA 19104, USA.
E-mail: abhich@mail.med.upenn.edu

Tianxi Cai
Department of Biostatistics
Harvard University
655 Huntington Avenue
Building 2, 4th Floor
Boston, MA 02115, USA.
E-mail: tcai@hsph.harvard.edu

