# Semi-Supervised linear regression

David Azriel1, Lawrence D. Brown2, Michael Sklar3, Richard Berk2, Andreas
Buja2 and Linda Zhao2 Corresponding author: David Azriel,
<EMAIL_ADDRESS>The work of Lawrence D. Brown and Linda Zhao was
supported in part by NSF grant DMS-1512084

###### Abstract

We study a regression problem where for some part of the data we observe both
the label variable ($Y$) and the predictors (${\bf X}$), while for other part
of the data only the predictors are given. Such a problem arises, for example,
when observations of the label variable are costly and may require a skilled
human agent. If the conditional expectation $E[Y|{\bf X}]$ is exactly linear
in ${\bf X}$ then typically the additional observations of the ${\bf X}$’s do
not contain useful information, but otherwise the unlabeled data can be
informative. In this case, our aim is at constructing the best linear
predictor. We suggest improved alternative estimates to the naive standard
procedures that depend only on the labeled data. Our estimation method can be
easily implemented and has simply described asymptotic properties. The new
estimates asymptotically dominate the usual standard procedures under certain
non-linearity condition of $E[Y|{\bf X}]$; otherwise, they are asymptotically
equivalent. The performance of the new estimator for small sample size is
investigated in an extensive simulation study. A real data example of
inferring homeless population is used to illustrate the new methodology.

Technion – Israel institute of Technology1, Wharton – University of
Pennsylvania2 and Stanford University3

##  1 Introduction

The term “semi-supervised learning” was coined in the machine learning
literature to describe a situation in which some of the data is labeled while
the rest of the data is unlabeled (Merz et al., 1992). Such situations occur
when the label variable is difficult to observe and may require a complicated
or expensive procedure. A typical example is web document classification,
where the classification is done by a human agent while there are many more
unlabeled on-line documents. Specifically, a sample of $n$ observations from
the joint distribution of $({\bf X},Y)$ is given, where $Y$ is a one-
dimensional label variable and ${\bf X}$ is a p-dimensional vector of
covariates or predictors. Also, an additional sample of size $m$ is observed
where only ${\bf X}$ is given. The purpose is to study procedures that make
use of the additional unlabeled data to better capture the shape of the
underlying joint distribution of the labeled data.

A large body of literature focuses on the case that $Y$ takes a small number
of values and the problem reduces to a classification task; see e.g., Zhu
(2005) and references therein. Wang et al. (2009) divide the different methods
into two approaches: distributional and margin-based. The distributional
approach relies on an assumption relating the conditional expectation
$E[Y|{\bf X}]$ to the marginal distribution of ${\bf X}$ and the margin-based
approach uses the extra information on $X$ for estimating the Bayes decision
boundary; see Liang et al. (2008) for a Bayesian perspective within this basic
approach. Other works consider continuous $Y$’s and use the unlabeled data to
learn the structure of the ${\bf X}$’s in order to better estimate a non-
parametric regression (Zhou, 2005; Lafferty and Wasserman, 2008; Johnson and
Zhang, 2008). These works could be very helpful in situations where non-
parametric regression is useful and unlabeled data are available.

Here we follow a third methodology. This methodology makes use of possible
unmodeled non-linearities in $E(Y|{\bf X})$ and adapts information from the
unlabeled data to partially adjust for these. In a recent manuscript,
Chakrabortty and Cai (2016) follow the same methodological plan but the
technical details within their approach are then quite different from ours.
Their estimates offer some advantages and some disadvantages in relation to
ours, which we intend to explore in a future study.

We aim at estimating the vector ${\boldsymbol{\beta}}$ composed of the
parameters of the best linear predictor of $Y$, but we do not necessarily
assume that $E[Y|{\bf X}]$ is exactly linear in ${\bf X}$. The estimates of
${\boldsymbol{\beta}}$ can then be used to provide predictions. The standard
approach in the statistical literature may be best summarized by the following
quote from Little (1992):

> The related problem of missing values in the outcome $Y$ was prominent in
> the early history of missing-data methods, but is less interesting in the
> following sense: If the $X$’s are complete and the missing values of $Y$ are
> missing at random, then the incomplete cases contribute no information to
> the regression of $Y$ on $X_{1},\ldots,X_{p}$.

But see Cochran (1977, Chapter 7) for a different view more closely in tune
with our development.

Buja et al. (2016) show that ${\boldsymbol{\beta}}$ does not depend on the
distribution of ${\bf X}$ if and only if $E[Y|{\bf X}]$ is linear in ${\bf
X}$. When the conditional expectation is linear, there is typically at most a
limited amount of additional information in the unlabeled data. Brown (1990)
shows that even in this case there may be some useful information, but it will
not provide an asymptotic advantage in the manner we suggest in the current
treatment when the conditional expectation is not linear. The assumption that
$E[Y|{\bf X}]$ is linear is unrealistic in many situations, and we show that
in the absence of such an assumption the unlabeled sample can be used to
provide useful information for the estimation of ${\boldsymbol{\beta}}$.

We consider two scenarios. In the first one the distribution of ${\bf X}$ is
known exactly. This is equivalent to having infinitely many unlabeled
observations ($m=\infty$). In the second, more frequently encountered,
scenario $m<\infty$ but it is still significant in the sense that $m/n$ is
bounded away from 0. We call these two scenarios total and partial
information, respectively. In both situations we provide an asymptotically
better estimate of ${\boldsymbol{\beta}}$ than the standard least squares
estimator, and hence also a better linear predictor. Our new procedure is
closely related to the semi-supervised estimations of means, which is the
topic of a recent manuscript (Zhang et al., 2016).

To understand our basic methodology, consider the simple linear regression
model,

$Y=\alpha+\beta X+\delta,$ (1)

where $\alpha$, $\beta$ are least squares coefficients and $\delta$ is a
remainder term; these will be defined in more detail shortly. When $E(X)$ and
$E(X^{2})$ are known, we construct an estimator of $\beta$ that asymptotically
dominates the least squares estimator. This is done by replacing model (1)
with a different model where $\beta$ is the intercept and also the expectation
of the newly defined label variable. For this latter model, Zhang et al.
(2016) show that the intercept estimator dominates the simple empirical
estimate in a semi-supervised setting. We use their ideas to show that the
intercept estimator in this new model is better than the least squares
estimator from the original model as explained below.

To present the intercept model and the methodology, begin by stating the
general form of the best linear approximation. The approximation is used in
two different models, and therefore we introduce now the notation ${\bf U}$
and $W$ instead of ${\bf X}$ and $Y$. Suppose that $W\in\mathbb{R}$ and ${\bf
U}\in\mathbb{R}^{d}$ are random variables with joint distribution $G$ and
finite second moments. The best linear predictor is

${\boldsymbol{\theta}}=({\theta}_{1},\ldots,{\theta}_{d})^{T}={\arg\min}_{\tilde{\boldsymbol{\theta}}\in{\mathbb{R}}^{d}}E(W-\sum_{j=1}^{d}U_{j}\tilde{\theta}_{j})^{2}=\left\\{E({\bf
U}{\bf U}^{T})\right\\}^{-1}E({\bf U}W).$ (2)

Notice that this is a population version of the least squares where
$\sum_{j=1}^{d}U_{j}\tilde{\theta}_{j}$ minimizes the $L_{2}$ distance in the
population from $W$ to any linear function of ${\bf U}$. It follows that

$W=\sum_{j=1}^{d}\theta_{j}U_{j}+r,$ (3)

where the remainder term $r=W-\sum_{j=1}^{d}\theta_{j}U_{j}$ is orthogonal to
${\bf U}$, i.e., $E(r{\bf U})=0$. Given a sample of $n$ observations from $G$,
the standard least squares estimate, $\hat{\boldsymbol{\theta}}_{LSE}$,
satisfies asymptotically (White, 1980)

$\sqrt{n}(\hat{\boldsymbol{\theta}}_{LSE}-{\boldsymbol{\theta}})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N\left(0,\left\\{E({\bf U}{\bf
U}^{T})\right\\}^{-1}E\left(r^{2}{\bf U}{\bf U}^{T}\right)\left\\{E({\bf
U}{\bf U}^{T})\right\\}^{-1}\right).$ (4)

Unlike the standard fixed ${\bf X}$ assumption, the asymptotic variance has a
“sandwich” form, $\left\\{E({\bf U}{\bf U}^{T})\right\\}^{-1}$ forming the
“bread” and $E\left(r^{2}{\bf U}{\bf U}^{T}\right)$ the “meat”. See Buja et
al. (2016) for further discussion of this form of the sandwich.

We now return to model (1) and assume that $E(X)$ and $E(X^{2})$ are known. We
can therefore assume w.l.o.g that $X$ is standardized (i.e., $E(X)=0$ and
$E(X^{2})=1$). In this case we can write model (1) with $\alpha=E(Y)$,
$\beta=E(XY)$ and $\delta=Y-(\alpha+\beta X)$; here we consider model (3) with
$W=Y$, $U_{1}=1$, $U_{2}=X$ and the remainder term is $r=\delta$. The standard
LSE is a consistent estimate for $\beta$ and satisfies, according to (4),

$\sqrt{n}(\hat{\beta}_{LSE}-\beta)\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N\left\\{0,E(\delta X)^{2}\right\\}.$

Our aim it to better estimate $\beta=E(XY)$. To this end, we multiply (1) by
$X$ setting $XY$ to be the labeled variable $W$. Furthermore, we also center
the varibale ${\bf U}$ setting $E(XY)=\beta$ to be the intercept.
Specifically,

$XY=\alpha X+\beta X^{2}+X\delta=\beta+aX+b(X^{2}-1)+\tilde{\delta}.$ (5)

Now we consider model (3) with $W=XY$, $U_{1}=1$, $U_{2}=X$, $U_{3}=X^{2}-1$
and $r=\tilde{\delta}$; here $a,b$ are $\theta_{2},\theta_{3}$ defined by (2).
In setting $\beta=E(XY)$ to be the intercept coefficient we used that
$EX=E(X^{2})-1=0$, i.e., that we know the first and second moment of $X$. We
define our new estimator, $\hat{\beta}_{TI}$ (TI for total information), to be
the intercept estimator of (5). This estimate has a simple closed formula,
which is given below in (14). The sandwich theorem (4) implies that

$\sqrt{n}(\hat{\beta}_{TI}-\beta)\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,E\tilde{\delta}^{2}).$

Since $\delta X-\tilde{\delta}=X(a-\alpha)+(X^{2}-1)(b-\beta)$ is a linear
function of $X,X^{2}$, then it is orthogonal to the remainder term of (5),
which is $\tilde{\delta}$, (i.e., $E(\delta
X-\tilde{\delta})\tilde{\delta}=0$) and therefore,

$E(\delta X)^{2}=E\\{\tilde{\delta}+(\delta
X-\tilde{\delta})\\}^{2}=E\tilde{\delta}^{2}+E(\delta X-\tilde{\delta})^{2}.$
(6)

This implies that the asymptotic variance of $\hat{\beta}_{TI}$ is smaller
than that of $\hat{\beta}_{LSE}$ with equality iff $\delta
X\equiv\tilde{\delta}$. The latter occurs iff $\alpha=a$ and $b=\beta$ or
equivalently that $\delta$ is uncorrelated with both $X^{2}$ and $X^{3}$. This
occurs when the non-linear part of $E[Y|X]$ (if exists) is uncorrelated with
$X^{2}$ and $X^{3}$. In this case, Models (1) and (5) are essentially the same
model and nothing is gained in the new methodology. On the other hand, when
$\delta$ is correlated with either $X^{2}$ or $X^{3}$, then Models (1) and (5)
are different and $\hat{\beta}_{TI}$ has smaller asymptotic variance than
$\hat{\beta}_{LSE}$. We further show that a similar decomposition to (6) holds
for the partial information case and we generalize the results to the
$p$-dimensional case.

The rest of the paper is organized as follows. Section 2 provides the basic
setting the and loss functions that we use. The main results are given in
Section 3. Section 4 describes an extensive simulation study and Section 5
discusses estimation of the asymptotic covariance matrix of the estimates. An
implementation of the new methodology to infer homeless population in Los
Angeles is discussed in Section 6. Section 7 concludes with final remarks. The
proofs are given in Section 8.

##  2 Preliminaries

Consider a sample of $n$ i.i.d observations $({\bf
X}^{(1)},Y^{(1)}),\ldots,({\bf X}^{(n)},Y^{(n)})$ from a joint distribution
$G$, where ${\bf X}\in\mathbb{R}^{p}$, $Y\in\mathbb{R}$, and an additional set
of $m$ independent observations $({\bf X}^{(n+1)}\ldots,{\bf X}^{(n+m)})$ from
the marginal distribution of ${\bf X}$. We use super-index to denote the
number of the observation, and sub-index to denote coordinates of ${\bf X}$.
The notation ${\bf X},Y$ without super-index denotes a random vector whose
distribution is $G$.

We write $\vec{\bf X}=(1,X_{1},\ldots,X_{p})^{T}$ to be a vector ${\bf X}$
with an additional constant 1 to accommodate an intercept term. Assume that
the second moments of $G$ exists and that the matrix $E\left(\vec{\bf
X}\vec{\bf X}^{T}\right)$ is invertible. Then, we can define

$(\alpha,{\boldsymbol{\beta}})={\arg\min}_{\tilde{\alpha}\in\mathbb{R},\tilde{\boldsymbol{\beta}}\in\mathbb{R}^{p}}E(Y-\tilde{\alpha}-\tilde{\boldsymbol{\beta}}^{T}X)^{2}=\left\\{E\left(\vec{\bf
X}\vec{\bf X}^{T}\right)\right\\}^{-1}E\left(\vec{\bf X}Y\right).$ (7)

In the presence of non-linearity, ${\boldsymbol{\beta}}$ is still a meaningful
parameter that describes the overall association between $Y$ and ${\bf X}$
(Buja et al., 2016). We have in mind two related purposes. The first purpose
is just to better estimate the parameters of interest, while the second
purpose pertains to prediction. The latter is formalized now. Suppose that an
independent observation $({\bf X}^{*},Y^{*})\sim G$ is given. The optimal
linear predictor is $\alpha+{\boldsymbol{\beta}}^{T}X^{*}$. We consider the
excess loss of an estimator $\tilde{\alpha},\tilde{\boldsymbol{\beta}}$

