# Unbiased Estimation of Structured Prediction Error

Kevin Fry  Department of Statistics
Stanford University
Stanford, CA
<EMAIL_ADDRESS>Jonathan E. Taylor  Department of Statistics
Stanford University
Stanford, CA
<EMAIL_ADDRESS>
###### Abstract

Many modern datasets, such as those in ecology and geology, are composed of
samples with spatial structure and dependence. With such data violating the
usual independent and identically distributed (IID) assumption in machine
learning and classical statistics, it is unclear a priori how one should
measure the performance and generalization of models. Several authors have
empirically investigated cross-validation (CV) methods in this setting,
reaching mixed conclusions. We provide a class of unbiased estimation methods
for general quadratic errors, correlated Gaussian response, and arbitrary
prediction function $g$, for a noise-elevated version of the error. Our
approach generalizes the coupled bootstrap (CB) of Oliveira et al., [19] from
the normal means problem to general normal data, allowing correlation both
within and between the training and test sets. CB relies on creating bootstrap
samples that are intelligently decoupled, in the sense of being statistically
independent. Specifically, the key to CB lies in generating two independent
“views” of our data and using them as stand-ins for the usual independent
training and test samples. Beginning with Mallows’ $C_{p}$ [15], we generalize
the estimator to develop our generalized $C_{p}$ estimators (GC). We show at
under only a moment condition on $g$, this noise-elevated error estimate
converges smoothly to the noiseless error estimate. We show that when Stein’s
unbiased risk estimator (SURE) applies, GC converges to SURE as in the normal
means problem. Further, we use these same tools to analyze CV and provide some
theoretical analysis to help understand when CV will provide good estimates of
error. Simulations align with our theoretical results, demonstrating the
effectiveness of GC and illustrating the behavior of CV methods. Lastly, we
apply our estimator to a model selection task on geothermal data in Nevada.

##  1 Background

Prediction error estimation is of central importance to modern statistics,
machine learning, and artifical intelligence. Understanding how well a model
will predict on unseen data is crucial to the determination of whether it
should be used, how it should be used, and to what extent humans can trust and
rely on its predictions. In the standard setting, it is assumed that data
points are IID or exchangeable, which is amenable to standard practices of
data splitting and CV. However, in settings where the data is not IID, the
answer becomes more complicated. Methods must account for the correlation
between sample points that bias usual error estimates. However, in practice we
see practitioners still employing CV for the purpose of prediction error
estimation of spatial models [22, 31].

###  1.1 Problems with CV

We begin by reviewing CV methods. Recall that a CV method uses some procedure
to split the dataset into $k$ folds, repeatedly trains the model on $k-1$
folds and gets an estimate of performance on the held out fold. The estimates
on the held-out folds are averaged to produce the final estimate of
performance. The most basic version of this is standard $k$-fold CV (KFCV),
whereby observations are randomly assigned to each fold, only requiring each
fold be the same size. However, this IID assumption breaks down in spatial
settings. Aiming to remedy this, spatial $k$-fold CV (SPCV) splits the data
into $k$ folds while keeping nearby points together, e.g. all observations
from the same county are put in the same fold. The idea here is to cut down
the correlation between folds by keeping spatially nearby points together.
Going further with this idea, one can consider a buffered-leave-one-out CV
(BLOOCV) approach, whereby each observation is its own fold, but we
additionally remove observations within a spatial buffer of the validation
observation from training.

However, these spatial CV approaches cannot be blindly applied to spatial data
problems. One may think this is due only to the varying correlation structure
between differing training locations and testing locations. This is not the
case. Even when the training and testing locations are the same CV has issues.
An example of a such a situation would be measuring temperature at the same
place at different times, while using time-invariant features for prediction.
That is we have a $X,Y$ and $Y^{*}$ with samples indexed to the same
locations. This sampling scheme, and others, are encapsulated by Section 1.2.
We demonstrate even here CV methods struggle with an ordinary least squares
(OLS) example with $n=100$ and $p=5$. We generate $Y^{*},Y$ according to the
two mechanisms in Section 1.2. In the first, they are independent draws
conditional on $X$ from $\mathcal{N}(X\beta,\Sigma_{Y})$, referred to as the
no shared noise model (NSN). In the second, they are correlated conditional on
$X$ and they are drawn according to Eq. 1 with non-zero $\Sigma_{Y,Y^{*}}$.
This is referred to as the shared structural noise model (SSN). The error
metric is the relative mean squared error, with (Monte-Carlo estimated) ground
truth (red) as the baseline. Full simulation details are given in Section 6.

Figure 1: Relative MSE of adaptive linear regression for independent
$Y^{*},Y$. ‘GenCp’ refers to our proposed method, specifically the variant
given in Eq. 4. ‘KFCV’ is for standard $k$-fold CV and ‘SPCV’ is spatial CV
with $k$-means learned groups. Full details given in Section 6.

It is clear both regular $k$-fold CV and spatial CV are biased estimators of
the prediction error. We primarily focus on the NSN and SSN settings, but do
discuss CV for general settings, including out-of-sample error estimation, in
Section 5.

###  1.2 Setting and notation

We consider this problem in the setting of normally distributed data with
arbitrary correlation where $X\in\mathbb{R}^{n\times p}$ is a common feature
matrix that is considered fixed throughout and shared between
$Y^{*},Y\in\mathbb{R}^{n}$. Explicitly,

$\begin{pmatrix}Y^{*}\\\ Y\end{pmatrix}\mid
X\sim\mathcal{N}\left(\mu,\Sigma\right),\qquad\Sigma=\begin{pmatrix}\Sigma_{Y^{*}}&\Sigma_{Y^{*},Y}\\\
\Sigma_{Y,Y^{*}}&\Sigma_{Y}\end{pmatrix}$ (1)

And for a prediction function $g$ fit on $Y$ we want to estimate the
generalized quadratic loss

$\mathbb{E}\|Y^{*}-g(Y)\|_{\Theta_{p}}^{2}=\mathbb{E}[(Y^{*}-g(Y))^{\top}\Theta_{p}(Y^{*}-g(Y))]$

Under what setting would we arrive at such data? A useful class of data-
generating mechanisms is the following. Consider a process $Y$, such that

$\displaystyle Y_{t}(X_{t})$ $\displaystyle=f(X_{t})+\epsilon_{t}$
$\displaystyle\epsilon_{t}$ $\displaystyle=\epsilon^{S}_{t}+\epsilon^{M}_{t}$

where $t\in\mathbb{R}^{d}$ is some spatial index, for instance $t$ could be
time, locations on Earth, or regions of the brain. The superscripts $S$ and
$M$ denote spatial noise (e.g. due to unmodeled features) and measurement
noise, respectively. We then consider a replicate $Y^{*}$, which could take on
either of the following forms

$\displaystyle Y^{*}_{t}(X_{t})$ $\displaystyle=f(X_{t})+\epsilon^{*}_{t}$
$\displaystyle\epsilon^{*}_{t}$
$\displaystyle=\epsilon^{S,*}_{t}+\epsilon^{M,*}_{t}$ (NSN) (2)
$\displaystyle\epsilon^{*}_{t}$
$\displaystyle=\epsilon^{S}_{t}+\epsilon^{M,*}_{t}$ (SSN)

Here we distinguish between the two models by referring to the former as the
no shared noise (NSN) model, and the latter as the shared structural noise
model (SSN) model. This is to emphasize NSN generates independent vectors $Y$
and $Y^{*}$ given $X$, while SSN generates vectors $Y,Y^{*}$ that share the
structural component of their covariance matrices. Note here that we do not
require $\epsilon_{t},\epsilon^{*}_{t}$, or even their spatial and measurement
components to be independent or identically distributed. Now for some sample
of $n$ locations and model $g$, we want to estimate the prediction error,
which for this paper we define as

$\mathrm{Err}(g)=\mathbb{E}\|Y^{*}-g(Y)\|_{\Theta_{p}}^{2}$

In words, we are interested in the MSE for a new replicate $Y^{*}$ with the
same features $X$ at the same locations we have observations of $Y$ for. This
can be viewed as the spatial analog to the in-sample error setting considered
by Mallows’ $C_{p}$ and other covariance penalty estimators [1, 6, 7, 11, 15].
The generalization being that we allow $Y^{*}$ to be correlated with $Y$,
which could arise due to unmodeled spatial variation, for instance.

###  1.3 Stein’s unbiased risk estimator

One cannot write about unbiased risk estimation in normal data without mention
of the seminal work of Charles Stein [26]. The most well-known and widely
applied version of the result provides an unbiased estimator of squared error
for IID normal data $Y_{i}\sim\mathcal{N}(\mu,\sigma^{2})$, and weakly
differentiable estimator $g$. This has been used to provide unbiased estimates
of risk and prove the degrees of freedom of notable prediction functions, such
as the LASSO in [30]. A version of SURE exists for general Gaussian data, but
unbiasedly estimates a generalized quadratic form matched to the covariance
structure of the data. This whitens the data and in effect reduces it to the
isotropic Gaussian case.

######  Theorem 1.1.

Let $Y\sim\mathcal{N}(\mu,\Sigma)$ with $\Theta=\Sigma^{-1}$. Let estimator
$g:\mathbb{R}^{n}\mapsto\mathbb{R}$ be s.t.
$f(z)=\Theta^{1/2}g(\Theta^{-1/2}z)$ is weakly differentiable where
$\Theta^{1/2}$ is the symmetric square root. Write $\nabla\cdot
g:=\sum_{i}\nabla_{i}g_{i}$ for the divergence of $g$ with $\nabla_{i}g_{j}$
being the weak partial derivative of the $j$-th component of $g$ with respect
to $y_{i}$. Assume that $\mathbb{E}\|g(Y)\|_{\Theta}^{2}<\infty$ and
$\mathbb{E}|\nabla_{i}g_{i}(Y)|<\infty$. Then,

$\mathrm{SURE}(g,\Theta)=\|Y-g(Y)\|_{\Theta}^{2}+2(\nabla\cdot g)(Y)$ (3)

is an unbiased estimator of $\mathrm{Err}(g)$.

This is helpful in many cases but again requires weak differentiability and
calculating a divergence term. Our proposed estimator does not make any
assumptions on $g$ and works with any generalized quadratic
$\|\cdot\|_{Q}^{2}$ for any $Q$. However, in Section 3 we show that our
estimator does converge to SURE when SURE is well defined.

SURE has been extended to discontinuous functions by [29]. However, this
extension relies on calculating the discontinuity sets of the function, which
can be very difficulat. Despite this, Tibshirani, [29] use a different
approach to derives degrees of freedom for best subsets and relaxed lasso for
IID normal response and orthogonal $X$. More recently, Oliveira et al., [19]
break from the SURE approach, proposing the coupled bootstrap (CB) to estimate
error for general $g$ and IID response.

###  1.4 Contributions

This paper has two main contributions. First, the most general form of
generalized $C_{p}$ provides an unbiased estimate for generalized quadratic
prediction error when the joint vector $(Y^{*},Y)$ has arbitrary covariance
structure $\Sigma$ and arbitrary prediction function $g$ fitted on noise-
elevated data. We present several GC variants depending on the specific
setting, but refer to them interchangeably as generalized $C_{p}$ or GC
estimators since they fall into the same general framework and can be seen as
generalizations of Mallows’ $C_{p}$. Further, we prove the noise-elevated
error smoothly converges to the noiseless error in Proposition 3.1. Empirical
results verify GC estimators with small added noise do indeed closely
approximate the noiseless error. We also show in Theorem 3.1 that when SURE
exists, GC converges to SURE. In addition, we develop the special case of
estimation of prediction error of bagged arbitrary prediction function $g$.
Second, we apply a similar analysis to CV methods to provide a theorem
describing when CV would provide consistent error estimates. This theoretical
result not only adds theoretical context explaining empirical results from
[22, 31], but also provides guidelines for how practitioners should choose a
CV method for their dataset and application.

The structure of the rest of the paper is as follows:

  * •

