# Cross Validation for Correlated Data

Assaf Rabinowicz
Department of Statistics, Tel-Aviv University, Tel-Aviv, Israel, 69978

Saharon Rosset
Department of Statistics, Tel-Aviv University, Tel-Aviv, Israel, 69978  The
authors gratefully acknowledge Israel Science Foundation, grant 1804/16

(April 05, 2020)

# Cross-Validation for Correlated Data

Assaf Rabinowicz
Department of Statistics, Tel-Aviv University, Tel-Aviv, Israel, 69978

Saharon Rosset
Department of Statistics, Tel-Aviv University, Tel-Aviv, Israel, 69978  The
authors gratefully acknowledge Israel Science Foundation, grant 1804/16

(April 05, 2020)

###### Abstract

K-fold cross-validation (CV) with squared error loss is widely used for
evaluating predictive models, especially when strong distributional
assumptions cannot be taken. However, CV with squared error loss is not free
from distributional assumptions, in particular in cases involving non-i.i.d.
data. This paper analyzes CV for correlated data. We present a criterion for
suitability of standard CV in presence of correlations. When this criterion
does not hold, we introduce a bias corrected cross-validation estimator which
we term $CV_{c},$ that yields an unbiased estimate of prediction error in many
settings where standard CV is invalid. We also demonstrate our results
numerically, and find that introducing our correction substantially improves
both, model evaluation and model selection in simulations and real data
studies.

Keywords: Prediction Error Estimation; Model Selection; Dependent Data; Linear
Mixed Model; Gaussian Process Regression

##  1 INTRODUCTION

Datasets with correlation structures are common in modern statistical
applications in various fields, such as geostatistics (Goovaerts 1999),
genetics (Maddison 1990) and ecology (Roberts et al. 2017). Different modeling
methods address the correlation structure differently. Some modeling methods,
such as Gaussian process regression (Rasmussen and Williams 2006, GPR) and
generalized least squares (Hansen 2007, GLS), utilize explicitly the
correlation structure for achieving better prediction accuracy. Other
predictive models, like random forest (Breiman 2001, RF), gradient boosting
machines (Friedman 2002, GBM) and other machine learning models, do not
consider explicitly the correlation structure but are still potentially able
to utilize the correlation implicitly. The analysis in this paper mainly
focuses on correlation that appears due to latent objects, such as random
effects and random fields as appear in generalized linear mixed models
(Verbeke 1997, GLMM) and generalized Gaussian process regression (Rasmussen
and Williams 2006, GGPR) in clustered, temporal and spatial datasets. A simple
example that demonstrates the way that latent variable realizations affect the
correlation structure is linear mixed models (Verbeke 1997, LMM) with random
intercept for clustered data:

$\displaystyle
y_{i,j}=\boldsymbol{\phi}_{i,j}^{t}\boldsymbol{\beta}+b_{i}+\epsilon_{i,j},\;i\in\\{1,...,I\\},\;j\in\\{1,...,n_{i}\\},$
(1)

where $y_{i,j}\in\mathbb{R}$ is the $j^{th}$ observation for cluster $i,$
$\boldsymbol{\phi}_{i,j}\in\mathbb{R}^{p}$ are the observed fixed effects
covariates, $\boldsymbol{\beta}\in\mathbb{R}^{p}$ is the fixed effects
coefficients vector, $b_{i}\sim N(0,\sigma_{b}^{2})$ are independent random
effects and $\epsilon_{i,j}\sim N(0,\sigma_{\epsilon}^{2})$ are the i.i.d.
errors. Since all the observations in cluster $i$ share the same random effect
realization, they are correlated. For more information see Verbeke (1997).

When it comes to prediction, the question of whether there is correlation
between the observations from the training set and the prediction set — the
sample that is used for the model’s parameter estimation, and the set of
points whose response is predicted based on the trained model, respectively —
plays an important role. For example, in Eq. (1), if the training and
prediction sets are sampled from the same clusters, then once the random
effect realizations are estimated in the model training, they can be utilized
for better prediction accuracy of the dependent variable in the prediction
set. Under some conditions, a predictor that uses the estimated random effect
realizations is the best linear unbiased predictor (BLUP), for more
information see Harville et al. (1976). Another scenario is when the training
and prediction sets are sampled from different clusters. In this scenario,
observations of the prediction set are not correlated with observations of the
training set, and therefore estimating the random effect realizations of the
training set cannot be utilized for achieving better prediction accuracy. Of
course, there are other correlation settings, e.g., when the random effect
realizations of the training and the prediction sets are not the same but
correlated. In addition to the distributional settings that are covered by
GLMM and GGPR, explicit and implicit utilization of the correlation between
the training and prediction sets is common in other distributional settings in
various applications, including applications involving spatial datasets (Ward
and Gleditsch 2018), longitudinal datasets (Hand 2017) and datasets with
hierarchical clustering structure (Raudenbush and Bryk 2002).

The correlation structure of the training and prediction sets, and in
particular the correlation between the training and prediction sets, can
affect the model’s prediction error and therefore should be carefully
addressed when estimating the prediction error. Since many model selection
procedures are based on prediction error estimation, ignoring the correlation
setting may also affect model selection decision.

Various prediction error measures have been used and analyzed under different
correlation settings, e.g., $AIC$-type (Vaida and Blanchard 2005) and
$C_{p}$-type (Hodges and Sargent 2001). This paper focuses on K-fold cross-
validation prediction error estimator (Stone 1974, CV) which is the most
widely used method for estimating prediction error (Hastie et al. 2009). We
introduce a new perspective on CV and propose an applicable framework for
analyzing how CV is affected by the correlation structure in various
distributional settings. In addition, a bias corrected cross-validation
measure, $CV_{c},$ is introduced. $CV_{c}$ is suitable for many scenarios
where CV is biased due to correlations.

Section 2 presents the setting of the problem, the theoretical results of the
proposed approach, and comparison with other methods. Section 3 presents
numerical analyses that support the theoretical results of Section 2.

##  2 CV FOR CORRELATED DATA

###  2.1 K-fold Estimator and Generalization Error

Let
$\\{\boldsymbol{x}_{i}\\}_{i=1}^{n},\;\boldsymbol{x}_{i}\in\mathbb{R}^{p},$ be
an independent and identically distributed (i.i.d.) sample form the
probability distribution function $P_{\boldsymbol{x}}$. Let
$\\{y_{i}\\}_{i=1}^{n},\;y_{i}\in\mathbb{R},$ be a sample that is drawn
independently from the distribution
$P_{y|\boldsymbol{x}=\boldsymbol{x}_{i},\boldsymbol{s}=\boldsymbol{s}_{0}},$
where $\boldsymbol{s}_{0}\in\mathbb{R}^{q}$ is the latent variable realization
that induces correlation structure between $y_{i}\,^{\prime}s.$

Denote $T=\\{y_{i},\boldsymbol{x}_{i}\\}_{i=1}^{n}$ as the training sample for
fitting a predictive model. Also denote
$X=\left[\begin{array}[]{c}-\;\boldsymbol{x}_{1}^{t}-\\\ \vdots\\\
-\;\boldsymbol{x}_{n}^{t}-\end{array}\right]$ and
$\boldsymbol{y}=\left[\begin{array}[]{c}y_{1}\\\ \vdots\\\
y_{n}\end{array}\right].$

One example of this setting is linear mixed models (LMM):

$\boldsymbol{y}=\Phi\boldsymbol{\beta}+Z\boldsymbol{s}+\boldsymbol{\epsilon},$

where $\Phi\in\mathbb{R}^{n\times(p-q)}$ contains the fixed effects
covariates, $Z\in\mathbb{R}^{n\times q}$ contains the random effects
covariates, $\boldsymbol{s}\in\mathbb{R}^{q}$ is a normally distributed random
effects vector with a general covariance matrix and
$\boldsymbol{\epsilon}\in\mathbb{R}^{n}$ is the normally distributed error
term with a general covariance matrix. In this setting $X=\\{\Phi,Z\\}.$ Eq.
(1) is a special case of this setting. Extensions of this model are
generalized linear mixed models (Breslow and Clayton 1993, GLMM), where
$\boldsymbol{y}=g^{-1}(\Phi\boldsymbol{\beta}+Z\boldsymbol{s}+\boldsymbol{\epsilon})$
and $g$ is the link function, as well as hierarchical generalized linear
models (Lee and Nelder 1996, HGLM) where $\boldsymbol{s}$ and
$\boldsymbol{\epsilon}$ do not necessarily follow the normal distribution.
Other examples of this setting that are commonly analyzed and represented
using graphical probabilistic models tools are hidden Markov models and
mixture models (Jordan 2004).

Once a model is fitted, it is natural to evaluate the prediction ability of
the model. In CV prediction error estimator, $T$ is randomly partitioned into
K folds, $T_{1}=\\{y_{i},\boldsymbol{x}_{i}\\}_{i=1}^{n_{1}},$
$T_{2}=\\{y_{i},\boldsymbol{x}_{i}\\}_{i=(n_{1}+1)}^{n_{2}},...,T_{K}=\\{y_{i},\boldsymbol{x}_{i}\\}_{i=(n_{K-1}+1)}^{n_{K}},$111We
assume the observations are randomly ordered, so the folds are exchangeable.
where $n_{k}-n_{k-1}$ is the sample size in fold $k$ and $n_{K}=n.$ For each
$k\in\\{1,...,K\\},$ the model is trained on the entire data except the
$k^{th}$ fold, denoted by $T_{-k}=\bigcup_{j\neq
k}T_{j}=\\{\boldsymbol{y}_{-k},X_{-k}\\}.$ The prediction error of the trained
model is then measured on the holdout fold,
$T_{k}=\\{\boldsymbol{y}_{k},X_{k}\\},$ with respect to some loss function
$L(\cdot\,,\cdot):\;\mathbb{R}\times\mathbb{R}\to\mathbb{R}.$ Thus, each fold
does not train the model used for predicting its outcome. The CV prediction
error estimator is calculated by averaging out the estimated prediction error
across all the folds, i.e.,

$CV=\frac{1}{n}\sum_{k=1}^{K}\sum_{i\in
k^{th}\,fold}L\big{(}y_{i},\hat{y}(\boldsymbol{x}_{i};T_{-k})\big{)},$

where $\hat{y}(\boldsymbol{x}_{i};T_{-k})$ is the predictor of $y_{i},$
constructed by training with $T_{-K}$ and predicting on $\boldsymbol{x}_{i}.$
A special case of CV is leave one out CV (Stone 1974, LOO), which is defined
by setting $K=n,$ i.e., each observation defines a fold and therefore the
model is trained on $n-1$ observations and its prediction error is evaluated
on a single observation. Under some conditions, LOO is superior to other CV
variants in terms of minimizing asymptotic bias and instability (for more
information see Burman (1989); Arlot et al. (2010)), but can be
computationally prohibitive in some settings. A comprehensive comparison of CV
variants in terms of prediction error estimation and models selection can be
found in Zhang and Yang (2015).

In case of i.i.d. sampling, CV is considered to be an estimator of the
generalization error, which is defined here in a wide sense that also covers
settings with correlated data:

$\displaystyle
generalization\;error=\mathbb{E}_{T_{te},T_{tr}}\frac{1}{n_{te}}\sum_{i=1}^{n_{te}}L\big{(}y_{te,i},\hat{y}(\boldsymbol{x}_{te,i};T_{tr})\big{)},$
(2)

where

  * •

$T_{tr}=\\{\boldsymbol{y}_{tr},X_{tr}\\}=\\{y_{tr,i},\boldsymbol{x}_{tr,i}\\}_{i=1}^{n_{tr}}$

    * –

$\boldsymbol{x}_{tr,i}$ is an i.i.d. sample from $P_{\boldsymbol{x}}$

    * –

$y_{tr,i}$ is sampled from
$P_{y|\boldsymbol{x}=\boldsymbol{x}_{tr,i},\boldsymbol{s}=\boldsymbol{s}_{0}}$

  * •

$T_{te}=\\{\boldsymbol{y}_{te},X_{te}\\}=\\{y_{te,i},\boldsymbol{x}_{te,i}\\}_{i=1}^{n_{te}}$

    * –