$L(\tilde{\alpha},\tilde{\boldsymbol{\beta}})=\left(Y^{*}-\tilde{\alpha}-\tilde{\boldsymbol{\beta}}^{T}{\bf
X}^{*}\right)^{2}-\left(Y^{*}-{\alpha}-{\boldsymbol{\beta}}^{T}{\bf
X}^{*}\right)^{2}.$

We have that (see Lemma 1)

$EL(\tilde{\alpha},\tilde{\boldsymbol{\beta}})=E\Big{[}\bar{Y}-\sum_{j=1}^{p}\tilde{\beta}_{j}\\{{\bar{X}}_{j}-E(X_{j})\\}-E(Y)\Big{]}^{2}+E\left\\{({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
M}({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})\right\\}=E\tilde{L}(\tilde{\alpha},\tilde{\boldsymbol{\beta}}),$

where $\tilde{L}$ is the expression inside the expectation and ${\bf M}$ is
the covariance matrix of ${\bf X}$. We further define the asymptotic risk as

$R(\tilde{\alpha},\tilde{\boldsymbol{\beta}})=\lim_{B\to\infty}\lim_{n\to\infty}E\min\\{n\tilde{L}(\tilde{\alpha},\tilde{\boldsymbol{\beta}}),B\\}.$

The loss is of order $1/n$ and therefore we consider expectation of
$n\tilde{L}$. This is truncated by an arbitrarily large number $B$, since when
the loss is large, it makes sense not to penalize any further. Also, the
truncation helps to avoid issues of uniform integrability. This is done for
example in Le Cam and Yang (1990), Chapter 5. The following proposition
provides a simple expression for the asymptotic risk.

######  Proposition 1.

Let $\tilde{\alpha},\tilde{\boldsymbol{\beta}}$ satisfies

$\tilde{\alpha}=\bar{Y}-\sum_{j=1}^{p}\tilde{\beta}_{j}{\bar{X}}_{j},$ (8)

and assume that $\tilde{\boldsymbol{\beta}}$ satisfies
$\sqrt{n}(\tilde{\boldsymbol{\beta}}-{\boldsymbol{\beta}})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,{\boldsymbol{\Sigma}})$ then,

${R}(\tilde{\alpha},\tilde{\boldsymbol{\beta}})=E\Big{[}Y-\sum_{j=1}^{p}{\beta}_{j}\\{{X}_{j}-E(X_{j})\\}-E(Y)\Big{]}^{2}+{\rm
Trace}({\bf M}{\boldsymbol{\Sigma}}),$ (9)

where ${\bf M}$ is the covariance matrix of ${\bf X}$.

The first term in (9) does not depend on the distribution of
$\tilde{\boldsymbol{\beta}}$. Hence, Proposition 1 shows that the excess risk
is minimized when $\text{Trace}({\bf M}{\boldsymbol{\Sigma}})$ is small. Thus,
we aim at estimators $\tilde{\boldsymbol{\beta}}$ with asymptotic distribution
$N(0,{\boldsymbol{\Sigma}})$ such that ${\boldsymbol{\Sigma}}$ is “smaller”
than the covariance matrix of LSE, in the sense that the difference is
positive semidefinite. Such an estimator asymptotically better estimates
${\boldsymbol{\beta}}$ and also has smaller asymptotic excess risk.

For two estimates $(\tilde{\alpha}^{(1)},\tilde{\bf\beta}^{(1)})$ and
$(\tilde{\alpha}^{(2)},\tilde{\bf\beta}^{(2)})$, Proposition 1 implies that

$EL(\tilde{\alpha}^{(1)},\tilde{\boldsymbol{\beta}}^{(1)})-EL(\tilde{\alpha}^{(2)},\tilde{\boldsymbol{\beta}}^{(2)})\approx\text{Trace}\left\\{{\bf
M}({\boldsymbol{\Sigma}}^{(1)}-{\boldsymbol{\Sigma}}^{(2)})\right\\}/n,$

and therefore the difference of the prediction errors is

$E\left\\{Y^{*}-\tilde{\alpha}^{(1)}-\left(\tilde{\boldsymbol{\beta}}^{(1)}\right)^{T}{\bf
X}^{*}\right\\}^{2}-E\left\\{Y^{*}-\tilde{\alpha}^{(2)}-\left(\tilde{\boldsymbol{\beta}}^{(2)}\right)^{T}{\bf
X}^{*}\right\\}^{2}\approx\text{Trace}\left\\{{\bf
M}({\boldsymbol{\Sigma}}^{(1)}-{\boldsymbol{\Sigma}}^{(2)})\right\\}/n.$ (10)

It follows that the difference of the errors of the prediction of the mean is
also

$E\left\\{E(Y^{*}|{\bf
X}^{*})-\tilde{\alpha}^{(1)}-\left(\tilde{\boldsymbol{\beta}}^{(1)}\right)^{T}{\bf
X}^{*}\right\\}^{2}-E\left\\{E(Y^{*}|{\bf
X}^{*})-\tilde{\alpha}^{(2)}-\left(\tilde{\boldsymbol{\beta}}^{(2)}\right)^{T}{\bf
X}^{*}\right\\}^{2}\approx\text{Trace}\left\\{{\bf
M}({\boldsymbol{\Sigma}}^{(1)}-{\boldsymbol{\Sigma}}^{(2)})\right\\}/n.$ (11)

Notice that we estimate the intercept, as in (8), using ${\bar{X}}_{j}$ rather
than the possibly known expectations $E(X_{j})$. This is in accordance with
the findings of Zhang et al. (2016).

##  3 Main results

In this section we provide the main theoretical results of the paper. Before
we deal with the p-dimensional case, we first introduce the results for one-
dimensional ${\bf X}$. The reason is twofold: first, the presentation in the
one-dimensional case is simpler and captures the main ideas, and second, our
results for the p-dimensional ${\bf X}$ are obtained by reducing the
regression problem to p problems, each of which is closely related to the one
dimensional ${\bf X}$ regression.

###  3.1 One dimensional ${X}$

In this section we study the one dimensional case. Summary of the notation
used here is presented in Table 1. When ${\bf G}$ has finite second moments we
can write

$Y=\alpha+\beta X+\delta,$ (12)

where $\alpha=E(Y)-\beta E(X)$, $\beta=\frac{{\rm Cov}(X,Y)}{Var(X)}$ and
$\delta=Y-\alpha-\beta X$. Equation (12) is the best linear approximation in
the population in the sense that $\alpha+\beta X$ minimizes
$E(\tilde{\alpha}+\tilde{\beta}X-Y)^{2}$ over all
$\tilde{\alpha},\tilde{\beta}\in\mathbb{R}$ and $\delta$ is orthogonal to
$1,X$, i.e., $E(\delta)=E(\delta X)=0$. The regular LSE is

$\hat{\beta}_{LSE}=\frac{\sum_{i=1}^{n}(X^{(i)}-\bar{X})Y^{(i)}}{\sum_{i=1}^{n}(X^{(i)}-\bar{X})^{2}}.$

This is a consistent and asymptotically unbiased and asymptotically normal
estimator for $\beta$. See Buja et al. (2016) for contemporary discussion.

Table 1: Summary of the notation used in Section 3.1.  Basic model | $Y=\alpha+\beta X+\delta$ |  $\alpha=E(Y)-\beta E(X)$, $\beta=\frac{{\rm Cov}(X,Y)}{Var(X)}$ and $\delta=Y-\alpha-\beta X$
---|---|---
Intercept model | $W=\beta+aU_{1}+bU_{2}+\tilde{\delta}$ |  $W=\frac{Y\\{X-E(X)\\}}{Var(X)}$, $U_{1}=\frac{X-E(X)}{Var(X)}$, $U_{2}=\frac{\\{X-E(X)\\}X}{Var(X)}-1$
|  |  $a,b$ are $\theta_{2},\theta_{3}$ defined by (2) for model (13)
|  | and $\tilde{\delta}$ is the remainder term
TI estimator | $\hat{\beta}_{TI}=\bar{W}-\hat{a}\bar{U}_{1}-\hat{b}\bar{U}_{2}$ |  $\bar{W}=\frac{1}{n}\sum_{i=1}^{n}W^{(i)}$, $\bar{U}_{1}=\frac{1}{n}\sum_{i=1}^{n}U_{1}^{(i)}$, $\bar{U}_{2}=\frac{1}{n}\sum_{i=1}^{n}U_{2}^{(i)}$
|  |  $\hat{a},\hat{b}$ are the LSE of model (13)
PI estimator | $\hat{\beta}_{PI}=\bar{\check{W}}-\hat{a}\bar{\check{U}}_{1}-\hat{b}\bar{\check{U}}_{2}$ |  $\bar{\check{W}},\bar{\check{U}}_{1},\bar{\check{U}}_{2}$ are the means over the labeled sample,
|  | $\check{W}^{(i)}=\frac{\\{X^{(i)}-\check{E}(X)\\}Y^{(i)}}{\widecheck{Var}(X)},\check{U}_{1}^{(i)}=\frac{\\{X^{(i)}-\check{E}(X)\\}}{\widecheck{Var}(X)}$
|  |  $\check{U}_{2}^{(i)}=\frac{\\{X^{(i)}-\check{E}(X)\\}X^{(i)}}{\widecheck{Var}(X)}-1$, $\hat{a},\hat{b}$ are the LSE
|  | of model (13) with $\check{W},\check{U}_{1},\check{U}_{2}$ instead of $W,U_{1},U_{2}$
Asymptotics | $\nu=\lim\frac{n}{n+m}$ |  $n$ ($m$) the labeled (unlabeled) sample size
|  $\sigma^{2}_{LSE}$, $\sigma^{2}_{PI}$, $\sigma^{2}_{TI}$ | asymptotic variance of the LSE, PI, TI estimators
|  | as in Theorem 1
Main results | $\sigma^{2}_{LSE}=\sigma^{2}_{TI}+\sigma^{2}_{\rm diff}$ | $\sigma^{2}_{\rm diff}=E\left[\frac{\delta\\{X-E(X)\\}}{Var(X)}-\tilde{\delta}\right]^{2}$
| $\sigma^{2}_{PI}=\sigma^{2}_{TI}+\nu\sigma^{2}_{\rm diff}$ |

For the total information estimator (TI) we consider the following regression
model, which is obtained by multiplying (12) by $\frac{X-E(X)}{Var(X)}$,

$W=\beta+aU_{1}+bU_{2}+\tilde{\delta},$ (13)

where $W=\frac{Y\\{X-E(X)\\}}{Var(X)}$, $U_{1}=\frac{X-E(X)}{Var(X)}$,
$U_{2}=\frac{\\{X-E(X)\\}X}{Var(X)}-1$, $a,b$ are the coefficients
$\theta_{2},\theta_{3}$ of the best linear approximation defined by (2) and
$\tilde{\delta}$ is the remainder term. Here again $\beta+aU_{1}+bU_{2}$ is
the best linear approximation of $W$ in the population. The multiplication
term makes the expectation of $W$ to be
$E\left[\frac{Y\\{X-E(X)\\}}{Var(X)}\right]=\frac{Cov(X,Y)}{Var(X)}=\beta$.
The covariates $U_{1},U_{2}$ are the predictors $1,X$ in (12) multiplied by
$\frac{X-E(X)}{Var(X)}$. We subtract 1 in $U_{2}$ in order to set
$E(U_{2})=0$. Since $E(U_{1})=E(U_{2})=0$, then the intercept is $E(W)=\beta$.
Thus, we define the total information estimator to be the intercept estimator
of (13), i.e.,

$\hat{\beta}_{TI}=\bar{W}-\hat{a}\bar{U}_{1}-\hat{b}\bar{U}_{2},$ (14)

where $\hat{a},\hat{b}$ are the regular LSE of model (13), and $\bar{\cdot}$
denotes the mean over the supervised sample with $n$ observations.

For the partial information (PI) estimator we estimate $W,U_{1},U_{2}$ as
follows:

$\check{W}^{(i)}=\frac{\\{X^{(i)}-\check{E}(X)\\}Y^{(i)}}{\widecheck{Var}(X)},\check{U}_{1}^{(i)}=\frac{\\{X^{(i)}-\check{E}(X)\\}}{\widecheck{Var}(X)}\text{
and
}\check{U}_{2}^{(i)}=\frac{\\{X^{(i)}-\check{E}(X)\\}X^{(i)}}{\widecheck{Var}(X)}-1~{},i=1,\ldots,n,$

where $\check{E}(X)=\frac{1}{n+m}\sum_{i=1}^{n+m}X^{(i)}$,
$\widecheck{Var}(X)=\check{E}(X^{2})-\\{\check{E}(X)\\}^{2}$, and
$\check{E}(X^{2})=\frac{1}{n+m}\sum_{i=1}^{n+m}\\{X^{(i)}\\}^{2}$. The partial
information estimator is

$\hat{\beta}_{PI}=\bar{\check{W}}-\hat{a}\bar{\check{U}}_{1}-\hat{b}\bar{\check{U}}_{2},$

where $\hat{a},\hat{b}$ are the regular LSE of the regression model (13) with
$\check{W},\check{U}_{1},\check{U}_{2}$ instead of $W,U_{1},U_{2}$. (The
estimates $\hat{a},\hat{b}$ are different in the total and partial information
cases but the same notation is presented for simplicity.) We use $\bar{\cdot}$
(respectively, $\check{\cdot}$) to denote empirical mean with respect to the
labeled $n$ (respectively, full $n+m$) sample. The following theorem states
the asymptotic distribution of $\hat{\beta}_{LSE},\hat{\beta}_{TI}$ and
$\hat{\beta}_{PI}$. The first part of the theorem is known and is stated here
for comparison purposes.

######  Theorem 1.

  1. (i)

Suppose that $Var(X)\in(0,\infty)$ and that $\delta\\{X-E(X)\\}$ has finite
second moment, then,