In Section 2 we develop GC methods for error (and therefore risk) estimation
in general Gaussian data. We begin with adaptive linear smoothers, build into
bagging these smoothers, and ultimately to general (bagged) models. We note
that GC even allows for correlation between the observed $Y$ and unobserved
$Y^{*}$, as well as $Y^{*}$ and $Y$ having different covariance structure. In
Section 2.4 we extend these results to bagged estimators, and discuss
estimating the covariance $\Sigma$ in Section 2.6.

  * •

In Section 3 we present theoretical results of GC estimators. We show that the
noise-elevated error converges to the noiseless error as we reduce the added
noise to 0. We also demonstrate GC recovers SURE int the limit when SURE
exists.

  * •

In Section 4 we extend the usual notion of optimism to our setting, and
analyze the natural nonparamteric bootstrap estimators that arise from it,
showing their insufficiency for the task.

  * •

In Section 5 we use our results to formally investigate what CV methods are
estimating in the spatial setting, and how that differs from what one may
think it is estimating. We use this insight to attempt to resolve the debate
between [22] and [31], and give advice on how practitioners should choose a CV
method for their particular application.

  * •

In Section 6 we present simulations validating the unbiasedness of our
estimator and comparing it to CV methods. In Section 7 we use GC for model
selection on geothermal data in Nevada. We also demonstrate that in practical
settings, CV estimators deviate from GC.

  * •

In Section 8 we conclude, discussing future steps such as predicting local
maxima, incorporating the approach in CV methods, and extending to other
distributions.

###  1.5 Related work

Estimating and understanding model performance is a classical problem in
statistics and machine learning. In the above exposition, we have shown that
CV struggles to accurately estimate error in the SSN setting as well as the
easier NSN setting. Further we have established a notion of error that is an
extension of in-sample error for the spatial setting. For IID data and in-
sample error, covariance penalty approaches such as the classical Mallows’
$C_{p}$ and AIC [1, 15], with extensions to the ideas of optimism and degrees
of freedom [6, 7, 11] have given sufficient answers to estimating prediction
error of a model. We frame our work as building off of Mallows’, moving past
the IID assumption to spatially correlated data.

In the standard IID case, the foundational result of [15] offers a
satisfactory answer for linear models, needing only an estimate of the noise
level $\sigma$. We develop our approach with Mallows’ $C_{p}$ as a starting
point, extending it to the spatial settings of NSN and SSN defined in Section
1.2 for adaptive linear smoothers, which gives general predictive models as a
corollary. Our approach can be seen as an elaboration on the use of the
parametric bootstrap to estimate prediction error. There have been proposed
methods using a parametric bootstrap to estimate risk in IID normal data [3,
7, 32]. However, these methods do not capture the full optimism, see Oliveira
et al., [19] for a succinct and clear discussion of this. We also show the
spatial extenstions of these parametric bootstrap methods are also not fit-to-
task in Section 4. Tian, [27] estimates prediction error for linear models
after feature selection. One can view this as a special case of Theorem 2.3,
for IID Gaussian data and a subset of adaptive linear smoothers. Finally, our
method addresses bagged estimators which previous works do not. Oliveira et
al., [19] generalizes Tian, [27] to arbitrary model choice $g$, but again only
in the IID setting. We discuss these works in more depth in Section 2.7. Both
[19] and [27], as well as earlier works [23, 28], leverage a trick for
Gaussian vectors to generate two independent views with the same mean as the
original vector. This has since then been named “data fission” and generalized
to other distributions by Leiner et al., [13], and the interested reader
should refer to their paper for an extensive enumeration of “fissions”. In the
present work, we leverage this popular Gaussian “fission” trick as well. The
authors of [19] extend their results to the IID Poisson case, where an
analogous trick applies, in [20]. We give a more in-depth discussion of [19,
27] in Section 2.7.

Other authors have looked at other error metrics for CV methods and other
covariance penalties [2, 24]. Bates et al., [2] investigate CV in the IID
setting, and find that the method does not estimate out-of-sample error, that
is the expected error over a new draw $(x_{0},y_{0})$ for a model trained on
the data at hand, but instead more closely estimates the average out-of-sample
error over training sets $(X,Y)$. Although the majority of our paper is for
the fixed-$X$ setting, in Section 5 we derive insights for how CV performs
both in-sample and out-of-sample in the spatial setting. These insights
provide heuristics for how one should construct their CV procedure depending
on the particulars of the dataset and problem. Rosset and Tibshirani, [24]
takes a random-$X$ approach to in-sample error estimation, presenting a
generalization to Mallows’ $C_{p}$ they term $RC_{p}$. Although we also start
with Mallows’ $C_{p}$, our approach keeps to the fixed-$X$ regime and instead
moves beyond the IID setting.

As mentioned in the introduction, CV methods have been used by practitioners
to estimate prediction error in spatial data. Parmentier et al., [21] propose
the BLOOCV approach, where observations within a certain buffer radius of the
held-out observation are not used in training. With a sufficiently large
buffer, one would expect the training set to be independent of the held out
observation, recovering the nice properties of the IID setting. However, if
the future data points will be correlated with the training data, then BLOOCV
will obviously overestimate error. As discussed above in Section 1.1, CV
approaches do not estimate prediction error well in the SSN setting. Ploton et
al., [22] applied CV methods, including BLOOCV, to spatial data, arguing
spatial CV methods are necessary to not overestimate model performance. Soon
after, Wadoux et al., [31] responded to this, disagreeing and apparently
demonstrating standard $k$-fold CV was a better estimate of performance than
its spatial counterparts. Both of these papers are concerned with out-of-
sample error. After presenting GC, we use the same tools to analyze CV
estimators and derive a result explaining when CV should provide good
estimates of out-of-sample error for correlated data in Section 5.

##  2 Methodology

###  2.1 Mallows’ $C_{p}$

Mallows’ $C_{p}$ [15] is a classical result for getting an unbiased estimate
of prediction error in linear regression, also known as ordinary least
squares, for IID data. Formally, if

$(Y^{*},Y)\mid X\overset{iid}{\sim}\mathcal{N}(\mu,\sigma^{2}I)$

and $P=X(X^{\top}X)^{-1}X^{\top}$ is the usual OLS estimate, then

$\displaystyle\mathbb{E}\|Y^{*}-PY\|_{2}^{2}$
$\displaystyle=\mathbb{E}\|Y-PY\|_{2}^{2}+2p\sigma^{2}$

Although a useful result on its own, it falls short for our setting. It
requires IID data, only gives you predictions for the same sites you use for
fitting the model, and it does not hold for general models.

It is easy enough to relax the IID assumption. If instead, $(Y^{*},Y)\mid
X\overset{ind}{\sim}\mathcal{N}(\mu,\Sigma_{Y})$, one can easily see the
Mallows’ unbiased estimator is then:

$\mathbb{E}\|Y^{*}-PY\|_{2}^{2}=\mathbb{E}\|Y-PY\|_{2}^{2}+\mathrm{tr}(2P\Sigma_{Y})$

Now if we consider a general linear smoother $S$, we can arrive at an
estimator of the same form:

$\displaystyle\mathbb{E}\|Y^{*}-SY\|_{2}^{2}$
$\displaystyle=\mathbb{E}\|Y^{*}-Y\|_{2}^{2}+\mathbb{E}\|Y-SY\|_{2}^{2}+2\mathbb{E}\left[(Y^{*}-Y)^{\top}(Y-SY)\right]$
$\displaystyle=2\mathrm{tr}(\Sigma_{Y})+\mathbb{E}\|Y-SY\|_{2}^{2}-2\mathbb{E}\left[(Y-\mu)^{\top}(I-S)(Y-\mu)\right]$
$\displaystyle=\mathbb{E}\|Y-SY\|_{2}^{2}+2\mathrm{tr}(S\Sigma_{Y})$

Now we can add estimation and prediction matrices $\Theta_{e},\Theta_{p}$
respectively. Here $\Theta_{e}$ is some augmentation matrix applied to the
data before feeding it into the model. $\Theta_{p}$ is used to generalize our
error matrix with a generalized quadratic $\|\cdot\|_{\Theta_{p}}$. A useful
example is when $\Theta_{e}$ is the selector matrix for a training set
$\mathcal{T}$. For example, if $\mathcal{T}=\\{1,3\\}$ and $n=4$, then

$\Theta_{e}=\begin{bmatrix}1&0&0&0\\\ 0&0&1&0\end{bmatrix}$

Then setting $\Theta_{p}=I-\Theta_{e}^{\top}\Theta_{e}$ we derive an estimator
for a train/test split of our data:

$\displaystyle\mathbb{E}\|Y^{*}-S_{e}Y\|_{\Theta_{p}}^{2}$
$\displaystyle=\mathbb{E}\|Y-S_{e}Y\|_{\Theta_{p}}^{2}+2\mathrm{tr}(\Theta_{p}S_{e}\Sigma_{Y})$
(4)

where $S_{e}$ denotes the smoother based on $\Theta_{e}X,\Theta_{e}Y$. For the
rest of the paper, we often suppress the subscript $e$ notation as the results
do not change for $S$ fit on the full training data or on the training subset.

Above, we have constructed a Mallows’ $C_{p}$-like estimator for structured
data test set error for linear smoothing models. This is of some interest on
its own, but still leaves out a vast range of models used by practitioners
today.

###  2.2 Data “fission” by randomization

To generalize beyond non-adaptive linear smoothers, that is smoothers that are
not dependent on $Y$, we need to introduce a simple but clever tool.

A well known trick allows you to “fission” a Gaussian vector
$Y\sim\mathcal{N}(\mu,\Sigma_{Y})$ to generate two independent views with the
same mean $\mu$. This approach has been used for various applications such as
selective inference [28] and nonparametric confidence sets [23]. This method
is expanded upon in [13], from which we borrow the term data fission we will
use throughout. By cleverly adding noise $\omega$:

$\begin{gathered}\omega\sim\mathcal{N}(\mathbf{0},\Sigma_{Y}),\qquad
W=Y+\sqrt{\alpha}\omega\sim\mathcal{N}(\mu,\Sigma_{W}),\qquad
W^{\perp}=Y-\omega/\sqrt{\alpha}\sim\mathcal{N}(\mu,\Sigma_{W^{\perp}})\\\
\Sigma_{W}=(1+\alpha)\Sigma_{Y},\qquad\Sigma_{W^{\perp}}=(1+1/\alpha)\Sigma_{Y}\end{gathered}$
(5)

We note that one cannot recover $Y$ from either view alone, but with both we
have the decomposition

$\begin{gathered}Y=AW+A^{\perp}W^{\perp}\\\
A=I-A^{\perp},A^{\perp}=(1+1/\alpha)^{-1}I\end{gathered}$ (6)

With this, we can derive a Mallows’-like unbiased estimate of
$\mathrm{Err}(S)=\mathbb{E}\|Y^{*}-S(\Theta_{e},W)Y\|_{\Theta_{p}}^{2}$ for a
data-dependent linear smoother $S(\Theta_{e},\cdot)$:

######  Theorem 2.1.

With the above setup, an unbiased estimate of the prediction error
$\mathbb{E}\|Y^{*}-S(W)Y\|_{\Theta_{p}}^{2}$ is given by

$\|W^{\perp}-S(W)Y\|_{\Theta_{p}}^{2}+\mathrm{tr}(\Theta_{p}\Sigma_{Y})-\mathrm{tr}(\Theta_{p}\Sigma_{W^{\perp}})+2\mathrm{tr}(\Theta_{p}S(W)\Sigma_{Y})$
(7)

and an unbiased estimate of $\mathbb{E}\|Y^{*}-S(W)W\|_{\Theta_{p}}^{2}$ is
given by

$\|W^{\perp}-S(W)W\|_{\Theta_{p}}^{2}+\mathrm{tr}(\Theta_{p}\Sigma_{Y})-\mathrm{tr}(\Theta_{p}\Sigma_{W^{\perp}})$
(8)

The proof is an exercise in linear algebra and deferred to Section A.1. This
has a nice interpretation. The first two terms are adjustments for the
differences in noise level between $Y^{*}$ and $W^{\perp}$, and the last term
is the usual Mallows’ correction, now with $S(W)$ random. We present refitting
on $W$ and $Y$ as precursors to our bagging results Theorems 2.5 and 2.9,
where $S(W)W$ is the conventional refit for each learner in the bag, and
$S(W)Y$ is an alternative that intuitively should provide better predictions.

