# On Measuring Model Complexity in Heteroscedastic
Linear Regression

Bo Luan  The Ohio State University  Yoonkyung Lee  The Ohio State University
Yunzhang Zhu  The Ohio State University

###### Abstract

Heteroscedasticity is common in real world applications and is often handled
by incorporating case weights into a modeling procedure. Intuitively, models
fitted with different weight schemes would have a different level of
complexity depending on how well the weights match the inverse of error
variances. However, existing statistical theories on model complexity, also
known as model degrees of freedom, were primarily established under the
assumption of equal error variances. In this work, we focus on linear
regression procedures and seek to extend the existing measures to a
heteroscedastic setting. Our analysis of the weighted least squares method
reveals some interesting properties of the extended measures. In particular,
we find that they depend on both the weights used for model fitting and those
for model evaluation. Moreover, modeling heteroscedastic data with optimal
weights generally results in fewer degrees of freedom than with equal weights,
and the size of reduction depends on the unevenness of error variance. This
provides additional insights into weighted modeling procedures that are useful
in risk estimation and model selection.

Keywords: Heteroscedasticity, Model complexity, Model degrees of freedom,
Weighted least squares

##  1 Introduction

Heteroscedasticity is common in real world data and needs careful
consideration for modeling. Mishandling it may lead to inefficient estimates
for model parameters and unreliable predictions. To deal with
heteroscedasticity, a commonly used strategy is to incorporate case weights
into a modeling procedure. If the error variances associated with cases were
known, ideally we would specify the weights inversely proportional to the
error variances according to the normal likelihood. However, such case-
specific variability is usually unknown a priori, and weights need to be
estimated. This generally results in various weight schemes depending on how
they are estimated. For models with different weight schemes, we are
interested in whether and how their complexity measures vary. A deeper
understanding of this issue will offer insights into the impact of weight
schemes on model performance in relation to the true error variance and
further provide guidance on weight selection and ensuing model selection.

In statistics, the term “model degrees of freedom” is used as a measure of the
complexity or flexibility of a model (Tibshirani and Hastie,, 1987; Ye,, 1998;
Efron,, 2004). For least squares linear regression as an example, the model
degrees of freedom is defined as the number of free parameters in the model,
which coincides with the number of predictors. More generally, the degrees of
freedom of a model can be defined through the expected downward bias of the
training error as an estimate of the prediction error. This expected bias is
also known as the expected optimism (Efron,, 1986). Intuitively, models fitted
with different weight schemes would differ in complexity, since both training
and prediction errors are affected by case weights. However, to the best of
our knowledge, no complexity measure has been formally established for
weighted modeling procedures.

In this paper, we examine the expected optimism for in-sample prediction as
well as out-of-sample prediction with case weights. The classical optimism
theory (Efron,, 1986) takes in-sample prediction error as the primary target.
For a standard regression setting with constant error variance, the model
degrees of freedom is defined through the covariances between observed
responses and their fitted values. A shortcoming of this notion of model
degrees of freedom is that it fails to account for the flexibility of a model
at new feature values, which could be a major issue with high dimensional
data. To overcome this issue, Luan et al., (2021) proposed a different version
of model degrees of freedom, called the predictive model degrees of freedom,
based on the out-of-sample prediction error. Both versions of model degrees of
freedom have been derived under the assumption of equal error variance for
regression models, and they do not directly apply to modeling procedures with
case weights.

We aim to extend the optimism theory and the two notions of model degrees of
freedom to a more general setting with heteroscedastic data for linear
regression procedures. First, we assume that the error variance is a function
of covariates rather than constant. This extension has also been considered in
Anscombe, (1961), Bickel, (1978), Box and Hill, (1974), Carroll and Ruppert,
(1981), for example. Second, unlike the existing theories that evaluate a
model with squared loss, we use a weighted squared loss as an evaluation
metric for a given weight scheme. Using these weights for evaluation, we
define analogous expected optimism and model degrees of freedom.

As an example of linear regression procedures, we examine the two extended
measures of model degrees of freedom for the weighted least squares method in
two scenarios. In the first scenario of subset regression, we study their
monotonicity as the subset size increases. According to these two measures,
the weighted least squares model has smaller degrees of freedom than the
ordinary least squares model of the same size, and the size of reduction in
the degrees of freedom depends on the unevenness of the error variances. In
the second scenario, we fix the set of covariates as well as weights for
evaluation and examine the impact of the weight scheme for fitting a weighted
least squares model to data, which we call weights for model fitting. From
sensitivity analyses of degrees of freedom around the optimal weights, we gain
some useful insights into the relationships between the two model complexity
measures and case weights for model fitting.

Using the two extended measures of model degrees of freedom, we propose
estimators for both in-sample and out-of-sample prediction errors in a
heteroscedastic setting. For the former, our estimator is a direct extension
of Mallows’s $C_{p}$ (Mallows,, 1973), whereas for the latter, the proposed
estimator is an adaptation of the risk estimator developed in Luan et al.,
(2021) to heteroscedastic data. In practice, the predictive model degrees of
freedom is usually unknown due to its dependency on the underlying data
distribution. We thereby propose methods that can estimate the predictive
model degrees of freedom efficiently and use them for estimation of the out-
of-sample prediction error.

These proposed estimators are then assessed through a series of simulation
studies and real data analysis. First, we evaluate our out-of-sample risk
estimator in terms of prediction error estimation and model selection. We find
that it has smaller variance than the leave-one-out and five-fold cross
validation strategies. In addition, it tends to favor slightly simpler models
without compromising the predictive performance too much. Second, we evaluate
the impact of weight misspecification on model selection for both in-sample
and out-of-sample predictions. In general, ignoring heteroscedasticity may
lead to misidenfication of the model size and inflation in prediction risk.

This paper is organized as follows. In Section 2, we first review the two
notions of model degrees of freedom defined under the assumption of equal
error variance. We then propose their extended versions in a heteroscedastic
setting in Section 3. We examine properties of the extended measures for the
weighted least squares method in Section 4, and propose approaches to
estimation of the predictive model degrees of freedom and out-of-sample
prediction error in Section 5. In Section 6, we evaluate the proposed
estimators through simulation studies and real data analysis. We then conclude
the paper in Section 7. For conciseness, all proofs will be deferred to
Appendix.

##  2 Review of Model Complexity Theories

In this section, we review the existing theories for the complexity of
regression models. To start with, we first introduce the notation that will be
used in this work. Consider a regression setting with a continuous response
$Y\in\mathbb{R}$. For a given set of covariates $X=(X_{1},\ldots,X_{d})$,
assume that $Y=\mu(X;\bm{\beta})+\varepsilon$, where $\mu$ is the mean
regression function with unknown parameter $\bm{\beta}$, and $\varepsilon$ is
some random error such that $\mathrm{E}(\varepsilon|X)=0$ and
$\mathrm{Var}(\varepsilon|X)=\sigma_{\varepsilon}^{2}>0$. Here, we allow
$X_{1},\ldots,X_{d}$ to be either numerical or categorical.

Assume that a dataset of size $n$ is used for model fitting. Let
$\mathbf{y}=(y_{1},\ldots,y_{n})^{\top}$ be the vector of responses and
$\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})^{\top}\in\mathbb{R}^{n\times
m}$ be the full design matrix that includes all $d$ covariates. We encode a
categorical variable with $L$ levels by $L-1$ binary columns in $\mathbf{X}$.
Thus, $m>d$ if there is at least one categorical variable with more than two
levels.

For brevity, we abbreviate $\mu(\mathbf{x}_{i};\bm{\beta})$ by $\mu_{i}$ and
write $\bm{\mu}=(\mu_{1},\ldots,\mu_{n})^{\top}$. For a given modeling
procedure, define $\hat{\mu}_{i}$ to be the fitted value of $y_{i}$ and write
$\hat{\bm{\mu}}=(\hat{\mu}_{i},\ldots,\hat{\mu}_{n})^{\top}$. With the
training data $(\mathbf{X},\mathbf{y})$, $\hat{\bm{\mu}}$ is obtained by
minimizing the (unweighted) training error

$\mathrm{ErrT}_{\mathbf{X},\mathbf{y}}=\frac{1}{n}\sum_{i=1}^{n}(y_{i}-\hat{\mu}_{i})^{2}=\frac{1}{n}\|\mathbf{y}-\hat{\bm{\mu}}\|^{2}.$

###  2.1 Classical Model Degrees of Freedom

We first review the classical optimism theory presented in Efron, (1986, 2004)
that focuses on in-sample prediction error. Let $\tilde{y}_{i}$ be an
independent copy of $y_{i}$ given $\mathbf{x}_{i}$. The (unweighted) in-sample
prediction error is defined as

$\mathrm{ErrF}_{\mathbf{X},\mathbf{y}}=\frac{1}{n}\mathrm{E}(\|\tilde{\mathbf{y}}-\hat{\bm{\mu}}\|^{2}|\mathbf{X},\mathbf{y}),$

where the expectation is with respect to
$\tilde{\mathbf{y}}=(\tilde{y}_{1},\ldots,\tilde{y}_{n})^{\top}$.

The expected optimism, denoted by $\mathrm{OptF}_{\mathbf{X}}$, is the
expected difference between the prediction error
$\mathrm{ErrF}_{\mathbf{X},\mathbf{y}}$ and training error
$\mathrm{ErrT}_{\mathbf{X},\mathbf{y}}$. The expectation is taken with respect
to $\mathbf{y}$ given $\mathbf{X}$. Efron, (1986) showed that

$\mathrm{OptF}_{\mathbf{X}}\coloneqq\mathrm{E}(\mathrm{ErrF}_{\mathbf{X},\mathbf{y}}-\mathrm{ErrT}_{\mathbf{X},\mathbf{y}}|\mathbf{X})=\frac{2}{n}\sum_{i=1}^{n}\mathrm{Cov}(y_{i},\hat{\mu}_{i}|\mathbf{X}).$
(2.1)

Following this covariance representation of the expected optimism, Efron,
(2004) then formally defined the model degrees of freedom of a fitted model
$\hat{\bm{\mu}}$ as

$\mathrm{df}_{\rm
F}=\sum_{i=1}^{n}\frac{\mathrm{Cov}(y_{i},\hat{\mu}_{i}|\mathbf{X})}{\sigma_{\varepsilon}^{2}}.$
(2.2)

In particular, for a linear modeling procedure
$\hat{\bm{\mu}}=\mathbf{H}\mathbf{y}$, where
$\mathbf{H}=(h_{ij})\in\mathbb{R}^{n\times n}$ depends only on $\mathbf{X}$,
it is easy to show that
$\mathrm{Cov}(y_{i},\hat{\mu}_{i}|\mathbf{X})=\sigma_{\varepsilon}^{2}h_{ii}$
and $\mathrm{df}_{\rm F}={\rm tr}(\mathbf{H})$. Here, $\mathbf{H}$ is usually
called the hat matrix.

When the input space is large, the in-sample prediction error may be limited
in measuring how well a modeling procedure will generalize at new feature
values that arise in the future as in many prediction-oriented applications.
As a result, the classical model degrees of freedom $\mathrm{df}_{\rm F}$
needs some adjustment when prediction beyond the training data is concerned.
This motivates the use of out-of-sample prediction error for an alternative
definition of model degrees of freedom, which we will discuss next.

###  2.2 Predictive Model Degrees of Freedom

In this section, we review the model complexity theory for out-of-sample
prediction proposed by Luan et al., (2021). Let
$(\mathbf{x}_{\ast},\varepsilon_{\ast})$ be an independent copy of
$(\mathbf{x}_{i},\varepsilon_{i})$ and define
$\mu_{\ast}=\mu(\mathbf{x}_{\ast};\bm{\beta})$ and
$y_{\ast}=\mu_{\ast}+\varepsilon_{\ast}$. Let $\hat{\mu}_{\ast}$ be a
prediction of $y_{\ast}$. The (unweighted) out-of-sample prediction error is
defined as

$\mathrm{ErrR}_{\mathbf{X},\mathbf{y}}=\mathrm{E}[(y_{\ast}-\hat{\mu}_{\ast})^{2}|\mathbf{X},\mathbf{y}],$

where the expectation is taken with respect to $(\mathbf{x}_{\ast},y_{\ast})$.
Following the framework for in-sample prediction, we define the expected
optimism as
$\mathrm{OptR}_{\mathbf{X}}=\mathrm{ErrR}_{\mathbf{X}}-\mathrm{ErrT}_{\mathbf{X}}$,
where $\mathrm{ErrR}_{\mathbf{X}}$ and $\mathrm{ErrT}_{\mathbf{X}}$ are the
expectations of $\mathrm{ErrR}_{\mathbf{X},\mathbf{y}}$ and
$\mathrm{ErrT}_{\mathbf{X},\mathbf{y}}$ with respect to $\mathbf{y}$
respectively.

Consider a linear procedure with hat matrix $\mathbf{H}$. For any new
covariate vector $\mathbf{x}_{\ast}$, the procedure predicts the response at
$\mathbf{x}_{\ast}$ using a linear combination of $y_{i}$’s, denoted by
$\hat{\mu}_{\ast}=\mathbf{h}_{\ast}^{\top}\mathbf{y}$, where
$\mathbf{h}_{\ast}\in\mathbb{R}^{n}$ depends only on $\mathbf{x}_{\ast}$ and
$\mathbf{X}$. This coefficient vector $\mathbf{h}_{\ast}$ is called the hat
vector at $\mathbf{x}_{\ast}$ reminiscent of the hat matrix. Then, we can
write the bias-variance decomposition for $\mathrm{ErrR}_{\mathbf{X}}$ and
$\mathrm{ErrT}_{\mathbf{X}}$ as

$\displaystyle\mathrm{ErrR}_{\mathbf{X}}=\sigma_{\varepsilon}^{2}+\mathrm{E}[(\mu_{\ast}-\mathbf{h}_{\ast}^{\top}\bm{\mu})^{2}|\mathbf{X}]+\sigma_{\varepsilon}^{2}\mathrm{E}(\|\mathbf{h}_{\ast}\|^{2}|\mathbf{X}),$
(2.3)
$\displaystyle\mathrm{ErrT}_{\mathbf{X}}=\sigma_{\varepsilon}^{2}+\frac{1}{n}\|\bm{\mu}-\mathbf{H}\bm{\mu}\|^{2}+\frac{1}{n}\sigma_{\varepsilon}^{2}{\rm
tr}(\mathbf{H}^{\top}\mathbf{H}-2\mathbf{H}),$ (2.4)

where the three terms on the right hand side are the irreducible error,
(squared) bias and variance respectively. The expected optimism
$\mathrm{OptR}_{\mathbf{X}}$ can then be expressed as

$\displaystyle\mathrm{OptR}_{\mathbf{X}}$ $\displaystyle=\Delta
B_{\mathbf{X}}+\frac{2}{n}\sigma_{\varepsilon}^{2}\left[\mathrm{df}_{\rm
F}+\frac{n}{2}\left(\mathrm{E}(\|\mathbf{h}_{\ast}\|^{2}|\mathbf{X})-\frac{1}{n}{\rm
tr}(\mathbf{H}\mathbf{H}^{\top})\right)\right],$ (2.5)

where $\Delta
B_{\mathbf{X}}=\mathrm{E}[(\mu_{\ast}-\mathbf{h}_{\ast}^{\top}\bm{\mu})^{2}|\mathbf{X}]-\frac{1}{n}\|\bm{\mu}-\mathbf{H}\bm{\mu}\|^{2}$
is the excess bias and $\mathrm{df}_{\rm F}$ is the classical model degrees of
freedom defined in (2.2). The second term in the right hand side of (2.5) is
the excess variance, and it doesn’t depend on the underlying true model $\mu$.
Further, the quantity in the square brackets depends only on the modeling
procedure through the hat matrix $\mathbf{H}$ and hat vector
$\mathbf{h}_{\ast}$. Since model complexity is an intrinsic property of a
modeling procedure, Luan et al., (2021) defined this quantity as a model
complexity measure under the out-of-sample prediction setting. Denoting it by
$\mathrm{df}_{\rm R}$:

$\mathrm{df}_{\rm R}=\mathrm{df}_{\rm
F}+\frac{n}{2}\left(\mathrm{E}(\|\mathbf{h}_{\ast}\|^{2}|\mathbf{X})-\frac{1}{n}{\rm
tr}(\mathbf{H}^{\top}\mathbf{H})\right),$ (2.6)

they call $\mathrm{df}_{\rm R}$ the predictive model degrees of freedom to
differentiate it from the classical model degrees of freedom and emphasize its
use in an out-of-sample prediction setting.

##  3 Model Complexity in a Heteroscedastic Setting

In this section, we extend the existing theories on model complexity to a
heteroscedastic setting. Our extensions include two aspects. First, we allow
heteroscedasticity in error $\varepsilon$. In particular, we assume that
$\mathrm{Var}(\varepsilon|X)=\sigma_{\varepsilon}^{2}\tau(X)$, where $\tau$ is
a function satisfying $\tau(\cdot)>0$ and $\mathrm{E}[\tau(X)]=1$. We write
$\tau_{i}=\tau(\mathbf{x}_{i})$ for simplicity. Second, for model evaluation,
we consider a general weighted squared error loss. Let
$w\colon\mathbb{R}^{m}\to\mathbb{R}^{+}$ be a deterministic function and write
$w_{i}=w(\mathbf{x}_{i})$. We define the weighted training error given
$(\mathbf{X},\mathbf{y})$ and $w$ as