$\sqrt{n}(\hat{\beta}_{LSE}-\beta)\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,\sigma^{2}_{LSE}),$

where
$\sigma^{2}_{LSE}=E\left[\frac{\delta^{2}\\{X-E(X)\\}^{2}}{\\{Var(X)\\}^{2}}\right]$.

  2. (ii)

Suppose further that the vector $(W,U_{1},U_{2})$ has finite second moments
and that the matrix $E({\bf U}{\bf U}^{T})$ is invertible for ${\bf
U}=\left(1,U_{1},U_{2}\right)^{T}$. Then,

$\sqrt{n}(\hat{\beta}_{TI}-\beta)\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,\sigma^{2}_{TI}),$

where $\sigma^{2}_{TI}=E(\tilde{\delta}^{2})$ and
$\sigma^{2}_{LSE}=\sigma^{2}_{TI}+\sigma^{2}_{\rm diff}$ where
$\sigma^{2}_{\rm
diff}=E\left[\frac{\delta\\{X-E(X)\\}}{Var(X)}-\tilde{\delta}\right]^{2}$.

  3. (iii)

Suppose further that $\lim\frac{n}{n+m}=\nu$, then,

$\sqrt{n}(\hat{\beta}_{PI}-\beta)\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,\sigma^{2}_{PI}),$

where $\sigma^{2}_{PI}=\sigma^{2}_{TI}+\nu\sigma^{2}_{\rm diff}$.

Therefore, if $\sigma^{2}_{\rm diff}>0$ then
$\sigma^{2}_{TI}<\sigma^{2}_{LSE}$ and if further $\nu<1$ then
$\sigma^{2}_{PI}<\sigma^{2}_{LSE}$.

######  Corollary 1.

Theorem 1 and Proposition 1 imply that
$R(\hat{\alpha}_{LSE},\hat{\beta}_{LSE})-R(\hat{\alpha}_{PI},\hat{\beta}_{PI})=(1-\nu)\sigma^{2}_{\rm
diff}Var(X)$.

The improvement of $\hat{\beta}_{TI}$ and $\hat{\beta}_{PI}$ over
$\hat{\beta}_{LSE}$ depends on the assumption that $\sigma^{2}_{\rm diff}>0$.
The quantity $\sigma^{2}_{\rm diff}$ measures the difference between the
original regression model (12) and the intercept model (13). If $E(\delta
X^{2})=E(\delta X^{3})=0$ then $a=\alpha$ and $b=\beta$ and the two models are
essentially the same, in which case $\sigma^{2}_{\rm diff}=0$, otherwise
$\sigma^{2}_{\rm diff}>0$. In other words, when $E[Y|X]$ is non-linear and the
non-linear part is correlated with either $X^{2}$ or $X^{3}$ then
$\sigma^{2}_{\rm diff}>0$.

###  3.2 Multidimensional ${\bf X}$

We now consider the general $p$-dimensional case as described in the beginning
of Section 2. The notation of this section is summarized in Table 2. The model
can be written as

$Y=\alpha+\beta_{1}X_{1}+\cdots+\beta_{p}X_{p}+\delta,$ (15)

where $\alpha,{\boldsymbol{\beta}}$ are the coefficients of the best linear
predictor in the population defined by (7). The remainder term
$\delta=Y-\sum_{j=1}^{p}\beta_{j}X_{j}$ satisfies $E(\delta)=E(\delta
X_{1})=\cdots=E(\delta X_{p})=0$. Our aim here is at estimating
${\boldsymbol{\beta}}$. We use the adjustment representation of Buja et al.
(2016) that reduces the p-dimensional estimation procedure to $p$ separate
simple regression problems. Correspondingly, we will define our new estimates
by solving $p$ mean-estimation-problems separately, one for each coordinate
$j$.

Table 2: Summary of the notation used in Section 3.2 when $j$ is the chosen coordinate. Basic model | $Y=\alpha+\sum_{j=1}^{p}\beta_{j}X_{j}+\delta$ |  $\alpha,{\boldsymbol{\beta}}$ are defined by (7)
---|---|---
Adjusted regressor | $X_{j\bullet}=X_{j}-\vec{\bf X}_{-j}^{T}{\boldsymbol{\beta}}_{-j\bullet}$ | ${\boldsymbol{\beta}}_{-j\bullet}=\left\\{E\left(\vec{\bf X}_{-j}\vec{\bf X}_{-j}^{T}\right)\right\\}^{-1}E\left(\vec{\bf X}_{-j}X_{j}\right)$
|  | $\vec{\bf X}_{-j}=(1,X_{1},\ldots,X_{j-1},X_{j+1},\ldots,X_{p})^{T}$
Intercept model | $W_{j}=\beta_{j}+aU_{1}+\sum_{j^{\prime}=1}^{p}b_{j^{\prime}}U_{j^{\prime}+1}+\tilde{\delta}_{j}$ | $W_{j}=\frac{YX_{j\bullet}}{E(X_{j\bullet}^{2})},U_{1}=\frac{X_{j\bullet}}{E(X_{j\bullet}^{2})}$
|  |  $U_{j^{\prime}+1}=\frac{X_{j^{\prime}}X_{j\bullet}}{E(X_{j\bullet}^{2})}\text{ for }j^{\prime}\neq j$, $U_{j+1}=\frac{X_{j}X_{j\bullet}}{E(X_{j\bullet}^{2})}-1$,
|  |  $a,{\bf b}$ are $\theta_{2},\cdots,\theta_{p+1}$ defined
|  | by (2) for model (17),
|  |  $\tilde{\delta}_{j}$ is the remainder term
TI estimator | $\\{\hat{\boldsymbol{\beta}}_{TI}\\}_{j}=\bar{W}_{j}-\hat{a}\bar{U}_{1}-\sum_{j^{\prime}=1}^{p}\hat{b}_{j^{\prime}}\bar{U}_{j^{\prime}+1}$ |  $\bar{\cdot}$ is the mean over the labeled sample
|  |  $\hat{a},\hat{\bf b}$ are the LSE of model (17)
PI estimator | $\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}=\bar{\check{W}}_{j}-\hat{a}\bar{\check{U}}_{1}-\sum_{j^{\prime}=1}^{p}\hat{b}_{j^{\prime}}\bar{\check{U}}_{j^{\prime}+1}$ |  $\bar{\cdot}$ is the mean over the labeled sample
|  | $\check{W}_{j}=\frac{YX_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})},~{}\check{U}_{1}=\frac{X_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})}$
|  | $\check{U}_{j^{\prime}+1}=\frac{X_{j^{\prime}}X_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})}\text{ for }j^{\prime}\neq j,~{}\check{U}_{j+1}=\frac{X_{j}X_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})}-1$
|  |  $\hat{a},\hat{\bf b}$ are the LSE of model (17) with
|  |  $(\check{W}_{j},\check{U}_{1},\ldots,\check{U}_{p+1})$ replacing $({W}_{j},U_{1},\ldots,U_{p+1})$
Asymptotics | $\nu=\lim\frac{n}{n+m}$ |  $n$ ($m$) the labeled (unlabeled) sample size
|  ${\boldsymbol{\Sigma}}_{LSE}$, ${\boldsymbol{\Sigma}}_{PI}$, ${\boldsymbol{\Sigma}}_{TI}$ | asymptotic covariance matrix of LSE, PI, TI
|  | as in Theorem 2
Main results | ${\boldsymbol{\Sigma}}_{LSE}={\boldsymbol{\Sigma}}_{TI}+{\boldsymbol{\Sigma}}_{\rm diff}$ | ${\boldsymbol{\Sigma}}_{\rm diff}=Cov(\delta{\bf X}_{\bullet}-\tilde{\boldsymbol{\delta}})$
| ${\boldsymbol{\Sigma}}_{PI}={\boldsymbol{\Sigma}}_{TI}+\nu{\boldsymbol{\Sigma}}_{\rm diff}$ |

Let $\vec{\bf X}_{-j}=(1,X_{1},\ldots,X_{j-1},X_{j+1},\ldots,X_{p})^{T}$ and
let

${\boldsymbol{\beta}}_{-j\bullet}=\arg\min_{\tilde{\boldsymbol{\beta}}}E\left(X_{j}-\tilde{\boldsymbol{\beta}}^{T}\vec{\bf
X}_{-j}\right)^{2}=\left\\{E\left(\vec{\bf X}_{-j}\vec{\bf
X}_{-j}^{T}\right)\right\\}^{-1}E\left(\vec{\bf X}_{-j}X_{j}\right).$ (16)

Now, define $X_{j\bullet}=X_{j}-\vec{\bf
X}_{-j}^{T}{\boldsymbol{\beta}}_{-j\bullet}$. Each $\beta_{j}$ can be written
in the one dimensional form
$\beta_{j}=\frac{E(YX_{j\bullet})}{E(X_{j\bullet}^{2})}.$

The standard LSE can be viewed in a similar manner. Let ${\bf
Y}=(Y^{(1)},\ldots,Y^{(n)})^{T},~{}{\bf
X}_{j}=(X_{j}^{(1)},\ldots,X_{j}^{(n)})^{T}$ and let ${\bf X}_{-j}=({\bf
1},{\bf X}_{1},\ldots,{\bf X}_{j-1},{\bf X}_{j+1},\ldots,{\bf X}_{p})$. Define
also $\hat{\boldsymbol{\beta}}_{-j\bar{\bullet}}=\left\\{{\bf X}_{-j}^{T}{\bf
X}_{-j}\right\\}^{-1}{\bf X}_{-j}^{T}{\bf X}_{j}$ and ${\bf
X}_{j\bar{\bullet}}={\bf X}_{j}-{\bf
X}_{-j}\hat{\boldsymbol{\beta}}_{-j\bar{\bullet}}$; then,
$\\{{\hat{\boldsymbol{\beta}}}_{LSE}\\}_{j}=\frac{\langle\bf{Y},{\bf
X}_{j\bar{\bullet}}\rangle}{||{\bf X}_{j\bar{\bullet}}||^{2}}$. Recall that
$\bar{\cdot}$ denotes mean over the labeled sample; thus, $\bar{\bullet}$
denotes adjustments over the labeled sample, whereas $\bullet$ denotes
adjustments over the population.

The total information estimator is the intercept estimator of the regression
model in the population obtained by multiplying (15) by
$\frac{X_{j\bullet}}{E(X_{j\bullet}^{2})}$, that is,

$W_{j}=\beta_{j}+aU_{1}+\sum_{j^{\prime}=1}^{p}b_{j^{\prime}}U_{j^{\prime}+1}+\tilde{\delta}_{j},$
(17)

where

$W_{j}=\frac{YX_{j\bullet}}{E(X_{j\bullet}^{2})},~{}U_{1}=\frac{X_{j\bullet}}{E(X_{j\bullet}^{2})},~{}U_{j^{\prime}+1}=\frac{X_{j^{\prime}}X_{j\bullet}}{E(X_{j\bullet}^{2})}\text{
for }j^{\prime}\neq j,\text{ and
}U_{j+1}=\frac{X_{j}X_{j\bullet}}{E(X_{j\bullet}^{2})}-1$ (18)

are $Y,1,X_{1},\cdots,X_{p}$ multiplied by
$\frac{X_{j\bullet}}{E(X_{j\bullet}^{2})}$ and $a,{\bf b}$ are
$\theta_{2},\ldots,\theta_{p+1}$ defined by (2) for $W=W_{j}$ and ${\bf
U}=(U_{1},\ldots,U_{p+1})$. Since $E(U_{1})=\cdots=E(U_{p+1})=0$, then the
intercept is $E(W_{j})=\beta_{j}$. In setting $E(W_{j})=\beta_{j}$ and also
the expectations of the $U$’s to be zero, we exploited the knowledge of the
moments of ${\bf X}$, yielding higher efficiency as shown below. The remainder
term in (17),
$\tilde{\delta}_{j}=W_{j}-\beta_{j}-aU_{1}-\sum_{j^{\prime}=1}^{p}b_{j}^{\prime}U_{j^{\prime}+1}$
is orthogonal to $(U_{1},\ldots,U_{p+1})$. Notice that the vector
$(U_{1},\ldots,U_{p+1})$ depends on $j$ but this is suppressed in the
notation.

Specifically, we define the total information (TI) estimator to be

$\\{\hat{\boldsymbol{\beta}}_{TI}\\}_{j}=\bar{W}_{j}-\hat{a}\bar{U}_{1}-\sum_{j^{\prime}=1}^{p}\hat{b}_{j^{\prime}}\bar{U}_{j^{\prime}+1},$

where $\hat{a},\hat{\bf b}$ are the regular LSE of (17), and $\bar{\cdot}$
denotes the mean over the labeled sample with $n$ observations.

To define the partial information estimator let

${\boldsymbol{\beta}}_{-j\check{\bullet}}=\left\\{\check{E}\left(\vec{\bf
X}_{-j}\vec{\bf X}_{-j}^{T}\right)\right\\}^{-1}\check{E}\left(\vec{\bf
X}_{-j}X_{j}\right),$ (19)

where $\check{E}$ is the the empirical mean based on the full ${\bf X}$ sample
of size $n+m$. Now define $X_{j\check{\bullet}}=X_{j}-\left\\{\vec{\bf
X}_{-j}\right\\}^{T}{\boldsymbol{\beta}}_{-j\check{\bullet}}$ and

$\check{W}_{j}=\frac{YX_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})},~{}\check{U}_{1}=\frac{X_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})},~{}\check{U}_{j^{\prime}+1}=\frac{X_{j^{\prime}}X_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})}\text{
for }j^{\prime}\neq j,\text{ and
}\check{U}_{j+1}=\frac{X_{j}X_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})}.$

Here ${\check{\bullet}}$ denotes adjustments over the full $n+m$ sample. The
partial information (PI) estimator is

$\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}=\bar{\check{W}}_{j}-\hat{a}\bar{\check{U}}_{1}-\sum_{j^{\prime}=1}^{p}\hat{b}_{j^{\prime}}\bar{\check{U}}_{j^{\prime}+1},$