We note that Theorem 2.1 we need not use the trace correction, but any
quantity with that expectation. A natural choice is to extend the choice of
Oliveira et al., [19] to our setting, which they showed has smaller variance
than other options:

$\mathrm{GC}^{\\!\perp\\!\\!\\!\perp}(Y):=\|W^{\perp}-S(W)Y\|_{\Theta_{p}}^{2}+\mathrm{tr}(\Theta_{p}\Sigma_{Y})-\mathrm{tr}(\Theta_{p}\Sigma_{Y^{*}})-\|\omega\|_{\Theta_{p}}^{2}/\alpha+2\mathrm{tr}(\Theta_{p}S(W)\Sigma_{Y})$
(9)
$\mathrm{GC}^{\\!\perp\\!\\!\\!\perp}(W):=\|W^{\perp}-S(W)W\|_{\Theta_{p}}^{2}+\mathrm{tr}(\Theta_{p}\Sigma_{Y})-\mathrm{tr}(\Theta_{p}\Sigma_{Y^{*}})-\|\omega\|_{\Theta_{p}}^{2}/\alpha$
(10)

We verify this random correction has smaller variance than the trace
correction in Fig. 10.

The above is for one realization of $W,W^{\perp}$. We can repeat this for many
such bootstrap samples $W_{b},W_{b}^{\perp}$ and average them to get a better
estimate. This is what has been called the coupled bootstrap (CB) by Oliveira
et al., [19]. Despite the name, this technique in fact cleverly creates
marginally decoupled samples, in the sense that they are statistically
independent. So we have the following corollary

######  Corollary 2.2.

, In the setting of Theorem 2.1

$\frac{1}{B}\sum_{b=1}^{B}\mathrm{GC}_{b}^{\\!\perp\\!\\!\\!\perp}(W_{b})$

is an unbiased estimate of $\mathbb{E}\|Y^{*}-S(W)W\|_{\Theta_{p}}^{2}$ and

$\frac{1}{B}\sum_{b=1}^{B}\mathrm{GC}_{b}^{\\!\perp\\!\\!\\!\perp}(Y)$

is an unbiased estimate of $\mathbb{E}\|Y^{*}-S(W)Y\|_{\Theta_{p}}^{2}$.

Thus we suggest generating $W,W^{\perp}$ for small $\alpha$, and using the
above bootstrapped estimator of the error of $S(W)W$ to get an estimate of the
estimators $S(Y)Y$. In all our applications and simulations, we employ this
bootstrap version of the estimator.

We claim the above applies to decision trees. To see a decision tree is an
adaptive linear smoother is in fact quite simple. Conditional on the data it
was fit on, a tree is simply a function $S(\cdot)$ that groups observations
based on their features, and predicts the mean of the observations in each
group. That is,

$S(W)=\begin{pmatrix}n_{1}^{-1}\mathbf{1}_{n_{1}}\mathbf{1}_{n_{1}}^{\top}&0&\cdots&\cdots&0\\\
0&\ddots&\ddots&&\vdots\\\
\vdots&\ddots&n_{\ell}^{-1}\mathbf{1}_{n_{\ell}}\mathbf{1}_{n_{\ell}}^{\top}&\ddots&\vdots\\\
\vdots&&\ddots&\ddots&0\\\
0&\cdots&\cdots&0&n_{L}^{-1}\mathbf{1}_{n_{L}}\mathbf{1}_{n_{L}}^{\top}\end{pmatrix}$

###  2.3 Correlated $Y,Y^{*}$

Our first main result Theorem 2.3 moves beyond the case of $Y^{*}$ and $Y$
being IID vectors conditional on $X$ and simply considers them to be a jointly
Gaussian vector.

It uses the fact that we can regress normal vector $Y^{*}$ on $Y$ to obtain
the decomposition:

$Y^{*}=N+\Gamma Y,\qquad\Gamma=\Sigma_{Y^{*},Y}\Sigma_{Y}^{-1},\qquad
N=Y^{*}-\Gamma Y,\qquad\Delta=(S(W)-\Gamma)Y,\qquad
N^{\perp}=(I-\Gamma)W^{\perp}$

Note that $N\\!\perp\\!\\!\\!\perp Y$. Similarly define

$\Gamma_{W}=\Sigma_{Y^{*},W}\Sigma_{W}^{-1},\qquad
N^{*}=Y^{*}-\Gamma_{W}W,\qquad
N^{W}=(I-\Gamma_{W})W^{\perp},\qquad\Delta^{W}=(S(W)-\Gamma_{W})W$

######  Theorem 2.3.

Suppose

$\begin{pmatrix}Y^{*}\\\ Y\end{pmatrix}\mid
X\sim\mathcal{N}\left(\mu,\Sigma\right),\qquad\Sigma=\begin{pmatrix}\Sigma_{Y^{*}}&\Sigma_{Y^{*},Y}\\\
\Sigma_{Y,Y^{*}}&\Sigma_{Y}\end{pmatrix}$

and denote randomized views $W,W^{\perp}$ of $Y$ defined as in Eq. 5. Then,

$\mathrm{GC}^{cor}(Y):=\|N^{\perp}-\Delta\|_{\Theta_{p}}^{2}+\mathrm{tr}(\Theta_{p}(\Sigma_{N}-\Sigma_{(I-\Gamma)Y}))-\|(I-\Gamma)\omega\|_{\Theta_{p}}^{2}/\alpha+2\mathrm{tr}\left((I-\Gamma)^{\top}\Theta_{p}(S(W)-\Gamma)\Sigma_{Y}\right)$
(11)

is an unbiased estimate of $\mathbb{E}\|Y^{*}-S(W)Y\|_{\Theta_{p}}^{2}$, and

$\mathrm{GC}^{cor}(W):=\|N^{W}-\Delta^{W}\|_{\Theta_{p}}^{2}+\mathrm{tr}(\Theta_{p}(\Sigma_{N^{*}}-\Sigma_{(I-\Gamma_{W})Y}))-\|(I-\Gamma_{W})\omega\|_{\Theta_{p}}^{2}/\alpha$
(12)

is an unbiased estimate of $\mathbb{E}\|Y^{*}-S(W)W\|_{\Theta_{p}}^{2}$.

The proof can be found in Section A.2. Essentially we regress out the
correlations between $Y^{*}$ and $Y$, and then apply the same approach as in
Theorem 2.1. And again we immediately have a bootstrapped version of the
estimator:

######  Corollary 2.4.

In the setting of Theorem 2.3

$\frac{1}{B}\sum_{b=1}^{B}\mathrm{GC}_{b}^{cor}(W_{b})$

is an unbiased estimate of $\mathbb{E}\|Y^{*}-S(W)W\|_{\Theta_{p}}^{2}$ and

$\frac{1}{B}\sum_{b=1}^{B}\mathrm{GC}_{b}^{cor}(Y)$

is an unbiased estimate of $\mathbb{E}\|Y^{*}-S(W)Y\|_{\Theta_{p}}^{2}$.

###  2.4 Bagged estimators

Bootstrap aggregating (bagging) is a common approach to create an ensemble
model from training many base models on different bootstrap samples of the
data, originally proposed by Breiman, [4]. Bagging as been noted for its
ability to smooth out discontinuous estimators, hence it also being known as
bootstrap smoothing [5, 8].

####  2.4.1 Uncorrelated $Y,Y^{*}$

In the special case where $\Sigma_{\omega}=\Sigma_{Y}$, one may view the
fissioned view $W$ of our data as a parametric bootstrap with maximally-
flexible or minimally-biased model as our mean, which is simply $Y$ itself.
Thus it is natural to consider fitting $K$ models to fissioned views $W_{k}$
and estimating the error of the resulting bagged estimator:

$\mathbb{E}\left[\left\|Y^{*}-\frac{1}{K}\sum_{k=1}^{K}S(W_{k})W_{k}\right\|_{\Theta_{p}}^{2}\
\middle|\ X\right],$

or a modified bagged estimator that does a “full refit”:

$\mathbb{E}\left[\left\|Y^{*}-\frac{1}{K}\sum_{k=1}^{K}S(W_{k})Y\right\|_{\Theta_{p}}^{2}\
\middle|\ X\right],$

It is in fact quite straight-forward to extend Theorem 2.1 to the bagged
setting. When viewed in the appropriate way, it is a simple linear algebra
exercise deferred to Section A.3.

######  Theorem 2.5.

Let $W_{k},W_{k}^{\perp}$ be fissioned views of $Y$ from independent
$\omega_{k}$’s, and consider the bagged estimator $K^{-1}\sum_{k}S(W_{k})Y$.
Then an unbiased estimate of
$\mathbb{E}\left\|Y^{*}-\frac{1}{K}\sum_{k=1}^{K}S(W_{k})Y\right\|_{\Theta_{p}}^{2}$
is

$K^{-1}\left[\sum_{k=1}^{K}\mathrm{GC}^{\\!\perp\\!\\!\\!\perp}_{k}(Y)-\mathrm{tr}\left((Z_{K}^{\perp})^{\top}\Theta_{p}Z_{K}^{\perp}\right)\right]$

and an unbiased estimate of
$\mathbb{E}\left\|Y^{*}-\frac{1}{K}\sum_{k=1}^{K}S(W_{k})W_{k}\right\|_{\Theta_{p}}^{2}$
is

$K^{-1}\left[\sum_{k=1}^{K}\mathrm{GC}^{\\!\perp\\!\\!\\!\perp}_{k}(W_{k})-\mathrm{tr}\left((Z_{K}^{\perp})^{\top}\Theta_{p}Z_{K}^{\perp}\right)\right]$

where the subscript $k$ in $\mathrm{GC}^{\\!\perp\\!\\!\\!\perp}_{k}$ denotes
using the $k$-th fissioned view $W_{k},W_{k}^{\perp}$ and
$Z_{K}^{\perp}=Z_{K}(I-K^{-1}\mathbf{1}_{K}\mathbf{1}_{K}^{\top})$ is the
vector of centered predictions.

####  2.4.2 General $Y,Y^{*}$

Just as in the case of traditional estimators, we also have a correction for
bagged estimators in the case where $Y$ and $Y^{*}$ are correlated. We simply
rewrite

$\mathbb{E}\left\|Y^{*}-\frac{1}{K}\sum_{k=1}^{K}S(W_{k})Y\right\|_{\Theta_{p}}^{2}=\mathbb{E}\left\|N-\frac{1}{K}\sum_{k=1}^{K}\Delta_{k}\right\|_{\Theta_{p}}^{2}$

where $\Delta_{k}=(S(W_{k})-\Gamma)Y$. Then the proof from Theorem 2.5 applies
using the estimator from Theorem 2.3:

######  Theorem 2.6.

Assume the setting of $Y\not\\!\perp\\!\\!\\!\perp Y^{*}\mid X$ as in Section
2.3. Let $W_{k},W_{k}^{\perp}$ be fissioned views of $Y$ from independent
$\omega_{k}$’s, and consider the bagged estimator $K^{-1}\sum_{k}S(W_{k})Y$.
Then an unbiased estimate of
$\mathbb{E}\left\|Y^{*}-\frac{1}{K}\sum_{i=1}^{K}S(W_{k})Y\right\|_{\Theta_{p}}^{2}$
is

$K^{-1}\left[\sum_{k=1}^{K}\mathrm{GC}^{cor}_{k}(Y)-\mathrm{tr}\left((Z_{K}^{\perp})^{\top}\Theta_{p}Z_{K}^{\perp}\right)\right]$

and an unbiased estimate of
$\mathbb{E}\left\|Y^{*}-\frac{1}{K}\sum_{i=1}^{K}S(W_{k})W_{k}\right\|_{\Theta_{p}}^{2}$
is

$K^{-1}\left[\sum_{k=1}^{K}\mathrm{GC}^{cor}_{k}(W_{k})-\mathrm{tr}\left((Z_{K}^{\perp})^{\top}\Theta_{p}Z_{K}^{\perp}\right)\right]$