$\boldsymbol{x}_{te,i}$ is an i.i.d. sample from $P_{\boldsymbol{x}}$

    * –

$y_{te,i}$ is sampled from
$P_{y|\boldsymbol{x}=\boldsymbol{x}_{te,i},\boldsymbol{s}=\boldsymbol{s}_{te}}$

The reason that CV is commonly considered to be an estimator of the
generalization error is the random mechanism that is embedded in the CV
procedure. In this perspective, $T_{k}$ and $T_{-k}$ are equivalent to
$T_{te}$ and $T_{tr},$ respectively. Then, averaging the prediction error,
$\sum_{i\in
k^{th}\,fold}L\big{(}y_{i},\hat{y}(\boldsymbol{x}_{i};T_{-k})\big{)}/\big{(}n_{k}-n_{k-1}\big{)},$
over the different folds, estimates Eq. (2).

It is important to emphasize that unlike $T,$ which is the available dataset
for modeling, $\\{T_{tr},T_{te}\\}$ are used in the paper for demonstrating
different prediction goals, rather than actual datasets.

######  Remark 2.1.

Note, it is typically assumed that $T_{tr}$ is distributed as $T$ and
therefore of size $n.$ In this case, the size of $T_{-k}$ and $T_{tr}$ is
obviously different, and additional bias in CV evaluation is introduced. This
is typically ignored, especially when considering LOO and assuming training
with $n-1$ or $n$ observations carries little difference. In what follows we
also ignore this and implicitly assume that $T_{tr}$ is of the same size as
$T_{-k}$ for all $k$.

When

$\displaystyle\\{\boldsymbol{y}_{k},X_{k}\\}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\\{\boldsymbol{y}_{-k},X_{-k}\\},\;\forall
k\in\\{1,...,K\\}$ (3)

and

$\\{\boldsymbol{y}_{te},X_{te}\\}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\\{\boldsymbol{y}_{tr},X_{tr}\\},$

it is clear that CV is an unbiased estimator of the generalization error,
however a careful analysis is required for the case when the folds are
dependent.

Next we investigate how deviating from the condition in Eq. (3) contributes a
bias to CV with respect to the generalization error. Based on it, a bias
corrected CV estimator, $CV_{c},$ will be presented.

In what follows we limit the discussion to squared error loss function. Also
note, from now on LOO setting will be assumed, and consequently
$n_{te}=1,n_{tr}=n-1,$ however the results are valid to other CV partitioning
settings. In particular, it can easily be seen that the size of the test set,
$n_{te},$ has no bearing on generalization error definition.

###  2.2 A General Formulation of CV Bias

Let

$\displaystyle
w_{cv}=\mathbb{E}_{T_{te},T_{tr}}\big{(}y_{te}-\hat{y}(\boldsymbol{x}_{te};T_{tr})\big{)}^{2}-\frac{1}{n}\mathbb{E}_{T}\|\boldsymbol{y}-\hat{\boldsymbol{y}}_{cv}(T)\|^{2},$

where
$\hat{\boldsymbol{y}}_{cv}(T)=\big{(}\hat{y}(\boldsymbol{x}_{1};T_{-1}),...,\hat{y}(\boldsymbol{x}_{n};T_{-n})\big{)}^{t}$
is the CV predictor of $\boldsymbol{y}.$

An unbiased estimator of the generalization error is:

$\displaystyle CV_{c}=CV+w_{cv}.$ (4)

Before analyzing different correlation settings, we derive a more explicit
expression for $w_{cv}.$ For simplicity, subscript notations in the
expectation operator are frequently omitted. Therefore, unless a specific
object is specified, $\mathbb{E}$ averages all the random variables that it
operates on.

Since

$\displaystyle\mathbb{E}\big{(}y_{te}-\hat{y}(\boldsymbol{x}_{te};T_{tr})\big{)}^{2}=$
$\displaystyle\mathbb{E}\big{(}y_{te}-\mathbb{E}\,y_{te}\big{)}^{2}+\big{(}\mathbb{E}\,y_{te}-\mathbb{E}\,\hat{y}(\boldsymbol{x}_{te};T_{tr})\big{)}^{2}$
$\displaystyle+\mathbb{E}\big{(}\mathbb{E}\,\hat{y}(\boldsymbol{x}_{te};T_{tr})-\hat{y}(\boldsymbol{x}_{te};T_{tr})\big{)}^{2}$
$\displaystyle+2\mathbb{E}\big{(}y_{te}-\mathbb{E}\,y_{te}\big{)}\big{(}\mathbb{E}\,\hat{y}(\boldsymbol{x}_{te};T_{tr})-\hat{y}(\boldsymbol{x}_{te};T_{tr})\big{)},$
$\displaystyle\mathbb{E}\|\boldsymbol{y}-\hat{\boldsymbol{y}}_{cv}(T)\|^{2}=$
$\displaystyle\mathbb{E}\|\boldsymbol{y}-\mathbb{E}\,\boldsymbol{y}\|^{2}+\|\mathbb{E}\,\boldsymbol{y}-\mathbb{E}\,\hat{\boldsymbol{y}}_{cv}(T)\|^{2}+\mathbb{E}\|\mathbb{E}\,\hat{\boldsymbol{y}}_{cv}(T)-\hat{\boldsymbol{y}}_{cv}(T)\|^{2}$
$\displaystyle+2\mathbb{E}\big{(}\boldsymbol{y}-\mathbb{E}\,\boldsymbol{y}\big{)}^{t}\big{(}\mathbb{E}\,\hat{\boldsymbol{y}}_{cv}(T)-\hat{\boldsymbol{y}}_{cv}(T)\big{)},$

and $y_{i}$ and $y_{te}$ have the same marginal distribution, which gives
$\mathbb{E}\|\boldsymbol{y}-\mathbb{E}\,\boldsymbol{y}\|^{2}/n-\mathbb{E}\big{(}y_{te}-\mathbb{E}\,y_{te}\big{)}^{2}=0,$
then:

$\displaystyle w_{cv}=$
$\displaystyle\big{(}\mathbb{E}\,\hat{y}(\boldsymbol{x}_{te};T_{tr})-\mathbb{E}\,y_{te}\big{)}^{2}-\frac{1}{n}\|\mathbb{E}\,\hat{\boldsymbol{y}}_{cv}(T)-\mathbb{E}\,\boldsymbol{y}\|^{2}$
(5)
$\displaystyle\mathbb{E}\big{(}\hat{y}(\boldsymbol{x}_{te};T_{tr})-\mathbb{E}\,\hat{y}(\boldsymbol{x}_{te};T_{tr})\big{)}^{2}-\frac{1}{n}\mathbb{E}\|\hat{\boldsymbol{y}}_{cv}(T)-\mathbb{E}\,\hat{\boldsymbol{y}}_{cv}(T)\|^{2}$
$\displaystyle-2\mathbb{E}\,\mathrm{Cov}\big{(}\hat{y}(\boldsymbol{x}_{te};T_{tr}),y_{te}\big{)}+\frac{2}{n}\mathbb{E}\,\mathrm{tr}\big{[}\mathrm{Cov}\big{(}\hat{\boldsymbol{y}}_{cv}(T),\boldsymbol{y}\big{)}\big{]},$

where $\mathrm{tr}$ is the trace operator and $\mathrm{Cov}$ is the covariance
operator which contains conditional expectation of the dependent variables
$\boldsymbol{y},\;y_{te}$ and $\boldsymbol{y}_{tr}$ given their covariates,
$X,\;\boldsymbol{x}_{te}$ and $X_{tr}$ e.g.,
$\mathbb{E}\,\mathrm{tr}\big{[}\mathrm{Cov}\big{(}\hat{\boldsymbol{y}}_{cv}(T),\boldsymbol{y}\big{)}\big{]}\coloneqq\mathbb{E}_{X}\mathrm{tr}\big{[}\mathrm{Cov}\big{(}\hat{\boldsymbol{y}}_{cv}(T),\boldsymbol{y}|X\big{)}\big{]}.$

The first two lines in Eq. (5) are the differences between the bias and the
variance of $\hat{\boldsymbol{y}}_{cv}(T)$ and
$\hat{y}(\boldsymbol{x}_{te};T_{tr}),$ where the expectation is taken also
over the covariates. The third line relates to the covariances between the
response and its predictor in each scheme — CV prediction error and
generalization error.

###  2.3 Criterion for CV Unbiasedness

Let $P_{T_{te},T_{tr}}$ and $P_{T_{k},T_{-k}}$ be the joint distributions of
$\\{T_{te},T_{tr}\\}$ and $\\{T_{k},T_{-k}\\},$ respectively. Theorem 2.1
describes a simple generic condition when no correction is required for CV.

######  Theorem 2.1.

If $P_{T_{te},T_{tr}}=P_{T_{k},T_{-k}}\;\forall k\in\\{1,..,n\\},$ then
$w_{cv}=0.$

###### Proof.

Since $\\{T_{k},T_{-k}\\}$ were drawn from the same distribution as
$\\{T_{te},T_{tr}\\},$ then an expectation over any transformation of them is
equal, in particular:

$\big{(}\mathbb{E}\,y_{te}-\mathbb{E}\,\hat{y}(\boldsymbol{x}_{te};T_{tr})\big{)}^{2}=\frac{1}{n}\|\mathbb{E}\,\boldsymbol{y}-\mathbb{E}\,\hat{\boldsymbol{y}}_{cv}(T)\|^{2}.$

Similarly

$\displaystyle\mathbb{E}\big{(}\hat{y}(\boldsymbol{x}_{te};T_{tr})-\mathbb{E}\,\hat{y}(\boldsymbol{x}_{te};T_{tr})\big{)}^{2}$
$\displaystyle=\frac{1}{n}\mathbb{E}\|\hat{\boldsymbol{y}}_{cv}(T)-\mathbb{E}\,\hat{\boldsymbol{y}}_{cv}(T)\|^{2}$
$\displaystyle
2\mathbb{E}\,\mathrm{Cov}\big{(}\hat{y}(\boldsymbol{x}_{te};T_{tr}),y_{te}\big{)}$
$\displaystyle=\frac{2}{n}\mathbb{E}\,\mathrm{tr}\big{[}\mathrm{Cov}\big{(}\hat{\boldsymbol{y}}_{cv}(T),\boldsymbol{y}\big{)}\big{]}.$

∎

Theorem 2.1 states a very basic and intuitive condition of CV unbiasedness —
When the CV partitioning preserves the distributional relation between the
prediction set to the training set, then CV is unbiased.

Let
$P_{y_{k}|\boldsymbol{x}_{k},T_{-k}},\;P_{y_{te}|\boldsymbol{x}_{te},T_{tr}}$
and $P_{y_{te}|\boldsymbol{x}_{te},T_{-k}}$ denote the respective conditional
distribution. The condition in Theorem 2.1 can be compacted to the following
one:

$P_{y_{k}|\boldsymbol{x}_{k},T_{-k}}=P_{y_{te}|\boldsymbol{x}_{te},T_{tr}}.$

Moreover, since $T_{tr}$ is assumed to be distributed as $T_{-k},$ then this
condition can be rewritten as follows:

$\displaystyle
P_{y_{k}|\boldsymbol{x}_{k},T_{-k}}=P_{y_{te}|\boldsymbol{x}_{te},T_{-k}}.$
(6)

Of course, when $w_{cv}=0$ the CV is unbiased and therefore suitable. Based on
Theorem 2.1, it is important to stress that CV is suitable not only for the
case that $\\{y_{k},\boldsymbol{x}_{k}\\}_{k=1}^{n}$ are independent, neither
only for the case when $\\{y_{k},\boldsymbol{x}_{k}\\}_{k=1}^{n}$ are
exchangeable, as is commonly mistakenly referred in the literature (Anderson
et al. 2018; Roberts et al. 2017). The biasedness of CV only relates to the
question whether $T_{-k}$ contributes more information for predicting $y_{k}$
than $T_{tr}$ contributes for predicting $y_{te}.$