$\mathrm{wErrT}_{\mathbf{X},\mathbf{y}}=\frac{1}{n}\sum_{i=1}^{n}w_{i}(y_{i}-\hat{\mu}_{i})^{2}=\frac{1}{n}\|\mathbf{y}-\hat{\bm{\mu}}\|_{\mathbf{W}}^{2},$

where $\mathbf{W}=\mathrm{diag}(w_{1},\ldots,w_{n})$ and
$\|\mathbf{a}\|_{\mathbf{W}}^{2}=\mathbf{a}^{\top}\mathbf{W}\mathbf{a}$ for
$\mathbf{a}\in\mathbb{R}^{n}$. In particular, if we assume
$y_{i}\sim\mathcal{N}(\mu_{i},\sigma_{\varepsilon}^{2}\tau_{i})$ given
$\mathbf{x}_{i}$ and take $w_{i}=1/\tau_{i}$, then minimizing the weighted
training error is equivalent to maximizing the log-likelihood of $\mathbf{y}$
given $\mathbf{X}$. For in-sample and out-of-sample prediction errors, we
incorporate case weights in a similar way and define

$\mathrm{wErrF}_{\mathbf{X},\mathbf{y}}=\frac{1}{n}\mathrm{E}(\|\tilde{\mathbf{y}}-\hat{\bm{\mu}}\|_{\mathbf{W}}^{2}|\mathbf{X},\mathbf{y})\text{
and
}\mathrm{wErrR}_{\mathbf{X},\mathbf{y}}=\mathrm{E}[w_{\ast}(y_{\ast}-\hat{\mu}_{\ast})^{2}|\mathbf{X},\mathbf{y}],$

where $\tilde{\mathbf{y}}$ is an independent copy of $\mathbf{y}$ given
$\mathbf{X}$ and $w_{\ast}=w(\mathbf{x}_{\ast})$.

###  3.1 Classical Model Degrees of Freedom

We first look at the classical model degrees of freedom. With the weighted
training error and in-sample prediction error, it is easy to adjust the
expected optimism (2.1) with weights as

$\mathrm{wOptF}_{\mathbf{X}}=\frac{2}{n}\sum_{i=1}^{n}w_{i}\mathrm{Cov}(y_{i},\hat{\mu}_{i}|\mathbf{X}).$

Based on the expected optimism, we define the extended model degrees of
freedom as

$\mathrm{df}_{\rm
F}=\sum_{i=1}^{n}\frac{w_{i}}{\bar{w}}\frac{\mathrm{Cov}(y_{i},\hat{\mu}_{i}|\mathbf{X})}{\sigma_{\varepsilon}^{2}},$
(3.1)

where $\bar{w}=\sum_{i=1}^{n}w_{i}/n$. The use of $w_{i}/\bar{w}$ as a
normalized weight for each covariance summand guarantees that
$\mathrm{df}_{\rm F}$ is invariant to scaling of the weight. In addition, the
normalizing constant $\bar{w}$ ensures that the normalized weights sum up to
the sample size $n$. In particular, when $w_{1}=\cdots=w_{n}$, the extended
definition just reduces to the model degrees of freedom (2.2) for equal error
variances.

Using the extended model degrees of freedom $\mathrm{df}_{\rm F}$, we can
immediately get an unbiased estimator of $\mathrm{wErrF}_{\mathbf{X}}$, which
is given by

$\widehat{\mathrm{wErrF}}=\mathrm{wErrT}_{\mathbf{X},\mathbf{y}}+\frac{2}{n}\bar{w}\sigma_{\varepsilon}^{2}\mathrm{df}_{\rm
F}.$

For a linear procedure with hat matrix $\mathbf{H}$, we have
$\mathrm{Cov}(y_{i},\hat{\mu}_{i}|\mathbf{X})=h_{ii}\sigma_{\varepsilon}^{2}\tau_{i}$.
Then (3.1) is simplified to

$\mathrm{df}_{\rm
F}=\sum_{i=1}^{n}\frac{w_{i}}{\bar{w}}h_{ii}\tau_{i}=\frac{1}{\bar{w}}{\rm
tr}(\mathbf{H}\mathbf{T}\mathbf{W}),$ (3.2)

where $\mathbf{T}=\mathrm{diag}(\tau_{1},\ldots,\tau_{n})$. When
$\mathbf{W}=\mathbf{T}=\mathbf{I}_{n}$, this further reduces to
$\mathrm{df}_{\rm F}={\rm tr}(\mathbf{H})$ in the homoscedastic setting. The
following remarks demonstrate some interesting properties of this extension
(3.2).

######  Remark 3.1 (Optimal weights).

Assume the optimal weights are used for model evaluation, i.e.,
$w(\cdot)\propto 1/\tau(\cdot)$. Then we have

$\mathrm{df}_{\rm F}=\zeta_{\tau}{\rm tr}(\mathbf{H}),$ (3.3)

where $\zeta_{\tau}=\left(\frac{1}{n}\sum_{i=1}^{n}\tau_{i}^{-1}\right)^{-1}$
is the harmonic mean of $\tau_{1},\ldots,\tau_{n}$. Since $\tau_{i}$’s are all
positive, their harmonic mean is always no greater than the arithmetic mean
$\bar{\tau}$, and the equality holds if and only if all $\tau_{i}$’s are
equal. This implies that

$\mathrm{df}_{\rm F}=\zeta_{\tau}{\rm tr}(\mathbf{H})\leq\bar{\tau}{\rm
tr}(\mathbf{H})\approx{\rm tr}(\mathbf{H}),$

where the approximation in the last step is due to
$\bar{\tau}\approx\mathrm{E}[\tau(\mathbf{x})]=1$. Note that ${\rm
tr}(\mathbf{H})$ is the classical model degrees of freedom of a linear
procedure under the homoscedastic setting. This suggests that the unevenness
in error variances leads to a deflation in model complexity.

######  Remark 3.2 (Equal error variances).

When error variances are all equal, i.e., $\tau_{1}=\cdots=\tau_{n}=1$, (3.2)
becomes $\mathrm{df}_{\rm F}=\sum_{i=1}^{n}(w_{i}/\bar{w})h_{ii}$, which is a
weighted average of the leverages $h_{ii}$. Assume $h_{ii}>0$ for
$i=1,\ldots,n$. Then, for any $w_{1},\ldots,w_{n}>0$, we have
$n\min_{i}h_{ii}\leq\mathrm{df}_{\rm F}\leq n\max_{i}h_{ii}$. In other words,
$\mathrm{df}_{\rm F}$ is maximized (minimized) when all weights are assigned
to the case with the largest (smallest) leverage. Since leverages are directly
related to the impact of observations on the model, putting more weights on
more (less) influential observations will make the model more (less) flexible.
The extended model degrees of freedom properly reflects this relationship.

######  Remark 3.3 (Duplicate data).

Assume $\mathbf{x}_{1},\ldots,\mathbf{x}_{k}$ are $k$ distinct covariate
vectors in the training data with $n_{1},\ldots,n_{k}$ replicates respectively
($n=\sum_{i=1}^{k}n_{i}$). For each $i=1,\ldots,k$, let $y_{i,1},\ldots
y_{i,n_{i}}$ be i.i.d. observations given $\mathbf{x}_{i}$ with constant
variance $\sigma_{\varepsilon}^{2}$, and define
$\bar{y}_{i}=\sum_{j=1}^{n_{i}}y_{i,j}/n_{i}$. Note that

$\sum_{i=1}^{k}\sum_{j=1}^{n_{i}}(y_{i,j}-\hat{\mu}_{i})^{2}=\sum_{i=1}^{k}\sum_{j=1}^{n_{i}}(y_{i,j}-\bar{y}_{i})^{2}+\sum_{i=1}^{k}n_{i}(\bar{y}_{i}-\hat{\mu}_{i})^{2},$

and the first term on the right hand side of the above equation doesn’t depend
on the model. Thus, a weighted least squares model based on the collapsed data
$\\{(\mathbf{x}_{i},\bar{y}_{i})|i=1,\ldots,k\\}$ with weights $w_{i}=n_{i}$
is identical to the ordinary least squares model based on the full data
$\\{(\mathbf{x}_{i},y_{i,j})|i=1,\ldots,k;\,j=1,\ldots,n_{i}\\}$. This dual
representation of the model brings a question about its degrees of freedom. If
the collapsed data were to be viewed as a collection of genuine $k$
independent observations with unequal weights, Remark 3.1 appears to suggest
that the weighted model with a linear procedure has fewer degrees of freedom
than the unweighted model fit to the same $k$ observations. However, this
difference in model complexity comes from the scale of weights tied to the
sample size we assume. If the weights are scaled to match the original sample
size $n$ rather than the distinct number of observations $k$, the dual
representation of the model fit to duplicate data leads to the same model
degrees of freedom.

###  3.2 Predictive Model Degrees of Freedom

We now extend the notion of model complexity for the heteroscedastic setting
to the out-of-sample prediction framework for all linear procedures. Let
$w_{\ast}=w(\mathbf{x}_{\ast})$ and $\tau_{\ast}=\tau(\mathbf{x}_{\ast})$, and
define the weighted out-of-sample prediction error as

$\mathrm{wErrR}_{\mathbf{X},\mathbf{y}}=\mathrm{E}[w_{\ast}(y_{\ast}-\hat{\mu}_{\ast})^{2}|\mathbf{X},\mathbf{y}],$

where the expectation is with respect to $(\mathbf{x}_{\ast},y_{\ast})$.
Consider a linear procedure with hat matrix $\mathbf{H}$ and hat vector
$\mathbf{h}_{\ast}$ at $\mathbf{x}_{\ast}$. Similar to (2.3) and (2.4) in the
homoscedastic setting, we have the following bias-variance decompositions for
$\mathrm{wErrR}_{\mathbf{X}}$ and $\mathrm{wErrT}_{\mathbf{X}}$:

$\displaystyle\mathrm{wErrR}_{\mathbf{X}}=\sigma_{\varepsilon}^{2}\mathrm{E}(w_{\ast}\tau_{\ast})+\mathrm{E}[w_{\ast}(\mu_{\ast}-\mathbf{h}_{\ast}^{\top}\bm{\mu})^{2}|\mathbf{X}]+\sigma_{\varepsilon}^{2}\mathrm{E}(w_{\ast}\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2}|\mathbf{X}),$
(3.4)
$\displaystyle\mathrm{wErrT}_{\mathbf{X}}=\frac{1}{n}\sigma_{\varepsilon}^{2}\sum_{i=1}^{n}w_{i}\tau_{i}+\frac{1}{n}\|\bm{\mu}-\mathbf{H}\bm{\mu}\|_{\mathbf{W}}^{2}+\frac{1}{n}\sigma_{\varepsilon}^{2}{\rm
tr}(\mathbf{H}^{\top}\mathbf{W}\mathbf{H}\mathbf{T}-2\mathbf{H}\mathbf{T}\mathbf{W}),$
(3.5)

where $\mathrm{wErrR}_{\mathbf{X}}$ is the expectation of
$\mathrm{wErrR}_{\mathbf{X},\mathbf{y}}$ with respect to $\mathbf{y}$. See
Appendix A.1 for detailed derivation of (3.4) and (3.5). As a consequence, the
expected optimism
$\mathrm{wOptR}_{\mathbf{X}}\coloneqq\mathrm{wErrR}_{\mathbf{X}}-\mathrm{wErrT}_{\mathbf{X}}$
is of the form

$\displaystyle\mathrm{wOptR}_{\mathbf{X}}$
$\displaystyle=\sigma_{\varepsilon}^{2}\left(\mathrm{E}(w_{\ast}\tau_{\ast})-\frac{1}{n}\sum_{i=1}^{n}w_{i}\tau_{i}\right)$
(3.6)
$\displaystyle\qquad+\left(\mathrm{E}[w_{\ast}(\mu_{\ast}-\mathbf{h}_{\ast}^{\top}\bm{\mu})^{2}|\mathbf{X}]-\frac{1}{n}\|\bm{\mu}-\mathbf{H}\bm{\mu}\|_{\mathbf{W}}^{2}\right)$
$\displaystyle\qquad+\frac{2}{n}\bar{w}\sigma_{\varepsilon}^{2}\left[\mathrm{df}_{\rm
F}+\frac{n}{2\bar{w}}\left(\mathrm{E}(w_{\ast}\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2}|\mathbf{X})-\frac{1}{n}{\rm
tr}(\mathbf{W}\mathbf{H}\mathbf{T}\mathbf{H}^{\top})\right)\right],$

where $\mathrm{df}_{\rm F}={\rm tr}(\mathbf{H}\mathbf{T}\mathbf{W})/\bar{w}$
is the model degrees of freedom for a linear procedure defined in (3.2). We
call the first term on the right the excess irreducible error, which equals
$0$ in the homoscedastic setting due to $\tau$ and $w$ both being constant. In
general, this term should be negligible as long as the sample size is large.
The second and third terms are respectively the excess bias and excess
variance adjusted by the weight function $w$. Following Luan et al., (2021),
we still use $\Delta B_{\mathbf{X}}$ to denote the excess bias. We then
propose the extended predictive model degrees of freedom based on the excess
variance. Formally, we define it as

$\mathrm{df}_{\rm R}=\mathrm{df}_{\rm
F}+\frac{n}{2\bar{w}}\left(\mathrm{E}(w_{\ast}\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2}|\mathbf{X})-\frac{1}{n}{\rm
tr}(\mathbf{W}\mathbf{H}\mathbf{T}\mathbf{H}^{\top})\right).$ (3.7)

In particular, when $w(\cdot)$ is constant, we have

$\mathrm{df}_{\rm R}=\mathrm{df}_{\rm
F}+\frac{n}{2}\left(\mathrm{E}(\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2}|\mathbf{X})-\frac{1}{n}{\rm
tr}(\mathbf{H}\mathbf{T}\mathbf{H}^{\top})\right).$

If further $\tau_{i}$’s are all equal, this reduces to (2.6).

The expression in (3.7) suggests that the predictive model degrees of freedom
$\mathrm{df}_{\rm R}$ adjusts the classical model degrees of freedom
$\mathrm{df}_{\rm F}$ with an extra term that accounts for out-of-sample
prediction. In fact, this adjustment can be expressed and interpreted in two
different ways, which we summarize in the following proposition.

######  Proposition 3.1.

For a linear procedure, the adjustment $\mathrm{df}_{\rm R}-\mathrm{df}_{\rm
F}$ can be expressed as

  1. (i)

Covariance penalty representation

$\mathrm{df}_{\rm R}-\mathrm{df}_{\rm
F}=\frac{n}{2\bar{w}}\sum_{i=1}^{n}\tau_{i}\left[\mathrm{E}\left(w_{\ast}\frac{\mathrm{Cov}^{2}(y_{i},\hat{\mu}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})}{(\sigma_{\varepsilon}^{2}\tau_{i})^{2}}\Bigg{|}\mathbf{X}\right)-\frac{1}{n}\sum_{j=1}^{n}w_{j}\frac{\mathrm{Cov}^{2}(y_{i},\hat{\mu}_{j}|\mathbf{X})}{(\sigma_{\varepsilon}^{2}\tau_{i})^{2}}\right];$

  2. (ii)

Generalized degrees of freedom (GDF) representation

$\mathrm{df}_{\rm R}-\mathrm{df}_{\rm
F}=\frac{n}{2\bar{w}}\sum_{i=1}^{n}\tau_{i}\left[\mathrm{E}\left(w_{\ast}\left(\frac{\partial\mathrm{E}(\hat{\mu}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})}{\partial\mu_{i}}\right)^{2}\Bigg{|}\mathbf{X}\right)-\sum_{j=1}^{n}w_{j}\left(\frac{\partial\mathrm{E}(\hat{\mu}_{j}|\mathbf{X})}{\partial\mu_{i}}\right)^{2}\right].$

The covariance penalty representation unifies the two measures of model
degrees of freedom $\mathrm{df}_{\rm F}$ and $\mathrm{df}_{\rm R}$ by
quantifying the level of dependency of the predicted values on the observed
values. More specifically, the adjustment $\mathrm{df}_{\rm
R}-\mathrm{df}_{\rm F}$ compares the overall level of such dependency for out-
of-sample predictions against that for in-sample ones. The GDF representation
gets its name from the generalized degrees of freedom (GDF) proposed by Ye,
(1998). In his work, he defined
$\sum_{i=1}^{n}\frac{\partial\mathrm{E}(\hat{\mu}_{i}|\mathbf{X})}{\partial\mu_{i}}$
as the GDF for a modeling procedure $\hat{\bm{\mu}}$. The GDF representation
of $\mathrm{df}_{\rm R}-\mathrm{df}_{\rm F}$ expresses this adjustment in
terms of the rate of change of the predicted values on the test data against
that on the training data with respect to the true mean components
$\mu_{i}$’s.

##  4 Weighted Least Squares

We study the two proposed measures of model degrees of freedom for the
weighted least squares method in this section. Assume
$\mathbf{X}\in\mathbb{R}^{n\times m}$ has full column rank. Let
$q\colon\mathbb{R}^{m}\to\mathbb{R}_{+}$ be a known weight function, and
define $q_{i}=q(\mathbf{x}_{i})$ and
$\mathbf{Q}=\mathrm{diag}(q_{1},\ldots,q_{n})$. For the weighted least squares
problem of
$\min_{\mathbf{b}\in\mathbb{R}^{p}}\|\mathbf{y}-\mathbf{X}\mathbf{b}\|_{\mathbf{Q}}^{2}$,
the solution has the explicit form of
$\hat{\bm{\beta}}=(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}\mathbf{y}$.
The corresponding hat matrix and hat vector are then given by
$\mathbf{H}=\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}$
and
$\mathbf{h}_{\ast}=\mathbf{Q}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{x}_{\ast}$.