where the subscript $k$ in $\mathrm{GC}^{cor}_{k}$ denotes using the $k$-th
fissioned view $W_{k},W_{k}^{\perp}$ and
$Z_{K}^{\perp}=Z_{K}(I-K^{-1}\mathbf{1}_{K}\mathbf{1}_{K}^{\top})$ is the
vector of centered predictions.

The proof is identical to that of Theorem 2.5 and so we will not repeat it.
Readers will notice this is all for parametric bagging, as opposed to the more
often used nonparametric bagging. One may wonder if these parametric bagged
models will perform similarly to the nonparametric bagged models. We
investigate this in the case of random forests (RF) in Appendix B. In our
simulations, we show that they do indeed perform similarly in terms of
expected prediction error.

It is also worth noting that there is not a version of the nonparametric
bootstrap developed for our setting, and so a parametric bootstrap is also
more fit-to-purpose. There are several existing works proposing modifications
to RF to make it better suited to spatial data [9, 12, 25]. Hengl et al., [12]
adds distance features to account for spatial correlation. Fox et al., [9]
sums an RF prediction and a regression kriging prediction on the residuals of
RF. Saha et al., [25] assumes an autoregressive structure on the covariance
matrix to allow GLS-style fitting of trees in RF. It would be interesting in
future to see how our parametric RF compares to these approaches.

###  2.5 General models

Our initial focus on linear adaptive smoothers is interesting for what it
allows us to do with random forests as discussed above, which is interesting
on its own. However, if one looks carefully at our proofs, you will see that
when we refit on $W$, we need not use the fact that our model is a linear
adaptive smoother. One can easily see this by noting that for any model $g$ we
can define $S(W)=\mathrm{diag}(g(W))\mathrm{diag}(W)^{-1}$, and thus
$S(W)W\overset{a.s.}{=}g(W)$. Therefore we can provide estimators under both
arbitrary $g$ and $Y\not\\!\perp\\!\\!\\!\perp Y^{*}\mid X$.

######  Corollary 2.7.

Assume the setting of $Y\not\\!\perp\\!\\!\\!\perp Y^{*}\mid X$ as in Section
2.3. Additionally let $\Delta_{g}=g(W)-\Gamma_{W}W$. Then an unbiased estimate
of $\mathrm{Err}(g)$ is given by

$\mathrm{GC}^{g}(W)=\mathbb{E}\|N^{\perp}-\Delta_{g}\|_{\Theta_{p}}^{2}+\mathrm{tr}(\Theta_{p}(\Sigma_{N^{*}}-\Sigma_{(I-\Gamma_{W})Y}))-\|(I-\Gamma_{W})\omega\|_{\Theta_{p}}^{2}/\alpha$

And so bootstrapping this also provides an unbiased estimator:

######  Corollary 2.8.

In the setting of Corollary 2.7

$\frac{1}{B}\sum_{b=1}^{B}\mathrm{GC}_{b}^{g}(W)$

is an unbiased estimate of $\mathrm{Err}(g)$.

Corollary 2.8 then gives us a direct generalization of the coupled bootstrap
approach of [19] to the SSN model with arbitrarily correlated Gaussian data.

Of course Corollary 2.7 also then means we can immediately extend our result
to bagging general models $g$ when $Y\not\\!\perp\\!\\!\\!\perp Y^{*}\mid X$.
Now letting

$\mathbb{R}^{n\times K}\ni Z_{K}=\begin{pmatrix}|&|&|\\\
\cdots&Y^{*}-g(W_{k})&\cdots\\\
|&|&|\end{pmatrix},\qquad\mathbb{R}^{n}\ni\bar{Z}_{K}=Y^{*}-\frac{1}{K}\sum_{k=1}^{K}g(W_{k})$

we have the following theorem.

######  Theorem 2.9.

Assume the setting of $Y\not\\!\perp\\!\\!\\!\perp Y^{*}\mid X$ as in Section
2.3. If $W_{k},W_{k}^{\perp}$ are fissioned views of $Y$ from independent
$\omega_{k}$’s, and consider the bagged estimator $K^{-1}\sum_{k}g(W_{k})$.
Then an unbiased estimate of
$\mathbb{E}\left\|Y^{*}-\frac{1}{K}\sum_{i=1}^{K}g(W_{k})\right\|_{\Theta_{p}}^{2}$
is

$K^{-1}\left[\sum_{k=1}^{K}\mathrm{GC}^{g}(W_{k})-\mathrm{tr}\left((Z_{K}^{\perp})^{\top}\Theta_{p}Z_{K}^{\perp}\right)\right]$

where $Z_{K}^{\perp}=Z_{K}(I-K^{-1}\mathbf{1}_{K}\mathbf{1}_{K}^{\top})$ is
the vector of centered predictions.

This lets our estimator extend prediction functions such as the lasso or even
CV-tuned lasso and ridge. We simulate these estimators for independent and
correlated data vectors in Fig. 7 and Fig. 8 respectively, achieving almost
unbiased results.

###  2.6 Estimation of covariance matrix

In practice we will often not know $\Sigma$ a priori and thus must estimate
it. The standard method for this problem is to assume the covariance function
$C(s,t)$ belongs to some parametrized model, and fit those parameters to the
residuals after fitting some prediction function. Typically such models depend
only on the distance between locations $C(d)$ for $d=\|s-t\|_{2}$. For
example, a widely used and effective model is the Matern covariance model. It
is defined as

$C(d)=b+\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{d}{\rho}\right)^{\nu}K_{\nu}\left(\sqrt{2\nu}\frac{d}{\rho}\right)$
(13)

where $b,\nu,\rho,\sigma^{2}$ are the parameters of the Matern model,
$\Gamma(\cdot)$ is the Gamma function, and $K_{\nu}(\cdot)$ is the modified
Bessel function of the second kind. An in-depth discussion of the Matern model
is given in Minasny and McBratney, [17]. One notable strength of this model is
depending on parameter values, one can go from an Exponential covariance
function to a Gaussian covariance function. Software packages exist for
fitting a Matern model, such as the Python package scikitgstat described in
[18], which is the package we use in this paper. They fit the Matern model
with a bounded parameter version of the Levenberg-Marquadt algorithm [14, 16].

Once we have fitted a Matern model, we can use it and the locations
$\\{\ell_{i}\\}$ of our data $Y$ to construct our covariance matrix $\Sigma$.
We simply set $\Sigma_{ij}=C(\|\ell_{i}-\ell_{j}\|_{2})$. If we are
considering SSN, then we need to separate out the measurement noise from the
structured noise. Fortunately, the $b$ in the Matern model is exactly the
measurement noise, known as the nugget in spatial statistics. We simply
subtract it out from the diagonal to get our structural noise covariance
matrix.

An important decision is which prediction function to fit to the data to
generate residuals in the first place. We do not offer theoretical analysis of
this choice here but have found an approach that works well empirically. We
find a good choice is to use the same prediction function to generate
residuals for covariance estimation that one aims to estimate prediction error
for.

Obviously we pay some price in doing this estimation, however Fig. 9 show the
price paid is not too large, at most 5%. We simulated independent data vectors
and used decision trees, ridge regression tuned by CV (RidgeCV), and RF as our
models.

###  2.7 Relationship to Tian and Oliveira

The present work is quite related to Oliveira et al., [19] and our estimators
can be seen as generalizing CB to arbitrarily distributed normal response with
common mean. Oliveira et al., [19] considered the IID normal case for almost
arbitrary prediction function $g$, building upon the earlier work of Tian,
[27] which only considered linear models with feature selection of the form
$S(W)Y$. Our most general result is Corollary 2.8, which provides a direct
generalization of Oliveira et al., [19] to general Gaussian data. We also
present Theorems 2.1, LABEL: and 2.3 which extend Tian, [27] to general
Gaussian data on our way to our bagging result Theorems 2.5, LABEL:, 2.6,
LABEL: and 2.9, and its utility compared to the IID CB presented in Oliveira
et al., [19] is worth discussing.

When considering a general model $g$, it is statistically better to choose
$\alpha$ small, minimizing the perturbation to the model fit, and then
shrinking the resulting large perturbations in $W^{\perp}$ through
bootstrapping. However, this is approach is suboptimal when $g$ is a bagged
model. By taking the approach of fitting the base estimator in a bagged model
on the bootstrap samples $W_{k}$ with a sizeable $\alpha$, we can greatly
increase computational efficiency, even in the IID case. Note that Theorems
2.5, LABEL: and 2.6 only require a single round of bootstrapping. If we
applied the standard CB estimator of Oliveira et al., [19] to a bagged model,
we would bootstrap once in computing CB estimator, and again in the fitting of
the bagged estimator $g$. Thus our results offer a way to generate unbiased
estimates for bagged models not only under a more general data generating
process, but also with a significant computational speedup (proportional to
the number of bootstrap samples).

Oliveira et al., [19] do include a small discussion about general covariance
$\Sigma$ and generalized quadratics $\|\cdot\|_{A}$. However, this is still
for the case where $Y^{*},Y$ are IID conditional on $X$, i.e. NSN in Section
1.2. Our work explicitly allows for correlation between $Y^{*},Y$ via the SSN
model thanks to Theorems 2.3, LABEL: and 2.7. This is crucial for spatial
settings, as it is unreasonable to assume the test points $Y^{*}$ are
independent of the training points $Y$.

##  3 Theoretical properties

Our theoretical results demonstrate our method unbiasedly estimates prediction
error for models fit on noise-elevated data $W$. Simulations show that
choosing $\alpha$ close to 0 results in estimates close to that fit on the
original data $Y$. Here we provide some theoretical results backing up these
empirically observed properties. These are all straight-forward extensions of
results from Oliveira et al., [19] to the non-IID setting.

First we can show the noise-elevated expected error

$\mathrm{Err}_{\alpha}(g)=\mathbb{E}\|Y^{*}-g(W_{\alpha})\|_{2}^{2}$

is smooth in $\alpha$ as $\alpha\to 0$, where
$W_{\alpha}=Y+\sqrt{\alpha}\omega$. In this section only, we use this
subscript notation to denote the scale of the noise elevation as it becomes
necessary to distinguish between different amounts of noise elevation.

######  Proposition 3.1.

If for some $\beta>0$ and integer $k\geq 0$,

$\mathbb{E}\left[\|g(W_{\beta})\|_{\Theta_{p}}^{2}\|W_{\beta}-\mu\|_{\Sigma_{Y}^{-1}}^{2m}\right]<\infty,\qquad
m=0,\dots,k$

then $\alpha\mapsto\mathrm{Err}_{\alpha}(g)$ has $k$ continuous derivatives on
$[0,\beta)$.

The proof is deferred to Appendix C. We repeat an important point made in
Oliveira et al., [19], which is that the above result is achieved without
assuming continuity of $g$. Intuitively, the additional Gaussian noise smooths
out the discontinuities in $g$. Further, in our version of the theorem one may
work with any generalized quadratic $\Theta_{p}$ for the measure of error.

We can also analyze what happens to the estimator in the case where we have
infinite bootstrap samples and take $\alpha\to 0$, that is look at
$\lim_{\alpha\to 0}\mathrm{GCP}_{\alpha}^{\infty}(g)$ for

$\mathrm{GCP}_{\alpha}^{\infty}(g)=\mathbb{E}\left[\|W_{\alpha}^{\perp}-g(W_{\alpha})\|_{\Theta_{p}}^{2}-\|\omega\|_{\Theta_{p}}^{2}/\alpha\mid
Y\right]$

In Oliveira et al., [19] they show that when SURE exists, the above quantity
converges to SURE as $\alpha\to 0$. As one would hope, this results extends to
our setting as well. The proof is a simple extension of the theorem in the IID
case, thus we do not provide a proof here but interested readers should refer
to Oliveira et al., [19].

######  Theorem 3.1.

Consider the setting of Theorem 1.1, but at the $\beta$ noise-elevated level.
That is $\mathbb{E}\|g(W_{\beta})\|_{\Theta}^{2}<\infty$ and
$\mathbb{E}|\nabla_{i}f_{i}(\Theta^{1/2}Y)|<\infty$. Then