We can demonstrate the use of Theorem 2.1 for a simple application — using LMM
for predicting new observations from the same clusters that appear in the
training set. In the case $\boldsymbol{s}_{te}=\boldsymbol{s}_{0},$ i.e.,

$\displaystyle\boldsymbol{y}=$
$\displaystyle\Phi\boldsymbol{\beta}+Z\boldsymbol{s}_{0}+\boldsymbol{\epsilon}$
(7) $\displaystyle\boldsymbol{y}_{tr}=$
$\displaystyle\Phi_{tr}\boldsymbol{\beta}+Z_{tr}\boldsymbol{s}_{0}+\boldsymbol{\epsilon}_{tr}$
$\displaystyle y_{te}=$
$\displaystyle\boldsymbol{\phi}_{te}\boldsymbol{\beta}+\boldsymbol{z}_{te}\boldsymbol{s}_{0}+\epsilon_{te},$

and

$\displaystyle X=$
$\displaystyle\\{\Phi,Z\\},\;\Phi\in\mathbb{R}^{n\times(p-q)},\;Z\in\mathbb{R}^{n\times
q}$ $\displaystyle X_{tr}=$
$\displaystyle\\{\Phi_{tr},Z_{tr}\\},\;\Phi_{tr}\in\mathbb{R}^{(n-1)\times(p-q)},\;Z_{tr}\in\mathbb{R}^{(n-1)\times
q}$ $\displaystyle\boldsymbol{x}_{te}=$
$\displaystyle\\{\boldsymbol{\phi}_{te},\boldsymbol{z}_{te}\\},\;\boldsymbol{\phi}_{te}\in\mathbb{R}^{1\times(p-q)},\;\boldsymbol{z}_{te}\in\mathbb{R}^{1\times
q},$

$\boldsymbol{s}_{0}\in\mathbb{R}^{q}$ is the random effect realization vector,
where each entry is a random effect realization for a different cluster and
$\boldsymbol{\epsilon}\in\mathbb{R}^{n},\;\boldsymbol{\epsilon}_{tr}\in\mathbb{R}^{(n-1)},\;\epsilon_{te}\in\mathbb{R}$
are i.i.d. normal errors terms.

As was mentioned previously, the observations in $X,\;X_{tr}$ and
$\boldsymbol{x}_{te}$ are i.i.d.. Also, in this example
$\boldsymbol{y},\;\boldsymbol{y}_{tr}$ and $y_{te}$ were drawn given the same
latent variable realization, $\boldsymbol{s}_{0},$ therefore Theorem 2.1’s
condition —
$P_{y_{k}|\boldsymbol{x}_{k},T_{-k}}=P_{y_{te}|\boldsymbol{x}_{te},T_{-k}}$ —
is satisfied and $w_{cv}=0.$ This use case of predicting new points from the
same clusters that were used in the training data is common, for example see
Gelman (2006).

The principle of CV suitability for the setting in Eq. (7) is discussed in the
LMM literature (Fang 2011; Little et al. 2017), however we did not find any
general mathematical formalization of it. Commonly, CV is avoided in
applications involving correlated data based on the wrong perception that CV
is always unsuitable for these cases (Anderson et al. 2018; Roberts et al.
2017). It is also important to stress that since the condition in Theorem 2.1
only relates to the distributional relation between $\\{T_{te},T_{tr}\\}$ and
$\\{T_{k},T_{-k}\\}$ rather than specifying a distribution, Theorem 2.1 can be
implemented in applications where the distributional settings are not fully
specified, as is common when implementing machine learning algorithms.

###  2.4 CV Correction

Now, consider the setting where Theorem 2.1’s condition is not satisfied,
i.e., when $P_{T_{te},T_{tr}}\neq P_{T_{k},T_{-k}}.$ A simple scenario for
this setting can again be taken from LMM for clustered data, where $T_{te}$
contains a new latent random effects realization i.e.,

$\displaystyle\boldsymbol{y}=$
$\displaystyle\Phi\boldsymbol{\beta}+Z\boldsymbol{s}_{0}+\boldsymbol{\epsilon}$
(8) $\displaystyle\boldsymbol{y}_{tr}=$
$\displaystyle\Phi_{tr}\boldsymbol{\beta}+Z_{tr}\boldsymbol{s}_{0}+\boldsymbol{\epsilon}_{tr}$
$\displaystyle y_{te}=$
$\displaystyle\boldsymbol{\phi}_{te}\boldsymbol{\beta}+\boldsymbol{z}_{te}\boldsymbol{s}_{te}+\epsilon_{te},$

where $\boldsymbol{s}_{te}\neq\boldsymbol{s}_{0}.$ Other relevant components
are defined in the same way as in Eq. (7). In this scenario the correlation
between $\boldsymbol{y}_{tr}$ and $y_{te}$ is different than the correlation
between $\boldsymbol{y}_{-k}$ and $y_{k},$ and a further analysis of $w_{cv}$
is required. This can occur for example when the random effects are intercepts
for clusters (e.g., cities), and the data was collected at one point of time,
but the prediction task is performed on the same clusters at another point of
time.

Here, first we will find an estimator of $w_{cv}$ when
$\hat{\boldsymbol{y}}_{cv}(T)$ is linear in $\boldsymbol{y},$ then $w_{cv}$
for nonlinear predictors will be discussed as well. However, before that, let
us demonstrate how $\hat{\boldsymbol{y}}_{cv}(T)$ can be formalized linearly
in $\boldsymbol{y}$ for some models.

######  Definition 2.1.

$\hat{\boldsymbol{y}}_{cv}(T)$ is linear in $\boldsymbol{y}$ if:

$\displaystyle\hat{\boldsymbol{y}}_{cv}(T)=H_{cv}\boldsymbol{y},$ (9)

where $H_{cv}\in\mathbb{R}^{n\times n}$ does not contain $\boldsymbol{y}$ and
is constructed as follows:

$H_{cv}=\left[\begin{array}[]{cccc}0&h_{1,2}&...&h_{1,n}\\\
h_{2,1}&0&&h_{2,n}\\\ ...\\\ h_{n,1}&h_{n,2}&...&0\end{array}\right],$

$h_{k,k^{\prime}}\in\mathbb{R}\;\forall k,k^{\prime}\in\\{1,...,n\\}.$

The cross validation’s principle that $y_{k}$ is not involved in predicting
itself is reflected by having zeros in the diagonal, therefore
$\hat{y}(\boldsymbol{x}_{k},T_{-k})=\boldsymbol{h}_{k}\boldsymbol{y}_{-k},$
where $\boldsymbol{h}_{k}\in\mathbb{R}^{1\times n-1}$ is the $k^{th}$ row of
$H_{cv}$ without the diagonal element, e.g.,

$\boldsymbol{h}_{1}=[h_{1,2},...,h_{1,n}].$

Of course, $H_{cv}$ can be defined for other K-fold CV settings by adjusting
the dimension of the blocks with the zeros on the diagonal (in LOO the
dimension of each block is $1\times 1,$ however for general K-fold CV is
$n/K\times n/K$).

Examples for linear models are ordinary least squares (OLS), generalized least
squares (GLS), ridge regression, smoothing splines, LMM, Gaussian process
regression (GPR) and kernel regression.

######  Definition 2.2.

Let $\boldsymbol{h}_{te}\in\mathbb{R}^{1\times n-1}$ be the hat vector of
$y_{te},$ constructed from
$\\{X_{tr},\mathrm{Cov}(y_{te},\boldsymbol{y}_{tr})\\}$ in the same way as
$\boldsymbol{h}_{k}$ is constructed from
$\\{X_{-k},\mathrm{Cov}(y_{k},\boldsymbol{y}_{-k})\\}.$

By definition 2.2, $\boldsymbol{h}_{te}$ is the equivalent of
$\boldsymbol{h}_{k}$ for the set $\\{T_{te},T_{tr}\\}$ and they are
distributed the same.

######  Theorem 2.2.

Let $\hat{\boldsymbol{y}}_{cv}(T)=H_{cv}\boldsymbol{y}$ be a linear predictor
of $\boldsymbol{y},$ and
$\hat{y}_{te}(\boldsymbol{x}_{te},T_{tr})=\boldsymbol{h}_{te}\boldsymbol{y}_{tr}$
is its corresponding predictor of $y_{te}.$ Then:

$w_{cv}=\frac{2}{n}\mathbb{E}\big{[}\mathrm{tr}\Big{(}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}-n\boldsymbol{h}_{te}\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}\big{]}.$

###### Proof.

Based on Eq. (5) and by assuming linear predictor,

$-2\mathbb{E}\,\mathrm{Cov}\big{(}\hat{y}(\boldsymbol{x}_{te};T_{tr}),y_{te}\big{)}+\frac{2}{n}\mathbb{E}\,\mathrm{tr}\big{[}\mathrm{Cov}\big{(}\hat{\boldsymbol{y}}_{cv}(T),\boldsymbol{y}\big{)}\big{]}=w_{cv}.$

Therefore, it is only left to show that

$\displaystyle\big{(}\mathbb{E}\,y_{te}-\mathbb{E}\,\boldsymbol{h}_{te}\boldsymbol{y}_{tr}\big{)}^{2}$
$\displaystyle=\frac{1}{n}\|\mathbb{E}\,\boldsymbol{y}-\mathbb{E}\,H_{cv}\boldsymbol{y}\|^{2}$
(10)
$\displaystyle\mathbb{E}\big{(}\boldsymbol{h}_{te}\boldsymbol{y}_{tr}-\mathbb{E}\,\boldsymbol{h}_{te}\boldsymbol{y}_{tr}\big{)}^{2}$
$\displaystyle=\frac{1}{n}\mathbb{E}\|H_{cv}\boldsymbol{y}-\mathbb{E}\,H_{cv}\boldsymbol{y}\|^{2}.$
(11)

Since, by definition $y_{k}$ and $y_{te}$ were drawn from the same distortion,
and $\boldsymbol{y}_{-k}$ and $\boldsymbol{y}_{tr}$ were drawn from the same
distribution, then although $P_{T_{te},T_{tr}}\neq P_{T_{k},T_{-k}},$ still

$\|\mathbb{E}\,\boldsymbol{y}-\mathbb{E}\,H_{cv}\boldsymbol{y}\|^{2}=n\big{(}\mathbb{E}\,y_{k}-\mathbb{E}\,\boldsymbol{h}_{k}\boldsymbol{y}_{-k}\big{)}^{2}=n\big{(}\mathbb{E}\,y_{te}-\mathbb{E}\,\boldsymbol{h}_{te}\boldsymbol{y}_{tr}\big{)}^{2}.$

Therefore Eq. (10) holds. Similarly with Eq. (11). ∎

Using Theorem 2.2, given
$\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}$ and
$\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)},$ an estimator of the
generalization error for a linear predictor is:

$\displaystyle\widehat{CV_{c}}=\frac{1}{n}\big{(}\boldsymbol{y}-H_{cv}\boldsymbol{y}\big{)}^{t}\big{(}\boldsymbol{y}-H_{cv}\boldsymbol{y}\big{)}+\frac{2}{n}\Big{[}\mathrm{tr}\Big{(}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}-n\boldsymbol{h}_{te}\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}\Big{]}.$
(12)

A special case is when
$\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}=0.$ For example, in the
clustered LMM setting that was given in Eq. (8),
$\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}=0$ when the latent
variable realizations of $\boldsymbol{y}_{tr}$ and $y_{te},$ i.e.,
$\boldsymbol{s}_{0}$ and $\boldsymbol{s}_{te},$ are independent. In this case:

$\displaystyle\widehat{CV_{c}}=\frac{1}{n}\big{(}\boldsymbol{y}-H_{cv}\boldsymbol{y}\big{)}^{t}\big{(}\boldsymbol{y}-H_{cv}\boldsymbol{y}\big{)}+\frac{2}{n}\mathrm{tr}\big{[}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\big{]}.$
(13)

In case $\hat{\boldsymbol{y}}_{cv}(T)$ is not a linear predictor and Theorem
2.1’s condition does not hold, then there is no closed form for $w_{cv}.$
However, a correction is required as CV is not an unbiased estimator of the
prediction error in this case. See Section 2.6 for some ad-hoc solutions that
may apply in specific cases.