Consider evaluating the model using a weighted squared error loss with weight
function $w$. Let
$\bm{\Sigma}=\mathrm{E}(w_{\ast}\mathbf{x}_{\ast}\mathbf{x}_{\ast}^{\top})/\bar{w}$.
Then it is easy to show that

$\displaystyle\mathrm{df}_{\rm F}$ $\displaystyle=\frac{1}{\bar{w}}{\rm
tr}[(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}\mathbf{T}\mathbf{W}\mathbf{X}],$
(4.1) $\displaystyle\mathrm{df}_{\rm R}$ $\displaystyle=\mathrm{df}_{\rm
F}+\frac{n}{2}{\rm
tr}\left[(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}\mathbf{T}\mathbf{Q}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\left(\bm{\Sigma}-\frac{1}{n\bar{w}}\mathbf{X}^{\top}\mathbf{W}\mathbf{X}\right)\right].$
(4.2)

###  4.1 Subset Regression

We first study $\mathrm{df}_{\rm F}$ and $\mathrm{df}_{\rm R}$ of the weighted
least squares model in the subset regression setting. For simplicity, we
assume the weight scheme used for modeling fitting is the same as that for
model evaluation, i.e., $q=w$. Let
$\mathcal{S}_{1}\subset\cdots\subset\mathcal{S}_{d}\coloneqq\\{1,\ldots,d\\}$
be a sequence of nested subsets of variables such that $|\mathcal{S}_{p}|=p$
for each $p=1,\ldots,d$. Let $\mathbf{X}_{p}$ and $\bm{\Sigma}_{p}$ be the
submatrices of $\mathbf{X}$ and $\bm{\Sigma}$ corresponding to
$\mathcal{S}_{p}$ respectively. Consider regressing $\mathbf{y}$ on
$\mathbf{X}_{p}$ using the weighted least squares method with weight function
$w$. Then (4.1) and (4.2) are simplified to

$\displaystyle\mathrm{df}_{\rm F}(p)$ $\displaystyle=\frac{1}{\bar{w}}{\rm
tr}[(\mathbf{X}_{p}^{\top}\mathbf{W}\mathbf{X}_{p})^{-1}\mathbf{X}_{p}^{\top}\mathbf{W}\mathbf{T}\mathbf{W}\mathbf{X}_{p}],$
$\displaystyle\mathrm{df}_{\rm R}(p)$
$\displaystyle=\frac{1}{2}\mathrm{df}_{\rm F}(p)+\frac{n}{2}{\rm
tr}[(\mathbf{X}_{p}^{\top}\mathbf{W}\mathbf{X}_{p})^{-1}\mathbf{X}_{p}^{\top}\mathbf{W}\mathbf{T}\mathbf{W}\mathbf{X}_{p}(\mathbf{X}_{p}^{\top}\mathbf{W}\mathbf{X}_{p})^{-1}\bm{\Sigma}_{p}].$

For each $j=1,\ldots,d$, let $\mathbf{x}_{(j)}$ be the column(s) in
$\mathbf{X}$ associated with variable $j$. If the $j$th variable is continuous
or binary, $\mathbf{x}_{(j)}\in\mathbb{R}^{n}$. If the variable is categorical
with $k\geq 3$ levels, $\mathbf{x}_{(j)}\in\mathbb{R}^{n\times(k-1)}$. For the
sake of convenience and conciseness, we assume
$\mathbf{x}_{(j)}\in\mathbb{R}^{n}$ in the rest of this section, though all
results can be tailored to categorical variables with more than two levels
accordingly.

We examine the size of change in the model degrees of freedom as a covariate
is added. Normalizing observations with weights, let
$\tilde{\mathbf{x}}_{(p+1)}=\mathbf{W}^{\frac{1}{2}}\mathbf{x}_{(p+1)}/\sqrt{\bar{w}}$
and
$\tilde{\mathbf{X}}_{p}=\mathbf{W}^{\frac{1}{2}}\mathbf{X}_{p}/\sqrt{\bar{w}}$.
Similarly, let
$\tilde{\mathbf{y}}=\mathbf{W}^{\frac{1}{2}}\mathbf{y}/\sqrt{\bar{w}}$ and
note that
$\mathrm{Var}(\tilde{\mathbf{y}}|\mathbf{X})=\sigma_{\varepsilon}^{2}\mathbf{W}^{\frac{1}{2}}\mathbf{T}\mathbf{W}^{\frac{1}{2}}$.
Define
$\tilde{\mathbf{H}}=\tilde{\mathbf{X}}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}$
and
$\tilde{\mathbf{r}}_{p+1}=(\mathbf{I}_{n}-\tilde{\mathbf{H}})\tilde{\mathbf{x}}_{(p+1)}$.
Here, $\tilde{\mathbf{r}}_{(p+1)}$ is the residual vector when the normalized
new variable vector $\tilde{\mathbf{x}}_{(p+1)}$ is regressed on the
normalized design matrix $\tilde{\mathbf{X}}_{p}$. The following theorem
describes the increment in $\mathrm{df}_{\rm F}$ in terms of
$\tilde{\mathbf{r}}_{(p+1)}$ as $p$ increases by 1.

######  Theorem 4.1.

Assume $p<d<n$. Then for the weighted least squares model with weight matrix
$\mathbf{W}$ and design matrix $\mathbf{X}_{p}$, $\mathrm{df}_{\rm F}$
increases in $p$, and

$\Delta\mathrm{df}_{\rm F}(p)\coloneqq\mathrm{df}_{\rm
F}(p+1)-\mathrm{df}_{\rm
F}(p)=\frac{\|\tilde{\mathbf{r}}_{(p+1)}\|_{\mathbf{U}}^{2}}{\|\tilde{\mathbf{r}}_{(p+1)}\|^{2}},$

where
$\mathbf{U}=\mathbf{W}^{\frac{1}{2}}\mathbf{T}\mathbf{W}^{\frac{1}{2}}/\bar{w}$.

Theorem 4.1 indicates that the increment in $\mathrm{df}_{\rm F}$ generally
depends on the amount of collinearity between the normalized new variable and
existing variables. A special case occurs when $w\propto 1/\tau$. In this
scenario, we have $\mathrm{df}_{\rm F}=\zeta_{\tau}p$ and
$\Delta\mathrm{df}_{\rm F}(p)=\zeta_{\tau}$ (the harmonic mean of
$\tau_{1},\ldots,\tau_{n}$). Hence, $\mathrm{df}_{\rm F}$ is a linear function
of $p$. Its increment $\zeta_{\tau}$ is constant for all $p$ and depends only
on the unevenness of the true error variances.

For the predictive model degrees of freedom $\mathrm{df}_{\rm R}$, assume that
$\bm{\Sigma}$ is positive definite. For $p=1,\ldots,d-1$, write

$\bm{\Sigma}_{p+1}=\begin{pmatrix}\bm{\Sigma}_{p}&\bm{\phi}_{p,p+1}\\\
\bm{\phi}_{p,p+1}^{\top}&\sigma_{p+1}^{2}\end{pmatrix},\quad p=1,\ldots,d-1.$

Let
$a_{p,p+1}=\sigma_{p+1}^{2}-\bm{\phi}_{p,p+1}^{\top}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}$.
Note that $a_{p,p+1}>0$ when $\bm{\Sigma}$ is positive definite. We derive the
increment $\Delta\mathrm{df}_{\rm R}(p)=\mathrm{df}_{\rm
R}(p+1)-\mathrm{df}_{\rm R}(p)$ in the following theorem.

######  Theorem 4.2.

Assume $p<d<n$. Then for the weighted least squares model with weight matrix
$\mathbf{W}$ and design matrix $\mathbf{X}_{p}$,

$\displaystyle\Delta\mathrm{df}_{\rm R}(p)=\frac{1}{2}\Delta\mathrm{df}_{\rm
F}(p)$ $\displaystyle+\frac{n}{2}\Delta\mathrm{df}_{\rm
F}(p)\frac{(\|\tilde{\mathbf{x}}_{(p+1)}-\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\|_{\mathbf{M}}^{2}+a_{p,p+1})}{\|\tilde{\mathbf{r}}_{(p+1)}\|^{2}}$
(4.3)
$\displaystyle-n\frac{\langle\tilde{\mathbf{r}}_{(p+1)},\tilde{\mathbf{x}}_{(p+1)}-\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\rangle_{\mathbf{U}\mathbf{M}}}{\|\tilde{\mathbf{r}}_{(p+1)}\|^{2}}$

where
$\mathbf{M}=\tilde{\mathbf{X}}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\bm{\Sigma}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}$,
and $\langle\cdot,\cdot\rangle_{\mathbf{A}}$ for a matrix $\mathbf{A}$ denotes
the inner product with respect to $\mathbf{A}$.

The component
$\tilde{\mathbf{x}}_{(p+1)}-\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}$
can be seen as the residual vector $\tilde{\mathbf{r}}_{(p+1)}$ at the
population level, in which we simply replace
$\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p}$ and
$\tilde{\mathbf{X}}^{\top}\tilde{\mathbf{x}}_{(p+1)}$ with their population
versions $n\bar{w}\Sigma_{p}$ and $n\bar{w}\phi_{p,p+1}$ respectively. In
general, $\Delta\mathrm{df}_{\rm R}(p)$ is not necessarily positive due to the
third term in the right hand side of (4.3). However, when $w\propto 1/\tau$,
this term vanishes due to the fact that
$(\mathbf{I}_{n}-\tilde{\mathbf{H}})\mathbf{U}\mathbf{M}=\mathbf{O}_{n\times
n}$. In this case, we have

$\mathrm{df}_{\rm R}(p)=\frac{\zeta_{\tau}}{2}(p+n\bar{w}{\rm
tr}[(\mathbf{X}_{p}^{\top}\mathbf{W}\mathbf{X}_{p})^{-1}\bm{\Sigma}_{p})]),$
(4.4)

which is strictly increasing in $p$. The following remark gives an
approximation of (4.4) under the assumption that $\mathbf{x}_{i}$’s are
multivariate normal.

######  Remark 4.1.

Assume $\mathbf{X}\in\mathbb{R}^{n\times p}$. We slightly modify the setting
by assuming $\tau$ is a function of a random variable $\mathbf{u}$ that is
independent of $\mathbf{x}$. Further, assume that $\mathbf{x}_{\ast}$ is
multivariate normal with $\mathrm{E}(\mathbf{x}_{\ast})=\mathbf{0}$. Then it
can be shown that

$\mathrm{df}_{\rm
R}(p)\approx\frac{p}{2}\zeta_{\tau}\left(1+\frac{n_{\tau}}{n_{\tau}-p-1}\right),$
(4.5)

where $n_{\tau}$ is the nearest integer to
$\left(\sum_{i=1}^{n}1/\tau_{i}\right)^{2}/\sum_{i=1}^{n}(1/\tau_{i})^{2}$.
For the derivation of (4.5), see Appendix A.5.

Figure 1 shows the approximation (4.5) under the assumption that
$\tau(u_{i})=(1+|u_{i}|)^{\lambda}$ ($\lambda\geq 0$) and
$\mu(\mathbf{x}_{i})=\sqrt{10/d}\sum_{j=1}^{d}x_{ij}$, where $x_{ij}$’s and
$u_{i}$ are i.i.d. $\mathcal{N}(0,1)$. We set $n=100$ and $d=50$ and consider
$\lambda=0.5,1,2\text{ and }3$. The distribution of $\tau(u)$ becomes more
skewed to the right as $\lambda$ increases. Consequently, the rate of increase
in both $\mathrm{df}_{\rm F}$ and $\mathrm{df}_{\rm R}$ drops as the skewness
level $\lambda$ increases. The approximation (4.5) works well when the
distribution of $\tau$ is less skewed, but tends to overestimate the true
predictive model degrees of freedom with more uneven $\tau_{i}$’s at a large
$p$ due to the quickly shrinking denominator $n_{\tau}-p-1$.

Figure 1: The approximated predictive model degrees of freedom
$\mathrm{df}_{\rm R}$ of the weighted least squares model with $w=1/\tau$ as
the weight function. The results are obtained with
$\mu(\mathbf{x}_{i})=\sqrt{10/d}\sum_{j=1}^{d}x_{ij}$ and
$\tau(u_{i})=(1+|u_{i}|)^{\lambda}$ ($\lambda\geq 0$), where $d=50$, $n=100$,
and $x_{ij}$’s and $u_{i}$ are i.i.d. $\mathcal{N}(0,1)$.

In fact, the size of $\zeta_{\tau}$ and $n_{\tau}$ are both directly related
to the unevenness of $\tau_{i}$’s. When $w=1/\tau=1$, we have $\zeta_{\tau}=1$
and $n_{\tau}=n$. In this case, (4.5) becomes

$\mathrm{df}_{\rm R}(p)\approx\frac{p}{2}\left(1+\frac{n}{n-p-1}\right),$
(4.6)

which is the approximation of $\mathrm{df}_{\rm R}$ for subset regression
given by Luan et al., (2021) under the homoscedastic setting. On the other
hand, when $\tau_{i}$’s are not all equal, $\zeta_{\tau}<1$ and $n_{\tau}\leq
n$. In particular, Figure 2 suggests that the more uneven $\tau_{i}$’s are,
the smaller $\zeta_{\tau}$ and $n_{\tau}$ will be. Since $n_{\tau}$ in (4.5)
plays a similar role as $n$ in (4.6), we call it the effective sample size.
When $\tau_{i}$’s are highly uneven, observations with large variance will
receive little weight. The resulting effect is as if the sample size was
reduced.

Figure 2: The coefficient $\zeta_{\tau}$ and effective sample size $n_{\tau}$
as a function of the skewness parameter $\lambda$ when $n=100$, $d=50$, and
$\tau(u)=(1+|u|)^{\lambda}$ ($\lambda\geq 0$) with $u\sim\mathcal{N}(0,1)$ and
independent of the covariates $\mathbf{x}$.

###  4.2 Effect of Weight Schemes

Since the optimal weight scheme is usually unknown a priori, misspecification
of weights could happen in practice. A natural question is then how the model
complexity is affected by the choice of weights. In this section, we try to
answer this question by fixing the subset of variables for regression and
varying the weight scheme $q$ for the weighted least squares model. Since
these models are trained under different weight schemes, it is important to
assess their complexity using a common criterion. Here, we set $w=1/\tau$ and
define model degrees of freedoms for all weighted least squares models using
this weight.

Assume that $\mathbf{X}\in\mathbb{R}^{n\times p}$ is the design matrix. To
avoid confusion, we use $\mathrm{wErrT}_{\mathbf{X}}(q)$,
$\mathrm{wErrF}_{\mathbf{X}}(q)$ and $\mathrm{wErrR}_{\mathbf{X}}(q)$ in this
section to denote the training and prediction errors of the weighted least
squares model $\hat{\bm{\mu}}(q)$ with weight scheme $q$. Similarly, we use
$\mathrm{df}_{\rm F}(q)$ and $\mathrm{df}_{\rm R}(q)$ here to denote the
corresponding measures of model degrees of freedom. By replacing $\mathbf{W}$
with $\mathbf{T}^{-1}$ in (4.1) and (4.2), we get $\mathrm{df}_{\rm
F}(q)=\zeta_{\tau}p$ and

$\mathrm{df}_{\rm R}(q)=\zeta_{\tau}p+\frac{n}{2}{\rm
tr}[(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}\mathbf{T}\mathbf{Q}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{R}],$
(4.7)

where $\mathbf{R}=\bm{\Sigma}-\mathbf{X}^{\top}\mathbf{T}^{-1}\mathbf{X}/n$.

For the classical model degrees of freedom $\mathrm{df}_{\rm F}$, we find that
$\mathrm{df}_{\rm F}$ is constant with respect to the weight scheme $q$. In
other words, the average covariance between the fitted values $\hat{\mu}_{i}$
and the observations $y_{i}$ is the same for all weighted least squares
models. Since
$\mathrm{wErrF}_{\mathbf{X}}(q)=\mathrm{wErrT}_{\mathbf{X}}(q)+(2/n)\bar{w}\sigma_{\varepsilon}^{2}\mathrm{df}_{\rm
F}(q)$, the best weighted least squares model should be the one that minimizes
the expected training error $\mathrm{wErrT}_{\mathbf{X}}(q)$. Apparently, this
minimizer is the model with the optimal weight scheme $q=1/\tau$, since it is
specifically determined to minimize the training error
$\mathrm{wErrT}_{\mathbf{X},\mathbf{y}}(q)=\|\mathbf{y}-\hat{\bm{\mu}}(q)\|_{\mathbf{T}^{-1}}^{2}/n$
for each $\mathbf{y}$.

For the predictive model degrees of freedom (4.7), the following theorem gives
its gradient with respect to $\mathbf{q}=(q_{1},\ldots,q_{n})^{\top}$.

######  Theorem 4.3.

The gradient of the predictive model degrees of freedom $\mathrm{df}_{\rm R}$
with respect to the weight vector $\mathbf{q}=\mathrm{diag}(\mathbf{Q})$ is
given by

$\frac{\partial\mathrm{df}_{\rm
R}}{\partial\mathbf{q}}=n\,\mathrm{diag}\left((\mathbf{I}_{n}-\mathbf{H})\mathbf{T}\mathbf{Q}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{R}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\right),$

where
$\mathbf{H}=\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}$.