$\lim_{\alpha\to
0}\mathrm{GCP}_{\alpha}^{\infty}(g)\overset{a.s.}{=}\mathrm{SURE}(g,\Theta)$

Thus the limit of $\mathrm{GCP}_{\alpha}^{\infty}(g)$ is unbiased for
$\mathrm{Err}(g)$.

In summary, our proposed estimators provide unbiased estimates of the noise-
elevated error $\mathrm{Err}_{\alpha}$, which converges to target error
$\mathrm{Err}$ as $\alpha\to 0$. Furthermore, when SURE exists, our limiting
estimator agrees with SURE.

##  4 Optimism and parametric bootstrap estimators

Although we present our results as stemming from Mallows’ $C_{p}$, we take
time here to relate them to more general notions in covariance penalty
approaches. Recall Efron’s optimism theorem [6, 7] in the IID setting:

$\mathbb{E}\|Y^{*}-g(Y)\|_{2}^{2}=\mathbb{E}\|Y-g(Y)\|_{2}^{2}+2\sum_{i=1}^{n}\mathrm{Cov}(Y_{i},g(Y_{i}))=\mathbb{E}\|Y-g(Y)\|_{2}^{2}+2\mathrm{tr}(\mathrm{Cov}(Y,g(Y)))$
(14)

Estimating this covariance term is how approaches such as Mallows’ and AIC
work. In the spatial setting where $Y^{*},Y\overset{ind}{\sim}P$, the
decomposition remains the same, but in the more general setting where
$Y^{*},Y$ are correlated with eachother and not necessarily IID vectors, the
error decomposees as:

$\mathbb{E}\|Y^{*}-g(Y)\|_{2}^{2}=\mathbb{E}\|Y-g(Y)\|_{2}^{2}+2\mathrm{tr}\left(\mathrm{Cov}(Y,g(Y))\right)-2\mathrm{tr}\left(\mathrm{Cov}(Y^{*},g(Y))\right)+\mathrm{tr}\Sigma_{Y^{*}}-\mathrm{tr}\Sigma_{Y}$
(15)

In light of the above, the approach of this paper can be seen as estimating
this generalized optimism term for Gaussian data. Importantly the optimism is
no longer a function of only the training data, but also of the new
realization $Y^{*}$ through its covariance with the predictions.

With this we can now examine the estimators for IID data in the literature [3,
7, 32] and extend them to our setting. A clear and more detailed explanation
for the IID case is given in Oliveira et al., [19], but we present a brief
summary here borrowing similar notation. These works estimate the optimism in
Eq. 14 with a parametric bootstrap. For a parametric bootstrap sample
$Y^{(b)}\mid Y\sim\mathcal{N}(h(Y),\alpha\sigma^{2}I_{n})$ for some initial
estimate $h(Y)$ and $\alpha\in(0,1]$. The works cited above choose $h(Y)=Y$.
Then using the bootstrap estimate of covariance:

$\widehat{\mathrm{Cov}}=\mathrm{tr}\left((B-1)^{-1}\sum_{b=1}^{B}g(Y^{(b)})(Y^{(b)}-\bar{Y}^{(B)})^{\top}\right)$

we can construct the Efron (Efr) and Breiman-Ye (BY) estimates of error:

$\mathrm{Efr}_{\alpha}(g)=\|Y-g(Y)\|_{2}^{2}+2\widehat{\mathrm{Cov}},\qquad\mathrm{BY}_{\alpha}(g)=\|Y-g(Y)\|_{2}^{2}+\frac{2}{\alpha}\widehat{\mathrm{Cov}}$

We note that Efron, [7] specifically proposes $\mathrm{Efr}_{\alpha}$ with
$\alpha=1$, but we present it in this more general form. With $\alpha=1$ we
have $\mathrm{Efr}_{\alpha}=\mathrm{BY}_{\alpha}$. Now
$\widehat{\mathrm{Cov}}$ is intuitively an unbiased estimate of
$\mathrm{Cov}(Y^{(b)},g(Y^{(b)}))$, and is argued for small $\alpha$ to also
be close to $\mathrm{Cov}(Y,g(Y))$. The essential shortcoming in the IID case
is that $\widehat{\mathrm{Cov}}$ is computed conditional on $Y$, and hence
only captures part of $\mathrm{Cov}(Y^{(b)},g(Y^{(b)}))$, in particular:

$\mathbb{E}[\widehat{\mathrm{Cov}}\mid
Y]=\mathrm{tr}\left(\mathbb{E}\left[\mathrm{Cov}(Y^{(b)},g(Y^{(b)})\mid
Y)\right]\right)$

and thus $\mathbb{E}[\widehat{\mathrm{Cov}}]$ is just the expected conditional
covariance, while by the law of total covariance:

$\mathrm{Cov}(Y^{(b)},g(Y^{(b)}))=\mathbb{E}[\mathrm{Cov}(Y^{(b)},g(Y^{(b)})\mid
Y)]+\mathrm{Cov}(\mathbb{E}[Y^{(b)}\mid Y],\mathbb{E}[g(Y^{(b)})\mid Y])$

Thus it is clear $\mathrm{Efr}_{\alpha}$ will be biased for noise-elevated
covariance for small $\alpha$. The behavior of $\mathrm{BY}_{\alpha}$ is nice
in special cases such as linear smoothers, but in general it is not clear what
it is estimating or how to pick $\alpha$. The same results hold moving to a
setting where $Y^{*}\\!\perp\\!\\!\\!\perp Y$ but
$Y\sim\mathcal{N}(\mu,\Sigma_{Y})$. Therefore, these estimators do not
estimate the noise-elevated error, and it is not clear what they do estimate.

Now if we move to the SSN model where $Y^{*}$ and $Y$ are dependent, we need
to use the parametric bootstrap to instead estimate Eq. 15. So we have to
additionally estimate $\mathrm{Cov}(Y^{*},g(Y))$. A straightforward approach
would be to use two sets of bootstrap samples $Y^{(b)},Y^{*(b)}$ such that

$\begin{pmatrix}Y^{*(b)}\\\ Y^{(b)}\end{pmatrix}\mid
Y\sim\mathcal{N}\begin{pmatrix}\begin{pmatrix}Y\\\
Y\end{pmatrix},\alpha\begin{pmatrix}\Sigma_{Y^{*}}&\Sigma_{Y^{*},Y}\\\
\Sigma_{Y,Y^{*}}&\Sigma_{Y}\end{pmatrix}\end{pmatrix}$

and use the estimate

$\widehat{\mathrm{Cov}}^{*}=\mathrm{tr}\left((B-1)^{-1}\sum_{b=1}^{B}g(Y^{(b)})(Y^{*(b)}-\bar{Y}^{*(B)})^{\top}\right)$

Then the analogous Efron and BY estimators would be

$\mathrm{Efr}_{\alpha}(g)=\|Y-g(Y)\|_{2}^{2}+2\widehat{\mathrm{Cov}}-2\widehat{\mathrm{Cov}}^{*}+\mathrm{tr}\Sigma_{Y^{*}}-\mathrm{tr}\Sigma_{Y}$
$\mathrm{BY}_{\alpha}(g)=\|Y-g(Y)\|_{2}^{2}+\frac{2}{\alpha}\widehat{\mathrm{Cov}}-\frac{2}{\alpha}\widehat{\mathrm{Cov}}^{*}+\mathrm{tr}\Sigma_{Y^{*}}-\mathrm{tr}\Sigma_{Y}$
(16)

For similar reasons, it is not clear these will estimate the optimism at the
noise-elevated level. We can see BY is successful for linear smoothers as

$\mathrm{Cov}(Y^{*},g(Y))=\mathrm{tr}(\Sigma_{Y^{*},Y}S^{\top})=\alpha^{-1}\mathbb{E}\left[\widehat{\mathrm{Cov}}^{*}\mid
Y\right]$

However for more complicated models it seems harder to understand its
behavior. Further, as noted by Oliveira et al., [19], BY can suffer higher
variance when the prediction function $g$ is unstable, whereas estimators that
smooth $g$ with added noise reduce this instability and therefore reduce
estimator variance.

Figure 2: Relative MSE of depth 5 decision tree for the SSN model. ‘Corr X’
refers to the usual way we gnerate $X$ described in Section 6. ‘IID X’ refers
to $X$ generated as IID draws from an isotropic Gaussian
$\mathcal{N}(0,I\sigma^{2})$. ‘GenCp .05’ refers to the estimator proposed in
Corollary 2.7 with $\alpha=.05$. Labels of the form ‘BY u’ refer to the BY
estimator of Eq. 16 with $\alpha=u$.

Fig. 2 illustrates the above described qualities of BY. We simulate as
described in Section 6, only varying how we generate $X$ and $\beta$. We
choose $p=s=200$. On the left is the usual construction detailed in Section 6,
on the right we generate an uncorrelated $X$. We see that when $X$ is
structured that BY underestimates both $\mathrm{Err}(g)$ and
$\mathrm{Err}_{\alpha}(g)$ for larger $\alpha$. On the other hand when $X$ is
IID, we see that although BY is unbiased for all values of $\alpha$, it also
exhibits much higher variance for all values of $\alpha$.

##  5 What is CV estimating under structured noise?

With our approach to the problem, we can also provide some insight into what a
particular CV method is estimating in spatial settings and thus when it may be
suitable to employ.

We begin with some simple insights from decomposing the MSE. Here we depart
from our usual setting of NSN or SSN. Instead we investigate the more typical
setting for CV methods, out-of-sample error. Examining an estimation fold $E$
and prediction fold $P$, we see that for a linear smoother $S$

$\mathbb{E}\|Y-S(X_{E},X_{P})Y\|_{\Theta_{P}}^{2}=\mathbb{E}\|\mu_{P}-S(X_{E},X_{P})Y\|_{\Theta_{P}}^{2}+\mathrm{tr}(\Sigma_{Y}\Theta_{P})-2\mathrm{tr}(S(X_{E},X_{P})\Sigma_{Y}\Theta_{P})$

If we assume our model is (approximately) unbiased, we see that

$\mathbb{E}\|Y-S(X_{E},X_{P})Y\|_{\Theta_{P}}^{2}=\mathrm{tr}(M_{P,E})$

for some $M_{P,E}\in\mathbb{R}^{n\times n}$ which is a function of only
observable quantities and the (assumed known) covariance. The bias of a single
fold from a CV method is due to its $\mathrm{tr}(M_{P,E})$ not matching that
of the true training and test sets.Thus we can simulate this for single splits
$E,P$ from various CV methods to get some insight into when their estimates
match up with the error they aim to estimate.

Fig. 3 illustrates two such comparisons in the case of OLS. We plot the
correction terms of one fold of several CV methods relative to some train/test
split of interest. Namely, we look at a standard random split as well as a
spatial split of the data into train and test sets. We then apply KFCV, SPCV,
and BLOOCV to the training set to get one CV estimation/prediction fold. The
average correction term with standard error bars, normalized by the correction
term of the true train/test split, is plotted. The data $(X,Y)$ are generated
as described in Section 6 with the following differences. We generate $n=2500$
samples placed in a grid on a $50\times 50$ space. We set $p=s=30$. For
smoothing the spikes in columns of $X$, we use a Matern kernel with
$\nu=.5,\rho=5.$

Figure 3: Relative correction term $M_{P,E}$ for various CV methods compared
to different train/test splits. ‘KFCV’ refers to one split of $k$-fold CV, and
‘SPCV’ denotes one split of spatial CV with groups assigned by $k$-means, both
with $k$=5. ‘BLOOCV’ refers to one split of BLOOCV with buffer radius 10.

These initial results point to CV methods providing good estimates of error
when the folds in CV emulate the training and testing sets. We provide the
following theorem to formalize this intuition in a more general setting. The
result is stated for one split, but directly extends to CV. The proof is
deferred to Appendix D.

######  Theorem 5.1.