where $\hat{a},\hat{b}_{j^{\prime}}$ are the regular LSE of model (17) with
$(\check{W}_{j},\check{U}_{1},\ldots,\check{U}_{p+1})$ replacing
$({W}_{j},U_{1},\ldots,U_{p+1})$.

The following theorem establishes the asymptotic distribution of the estimates
and states conditions under which $\hat{\boldsymbol{\beta}}_{TI}$ and
$\hat{\boldsymbol{\beta}}_{PI}$, asymptotically dominate
$\hat{\boldsymbol{\beta}}_{LSE}$. The asymptotic distribution of
$\hat{\boldsymbol{\beta}}_{LSE}$ is already known (White, 1980) and it is
presented here so that the comparison to $\hat{\boldsymbol{\beta}}_{TI}$ and
$\hat{\boldsymbol{\beta}}_{PI}$ can be made.

######  Theorem 2.

  1. (i)

Suppose that the vector ${\bf
X}_{\bullet}=\left\\{\frac{X_{1\bullet}}{E(X_{1\bullet}^{2})},\ldots,\frac{X_{p\bullet}}{E(X_{p\bullet}^{2})}\right\\}$
is well defined (i.e., the projections in (16) exist and $E(X_{j\bullet})^{2}$
is positive and finite for $j=1,\ldots,p$) and that $\delta{\bf X}_{\bullet}$
has finite second moments, then,

$\sqrt{n}(\hat{\boldsymbol{\beta}}_{LSE}-{\boldsymbol{\beta}})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,{\boldsymbol{\Sigma}}_{LSE}),$

where ${\boldsymbol{\Sigma}}_{LSE}=Cov(\delta{\bf X}_{\bullet})$.

  2. (ii)

Suppose further that for each $j=1,\ldots,p$, the vector
$(W_{j},U_{1},U_{2},\ldots,U_{p+1})$ has finite second moments and that the
matrix $E({\bf U}_{j}{\bf U}_{j}^{T})$ is invertible for ${\bf
U}_{j}=\left(1,U_{1},U_{2},\ldots,U_{p+1}\right)^{T}$; then,

$\sqrt{n}(\hat{\boldsymbol{\beta}}_{TI}-{\boldsymbol{\beta}})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,{\boldsymbol{\Sigma}}_{TI}),$

where ${\boldsymbol{\Sigma}}_{TI}=Cov(\tilde{\boldsymbol{\delta}})$ for
$\tilde{\boldsymbol{\delta}}=(\tilde{\delta}_{1},\ldots,\tilde{\delta}_{p})$.
Furthermore, we have that
${\boldsymbol{\Sigma}}_{LSE}={\boldsymbol{\Sigma}}_{TI}+{\boldsymbol{\Sigma}}_{\rm
diff}$ where ${\boldsymbol{\Sigma}}_{\rm diff}=Cov({\bf V})$ for ${\bf
V}=\delta{\bf X}_{\bullet}-\tilde{\boldsymbol{\delta}}$.

  3. (iii)

Suppose further that $\lim\frac{n}{n+m}=\nu$, then,

$\sqrt{n}(\hat{\boldsymbol{\beta}}_{PI}-{\boldsymbol{\beta}})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,{\boldsymbol{\Sigma}_{PI}}),$

where
${\boldsymbol{\Sigma}}_{PI}={\boldsymbol{\Sigma}}_{TI}+\nu{\boldsymbol{\Sigma}}_{\rm
diff}$.

Therefore, if ${\boldsymbol{\Sigma}}_{\rm diff}$ is not the zero matrix, then
${\boldsymbol{\Sigma}}_{LSE}-{\boldsymbol{\Sigma}}_{TI}$ is positive definite,
and if further $\nu<1$, then
${\boldsymbol{\Sigma}}_{LSE}-{\boldsymbol{\Sigma}}_{PI}$ is also positive
definite.

######  Corollary 2.

Theorem 2 and Proposition 1 imply that
$R(\hat{\alpha}_{LSE},\hat{\boldsymbol{\beta}}_{LSE})-R(\hat{\alpha}_{PI},\hat{\boldsymbol{\beta}}_{PI})=\text{Trace}\\{(1-\nu){\boldsymbol{\Sigma}}_{\rm
diff}{\bf M}\\}$, where ${\bf M}=Cov({\bf X})$.

In short, Theorem 2 states that if the regression models (15) and (17) are
well defined and the residuals have finite second moments, then
${\boldsymbol{\Sigma}}_{LSE}-{\boldsymbol{\Sigma}}_{TI}$ is positive semi-
definite and it is strictly positive unless ${\bf V}\equiv 0$. If further, the
unlabeled sample size is not negligible, i.e. $\lim m/n>0$, then
${\boldsymbol{\Sigma}}_{LSE}-{\boldsymbol{\Sigma}}_{PI}$ is also positive
definite.

As in the one-dimensional case, $\hat{\boldsymbol{\beta}}_{TI}$ and
$\hat{\boldsymbol{\beta}}_{PI}$ improve over $\hat{\boldsymbol{\beta}}_{LSE}$
only when ${\boldsymbol{\Sigma}}_{\rm diff}$ is not the zero matrix or
equivalently that ${\bf V}$ is not zero. For each $j$, $Var(V_{j})$ measures
the difference between the original model (15) and the intercept model (17).
When

$E(X_{j\bullet}^{2}\delta)=E(X_{j\bullet}^{2}X_{1}\delta)=\cdots=E(X_{j\bullet}^{2}X_{p}\delta)=0$
(20)

then $\delta X_{j\bullet}\equiv\delta_{j}$ and the two models are essentially
the same, in which case $V_{j}\equiv 0$. Otherwise, if for some $j$ (20) does
not hold, then $\hat{\boldsymbol{\beta}}_{PI}$ improves over
$\hat{\boldsymbol{\beta}}_{LSE}$ when $\nu<1$. In the one-dimensional case,
(20) is equivalent to $E(X^{2}\delta)=E(X^{3}\delta)=0$, i.e., that $X^{2}$
and $X^{3}$ are uncorrelated with $\delta$. Generally, (20) implies that
certain linear combinations of $E(X_{j}X_{j^{\prime}}\delta)$ and
$E(X_{j}X_{j^{\prime}}X_{j^{\prime\prime}}\delta)$ are 0. Thus, roughly
speaking, when $E[Y|{\bf X}]$ is non-linear and the non-linear part is
correlated with several second or third moments of ${\bf X}$, then we expect
improvement of $\hat{\boldsymbol{\beta}}_{PI}$.

###  3.3 New methodology (summary)

We now provide a step-by-step description of the new methodology for
estimation of ${\boldsymbol{\beta}}$.

For each $j=1,\ldots,p$:

  1. 1.

(Adjust the regressors) Let ${\boldsymbol{\beta}}_{-j\check{\bullet}}$ as
defined in (19) and define
$X_{j\check{\bullet}}^{(i)}=X_{j}^{(i)}-\left\\{\vec{\bf
X}_{-j}^{(i)}\right\\}^{T}{\boldsymbol{\beta}}_{-j\check{\bullet}}$ for
$i=1,\ldots,n$.

  2. 2.

(Define the intercept model) Define for $i=1,\ldots,n$

$\check{W}_{j}^{(i)}=\frac{Y^{(i)}X_{j\check{\bullet}}^{(i)}}{\check{E}(X_{j\check{\bullet}}^{2})},~{}\check{U}_{1}^{(i)}=\frac{X_{j\check{\bullet}}^{(i)}}{\check{E}(X_{j\check{\bullet}}^{2})},~{}\check{U}_{j^{\prime}+1}^{(i)}=\frac{X_{j^{\prime}}^{(i)}X_{j\check{\bullet}}^{(i)}}{\check{E}(X_{j\check{\bullet}}^{2})}\text{
for }j^{\prime}\neq j,\text{ and
}\check{U}_{j+1}^{(i)}=\frac{X_{j}^{(i)}X_{j\check{\bullet}}^{(i)}}{\check{E}(X_{j\check{\bullet}}^{2})}-1.$

  3. 3.

(Define the intercept estimator) The partial information (PI) estimator is

$\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}=\bar{\check{W}}_{j}-\hat{a}\bar{\check{U}}_{1}-\sum_{j^{\prime}=1}^{p}\hat{b}_{j^{\prime}}\bar{\check{U}}_{j^{\prime}+1},$

where $\hat{a},\hat{b}_{j^{\prime}}$ are the regular LSE of the regression
model

$W_{j}^{(i)}=\beta_{j}+a\check{U}_{1}^{(i)}+\sum_{j^{\prime}=1}^{p}b_{j^{\prime}}\check{U}_{j^{\prime}+1}^{(i)}+\tilde{\delta}_{j}^{(i)},\text{for
}i=1,\ldots,n.$

The methodology is built from standard least squares procedures. An R code
that implements the algorithm and also computes estimates of the variance as
in Section 5.1 below, is available at the homepage of the first author.

##  4 A simulation study

We compare the performance of the Partial Information (PI) and Total
Information (TI) estimators against the Least Squares Estimate (LSE) across a
wide range of settings.

###  4.1 Toy example

We start by studying the following toy model

$Y=\alpha X^{2}+\beta X+\varepsilon=\alpha+\beta
X+\underbrace{\alpha(X^{2}-1)+\varepsilon}_{=\delta},$ (21)

where $X$ and $\varepsilon$ are i.i.d N(0,1). Under model (21), the linear and
non-linear part are determined by $\beta$ and $\alpha$ separately and the
linear coefficient $\beta$ does not affect the residual $\delta$. In this
case, it easy to calculate the asymptotic variance of the estimates explicitly

$\sigma^{2}_{LSE}=10\alpha^{2}+1,~{}\sigma^{2}_{TI}=6\alpha^{2}+1,~{}\sigma^{2}_{\rm
diff}=4\alpha^{2}.$

Considering the excess risk, then by Proposition 1 the ratio of excess risks
of $\hat{\beta}_{PI}$ ($\hat{\beta}_{TI}$, respectively) and
$\hat{\beta}_{LSE}$ converges to
$\frac{2\alpha^{2}+1+\sigma^{2}_{PI}}{2\alpha^{2}+1+\sigma^{2}_{LSE}}$
($\frac{2\alpha^{2}+1+\sigma^{2}_{TI}}{2\alpha^{2}+1+\sigma^{2}_{LSE}}$,
respectively). The limits equal 1 in the linear case ($\alpha=0$) and approach
0 when the non-linear part is dominant ($\alpha\to\infty$).

Figure 1 summarizes 10,000 simulations of model (21) with $\beta=1$ and
different values of $\alpha$, $n$; for the PI case we used $m=2n$. For all the
scenarios we found that $\hat{\beta}_{TI}$, $\hat{\beta}_{PI}$ have smaller
excess risk when $n$ is large enough. As $\alpha$ increases, the departure
from linearity is more significant and the ratio of the excess risks is
smaller. For small $n$, LSE is superior for all scenarios as it is a simpler
estimate. When the model is close to linear ($\alpha=1/4$), the new estimates
are better for $n\geq 200$ and for the other values of $\alpha$ the new
estimates are better for $n\geq 70$. For small values of $\alpha$ the limiting
excess risk is close to the actual risk even when $n\approx 300$ but this does
not hold true for larger values of $\alpha$. In short, we found that the new
estimates are better for large $n$, and are much better when the non-linear
part is significant.

(a) TI  (b) PI

Figure 1: Estimates of the excess risk ratio for the TI and PI estimates based
on 10,000 simulations of model (21) with $\beta=1$ and different values of
$\alpha$, $n$; for the PI case we used $m=2n$. A confidence interval based on
two standard deviations is also plotted. The horizontal lines represents the
limiting excess risk.

###  4.2 Simulation Parameters

We now compare our new estimates to LSE for a broader variety of scenarios.
Estimates of the Excess Risk Ratio (ERR) for PI and TI were obtained for every
combination of the following parameter choices:

$n$ (size of labeled dataset) | 12, 25, 50, 100, 250, 500
---|---
$m$ (size of unlabeled dataset) | $n,2n$
$p$ (number of predictors) | 1, 4
Distribution of $X$ | Gaussian, Lognormal, Exponential, Cubed Gaussian
Errors | Gaussian, $N(0,e^{2\|X\|})$
$E[Y|X]$ |  $X$, $e^{X}$, $X^{3}$, $\sqrt{X}$

Fixing each of the 768 parameters settings, many sample datasets were
generated, along with a large test dataset (n = 100,000). For each sample
dataset, PI, TI, and LSE were fit, (calculating TI using the PI method with
$m\geq 500n$.) Next, an estimate of expected MSE was obtained using the test
dataset. Along with the best linear fit for the test dataset (given by least
squares regression), these determine the Excess Risk Ratio for PI ($ERR_{PI}$)
as follows:

$ERR_{PI}=\frac{MSE_{PI}-MSE_{BLF}}{MSE_{LSE}-MSE_{BLF}},$

where
$MSE_{PI}$ is the expected MSE of the PI estimator (estimated from mean
performance on the test dataset),
$MSE_{LSE}$ is the expected MSE of the LSE estimator, and
$MSE_{BLF}$ is the MSE of the best linear fit of the test dataset.

An analagous calculation was performed to calculate $ERR_{TI}$. Additional
datasets were sampled to improve estimates of $MSE_{PI}$ and $MSE_{LSE}$,
until standard errors for $ERR_{PI}$ and $ERR_{TI}$ fell below $1\%$ (as
estimated by the delta method), or a maximum of $100,000$ sample datasets was
generated.

###  4.3 Simulation Results

Table 3 provides the proportions of scenarios where the PI estimate yields a
statistically significant smaller excess risk than the LSE and the proportion
where the opposite holds true, for different $p$ and $n$. Statistical
significance is measured by two standard errors. The results demonstrate that
the PI estimate outperforms the LSE across a wide range of scenarios for large
$n$, and the proportion increases with $n$. When $n=500$, LSE is significantly
better only for about 10% of the scenarios we studied. More detailed comments
follow:

Table 3: Proportions of scenarios where $ERR_{PI}$ is significantly smaller than 1 and where it is significantly larger than 1.  |  p n | 12 | 25 | 50 | 100 | 250 | 500
---|---|---|---|---|---|---|---
$ERR_{PI}<1$ significantly | 1 | 0.297 | 0.344 | 0.469 | 0.594 | 0.656 | 0.609
4 | 0.047 | 0.234 | 0.391 | 0.438 | 0.594 | 0.547
$ERR_{PI}>1$ significantly | 1 | 0.344 | 0.234 | 0.172 | 0.125 | 0.125 | 0.109
4 | 0.578 | 0.391 | 0.234 | 0.188 | 0.188 | 0.125

  1. 1.

A proposed rule of thumb when $p=1$: As indicated by Table 3 across a wide
range of parameter settings, PI performs almost as well as LSE when $n\geq
100$, $m\geq 100$ and linear $E[Y|X]$, and better than the LSE for nonlinear
$E[Y|X]$ (that is, $ERR<1$). Furthermore, PI’s underperformance is mild when
$E[Y|X]$ is linear. Thus, assuming that the scenarios we studied are
representative sample of reality, we recommend using the PI method when the
following conditions are satisfied, as a rule of thumb:

    1. (a)

$n\geq 100$

    2. (b)

$m\geq 100$

    3. (c)

$E[Y|X]$ could be non-linear

  2. 2.

If $p=1$, $n>100$, $m>100$, and $E[Y|X]$ is one of the non-linear functions
tested, then both PI and TI appear to outperform LSE; furthermore, the margin
of outperformance widens as $n$ increases (See Figure 2 (a), and Figure 2 (b),
respectively). In particular, PI does much better for $E[Y|{X}]=X^{3}$ and
$e^{X}$, as compared to $\sqrt{X}$. Intuitively, this results from the larger
non-linear moments of $X^{3}$ and $e^{X}$, causing ${\boldsymbol{\Sigma}}_{\rm
diff}$ to be large (as defined in Theorem 2, part ii)

  3. 3.

Increasing $m$, the size of the unlabeled dataset, improves performance of PI.
Indeed, comparing PI where $m=2n$ in Graph X against TI ($m\geq 500n$) in
Graph Y, one observes that greatly increasing the pool of unlabeled samples
can reduce Excess Risk, in some cases on the order of $5-10$ per cent (Compare
Figures 2 (a) and (b)). But, a word of caution: if $n$ is small, then LSE may
still outperform both PI and TI, even when $m$ is large and $E[Y|X]$ is non-
linear.

(a) PI  (b) TI

Figure 2: Estimates of the excess risk ratio for the TI and PI estimates, for
several possible functions of $E[Y|X]$. In PI estimates (a), $m=2n$. In TI
estimates (b), the PI method is used with $m\geq 500n$. In all cases, X is
Gaussian, errors are Gaussian, and $p=1$. Confidence intervals are $\pm 2$
standard errors.

  4. 4.

If $n$ is small, or if $E[Y|X]$ is linear, then LSE performs better than TI
and PI. However, $ERR\to 1$ as $n\to\infty$, holding other parameters fixed
(See Figures 2 (a) and (b)). This result is consistent with the asymptotic
agreement of PI, TI, and LSE when $E[Y|X]$ is linear (See Theorem 2)

  5. 5.

Extra heteroskedastic noise appears to decrease the relative advantage of PI.
See Figure 3 (a).

  6. 6.

Somewhat similarly, increasing the number of parameters to $4$ (where Y is the
sum of its relationships with each $X_{i}$) also decreases the relative
advantage of PI. See Figure 3 (b).

  7. 7.

If $X$ is not Gaussian, PI still does well, especially for large $n$. See
Figure 3 (c) for an example where $X$ is exponential. (Note that standard
errors may be understated, especially in the case $E[Y|X]=\exp(X)$, partly due
to the large moments of $Y$.)

(a) Heteroskedastic Errors (b) $p=4$ (c) $X$ Exponential

Figure 3: Estimates of the excess risk ratio of PI estimates, across several
possible functions of $E[Y|X]$, with tweaks to the previous settings of Figure
2 (a): these were $m=2n$, $X$ Gaussian, $p=1$, and Guassian errors. The
following is changed in each plot: Figure (a): Errors are heteroscedastic,
distributed $N(0,e^{2\|x\|})$. Figure (b): $p=4$, where $E[Y|X]$ is the sum of
the given relationship in each $X$-variable. Figure (c): $X$ is distributed
exponentially. Confidence intervals are $\pm 2$ standard errors.

##  5 Estimating the standard errors of the estimates

In this section we study estimation of the asymptotic covariance matrix. We
suggest three estimates for ${\boldsymbol{\beta}}_{PI}$ in Section 5.1 and
study their performance by simulation (Section 5.2). Section 5.3 deals with
the estimation of the prediction error (excess risk).

###  5.1 Three estimates of ${\boldsymbol{\Sigma}}_{PI}$

We proved that under certain conditions,

$\sqrt{n}(\hat{\boldsymbol{\beta}}_{LSE}-{\boldsymbol{\beta}})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,{\boldsymbol{\Sigma}}_{LSE}),~{}\sqrt{n}(\hat{\boldsymbol{\beta}}_{PI}-{\boldsymbol{\beta}})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,{\boldsymbol{\Sigma}}_{PI})\text{ and,
}\sqrt{n}(\hat{\boldsymbol{\beta}}_{TI}-{\boldsymbol{\beta}})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,{\boldsymbol{\Sigma}}_{TI}).$

In these expressions the (asymptotic) variability of the estimates are
determined by the matrices ${\boldsymbol{\Sigma}}_{LSE}$,
${\boldsymbol{\Sigma}}_{PI}$ and ${\boldsymbol{\Sigma}}_{TI}$. These matrices
are unknown and therefore need to be estimated. In this section we consider
only estimation of ${\boldsymbol{\Sigma}}_{LSE}$ and
${\boldsymbol{\Sigma}}_{PI}$. Estimation of the asymptotic variance yields
prediction error estimates or $L_{2}$ risk in a standard fashion. In section
5.3 we discuss estimation of the excess risk.

For ${\boldsymbol{\Sigma}}_{LSE}$ we follow Buja et al. (2016) who show that
the standard estimates are inconsistent when non-linearity is present and
suggest two alternative estimates for the variance. The parametric sandwich
estimator for the asymptotic variance of ${\hat{\boldsymbol{\beta}}}_{LSE}$ is

$\hat{\rm
AV}_{parm}({\hat{\boldsymbol{\beta}}}_{LSE})=\left(\mathbb{X}^{T}{\mathbb{X}}\right)^{-1}\mathbb{X}^{T}\hat{D}{\mathbb{X}}\left(\mathbb{X}^{T}{\mathbb{X}}\right)^{-1},$

where $\mathbb{X}$ is the design matrix and $\hat{D}$ is a diagonal matrix
with $\hat{D}_{i}=\hat{\delta}_{i}^{2}$ the standard residual estimator. The
second estimate is derived from a pairs bootstrap where a pair $({\bf
X}^{*},Y^{*})$ is sampled $N_{BS}$ times from the empirical joint distribution
of $({\bf X},Y)$, yielding a sample (under the empirical measure) of size
$N_{BS}$ of ${\hat{\boldsymbol{\beta}}}_{LSE}$. The resulting estimate of the
variance is denoted by $\hat{\rm AV}_{BS}({\hat{\boldsymbol{\beta}}}_{LSE})$.

For estimating the variance of ${\boldsymbol{\Sigma}}_{PI}$, we consider three
estimates described as follows:

  1. 1.

Parametric: The asymptotic variances of ${\hat{\boldsymbol{\beta}}}_{LSE}$ and
${\hat{\boldsymbol{\beta}}}_{PI}$ are

${\rm
AV}({\hat{\boldsymbol{\beta}}}_{LSE})=Cov(\tilde{\boldsymbol{\delta}})+Cov({\bf
V}),~{}{\rm
AV}({\hat{\boldsymbol{\beta}}}_{PI})=Cov(\tilde{\boldsymbol{\delta}})+\nu
Cov({\bf V}).$ (22)

To estimate the variance of $Cov(\tilde{\boldsymbol{\delta}})$ notice that for
$j=1,\ldots,p$ and $i=1,\ldots,n$,

$\check{\tilde{\delta}}_{j}^{(i)}=\check{W}^{(i)}-\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}-\hat{a}\check{U}_{1}^{(i)}-\sum_{j^{\prime}=1}^{p}\hat{b}_{j^{\prime}}\check{U}_{j^{\prime}+1}^{(i)}=\tilde{\delta}_{j}^{(i)}+o_{p}(1/\sqrt{n}),$

and hence,
$\overline{Cov}(\check{\tilde{\boldsymbol{\delta}}})\stackrel{{\scriptstyle
p}}{{\longrightarrow}}Cov(\tilde{\boldsymbol{\delta}})$, where
$\overline{Cov}$ is the empirical covariance based on the labeled $n$ sample.
Therefore, using (22) we obtain that a consistent estimator to ${\rm
AV}({\hat{\boldsymbol{\beta}}}_{PI})$ is

$\left\\{\hat{\rm
AV}_{parm}({\hat{\boldsymbol{\beta}}}_{PI})\right\\}_{j,j^{\prime}}=\overline{Cov}(\check{\tilde{\boldsymbol{\delta}}}_{j},\check{\tilde{\boldsymbol{\delta}}}_{j^{\prime}})+\nu\left(\left\\{\hat{\rm
AV}_{parm}(\\{{\hat{\boldsymbol{\beta}}}_{LSE}\\})\right\\}_{j,j^{\prime}}-\overline{Cov}(\check{\tilde{\boldsymbol{\delta}}}_{j},\check{\tilde{\boldsymbol{\delta}}}_{j^{\prime}})\right),$

for $j\neq j^{\prime}$ and for $j=1,\ldots,p$,

$\left\\{\hat{\rm
AV}_{parm}({\hat{\boldsymbol{\beta}}}_{PI})\right\\}_{j,j}=\frac{1}{n}\sum_{i=1}^{n}\left\\{\check{\tilde{\delta}}_{j}^{(i)}\right\\}^{2}+\nu\max\left(\left\\{\hat{\rm
AV}_{parm}(\\{{\hat{\boldsymbol{\beta}}}_{LSE}\\})\right\\}_{j,j}-\frac{1}{n}\sum_{i=1}^{n}\left\\{\check{\tilde{\delta}}_{j}^{(i)}\right\\}^{2},0\right).$

  2. 2.

Bootstrap: Let $\left\\{({\bf X}^{(i)}_{*},Y^{(i)}_{*})\right\\}_{i=1}^{n}$ be
a sample from the empirical distribution of $({\bf X},Y)$ (based on the
labeled observations) and $\\{{\bf X}^{(i)}_{*}\\}_{i=n+1}^{n+m}$ be a sample
from the empirical distribution of ${\bf X}$ (based on the unlabeled
observations). Thus,
${\hat{\boldsymbol{\beta}}}_{PI}^{*}={\hat{\boldsymbol{\beta}}}_{PI}^{*}\left(\left\\{({\bf
X}^{(i)}_{*},Y^{(i)}_{*})\right\\}_{i=1}^{n},\\{{\bf
X}^{(i)}_{*}\\}_{i=n+1}^{n+m}\right)$ is a sample under the empirical measure
of ${\hat{\boldsymbol{\beta}}}_{PI}$, which can yield an estimate of the
variance of the estimator, denoted by $\hat{\rm
AV}_{BS}({\hat{\boldsymbol{\beta}}}_{PI})$.

  3. 3.

Variance bootstrap: In the proof of Theorem 2 we showed that for
$j=1,\ldots,p$,

$\\{\hat{\boldsymbol{\beta}}_{LSE}\\}_{j}=\beta_{j}+\frac{\frac{1}{n}\sum_{i=1}^{n}{X}_{j{\bullet}}^{(i)}\delta^{(i)}}{E(X_{j\bullet}^{2})}+o_{p}(1/\sqrt{n})\text{
and
}\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}=\beta_{j}+\frac{1}{n}\sum_{i=1}^{n}\tilde{\delta}^{(i)}+\frac{1}{n+m}\sum_{i=1}^{n+m}V_{j}^{(i)}+o_{p}(1/\sqrt{n}).$

Hence, the definition of the vector ${\bf V}$ implies that

$\hat{\boldsymbol{\beta}}_{LSE}-\hat{\boldsymbol{\beta}}_{PI}=\frac{1}{n}\sum_{i=1}^{n}-{\bf
V}^{(i)}+\frac{1}{n+m}\sum_{i=1}^{n+m}{\bf V}^{(i)}+o_{p}(1/\sqrt{n})\\\
=\frac{1}{n}\sum_{i=1}^{n}{\bf
V}^{(i)}\left(\frac{n}{n+m}-1\right)+\frac{1}{n+m}\sum_{i=n+1}^{n+m}{\bf
V}^{(i)}+o_{p}(1/\sqrt{n}).$

Therefore, dropping the smaller order terms we obtain

$nCov(\hat{\boldsymbol{\beta}}_{LSE}-\hat{\boldsymbol{\beta}}_{PI})\approx
Cov({\bf V})\left(1-\frac{n}{n+m}\right)^{2}+\frac{nm}{(n+m)^{2}}Cov({\bf
V})\\\ =Cov({\bf V})\left(1-\frac{n}{n+m}\right)\to Cov({\bf
V})\left(1-\nu\right)={\boldsymbol{\Sigma}}_{LSE}-{\boldsymbol{\Sigma}}_{PI}.$

Thus, a a consistent estimator of the difference
${\boldsymbol{\Sigma}}_{LSE}-{\boldsymbol{\Sigma}}_{PI}=:{\boldsymbol{\Delta}}$
can be obtained by bootstrapping
$nCov({\hat{\boldsymbol{\beta}}}_{PI}^{*}-{\hat{\boldsymbol{\beta}}}_{LSE}^{*})$.
We denote this estimator by $\widehat{\boldsymbol{\Delta}}$. A consistent
estimate to $AV({\hat{\boldsymbol{\beta}}}_{PI})$ is