For a general weight scheme $q$, it is usually not easy to see directly from
$\partial\mathrm{df}_{\rm R}/\partial\mathbf{q}$ how $\mathrm{df}_{\rm R}$
behaves when $q$ is perturbed. However, when $q=1/\tau$, we have
$(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}\mathbf{Q}\mathbf{X}=(\mathbf{I}_{n}-\mathbf{H})\mathbf{X}=\mathbf{O}_{n\times
p}$, which further implies that $(\partial\mathrm{df}_{\rm
R}/\partial\mathbf{q})|_{q=1/\tau}=\mathbf{0}_{n}$. This indicates that
$\mathrm{df}_{\rm R}$, as a function of the weight scheme $q$, is relatively
flat around the optimal weight $1/\tau$. In practice, we usually need to
estimate $\tau$ with an estimator $\hat{\tau}$ and use $\hat{q}=1/\hat{\tau}$
as a weight function for model fitting. A good estimate of $\tau$ guarantees
that $\hat{q}$ is close to the optimal weight $1/\tau$. Therefore, it is
desirable for $\mathrm{df}_{\rm R}$ to be flat around the optimal weight in
the sense that the discrepancy between the optimal weight scheme $1/\tau$ and
a reasonably good estimate of it would make little change in $\mathrm{df}_{\rm
R}$.

A more quantitative study of $\mathrm{df}_{\rm R}$ is given as follows. Assume
that $q$ is of the form

$q(\alpha;\mathbf{x})=\frac{1}{\alpha+(1-\alpha)\tau(\mathbf{x})},\text{ for
}0\leq\alpha\leq 1.$ (4.8)

Then $\mathrm{E}[1/{q(\alpha;\mathbf{x})}]=1=\mathrm{E}[\tau(\mathbf{x})]$ for
any $\alpha$. Note that $q=1/\tau$ when $\alpha=0$ and $q=1$ when $\alpha=1$.
Thus, (4.8) constructs a path from the optimal to the equal weight scheme. To
examine $\mathrm{df}_{\rm R}$ as a function of $\alpha$ more closely, we
consider a simulation study below. We set $n=100$ and $d=20$, and assume that
$\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{d})$.
Also assume that
$\tau(\mathbf{x})\propto(1+|\mathbf{x}^{\top}\bm{\beta}|)^{2}$ with
$\beta_{j}\propto(1-j/d)^{5}$ and $\|\bm{\beta}\|^{2}=10$. We generate 100
replicates of $\mathbf{X}$. Figure 3 demonstrates the predictive model degrees
of freedom of models with the first $p$ variables for $p=1,\,5,\,10$ and $20$.
We can see that $\mathrm{df}_{\rm R}$ can be either increasing or decreasing
in $\alpha$. However, as the number of variables increases, $\mathrm{df}_{\rm
R}$ is more likely to be an increasing function of $\alpha$. Note that, as
more relevant variables are added, the weighted prediction errors
$(y_{\ast}-\hat{\mu}_{\ast})/\sqrt{\tau_{\ast}}$ of the model with the optimal
weight scheme are generally homoscedastic, while those of the other models are
not. Consequently, the variance of the predicted values produced by the model
with optimal weights is the lowest on average. Since the predictive model
degrees of freedom is defined through the excess variance, the optimal weight
scheme generally leads to the simplest model.

Figure 3: The predictive model degrees of freedom $\mathrm{df}_{\rm R}$ of the
weighted least squares method as a function of the coefficient $\alpha$ when
$\mathbf{Q}=(\alpha\mathbf{I}_{n}+(1-\alpha)\mathbf{T})^{-1}$
($0\leq\alpha\leq 1$). We assume that
$\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{d})$,
$\tau(\mathbf{x})\propto(1+|\mathbf{x}^{\top}\bm{\beta}|)^{2}$,
$\beta_{j}\propto(1-j/d)^{5}$ and $\|\bm{\beta}\|^{2}=10$ with $n=100$ and
$d=20$. The gray curves represent $\mathrm{df}_{\rm R}$ based on 100 different
realizations of $\mathbf{X}$, whereas the red curve is the average over these
replicates.

On the other hand, though the weight scheme $q$ for model fitting affects the
predictive model degrees of freedom, our results in Figure 3 suggest that its
impact may not be big enough to produce appreciable range of difference. When
the model size is not too large, which is usually the case in many real
applications, the difference in $\mathrm{df}_{\rm R}$ due to the weight scheme
$q$ for model fitting is fairly small or even trivial, as we can see from the
figure. In addition, the expected optimism (3.6) implies that

$\mathrm{wErrR}_{\mathbf{X}}=\mathrm{wErrT}_{\mathbf{X}}+\sigma_{\varepsilon}^{2}\left(\mathrm{E}(w_{\ast}\tau_{\ast})-\frac{1}{n}\sum_{i=1}^{n}w_{i}\tau_{i}\right)+\Delta
B_{\mathbf{X}}+\frac{2}{n}\bar{w}\sigma_{\varepsilon}^{2}\mathrm{df}_{\rm R}.$
(4.9)

Since $\mathrm{df}_{\rm R}$ contributes to the prediction error
$\mathrm{wErrR}_{\mathbf{X}}$ through
$(2/n)\bar{w}\sigma_{\varepsilon}^{2}\mathrm{df}_{\rm R}$, the impact of the
difference in $\mathrm{df}_{\rm R}$ for different weighted least squares
models will be reduced when the sample size $n$ is large. By contrast, the
choice of $q$ may have a much larger impact on the training error
$\mathrm{wErrT}_{\mathbf{X}}$ or excess bias $\Delta B_{\mathbf{X}}$. In
Section 5, we will take these terms into account and provide a risk estimator
for model evaluation and comparison.

##  5 Prediction Error Estimation

In this section, we aim to construct an estimator for the out-of-sample
prediction error based on (4.9). For the excess irreducible error, the law of
large numbers guarantees that
$\mathrm{E}(w_{\ast}\tau_{\ast})-\sum_{i=1}^{n}w_{i}\tau_{i}/n$ is generally
negligible when $n$ is large. The excess bias $\Delta B_{\mathbf{X}}$ depends
on the unknown true mean function $\mu(\cdot;\bm{\beta})$, and thus needs to
be estimated. The predictive model degrees of freedom depends on the
distribution of $\mathbf{x}_{\ast}$ through
$\mathrm{E}(w_{\ast}\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2}|\mathbf{X})$. In
practice, this expectation also needs to be estimated.

For simplicity, we assume that $\tau$ is either known a priori or can be
estimated by $\hat{\tau}$ using some independent data (e.g., historical data)
throughout this section. Further, we assume that $q=w\propto 1/\tau$ (or
$1/\hat{\tau}$).

###  5.1 Estimation of $\Delta B_{\mathbf{X}}$

We first consider estimating the excess bias $\Delta B_{\mathbf{X}}$. Here, we
adjust with weights the method proposed in Luan et al., (2021) that uses the
leave-one-out cross validation (LOOCV) trick. Let $\mathbf{X}^{-i}$,
$\mathbf{y}^{-i}$ and $\bm{\mu}^{-i}$ be the corresponding terms with the
$i$th case deleted. Write
$\mathbf{H}=(\mathbf{h}_{1},\ldots,\mathbf{h}_{n})^{\top}$ and define
$\mathbf{h}^{-i}_{i}$ to be the hat vector at $\mathbf{x}_{i}$ based on
$\mathbf{X}^{-i}$. Conceptually, when $n$ is large, we have

$\mathrm{E}[w_{\ast}(\mu_{\ast}-\mathbf{h}_{\ast}^{\top}\bm{\mu})^{2}|\mathbf{X}]\approx\frac{1}{n}\sum_{i=1}^{n}w_{i}(\mu_{i}-(\mathbf{h}^{-i}_{i})^{\top}\bm{\mu}^{-i})^{2}.$

For linear procedures, the right hand side of the above equation can be
evaluated efficiently due to the LOOCV identity
$y_{i}-(\mathbf{h}_{i}^{-i})^{\top}\mathbf{y}^{-i}=(y_{i}-\mathbf{h}_{i}^{\top}\mathbf{y})/(1-h_{ii})$
(Craven and Wahba,, 1978). Since $\mathbf{h}_{i}$ doesn’t depend on either
$\bm{\mu}$ or $\mathbf{y}$, the identity still holds if we replace $y_{i}$
with $\mu_{i}$ and $\mathbf{y}^{-i}$ with $\bm{\mu}^{-i}$ respectively. Then,
we have

$\displaystyle\Delta B_{\mathbf{X}}$
$\displaystyle\approx\frac{1}{n}\sum_{i=1}^{n}w_{i}[(\mu_{i}-(\mathbf{h}^{-i}_{i})^{\top}\bm{\mu}^{-i})^{2}-(\mu_{i}-\mathbf{h}_{i}^{\top}\bm{\mu})^{2}]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}w_{i}\left[\frac{(\mu_{i}-\mathbf{h}_{i}^{\top}\bm{\mu})^{2}}{(1-h_{ii})^{2}}-(\mu_{i}-\mathbf{h}_{i}^{\top}\bm{\mu})^{2}\right]$
$\displaystyle=\frac{1}{n}\bm{\mu}^{\top}\mathbf{A}\bm{\mu}$
$\displaystyle=\frac{1}{n}\mathbf{y}^{\top}\mathbf{A}\mathbf{y}-\frac{2}{n}\bm{\mu}^{\top}\mathbf{A}\bm{\varepsilon}-\frac{1}{n}\bm{\varepsilon}^{\top}\mathbf{A}\bm{\varepsilon},$

where
$\mathbf{A}=(\mathbf{I}_{n}-\mathbf{H})^{\top}\mathbf{D}(\mathbf{I}_{n}-\mathbf{H})$
and $\mathbf{D}=\mathrm{diag}\left(w_{i}/(1-h_{ii})^{2}-w_{i}\right)$. Since
$\mathrm{E}(\bm{\mu}^{\top}\mathbf{A}\bm{\varepsilon}|\mathbf{X})=0$ and
$\mathrm{E}(\bm{\varepsilon}^{\top}\mathbf{A}\bm{\varepsilon}|\mathbf{X})=\sigma_{\varepsilon}^{2}{\rm
tr}(\mathbf{A}\mathbf{T})$, we can estimate $\Delta B_{\mathbf{X}}$ by

$\hat{\delta}=\frac{1}{n}\mathbf{y}^{\top}\mathbf{A}\mathbf{y}-\frac{1}{n}\sigma_{\varepsilon}^{2}{\rm
tr}(\mathbf{A}\mathbf{T}).$

Assume $\mathrm{df}_{\rm R}$ is known. Then (4.9) implies that an estimator of
$\mathrm{wErrR}_{\mathbf{X}}$ is given by

$\widehat{\mathrm{wErrR}}=\mathrm{wErrT}_{\mathbf{X},\mathbf{y}}+\hat{\delta}+\frac{2}{n}\bar{w}\sigma_{\varepsilon}^{2}\mathrm{df}_{\rm
R}=\widehat{\mathrm{wErrR}}_{\rm
loocv}+\frac{1}{n}\sigma_{\varepsilon}^{2}\xi_{\mathbf{X}},$ (5.1)

where $\widehat{\mathrm{wErrR}}_{\rm
loocv}=\sum_{i=1}^{n}w_{i}(y_{i}-(\mathbf{h}_{i}^{-i})^{\top}\mathbf{y}^{-i})^{2}/n$
is the LOOCV error and $\xi_{\mathbf{X}}=2\bar{w}\mathrm{df}_{\rm R}-{\rm
tr}(\mathbf{A}\mathbf{T})$. Thus, $\widehat{\mathrm{wErrR}}$ is an adjusted
version of the LOOCV error. In fact, this adjustment has a close connection
with the variance of prediction errors in the full model and leave-one-out
models, which we summarize in the following theorem.

######  Theorem 5.1.

The adjustment term $\frac{1}{n}\sigma_{\varepsilon}^{2}\xi_{\mathbf{X}}$ in
(5.1) can be expressed as

$\frac{1}{n}\sigma_{\varepsilon}^{2}\xi_{\mathbf{X}}=\mathrm{E}[w_{\ast}\mathrm{Var}(\hat{\varepsilon}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})|\mathbf{X}]-\frac{1}{n}\sum_{i=1}^{n}w_{i}\mathrm{Var}(\hat{\varepsilon}^{-i}_{i}|\mathbf{X})-\sigma_{\varepsilon}^{2}\left(\mathrm{E}(w_{\ast}\tau_{\ast})-\frac{1}{n}\sum_{i=1}^{n}w_{i}\tau_{i}\right),$

where $\hat{\varepsilon}_{\ast}$ is the prediction error of the full model at
$\mathbf{x}_{\ast}$, and $\hat{\varepsilon}_{i}^{-i}$ is that of the leave-
the-$i$th-out model at $\mathbf{x}_{i}$.

Note that the third term on the right hand side of the expression above is the
excess irreducible error, which is approximately 0 when $n$ is large.
Therefore, $\frac{1}{n}\sigma_{\varepsilon}^{2}\xi_{\mathbf{x}}$ can be
approximated by

$\frac{1}{n}\sigma_{\varepsilon}^{2}\xi_{\mathbf{X}}\approx\mathrm{E}[w_{\ast}\mathrm{Var}(\hat{\varepsilon}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})|\mathbf{X}]-\frac{1}{n}\sum_{i=1}^{n}w_{i}\mathrm{Var}(\hat{\varepsilon}^{-i}_{i}|\mathbf{X}).$

This indicates that the adjustment term
$\frac{1}{n}\sigma_{\varepsilon}^{2}\xi_{\mathbf{X}}$ is a measure of the
difference between the full model and leave-one-out models in the variance of
prediction errors weighted by the case weights $w_{1},\ldots,w_{n}$. In
addition, as shown in Luan et al., (2021), the adjustment term is generally
negligible when $n$ is large and $p/n$ is relatively small. But when $p$ is
close to $n$, ${\rm tr}(\mathbf{A}\mathbf{T})$ could be very large and
variable due to the leverages $h_{ii}$ approaching 1. This may make the
adjustment negative and large in magnitude, and further cause
$\widehat{\mathrm{wErrR}}$ to be negative. Since the term ${\rm
tr}(\mathbf{A}\mathbf{T})$ in $\xi_{\mathbf{X}}$ comes from the excess bias
estimator $\hat{\delta}$, to circumvent this issue, they suggested to use
$\hat{\delta}_{+}=\max(\hat{\delta},0)$ in place of $\hat{\delta}$ in (5.1).
We will follow this adjustment in the rest of the paper.

###  5.2 Estimation of $\mathrm{df}_{\rm R}$

In practice, the distribution of $\mathbf{x}_{\ast}$ is usually unknown, which
makes $\mathrm{E}(w_{\ast}\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2}|\mathbf{X})$
and $\mathrm{df}_{\rm R}$ not directly available. To estimate
$\mathrm{E}(w_{\ast}\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2}|\mathbf{X})$, one
can again use the LOOCV trick and approximate the expectation by
$\sum_{i=1}^{n}w_{i}\|\mathbf{h}_{i}^{-i}\|_{\mathbf{T}}^{2}/n$.
Unfortunately, unlike estimating the excess bias $\Delta B_{\mathbf{X}}$, this
approach inevitably requires fitting $n$ leave-one-out models, which could be
highly inefficient when $n$ is large.

Another approach is to estimate the expectation
$\mathrm{E}(w_{\ast}\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2}|\mathbf{X})$ using
some synthetic data generated based on the training data $\mathbf{X}$ that
might mimic genuine random samples from the joint distribution of
$\mathbf{x}$. Assume that $B$ synthetic examples are generated, denoted by
$\mathbf{x}_{\ast}^{(1)},\ldots,\mathbf{x}_{\ast}^{(B)}$. Then we can estimate
the predictive model degrees of freedom by

$\widehat{\mathrm{df}}_{\mathrm{R}}=\mathrm{df}_{\rm
F}+\frac{n}{2\bar{w}}\left(\frac{1}{B}\sum_{b=1}^{B}w_{\ast}^{(b)}\|\mathbf{h}_{\ast}^{(b)}\|_{\mathbf{T}}^{2}-\frac{1}{n}{\rm
tr}(\mathbf{W}\mathbf{H}\mathbf{T}\mathbf{H}^{\top})\right),$

where $w_{\ast}^{(b)}=w(\mathbf{x}_{\ast}^{(b)})$ and
$\mathbf{h}_{\ast}^{(b)}$ is the hat vector at $\mathbf{x}_{\ast}^{(b)}$
($b=1,\ldots,B$).

In this work, we consider two different procedures to generate synthetic data:
the Classification and Regression Trees (CART) data synthesis method (Reiter,,
2005) and the Naive Bayes Estimation (NBE) learning algorithm (Lowd and
Domingos,, 2005). One big advantage of these two methods is that they are
efficient and readily available without any extra model tuning. For a given
set of variables $\mathbf{x}=(x_{1},\ldots,x_{d})$, the former seeks to
estimate the full conditional distributions
$p(x_{j}|x_{1},\ldots,x_{j-1},x_{j+1},x_{d})$ using CART (Breiman et al.,,
2017) and sample sequentially from these distributions. The NBE algorithm, on
the other hand, takes a parametric approach by assuming conditional
independence among $X_{j}$’s given a latent variable $Z$ and estimating the
conditional joint distribution $p(\mathbf{x}|z)$ by a mixture of models
(Gaussian for continuous variables and multinomial for categorical ones) using
the EM algorithm (Dempster et al.,, 1977). In Section 6, we compare the
performance of these two methods in estimation of $\mathrm{df}_{\rm R}$
through some simulation studies.