Consider some $d$-dimensional study area of interest
$R\subseteq\mathbb{R}^{d}$ that has features $X_{r}\in\mathbb{R}^{p}$ and
response $Y_{r}\in\mathbb{R}$ associated with each point $r\in\mathbb{R}$.
Suppose we have a mechanism of generating sample data
$F_{\text{sample}}:\Omega\times\mathbb{N}\times\mathbb{N}\to 2^{R}\times
2^{R}$ to sample points from $R$, and a mechanism specifying how a user would
sub-sample data $F_{\text{user}}:2^{R}\times\mathbb{N}\times\mathbb{N}\to
2^{R}\times 2^{R}$. Consider sampled sets

$\displaystyle
T_{tr}(\omega),T_{ts}(\omega)=F_{\text{sample}}(\omega,n_{tr},n_{ts}),\quad
E(\omega),P(\omega)=F_{\text{user}}(T_{tr}(\omega),n_{e},n_{p}),\quad$

Assume that $F_{\text{sample}},F_{\text{user}}$ are such that, for some
$\omega^{\prime}\overset{d}{=}\omega$,
$T_{tr}(\omega^{\prime}),T_{ts}(\omega^{\prime})=F_{\text{sample}}(\omega^{\prime},n_{e},n_{p})$
we have $E(\omega)\overset{d}{=}T_{tr}(\omega^{\prime})$ for all
$n_{e},n_{p}$. Further, define $r_{ts}$ as a simple random draw from
$T_{ts}(\omega)$ and similarly let $r_{p}$ be a simple random draw from
$P(\omega)$. Assume $r_{ts}\overset{d}{=}r_{p}$. Let $Y,r_{ts},r_{p}$ be such
that $Y_{r_{ts}},Y_{r_{p}}$ have finite second moments. Consider some general
model $g$ such that
$g(X_{tr},Y_{tr},X_{r_{ts}})\overset{p}{\to}g_{\infty}(X_{r_{ts}})$ as
$n_{tr}\to\infty$ for some $g_{\infty}$ such that
$g(X_{E},Y_{E},X_{r_{p}})\overset{p}{\to}g_{\infty}(X_{r_{p}})$. Additionally
let $g$ satisfy

$\displaystyle\mathbb{E}[g(X_{tr},Y_{tr},X_{r_{ts}})^{2}]\to\mathbb{E}[g_{\infty}(X_{r_{ts}})^{2}],\quad\mathbb{E}[g(X_{E},Y_{E},X_{r_{p}})^{2}]\to\mathbb{E}[g_{\infty}(X_{r_{p}})^{2}]$

Then,

$\lim_{n_{tr}\to\infty}\mathbb{E}\left[n_{ts}^{-1}\|Y_{ts}-g(X_{tr},Y_{tr},X_{ts})\|_{2}^{2}\right]=\lim_{n_{e}\to\infty}\mathbb{E}\left[n_{p}^{-1}\|Y_{P}-g(X_{E},Y_{E},X_{P})\|_{2}^{2}\right]$

Intuitively this says that a practitioner should choose a CV method that
splits the training set into folds such that the distribution of the folds
matches the distribution of your training and test sets.

This theorem can also help us explain the results presented in [31]. They
found regular $k$-fold CV to provide unbiased estimates in spatial settings
when the training set comes from a simple random sample of locations. Since a
simple random sample of a simple random sample is still a simple random
sample, the unbiasedness of $k$-fold CV in this setting is unsurprising in
light of the above theorem. However, when Wadoux et al., [31] simulated the
training set from a clustered sampling of locations, none of the CV methods
unbiasedly estimated prediction error. This is obvious in light of the
analysis and discussion in this section. When the original sample is a
clustered sample, none of the CV methods subsample points in a way that
replicates the distribution of the training set or test set locations, and
thus don’t estimate the error accurately.

##  6 Simulations

We now detail the setup for our simulation examples. Unless otherwise
specified, simulation dataset $(X,y)$ is generated in the following way. For
sample size $n=100$, features $p=200$, and sparsity $s=5$. We begin by
generating $X$ and $\beta$ by the following procedure:

  1. 1.

Place samples equally spaced on a $10\times 10$ square with locations
$(c_{x},c_{y})\in[0,10]\times[0,10]$.

  2. 2.

Independently generate columns $X_{i}$ by a two-step process:

    1. (a)

Randomly select $2\log(n)$ locations and place spikes at those locations drawn
from a $\mathrm{Unif}([-3,-1]\cup[1,3])$ distribution.

    2. (b)

Smooth spikes to interpolate values at all other locations using a Matern
kernel. We use Euclidean distance and parameters $\nu=.5,\rho=1.$

  3. 3.

Draw $\beta_{i}\sim\mathrm{Unif}[-1,1]$ for $i=i_{1},\dots,i_{s}$ for randomly
picked indices $i_{j}$. The remaining elements of $\beta$ are set to 0.

From this point forward, $X,\beta$ are held fixed. Then for each iteration of
a simulation, we generate

$\begin{pmatrix}Y\\\ Y^{*}\end{pmatrix}=\begin{pmatrix}X\beta\\\
X\beta\end{pmatrix}+\epsilon$

for noise $\epsilon\sim\mathcal{N}(0,\Sigma)$. $\Sigma$ will be generated
differently depending on the setting: NSN or SSN.

For NSN, that is where $Y\\!\perp\\!\\!\\!\perp Y^{*}\mid X$, $\Sigma$ is
block diagonal where both blocks are equal to
$\Sigma_{Y}=\delta\Sigma_{S}+(1-\delta)I$ for $\delta=.75$.

  * •

$\delta$ controls the relative amounts of structured and measurement noise.
Usually we choose $\delta=.75$.

  * •

$\Sigma_{S}$ is generated from a Matern kernel using Euclidean distance and
parameters $\nu=2.5,\rho=5.$

For SSN, that is where $Y\not\\!\perp\\!\\!\\!\perp Y^{*}\mid X$, then
$\Sigma$ has the same diagonal blocks, and the covariance between $Y^{*},Y$ is
simply $\delta\Sigma_{S}$.

In either case, the covariance matrix is then scaled to achieve a signal-to-
noise ratio (SNR) of $0.4$. That is we scale $\Sigma$ so that
$\mathrm{Var}_{n}(X\beta)/\Sigma_{ii}=0.4$ where $\mathrm{Var}_{n}$ denotes
the empirical variance of the vector $X\beta$.

In some simulations we employ train/test splits. What is referred to as a
random split is simple random sampling of locations. What we refer to as
spatial splits generated by a two-step clustered-sampling approach:

  1. 1.

Split the space into a $5\times 5$ grid of sub-squares.

  2. 2.

For training proportion $p_{tr}$, sample $2p_{tr}$ fraction of sub-squares.

  3. 3.

Sample $np_{tr}$ locations uniformly from the selected sub-squares to create
the training set. All unsampled points become the test set.

This creates clustered samples, replicating what is often found in spatial
datasets.For such simulations, we use the corresponding
$\Theta_{e},\Theta_{p}$ constructions described in Section 2.1. That is, our
error metric is the MSE on the test locations only and for $g$ trained only on
the training locations.

Figures include:

  * •

A red line for the ground truth target $\mathrm{Err}$ estimated by Monte-
Carlo. We plot MSE relative to the ground-truth, so the red line will always
be at 1.

  * •

A dashed gray line for the noise-elevated target $\mathrm{Err}_{\alpha}$ at
$\alpha=.05$, also estimated by Monte-Carlo.

  * •

Error bars, where present, to show standard errors.

Common labels for different estimation methods in figures are:

  * •

‘GenCp’ refers to the variation of GC estimator appropriate for the model in
question.

  * •

‘KFCV’ refers to standard $k$-fold CV with $k=5$.

  * •

‘SPCV’ refers to spatial CV based on groups determined by $k$-means clustering
on locations with $k=5$.

  * •

‘BY’, where present, refers to the generalization of the Breiman-Ye estimator
laid out in Section 4.

  * •

‘Split’, where present, refers to using the train/test split (on $Y$) as an
estimate of $\mathrm{Err}(g)$, which is test locations MSE on $Y^{*}$ for $g$
trained on the training locations of $Y$.

We use a variety of prediction models in our simulations. We detail their
parameter settings here:

  * •

The relaxed lasso is fit with $\lambda=.31$. This $\lambda$ was chosen because
it is close to the middle of the solution path.

  * •

Decision trees are fit to depth 3.

  * •

RF models are fit with 100 trees of depth 3.

  * •

The lasso is fit with $\lambda=.31$. This $\lambda$ was chosen because it is
close to the middle of the solution path.

  * •

LassoCV is the lasso with with $\lambda$ tuned by CV. We tune on a grid of 10
values logarithmically spaced between .01 and 10.

  * •

RidgeCV is the ridge with with $\lambda$ tuned by CV. We tune on a grid of 10
values logarithmically spaced between .01 and 10.

The overarching takeaway from these various simulations is that GC is always
unbiased when the covariance is known, and appear to only suffer at most a
small bias when the covariance is estimated ($\sim 5\%$).

###  6.1 Adaptive linear smoothers

Figure 4: Relative MSE of adaptive linear smoother prediction functions for
the NSN model. We see KFCV underestimates the MSE, while SPCV overestimates
the MSE, while GenCp appears unbiased. Figure 5: Relative MSE of adaptive
linear smoother prediction functions for the SSN model. In this correlated
case, we see both KFCV and SPCV overestimate MSE, while GenCp appears
unbiased.

###  6.2 Bagged adaptive linear smoothers

For these simulations we chose $n=400,p=s=5$ and the following modifications.
We used a non-linear mean function, the Friedman function $\mu(x)=(10\sin(\pi
x_{1}x_{2})+20(x_{3}-0.5)^{2}+10x_{4}+5x_{5})/6$ from [10]. We also applied
random and spatial 25/75 train/test splits.

Figure 6: Relative MSE of bagged prediction functions. Here ‘Split’ refers to
the naive MSE from the train/test split on original response $Y$. The red line
is the MSE on the new sample $Y^{*}$.

###  6.3 General models

Figure 7: Relative MSE of general prediction functions for the NSN model.
Figure 8: Relative MSE of general prediction functions for the SSN model.

###  6.4 Estimating covariance

We set $n=400,p=s=5$. We used a non-linear mean function $\mu$. In particular,
we used the Friedman function $\mu(x)=(10\sin(\pi
x_{1}x_{2})+20(x_{3}-0.5)^{2}+10x_{4}+5x_{5})/6$ [10]. We see that although we
do lose estimation accuracy when we estimate the covariance, it is limited to
about 5%.

Figure 9: Comparing relative MSE using an estimated covariance matrix and an
oracle. We plot results for a depth 3 decision tree as well as ridge tuned by
cross-validation and random forest. EstCov (blue) refers to the estimates
generated using a covariance matrix estimated from the data, while Oracle
(red) refers to estimates generated with the ground-truth covariance matrix.

###  6.5 Random vs trace correction

We set $n=400,p=200,s=5$.

Figure 10: Comparing variance of using random correction vs trace correction
for lasso and ridge tuned by cross-validation.

Fig. 10 shows that the random correction has smaller standard error bars in
both cases while maintaining unbiasedness. Therefore it is more desirable than
the trace correction.

##  7 Application: Nevada geothermal data

In this application, we make use of data from one of our industry
collaborators, Zanskar Geothermal and Minerals. They provided us with
geological and temperature data for the state of Nevada in the USA. In
particular, the dataset has 1,009 sampled locations with temperature data, and
each location has 21 geological features measured for it as well. Examples of
features are the depth at which the measurement was taken, gravity density,
magnetic density, geodetic strain rate, quaternary fault distance. We elide
the details of these features, trusting that our geology colleagues have
measured features predictive of the geological temperature.

We have two goals. The first is to replicate, on real data, simulation results
showing CV methods do disagree with GC, suggesting the inefficacy of CV
methods in spatial settings. The second goal is to demonstrate how one could
use GC estimates as a criteria for model selection.

We fit a RF model using 100 trees of depth 2. We then produce error estimates
using Gen $C_{p}$, KFCV, and SPCV with standard error bars under the SSN
model. Given our theoretical and simulation results, this suggests CV methods
overestimate MSE on this Nevada dataset.

Figure 11: Comparing estimated MSE of various estimation methods. We see that
KFCV and SPCV estimate the error to be significantly higher than our GenCp
method.

For model comparison, we compare various models and parameter settings:

  * •