$\hat{\rm AV}_{VBS}({\hat{\boldsymbol{\beta}}}_{PI})=\hat{\rm
AV}_{BS}({\hat{\boldsymbol{\beta}}}_{LSE})-\widehat{\boldsymbol{\Delta}}.$

This estimator has the advantage over the previous suggestions that it is
always smaller than the estimated variance of the LSE (in the sense that the
difference is positive definite), in a similar fashion to the asymptotic
distribution.

In the next section we compare the different estimates through a simulation
study.

###  5.2 Simulations for the estimates of the variance

We now investigate the estimates of the variance under the toy model (21). We
study the performance of the estimates for $\alpha=1,1/2,1/4,1/8$, where small
$\alpha$ corresponds to little non-linearity. We considered $n=250$ and
$m=1500$ similar to the numbers of Section 6. For each $\alpha$ we repeated
the simulation 1000 times and also $N_{BS}=1000$.

The simulations are summarized in Table 4 and in Figure 4. We find that the
bootstrap estimate (henceforth BS) is more variable than the parametric
(henceforth PARM) and the variance bootstrap (henceforth VBS). PARM and VBS
are comparable but for small $\alpha$’s for estimating the difference, the
histogram of VBS is narrower around the true value. When $\alpha=1/8$, the
estimated difference under BS (PARM) is negative 32.9% (34.7%) of the
simulations, while for VBS it is never negative. However, since VBS cannot be
negative, the mean of the simulations is upwards biased. In short, we find
that PARM and VBS outperform BS and VBS has the advantage of never being
smaller than $\hat{AV}(\hat{\beta}_{LSE})$.

Table 4: The mean (std) of the estimates of the asymptotic variance $AV({\hat{\beta}}_{PI})$ and the difference $AV({\hat{\beta}}_{LSE})-AV({\hat{\beta}}_{PI})$.  | Variance | Difference
---|---|---
$\alpha$ | Bootstrap | Parametric | Variance BS | true | Bootstrap | Parametric | Variance BS | true
1 | 7.65 (3.71) | 6.86 (3.11) | 6.58 (2.40) | 7 | 2.63 (2.53) | 3.51 (1.53) | 3.70 (2.13) | 3.43
1/2 | 2.65 (0.97) | 2.42 (0.82) | 2.36 (0.70) | 2.5 | 0.64 (0.53) | 0.88 (0.45) | 0.93 (0.46) | 0.86
1/4 | 1.43 (0.40) | 1.33 (0.33) | 1.31 (0.31) | 1.37 | 0.15 (0.18) | 0.24 (0.22) | 0.27 (0.16) | 0.21
1/8 | 1.11 (0.23) | 1.06 (0.20) | 1.03 (0.19) | 1.09 | 0.02 (0.07) | 0.06 (0.17) | 0.10 (0.06) | 0.05

(a) $\alpha=1/4$

(b) $\alpha=1/8$

Figure 4: Histograms of the estimates (Bootstrap, Parametric, Variance BS) of
the asymptotic variance $AV({\hat{\beta}}_{PI})$ and the difference
$AV({\hat{\beta}}_{LSE})-AV({\hat{\beta}}_{PI})$. The true asymptotic value is
illustrated by the vertical gray line.

###  5.3 Estimating the prediction error

The three estimates of ${\mathbb{\Sigma}}_{PI}$ can be used to estimate the
prediction error in a standard fashion. Here, we describe how the variance
bootstrap estimate together with the result of Proposition 1 provide an
estimate of the excess risk.

By Proposition 1, the excess risk of $\hat{\boldsymbol{\beta}}_{LSE}$ can be
approximated by

$E\Big{[}Y-E(Y)-\sum_{j=1}^{p}{\beta}_{j}\\{{X}_{j}-E(X_{j})\\}\Big{]}^{2}+{\rm
Trace}({\bf M}{\boldsymbol{\Sigma}}_{LSE})=Var(Y)+{\boldsymbol{\beta}}^{T}{\bf
M}{\boldsymbol{\beta}}-2\sum_{j=1}^{p}\beta_{j}Cov(Y,X_{j})+{\rm Trace}({\bf
M}{\boldsymbol{\Sigma}}_{LSE})$

and hence can be consistently estimated by

$\hat{\sigma}^{2}_{Y}+\hat{\boldsymbol{\beta}}_{LSE}^{T}\hat{\bf
M}{\boldsymbol{\beta}}_{LSE}-2\sum_{j=1}^{p}\\{\hat{\boldsymbol{\beta}}_{LSE}\\}_{j}\overline{Cov}(Y,X_{j})+{\rm
Trace}\\{\hat{\bf M}\hat{AV}_{BS}({\boldsymbol{\Sigma}}_{LSE})\\},$

where $\hat{\sigma}^{2}_{Y}=\frac{1}{n}\sum_{i=1}^{n}(Y^{(i)}-\bar{Y})^{2}$,
$\hat{\bf M}=\frac{1}{n+m}\sum_{i=1}^{n+m}({\bf X}^{(i)}-\bar{\bf
X})\left\\{{\bf X}^{(i)}-\bar{\bf X}\right\\}^{T}$ and
$\overline{Cov}(Y,X_{j})=\frac{1}{n}\sum_{i=1}^{n}(Y^{(i)}-\bar{Y})(X^{(i)}_{j}-\bar{X}_{j})$
is the empirical covariance. Therefore, a consistent estimate to the excess
risk ratio is

$\widehat{ERR}_{PI}=\frac{\hat{\sigma}^{2}_{Y}+\hat{\boldsymbol{\beta}}_{LSE}^{T}\hat{\bf
M}{\boldsymbol{\beta}}_{LSE}-2\sum_{j=1}^{p}\\{\hat{\boldsymbol{\beta}}_{LSE}\\}_{j}\overline{Cov}(Y,X_{j})+{\rm
Trace}\\{\hat{\bf
M}(\hat{AV}_{BS}({\boldsymbol{\Sigma}}_{LSE})-\hat{\boldsymbol{\Delta}})\\}}{\hat{\sigma}^{2}_{Y}+\hat{\boldsymbol{\beta}}_{LSE}^{T}\hat{\bf
M}{\boldsymbol{\beta}}_{LSE}-2\sum_{j=1}^{p}\\{\hat{\boldsymbol{\beta}}_{LSE}\\}_{j}\overline{Cov}(Y,X_{j})+{\rm
Trace}\\{\hat{\bf M}\hat{AV}_{BS}({\boldsymbol{\Sigma}}_{LSE})\\}},$ (23)

where $\hat{\boldsymbol{\Delta}}$ is defined by the Variance bootstrap.

##  6 Inferring homeless population

We now consider the Los Angeles homeless data set (Kriegler and Berk, 2010).
Our purpose here is not to carefully analyze this data set but rather to
demonstrate our new method and to compare it to the standard LSE. For our
analysis, we used the same covariates as in Kriegler and Berk (2010); see this
reference for more details about the data.

There are about 2000 census tracts in the Los Angeles county, and a sample was
conducted in order to infer the homeless population. The sample consisted of
two parts. First, tracts believed to have large numbers of homeless people
were pre-selected and visited; there are about 240 such tracts. Second, a
sample of about 260 tracts was chosen from the remainder by a random sampling
technique and the homeless population was counted, leaving about 1500 tracts
to be imputed.

Leaving aside the preselected tracts this is almost exactly our semi-
supervised setting among the remaining population. We have $n=261$ and
$m=1536$. The major difference from the setting in our introduction here is
that sampling from the remaining population is without replacement and hence
is not an i.i.d random sample. Since we have not discussed this type of
sampling in detail, we will ignore this issue in the subsequent discussion
which is for illustrative purposes. (The difference in sampling schemes does
not affect the values of the various estimates, but does impact their standard
errors.)

We compare the LSE of the random sample and the new estimate (PI) where the
supervised part is the random sample (of size 261) and the unsupervised part
consists of the 1536 tracts to be imputed. The resulting estimates, as well as
their standard errors (SE) are given in Table 5. The SE of the LSE was
computed using pairs bootstrap ($N_{BS}=10^{5}$), and the SE of PI was
computed using the above variance bootstrap method (the SE of the intercept of
PI was estimated by pairs bootstrap).

We find that the SE’s of PI are smaller by 15%-40% than LSE in the predictors
MedianInc, PercVacant, PercCommercial, PercIndustrial. These predictors also
demonstrated a discrepancy between the SE from linear models to bootstrap SE
as reported in Buja et al. (2016) Section 2, indicating non-linearity in these
predictors (a discrepancy was observed also in the predictor PercResidential,
but the improvement in this predictor is relatively small). In summary, we
found that the additional information on the distribution of the predictors in
the unsupervised data provides estimates that are more accurate (smaller
variance).

Consider the prediction problem of estimating the homeless population for the
1500 tracts. An estimate of the excess risk ratio, using (23), is
$\widehat{ERR}_{PI}=0.904$, i.e., improvement of about 10%. An estimate of the
differences of prediction errors as in (10) and (11) is
$\text{Trace}(\hat{\boldsymbol{\Delta}}\hat{\bf M})/n=15.2$, while an estimate
of the prediction error of the mean is

$\frac{1}{m}\sum_{i=1}^{m}\vec{\bf
X}_{n+i}^{T}\hat{\boldsymbol{\Sigma}}_{LSE}\vec{\bf X}_{n+i}/n=42.8.$

Since 15.2/42.8 =0.355 the improvement in prediction risk is about 36%. When
considering the prediction of $Y$ itself as in (11) then the improvement is
much smaller, about 1.5%. The reason for this difference is that the
predictors we considered are all quite weak, and hence most of prediction
error comes from the variability in the distribution of $Y$ given ${\bf X}$
rather than the variability in estimating the $\beta$’s.

Table 5: The estimates $\hat{\boldsymbol{\beta}}_{LSE}$, $\hat{\boldsymbol{\beta}}_{PI}$ and their standard errors.  | $\hat{\beta}_{LSE}$ | $\hat{SE}_{LSE}$ | $\hat{\beta}_{PI}$ | $\hat{SE}_{PI}$ | $\frac{{SE}_{LSE}}{{SE}_{PI}}$
---|---|---|---|---|---
Intercept | 13.138 | 11.822 | 14.758 | 11.184 | 0.927
MedianInc ($K) | -0.065 | 0.056 | -0.080 | 0.044 | 1.284
PercVacant | 1.449 | 0.707 | 1.583 | 0.514 | 1.374
PercMinority | 0.060 | 0.105 | 0.063 | 0.099 | 1.058
PercResidential | -0.070 | 0.100 | -0.078 | 0.094 | 1.072
PercCommercial | 0.446 | 0.354 | 0.335 | 0.300 | 1.179
PercIndustrial | 0.101 | 0.188 | 0.202 | 0.143 | 1.317

##  7 Discussion

In this work we showed that ${\boldsymbol{\beta}}$ can be better estimated
when additional information on the predictors ${\bf X}$ is available. The key
idea is to replace the regression model

$Y=\alpha+\beta_{1}X_{1}+\cdots+\beta_{p}X_{p}+\delta,$ (24)

with $p$ regressions of the form

$W_{j}=\beta_{j}+aU_{1}+\sum_{j^{\prime}=1}^{p}b_{j^{\prime}}U_{j^{\prime}+1}+\tilde{\delta}_{j},$
(25)

for $j=1,\ldots,p$, where $W_{j},U_{1},\ldots,U_{p+1}$ is defined in (18).
Regression (25) can be used only when the first and second moments of ${\bf
X}$ are known. We showed that the intercept estimator of (25) is a better
estimate of $\beta_{j}$, in terms of smaller asymptotic variance, than the
standrd LSE of (24). Furthermore, we showed that even if the second moments of
${\bf X}$ are not known exactly then improvement can be made over the standard
LSE provided that the unlabeled sample is non-negligible with respect to the
labeled sample, i.e., $m/n$ is bounded away from 0.

Here we considered the classical framework where $p$ is fixed and
$n\to\infty$. Zhang et al. (2016) study the semi-supervised mean estimation
problem where $p,n$ both go to infinity but $p$ grows relatively slowly, i.e.,
$p=O(\sqrt{n})$. The findings of Zhang et al. (2016) indicate that our results
can be extended to the latter case. However, our method obviously requires
that $p<<n$ and for high-dimensional regression models where $p>n$, a
different approach is needed. High-dimensional missspecified regression models
are the topic of a recent manuscript by Bühlmann and van de Geer (2015).

In this work we did not study variable selection and considered the covariates
as given. However, improved asymptotic performance may be achieved by
including functions of the existing covariates as additional new covariates.
Section 3 of Zhang et al. (2016) discusses how this can be done for the
problem of semi-supervised estimation of a mean. On the other hand, when the
labeled sample, $n$, is small, then it may be better not to include all
covariates. We hope to study the problem of variable selection in this context
in a future research.

A possible extension of these results is to generalized linear models, where
the label variable ($Y$) takes a small number of values and the regression
model reduces to a classification problem. Our hope is that our method can be
extended to these cases and improvement can be made over naive classifiers
that consider only the labeled data.

In summary, in this work we considered the framework where the best linear
predictor is of interest even if $E(Y|{\bf X})$ is not linear. Under this
framework, additional information on the distribution of ${\bf X}$ is useful
to construct better estimates than the standard estimates. We believe that
this framework is of practical importance and also can lead to interesting
statistical issues, which are yet to be studied.

##  8 Proofs

######  Lemma 1.

We have that