Here, one question remains to be answered: Can we synthesize
$\\{(\mathbf{x}_{\ast},y_{\ast})\\}_{b=1,\ldots,B}$ directly and evaluate out-
of-sample prediction error based on the synthetic data? To study its
feasibility, we consider a simple example below. Assume
$y_{i}=\sum_{j=1}^{5}x_{ij}+\varepsilon_{i}$, where $x_{ij}$’s and
$\varepsilon_{i}$ are i.i.d. $\mathcal{N}(0,1)$. We set $n=200$ and
$B=\text{1,000}$, and generate 500 synthetic datasets based on a fixed
training dataset $(\mathbf{X},\mathbf{y})$ using both the NBE and CART
algorithms. We apply the ordinary least squares method with all variables
included to fit a linear model. It is easy to show that the true out-of-sample
prediction error is given by $\mathrm{ErrR}_{\mathbf{X}}=1+{\rm
tr}[(\mathbf{X}^{\top}\mathbf{X})^{-1}]\approx 1.027$. However, we find that
both data synthesizing methods overestimate the true prediction error, with an
average of 3.872 (SD = 0.185) for the CART method and 5.878 (SD = 0.270) for
the NBE algorithm. For each method, a scatterplot of the synthesized
$y_{\ast}$’s against $x_{\ast}$’s is given in Figure 4. While the linear
pattern between $x_{\ast}$ and $y_{\ast}$ still remains as expected,
$y_{\ast}$ clearly has variance greater than 1. This explains the inflated
prediction error estimates. Our numerical results indicate that such over-
dispersion generally becomes more notable as the signal-to-noise ratio
increases. Structurally, the full conditional distributions estimated using
CART and the joint distribution estimated through the NBE algorithm under the
assumption of conditional independence may not provide a good approximation of
the true joint distribution of $(\mathbf{x},y)$. Therefore, it is not
recommended to estimate the prediction error based on synthetic data directly.

Figure 4: Response $y$ against its mean $\mu=\mu(\mathbf{x})$ from synthetic
data generated based on a fixed training dataset of size 200. The true data
generating mechanism is $y_{i}=\sum_{j=1}^{5}x_{ij}+\varepsilon_{i}$, where
$x_{ij}$’s and $\varepsilon_{i}$ are i.i.d. $\mathcal{N}(0,1)$. 1,000
synthetic examples are generated for each method.

On the other hand, this result doesn’t invalidate our approach to estimating
$\mathrm{df}_{\rm R}$ using synthetic data. In fact, our numerical results in
Section 6.1 demonstrate great robustness of $\widehat{\mathrm{df}}_{\rm R}$ to
synthetic data even when all features are moderately correlated. Also, note
that the estimated predictive model degrees of freedom
$\widehat{\mathrm{df}}_{\rm R}$ contributes to the prediction error estimate
$\widehat{\mathrm{wErrR}}$ through the term
$(2/n)\bar{w}\sigma_{\varepsilon}^{2}\widehat{\mathrm{df}}_{\rm R}$. When $n$
is large, the factor $2/n$ in front can effectively alleviate the impact of
discrepancy between the true and estimated joint covariate distribution.
Besides, we can further reduce this impact by wisely selecting variables for
data synthesis. For example, when multiple functions of a certain variable
$X_{j}$ are used as predictors in a model (e.g., $X_{j},\ldots,X_{j}^{k}$ in
polynomial regression), it is good practice to synthesize $X_{j}$ first and
propagate it to other transformations.

##  6 Numerical Studies

In this section, we evaluate the performance of the estimators developed in
Section 5 through numerical studies.

###  6.1 Simulations

Throughout our experiments, we set $n=60$ and generate 20 continuous variables
and 20 categorical variables in total ($d=40$). Let
$\mathbf{V}=(V_{1},\ldots,V_{20})^{\top}\in\mathbb{R}^{20}$ and
$\mathbf{C}=(C_{1},\ldots,C_{20})^{\top}\in\\{-1,+1\\}^{20}$ be the sets of
continuous and binary categorical variables respectively. To specify
relationship among features with different level of complexity, we consider
the following two settings.

  * •

Single component. Assume that
$\mathbf{V}\sim\mathcal{N}(\mathbf{0},\bm{\Sigma}_{\mathbf{V}})$ is
independent of $\mathbf{C}$, and $C_{1},\ldots,C_{20}$ are i.i.d. with
$\mathrm{P}(C_{1}=+1)=\mathrm{P}(C_{1}=-1)=1/2$.

  * •

Mixture of components. Let $Z\in\\{-1,+1\\}$ be a latent random variable such
that $\mathrm{P}(Z=+1)=\mathrm{P}(Z=-1)=1/2$. Assume that
$\mathbf{V}|Z\sim\mathcal{N}(Z(\log 1,\ldots,\log
20)^{\top},\bm{\Sigma}_{\mathbf{V}})$, $\mathrm{P}(C_{1}=+1)=1/2$ and
$\mathrm{P}(C_{j}=+1|C_{j-1},Z)=1/2+ZC_{j-1}/(4+\log j)$ for $j=2,\ldots,20$.

For the variance-covariance matrix $\bm{\Sigma}_{\mathbf{V}}$, we consider
three different cases: 1) I.I.D.: $\bm{\Sigma}_{\mathbf{V}}=\mathbf{I}_{20}$;
2) Equicorrelation:
$\bm{\Sigma}_{\mathbf{V}}=\frac{1}{2}\mathbf{I}_{20}+\frac{1}{2}\mathbf{1}_{20}\mathbf{1}_{20}^{\top}$;
and 3) Autoregressive correlation:
$\bm{\Sigma}_{\mathbf{V}}=(2^{-|i-j|})_{i,j=1,\ldots,20}$. Let
$\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})^{\top}$ be the design
matrix. We assume that the true data generating mechanism
$\mu(\cdot;\bm{\beta})$ is linear, i.e.,
$\mu(\mathbf{x}_{i};\bm{\beta})=\sum_{j=1}^{d}\beta_{j}x_{ij}$. For
convenience, we use $\beta^{\rm v}_{j}$ and $\beta^{\rm c}_{j}$ to denote the
coefficient $\beta_{j}$’s for the corresponding continuous and categorical
variables respectively. We set $\beta^{\rm v}_{j}=\beta^{\rm
c}_{j}/2=a\left(1-j/20\right)^{5}$ and choose $a$ such that
$\sum_{j=1}^{20}(\beta^{\rm v}_{j})^{2}=5$. In subset regression, we add
continuous and categorical variables alternatingly in descending order of
their coefficients, i.e., in the order $V_{1},C_{1},\ldots,V_{20},C_{20}$.

For the error variances, we assume $\sigma_{\varepsilon}^{2}=1$ and
$\tau(\mathbf{x})\propto(1+|\mu(\mathbf{x};\bm{\beta})|)^{2}$. We assume
$\tau(\mathbf{x})$ can be obtained for each $\mathbf{x}$ without knowing the
mean function $\mu$. We use $w=1/\tau$ as the weight function when fitting
weighted least squares models.

####  6.1.1 Estimation of $\mathrm{df}_{\rm R}$

We consider the CART method and the NBE learning algorithm described in
Section 5.2 to generate synthetic data. Figure 5 illustrates and compares the
results using these two procedures based on a random realization of
$\mathbf{X}$.

Figure 5: Estimated predicted model degrees of freedom
$\widehat{\mathrm{df}}_{\rm R}$ based on synthetic data generated using CART
and NBE algorithm. For each method, 100 synthetic datasets of size 1,000 are
generated based on a single replicate of $\mathbf{X}$ of size $n=60$ and
$d=40$ (20 continuous and 20 categorical). The true mean function is linear in
$\mathbf{x}$, and continuous and categorical variables are added alternatingly
in the descending order of their coefficients.

Each method is applied 100 times to generate different synthetic datasets used
for estimation of $\mathrm{df}_{\rm R}$. We can see a consistent pattern
across all 6 examined cases. Both methods exhibit stability with respect to
the choice of synthetic data, though the variance slightly inflates as the
subset size $p$ increases. The CART method tends to underestimate the true
predictive model degrees of freedom, especially under the I.I.D. setting.
These estimates lie between the predictive model degrees of freedom
$\mathrm{df}_{\rm R}$ and the classical model degrees of freedom
$\mathrm{df}_{\rm F}$. This suggests that the synthetic data generated using
this method do not fully mimic out-of-sample examples but act like a “hybrid”
of out-of-sample and in-sample data. This is not surprising since the method
samples from a set of full conditional distributions as mentioned in Section
5.2, which does not guarantee that the resulting synthetic data follow the
same joint distribution as the training data. On the other hand, despite a
slight tendency of overestimation in the equicorrelation setting, the NBE
algorithm performs remarkably well in estimating the true predictive model
degrees of freedom, even when the data do not meet the “naive” assumption that
all features are conditionally independent given the latent cluster. Lowd and
Domingos, (2005) pointed out that the naive Bayes mixture model is a special
case of Bayesian networks, which are commonly used in general high-dimensional
probability representation and estimation. In addition, the naive Bayes
mixture model is readily available and orders of magnitude faster than
Bayesian networks for sampling and inference. This enables us to estimate the
predictive model degrees of freedom in an accurate and efficient fashion.

####  6.1.2 Estimation of $\Delta B_{\mathbf{X}}$ and
$\mathrm{wErrR}_{\mathbf{X}}$

For simplicity, we only consider the setting of a mixture of components with
i.i.d. continuous covariates ($\bm{\Sigma}_{\mathbf{V}}=\mathbf{I}_{20}$). We
generate 500 replicates of $(\mathbf{X},\mathbf{y})$. Based on each replicate,
the estimators $\hat{\delta}$ and $\widehat{\mathrm{wErrR}}$ are calculated
along with the true excess bias $\Delta B_{\mathbf{X}}$ and prediction error
$\mathrm{wErrR}_{\mathbf{X}}$. In general, the excess bias estimator is
unbiased when $p$ is small, but becomes slightly biased upward as $p$
increases due to the truncation of $\hat{\delta}$,
$\hat{\delta}_{+}=\max(\hat{\delta},0)$ when $\hat{\delta}<0$.

Figure 6: Comparison of the four estimators $\widehat{\mathrm{wErrR}}$,
$\widehat{\mathrm{wErrR}}_{\rm loocv}$, $\widehat{\mathrm{wErrR}}_{\rm cv(5)}$
and $\widehat{\mathrm{wErrR}}_{\rm syn}$ in prediction error estimation. The
gray lines are the risk estimates based on each of the 500 random replicates
of $(\mathbf{X},\mathbf{y})$ generated under the mixture components setting
with $\bm{\Sigma}_{\mathbf{V}}=\mathbf{I}_{20}$. The red line marks the
average of these estimates, whereas the black line shows the average of the
true prediction ever over all replicates.

With the corrected excess bias estimator $\hat{\delta}_{+}$ and the estimated
predictive model degrees $\widehat{\mathrm{df}}_{\rm R}$ based on the NBE
algorithm, we form the prediction error estimator $\widehat{\mathrm{wErrR}}$.
For comparison, we also include three other estimators: the leave-one-out
cross validation error $\widehat{\mathrm{wErrR}}_{\rm loocv}$, the 5-fold
cross validation error $\widehat{\mathrm{wErrR}}_{\rm cv(5)}$, and the
prediction error estimator based on synthetic data, denoted by
$\widehat{\mathrm{wErrR}}_{\rm syn}$. The results are illustrated in Figure 6.
We can see that our estimator is very similar to the LOOCV error when $p$ is
relatively small, but has much lower variance when $p$ gets large. This
results from the use of $\hat{\delta}_{+}$ that trades off a small increase in
bias for a large reduction in variance. Note that 5-fold cross validation
doesn’t work very well in our example, since the training folds contain only
48 observations but a maximum of 40 covariates, which inevitably inflates bias
and variance simultaneously. $\widehat{\mathrm{wErrR}}_{\rm syn}$ performs the
worst due to the over-dispersion issue in the synthetic data.

Figure 7: Comparison of the four prediction error estimators
$\widehat{\mathrm{wErrR}}$, $\widehat{\mathrm{wErrR}}_{\rm loocv}$,
$\widehat{\mathrm{wErrR}}_{\rm cv(5)}$ and $\widehat{\mathrm{wErrR}}_{\rm
syn}$ in terms of mean squared error (left) and model selection (right) based
on 500 replicates of $(\mathbf{X},\mathbf{y})$ generated under the mixture
components setting with $\bm{\Sigma}_{\mathbf{V}}=\mathbf{I}_{20}$. For each
$p$, the mean squared error of an estimator is the average squared difference
between the estimate and true prediction error over all replicates. For each
replicate, we define $p_{\ast}$ and $\hat{p}_{\ast}$, respectively, as the
minimizer of the true and estimated prediction errors.

A more quantitative comparison of the four prediction error estimators is
given in Figure 7. The left panel shows the mean squared error (on log scale)
of each estimator as a function of $p$. Our estimator
$\widehat{\mathrm{wErrR}}$ beats the other three due to its much smaller
variance. In terms of model selection, let $p_{\ast}$ and $\hat{p}_{\ast}$ be
the minimizers of the true and estimated prediction errors respectively for
each replicate. The right panel shows the histogram of
$\hat{p}_{\ast}-p_{\ast}$. We can see that $\widehat{\mathrm{wErrR}}$ spreads
the least, though its center is a bit off to the left of 0, meaning that it
tends to select a slightly simpler model than the true model. Since out-of-
sample prediction generally involves more uncertainty than in-sample
prediction, it is reasonable to be slightly conservative in model selection.
On the other hand, the other three estimators are more likely to select larger
models.

Figure 8: Comparison of models selected by the four prediction error
estimators $\widehat{\mathrm{wErrR}}$, $\widehat{\mathrm{wErrR}}_{\rm loocv}$,
$\widehat{\mathrm{wErrR}}_{\rm cv(5)}$ and $\widehat{\mathrm{wErrR}}_{\rm
syn}$ based on 500 replicates of $(\mathbf{X},\mathbf{y})$ generated under the
mixture components setting with $\bm{\Sigma}_{\mathbf{V}}=\mathbf{I}_{20}$.
The red dashed lines show the average ratio
$\mathrm{wErrR}_{\mathbf{X}}(\hat{p}_{\ast})/\mathrm{wErrR}_{\mathbf{X}}(p_{\ast})$.

A simpler model is desirable in the sense that it is usually computationally
more efficient and has better interpretability. But if prediction is a top
priority, how much predictive accuracy do we lose in pursuit of such model
parsimony? To answer this question, we look at the ratio of the true
prediction error of the selected model to that of the optimal one, defined as
$\mathrm{wErrR}_{\mathbf{X}}(\hat{p}_{\ast})/\mathrm{wErrR}_{\mathbf{X}}(p_{\ast})$.
For a good model selection criterion, this ratio is expected to be close to 1.
As illustrated in Figure 8, our risk estimator $\widehat{\mathrm{wErrR}}$
exhibits the best performance in both the mean and standard deviation of the
ratio. The 5-fold cross-validation error $\widehat{\mathrm{wErrR}}_{\rm
cv(5)}$ also works comparably well, but its variance is slightly larger. By
contrast, the risk estimator based on synthetic data is very unstable and can
easily result in poor models.

####  6.1.3 WLS vs. OLS

So far in our numerical studies, we have assumed that the the weight function
is correctly specified. In practice, however, since error variances are
usually unknown, case weights need to be estimated and could be misspecified,
which may consequently affect statistical inference and decision. In this
section, we examine the effect of weight misspecification for the least
squares method. In particular, we focus on model selection and compare two
extreme weight schemes: the optimal weight and equal weight schemes under both
in-sample and out-of-sample prediction settings. For simplicity, we still
assume that the data are generated under the mixture components setting with
$\bm{\Sigma}_{\mathbf{V}}=\mathbf{I}_{20}$, and that
$\tau(\mathbf{x})\propto(1+|\mu(\mathbf{x};\bm{\beta})|)^{\lambda}$ with
$\lambda=2$.

We use $\widehat{\mathrm{wErrF}}$ and $\widehat{\mathrm{wErrR}}$ as a model
selection criterion in the in-sample and out-of-sample prediction settings
respectively. As shown in Figure 9, the ordinary least squares model tends to
select fewer variables than the weighted least squares model with the optimal
weight scheme in both scenarios. This difference in the size of the selected
model is generally the same for in-sample and out-of-sample predictions. In
terms of the variance of the estimated model size $\hat{p}_{\ast}$, the two
methods are quite close, with the weighted least squares model being slightly
more variable under the out-of-sample prediction setting.

Figure 9: Comparison of the selected model size $\hat{p}_{\ast}$ between the
ordinary least squares model and the weighted least squares model with the
optimal weight scheme based on 500 replicates of $(\mathbf{X},\mathbf{y})$
generated under the mixture components setting with
$\bm{\Sigma}_{\mathbf{V}}=\mathbf{I}_{20}$ and
$\tau(\mathbf{x})\propto(1+|\mu(\mathbf{x};\bm{\beta})|)^{2}$.
$\widehat{\mathrm{wErrF}}$ and $\widehat{\mathrm{wErrR}}$ are used as the
model selection criterion in the in-sample and out-of-sample prediction
settings respectively.

Since $\widehat{\mathrm{wErrR}}$ tends to favor a slightly more parsimonious
model as we have shown in Section 6.1.2, misspecification of weights could
lead to an oversimplified model. This could further result in a nontrivial
increase in the prediction error if the true risk curve to the left of the
optimal model size is steep. In our example, the in-sample and out-of-sample
prediction errors could be inflated by 12% and 50% respectively on average due
to weight misspecification. This issue will become more substantial with more
uneven error variances.

###  6.2 Real Data Analysis