It is important to note that:

  * •

Theorem 2.2, as well as Eqs. (12) and (13), are valid for any K-fold setting
and not only for LOO. Of course, $H_{cv}$ is based on $K$ and should be
adjusted correspondingly.

  * •

The correction term in $\widehat{CV_{c}}$ depends on covariance matrices,
which are commonly estimated. The effect of using estimated covariance
matrices instead of the true covariance matrices is discussed and analyzed in
Sections 2.7 and 3.

In case there are several alternative models, it is common to use their
estimated prediction errors for selecting the best model.

######  Definition 2.3.

Given a set of models $\mathcal{H},$ the best $\widehat{CV_{c}}$ model is:

$\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV_{c}}(h),$

where $\widehat{CV_{c}}(h)$ is $\widehat{CV_{c}}$ for model $h.$

####  2.4.1 Specifying $\boldsymbol{h}_{te}$ and
$\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}$

The term
$n\boldsymbol{h}_{te}\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}$ in
Eq. (12) is relevant when $\boldsymbol{y}_{tr}$ and $y_{te}$ are correlated
(but not in the same way as $y_{k}$ and $\boldsymbol{y}_{-k})$ e.g., assuming
$T=\\{X,\boldsymbol{y}\\}$ has a hierarchical clustering structure:

$\displaystyle
y_{i,j,r}=\boldsymbol{x}_{i,j,r}^{t}\boldsymbol{\beta}+u_{i}+b_{i,j}+\epsilon_{i,j,r},$
$\displaystyle u_{i}\sim N(0,\sigma^{2}_{u}),\;b_{i,j}\sim
N(0,\sigma^{2}_{b}),\ \epsilon_{i,j,r}\sim N(0,\sigma^{2}_{\epsilon}),$
$\displaystyle\forall\;i\in\\{1,...,I\\},\;j\in\\{1,...,J\\},r\in\\{1,...,R\\}.$

Assume that the prediction goal is to estimate:

$y_{te}=\boldsymbol{x}_{te}^{t}\boldsymbol{\beta}+u_{te}+b_{te}+\epsilon_{te},$

where $u_{te}\in[u_{1},...,u_{I}],$ however $b_{te}$ does not depend on the
realizations of $\\{b_{i,j}\\}_{i\in\\{1,...I\\},\;j\in\\{1,...,J\\}}.$ In
words, the prediction goal is to predict a new observation that relates to one
of the high-level clusters (indexed by $i$), but from a new low-level cluster
(indexed by $j$).

In this setting, since the realization of $u_{te}$ appears in
$\boldsymbol{y},$ then $\mathrm{Cov}(\boldsymbol{y}_{tr},y_{te})\neq 0$ and
$n\boldsymbol{h}_{te}\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}$ is
required for calculating $\widehat{CV_{c}}.$

Although $T_{te}$ is an abstract object, $\boldsymbol{h}_{te}$ and
$\mathrm{Cov}(\boldsymbol{y}_{tr},y_{te})$ can be extracted from the data.
Since $\boldsymbol{h}_{te}$ is distributed as $\boldsymbol{h}_{k}$ (for random
k), then $\boldsymbol{h}_{te}$ can be extracted by selecting a random
$\boldsymbol{h}_{k}$ from $H_{cv}.$ $\mathrm{Cov}(\boldsymbol{y}_{tr},y_{te})$
is extracted correspondingly as
$\mathrm{Cov}(\boldsymbol{y}_{-k},y_{k}|\boldsymbol{b},\boldsymbol{\epsilon}),$
where $\boldsymbol{b}$ and $\boldsymbol{\epsilon}$ are vectors containing
$\\{b_{i,j}\\}_{i\in\\{1,...I\\},\;j\in\\{1,...,J\\}}$ and
$\\{\epsilon_{i,j,r}\\}_{i\in\\{1,...I\\},\;j\in\\{1,...,J\\},\;r\in\\{1,...,R\\}},$
respectively. By conditioning on these vectors, the remaining covariance stems
from $\boldsymbol{u}=[u_{1},...,u_{I}]$ only, corresponding to the covariance
between $\boldsymbol{y}_{tr}$ and $y_{te}.$

Simpler,
$n\boldsymbol{h}_{te}\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}$ can
be replaced by
$\mathrm{tr}\big{(}H_{cv}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y}|\boldsymbol{b},\boldsymbol{\epsilon})\big{)},$
such that:

$\displaystyle\widehat{CV_{c}}=$
$\displaystyle\frac{1}{n}\big{(}\boldsymbol{y}-H_{cv}\boldsymbol{y}\big{)}^{t}\big{(}\boldsymbol{y}-H_{cv}\boldsymbol{y}\big{)}+\frac{2}{n}\Big{[}\mathrm{tr}\Big{(}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}-\mathrm{tr}\Big{(}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}|\boldsymbol{b},\boldsymbol{\epsilon}\big{)}\Big{)}\Big{]}$
$\displaystyle=$
$\displaystyle\frac{1}{n}\big{(}\boldsymbol{y}-H_{cv}\boldsymbol{y}\big{)}^{t}\big{(}\boldsymbol{y}-H_{cv}\boldsymbol{y}\big{)}+\frac{2}{n}\mathrm{tr}\Big{(}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{|}\boldsymbol{u})\Big{)}.$

####  2.4.2 Interpretation of the Results

The correction
$2/n\times\Big{[}\mathrm{tr}\Big{(}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}-n\boldsymbol{h}_{te}\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}\Big{]}$
is intuitive since it expresses the difference between the correlation
structure of the target prediction problem and the correlation structure in
the available dataset, $T.$

Theorems 2.1 and 2.2 emphasize that the question whether CV is biased relates
to the distributional setting and it is indifferent to the implemented
algorithm. The implemented algorithm is expressed only by the values of
$H_{cv}$ in $\widehat{CV_{c}}.$

An interesting example that stresses this understanding is when generalized
least squares (GLS) is implemented in a use case with a correlation setting of
$\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}=0.$ Since GLS estimates
$\mathbb{E}(\boldsymbol{y}|X)$ and does not utilize explicitly the random
effects for achieving better prediction accuracy — as LMM does by estimating
$\mathbb{E}(\boldsymbol{y}|X,\boldsymbol{s})$ — then one may think that CV is
unbiased when GLS is implemented, regardless of the correlation structure of
$\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}.$ However, as was
mentioned, this is wrong since the CV biasedness relates to the distributional
setting rather than to the implemented algorithm. The bias in this case is:

$w_{cv}(GLS)=\mathbb{E}\,\sum_{k=1}^{n}\boldsymbol{x}_{k}^{t}\Big{(}X_{-k}^{t}\mathrm{Cov}(\boldsymbol{y}_{-k},\boldsymbol{y}_{-k}\big{)}^{-1}X_{-k}\Big{)}^{-1}X_{-k}^{t}\mathrm{Cov}\big{(}\boldsymbol{y}_{-k},\boldsymbol{y}_{-k}\big{)}^{-1}\mathrm{Cov}\big{(}\boldsymbol{y}_{-k},y_{k}\big{)}.$

If
$\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})=\sigma^{2}_{\epsilon}I+\rho\mathbbm{1}_{n\times
n},$ where $\mathbbm{1}_{n\times n}$ is a $n$ by $n$ matrix of ones and
$\rho\in\mathbb{R}^{+},$ then:

$\displaystyle w_{cv}$
$\displaystyle=\frac{\rho}{\sigma^{2}_{\epsilon}+\rho\big{(}n-1\big{)}}\sum_{k=1}^{n}\mathbb{E}\,\boldsymbol{x}_{k}\Big{(}X_{-k}^{t}\mathrm{Cov}\big{(}\boldsymbol{y}_{-k},\boldsymbol{y}_{-k}\big{)}^{-1}X_{-k}\Big{)}^{-1}X_{-k}^{t}\mathbbm{1}_{n-1},$

where the identity
$\mathrm{Cov}\big{(}\boldsymbol{y}_{-k},\boldsymbol{y}_{-k}\big{)}^{-1}\mathrm{Cov}\big{(}\boldsymbol{y}_{-k},y_{k}\big{)}=\mathbbm{1}_{n-1}\rho/\Big{(}\sigma^{2}_{\epsilon}+\rho\big{(}n-1\big{)}\Big{)}$
is based on the result by Miller (1981).

The intuition behind the biasedness of CV in this setting is that although GLS
does not utilizes explicitly estimated random effect realizations, its
estimated model coefficients are still affected by the random effect
realizations in $T.$ Similarly, if GLS was replaced by OLS, CV would still be
biased under this correlation setting.

It is also important to emphasize that Theorem 2.1 derives $w_{cv}$ explicitly
($w_{cv}=0$) for any model under its assumed conditions, however
$\widehat{CV_{c}}$ uses an approximated $w_{cv},$ which is only relevant for
linear models.

####  2.4.3 Comparison with Expected Optimism

Below, a comparison between the correction in $\widehat{CV_{c}}$ and the
expected optimism correction (Efron 1986) is presented.

Expected optimism correction was developed in a context of in-sample
prediction error measure:

$in{-}sample\;error=\mathbb{E}_{\boldsymbol{y}^{*},\boldsymbol{y}}\frac{1}{n}\|\boldsymbol{y}^{*}-\hat{\boldsymbol{y}}(X;T)\|^{2},$

where $\boldsymbol{y}^{*}\in\mathbb{R}^{n}$ is identically distributed but
independent copy of $\boldsymbol{y}.$ The in-sample prediction error is
estimated by

$\frac{1}{n}\|\boldsymbol{y}-\hat{\boldsymbol{y}}(X;T)\|^{2}+w,$

where $w$ is the expected optimism:

$w=\mathbb{E}_{\boldsymbol{y}^{*},\boldsymbol{y}}\big{(}\frac{1}{n}\|\boldsymbol{y}^{*}-\hat{\boldsymbol{y}}(X;T)\|^{2}-\frac{1}{n}\|\boldsymbol{y}-\hat{\boldsymbol{y}}(X;T)\|^{2}\big{)}.$

If $\hat{\boldsymbol{y}}(X;T)=H\boldsymbol{y},$ for some hat matrix $H,$ then

$w=\frac{2}{n}\mathrm{tr}\Big{(}H\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}.$

The similarity between $w$ and the correction in $\widehat{CV_{c}}$ in case
$\mathrm{Cov}\big{(}\boldsymbol{y}_{tr},y_{te}\big{)}=0,$ i.e., in case
$w_{cv}=\mathbb{E}\Big{(}2/n\times\mathrm{tr}\big{(}H_{cv}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})\big{)}\Big{)},$
is interesting since it reflects the relation between generalization error and
in-sample error and emphasizes the role of the linearity in this relation.

The fundamental difference between in-sample error and generalization error is
that in the latter, the covariates matrices, $X_{tr},\;\boldsymbol{x}_{te}$
are assumed to be random variables and therefore in generalization error,
unlike in in-sample prediction error:

  1. 1.

$X_{tr}$ and $\boldsymbol{x}_{te}$ are not identical.

  2. 2.

An expectation is taken also over $\\{X_{tr},\;\boldsymbol{x}_{te}\\}.$

As was mentioned in the previous sections, the inner-sampling mechanism of
$\\{T_{k},T_{-k}\\}$ in CV that emulates repeated sampling of
$\\{T_{te},T_{tr}\\}$ from
$\\{P_{\boldsymbol{x}},P_{y|\boldsymbol{x},\boldsymbol{s}}\\}$ expresses these
properties. These properties are also reflected in the correction. Since

$\frac{2}{n}\mathrm{tr}\Big{(}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}=\frac{2}{n}\sum_{k=1}^{n}\boldsymbol{h}_{k}\mathrm{Cov}(y_{k},\boldsymbol{y}_{-k}),$