‘XGB2’: An XGBoost model with 100 estimators of depth 2.

  * •

‘XGB6’: An XGBoost model with 100 estimators of depth 6.

  * •

‘RF2’: An RF model with 100 estimators of depth 2.

  * •

‘RF6’: An RF model with 100 estimators of depth 6.

  * •

‘ENetCV’: An elastic net model tuned by CV.

For each model, we compute in-sample MSE on the full Nevada dataset.

Figure 12: Comparing estimated in-sample MSE of various prediction models.

Fig. 12 shows the results on the Nevada dataset. We see that XGB2 and RF2
performed the best and second best, respectively. Thus GC suggests XGB2 would
be the best predictor of new temperature measurements at these locations.

##  8 Conclusion

We have proposed a method for unbiasedly estimating the error of arbitrary
predictive models on spatially correlated Gaussian data. We also present an
approach to bag any predictive model and provide unbiased estimates of its
prediction error. This approach is appealing because it allows one to
bootstrap in a way that respects the covariance structure of the data, and is
computationally less expensive than applying Corollary 2.7 with $g$ being the
bagged model.

After this analysis, we use our setup to provide some insight into an ongoing
debate around the use of various CV methods for prediction error estimation.
We see that train/test error estimates depend on the covariance of both the
training samples with the predictions as well as the testing samples with the
predictions. This suggests how to use and interpret CV estimates is heavily
dependent on the sampling procedure and ensuring that aligns with the CV
procedure used. For example, if future data are going to be sampled at random
from the study area, then traditional k-fold CV better emulates this sampling
process and should provide better estimates than a spatial approach.

###  8.1 Predicting local maxima

In many settings, the practitioner may not want to know the error on average
over new locations drawn from some distribution, but instead want to know the
model performance at a local maxima of its prediction. Now the choice of
prediction locations depends on the training locations via the model
predictions. GC treats the training and test locations as fixed a priori, and
thus not immediately applicable. Devising a CV splitting strategy that would
replicate this procedure, then perhaps Theorem 5.1 would apply, but such a
strategy does not seem immediately obvious.

###  8.2 Incorporating into CV

Essentially our method allows one to properly do train/test splits for spatial
Gaussian data. One could imagine incorporating this into a CV scheme, where
you use this method to compute multiple estimates of prediction error on
different folds and averaging them. What exactly this CV scheme would be
estimating is not immediately clear and would depend on the CV scheme used.
However, if done such that each train/test split is identically distributed,
then as in the IID case, it would offer a better estimate of predictive error
than a single split. A more thorough investigation of such an approach would
be interesting.

###  8.3 Estimating covariance

As noted in Section 2.6, we find that using the same prediction function for
covariance estimation and prediction error estimation works well in our
simulations. However, it would be interesting to investigate this phenomena
and derive some theory justifying this choice.

###  8.4 Other distributions

Our focus has been correlated Gaussian data, however it would be impactful to
be able to derive analogous results for other distributions and models, such
as common exponential families and their corresponding canonical generalized
linear models (GLMs). Oliveira et al., [19] note that you need only two
independent views of your data $Y$. In Oliveira et al., [20] they work out the
IID Poisson case. In most cases we cannot fission our data as we can with
Gaussians or Poissons, so a different approach may be necessary.

### Acknowledgments

KF is thankful to James Yang and Chen Cheng for useful discussions and
feedback. KF was supported by the National Science Foundation Graduate
Research Fellowship under Grant No. DGE-1656518.

## References

  * Akaike, [1998] Akaike, H. (1998).  Information Theory and an Extension of the Maximum Likelihood Principle, pages 199–213.  Springer New York, New York, NY.
  * Bates et al., [2021] Bates, S., Hastie, T., and Tibshirani, R. (2021).  Cross-validation: what does it estimate and how well does it do it?
  * Breiman, [1992] Breiman, L. (1992).  The little bootstrap and other methods for dimensionality selection in regression: X-fixed prediction error.  Journal of the American Statistical Association, 87(419):738–754.
  * Breiman, [1996] Breiman, L. (1996).  Bagging predictors.  Machine Learning, 24(2):123–140.
  * Buja and Stuetzle, [2016] Buja, A. and Stuetzle, W. (2016).  Smoothing effects of bagging: Von mises expansions of bagged statistical functionals.
  * 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.  Journal of the American Statistical Association, 99(467):619–632.
  * Efron, [2014] Efron, B. (2014).  Estimation and accuracy after model selection.  Journal of the American Statistical Association, 109(507):991–1007.  25346558[pmid], PMC4207812[pmcid].
  * Fox et al., [2020] Fox, E. W., Ver Hoef, J. M., and Olsen, A. R. (2020).  Comparing spatial regression to random forests for large environmental data sets.  PLOS ONE, 15(3):1–22.
  * Friedman, [1991] Friedman, J. H. (1991).  Multivariate Adaptive Regression Splines.  The Annals of Statistics, 19(1):1 – 67.
  * Hastie et al., [2009] Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The elements of statistical learning: data mining, inference and prediction.  Springer, 2 edition.
  * Hengl et al., [2018] Hengl, T., Nussbaum, M., Wright, M., Heuvelink, G., and Gräler, B. (2018).  Random forest as a generic framework for predictive modeling of spatial and spatio-temporal variables.  PeerJ.
  * Leiner et al., [2021] Leiner, J., Duan, B., Wasserman, L., and Ramdas, A. (2021).  Data blurring: sample splitting a single sample.
  * Levenberg, [1944] Levenberg, K. (1944).  A method for the solution of certain non-linear problems in least squares.  Quarterly of Applied Mathematics, 2(2):164–168.
  * Mallows, [1973] Mallows, C. L. (1973).  Some comments on ${C}_{p}$.  Technometrics, 15(4):661–675.
  * Marquardt, [1963] Marquardt, D. W. (1963).  An algorithm for least-squares estimation of nonlinear parameters.  Journal of the Society for Industrial and Applied Mathematics, 11(2):431–441.
  * Minasny and McBratney, [2005] Minasny, B. and McBratney, A. B. (2005).  The matérn function as a general model for soil variograms.  Geoderma, 128(3):192–207.  Pedometrics 2003.
  * Mälicke et al., [2021] Mälicke, M., Möller, E., Schneider, H. D., and Müller, S. (2021).  Scikit-gstat: A scipy flavoured geostatistical variogram analysis toolbox.
  * Oliveira et al., [2021] Oliveira, N. L., Lei, J., and Tibshirani, R. J. (2021).  Unbiased risk estimation in the normal means problem via coupled bootstrap techniques.
  * Oliveira et al., [2022] Oliveira, N. L., Lei, J., and Tibshirani, R. J. (2022).  Coupled bootstrap test error estimation for Poisson variables.
  * Parmentier et al., [2011] Parmentier, I., Harrigan, R. J., Buermann, W., Mitchard, E. T. A., Saatchi, S., Malhi, Y., Bongers, F., Hawthorne, W. D., Leal, M. E., Lewis, S. L., Nusbaumer, L., Sheil, D., Sosef, M. S. M., Affum-Baffoe, K., Bakayoko, A., Chuyong, G. B., Chatelain, C., Comiskey, J. A., Dauby, G., Doucet, J.-L., Fauset, S., Gautier, L., Gillet, J.-F., Kenfack, D., Kouamé, F. N., Kouassi, E. K., Kouka, L. A., Parren, M. P. E., Peh, K. S.-H., Reitsma, J. M., Senterre, B., Sonké, B., Sunderland, T. C. H., Swaine, M. D., Tchouto, M. G. P., Thomas, D., Van Valkenburg, J. L. C. H., and Hardy, O. J. (2011).  Predicting alpha diversity of african rain forests: models based on climate and satellite-derived data do not perform better than a purely spatial model.  Journal of Biogeography, 38(6):1164–1176.
  * Ploton et al., [2020] Ploton, P., Mortier, F., Réjou-Méchain, M., Barbier, N., Picard, N., Rossi, V., Dormann, C., Cornu, G., Viennois, G., Bayol, N., Lyapustin, A., Gourlet-Fleury, S., and Pélissier, R. (2020).  Spatial validation reveals poor predictive performance of large-scale ecological mapping models.  Nature Communications, 11(1):4540.
  * Robins and van der Vaart, [2006] Robins, J. and van der Vaart, A. (2006).  Adaptive nonparametric confidence sets.  The Annals of Statistics, 34(1):229 – 253.
  * Rosset and Tibshirani, [2020] Rosset, S. and Tibshirani, R. J. (2020).  From fixed-x to random-x regression: Bias-variance decompositions, covariance penalties, and prediction error estimation.  Journal of the American Statistical Association, 115(529):138–151.
  * Saha et al., [2020] Saha, A., Basu, S., and Datta, A. (2020).  Random forests for dependent data.
  * Stein, [1981] Stein, C. M. (1981).  Estimation of the Mean of a Multivariate Normal Distribution.  The Annals of Statistics, 9(6):1135 – 1151.
  * Tian, [2020] Tian, X. (2020).  Prediction error after model search.  The Annals of Statistics, 48(2):763 – 784.
  * Tian and Taylor, [2016] Tian, X. and Taylor, J. E. (2016).  Selective inference with a randomized response.
  * Tibshirani, [2015] Tibshirani, R. J. (2015).  Degrees of freedom and Model Search.  Statistica Sinica, 25(3):1265 – 1296.
  * Tibshirani and Taylor, [2012] Tibshirani, R. J. and Taylor, J. (2012).  Degrees of freedom in lasso problems.  The Annals of Statistics, 40(2):1198 – 1232.
  * Wadoux et al., [2021] Wadoux, A. M.-C., Heuvelink, G. B., de Bruin, S., and Brus, D. J. (2021).  Spatial cross-validation is not the right way to evaluate map accuracy.  Ecological Modelling, 457:109692.
  * 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 of unbiasedness of estimators

###  A.1 Proof of Theorem 2.1

The proof is a simple algebra exercise leveraging the independence of
$W^{\perp}$ and $W$

$\displaystyle\ \mathbb{E}\left[\|Y^{*}-SY\|_{\Theta_{p}}^{2}\mid
S(W)=S,X\right]-\mathbb{E}\left[\|W^{\perp}-SY\|_{\Theta_{p}}^{2}\mid
S(W)=S,X\right]$ $\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}\Sigma_{Y})-\mathrm{tr}(\Theta_{p}\Sigma_{W^{\perp}})-2\mathbb{E}\left[(W^{\perp}-\mu)^{\top}\Theta_{p}(\mu-
SY)\mid S(W)=S,X\right]$ $\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}\Sigma_{Y})-\mathrm{tr}(\Theta_{p}\Sigma_{W^{\perp}})+2\mathbb{E}\left[(W^{\perp}-\mu)^{\top}\Theta_{p}SY\mid
S(W)=S,X\right]$ $\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}\Sigma_{Y})-\mathrm{tr}(\Theta_{p}\Sigma_{W^{\perp}})+2\mathbb{E}\left[(W^{\perp}-\mu)^{\top}\Theta_{p}SA^{\perp}(W^{\perp}-\mu)\mid
S(W)=S,X\right]$ $\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}\Sigma_{Y})-\mathrm{tr}(\Theta_{p}\Sigma_{W^{\perp}})+2\mathrm{tr}(\Theta_{p}SA^{\perp}\Sigma_{W^{\perp}})$
$\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}\Sigma_{Y})-\mathrm{tr}(\Theta_{p}\Sigma_{W^{\perp}})+2\mathrm{tr}(\Theta_{p}S\Sigma_{Y})$

The proof is analogous for $S(W)W$ but in that case the last term is 0 since
$W^{\perp}\\!\perp\\!\\!\\!\perp W$ by construction.

###  A.2 Proof of Theorem 2.3

Noting that $\mathbb{E}[N^{\perp}]=\mathbb{E}[N]$, we have