We consider two real data problems in this section.

####  6.2.1 Cancer Mortality Data

We revisit the county-level cancer mortality data discussed in Luan et al.,
(2021). In addition to the 22 continuous features considered in the previous
analysis, we also take geographic variations into account by grouping all
states into 4 regions, namely Midwest, Northeast, South and West, defined in
Chapter 6 of the Geographic Areas Reference Manual (GARM) (Bureau of the
Census,, 1994). To restrict the scope of the analysis, we focus on the 19
states around the Great Lakes area illustrated in the left panel of Figure 11,
which amount to 1,148 observations (counties) in total. Minnesota is removed
due to data quality issues. Stratified random sampling based on each region is
used to create two datasets of size 150 and 998 for training and testing
respectively.

Figure 10: County-level cancer death rate of the states around the Great Lakes
area (left) and error variance as a function of $\log_{10}(\text{population})$
(right).

Figure 11: Standard deviation of the weighted prediction errors
$\sqrt{w_{i}}(y_{i}-\hat{y}_{i})$ by population deciles. For the OLS model,
$w_{i}=1$, whereas for the WLS model, $w_{i}$’s are calculated using the
smoothing spline model shown in the right panel of Figure 11.

We use observations from the remaining states to estimate the variance
function $\tau$ and determine the variable orderings in subset regression. For
estimation of $\tau$, we first fit an OLS model to the death rate with all
features included. We then split the residuals of the model into 10 quantile
groups by the the population size and calculate the variances of the residuals
within each group. We then estimate the variance function by fitting a
smoothing spline model to these sample variances as a function of the group
median of the logarithm of the population size, as illustrated in the right
panel of Figure 11. For ordering of variables, we bootstrap data for 500
times. With each bootstrap sample, we construct a variable sequence by adding
the variable that reduces the weighted residual sum of squares most until all
variables are included. The variables are then ordered by their average ranks
over all 500 bootstrap samples.

Figure 12: Estimated predictive model degrees of freedom $\mathrm{df}_{\rm R}$
and the excess bias $\Delta B_{\mathbf{X}}$ of the ordinary and weighted least
squares models averaged over 500 random splits of the data into training and
test sets. The shaded bands show the 2.5th and 97.5th percentiles of the
estimates for each $p$.

We use $w=1/\tau$ as the weight function in the weighted least squares model.
We first compare the OLS and WLS models based on 500 random splits of the data
into training and test sets. Figure 11 shows the standard deviations of the
weighted prediction errors $\sqrt{w_{i}}(y_{i}-\hat{y}_{i})$ by population
deciles. For the OLS model, these errors are the unweighted ones due to
$w_{i}=1$. It is clear that error variances become much more homogeneous when
case weights are incorporated. In Figure 12, we demonstrate the average
estimated predictive model degrees of freedom and excess bias of the two
models. We can see that the WLS model has both smaller model complexity and
excess bias than the OLS model, though the difference in the former is rather
negligible. In terms of model selection, the two models are fairly close, with
the OLS model choosing 9.07 (SD = 1.65) variables on average and 9.08 (SD =
1.78) for the WLS model. For a particular random split of the data, a
comparison between the WLS and OLS models in coefficient estimation is given
in Appendix B.1.

Table 1: Average performance of different variable selection methods based on
500 random data partitions. $p_{\ast}$ and $\hat{p}_{\ast}$ are the optimal
model size indicated by the test data and variable selection criterion
respectively. Numbers in parentheses are the standard deviations. $n=150$
---
Method | $\widehat{\mathrm{wErrR}}$ | $C_{p}$ | LOOCV | 5-fold CV | AIC | BIC
$\hat{p}_{\ast}/p_{\ast}$ |  $0.968$ ($0.450$) |  $1.191$ ($0.453$) |  $1.033$ ($0.458$) |  $1.062$ ($0.460$) |  $1.247$ ($0.469$) |  $0.703$ ($0.317$)
$\mathrm{wErrR}_{\rm te}(\hat{p}_{\ast})/\mathrm{wErrR}_{\rm te}(p_{\ast})$ |  $1.071$ ($0.123$) |  $1.071$ ($0.117$) |  $1.107$ ($0.123$) |  $1.064$ ($0.109$) |  $1.073$ ($0.123$) |  $1.085$ ($0.125$)
$n=40$
Method | $\widehat{\mathrm{wErrR}}$ | $C_{p}$ | LOOCV | 5-fold CV | AIC | BIC
$\hat{p}_{\ast}/p_{\ast}$ |  $0.663$ ($0.319$) |  $1.492$ ($0.707$) |  $1.230$ ($0.616$) |  $1.155$ ($0.554$) |  $2.146$ ($0.840$) |  $1.192$ ($0.612$)
$\mathrm{wErrR}_{\rm te}(\hat{p}_{\ast})/\mathrm{wErrR}_{\rm te}(p_{\ast})$ |  $1.370$ ($0.530$) |  $1.936$ ($1.719$) |  $1.689$ ($1.440$) |  $1.604$ ($1.105$) |  $2.589$ ($2.529$) |  $1.703$ ($1.386$)

In Table 1, we compare the average performance of $\widehat{\mathrm{wErrR}}$
in variable selection with five other commonly used criteria using the
estimated weight scheme. Here, we consider two scenarios. In the first one, we
set the sample size of the training data equal to 150. On average, models
selected by $\widehat{\mathrm{wErrR}}$ and BIC contain fewer variables than
the optimal ones on the test data. However, since the sample size is much
larger than the total number of variables in this case, the model sizes
identified by each method are not too different. In terms of out-of-sample
prediction, all criteria have similar performance.

If we reduce the sample size of training data to 40, more differences can be
observed among these six criteria. In general, $\widehat{\mathrm{wErrR}}$
outperforms others since it selects the least number of variables while having
the smallest prediction error and lowest variability. The model parsimony from
$\widehat{\mathrm{wErrR}}$ is due to the fact that model complexity increases
quickly when the number of features gets close to the sample size. Since out-
of-sample prediction in this scenario usually involves much higher
uncertainty, it is reasonable to only include the most relevant variables in
the model.

####  6.2.2 King County House Sales Data

We continue the discussion in Section 6.1.3 on the impact of weight scheme on
model selection using the King County house sales data, which is publicly
available at https://www.kaggle.com/harlfoxem/housesalesprediction. It
contains sale prices and property features of 21,597 houses sold in King
County, WA between May 2014 and May 2015. We split the data by year (2014 and
2015) for exploration and modeling purposes. We take square root of the
housing price as the response variable. After selecting 12 initial predictors,
we apply appropriate transformations. Further, using the same data, we
estimate the variance function $\tau$ based on the residuals from the OLS
model with all features included, as shown by the red curve in Figure 13.
Details of the procedure are provided in Appendix B.2.

Figure 13: Variances of residuals within each percentile group of the fitted
values of the OLS model with the square root of the housing price as a
response and all 12 predictors for King County house sales data.

For the modeling part, we randomly select a subset of size $n=200$ to fit a
set of WLS models using the estimated variance function $\tau$ along the
variable sequence. We repeat the procedure 500 times. For comparison, we also
fit an OLS model with the same set of variables for each subset size $p$. We
find that, under both in-sample and out-of-sample prediction settings, the two
models identify the same subset size for more than half of the time, which is
apparently different from what we have seen from the simulation study in
Section 6.1.3. This suggests that the training error and excess bias in this
example may have larger influence on model selection. What is consistent with
the simulation is that, when the difference does show up, it is most likely
that the WLS model selects more variables than the OLS model. In a few extreme
cases, this difference in the size of the selected model can reach up to 6
variables.

For a close-up examination of the above findings, we decompose
$\widehat{\mathrm{wErrF}}$ and $\widehat{\mathrm{wErrR}}$ into the training
error, excess bias and excess variance, and study the median percentages of
these components of the selected models. A figure of the decomposition is
provided in Appendix B.2. We find that the training error is the dominant
component in the risk estimate under both in-sample and out-of-sample
prediction settings, even when all variables are included. As a result, model
selection is greatly impacted by the training error, especially in the in-
sample scenario where the term explains more than 90% of the risk. Since the
training error is decreasing in the subset size, both WLS and OLS models tend
to favor the full model in this scenario. On the other hand, for out-of-sample
prediction, the excess variance plays a relatively more important role in
model selection (around 20-25% of the risk when all variables are included
against 9-12% in the in-sample prediction setting), which leads to a reduction
in the size of the selected model. Comparing the WLS and OLS models, the
former has a constantly higher percentage of training error and lower
percentage of excess variance than the latter. This explains why the WLS model
has a slight tendency of selecting a few more variables. In addition, we also
see a lower percentage of the excess bias in $\widehat{\mathrm{wErrR}}$ of the
WLS model. Since the overall percentage of this component is quite small in
this example, we expect its impact on model selection would be small as well.

##  7 Conclusion

In this work, we extend the existing theories on model complexity to a
heteroscedastic setting for linear regression procedures. In particular, we
generalize the classical model degrees of freedom and predictive model degrees
of freedom so that models that incorporate different case weights can be
distinguished appropriately. As a typical example, we explore the two extended
measures of model degrees of freedom of the weighted least squares method in
depth. Interestingly, we find that, in the subset regression setting, the two
measures are a shrunken version of those of the ordinary least squares method,
and the amount of shrinkage depends on how uneven the error variances are.
Furthermore, we set up connections between the two measures of model degrees
of freedom and weight schemes both for model fitting and model evaluation. Our
results provide valuable insights on how model complexity is affected by case
weights, and further provide useful guidance in risk estimation and model
selection.

For future study in model complexity, it will be interesting to extend the
current framework to a generalized linear model setting. Note that for
exponential family distributions other than normal, variance usually shares
the same set of parameters with the mean and varies from sample to sample. As
a result, our work on heteroscedastic linear regression models holds promise
for the possibility of the extension and thus deserves further investigation.

#### Acknowledgements

We thank Steven MacEachern for comments on the earlier version of this work.
This research was supported in part by the National Science Foundation Grants
DMS-17-12580, DMS-17-21445 and DMS-20-15490.

## References

  * Anscombe, (1961) Anscombe, F. J. (1961).  Examination of residuals.  In Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics, pages 1–36. University of California Press.
  * Bickel, (1978) Bickel, P. J. (1978).  Using residuals robustly I: Tests for heteroscedasticity, nonlinearity.  Annals of Statistics, 6(2):266–291.
  * Box and Hill, (1974) Box, G. E. and Hill, W. J. (1974).  Correcting inhomogeneity of variance with power transformation weighting.  Technometrics, 16(3):385–389.
  * Breiman et al., (2017) Breiman, L., Friedman, J. H., Olshen, R. A., and Stone, C. J. (2017).  Classification and Regression Trees.  Routledge.
  * Bureau of the Census, (1994) Bureau of the Census (1994).  Geographic Areas Reference Manual.  U.S. Department of Commerce, Economics and Statistics Administration.
  * Carroll and Ruppert, (1981) Carroll, R. J. and Ruppert, D. (1981).  On robust tests for heteroscedasticity.  Annals of Statistics, 9(1):206–210.
  * Craven and Wahba, (1978) Craven, P. and Wahba, G. (1978).  Smoothing noisy data with spline functions.  Numerische Mathematik, 31(4):377–403.
  * Dempster et al., (1977) Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).  Maximum likelihood from incomplete data via the EM algorithm.  Journal of the Royal Statistical Society: Series B (Methodological), 39(1):1–22.
  * Efron, (1986) Efron, B. (1986).  How biased is the apparent error rate of a prediction rule?  Journal of the American Statistical Association, 81(394):461–470.
  * Efron, (2004) Efron, B. (2004).  The estimation of prediction error: Covariance penalties and cross-validation.  Journal of the American Statistical Association, 99(467):619–632.
  * Lowd and Domingos, (2005) Lowd, D. and Domingos, P. (2005).  Naive Bayes models for probability estimation.  In Proceedings of the 22nd International Conference on Machine Learning, pages 529–536.
  * Luan et al., (2021) Luan, B., Lee, Y., and Zhu, Y. (2021).  Predictive model degrees of freedom in linear regression.  arXiv preprint arXiv:2106.15682.
  * Mallows, (1973) Mallows, C. (1973).  Some comments on $C_{p}$.  Technometrics, 15(4):661–675.
  * Mardia et al., (1979) Mardia, K., Kent, J., and Bibby, J. (1979).  Multivariate Analysis.  Academic Press.
  * Pivaro et al., (2017) Pivaro, G. F., Kumar, S., Fraidenraich, G., and Dias, C. F. (2017).  On the exact and approximate eigenvalue distribution for sum of Wishart matrices.  IEEE Transactions on Vehicular Technology, 66(11):10537–10541.
  * Reiter, (2005) Reiter, J. P. (2005).  Using CART to generate partially synthetic public use microdata.  Journal of Official Statistics, 21(3):441.
  * Tibshirani and Hastie, (1987) Tibshirani, R. and Hastie, T. (1987).  Local likelihood estimation.  Journal of the American Statistical Association, 82(398):559–567.
  * Ye, (1998) Ye, J. (1998).  On measuring and correcting the effects of data mining and model selection.  Journal of the American Statistical Association, 93(441):120–131.

##  Appendix A Proofs

###  A.1 Derivation of Equations (3.4) and (3.5)

Consider a linear procedure with hat matrix $\mathbf{H}$ and hat vector
$\mathbf{h}_{\ast}$ at $\mathbf{x}_{\ast}$. We first look at the bias-variance
decomposition for $\mathrm{wErrR}_{\mathbf{X}}$. Note that

$\mathrm{wErrR}_{\mathbf{X}}=\mathrm{E}[w_{\ast}(y_{\ast}-\mathbf{h}_{\ast}^{\top}\mathbf{y})^{2}|\mathbf{X}]=\mathrm{E}[\mathrm{E}[(w_{\ast}(y_{\ast}-\mathbf{h}_{\ast}^{\top}\mathbf{y})^{2}|\mathbf{x}_{\ast},\mathbf{X}]|\mathbf{X}],$

where the inner and outer expectations are taken with respect to
$(y_{\ast},\mathbf{y})$ and $\mathbf{x}_{\ast}$ respectively. For the inner
expectation, using the fact that $(\mathbf{x}_{\ast},y_{\ast})$ is independent
of $(\mathbf{X},\mathbf{y})$, we have

$\displaystyle\mathrm{E}[(w_{\ast}(y_{\ast}-\mathbf{h}_{\ast}^{\top}\mathbf{y})^{2}|\mathbf{x}_{\ast},\mathbf{X}]$
$\displaystyle=w_{\ast}\mathrm{E}[(y_{\ast}-\mu_{\ast}+\mu_{\ast}-\mathbf{h}_{\ast}^{\top}\bm{\mu}+\mathbf{h}_{\ast}^{\top}\bm{\mu}-\mathbf{h}_{\ast}^{\top}\mathbf{y})^{2}|\mathbf{x}_{\ast},\mathbf{X}]$
$\displaystyle=w_{\ast}(\mathrm{E}[(y_{\ast}-\mu_{\ast})^{2}|\mathbf{x}_{\ast}]+(\mu_{\ast}-\mathbf{h}_{\ast}^{\top}\bm{\mu})^{2}+\mathrm{E}[(\mathbf{h}_{\ast}^{\top}\mathbf{y}-\mathbf{h}_{\ast}^{\top}\bm{\mu})^{2}|\mathbf{x}_{\ast},\mathbf{X}])$
$\displaystyle=\sigma_{\varepsilon}^{2}w_{\ast}\tau_{\ast}+w_{\ast}(\mu_{\ast}-\mathbf{h}_{\ast}^{\top}\bm{\mu})^{2}+\sigma_{\varepsilon}^{2}w_{\ast}\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2},$

where $\mathbf{T}=\mathrm{diag}(\tau_{1},\ldots,\tau_{n})$. Then taking
another expectation of the above result with respect to $\mathbf{x}_{\ast}$
yields (3.4) immediately.

For $\mathrm{ErrT}_{\mathbf{X}}$, it follows that

$\displaystyle\mathrm{wErrT}_{\mathbf{X}}$
$\displaystyle=\frac{1}{n}\mathrm{E}(\|\mathbf{y}-\mathbf{H}\mathbf{y}\|_{\mathbf{W}}^{2}|\mathbf{X})$
$\displaystyle=\frac{1}{n}\|\bm{\mu}-\mathbf{H}\bm{\mu}\|_{\mathbf{W}}^{2}+\frac{1}{n}\sigma_{\varepsilon}^{2}{\rm
tr}[(\mathbf{I}_{n}-\mathbf{H})^{\top}\mathbf{W}(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}]$
$\displaystyle=\frac{1}{n}\sigma_{\varepsilon}^{2}{\rm
tr}(\mathbf{W}\mathbf{T})+\frac{1}{n}\|\bm{\mu}-\mathbf{H}\bm{\mu}\|_{\mathbf{W}}^{2}+\frac{1}{n}\sigma_{\varepsilon}^{2}{\rm
tr}(\mathbf{H}^{\top}\mathbf{W}\mathbf{H}\mathbf{T}-\mathbf{H}^{\top}\mathbf{W}\mathbf{T}-\mathbf{W}\mathbf{H}\mathbf{T})$
$\displaystyle=\frac{1}{n}\sigma_{\varepsilon}^{2}\sum_{i=1}^{n}{w_{i}\tau_{i}}+\frac{1}{n}\|\bm{\mu}-\mathbf{H}\bm{\mu}\|_{\mathbf{W}}^{2}+\frac{1}{n}\sigma_{\varepsilon}^{2}{\rm
tr}(\mathbf{H}^{\top}\mathbf{W}\mathbf{H}\mathbf{T}-2\mathbf{H}\mathbf{T}\mathbf{W}).$

This completes the proof of (3.5). In particular, when
$w(\cdot)=\tau(\cdot)=1$, (3.4) and (3.5) reduce to (2.3) and (2.4)
respectively.

###  A.2 Proof of Proposition 3.1

###### Proof.

Since $\hat{\mu}_{\ast}=\mathbf{h}_{\ast}^{\top}\mathbf{y}$ and
$\hat{\mu}_{j}=\mathbf{h}_{j}^{\top}\mathbf{y}$, we have

$\frac{\mathrm{Cov}(y_{i},\hat{\mu}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})}{\sigma_{\varepsilon}^{2}\tau_{i}}=\frac{\partial\mathrm{E}(\hat{\mu}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})}{\partial\mu_{i}}=h_{\ast,i}\text{
and
}\frac{\mathrm{Cov}(y_{i},\hat{\mu}_{j}|\mathbf{x}_{\ast},\mathbf{X})}{\sigma_{\varepsilon}^{2}\tau_{i}}=\frac{\partial\mathrm{E}(\hat{\mu}_{j}|\mathbf{X})}{\partial\mu_{i}}=h_{ji}.$