then
$2/n\times\mathrm{tr}\Big{(}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}$
averages $n$ identically distributed atoms,
$\boldsymbol{h}_{k}\mathrm{Cov}(y_{k},\boldsymbol{y}_{-k}),$ where each one of
them is an unbiased estimator of $w_{cv}.$ Unlike in $w,$ which relates to a
specific covariates matrix realization, $X,$ the atoms in
$2/n\times\mathrm{tr}\Big{(}H_{cv}\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}$
relate to different covariates realizations,
$\\{X_{-k},\boldsymbol{x}_{k}\\}_{k=1}^{n}.$

In addition, when
$\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})=\sigma^{2}_{\epsilon}I,$ while
$w_{cv}=0,$ $w=2\sigma^{2}_{\epsilon}/n\times\mathrm{tr}(H).$ This is since in
CV the sample is partitioned into training and test, however in expected
optimism approach the whole sample is used for both tasks — training and test.

###  2.5 Advanced Correlation Settings

####  2.5.1 Kriging

Many applications with spatial and temporal data are analyzed using random
functions framework, rather than multivariate random variable framework. A
comprehensive review about random functions data analysis can be found in Wang
et al. (2016). Here we focus on interpreting the results from Sections 2.3 and
2.4 in the context of a specific use case in the random functions framework —
Kriging with Gaussian process regression (Rasmussen and Williams 2006, GPR).
Numerical analysis of $\widehat{CV_{c}}$ implementation in real spatial
dataset is presented in Section 3.2.2.

In Kriging (Goovaerts 1999) the goal is to create a climate map on some
surface, $\mathbb{A},$ using climate predictions at a high-resolution grid of
the surface. The predictions at the grid points are based on a predictive
model that was fitted to a sample, $T,$ that was drawn from this surface, but
covers the surface sparsely. In many cases, GPR is the predictive modeling
method that is used for Kriging. In this method, as well as in other
functional data analysis methods, the mean and the covariance of the predicted
variable are formulated as functions. The mean function typically depends on
fixed effects of some covariates (like elevation). The estimated mean function
in GPR is a linear function of $\boldsymbol{y}.$ The covariance function,
which is termed the kernel function,
$\mathcal{K_{\mathbb{A}}}:(\mathbb{A}\times\mathbb{A})\to\mathbb{R},$ measures
the covariance between each two points in $\mathbb{A},$ whether the points are
in the sample $T$ or not. Unlike in the multivariate approach, where the
correlation is induced by a latent random variable realization, in the
functional approach the correlation structure of the surface $\mathbb{A},$ as
it is expressed by $\mathcal{K}_{\mathbb{A}},$ is induced by realization of a
stochastic process instance — latent random function, $\boldsymbol{s}.$ Since
Theorems 2.1 and 2.2 are based on the relation between
$P_{y_{k}|\boldsymbol{x}_{k},T_{-k}}$ and
$P_{y_{te}|\boldsymbol{x}_{te},T_{-k}},$ rather than whether the source of the
correlation between the observations is a latent random variable or latent
random function, these theorems can also be applied here.

Let us consider three scenarios. The first scenario is the classical Kriging
use case, where the observations of both samples, $T_{te}$ and $T,$ are
randomly sampled from the same surface, $\mathbb{A},$ i.e., observations of
both samples are drawn independently from
$\\{P_{\boldsymbol{x}},\;P_{y|\boldsymbol{x},\boldsymbol{s}=\boldsymbol{s}_{0}}\\},$
where $\boldsymbol{s}_{0}$ is the realization of the latent random function
$\boldsymbol{s}$ in the surface $\mathbb{A}.$ In this case Theorem 2.1’s
condition is satisfied and therefore $w_{cv}=0$ and CV is suitable.

The second scenario is when the realization of $\boldsymbol{s}$ is not the
same in $T_{te}$ and $T,$ and therefore while the observations of $T$ follow
$\\{P_{\boldsymbol{x}},\;P_{y|\boldsymbol{x},\boldsymbol{s}=\boldsymbol{s}_{0}}\\},$
the observation in $T_{te}$ follows
$\\{P_{\boldsymbol{x}},\;P_{y|\boldsymbol{x},\boldsymbol{s}=\boldsymbol{s}_{te}}\\}.$
In this case Theorem 2.1’s condition is not satisfied. An example for this
scenario is when $T_{te}$ is sampled from the same surface as $T,$
$\mathbb{A},$ however at a future time-point (e.g., when the goal it to create
a climate map for the next year based on this year’s data). In this case,
$y_{k}$ and $\boldsymbol{y}_{-k},$ which are sampled at the same time point,
are more correlated than $y_{te}$ and $\boldsymbol{y}_{-k},$ which are sampled
at different time-points. Therefore, CV is biased. As we saw in Section 2.4,
if a linear model (such as GPR) is used, then $\widehat{CV_{c}}$ is an
unbiased estimator of the generalization error and therefore should be used
instead of CV.

Another spatial application in this scenario is when $T_{te}$ is sampled from
the surface $\mathbb{A}^{\prime},$ which is different than $\mathbb{A}.$ Since
the surfaces are different, then their latent random function realizations are
different. Therefore, assuming observations in both samples, $T$ and $T_{te},$
were drawn from the same marginal distribution —
$\\{P_{\boldsymbol{x}},\;P_{y|\boldsymbol{x}}\\}$ —, then $\widehat{CV_{c}}$
should be used instead of CV.

Another interesting scenario that is not covered either by Theorem 2.1 or by
Theorem 2.2 is when $\\{y_{k},\boldsymbol{x}_{k}\\}_{k=1}^{n}$ and
$\\{y_{te},\boldsymbol{x}_{te}\\}$ are drawn from different marginal
distributions — $\\{P_{\boldsymbol{x}},\;P_{y|\boldsymbol{x}}\\}$ and
$\\{P^{te}_{\boldsymbol{x}},\;P^{te}_{y|\boldsymbol{x}}\\},$ respectively. An
example for this scenario is when Kriging is used for predicting extrapolated
spatial points with respect to the sample $T.$ It may happen due to sampling
challenges, such as sampling from mountainous and deep marine regions
(Rabinowicz and Rosset 2020). This scenario, which violates the setting that
is assumed in Theorems 2.1 and 2.2 requires further research.

####  2.5.2 Longitudinal Data

Another common setting with correlation structure is longitudinal data — where
there are several subjects that are repeatedly observed over time. In this
setting, due to the temporal orientation, the correlation structure is more
complicated than in a simple clustered data.

Let us consider three scenarios that are equivalent to the three scenarios
given in Section 2.5.1. However, unlike in Section 2.5.1, the multivariate
framework would be considered (rather than the functional framework).

The first scenario is when the prediction goal is to predict a new
observation, $T_{te},$ of one of the subjects in $T,$ sampled at a random
time-point from the same distribution that the time-points of the observations
in $T$ follow. Since $T_{te}$ and $T$ are sampled from the same subjects, then
$\boldsymbol{s}_{te}=\boldsymbol{s}_{0}$ and therefore
$P_{y_{te}|\boldsymbol{x}_{te},T_{-k}}=P_{y_{k}|\boldsymbol{x}_{k},T_{-k}}.$
By Theorem 2.1 this gives $w_{cv}=0,$ and therefore CV is suitable in this
case.

The second scenario is when the prediction goal is predicting a new
observation that relates to a subject that is not in $T$ and therefore
$\boldsymbol{s}_{te}\neq\boldsymbol{s}_{0}$ and
$P_{y_{te}|\boldsymbol{x}_{te},T_{-k}}\neq
P_{y_{k}|\boldsymbol{x}_{k},T_{-k}}.$ For this scenario, given that
$\\{y_{k},\boldsymbol{x}_{k}\\}$ and $\\{y_{te},\boldsymbol{x}_{te}\\}$ were
sampled from the same marginal distribution,
$\\{P_{\boldsymbol{x}},\;P_{y|\boldsymbol{x}}\\},$ and a linear model is
implemented, then by Theorem 2.2, $\widehat{CV_{c}}$ should be used instead of
CV. Numerical analysis of $\widehat{CV_{c}}$ implementation in this scenario
is presented in Section 3.1.

Another scenario is when $\\{y_{k},\boldsymbol{x}_{k}\\}$ and
$\\{y_{te},\boldsymbol{x}_{te}\\}$ are sampled from different marginal
distributions, $\\{P_{\boldsymbol{x}},P_{y|\boldsymbol{x}}\\}$ and
$\\{P^{te}_{\boldsymbol{x}},\;P^{te}_{y|\boldsymbol{x}}\\},$ respectively. The
marginal distributions can be different due to various reasonable prediction
goals, such as forcing the data point in $T_{te}$ to extrapolates the data
points in $T$ with respect to the time variable, which results in
$\boldsymbol{x}_{te}$ and $\boldsymbol{x}_{k}$ being nonidentically
distributed. This scenario violates the assumptions in Theorems 2.1 and 2.2
and therefore requires further research.

###  2.6 Comparison with Other Methods

Several cross-validation variants were proposed for settings involving
correlated data. Some of them were proposed from a perspective that
correlation between the folds causes K-fold CV to underestimate the
generalization error. As was shown above, this perception is wrong in many
scenarios. Other variants are relevant for very specific applications under
various sampling restrictions. Below, several cross-validation variants are
described and compared to $\widehat{CV_{c}}.$

One method is h-blocking (Burman et al. 1994), which is mainly relevant for
spatial data. In h-blocking, in order to reduce the correlation between the
folds, the analyzed surface is partitioned into blocks (folds) that are
separated from each other by some distance, $h.$ As was described above, many
use cases do not require any correlation reduction between the folds and
K-fold CV is suitable, however the use of h-blocking is often recommended
regardless of the predictive problems setup (Roberts et al. 2017). Let us
focus on a scenario when $P_{y_{te}|\boldsymbol{x}_{te},T_{-k}}\neq
P_{y_{k}|\boldsymbol{x}_{k},T_{-k}},$ and therefore the condition in Theorem
2.1 is not satisfied, causing K-fold CV to be biased. In this scenario,
although the h-blocking approach may seem reasonable, in fact, it suffers from
several issues that do not affect $\widehat{CV_{c}}.$ For example, frequently,
creating the separation between the folds requires omitting observations from
the training sample. In addition, the folds that are generated by h-blocking
have different distributions, in particular their distributions are different
than $P_{T_{te}}.$ Therefore, h-blocking may provide a biased prediction error
estimator with respect to the generalization error. Moreover, since some of
the blocks are at the edge of the surface, then the prediction of those blocks
becomes predicting spatial extrapolation, which might be inaccurate and does
not reflect the planned prediction problem. This implication can affect
dramatically the prediction error estimate (Roberts et al. 2017).

Another method is leave cluster out (Rice and Silverman 1991, LCO). This
method is relevant for the case when
$\mathrm{Cov}(\boldsymbol{y}_{tr},y_{te})=0$ and the training set, $T,$ has a
clustered correlation structure. LCO eliminates the correlation between the
folds by defining each cluster as a fold. This method suffers from several
challenges that do not appear in $\widehat{CV_{c}}.$ First, using LCO forces
the number of folds to be equal to the number of clusters. Another issue is
when different clusters contain a substantially different number of
observations, in which case LCO prediction error estimator can be biased with
respect to the generalization error. In addition, validity of LCO is
challenged when some clusters have different distribution than other clusters.
In this case, as opposed to the generalization error definition, the
observations in $X_{k}$ and $X_{-k}$ are nonidentically distributed.

Another cross-validation variant that is relevant for a balanced longitudinal
data setting is leave observation from each cluster out (Wu and Zhang 2002,
LOFCO). This method is relevant for the case when the goal is to predict a new
observation for each one of the subjects that appear in the training set. The
folds partitioning mechanism in LOFCO is that different folds refer to
different time-points, such that each fold contains the observations that were
collocated at the same time-point across all the subjects. This partitioning
is feasible due to the balanced data design assumption. The challenges in this
method are similar to the challenges mentioned above: it requires a balanced
data design, the number of folds are forced by the data structure, $X_{k}$ and
$X_{-k}$ are not identically distributed — in particular their time-points
covariate is nonidentically distributed.