$\displaystyle\ \mathbb{E}\left[\|Y^{*}-SY\|_{\Theta_{p}}^{2}\mid
S(W)=S,X\right]-\mathbb{E}\left[\|N^{\perp}-\Delta\|_{\Theta_{p}}^{2}\mid
S(W)=S,X\right]$ $\displaystyle=$ $\displaystyle\
\mathbb{E}\left[\|N-\Delta\|_{\Theta_{p}}^{2}\mid
S(W)=S,X\right]-\mathbb{E}\left[\|N^{\perp}-\Delta\|_{\Theta_{p}}^{2}\mid
S(W)=S,X\right]$ $\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}(\Sigma_{N}-\Sigma_{N^{\perp}}))-2\mathbb{E}\left[(N^{\perp}-\mathbb{E}[N])^{\top}\Theta_{p}(\mathbb{E}[N]-\Delta)\mid
S(W)=S,X\right]$ $\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}(\Sigma_{N}-\Sigma_{N^{\perp}}))+2\mathbb{E}\left[(N^{\perp}-\mathbb{E}[N])^{\top}\Theta_{p}(S-\Gamma)Y\mid
S(W)=S,X\right]$ $\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}(\Sigma_{N}-\Sigma_{N^{\perp}}))+2\mathbb{E}\left[(W^{\perp}-\mu)^{\top}(I-\Gamma)^{\top}\Theta_{p}(S-\Gamma)A^{\perp}(W^{\perp}-\mu)\mid
S(W)=S,X\right]$ $\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}(\Sigma_{N}-\Sigma_{N^{\perp}}))+2\mathrm{tr}\left((I-\Gamma)^{\top}\Theta_{p}(S-\Gamma)A^{\perp}\Sigma_{W^{\perp}}\right)$
$\displaystyle=$ $\displaystyle\
\mathrm{tr}(\Theta_{p}(\Sigma_{N}-\Sigma_{N^{\perp}}))+2\mathrm{tr}\left((I-\Gamma)^{\top}\Theta_{p}(S-\Gamma)\Sigma_{Y}\right)$

Then noting the first trace term on the rhs can be replaced with

$\mathrm{tr}(\Sigma_{N}-\Sigma_{(I-\Gamma)Y})-\|(I-\Gamma)\omega\|_{\Theta_{p}}^{2}/\alpha$

gives the desired result.

Analogous steps give the result for refitting with $W$, the only difference
being the last term is 0.

###  A.3 Proof of Theorem 2.5

Let $Z_{K}$ denote the matrix of residuals from our $K$ estimators, and
$\bar{Z}_{K}$ be the vector of average residuals, averaged across estimators.
There is a dependence here on the choice of refit vector, but it should be
obvious from context and we suppress that in our notation. In the case of
refitting with $Y$, these quantities would be

$\mathbb{R}^{n\times K}\ni Z_{K}=\begin{pmatrix}|&|&|\\\
\cdots&Y^{*}-S(W_{k})Y&\cdots\\\
|&|&|\end{pmatrix},\qquad\mathbb{R}^{n}\ni\bar{Z}_{K}=Y^{*}-\frac{1}{K}\sum_{k=1}^{K}S(W_{k})Y$

Now we can state our result. In matrix notation we can write

$\left\|Y^{*}-\frac{1}{K}\sum_{i=1}^{K}S(W_{k})Y\right\|_{\Theta_{p}}^{2}=\mathrm{tr}(\bar{Z}_{K}^{\top}\Theta_{p}\bar{Z}_{K})$

It is a simple algebra exercise in expanding out
$\mathrm{tr}\left((Z_{K}^{\perp})^{\top}\Theta_{p}Z_{K}^{\perp}\right)$ to see
that

$K\mathrm{tr}(\bar{Z}_{K}^{\top}\Theta_{p}\bar{Z}_{K})=\mathrm{tr}(Z_{K}^{\top}\Theta_{p}Z_{K})-\mathrm{tr}\left((Z_{K}^{\perp})^{\top}\Theta_{p}Z_{K}^{\perp}\right)$
(17)

Inspecting the right hand side of Eq. 17, we see the first term is the sum of
prediction errors for each individual estimator, and thus estimable by Theorem
2.3. The second term is just a function of the centered predictions, which are
observable. Taking expectations on both sides gives the result.

To see the analogous result for refitting with $W_{k}$, simply replace
$S(W_{k})Y$ with $S(W)W_{k}$ in the definitions of $Z_{K}$, $\bar{Z}_{K}$, and
$Z_{K}^{\perp}$.

##  Appendix B Parametrics v. nonparametric random forest

Here we present an additional simulation comparing the MSE of the usual
nonparametric bootstrapped RF and the parametric bootstrapped RF. We set
$n=400,p=s=30$ with locations placed in a grid on a $20\times 20$ square.

Figure 13: Comparing estimated MSE of parametric vs nonparametric random
forests (RF). ‘NPRF’ is ground-truth error of the nonparametrically
bootstrapped RF, that is the standard RF. ‘PRF’ is ground-truth error of the
parametrically bootstrapped RF. ‘GenCpPRF’ is the estimated error of the
parametrically bootstrapped RF using the estimator proposed in Corollary 2.7.

Fig. 13 shows that the parametric bootstrap RF has very similar MSE to the
usual nonparametric bootstrap RF. Verifying our theory, we see our generalized
$C_{p}$ estimator for parametric bagged models also agrees with both. Thus we
have some confidence in using the parametric bagged error estimates for the
nonparametric RF to achieve significant computational savings compared to the
fully general version of the estimator.

##  Appendix C Proof of Proposition 3.1

This proof is almost a one-for-one replication of the proof given in [19] in
the IID case, merely changing the quadratic form. However we provide the proof
for the non-IID case here for completeness. The proposition is a direct
consequence of the following lemma taking $f(w)=\|\mu-g(w)\|_{\Theta_{p}}^{2}$
and noting that
$\mathbb{E}\|Y^{*}-g(W)\|_{\Theta_{p}}^{2}=\mathbb{E}\|Y^{*}-\mu\|_{\Theta_{p}}^{2}+\mathbb{E}\|\mu-g(W)\|_{\Theta_{p}}^{2}$
since $Y^{*}\\!\perp\\!\\!\\!\perp W$.

######  Lemma C.1.

For $\alpha\geq 0$ let $W_{\alpha}\sim\mathcal{N}(\mu,\Sigma_{W_{\alpha}})$
where

$\Sigma_{W_{\alpha}}=(1+\alpha)\Sigma_{Y}$

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ such that for some $\beta>0$ and integer
$k\geq 0$,

$\mathbb{E}\left[f(W_{\beta})\|W_{\beta}-\mu\|_{\Sigma_{Y}^{-1}}^{2m}\right]<\infty,\qquad
m=0,\dots,k$

then $\alpha\mapsto\mathbb{E}[f(W_{\alpha})]$ has $k$ continuous derivatives
on $[0,\beta)$.

###### Proof.

We begin by proving continuity. Fix $\alpha\in[0,\beta)$. Note that

$\displaystyle\lim_{t\to\alpha}\mathbb{E}[f(W_{t})]$
$\displaystyle=\lim_{t\to\alpha}\int\frac{f(w)}{(2\pi(1+t))^{n/2}|\Sigma_{Y}|^{1/2}}\exp\left\\{-\frac{1}{2(1+t)}\|w-\mu\|_{\Sigma_{Y}^{-1}}^{2}\right\\}dw$
$\displaystyle=\int\lim_{t\to\alpha}\frac{f(w)}{(2\pi(1+t))^{n/2}|\Sigma_{Y}|^{1/2}}\exp\left\\{-\frac{1}{2(1+t)}\|w-\mu\|_{\Sigma_{Y}^{-1}}^{2}\right\\}dw$
$\displaystyle=\mathbb{E}[f(W_{\alpha})]$

where in the second line we use dominated convergence theorem (DCT). This is
justified as

$\displaystyle\frac{f(w)}{(2\pi(1+t))^{n/2}|\Sigma_{Y}|^{1/2}}\exp\left\\{-\frac{1}{2(1+t)}\|w-\mu\|_{\Sigma_{Y}^{-1}}^{2}\right\\}$
$\displaystyle\leq\frac{f(w)}{(2\pi)^{n/2}|\Sigma_{Y}|^{1/2}}\exp\left\\{-\frac{1}{2(1+\alpha)}\|w-\mu\|_{\Sigma_{Y}^{-1}}^{2}\right\\}$
$\displaystyle\leq\frac{f(w)}{(2\pi)^{n/2}|\Sigma_{Y}|^{1/2}}\exp\left\\{-\frac{1}{2(1+\beta)}\|w-\mu\|_{\Sigma_{Y}^{-1}}^{2}\right\\}$

with the last line being integrable by assumption. The proof the first
derivative, which generalizes to any higher derivative, relies on the Leibniz
integral rule for the derivative

$\frac{\partial}{\partial\alpha}\mathbb{E}[f(W_{\alpha})]=-\frac{n}{2(1+\alpha)}\mathbb{E}[f(W_{\alpha})]+\frac{1}{2(1+\alpha)^{2}}\mathbb{E}\left[f(W_{\alpha})\|W_{\alpha}-\mu\|_{\Sigma_{Y}^{-1}}^{2}\right]$

Then the integrands of these expectations are bounded by

$\frac{f(w)}{(2\pi)^{n/2}|\Sigma_{Y}|^{1/2}}\|w-\mu\|_{\Sigma_{Y}^{-1}}^{2m}\exp\left\\{-\frac{1}{2(1+\beta)}\|w-\mu\|_{\Sigma_{Y}^{-1}}^{2}\right\\}$

for $m=0,1$ respectively. Thus by DCT we have continuity of the derivative on
$[0,\beta)$. This general approach applies for higher derivatives using the
Leibniz integral rule and DCT. ∎

##  Appendix D Proof of Theorem 5.1

This is a simple result from swapping limits using Vitali’s convergence
theorem. Recall that Vitali’s convergence theorem states that if
$f_{n}\overset{p}{\to}f$ and $\mathbb{E}[f_{n}^{2}]\to\mathbb{E}[f^{2}]$, then
$f_{n}\overset{\mathcal{L}^{2}}{\to}f$. For
$f_{n}=Y_{r_{ts}}-g(X_{tr},Y_{tr},X_{r_{ts}})$ and
$f=Y_{r_{ts}}-g_{\infty}(X_{r_{ts}})$, they satisfy the conditions of the
theorem by our assumptions on $g$, thus we can conclude

$\lim_{n_{tr}\to\infty}\mathbb{E}[\|Y_{r_{ts}}-g(X_{tr},Y_{tr},X_{r_{ts}})\|_{2}^{2}]=\mathbb{E}[\|Y_{r_{ts}}-g_{\infty}(X_{r_{ts}})\|_{2}^{2}]$
$\lim_{n_{e}\to\infty}\mathbb{E}[\|Y_{r_{p}}-g(X_{E},Y_{E},X_{r_{p}})\|_{2}^{2}]=\mathbb{E}[\|Y_{r_{p}}-g_{\infty}(X_{r_{p}})\|_{2}^{2}]$

Thus we have

$\displaystyle\lim_{n_{tr}\to\infty}\mathbb{E}\left[n_{ts}^{-1}\|Y_{ts}-g(X_{tr},Y_{tr},X_{ts})\|_{2}^{2}\right]$
$\displaystyle=\lim_{n_{tr}\to\infty}\mathbb{E}\left[\|Y_{r_{ts}}-g(X_{tr},Y_{tr},X_{r_{ts}})\|_{2}^{2}\right]$
$\displaystyle=\mathbb{E}\left[\|Y_{r_{ts}}-g_{\infty}(X_{r_{ts}})\|_{2}^{2}\right]$
$\displaystyle=\mathbb{E}\left[\|Y_{r_{p}}-g_{\infty}(X_{r_{p}})\|_{2}^{2}\right]$
$\displaystyle=\lim_{n_{e}\to\infty}\mathbb{E}\left[n_{p}^{-1}\|Y_{P}-g(X_{E},Y_{E},X_{P})\|_{2}^{2}\right]$

where the second to last equality follows from the fact that
$r_{ts}\overset{d}{=}r_{p}$ by assumption, and thus
$Y_{r_{ts}}-g_{\infty}(X_{r_{ts}})\overset{d}{=}Y_{r_{p}}-g_{\infty}(X_{r_{p}})$.