$EL(\tilde{\alpha},\tilde{\boldsymbol{\beta}})=E\Big{[}\bar{Y}-\sum_{j=1}^{p}\tilde{\beta}_{j}\\{{\bar{X}}_{j}-E(X_{j})\\}-E(Y)\Big{]}^{2}+E\left\\{({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
M}({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})\right\\}.$ (26)

Proof.
First,

$R(\tilde{\alpha},\tilde{\boldsymbol{\beta}})=E\left(Y^{*}-\tilde{\alpha}-\tilde{\boldsymbol{\beta}}^{T}{\bf
X}^{*}\right)^{2}-E\left(Y^{*}-{\alpha}-{\boldsymbol{\beta}}^{T}{\bf
X}^{*}\right)^{2}\\\
=E\left\\{(\alpha-\tilde{\alpha})+({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
X}^{*}\right\\}\left\\{2Y^{*}-(\alpha+\hat{\alpha})-({\boldsymbol{\beta}}+\tilde{\boldsymbol{\beta}})^{T}{\bf
X}^{*}\right\\}\\\
=E\left\\{(\alpha-\tilde{\alpha})+({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
X}^{*}\right\\}^{2}+2E\left\\{(\alpha-\hat{\alpha})+({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
X}^{*}\right\\}\left\\{Y^{*}-\alpha-{\boldsymbol{\beta}}^{T}{\bf
X}^{*}\right\\}\\\
=E\left\\{(\alpha-\tilde{\alpha})+({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
X}^{*}\right\\}^{2},$

where the last equality holds true since
$Y^{*}-\alpha-{\boldsymbol{\beta}}^{T}{\bf X}^{*}$ is orthogonal to $X^{*},1$.
We can further write

$E\left\\{(\alpha-\tilde{\alpha})+({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
X}^{*}\right\\}^{2}=E\left[(\alpha-\tilde{\alpha})+({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}E({\bf
X})+({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}\\{{\bf X}^{*}-E({\bf
X})\\}\right]^{2}\\\
=E\left\\{(\alpha-\tilde{\alpha})+({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}E({\bf
X})\right\\}^{2}+E\\{({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
M}({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})\\},$ (27)

For the first term in (27) notice that
$\alpha=E(Y)-\sum_{j=1}^{p}\beta_{j}E(X_{j})$ and
$\tilde{\alpha}=\bar{Y}-\sum_{j=1}^{p}\tilde{\beta}_{j}{\bar{X}}_{j}$;
therefore,

$E\left\\{(\alpha-\tilde{\alpha})+({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}E({\bf
X})\right\\}^{2}=E\Big{[}\bar{Y}-\sum_{j=1}^{p}\tilde{\beta}_{j}\\{{\bar{X}}_{j}-E(X_{j})\\}-E(Y)\Big{]}^{2},$

and (26) is established. ∎

Proof of Proposition 1.
Since for any numbers $c,d,B$ we have that
$\min(c+d,B)\leq\min(c,B)+\min(d,B)$ and
$\min(c+d,B)\geq\min(c,B/2)+\min(d,B/2)$, then

${R}(\tilde{\alpha},\tilde{\boldsymbol{\beta}})=\lim_{B\to\infty}\lim_{n\to\infty}E\min\\{n\tilde{L}(\tilde{\alpha},\tilde{\boldsymbol{\beta}}),B\\}=\lim_{B\to\infty}\lim_{n\to\infty}E\min\left(n\Big{[}\bar{Y}-\sum_{j=1}^{p}\tilde{\beta}_{j}\\{{\bar{X}}_{j}-E(X_{j})\\}-E(Y)\Big{]}^{2},B\right)\\\
+\lim_{B\to\infty}\lim_{n\to\infty}E\min\left(n\left\\{({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
M}({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})\right\\},B\right).$ (28)

We start with the first summand in (28). Since
$\sqrt{n}(\tilde{\boldsymbol{\beta}}-{\boldsymbol{\beta}})$ is asymptotically
normal then

$\bar{Y}-\sum_{j=1}^{p}\tilde{\beta}_{j}\\{\bar{X}_{j}-E(X_{j})\\}-E(Y)-\Big{[}\bar{Y}-\sum_{j=1}^{p}{\beta}_{j}\\{\bar{X}_{j}-E(X_{j})\\}-E(Y)\Big{]}=\sum_{j=1}^{p}({\beta}_{j}-\tilde{\beta}_{j})\\{\bar{X}_{j}-E(X_{j})\\}=o_{p}(1/\sqrt{n}).$

Therefore, by the continuous mapping theorem and since the random variables
are bounded, then for any $B$,

$\lim_{n\to\infty}E\min\left(n\Big{[}\bar{Y}-\sum_{j=1}^{p}\tilde{\beta}_{j}\\{{\bar{X}}_{j}-E(X_{j})\\}-E(Y)\Big{]}^{2},B\right)\\\
=\lim_{n\to\infty}E\min\left(n\Big{[}\bar{Y}-\sum_{j=1}^{p}{\beta}_{j}\\{{\bar{X}}_{j}-E(X_{j})\\}-E(Y)\Big{]}^{2},B\right)\\\
=E\min\left(Z_{1}^{2},B\right),$ (29)

where $Z_{1}$ is normal with mean zero and variance
$\Big{[}{Y}-\sum_{j=1}^{p}{\beta}_{j}\\{{X}_{j}-E(X_{j})\\}-E(Y)\Big{]}^{2}$.

For the second summand in (28), the asymptotic normality
$\sqrt{n}(\tilde{\boldsymbol{\beta}}-{\boldsymbol{\beta}})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,{\boldsymbol{\Sigma}})$ implies that for any $B$,

$\lim_{n\to\infty}E\min\left(n\left\\{({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})^{T}{\bf
M}({\boldsymbol{\beta}}-\tilde{\boldsymbol{\beta}})\right\\},B\right)=E\min\left\\{{\bf
Z}_{2}^{T}{\bf Z}_{2},B\right\\},$ (30)

where ${\bf Z}_{2}$ is a normal vector with mean zero and variance matrix
${\bf M}^{1/2}{\boldsymbol{\Sigma}}{\bf M}^{1/2}$. Taking limits as
$B\to\infty$ in (29) and (30) completes the proof of the proposition since
${\rm Trace}({\bf M}^{1/2}{\boldsymbol{\Sigma}}{\bf M}^{1/2})={\rm Trace}({\bf
M}{\boldsymbol{\Sigma}})$. ∎

We will not prove Theorem 1 since it is a special case of Theorem 2.

Proof of Theorem 2.
Part (i) We have that (Buja et al., 2016)
$\\{\hat{\boldsymbol{\beta}}_{LSE}\\}_{j}=\beta_{j}+\frac{\langle{\boldsymbol{\delta}},{\bf
X}_{j\bar{\bullet}}\rangle}{||{\bf X}_{j\bar{\bullet}}||^{2}}$. Furthermore,

$\\{\hat{\boldsymbol{\beta}}_{LSE}\\}_{j}=\beta_{j}+\frac{\langle{\boldsymbol{\delta}},{\bf
X}_{j\bar{\bullet}}\rangle}{||{\bf
X}_{j\bar{\bullet}}||^{2}}=\beta_{j}+\frac{\frac{1}{n}\sum_{i=1}^{n}{X}_{j\bar{\bullet}}^{(i)}\delta^{(i)}}{\frac{1}{n}\sum_{i=1}^{n}\left\\{{X}_{j\bar{\bullet}}^{(i)}\right\\}^{2}}=\beta_{j}+\frac{\frac{1}{n}\sum_{i=1}^{n}{X}_{j{\bullet}}^{(i)}\delta^{(i)}}{E(X_{j\bullet}^{2})}+o_{p}(1/\sqrt{n}),$

since

$\frac{\sqrt{n}}{n}\sum_{i=1}^{n}\delta^{(i)}({X}_{j\bar{\bullet}}^{(i)}-{X}_{j\bullet}^{(i)})=\frac{\sqrt{n}}{n}\sum_{i=1}^{n}\delta^{(i)}\left\\{\vec{\bf
X}_{-j}^{(i)}\right\\}^{T}\left(\hat{\boldsymbol{\beta}}_{-j\bar{\bullet}}-{\boldsymbol{\beta}}_{-j{\bullet}}\right)\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0,\text{ and }\frac{1}{n}||{\bf
X}_{j\bar{\bullet}}||^{2}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}E(X_{j\bullet}^{2}).$

The results now follows from the multivariate CLT and Slutsky’s theorem.
Part (ii) For each $j$,

$\\{\hat{\boldsymbol{\beta}}_{TI}\\}_{j}=\bar{W}_{j}-\hat{a}\bar{U}_{1}-\sum_{j^{\prime}=1}^{p}\hat{b}_{j^{\prime}}\bar{U}_{j^{\prime}+1}=\bar{W}_{j}-{a}\bar{U}_{1}-\sum_{j^{\prime}=1}^{p}{b}_{j^{\prime}}\bar{U}_{j^{\prime}+1}+o_{p}(1/\sqrt{n})\\\
=\beta_{j}+\frac{1}{n}\sum_{i=1}^{n}\tilde{\delta}_{j}^{(i)}+o_{p}(1/\sqrt{n}).$
(31)

Again, the multivariate CLT and Slutsky’s theorem imply the result.

To prove the furthermore part, notice that for each $j$,

$V_{j}=\frac{\delta
X_{j\bullet}}{E(X_{j\bullet}^{2})}-\tilde{\delta}_{j}=\frac{\left(Y-\alpha-\sum_{j^{\prime}=1}^{p}\beta_{j^{\prime}}X_{j^{\prime}}\right)X_{j\bullet}}{E(X_{j\bullet}^{2})}-\frac{YX_{j\bullet}}{E(X_{j\bullet}^{2})}+\beta_{j}-b_{j}+a\frac{X_{j\bullet}}{E(X_{j\bullet}^{2})}+\sum_{j^{\prime}=1}^{p}b_{j^{\prime}}\frac{X_{j^{\prime}}X_{j\bullet}}{E(X_{j\bullet}^{2})}\\\
=\beta_{j}-{b_{j}}+(a-\alpha)\frac{X_{j\bullet}}{E(X_{j\bullet}^{2})}+\sum_{j^{\prime}=1}^{p}(b_{j^{\prime}}-\beta_{j^{\prime}})\frac{X_{j^{\prime}}X_{j\bullet}}{E(X_{j\bullet}^{2})}=(a-\alpha)U_{1}+\sum_{j^{\prime}=1}^{p}(b_{j^{\prime}}-\beta_{j^{\prime}})U_{j^{\prime}+1}$
(32)

is a linear function of $U_{1},\ldots,U_{p+1}$ and hence is orthogonal to
$\tilde{\delta}_{j}$, i.e., $EV_{j}\tilde{\delta}_{j}=0$ for all $j$.
Therefore,

${\boldsymbol{\Sigma}}_{LSE}=Cov({\bf
X}_{\bullet}\delta)=Cov\left\\{\tilde{\boldsymbol{\delta}}+({\bf
X}_{\bullet}\delta-\tilde{\boldsymbol{\delta}})\right\\}=Cov(\tilde{\boldsymbol{\delta}}+{\bf
V})=Cov(\tilde{\boldsymbol{\delta}})+Cov({\bf
V})={\boldsymbol{\Sigma}}_{TI}+{\boldsymbol{\Sigma}}_{\rm diff}.$

Part (iii) We show that

$\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}-\\{\hat{\boldsymbol{\beta}}_{TI}\\}_{j}=\frac{1}{n+m}\sum_{i=1}^{n+m}V_{j}^{(i)}+o_{p}(1/\sqrt{n}).$
(33)

Hence,
$\sqrt{n}(\hat{\boldsymbol{\beta}}_{PI}-\hat{\boldsymbol{\beta}}_{TI})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,\nu{\boldsymbol{\Sigma}}_{\rm diff})$. Since ${\bf
V}$ is orthogonal to $\tilde{\boldsymbol{\delta}}$, then

$\sqrt{n}(\hat{\boldsymbol{\beta}}_{PI}-{\boldsymbol{\beta}})=\sqrt{n}(\hat{\boldsymbol{\beta}}_{TI}-{\boldsymbol{\beta}})+\sqrt{n}(\hat{\boldsymbol{\beta}}_{PI}-\hat{\boldsymbol{\beta}}_{TI})\stackrel{{\scriptstyle{\cal
D}}}{{\longrightarrow}}N(0,{\boldsymbol{\Sigma}}_{TI}+\nu{\boldsymbol{\Sigma}}_{\rm
diff}).$

We now prove (33) by a somewhat lengthy calculation. We have that

$\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}=\bar{\check{W}}_{j}-\hat{a}\bar{\check{U}}_{1}-\sum_{j^{\prime}=1}^{p}\hat{b}_{j^{\prime}}\bar{\check{U}}_{j^{\prime}+1}=\bar{\hat{W}}_{j}-{a}\bar{\check{U}}_{1}-\sum_{j^{\prime}=1}^{p}{b}_{j^{\prime}}\bar{\check{U}}_{j^{\prime}+1}+o_{p}(1/\sqrt{n}).$
(34)

Therefore, (34) and (31) yield

$\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}-\\{\hat{\boldsymbol{\beta}}_{TI}\\}_{j}=\bar{\check{W}}_{j}-\bar{W}_{j}-{a}(\bar{\check{U}}_{1}-\bar{U}_{1})-\sum_{j^{\prime}=1}^{p}{b}_{j^{\prime}}(\bar{\check{U}}_{j^{\prime}+1}-\bar{U}_{j^{\prime}+1})+o_{p}(1/\sqrt{n}).$

We have that

$\bar{\check{W}}_{j}-\bar{W}_{j}={\alpha}(\bar{\check{U}}_{1}-\bar{U}_{1})+\sum_{j^{\prime}=1}^{p}{\beta}_{j^{\prime}}(\bar{\check{U}}_{j^{\prime}+1}-\bar{U}_{j^{\prime}+1})+o_{p}(1/\sqrt{n}),$

and hence,

$\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}-\\{\hat{\boldsymbol{\beta}}_{TI}\\}_{j}=(\alpha-a)(\bar{\check{U}}_{1}-\bar{U}_{1})+\sum_{j^{\prime}=1}^{p}({\beta}_{j^{\prime}}-b_{j^{\prime}})(\bar{\check{U}}_{j^{\prime}+1}-\bar{U}_{j^{\prime}+1})+o_{p}(1/\sqrt{n}).$
(35)

We now consider the summands in (35). We start with the first summand

$\bar{\check{U}}_{1}-\bar{U}=\frac{\frac{1}{n}\sum_{i=1}^{n}X_{j\check{\bullet}}^{(i)}}{\check{E}(X_{j\check{\bullet}}^{2})}-\frac{\frac{1}{n}\sum_{i=1}^{n}X_{j{\bullet}}^{(i)}}{{E}(X_{j{\bullet}})^{2}}=\frac{1}{n}\sum_{i=1}^{n}X_{j{\bullet}}^{(i)}\left(\frac{1}{\check{E}(X_{j\check{\bullet}}^{2})}-\frac{1}{{E}(X_{j{\bullet}})^{2}}\right)+\frac{({\boldsymbol{\beta}}_{-j{\bullet}}-{\boldsymbol{\beta}}_{-j\check{\bullet}})^{T}\frac{1}{n}\sum_{i=1}^{n}\vec{\bf
X}_{-j}^{(i)}}{\check{E}(X_{j\check{\bullet}}^{2})}\\\
=\frac{({\boldsymbol{\beta}}_{-j{\bullet}}-{\boldsymbol{\beta}}_{-j\check{\bullet}})^{T}\frac{1}{n}\sum_{i=1}^{n}\vec{\bf
X}_{-j}^{(i)}}{{E}(X_{j{\bullet}})^{2}}+o_{p}(1/\sqrt{n})\\\
=\frac{({\boldsymbol{\beta}}_{-j{\bullet}}-{\boldsymbol{\beta}}_{-j\check{\bullet}})^{T}\frac{1}{n+m}\sum_{i=1}^{n+m}\vec{\bf
X}_{-j}^{(i)}}{{E}(X_{j{\bullet}})^{2}}+\frac{({\boldsymbol{\beta}}_{-j{\bullet}}-{\boldsymbol{\beta}}_{-j\check{\bullet}})^{T}\left(\frac{1}{n}\sum_{i=1}^{n}\vec{\bf
X}_{-j}^{(i)}-\frac{1}{n+m}\sum_{i=1}^{n+m}\vec{\bf
X}_{-j}^{(i)}\right)}{{E}(X_{j{\bullet}})^{2}}+o_{p}(1/\sqrt{n})\\\
=\frac{({\boldsymbol{\beta}}_{-j{\bullet}}-{\boldsymbol{\beta}}_{-j\check{\bullet}})^{T}\frac{1}{n+m}\sum_{i=1}^{n+m}\vec{\bf
X}_{-j}^{(i)}}{{E}(X_{j{\bullet}})^{2}}+o_{p}(1/\sqrt{n}).$

Since

$0=\check{E}(X_{j\check{\bullet}})=\check{E}(X_{j}-{\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j})\Longrightarrow\check{E}(X_{j})=\check{E}({\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j})$

then,

$\check{E}(X_{j{\bullet}})=\check{E}(X_{j}-{\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j})=\check{E}({\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j})-\check{E}({\boldsymbol{\beta}}_{-j{\bullet}}^{T}\vec{\bf
X}_{-j})={({\boldsymbol{\beta}}_{-j\check{\bullet}}-{\boldsymbol{\beta}}_{-j{\bullet}})^{T}\frac{1}{n+m}\sum_{i=1}^{n+m}\vec{\bf
X}_{-j}^{(i)}},$

and therefore,

$\bar{\check{U}}_{1}-\bar{U}_{1}=\frac{-\frac{1}{n+m}\sum_{i=1}^{n+m}X_{j{\bullet}}^{(i)}}{{E}(X_{j{\bullet}})^{2}}+o_{p}(1/\sqrt{n}).$
(36)

We now consider the $j^{\prime}$-th term in (35) for $j^{\prime}\neq j$; then,
$E(X_{j^{\prime}}X_{j{\bullet}})=0$. We have

$\bar{\check{U}}_{j^{\prime}+1}-\bar{U}_{j^{\prime}+1}=\frac{\frac{1}{n}\sum_{i=1}^{n}X^{(i)}_{j^{\prime}}X^{(i)}_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})}-\frac{\frac{1}{n}\sum_{i=1}^{n}X^{(i)}_{j^{\prime}}X^{(i)}_{j{\bullet}}}{{E}(X_{j{\bullet}})^{2}}\\\
=\frac{1}{n}\sum_{i=1}^{n}X^{(i)}_{j^{\prime}}X^{(i)}_{j{\bullet}}\left\\{\frac{1}{\check{E}(X_{j\check{\bullet}}^{2})}-\frac{1}{{E}(X_{j{\bullet}})^{2}}\right\\}+\frac{({\boldsymbol{\beta}}_{-j{\bullet}}-{\boldsymbol{\beta}}_{-j\check{\bullet}})^{T}\frac{1}{n}\sum_{i=1}^{n}X_{j^{\prime}}^{(i)}\vec{\bf
X}_{-j}^{(i)}}{\check{E}(X_{j\check{\bullet}}^{2})}\\\
=\frac{({\boldsymbol{\beta}}_{-j{\bullet}}-{\boldsymbol{\beta}}_{-j\check{\bullet}})^{T}\frac{1}{n+m}\sum_{i=1}^{n+m}X_{j^{\prime}}^{(i)}\vec{\bf
X}_{-j}^{(i)}}{{E}(X_{j{\bullet}})^{2}}+o_{p}(1/\sqrt{n}).$

Since

$0=\check{E}(X_{j\check{\bullet}}X_{j^{\prime}})=\check{E}\\{(X_{j}-{\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j})X_{j^{\prime}}\\}\Longrightarrow\check{E}(X_{j}X_{j^{\prime}})=\check{E}({\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j}X_{j^{\prime}})$

then,

$\check{E}(X_{j{\bullet}}X_{j^{\prime}})=\check{E}\\{(X_{j}-{\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j})X_{j^{\prime}}\\}=\check{E}({\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j}X_{j^{\prime}})-\check{E}({\boldsymbol{\beta}}_{-j{\bullet}}^{T}\vec{\bf
X}_{-j}X_{j^{\prime}})={({\boldsymbol{\beta}}_{-j\check{\bullet}}-{\boldsymbol{\beta}}_{-j{\bullet}})^{T}\frac{1}{n+m}\sum_{i=1}^{n+m}X_{j^{\prime}}^{(i)}\vec{\bf
X}_{-j}^{(i)}},$

and therefore,

$\bar{\check{U}}_{j^{\prime}+1}-\bar{U}_{j^{\prime}+1}=\frac{-\frac{1}{n+m}\sum_{i=1}^{n+m}X_{j{\bullet}}^{(i)}X_{j^{\prime}}^{(i)}}{{E}(X_{j{\bullet}})^{2}}+o_{p}(1/\sqrt{n}).$
(37)

For the $j$-th term in (35), we have that
$E(X_{j}X_{j{\bullet}})=E(X_{j{\bullet}})^{2}$; hence,

$\bar{\check{U}}_{j+1}-\bar{U}_{j+1}=\frac{\frac{1}{n}\sum_{i=1}^{n}X^{(i)}_{j}X^{(i)}_{j\check{\bullet}}}{\check{E}(X_{j\check{\bullet}}^{2})}-\frac{\frac{1}{n}\sum_{i=1}^{n}X^{(i)}_{j}X^{(i)}_{j{\bullet}}}{{E}(X_{j{\bullet}})^{2}}\\\
=\left\\{\frac{1}{n}\sum_{i=1}^{n}X^{(i)}_{j^{\prime}}X^{(i)}_{j{\bullet}}-E(X_{j{\bullet}})^{2}\right\\}\left\\{\frac{1}{\check{E}(X_{j\check{\bullet}}^{2})}-\frac{1}{{E}(X_{j{\bullet}})^{2}}\right\\}+\frac{({\boldsymbol{\beta}}_{-j{\bullet}}-{\boldsymbol{\beta}}_{-j\check{\bullet}})^{T}\frac{1}{n}\sum_{i=1}^{n}X_{j}^{(i)}\vec{\bf
X}_{-j}^{(i)}+E(X_{j{\bullet}})^{2}-\check{E}(X_{j\check{\bullet}}^{2})}{\check{E}(X_{j\check{\bullet}}^{2})}\\\
=\frac{({\boldsymbol{\beta}}_{-j{\bullet}}-{\boldsymbol{\beta}}_{-j\check{\bullet}})^{T}\frac{1}{n}\sum_{i=1}^{n}X_{j}^{(i)}\vec{\bf
X}_{-j}^{(i)}+E(X_{j{\bullet}})^{2}-\check{E}(X_{j\check{\bullet}}^{2})}{{E}(X_{j{\bullet}})^{2}}+o_{p}(1/\sqrt{n}).$

Since

$\check{E}(X_{j\check{\bullet}}^{2})=\check{E}(X_{j\check{\bullet}}X_{j})=\check{E}\\{(X_{j}-{\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j})X_{j}\\}\Longrightarrow\check{E}(X_{j}^{2})=\check{E}({\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j}X_{j})+\check{E}(X_{j\check{\bullet}}^{2})$

then,

$\check{E}(X_{j{\bullet}}X_{j})=\check{E}\\{(X_{j}-{\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j})X_{j}\\}=\check{E}({\boldsymbol{\beta}}_{-j\check{\bullet}}^{T}\vec{\bf
X}_{-j}X_{j})+\check{E}(X_{j\check{\bullet}}^{2})-\check{E}({\boldsymbol{\beta}}_{-j{\bullet}}^{T}\vec{\bf
X}_{-j}X_{j})\\\
={({\boldsymbol{\beta}}_{-j\check{\bullet}}-{\boldsymbol{\beta}}_{-j{\bullet}})^{T}\frac{1}{n+m}\sum_{i=1}^{n+m}X_{j}^{(i)}\vec{\bf
X}_{-j}^{(i)}}+\check{E}(X_{j\check{\bullet}}^{2}).$

Thus, for the $j$-th term we obtain,

$\bar{\check{U}}_{j+1}-\bar{U}_{j+1}=\frac{-\frac{1}{n+m}\sum_{i=1}^{n+m}X_{j{\bullet}}^{(i)}X_{j}^{(i)}}{{E}(X_{j{\bullet}})^{2}}+1+o_{p}(1/\sqrt{n}).$
(38)

Tracking back to (35) through (36),(37), and (38), we obtain

$\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}-\\{\hat{\boldsymbol{\beta}}_{TI}\\}_{j}=(\beta_{j}-{b_{j}})+(a-\alpha)\frac{\frac{1}{n+m}\sum_{i=1}^{n+m}X_{j{\bullet}}^{(i)}}{E(X_{j\bullet}^{2})}+\sum_{j^{\prime}=1}^{p}(b_{j^{\prime}}-\beta_{j^{\prime}})\frac{\frac{1}{n+m}\sum_{i=1}^{n+m}X_{j{\bullet}}^{(i)}X_{j^{\prime}}^{(i)}}{E(X_{j\bullet}^{2})}+o_{p}(1/\sqrt{n}).$

Hence, (32) implies that

$\\{\hat{\boldsymbol{\beta}}_{PI}\\}_{j}-\\{\hat{\boldsymbol{\beta}}_{TI}\\}_{j}=\frac{1}{n+m}\sum_{i=1}^{n+m}V_{j}^{(i)}+o_{p}(1/\sqrt{n}),$

and (33) is established. ∎

## References

  * Brown (1990) Brown, L.D. (1990). An Ancillarity Paradox Which Appears in Multiple Linear Regression. Annals of Statistics, 18, 471–493.
  * Bühlmann and van de Geer (2015) Bühlmann, P., van de Geer, S. (2015). High-dimensional inference in misspecified linear models. Electronic Journal of Statistics, 9, 1449–1473.
  * Buja et al. (2016) Buja, A., Berk, R., Brown, L., George, E., Pitkin, E., Traskin, M., Zhao, L., Zhang, K. (2016). Models as Approximations: How Random Predictors and Model Violations Invalidate Classical Inference in Regression. Submitted; arXiv:1404.1578.
  * Chakrabortty and Cai (2016) Chakrabortty A., Cai T. (2016). Efficient and Adaptive Linear Regression in Semi-Supervised settings. Submitted.
  * Cochran (1977) Cochran W.G. (1977). Sampling Techniques (3rd ed.), New York: John Wiley.
  * Johnson and Zhang (2008) Johnson, R., Zhang, T. (2008). Graph-based semi-supervised learning and spectral kernel design. IEEE Transactions on Information Theory. 54, 275–288.
  * Kriegler and Berk (2010) Kriegler, B., Berk, R. (2010). Small area estimation of the homeless in Los Angeles: An application of cost-sensitive stochastic gradient boosting. Annals of Applied Statistics, 4, 1234–1255.
  * Lafferty and Wasserman (2008) Lafferty, J., Wasserman, L. (2008). Statistical analysis of semi-supervised regression. Advances in Neural Information Processing Systems, 20, 801–-808.
  * Liang et al. (2008) Liang, F., Mukherjee, S., West, M. (2007). The use of unlabeled data in predictive modeling. Statistical Science, 22, 189-205.
  * Le Cam and Yang (1990) Le Cam, L., Yang, G.L. (1990). Asymptotics in Statistics : Some Basic Concepts. New York: Springer.
  * Little (1992) Little, R.J.A. (1992). Regression With Missing X’s: A Review. Journal of the American Statistical Association, 87, 1227–1237.
  * Merz et al. (1992) Merz, C.J., Clair, D.C.St., Bond, W.E (1992). Semi-supervised adaptive resonance theory (smart2). International Joint Conference on Neural Networks, 3, 851–-856.
  * Wang et al. (2009) Wang, J., Shen, X., Pan, W. (2009). On efficient large margin semisupervised learning: Method and theory. The Journal of Machine Learning Research, 10, 719-742.
  * White (1980) White, H. (1980). Using Least Squares to Approximate Unknown Regression Functions. International Economic Review, 21, 149–170.
  * Zhang et al. (2016) Zhang, A., Brown L.D., Cai, T. (2016). Semi-supervised Inference: General Theory and estimation of means. Submitted.
  * Zhou (2005) Zhou, Z.H., Li, M. (2005). Semi-supervised regression with co-training. International Joint Conference on Artificial Intelligence, 5, 908–913.
  * Zhu (2005) Zhu, X., (2005). Semi-Supervised Learning Literature Survey. Computer Sciences, University of Wisconsin-Madison, 1530.