LCO, LOFCO and h-blocking reflect the understanding that nonstandard CV may be
needed in presence of correlations, and offer solutions for very specific
types of datasets and correlation settings that apply to any modeling
technique. In contrast, $\widehat{CV_{c}}$ can be applied in a wide range of
types of datasets and correlation settings, however it is limited to linear
models.

One last estimator that should be compared to $\widehat{CV_{c}}$ is

$\displaystyle\frac{\frac{1}{n}\|\boldsymbol{y}-H\boldsymbol{y}\|^{2}}{\Big{[}1-\frac{1}{n}\mathrm{tr}\Big{(}H\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}\Big{]}^{2}},$
(14)

where $H$ satisfies, $\hat{\boldsymbol{y}}(X;T)=H\boldsymbol{y}.$ This
estimator has the same spirit as generalized CV (Craven and Wahba 1978, GCV)
which approximates LOO for i.i.d. data as follows:

$GCV=\frac{\frac{1}{n}\|\boldsymbol{y}-\hat{\boldsymbol{y}}(X;T)\|^{2}}{\Big{[}1-\frac{1}{n}\mathrm{tr}\big{(}H\big{)}\Big{]}^{2}}.$

Eq. (14) was first proposed by Altman (1990) in context of time-series and
later by Opsomer et al. (2001) in context of spatial data analysis. Altman’s
motivating application was selecting the best bandwidth for modeling the trend
in the data. He argues that replacing $\mathrm{tr}\big{(}H\big{)}$ in GCV by
$\mathrm{tr}\Big{(}H\mathrm{Cov}\big{(}\boldsymbol{y},\boldsymbol{y}\big{)}\Big{)}$
denoises the correlation between observations. Although Eq. (14) and
$\widehat{CV_{c}}$ are different and developed for different applications,
they share some similarity in suggesting a corrected version for CV rather
than controlling the partitioning scheme. Analyzing numerically Eq. (14) did
not yield promising results, in fact its approximated expectation did not
converge to the same scale as the approximated generalization error, CV and
$\widehat{CV_{c}}.$ We conclude that the heuristic argument behind the
derivation of Eq. (14) as an estimator of the generalization error does not
hold in practical cases where the rigorous $\widehat{CV_{c}}$ gives valid
results.

###  2.7 Estimating $\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})$

In most applications $\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})$ is unknown
in advance and therefore in order to implement $\widehat{CV_{c}},$ it should
be estimated. In many cases, such as when GLS is implemented, the covariance
matrix is already estimated in the model fitting part and can be also used in
$\widehat{CV_{c}}.$ The use of estimated covariance matrix instead of the true
covariance matrix can add variance to $\widehat{CV_{c}}$ and even can make
$\widehat{CV_{c}}$ biased, especially when the estimation of the covariance
matrix is imprecise, due to small sample size or other reasons. Several papers
analyze the effect of plugging in the estimated covariance matrix instead of
the true one in context of prediction error estimation, in various CV versions
(Altman 1990; Francisco-Fernandez and Opsomer 2005) or AIC-type versions
(Vaida and Blanchard 2005; Liang et al. 2008; Greven and Kneib 2010;
Rabinowicz and Rosset 2020). Some papers derive extra corrections for their
prediction error estimators which are theoretically valid in certain use
cases, but are not applicable due to computational complexity (Liang et al.
2008), others show numerically that the prediction error estimator that are
based on the estimated covariance matrices instead of the true covariance
matrices performs well. Using simulations and real data analysis, we show in
Section 3, that the effect of using the estimated covariance matrices in
$\widehat{CV_{c}}$ on its bias and variance is small, especially when the
sample size of the training set is not very small. These results support the
results of Altman (1990); Vaida and Blanchard (2005); Francisco-Fernandez and
Opsomer (2005).

##  3 NUMERICAL RESULTS

This section compares $\widehat{CV_{c}}$ and CV, with respect to the
approximated generalization error, using simulation and real datasets
analyses. Different prediction goals and correlation structures are analyzed.
Relevant datasets and code can be found in https://github.com
/AssafRab/CVc.

###  3.1 Simulation

The dependent variable, $\boldsymbol{y}\in\mathbb{R}^{n},$ was sampled from
the following model:

$\displaystyle
y_{i,j,k}=0.1\sum_{r=1}^{9}x_{i,j,k,r}+u_{i}+\sum_{r=1}^{2}z_{i,j,k,r}b_{i,j,r}+\epsilon_{i,j,k},$
$\displaystyle i\in\\{1,...,I\\},\;j\in\\{1,...,5\\},\;k\in\\{1,...,10\\},$

where

  * •

the random effects, $u_{i},\;\boldsymbol{b}_{i,j}=[b_{1},b_{2}]$ and
$\epsilon_{i,j,k}$ are independent and distributed as follows:

$\displaystyle
u_{i}\overset{ind}{\sim}N\big{(}0,3^{2}\big{)},\;\boldsymbol{b}_{i,j}=[b_{i,j,1},b_{i,j,2}]\in\mathbb{R}^{2}\overset{ind}{\sim}N(0,\begin{bmatrix}3^{2}&0\\\
0&1\end{bmatrix}),\;\epsilon_{i,j,k}\overset{ind}{\sim}N\big{(}0,1\big{)},$

  * •

the covariates, $\\{\boldsymbol{x}_{r}\\}_{r=1}^{9}$ are:

    * –

$x_{i,j,k,1}=1,\;x_{i,j,k,2}=k,\;\forall i,j,k$ are the intercept and the time
covariates,

    * –

$x_{i,j,k,r}=\eta_{i}+\delta_{i,j,k},\;\forall i,j,k$ and $\forall
r\in\\{3,...,9\\},$ where
$\eta_{i}\overset{ind}{\sim}N(0,1),\;\delta_{i,j,k}\overset{ind}{\sim}N(0,1)$
are independent,

  * •

$z_{i,j,k,1}=1,\;z_{i,j,k,2}=k\;\forall i,j,k$ are the covariates for random
intercept and random slope.

This settings was simulated for three different sample sizes, $n=300/400/500$
(where $I$ varies respectively to $6/8/10$).

This is a hierarchical clustered structured with $I$ clusters of $50$
observations each, where within each cluster there are five subclusters of ten
observations each. There is a random intercept for the high-level clusters and
random intercept and slope for the subclusters. Therefore the covariance of
$\boldsymbol{y}|\boldsymbol{x}_{1},...,\boldsymbol{x}_{9}$ is:

$\mathrm{Cov}(y_{i,j,k},y_{i^{\prime},j^{\prime},k^{\prime}})=9\times
1_{[i=i^{\prime}]}+1_{[i=i^{\prime},j=j^{\prime}]}(9+1\times k\times
k^{\prime})+1_{[i=i^{\prime},j=j^{\prime},k=k^{\prime}]}.$

GLS was fitted using LOO for eight nested models. Model $1$ contains the
intercept and time, Model $2$ also contains $\boldsymbol{x}_{3},$ and so on.
Model $8$ contains all the covariates.

####  3.1.1 Estimating Prediction Error

$T_{tr}$ is a random subset of $T$ of size $n-1.$ $T_{te}$ is a single
observation, drawn from the same marginal distribution, however with new
independent realizations of all the random effects. Therefore, when
implementing CV, while in $T:$

$\mathrm{Cov}(\boldsymbol{y}_{-k},y_{k})\neq 0\;\;\forall\;k\in\\{1,...,n\\},$

in $\\{T_{tr},T_{te}\\}:$

$\mathrm{Cov}(\boldsymbol{y}_{tr},y_{te})=0.$

The generalization error was approximated by averaging
$\big{(}y_{te}-\hat{y}(\boldsymbol{x}_{te};T_{tr})\big{)}^{2},$ based on
$1000\times n$ samples of $\\{T_{te},T_{tr}\\}.$ The densities of CV and
$\widehat{CV_{c}}$ were approximated based on $1000$ samples of $T.$ Since in
this setting $\mathrm{Cov}(\boldsymbol{y}_{tr},y_{te})=0,$ the correction in
this case is
$2/n\times\mathrm{tr}\big{(}H_{cv}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})\big{)}.$

Figure 1 shows the distribution of CV and $\widehat{CV_{c}}$ including their
means for the saturated model when $n=400,$ compared to the generalization
error. Two versions of CV and $\widehat{CV_{c}},$ are presented — when the
variance parameters are known, and when they are estimated.

Figure 1: Densities of CV and $\widehat{CV_{c}}$ and their means compared to
the generalization error (in black). Two scenarios are presented: when the
variance parameters are known (solid line) and when they are estimated (dashed
line).

As can be seen, $\widehat{CV_{c}}$ is an unbiased estimator of the
generalization error, while CV estimator is biased. This is true for both
versions, when the variance parameters are known and when they are estimated.
Also, as expected, estimating the variance parameters increases the variance
of $\widehat{CV_{c}}.$ Still, the density of the $\widehat{CV_{c}}$ version
with the estimated variance parameters is similar to the version with the
known variance parameters — their averages are $60.14$ and $60.09$ (compared
to $60.00$ for the approximated generalization error), their standard
deviations are $12.07$ and $13.02,$ respectively.

In order to asses the performance of $\widehat{CV_{c}}$ version with the
estimated variance parameters, compared to the version with the known variance
parameters, a two sample Anderson-Darling test (Anderson and Darling 1952) was
used. The tested statistic is:

$\widehat{CV_{c}}-generalization\;error,$

where one sample uses $\widehat{CV_{c}}$ version with the estimated variance
parameters, and the other sample uses $\widehat{CV_{c}}$ version with the true
variance parameters. Implementing the function ad.test of the package kSamples
in R software, the p-value of the test is $0.33.$ The result indicates that in
this setting there is no evidence for significant difference between the
distribution of $\widehat{CV_{c}}-generalization\;error$ when the variance
parameters are known in advance or estimated.

Figure 2 shows the same analysis as in Figure 1, however for various sample
sizes.

Figure 2: Densities of CV and $\widehat{CV_{c}}$ and their means compared to
the generalization error (in black) for different sample sizes. Two scenarios
are presented: when the variance parameters are known (solid line) and when
they are estimated (dashed line).

As can be seen in Figure 2, the variance of the two $\widehat{CV_{c}}$
versions becomes similar as the sample size increases. While for $n=300$ the
standard deviations of the version with the true variance parameters is
$14.43$ and for the estimated is $16.04,$ when $n=500$ they are $11.25$ and
$11.93,$ respectively. This is expected, since larger sample sizes provide
more accurate variance parameters estimators. Also, as the sample size
increases the CV bias, $w_{cv},$ decreases. This is specific to our setting,
where $\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})$ becomes more sparse as the
sample size increases and therefore
$\mathbb{E}\big{[}2/n\times\mathrm{tr}\big{(}H_{cv}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})\big{)}\big{]}$
decreases.

####  3.1.2 $\mathrm{Cov}(\boldsymbol{y}_{tr},y_{te})\neq 0$ and Model
Selection

In this simulation $\widehat{CV_{c}}$ is analyzed under a different prediction
goal than in Section 3.1.1: estimating the generalization error of new
observations from the same high-level clusters, but from a different
subcluster, i.e.,

$y_{te}=0.1\sum_{r=1}^{9}x_{te,r}+u_{i}+\sum_{r=1}^{2}z_{te,r}b_{te,r}+\epsilon_{te},$

where $u_{i}$ is a random effect realization that appears in the training
sample, and $\boldsymbol{b}_{te}=[b_{te,1},b_{te,2}]^{t}$ is a new random
effect realization vector. The new prediction goal derives a different
correction — since in this setting $y_{te}$ is correlated with
$\boldsymbol{y}_{tr}$ through the random effect $u,$ then the correction in
this case is
$2/n\times\mathrm{tr}\big{(}H_{cv}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y}|u_{1},u_{2}...,u_{I})\big{)},$
for more details see Section 2.4.1.

Two different predictive algorithms were implemented:

  * •

LMM, which is the model that should be used in this setting under the
normality assumption, since it utilizes explicitly the correlation between
$y_{te}$ and $\boldsymbol{y},$ using the relevant estimated random effects
realizations, $\hat{u}_{1},...,\hat{u}_{I}.$

  * •