As a result,

$\sum_{i=1}^{n}\tau_{i}\mathrm{E}\left(w_{\ast}\frac{\mathrm{Cov}^{2}(y_{i},\hat{\mu}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})}{(\sigma_{\varepsilon}^{2}\tau_{i})^{2}}\Bigg{|}\mathbf{X}\right)=\sum_{i=1}^{n}\tau_{i}\mathrm{E}\left(w_{\ast}\left(\frac{\partial\mathrm{E}(\hat{\mu}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})}{\partial\mu_{i}}\right)^{2}\Bigg{|}\mathbf{X}\right)=\mathrm{E}(w_{\ast}\|\mathbf{h}_{\ast}\|_{\mathbf{T}}^{2}|\mathbf{X}),$

and

$\sum_{i=1}^{n}\tau_{i}\sum_{j=1}^{n}w_{j}\frac{\mathrm{Cov}^{2}(y_{i},\hat{\mu}_{j}|\mathbf{X})}{(\sigma_{\varepsilon}^{2}\tau_{i})^{2}}=\sum_{i=1}^{n}\tau_{i}\sum_{j=1}^{n}w_{j}\left(\frac{\partial\mathrm{E}(\hat{\mu}_{j}|\mathbf{X})}{\partial\mu_{i}}\right)^{2}={\rm
tr}(\mathbf{W}\mathbf{H}\mathbf{T}\mathbf{H}^{\top}).$

The theorem then follows. ∎

###  A.3 Proof of Theorem 4.1

###### Proof.

Let
$\mathbf{U}=\frac{1}{\bar{w}}\mathbf{W}^{\frac{1}{2}}\mathbf{T}\mathbf{W}^{\frac{1}{2}}$.
We can express $\mathrm{df}_{\rm F}(p)$ in terms of $\tilde{\mathbf{X}}_{p}$
and $\mathbf{U}$ as

$\mathrm{df}_{\rm F}(p)={\rm
tr}[(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}\mathbf{U}\tilde{\mathbf{X}}_{p}].$

Note that
$\tilde{\mathbf{X}}_{p+1}=(\tilde{\mathbf{X}}_{p},\tilde{\mathbf{x}}_{(p+1)})$.
By the inversion formula for a 2 by 2 partitioned matrix,

$(\tilde{\mathbf{X}}_{p+1}^{\top}\tilde{\mathbf{X}}_{p+1})^{-1}=\begin{pmatrix}\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p}&\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{x}}_{(p+1)}\\\
\tilde{\mathbf{x}}_{(p+1)}^{\top}\tilde{\mathbf{X}}_{p}&\tilde{\mathbf{x}}_{(p+1)}^{\top}\tilde{\mathbf{x}}_{(p+1)}\end{pmatrix}^{-1}=\left(\begin{array}[]{cc}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}+\frac{1}{b}\mathbf{s}\mathbf{s}^{\top}&-\frac{1}{b}\mathbf{s}\\\
-\frac{1}{b}\mathbf{s}^{\top}&\frac{1}{b}\end{array}\right),$

where
$\mathbf{s}=(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{x}}_{(p+1)}$
and
$b=\tilde{\mathbf{x}}_{(p+1)}^{\top}(\mathbf{I}_{n}-\tilde{\mathbf{X}}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top})\tilde{\mathbf{x}}_{(p+1)}$.
Letting
$\tilde{\mathbf{H}}=\tilde{\mathbf{X}}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}$
and
$\tilde{\mathbf{r}}_{(p+1)}=(\mathbf{I}_{n}-\tilde{\mathbf{H}})\tilde{\mathbf{x}}_{(p+1)}$,
we have $b=\|\tilde{\mathbf{r}}_{(p+1)}\|^{2}$. Then

$\displaystyle\mathrm{df}_{\rm F}(p+1)$ $\displaystyle={\rm
tr}[(\tilde{\mathbf{X}}_{p+1}^{\top}\tilde{\mathbf{X}}_{p+1})^{-1}\tilde{\mathbf{X}}_{p+1}^{\top}\mathbf{U}\tilde{\mathbf{X}}_{p+1}]$
$\displaystyle={\rm
tr}\left[\left((\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}+\frac{1}{b}\mathbf{s}\mathbf{s}^{\top}\right)\tilde{\mathbf{X}}_{p}^{\top}\mathbf{U}\tilde{\mathbf{X}}_{p}-\frac{1}{b}\mathbf{s}\tilde{\mathbf{x}}_{(p+1)}^{\top}\mathbf{U}\tilde{\mathbf{X}}_{p}\right]$
$\displaystyle\qquad-\frac{1}{b}\mathbf{s}^{\top}\tilde{\mathbf{X}}_{p}^{\top}\mathbf{U}\tilde{\mathbf{x}}_{(p+1)}+\frac{1}{b}\tilde{\mathbf{x}}_{(p+1)}^{\top}\mathbf{U}\tilde{\mathbf{x}}_{(p+1)}$
$\displaystyle=\mathrm{df}_{\rm
F}(p)+\frac{1}{b}\tilde{\mathbf{x}}_{(p+1)}^{\top}\tilde{\mathbf{X}}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}\mathbf{U}\tilde{\mathbf{X}}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{x}}_{(p+1)}$
$\displaystyle\qquad\qquad\,-\frac{2}{b}\tilde{\mathbf{x}}_{(p+1)}^{\top}\tilde{\mathbf{X}}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}\mathbf{U}\tilde{\mathbf{x}}_{(p+1)}$
$\displaystyle\qquad\qquad\,+\frac{1}{b}\tilde{\mathbf{x}}_{(p+1)}^{\top}\mathbf{U}\tilde{\mathbf{x}}_{(p+1)}$
$\displaystyle=\mathrm{df}_{\rm
F}(p)+\frac{1}{b}\tilde{\mathbf{x}}_{(p+1)}^{\top}(\mathbf{I}_{n}-\tilde{\mathbf{H}})\mathbf{U}(\mathbf{I}_{n}-\tilde{\mathbf{H}})\tilde{\mathbf{x}}_{(p+1)}$
$\displaystyle=\mathrm{df}_{\rm
F}(p)+\frac{1}{b}\tilde{\mathbf{r}}_{(p+1)}^{\top}\mathbf{U}\tilde{\mathbf{r}}_{(p+1)}$
$\displaystyle=\mathrm{df}_{\rm
F}(p)+\frac{\|\tilde{\mathbf{r}}_{(p+1)}\|_{\mathbf{U}}^{2}}{\|\tilde{\mathbf{r}}_{(p+1)}\|^{2}}.$

∎

###  A.4 Proof of Theorem 4.2

###### Proof.

Following the notation in A.3, define
$\mathbf{G}_{p}=(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}\mathbf{U}\tilde{\mathbf{X}}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}$.
Then $\mathrm{df}_{\rm R}$ can be written as

$\mathrm{df}_{\rm R}(p)=\frac{1}{2}\mathrm{df}_{\rm F}(p)+\frac{n}{2}{\rm
tr}(\mathbf{G}_{p}\bm{\Sigma}_{p}).$

It can be checked that
$\mathbf{G}_{p+1}=(\tilde{\mathbf{X}}_{p+1}^{\top}\tilde{\mathbf{X}}_{p+1})^{-1}\tilde{\mathbf{X}}_{p+1}^{\top}\mathbf{U}\tilde{\mathbf{X}}_{p+1}(\tilde{\mathbf{X}}_{p+1}^{\top}\tilde{\mathbf{X}}_{p+1})^{-1}$
has the form

$\mathbf{G}_{p+1}=\begin{pmatrix}\mathbf{F}&\mathbf{g}\\\
\mathbf{g}^{\top}&f\end{pmatrix},$

where

$\displaystyle\mathbf{F}$
$\displaystyle=\mathbf{G}_{p}-\frac{1}{b}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}\mathbf{U}\tilde{\mathbf{r}}_{(p+1)}\mathbf{s}^{\top}-\frac{1}{b}\mathbf{s}\tilde{\mathbf{r}}_{(p+1)}^{\top}\mathbf{U}\tilde{\mathbf{X}}_{p}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}+\frac{\|\tilde{\mathbf{r}}_{(p+1)}\|_{\mathbf{U}}^{2}}{b^{2}}\mathbf{s}\mathbf{s}^{\top},$
$\displaystyle\mathbf{g}$
$\displaystyle=\frac{1}{b}(\tilde{\mathbf{X}}_{p}^{\top}\tilde{\mathbf{X}}_{p})^{-1}\tilde{\mathbf{X}}_{p}^{\top}\mathbf{U}\tilde{\mathbf{r}}_{(p+1)}-\frac{\|\tilde{\mathbf{r}}_{(p+1)}\|_{\mathbf{U}}^{2}}{b^{2}}\mathbf{s},$
$\displaystyle f$
$\displaystyle=\frac{\|\tilde{\mathbf{r}}_{(p+1)}\|_{\mathbf{U}}^{2}}{b^{2}}.$

Note that

$\mathbf{F}=\mathbf{G}_{p}-(\mathbf{g}+f\mathbf{s})\mathbf{s}^{\top}-\mathbf{s}(\mathbf{g}+f\mathbf{s})^{\top}+f\mathbf{s}\mathbf{s}^{\top}=\mathbf{G}_{p}-\mathbf{g}\mathbf{s}^{\top}-\mathbf{s}\mathbf{g}^{\top}-f\mathbf{s}\mathbf{s}^{\top}.$

Then

$\displaystyle\mathrm{df}_{\rm R}(p+1)$
$\displaystyle=\frac{1}{2}\mathrm{df}_{\rm F}(p+1)+\frac{n}{2}{\rm
tr}(\mathbf{G}_{p+1}\bm{\Sigma}_{p+1})$
$\displaystyle=\frac{1}{2}\mathrm{df}_{\rm F}(p+1)+\frac{n}{2}[{\rm
tr}(\mathbf{F}\bm{\Sigma}_{p})+2\mathbf{g}^{\top}\bm{\phi}_{p,p+1}+f\sigma_{p+1}^{2}]$
$\displaystyle=\mathrm{df}_{\rm R}(p)+\frac{1}{2}\Delta\mathrm{df}_{\rm
F}(p)+\frac{n}{2}(f\sigma_{p+1}^{2}-f\mathbf{s}^{\top}\bm{\Sigma}_{p}\mathbf{s}-2\mathbf{g}^{\top}\bm{\Sigma}_{p}\mathbf{s}+2\mathbf{g}^{\top}\bm{\phi}_{p,p+1}).$

Now we show that

$\displaystyle\quad\;\frac{n}{2}(f\sigma_{p+1}^{2}-f\mathbf{s}^{\top}\bm{\Sigma}_{p}\mathbf{s}-2\mathbf{g}^{\top}\bm{\Sigma}_{p}\mathbf{s}+2\mathbf{g}^{\top}\bm{\phi}_{p,p+1})$
(A.1) $\displaystyle=\frac{n}{2}\Delta\mathrm{df}_{\rm
F}(p)\frac{(\|\tilde{\mathbf{x}}_{(p+1)}-\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\|_{\mathbf{M}}^{2}+a_{p,p+1})}{\|\tilde{\mathbf{r}}_{(p+1)}\|^{2}}$
$\displaystyle\qquad\qquad-n\frac{\langle\tilde{\mathbf{r}}_{(p+1)},\tilde{\mathbf{x}}_{(p+1)}-\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\rangle_{\mathbf{U}\mathbf{M}}}{\|\tilde{\mathbf{r}}_{(p+1)}\|^{2}}.$

Consider expressing the right hand side of (A.1) using $f$, $\mathbf{g}$ and
$\mathbf{s}$. For the first term, note that

$\displaystyle\|\tilde{\mathbf{x}}_{(p+1)}\|_{\mathbf{M}}^{2}$
$\displaystyle=\mathbf{s}^{\top}\bm{\Sigma}_{p}\mathbf{s},$
$\displaystyle\|\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\|_{\mathbf{M}}^{2}$
$\displaystyle=\bm{\phi}_{p,p+1}^{\top}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1},$
$\displaystyle\langle\tilde{\mathbf{x}}_{(p+1)},\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\rangle_{\mathbf{M}}$
$\displaystyle=\mathbf{s}^{\top}\bm{\phi}_{p,p+1}.$

Then

$\displaystyle\quad\;\frac{n}{2}\Delta\mathrm{df}_{\rm
F}(p)\frac{(\|\tilde{\mathbf{x}}_{(p+1)}-\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\|_{\mathbf{M}}^{2}+a_{p,p+1})}{\|\tilde{\mathbf{r}}_{(p+1)}\|^{2}}$
(A.2)
$\displaystyle=\frac{nf}{2}(\|\tilde{\mathbf{x}}_{(p+1)}\|_{\mathbf{M}}^{2}-2\langle\tilde{\mathbf{x}}_{(p+1)},\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\rangle_{\mathbf{M}}+\|\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\|_{\mathbf{M}}^{2}+a_{p,p+1})$
$\displaystyle=\frac{nf}{2}(\mathbf{s}^{\top}\bm{\Sigma}_{p}\mathbf{s}-2\mathbf{s}^{\top}\bm{\phi}_{p,p+1}+\sigma_{p+1}^{2}).$

For the second term, it can be similarly shown that

$n\frac{\langle\tilde{\mathbf{r}}_{(p+1)},\tilde{\mathbf{x}}_{(p+1)}-\tilde{\mathbf{X}}_{p}\bm{\Sigma}_{p}^{-1}\bm{\phi}_{p,p+1}\rangle_{\mathbf{U}\mathbf{M}}}{\|\tilde{\mathbf{r}}_{(p+1)}\|^{2}}=n(\mathbf{g}+f\mathbf{s})^{\top}(\bm{\Sigma}_{p}\mathbf{s}-\bm{\phi}_{p,p+1})$
(A.3)

Subtracting (A.3) from (A.2) yields (A.1) immediately. ∎

###  A.5 Derivation of Approximation (4.5)

Assume that $\mathbf{x}_{\ast}$ is multivariate normal with
$\mathrm{E}(\mathbf{x}_{\ast})=\mathbf{0}$ and
$\mathrm{Var}(\mathbf{x}_{\ast})=\bm{\Sigma}_{\mathbf{x}}$. Let
$\mathbf{u}_{1},\ldots,\mathbf{u}_{n}$ be i.i.d. random variables that are
independent of $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$. Define
$w_{i}=1/\tau(\mathbf{u}_{i})$ and $w_{\ast}=1/\tau(\mathbf{u}_{\ast})$. Then,
we have

$\bm{\Sigma}_{\mathbf{x}}=\frac{\bar{w}}{\mathrm{E}(w_{\ast})}\frac{1}{\bar{w}}\mathrm{E}(w_{\ast}\mathbf{x}_{\ast}\mathbf{x}_{\ast}^{\top})=\frac{\bar{w}}{\mathrm{E}(w_{\ast})}\bm{\Sigma}_{p}.$

Let $\mathbf{z}_{i}=\bm{\Sigma}_{\mathbf{x}}^{-\frac{1}{2}}\mathbf{x}_{i}$ and
$\mathbf{Z}=(\mathbf{z}_{1},\ldots,\mathbf{z}_{n})^{\top}$. Note that
$\mathbf{z}_{i}$’s are i.i.d. $\mathcal{N}(\mathbf{0},\mathbf{I}_{p})$. Thus,
$\mathbf{z}_{i}\mathbf{z}_{i}^{\top}$’s are independent Wishart random
matrices with degrees of freedom 1 and scale matrix $\mathbf{I}_{p}$, denoted
by $\mathcal{W}(1,\mathbf{I}_{p})$. Note that, when $w\propto 1/\tau$,

$\mathrm{df}_{\rm R}(p)=\frac{\zeta_{\tau}}{2}(p+n\bar{w}{\rm
tr}[(\mathbf{X}_{p}^{\top}\mathbf{W}\mathbf{X}_{p})^{-1}\bm{\Sigma}_{p})])$

and

$n\bar{w}{\rm
tr}[(\mathbf{X}^{\top}\mathbf{W}\mathbf{X})^{-1}\bm{\Sigma}_{p})]=\frac{\mathrm{E}(w_{\ast})}{\bar{w}}n\bar{w}{\rm
tr}[(\mathbf{Z}^{\top}\mathbf{W}\mathbf{Z})^{-1}]=\frac{\mathrm{E}(w_{\ast})}{\bar{w}}{\rm
tr}\left[\left(\sum_{i=1}^{n}\frac{w_{i}}{n\bar{w}}\mathbf{z}_{i}\mathbf{z}_{i}^{\top}\right)^{-1}\right].$