GLS, which is inferior to LMM under the normality assumption, however unlike
LMM, GLS is relevant for use cases where distributional assumptions cannot be
taken. This is since GLS can be interpreted as an extension of the least
squares algorithm, which does not rely on distributional assumptions.

Figure 3 shows the model density estimation of $\widehat{CV_{c}}$ and CV
compared to the approximated generalization error in the new setting when LMM
and GLS are implemented.

(a)

(b)

Figure 3: (3(a)): Using LMM, each curve refers to the density of the relevant
estimator, CV and $\widehat{CV_{c}}.$ Two scenarios are presented: when the
variance parameters are known (solid line) and when they are estimated (dashed
line). Their means are compared to the generalization error (in black).
(3(b)): Similarly, the analysis is presented when GLS is implemented instead
of LMM.

As can be seen in Figure 3, $\widehat{CV_{c}}$ estimates the generalization
error unbiasedly, while CV underestimates it. The CV bias in Figure 3(b) is
relatively small with respect to the biases in Figure 3(a) and Figure 2:

  * •

The reason that the bias in Figure 3(b) is smaller than the bias in Figure 2
although in both settings GLS is implemented, is that the deviation of
$P_{y_{te}|\boldsymbol{x}_{te},T_{tr}}$ from
$P_{y_{k}|\boldsymbol{x}_{k},T_{-k}}$ is smaller in the setting of Figure
3(b), as expressed in the corrections —
$2/n\times\mathrm{tr}\big{(}H_{cv}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y}|u_{1},u_{2}...,u_{I})\big{)}$
in the setting of Figure 3(b) and
$2/n\times\mathrm{tr}\big{(}H_{cv}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y})\big{)}$
in the setting of Figure 2.

  * •

The reason that the bias in Figure 3(b) is smaller than the bias in Figure
3(a) although the deviation of $P_{y_{te}|\boldsymbol{x}_{te},T_{tr}}$ from
$P_{y_{k}|\boldsymbol{x}_{k},T_{-k}}$ is the same in both settings is that in
Figure 3(a) LMM is implemented and therefore the correlation setting is
utilized more than in Figure 3(b). While in GLS the realization of
$\boldsymbol{b}$ affects the fixed effects estimates (see Section 2.4.2), in
LMM the realization of $\boldsymbol{b}$ affects both, the fixed effects
estimates and the estimates of $u_{1},...,u_{I}.$

For demonstrating the model selection performance, Figure 4 presents the
agreement rates of $\widehat{CV_{c}}$ and CV with the approximated
generalization error, over the repeated simulation runs. The results are
presented for different sample sizes and for different predictive models (LMM
and GLS).

(a)

(b)

Figure 4: (4(a)): Each bar refers to the agreement rate of the relevant
criterion, $\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV_{c}}(h)$
and $\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;CV(h),$ with
$\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;generalization\;error(h)$ for
different sample sizes. Also, two scenarios are presented: when the variance
parameters are known (in dark color) and when they are estimated (in light
color). (4(b)): Similarly, the analysis is presented when GLS is implemented
instead of LMM.

As can be seen from Figures 4, when LMM is implemented
$\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV_{c}}(h)$ performs
better than $\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV}(h),$
however when GLS is implemented both criteria perform the same. It is
important to note that in order to have a high agreement rate with the oracle,
estimating unbiasedly the generalization error, as $\widehat{CV_{c}}$ does, is
not enough since the variance of the estimator can mix the ranks of the models
relative to the oracle. Therefore, there is no guarantee that
$\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV_{c}}(h)$ performs
better than $\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV}(h).$
This is also applies to other model selection criteria that are based on
comparing the models’ estimated prediction errors. However, in practice, as we
see here and in the real dataset examples in the next section,
$\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV_{c}}(h)$ performs
better than $\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV}(h)$
when the CV bias is large.

In addition, for LMM, both
$\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV_{c}}(h)$ and
$\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV}(h)$ with the true
variance parameters perform better than the versions with the estimated
variance parameters, while for GLS both perform similarly. This can be
explained by the key role of the variance in predicting using LMM, compared to
predicting using GLS. Also, when the sample size is small ($n=300$), it is
more difficult to estimate accurately the variance parameters, and the
difference between the two
$\underset{h\in\mathcal{H}}{\mathrm{argmin}}\;\widehat{CV_{c}}(h)$ versions is
substantial.

###  3.2 Real Data Analysis

This section presents analysis of two real datasets, a dataset with a
clustered correlation structure, and a dataset with hierarchical spatial
correlation structure.

####  3.2.1 Black Friday Dataset — Clustered Correlation Structure

The Black Friday dataset222Presented by Analytics Vidhya. contains information
of $737,577$ transactions made in a retail store on Black Friday by $5,891$
customers. The median number of transactions by a customer is $53,$ and the
range is $[5,1025].$ The dataset is available in https://github.com/
AssafRab/CVc.

A training sample of all transactions from $100$ random customers was drawn
and the goal is to fit a predictive model that predicts the purchase amount of
a new transaction of a new random customer. Therefore, while in CV setting
$\mathrm{Cov}(\boldsymbol{y}_{-k},y_{k})\neq 0,$ in the prediction goal
setting $\mathrm{Cov}(\boldsymbol{y}_{tr},y_{te})=0.$ Therefore the condition
of Theorem 2.1 is not satisfied $(P_{T_{te},T_{tr}}\neq P_{T_{k},T_{-k}})$ and
K-fold CV is unsuitable.

Three GLS models were fitted using different covariates (see Table 6) with the
covariance matrix
$\mathrm{Cov}(y_{i,j},y_{i^{\prime},j^{\prime}})=\sigma^{2}_{customer}1_{[i=i^{\prime}]}+\sigma^{2}_{\epsilon}1_{[i=i^{\prime},j=j^{\prime}]},$
where $i$ is customer index and $j$ is the index for an observation in a
customer’s set of observations. $\sigma^{2}_{customer}$ and
$\sigma^{2}_{\epsilon}$ were estimated using restricted maximum likelihood of
normal distribution.

CV and $\widehat{CV_{c}}$ were calculated using the training sample. Test
error was calculated by averaging the prediction error of observations of
other random $900$ customers using GLS models that were fitted by $T.$ This
procedure was repeated $50$ times with different random training and test
sets, and the generalization error was approximated by averaging the test
errors of the $50$ runs.

Figure 6 presents the means of CV, $\widehat{CV_{c}},$ and the test error over
the $50$ runs, as well as the confidence intervals of the means.

Figure 5: Prediction error estimation. For each model, the purple horizontal
line in the box is the approximated generalization error, which was calculated
by averaging the test error over the $50$ runs. The box range is the two
standard deviations confidence interval of the mean. Similarly, the mean
estimated prediction error and the two standard deviations confidence interval
of CV and $\widehat{CV_{c}}$ are presented.

Covariate\Model | Model 1 | Model 2 | Model 3
---|---|---|---
Intercept | $\checkmark$ | $\checkmark$ | $\checkmark$
Product Category 1 | $\checkmark$ | $\checkmark$ | $\checkmark$
Marital Status |  | $\checkmark$ | $\checkmark$
Gender |  | $\checkmark$ | $\checkmark$
Age |  | $\checkmark$ | $\checkmark$
Occupation |  |  | $\checkmark$
City Category |  |  | $\checkmark$
Stay in Current City Years |  |  | $\checkmark$

Figure 6: Model 1 contains two covariates, model 2 contains five covariates
and model 3 contains eight covariates.

As can be seen in Figure 6, $\widehat{CV_{c}}$ estimates the generalization
error better than CV for all the three models. Also, unlike CV,
$\widehat{CV_{c}}$ selects the same model as the oracle, Model 1, and follows
properly the prediction error estimation trend across the models. Analyzing
the business aspect of this example, leads to the surprising conclusion that
including only the product category as a fixed effect gives the best
prediction model, and adding covariates like age and gender does not improve
predictive performance, when the data are properly analyzed using
$\widehat{CV_{c}}.$

####  3.2.2 California Housing Dataset — Clustered Random Field Correlation
Structure

Another dataset that was analyzed is the California housing dataset, which
contains house values and other housing parameters in California. The dataset
has a hierarchical spatial correlation structure — each observation has
spatial coordinates value, where some of the observations share the same
coordinates value and therefore define a cluster. This correlation structure
is similar to the correlation structure in Sections 3.1, however with spatial
data rather than longitudinal data. The dataset is available in the python
scikit-learn package and the code is available in
https://github.com/AssafRab/CVc.

To emphasize the hierarchical clustering structure, all the clusters
containing only a single observation were excluded from the analysis ($8,237$
observations out of $20,640$), so the analyzed dataset contains $12,403$
observations from $4,353$ clusters. The median number of observations in a
cluster is $2$ and the range is $[2,15].$

A training sample, $T,$ of $700$ independent clusters was randomly drawn. The
prediction goal is to predict the ’median house value in a block’ of a new
observation from a different cluster than the clusters in $T.$ Three GPR
models with different covariates were fitted (see Table 8) with the following
covariance function:

$\mathrm{Cov}(y_{i},y_{i^{\prime}})=K_{exp}(\boldsymbol{z}_{i}-\boldsymbol{z}_{i^{\prime}})+\sigma^{2}\times
1_{[i=i^{\prime}]},$

where $\boldsymbol{z}_{i}\in\mathbb{R}^{2}$ contains the latitude and
longitude of observation $i$ and $K_{exp}(\cdot)$ is the exponential kernel
function. The parameters of $K_{exp}(\cdot)$ and $\sigma^{2}$ were fitted
using maximum likelihood of normal distribution.

Since the prediction goal is to predict a new observation from a different
cluster than the clusters in $T,$ then
$\mathrm{Cov}(\boldsymbol{y}_{-k},y_{k})\neq\mathrm{Cov}(\boldsymbol{y}_{-k},y_{te}).$
Therefore the condition of Theorem 2.1 is not satisfied and standard K-fold CV
is unsuitable. It is important to note that unlike the Black Friday example,
here $\mathrm{Cov}(\boldsymbol{y}_{-k},y_{te})\neq 0$ although $y_{te}$ is
from a new cluster. This is due to the spatial correlation structure of the
dataset, i.e., the latent random function $\boldsymbol{s}.$ The
$\widehat{CV_{c}}$ correction in this case is
$2/n\times\mathrm{tr}\big{(}H_{cv}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y}|\boldsymbol{s})\big{)},$
where $\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y}|\boldsymbol{s})$, which
expresses the extra correlation that the clusters contribute over the spatial
correlation of the data (see Section 2.4.1), is estimated as follows:

$\widehat{\mathrm{Cov}}(\boldsymbol{y},\boldsymbol{y}|\boldsymbol{s})[i,i^{\prime}]=\begin{cases}\begin{array}[]{c}K_{exp}(0)-\frac{1}{|j,j^{\prime}\;\forall\boldsymbol{z}_{j}\neq\boldsymbol{z}_{j^{\prime}}|}\sum_{j,j^{\prime}\;\forall\boldsymbol{z}_{j}\neq\boldsymbol{z}_{j^{\prime}}}K_{exp}(\boldsymbol{z}_{j}-\boldsymbol{z}_{j^{\prime}})\\\
0\end{array}&\begin{array}[]{c}\boldsymbol{z}_{i}=\boldsymbol{z}_{i^{\prime}}\\\
\boldsymbol{z}_{i}\neq\boldsymbol{z}_{i^{\prime}}\\\ \end{array}\end{cases}.$

The analysis was repeated $50$ times with different random $T,$ at each run CV
and $\widehat{CV_{c}}$ were calculated and the test error was calculated by
averaging the prediction error of the remaining observations of the $3,653$
clusters that are not in $T$ using GPR models that were fitted by $T.$
Averaging the test error over the $50$ runs approximates the generalization
error.

Figure 8 presents the means of $CV,$ and $\widehat{CV_{c}},$ and the test
error over the $50$ runs, as well as the confidence intervals of the means.

Figure 7: Prediction error estimation. For each model, the purple horizontal
line in the box is the approximated generalization error, which was calculated
by averaging the test error over the $50$ runs. The box range is the two
standard deviations confidence interval of the mean. Similarly, the mean
estimated prediction error and the two standard deviations confidence interval
of CV and $\widehat{CV_{c}}$ are presented. The unit of the dependent variable
is in $\$100,000.$

Covariate\Model | Model 1 | Model 2 | Model 3
---|---|---|---
Intercept | $\checkmark$ | $\checkmark$ | $\checkmark$
Median income in block | $\checkmark$ | $\checkmark$ | $\checkmark$
Median house age in block | $\checkmark$ | $\checkmark$ | $\checkmark$
Average number of rooms |  | $\checkmark$ | $\checkmark$
Average number of bedrooms |  | $\checkmark$ | $\checkmark$
Block population |  |  | $\checkmark$
Average house occupancy |  |  | $\checkmark$

Figure 8: Model 1 contains three covariates, model 2 contains five covariates
and model 3 contains seven covariates.

Similarly to Figure 6, $\widehat{CV_{c}}$ estimates the generalization error
better than CV for all the models, selects the same model as the oracle (Model
3), and follows properly the prediction error estimation trend across the
models. Also, while in Figure 6 CV underestimates the prediction error of the
saturated model more than the smaller models, here CV underestimates the
saturated model less than the smaller models. Therefore, biasedness of CV
prediction error can be expressed in model selection in different ways —
sometimes by favoring over-parameterized models and sometimes by favoring
under-parameterized models.

##  4 Conclusion and Discussion

In this paper we tackle the problem of applying CV as an estimate of
generalization error in non-i.i.d. situations. While the fundamental concerns
that this presents are widely acknowledged, a clear understanding of when
adjustments are needed, and what adjustments are appropriate, seems lacking in
much of the literature (Roberts et al. 2017; Saeb et al. 2017; Anderson et al.
2018).

We first present a general formulation of the bias in using CV in presence of
correlations, which leads to a clear general definition of settings where no
correction to CV is needed (Theorem 2.1). It shows that non-i.i.d. situations
can still facilitate correctness of regular CV, as long as the dependence
structure between training and prediction points is consistent in CV and
actual prediction. This simple result appears to contradict some previous
claims in the literature. An example can be found in Roberts et al. (2017),
where mechanisms of controlling folds partitioning are suggested for Kriging
tasks, based on the conception that CV is always over-optimistic for
non-i.i.d. data.

We then present a derivation of a bias correction for linear models under
general correlation structures (Theorem 2.2), which is used in establishing a
bias corrected cross-validation version, $\widehat{CV_{c}}$. To implement our
correction, it is necessary to specify the covariance structures within the
training set and between the training and prediction sets. This is typically
also required to choose a modeling approach. For example, in a simple linear
mixed model with normal assumptions, if one assumes that the random effect
realizations are the same when predicting, then LMM prediction is appropriate,
while if random effect realizations are new, then using GLS for estimation of
fixed effects only for prediction is more appropriate (Verbeke 1997). However,
it is important to emphasize that the validity of the correction
$\widehat{CV_{c}}$ does not depend on selection of an appropriate modeling
approach. In other words, if one mistakenly uses GLS where LMM is appropriate,
$\widehat{CV_{c}}$ still gives an unbiased estimate of generalization error
for the resulting model, as long as the covariance structure is correctly
specified.

In practice, the covariance matrices are typically not fully known, but
partially estimated from the data (for example, variance parameters in LMM can
be estimated using restricted maximum likelihood, Verbeke (1997)), and this is
also required for applying $\widehat{CV_{c}}$ to correct CV results. This
could potentially add uncertainty to the estimates, but as demonstrated in
Section 3, it does not tend to affect their expectation.

A fundamental assumption that is taken throughout this paper is that the
marginal distributions of the training and prediction sets are the same. In
case the marginal distributions are different, i.e., the observations in the
training set are drawn from a different marginal distribution than the
prediction data, then $\widehat{CV_{c}}$ is unsuitable. This scenario requires
further research and relevant applications are given in Sections 2.5.1 and
2.5.2. Other important use cases of different marginal distributions for
training and prediction sets are discussed in Rabinowicz and Rosset (2020) and
solutions for these use cases are proposed in context of transductive
prediction error (rather than generalization error approach taken here).

An interesting situation also covered by $\widehat{CV_{c}}$ and not often
discussed, is when the training set contains i.i.d. observations, but new data
points where predictions are made are actually correlated with the training
set (for example, new observations in the same set of cities). In this case,
if the correlation structure is known, then $\widehat{CV_{c}}$ can still be
used to estimate the prediction error.

## References

  * Altman (1990) Altman, N. S. (1990), “Kernel smoothing of data with correlated errors,” Journal of the American Statistical Association, 85, 749–759.
  * Anderson et al. (2018) Anderson, C., Hafen, R., Sofrygin, O., Ryan, L., and HBGDki Community (2018), “Comparing predictive abilities of longitudinal child growth models,” Statistics in Medicine, 38, 3555–3570.
  * Anderson and Darling (1952) Anderson, T. W. and Darling, D. A. (1952), “Asymptotic theory of certain ”goodness of fit” criteria based on stochastic processes,” The Annals of Mathematical Statistics, 23, 193–212.
  * Arlot et al. (2010) Arlot, S., Celisse, A., et al. (2010), “A survey of cross-validation procedures for model selection,” Statistics Surveys, 4, 40–79.
  * Breiman (2001) Breiman, L. (2001), “Random forests,” Machine Learning, 45, 5–32.
  * Breslow and Clayton (1993) Breslow, N. E. and Clayton, D. G. (1993), “Approximate inference in generalized linear mixed models,” Journal of the American statistical Association, 88, 9–25.
  * Burman (1989) Burman, P. (1989), “A comparative study of ordinary cross-validation, v-fold cross-validation and the repeated learning-testing methods,” Biometrika, 76, 503–514.
  * Burman et al. (1994) Burman, P., Chow, E., and Nolan, D. (1994), “A cross-validatory method for dependent data,” Biometrika, 81, 351–358.
  * Craven and Wahba (1978) Craven, P. and Wahba, G. (1978), “Smoothing noisy data with spline functions,” Numerische Mathematik, 31, 377–403.
  * Efron (1986) Efron, B. (1986), “How biased is the apparent error rate of a prediction rule?” Journal of the American statistical Association, 81, 461–470.
  * Fang (2011) Fang, Y. (2011), “Asymptotic equivalence between cross-validations and Akaike information criteria in mixed-effects models,” Journal of Data Science, 9, 15–21.
  * Francisco-Fernandez and Opsomer (2005) Francisco-Fernandez, M. and Opsomer, J. D. (2005), “Smoothing parameter selection methods for nonparametric regression with spatially correlated errors,” Canadian Journal of Statistics, 33, 279–295.
  * Friedman (2002) Friedman, J. H. (2002), “Stochastic gradient boosting,” Computational Statistics & Data Analysis, 38, 367–378.
  * Gelman (2006) Gelman, A. (2006), “Multilevel (hierarchical) modeling: what it can and cannot do,” Technometrics, 48, 432–435.
  * Goovaerts (1999) Goovaerts, P. (1999), “Geostatistics in soil science: state-of-the-art and perspectives,” Geoderma, 89, 1–45.
  * Greven and Kneib (2010) Greven, S. and Kneib, T. (2010), “On the behaviour of marginal and conditional AIC in linear mixed models,” Biometrika, 97, 773–789.
  * Hand (2017) Hand, D. J. (2017), Practical longitudinal data analysis, Routledge.
  * Hansen (2007) Hansen, C. B. (2007), “Generalized least squares inference in panel and multilevel models with serial correlation and fixed effects,” Journal of Econometrics, 140, 670–694.
  * Harville et al. (1976) Harville, D. et al. (1976), “Extension of the Gauss-Markov theorem to include the estimation of random effects,” The Annals of Statistics, 4, 384–395.
  * Hastie et al. (2009) Hastie, T., Tibshirani, R., and Friedman, J. (2009), The elements of statistical learning: data mining, inference, and prediction, New York, NY: Springer.
  * Hodges and Sargent (2001) Hodges, J. S. and Sargent, D. J. (2001), “Counting degrees of freedom in hierarchical and other richly-parameterised models,” Biometrika, 88, 367–379.
  * Jordan (2004) Jordan, M. (2004), “Graphical models,” Statistical Science, 19, 140–155.
  * Lee and Nelder (1996) Lee, Y. and Nelder, J. A. (1996), “Hierarchical generalized linear models,” Journal of the Royal Statistical Society: Series B (Methodological), 58, 619–656.
  * Liang et al. (2008) Liang, H., Wu, H., and Zou, G. (2008), “A note on conditional AIC for linear mixed-effects models,” Biometrika, 95, 773–778.
  * Little et al. (2017) Little, M. A., Varoquaux, G., Saeb, S., Lonini, L., Jayaraman, A., Mohr, D. C., and Kording, K. P. (2017), “Using and understanding cross-validation strategies. Perspectives on Saeb et al.” GigaScience, 6, 1–6.
  * Maddison (1990) Maddison, W. P. (1990), “A method for testing the correlated evolution of two binary characters: are gains or losses concentrated on certain branches of a phylogenetic tree?” Evolution, 44, 539–557.
  * Miller (1981) Miller, K. S. (1981), “On the inverse of the sum of matrices,” Mathematics Magazine, 54, 67–72.
  * Opsomer et al. (2001) Opsomer, J., Wang, Y., and Yang, Y. (2001), “Nonparametric regression with correlated errors,” Statistical Science, 16, 134–153.
  * Rabinowicz and Rosset (2020) Rabinowicz, A. and Rosset, S. (2020), “Assessing Prediction Error at Interpolation and Extrapolation Points,” Electronic Journal of Statistics, 14, 272–301.
  * Rasmussen and Williams (2006) Rasmussen, C. E. and Williams, C. K. (2006), Gaussian processes for machine learning, MIT press Cambridge.
  * Raudenbush and Bryk (2002) Raudenbush, S. W. and Bryk, A. S. (2002), Hierarchical linear models: Applications and data analysis methods, vol. 1, Sage.
  * Rice and Silverman (1991) Rice, J. A. and Silverman, B. W. (1991), “Estimating the mean and covariance structure nonparametrically when the data are curves,” Journal of the Royal Statistical Society: Series B (Methodological), 53, 233–243.
  * Roberts et al. (2017) Roberts, D. R., Bahn, V., Ciuti, S., Boyce, M. S., Elith, J., Guillera-Arroita, G., Hauenstein, S., Lahoz-Monfort, J. J., Schröder, B., Thuiller, W., et al. (2017), “Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure,” Ecography, 40, 913–929.
  * Saeb et al. (2017) Saeb, S., Lonini, L., Jayaraman, A., Mohr, D. C., and Kording, K. P. (2017), “The need to approximate the use-case in clinical machine learning,” GigaScience, 6, 1–9.
  * Stone (1974) Stone, M. (1974), “Cross-validatory choice and assessment of statistical predictions,” Journal of the Royal Statistical Society: Series B (Methodological), 36, 111–133.
  * Vaida and Blanchard (2005) Vaida, F. and Blanchard, S. (2005), “Conditional Akaike information for mixed-effects models,” Biometrika, 92, 351–370.
  * Verbeke (1997) Verbeke, G. (1997), “Linear mixed models for longitudinal data,” in Linear mixed models in practice, Springer, pp. 63–153.
  * Wang et al. (2016) Wang, J.-L., Chiou, J.-M., and Müller, H.-G. (2016), “Functional data analysis,” Annual Review of Statistics and Its Application, 3, 257–295.
  * Ward and Gleditsch (2018) Ward, M. D. and Gleditsch, K. S. (2018), Spatial regression models, vol. 155, Sage Publications.
  * Wu and Zhang (2002) Wu, H. and Zhang, J.-T. (2002), “Local polynomial mixed-effects models for longitudinal data,” Journal of the American Statistical Association, 97, 883–897.
  * Zhang and Yang (2015) Zhang, Y. and Yang, Y. (2015), “Cross-validation for selecting a model selection procedure,” Journal of Econometrics, 187, 95–112.