When $n$ is large, Pivaro et al., (2017) suggest to approximate
$\sum_{i=1}^{n}\frac{w_{i}}{n\bar{w}}\mathbf{z}_{i}\mathbf{z}_{i}^{\top}$ by

$\sum_{i=1}^{n}\frac{w_{i}}{n\bar{w}}\mathbf{z}_{i}\mathbf{z}_{i}^{\top}\stackrel{{\scriptstyle\rm
d}}{{\approx}}\frac{\mathbf{V}}{n_{\tau}},\quad
n_{\tau}=\left[\frac{1}{\sum_{i=1}^{n}\left(\frac{w_{i}}{n\bar{w}}\right)^{2}}\right]=\left[\frac{\left(\sum_{i=1}^{n}1/\tau_{i}\right)^{2}}{\sum_{i=1}^{n}(1/\tau_{i})^{2}}\right],$

where $\mathbf{V}\sim\mathcal{W}(n_{\tau},\mathbf{I}_{p})$. Since
$\frac{\mathrm{E}(w_{\ast})}{\bar{w}}\approx 1$, we then have

$n\bar{w}\left(\sum_{i=1}^{n}w_{i}\mathbf{z}_{i}\mathbf{z}_{i}^{\top}\right)^{-1}\stackrel{{\scriptstyle\rm
d}}{{\approx}}n_{\tau}\mathbf{V}^{-1}.$

Assume $n_{\tau}>p+1$. Note that $\mathbf{V}^{-1}$ follows an inverse-Wishart
distribution with degrees of freedom $n_{\tau}$ and scale matrix
$\mathbf{I}_{p}$. It then follows that $\mathrm{E}[{\rm
tr}(\mathbf{V}^{-1})]=\frac{p}{n_{\tau}-p-1}$ (Mardia et al.,, 1979). Using
$\frac{n_{\tau}p}{n_{\tau}-p-1}$ as an approximation of $n\bar{w}{\rm
tr}[(\mathbf{X}^{\top}\mathbf{W}\mathbf{X})^{-1}\bm{\Sigma}_{p})]$, we then
have

$\mathrm{df}_{\rm
R}(p)\approx\frac{p}{2}\zeta_{\tau}\left(1+\frac{n_{\tau}}{n_{\tau}-p-1}\right).$

∎

###  A.6 Proof of Theorem 4.3

###### Proof.

We first assume that $\mathbf{Q}$ is an arbitrary $n\times n$ matrix such that
$(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}$ exists. Define
$\bm{\Pi}\colon\mathbb{R}^{n\times n}\to\mathbb{R}^{p\times n}$ as

$\bm{\Pi}(\mathbf{Q})=(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}.$

Then $\mathrm{df}_{\rm R}$ can be written as

$\mathrm{df}_{\rm R}=\zeta_{\tau}p+\frac{n}{2}{\rm
tr}[\bm{\Pi}(\mathbf{Q})\mathbf{T}\bm{\Pi}(\mathbf{Q})^{\top}\mathbf{R}].$

Note that

$\displaystyle\partial\bm{\Pi}(\mathbf{Q})$
$\displaystyle=-(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}(\partial\mathbf{Q})\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}+(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}(\partial\mathbf{Q})$
$\displaystyle=(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}(\partial\mathbf{Q})(\mathbf{I}_{n}-\mathbf{H}).$

Then we have

$\displaystyle\partial(\bm{\Pi}(\mathbf{Q})\mathbf{T}\bm{\Pi}(\mathbf{Q})^{\top}\mathbf{R})$
$\displaystyle=(\partial\bm{\Pi}(\mathbf{Q}))\mathbf{T}\bm{\Pi}(\mathbf{Q})^{\top}\mathbf{R}+\bm{\Pi}(\mathbf{Q})\mathbf{T}(\partial\bm{\Pi}(\mathbf{Q}))^{\top}\mathbf{R}$
$\displaystyle=(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}(\partial\mathbf{Q})(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}\mathbf{Q}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{R}$
$\displaystyle\qquad+(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}\mathbf{T}(\mathbf{I}_{n}-\mathbf{H})^{\top}(\partial\mathbf{Q})\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{R}.$

By the chain rule, it follows that

$\displaystyle\partial\mathrm{df}_{\rm R}$
$\displaystyle=\frac{n}{2}\partial{\rm
tr}[\bm{\Pi}(\mathbf{Q})\mathbf{T}\bm{\Pi}(\mathbf{Q})^{\top}\mathbf{R}]$
$\displaystyle=\frac{n}{2}{\rm
tr}[\partial(\bm{\Pi}(\mathbf{Q})\mathbf{T}\bm{\Pi}(\mathbf{Q})^{\top}\mathbf{R})]$
$\displaystyle=\frac{n}{2}{\rm
tr}[(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}(\partial\mathbf{Q})(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}\mathbf{Q}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{R}]$
$\displaystyle\qquad+\frac{n}{2}{\rm
tr}[(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{Q}\mathbf{T}(\mathbf{I}_{n}-\mathbf{H})^{\top}(\partial\mathbf{Q})\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{R}]$
$\displaystyle=n{\rm
tr}[(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}\mathbf{Q}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{R}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\partial\mathbf{Q}].$

Thus,

$\frac{\partial\mathrm{df}_{\rm
R}}{\partial\mathbf{Q}}=n(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}\mathbf{Q}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{R}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}.$

Since $\mathbf{Q}$ is a diagonal matrix, we further have

$\frac{\partial\mathrm{df}_{\rm
R}}{\partial\mathbf{q}}=n\,\mathrm{diag}\left((\mathbf{I}_{n}-\mathbf{H})\mathbf{T}\mathbf{Q}\mathbf{X}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{R}(\mathbf{X}^{\top}\mathbf{Q}\mathbf{X})^{-1}\mathbf{X}^{\top}\right),$

where $\mathbf{q}=\mathrm{diag}(\mathbf{Q})$. ∎

###  A.7 Proof of Theorem 5.1

###### Proof.

Let $\mathbf{H}$ and
$\hat{\bm{\varepsilon}}=(\hat{\varepsilon}_{1},\ldots,\hat{\varepsilon}_{n})^{\top}$
be the hat matrix and residual vector of the full model respectively. Write
$\hat{\bm{\varepsilon}}_{\rm
loocv}=(\hat{\varepsilon}_{1}^{-1},\ldots,\hat{\varepsilon}_{n}^{-n})^{\top}$.
Note that $\hat{\bm{\varepsilon}}=(\mathbf{I}_{n}-\mathbf{H})\mathbf{y}$,
$\hat{\bm{\varepsilon}}_{\rm
loocv}=\mathrm{diag}((1-h_{ii})^{-1})(\mathbf{I}_{n}-\mathbf{H})\mathbf{y}$,
and $\mathrm{Var}(\mathbf{y}|\mathbf{X})=\sigma_{\varepsilon}^{2}\mathbf{T}$.
Then, we have

$\mathrm{Var}(\mathbf{W}^{\frac{1}{2}}\hat{\bm{\varepsilon}}|\mathbf{X})=\sigma_{\varepsilon}^{2}\mathbf{W}^{\frac{1}{2}}(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}(\mathbf{I}_{n}-\mathbf{H})^{\top}\mathbf{W}^{\frac{1}{2}},$

and

$\mathrm{Var}(\mathbf{W}^{\frac{1}{2}}\hat{\bm{\varepsilon}}_{\rm
loocv}|\mathbf{X})=\sigma_{\varepsilon}^{2}\mathrm{diag}\left(\frac{\sqrt{w_{i}}}{1-h_{ii}}\right)(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}(\mathbf{I}_{n}-\mathbf{H})^{\top}\mathrm{diag}\left(\frac{\sqrt{w_{i}}}{1-h_{ii}}\right).$

Thus,

$\displaystyle{\rm tr}(\mathbf{A}\mathbf{T})$ $\displaystyle={\rm
tr}[(\mathbf{I}_{n}-\mathbf{H})^{\top}\mathrm{diag}(w_{i}(1-h_{ii})^{-2})(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}]-{\rm
tr}[(\mathbf{I}_{n}-\mathbf{H})^{\top}\mathrm{diag}(w_{i})(\mathbf{I}_{n}-\mathbf{H})\mathbf{T}]$
$\displaystyle=\frac{1}{\sigma_{\varepsilon}^{2}}({\rm
tr}[\mathrm{Var}(\mathbf{W}^{\frac{1}{2}}\hat{\bm{\varepsilon}}_{\rm
loocv}|\mathbf{X})]-{\rm
tr}[\mathrm{Var}(\mathbf{W}^{\frac{1}{2}}\hat{\bm{\varepsilon}}|\mathbf{X})])$
$\displaystyle=\frac{1}{\sigma_{\varepsilon}^{2}}\sum_{i=1}^{n}w_{i}[\mathrm{Var}(\hat{\varepsilon}^{-i}_{i}|\mathbf{X})-\mathrm{Var}(\hat{\varepsilon}_{i}|\mathbf{X})].$

On the other hand, since $\mathrm{df}_{\rm R}$ is defined through the excess
variance of the predictions on the test data against those on the training
data, it can be rewritten as

$\displaystyle\mathrm{df}_{\rm R}$
$\displaystyle=\frac{1}{2\bar{w}\sigma_{\varepsilon}^{2}}\left(n\mathrm{E}[w_{\ast}\mathrm{Var}(\hat{\mu}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})|\mathbf{X}]-\sum_{i=1}^{n}w_{i}[\mathrm{Var}(\hat{\mu}_{i}|\mathbf{X})-2\mathrm{Cov}(y_{i},\hat{\mu}_{i}|\mathbf{X})]\right)$
$\displaystyle=\frac{1}{2\bar{w}\sigma_{\varepsilon}^{2}}\left(n\mathrm{E}[w_{\ast}\mathrm{Var}(\hat{\varepsilon}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})|\mathbf{X}]-\sum_{i=1}^{n}w_{i}\mathrm{Var}(\hat{\varepsilon}_{i}|\mathbf{X})-n\sigma_{\varepsilon}^{2}\mathrm{E}(w_{\ast}\tau_{\ast})+\sigma_{\varepsilon}^{2}\sum_{i=1}^{n}w_{i}\tau_{i}\right).$

Therefore, we have

$\frac{1}{n}\sigma_{\varepsilon}^{2}\xi_{\mathbf{X}}=\mathrm{E}[w_{\ast}\mathrm{Var}(\hat{\varepsilon}_{\ast}|\mathbf{x}_{\ast},\mathbf{X})|\mathbf{X}]-\frac{1}{n}\sum_{i=1}^{n}w_{i}\mathrm{Var}(\hat{\varepsilon}^{-i}_{i}|\mathbf{X})-\sigma_{\varepsilon}^{2}\left(\mathrm{E}(w_{\ast}\tau_{\ast})-\frac{1}{n}\sum_{i=1}^{n}w_{i}\tau_{i}\right).$

∎

##  Appendix B Data Analysis

###  B.1 Cancer Mortality Data

Table 2 shows the variables and their estimated coefficients of the models
selected by $\widehat{\mathrm{wErrR}}$ for both the OLS and WLS methods based
on a particular partition of the data into the training and test sets. In this
example, both methods select the first 10 variables, including the categorical
variable Region that consists of 3 categories: Midwest, Northeast and South.
We see that the estimated coefficients for these 3 regions have relatively
large change when case weights are incorporated into the model. Since the
weights are estimated based on the population size, and county population
varies a lot across different regions, the change in these coefficients is a
direct result of the weight scheme. On the other hand, the contrasts of these
coefficients as well as the coefficients of other numerical variables remain
relatively stable from the OLS to the WLS model, suggesting that the weight
scheme affects more on the variance structure of the model rather than the
regression mean.

Table 2: Estimated coefficients in the optimal models identified by
$\widehat{\mathrm{wErrR}}$ for the OLS and WLS methods.

Variable |  Description |  Transformation |  OLS |  WLS
---|---|---|---|---
RegionMidwest |  Whether the county is in the Midwest |  - |  $43.222$ * |  $49.778$ *
RegionNortheast |  Whether the county is in the Northeast |  - |  $32.278$ |  $38.988$ *
RegionSouth |  Whether the county is in the South |  - |  $51.235$ * |  $57.468$ *
avgDeathRateEst2015 |  Estimated death rate based on 2015 population estimates and average number of reported cancer mortalities from 2010 to 2016 |  - |  $0.243$ *** |  $0.237$ ***
incidenceRate |  Mean per capita (100,000) cancer diagnoses |  - |  $0.209$ *** |  $0.195$ ***
PctPublicCoverage |  Percent of county residents with government provided health coverage |  logit |  $-41.652$ *** |  $-38.457$ ***
PctBachDeg25_Over |  Percent of county residents aged 25 and over with bachelor’s degree as the highest education attained |  logit |  $-9.127$ * |  $-10.532$ **
PctMarried |  Percent of county residents who are married |  logit |  $-9.351$ |  $-11.408$ $\cdot$
PctPrivateCoverage |  Percent of county residents with private health coverage |  logit |  $-27.979$ *** |  $-23.903$ ***
AvgHouseholdSize |  Average number of people in a household |  - |  $-0.148$ |  $-0.833$
AvgAnnCount |  Average number of reported cases of cancer diagnosed annually |  log |  $-1.078$ |  $-0.851$
PctUnemployed16_over |  Unemployment rate for county residents aged 16 and over |  logit |  $11.463$ *** |  $12.204$ ***

  * •

Significance codes: ‘***’ $<0.001$, ‘**’ $0.001\sim 0.01$, ‘*’ $0.01\sim
0.05$, ‘.’ $0.05\sim 0.1$, ‘ ’ $>0.1$.

###  B.2 King County House Sales Data

In this section, we provide more details on the King County housing data
analysis described in Section 6.2 of the paper. Table 3 lists all 12 variables
considered in the regression analysis. To simplify the analysis, some of these
variables are redefined based on variables in the original dataset, including

  * •

bedrooms: We convert the original integer-valued variable to a categorical
variable of 3 levels.

  * •

condition: We map the original condition indexed from 1 (worst) to 5 (best) to
a categorical variable of 3 levels: poor (1 and 2), average (3) and good (4
and 5).

  * •

yr_built: The variable is initially the year when the house was built. We
group its values into 6 time periods.

  * •

zone: This variable summarizes zipcode (the zip code of the area in which the
house is located) in the original dataset into 3 zones. To determine the
grouping, we fit an OLS model with the first 11 variables in Table 3 using all
observations in 2014. The residuals are then grouped by the zip code and the
median residuals for each zip code are studied through a boxplot shown in the
left panel of Figure 15. Two outlying zip codes with much higher median
residuals are identified. We define these two zip codes as zone 1 (98039) and
zone 2 (98004) respectively. The remaining zip codes are merged into one
common zone called zone 3.

Table 3: List of independent variables used in King County house sales data analysis Variable | Type |  Description | Transformation
---|---|---|---
basement | categorical |  Whether the house has basement or not: yes/no | -
bathrooms | numeric |  Number of bathrooms in the house | -
bedrooms | categorical |  Number of bedrooms in the house: 1-2/3-5/6 or more | -
condition | categorical |  Condition of the house: poor/average/good |
grade | numeric |  Overall grade of the house regarding the types of materials used and the quality of workmanship | -
grade2 | numeric |  Squared of grade | -
sqft_living | numeric |  Interior living space of the houses in square feet | square root
sqft_living15 | numeric |  Average interior living space for the closest 15 houses in square feet | square root
view | numeric |  Grade of the view around the house | -
waterfront | categorical |  Whether the house has a waterfront or not: yes/no | -
yr_built | categorical |  Year when the house was built: 1920 or before/1921-1940/1941-1960/1961-1980/1981-2000/2001 or after | -
zone | categorical |  Location of the house: zone 1/zone 2/zone 3 | -

Figure 14: Boxplot of the median residuals of the OLS model for each zip code
(left) and the map of the three zones defined through the boxplot (right).

Figure 15: Comparison of the selected model size $\hat{p}_{\ast}$ between the
OLS model and WLS model based on 500 random subsets of size $n=200$ from the
King County house sales data in 2015. The upper panels compare the model size
directly whereas the lower panels show the difference in model size.

To estimate the variance function $\tau$, we first fit an OLS model with all
variables included using the samples in 2014. The residuals are then bucketed
by 100 percentiles of the fitted values, and their variances are calculated
and plotted against the median fitted values for each bucket, as shown in
Figure 10 in the paper. It turns out that a quadratic function of the fitted
value fits the error variances well. This is approximately equivalent to
assuming that the variance function $\tau$ is of the form

$\tau(\mathbf{x})=a(1+b|\mu|)^{2},\quad a>0,\,b\neq 0.$

With the estimated variance function, the variable sequence is then determined
using a similar approach described in Section 6.2.1 based on all samples in
2014.

Figure 15 shows the size of the selected model for the two procedures based on
500 random training subsets, whereas Figure 16 presents the decompositions of
$\widehat{\mathrm{wErrF}}$ and $\widehat{\mathrm{wErrR}}$ into the training
error, excess bias and excess variance for the same training subsets.

Figure 16: Median percentages of the training error, excess bias and excess
variance in the prediction error estimates $\widehat{\mathrm{wErrF}}$ and
$\widehat{\mathrm{wErrR}}$ of the selected model as the model size $p$ varies.
The shaded bands represent the 2.5th and 97.5th sample percentiles for each
$p$.