# Assessing Prediction Error at Interpolation and Extrapolation Points

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 Israeli Science Foundation, grant 1804/16

###### Abstract

Common model selection criteria, such as $AIC$ and its variants, are based on
in-sample prediction error estimators. However, in many applications involving
predicting at interpolation and extrapolation points, in-sample error cannot
be used for estimating the prediction error. In this paper new prediction
error estimators, $tAI$ and $Loss(w_{t})$ are introduced. These estimators
generalize previous error estimators, however are also applicable for
assessing prediction error in cases involving interpolation and extrapolation.
Based on the prediction error estimators, two model selection criteria with
the same spirit as $AIC$ are suggested. The advantages of our suggested
methods are demonstrated in simulation and real data analysis of studies
involving interpolation and extrapolation in a linear mixed model framework.

Keywords: Model selection; $AIC$; Linear Mixed Model; Kriging

##  1 Introduction

Predicting a phenomenon at different points than the points appearing in the
training sample plays an important role across many research fields such as in
Geostatistics (Li and Heap 2014; Kyriakidis and Journel 1999), health (Manton
et al. 2012) and Econometrics (Baltagi 2008). In many of these use cases the
new predicted points are interpolation or extrapolation points with respect to
space or to time. For example, Brown and Comrie (2002) interpolated climate
values in Southwestern U.S., where the coverage of climate information is
sparse. By predicting at interpolation points, they created a high-resolution
map of seasonal temperature and precipitation in this area. Another example
given by Stewart et al. (2009) is forecasting the effects of obesity and
smoking on U.S. life expectancy in 2020 by using a data set for the years 2003
through 2006.

Modeling approaches involving prediction at interpolation and extrapolation
points were studied in Machine Learning, mainly in the context of transductive
Support Vector Machine (Joachims 1999), however also in regression (Le et al.
2006).

Assessing prediction error at interpolation and extrapolation points, or more
generally at transduction points, cannot be done using traditional in-sample
prediction error estimators as is used in $AIC$ (Akaike 1974) and its
variants. Similarly, K-fold Cross-Validation, which estimates the
generalization error, is also unsuitable in these cases, where prediction
points are specified.

This paper introduces a prediction error estimator, $tAI,$ which generalizes
previous in-sample prediction error estimators like $mAI$ (Vaida and Blanchard
2005) and $cAI$ (Vaida and Blanchard 2005), however, it doesn’t assume that
the predicted points are the same as the points appearing in the training
sample and therefore is applicable to a wider range of use cases, such as
cases involving prediction at interpolation and extrapolation points. Since
prediction error assessment is highly related to model selection, a new model
selection criterion, $tAIC,$ which is based on $tAI,$ is proposed as well.
$tAI$ is suitable when the observations are normally distributed, whether they
are correlated or not and therefore is applicable for various parametric
models with different variance structure assumptions such as Linear Mixed
Model (LMM), Gaussian Process Regression (GPR), Generalized Least Squares
(GLS) and Linear Regression. Relaxing the normality requirement of $tAI$, we
also propose in Section 5 an approach for inference on interpolation and
extrapolation that is based on squared error loss rather than likelihood, and
hence generalizes the Optimism approach in model selection (Efron 1986).

In many use cases involving predicting at interpolation and extrapolation
points, the dependent variable has a correlation structure (Li and Heap 2014;
Kyriakidis and Journel 1999). For example, in the use case that is given by
Brown and Comrie (2002), it is natural to assume a spatial correlation
structure on the Southwestern U.S. area. Similarly, in repeated measures
studies that forecast long-term treatment effects, a correlation structure
with respect to time is commonly assumed (Ho et al. 2011). Therefore, use
cases involving correlated data and models that are implemented on correlated
data, such as LMM, GPR and GLS, are good platforms for analyzing how
predicting at interpolation and extrapolation points influences prediction
error estimation and model selection. Before introducing $tAI,$ a setup which
puts LMM, GPR and GLS under a unified framework, will be defined:

Let $\boldsymbol{y}\in\mathbb{R}^{n}$ and the fixed matrices
$\\{X\in\mathbb{R}^{n\times p},\,Z\in\mathbb{R}^{n\times q}\\}$ be a training
sample, $\boldsymbol{y^{*}}\in\mathbb{R}^{n^{*}}$ and the fixed matrices
$\\{X^{*}\in\mathbb{R}^{n^{*}\times p},\,Z^{*}\in\mathbb{R}^{n^{*}\times
q}\\}$ be a prediction set, where

$\displaystyle\boldsymbol{y}\sim
N(\boldsymbol{\mu},V),\,\boldsymbol{y^{*}}\sim N(\boldsymbol{\mu^{*}},V^{*}),$
(1)

$\boldsymbol{\mu}=X\boldsymbol{\beta},\,\boldsymbol{\mu^{*}}=X^{*}\boldsymbol{\beta},$
$V$ is a function of $Z$ and $V^{*}$ is a function of $Z^{*}.$ For example, in
LMM it is typically assumed that the columns of $Z,\,Z^{*}$ are associated
with normally distributed random effects with covariance matrix
$G\in\mathbb{R}^{q\times q}$ such that

$\displaystyle V$
$\displaystyle=ZGZ^{{}^{\prime}}+\sigma^{2}I_{n},\,\,\sigma\in\mathbb{R}^{+}$
$\displaystyle V^{*}$
$\displaystyle=Z^{*}GZ^{*^{\prime}}+\sigma^{2}I_{n^{*}},$

where $I_{n},\,I_{n}^{*}$ are the identity matrices with dimensions $n$ and
$n^{*}$ respectively. In GPR it is often assumed that

$\displaystyle V$ $\displaystyle=K(Z,Z)+\sigma^{2}I_{n}$ $\displaystyle V^{*}$
$\displaystyle=K(Z^{*},Z^{*})+\sigma^{2}I_{n}^{*}$

where $K$ is some kernel function.

In addition denote

$\displaystyle R^{*}$
$\displaystyle=\mathrm{Var}(\boldsymbol{y^{*}}|\boldsymbol{y})$ $\displaystyle
R$ $\displaystyle=\mathrm{Var}(\boldsymbol{y^{new}}|\boldsymbol{y}),$

where $\boldsymbol{y^{new}}\sim N(\boldsymbol{\mu},V)$ is an i.i.d copy of
$\boldsymbol{y}.$

By normality of $\boldsymbol{y}$ and $\boldsymbol{y^{*}},$

$\mathbb{E}(\boldsymbol{y^{*}}|\boldsymbol{y})=X^{*}\boldsymbol{\beta}+\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})V^{-1}(\boldsymbol{y}-X\boldsymbol{\beta}).$

Given $V,\,\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})$ and the ML
estimator of $\boldsymbol{\beta},$

$\boldsymbol{\hat{\beta}}=(X^{{}^{\prime}}V^{-1}X)^{-1}X^{{}^{\prime}}V^{-1}\boldsymbol{y},$

$\mathbb{E}(\boldsymbol{y^{*}}|\boldsymbol{y})$ can be used for predicting
$\boldsymbol{y^{*}}$ as follows

$\displaystyle\boldsymbol{\hat{f}^{*}}$
$\displaystyle=\hat{\mathbb{E}}(\boldsymbol{y^{*}}|\boldsymbol{y})$ (2)
$\displaystyle=X^{*}(X^{{}^{\prime}}V^{-1}X)^{-1}X^{{}^{\prime}}V^{-1}\boldsymbol{y}+\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})V^{-1}\left\\{I_{n}-X(X^{{}^{\prime}}V^{-1}X)^{-1}X^{{}^{\prime}}V^{-1}\right\\}\boldsymbol{y}.$

This procedure generalizes standard prediction practices in LMM, GPR and GLS.
In addition, $\boldsymbol{\hat{f}^{*}}$ is the Best Linear Unbiased Predictor
(BLUP) (Harville 1976).

$tAI$ is an estimator of the following prediction error,

$\displaystyle-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})=-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}\log\left[\frac{\exp{\left\\{\frac{-1}{2}(\boldsymbol{y^{*}}-\boldsymbol{\hat{f}^{*}})^{{}^{\prime}}R^{*^{-1}}(\boldsymbol{y^{*}}-\boldsymbol{\hat{f}^{*}})\right\\}}}{\sqrt{(2\pi)^{n^{*}}|R^{*}|}}\right].$
(3)

Correspondingly, given a set of candidate models, $tAIC$ would be defined as a
model selection criterion selecting a model with the minimal $tAI.$ This
methodology of estimating the prediction errors for different models and then
selecting the model with the minimal prediction error, is the same as is
implemented in $AIC$ and its variants.

$\\{X^{*},Z^{*},R^{*}\\}=\\{X,Z,R\\}$ is not assumed in the setup above and in
its associated prediction error measure, eq. (3). Therefore, $tAI$ is
applicable in various use cases that require flexibility in defining
$\\{X^{*},Z^{*},R^{*}\\}.$ For example, in the use case mentioned above of
Brown and Comrie (2002), where GPR is used for predicting interpolated climate
values (Kriging), it is reasonable to define $\\{X^{*},Z^{*}\\}$ as the data
points at the high-resolution spatial array rather than as the data points at
the training sample, $\\{X,Z\\},$ which cover the area sparsely. Therefore,
while prediction error estimators that are based on in-sample error estimation
and generalization error estimation are unsuitable to this case, $tAI$ is
suitable. For similar considerations, $tAIC$ is required in repeated measures
studies in health and Biomedicine, when the main interest is to select LMM
model minimizing the prediction error at long-term points,
$\\{X^{*},Z^{*},R^{*}\\}$, which are different than the points that are used
for model building, $\\{X,Z,R\\}$ (Pope III et al. 2002; Li et al. 2008).

Beside downscaling of climate maps and estimating long-term effect in clinical
studies, interpolation and extrapolation using LMM and Kriging are important
tools for many research topics in mining engineering, agriculture,
environmental sciences, especially when sampling is difficult and expensive
like in mountainous and deep marine regions (Li and Heap 2011; Stahl et al.
2006; Vicente-Serrano et al. 2003). $tAI$ and $tAIC$ are relevant for all
these research topics as well as for others which don’t involve interpolation
and extrapolation but still don’t satisfy
$\\{X^{*},Z^{*},R^{*}\\}=\\{X,Z,R\\}.$ Various use cases will be presented and
analyzed in Sections 3 and 4.

##  2 $tAI$ and $tAIC$

$tAI$ is derived by estimating
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}$
by the averaged log-likelihood of the training sample,

$-\frac{1}{n}l(\boldsymbol{y})=-\frac{1}{n}\log\left[\frac{\exp{\left\\{\frac{-1}{2}(\boldsymbol{y}-\boldsymbol{\hat{f}})^{{}^{\prime}}R^{-1}(\boldsymbol{y}-\boldsymbol{\hat{f}})\right\\}}}{\sqrt{(2\pi)^{n}|R|}}\right],$

plus a penalty correction

$C_{tAI}=\mathbb{E}_{\boldsymbol{y}}\left[-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})-\left\\{-\frac{1}{n}l(\boldsymbol{y})\right\\}\right],$

where

$\boldsymbol{\hat{f}}=X\boldsymbol{\hat{\beta}}+(V-R)V^{-1}(\boldsymbol{y}-X\boldsymbol{\hat{\beta}})$

is the estimated conditional expectation,
$\hat{\mathbb{E}}(\boldsymbol{y^{*}}|\boldsymbol{y}),$ when
$\\{X^{*},Z^{*},R^{*}\\}=\\{X,Z,R\\}.$ This approach of estimating prediction
error by deriving the bias of the training error is also used in $AIC$ and its
variants (Akaike 1974). Consequently, the estimator

$tAI=-\frac{1}{n}l(\boldsymbol{y})+C_{tAI}$

doesn’t contain $\boldsymbol{y^{*}}$ but still satisfies

$E_{\boldsymbol{y}}tAI=\mathbb{E}_{\boldsymbol{y}}\left\\{-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})\right\\},$

and $tAI$ can be seen either as an estimator of
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}$
or of its expectation
$-\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}=-\mathbb{E}_{\boldsymbol{y},\boldsymbol{y^{*}}}l(\boldsymbol{y^{*}})/n^{*}.$

The following theorem and corollary introduce a general expression for
$C_{tAI}$ and therefore also for $tAI$ and $tAIC$.

######  Theorem 1.

Consider the setup given in eq. (1). In addition, let
$H\boldsymbol{y}\in\mathbb{R}^{n}$ and
$H^{*}\boldsymbol{y}\in\mathbb{R}^{n^{*}}$ be predictors of $\boldsymbol{y}$
and $\boldsymbol{y^{*}}$ respectively, where $H$ and $H^{*}$ don’t contain
$\boldsymbol{y},\,\boldsymbol{y^{*}}$ and satisfy
$H\boldsymbol{\mu}=\boldsymbol{\mu},\,H^{*}\boldsymbol{\mu}=\boldsymbol{\mu^{*}}.$
Then

$\displaystyle\mathbb{E}_{\boldsymbol{y}}\left\\{-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})+\frac{1}{n}l(\boldsymbol{y})\right\\}=$
$\displaystyle\frac{1}{n}tr\left(R^{-1}HV\right)-\frac{1}{n^{*}}tr\left(R^{*^{-1}}H^{*}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})\right)$
$\displaystyle+\frac{1}{2}\left\\{\log\left(\frac{|R^{*}|^{\frac{1}{n^{*}}}}{|R|^{\frac{1}{n}}}\right)+\frac{1}{n^{*}}tr\left(R^{*^{-1}}V^{*}\right)-\frac{1}{n}tr\left(R^{-1}V\right)\right\\}$
$\displaystyle+\frac{1}{2}\left\\{\frac{1}{n^{*}}tr\left(R^{*^{-1}}H^{*}VH^{*^{\prime}}\right)-\frac{1}{n}tr\left(R^{-1}HVH^{{}^{\prime}}\right)\right\\}.$

The proof is attached in Appendix A.

######  Corollary 2.

Under the set-up described in eq. (1), the conditions in Theorem 1 are
satisfied by the BLUPs $\boldsymbol{\hat{f}}$ and $\boldsymbol{\hat{f}^{*}},$
where

$\displaystyle\boldsymbol{\hat{f}}$ $\displaystyle=H\boldsymbol{y}$
$\displaystyle H$
$\displaystyle=X(X^{{}^{\prime}}V^{-1}X)^{-1}X^{{}^{\prime}}V^{-1}+(V-R)V^{-1}\left\\{I_{n}-X(X^{{}^{\prime}}V^{-1}X)^{-1}X^{{}^{\prime}}V^{-1}\right\\},$

and

$\displaystyle\boldsymbol{\hat{f}^{*}}$ $\displaystyle=H^{*}\boldsymbol{y}$
$\displaystyle H^{*}$
$\displaystyle=X^{*}(X^{{}^{\prime}}V^{-1}X)^{-1}X^{{}^{\prime}}V^{-1}+\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})V^{-1}\left\\{I_{n}-X(X^{{}^{\prime}}V^{-1}X)^{-1}X^{{}^{\prime}}V^{-1}\right\\}.$

By Corollary 2 and Theorem 1, $C_{tAI}$ can be calculated under the setup that
is described in eq. (1). Therefore, $tAI$ can be implemented in LMM, GPR and
other related models.

Beside prediction error estimation, these results can be used for defining the
following model selection criterion.

######  Definition 3.

Given set of models $\mathcal{H},$ satisfying the conditions in Theorem 1,
$tAIC$ is the following criterion

$h_{best}=\underset{h\in\mathcal{H}}{\mathrm{argmin}}\,tAI_{h},$

where $tAI_{h}$ is $tAI$ for model $h.$

###  2.1 Comparison with other prediction error estimators

The prediction error estimators that appear in $cAIC$ and $mAIC$ (Vaida and
Blanchard 2005) were developed for normal linear models under different
restrictions on the variance structure, but assuming
$\\{X,Z,R\\}=\\{X^{*},Z^{*},R^{*}\\}$. $cAIC$ is aimed at the LMM and GPR
case, where $\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})\neq 0$, while
$mAIC$ considers the GLS case where
$\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})=0$. For $cAIC$ the prediction
error estimate is:

$cAI=-\frac{1}{n}l(\boldsymbol{y})+\frac{tr(H)}{n},\,\,\text{}$

while for $mAIC$ it is:

$mAI=-\frac{1}{n}l(\boldsymbol{y})+\frac{p}{n}.$

$cAIC,\;mAIC$ are defined from $cAI,\;mAI$ similarly to $tAIC.$

It is easy to confirm that when $\\{X,Z,R\\}=\\{X^{*},Z^{*},R^{*}\\}$ our
$tAI$ formula indeed reduce to the $cAI$ and $mAI$ formulas.

In addition, for GLS, we can also show an interesting interpretation for the
difference between the $mAI$ and $tAI$ expressions. With a little algebra we
get:

$\displaystyle C_{tAI}(GLS)=$
$\displaystyle\frac{p}{n}+\frac{1}{2}\log\left(\frac{|V^{*}|^{\frac{1}{n^{*}}}}{|V|^{\frac{1}{n}}}\right)$
$\displaystyle+\frac{1}{2}tr\left\\{(X^{{}^{\prime}}V^{-1}X)^{-1}\left(\frac{1}{n^{*}}X^{*^{\prime}}V^{*^{-1}}X^{*}-\frac{1}{n}X^{{}^{\prime}}V^{-1}X\right)\right\\}$
$\displaystyle=$
$\displaystyle\frac{p}{n}+\frac{1}{2}\log\left(\frac{|V^{*}|^{\frac{1}{n^{*}}}}{|V|^{\frac{1}{n}}}\right)+\frac{1}{2}tr\left[\mathrm{Var}(\boldsymbol{\hat{\beta}})\left\\{\frac{1}{n^{*}}\mathrm{Var}(\boldsymbol{\hat{\beta}^{*}})^{-1}-\frac{1}{n}\mathrm{Var}(\boldsymbol{\hat{\beta}})^{-1}\right\\}\right],$

where

$\boldsymbol{\hat{\beta}^{*}}=(X^{{}^{\prime}*}V^{*^{-1}}X^{*})^{-1}X^{{}^{\prime}*}V^{*^{-1}}\boldsymbol{y^{*}}.$

Since $\mathrm{Var}(\boldsymbol{\hat{\beta}})$ achieves the Cramer-Rao bound:

$\displaystyle C_{tAI}(GLS)$
$\displaystyle=\frac{p}{n}+\frac{1}{2}\log\left(\frac{|V^{*}|^{\frac{1}{n^{*}}}}{|V|^{\frac{1}{n}}}\right)+\frac{1}{2}tr\left[\mathcal{I}(\boldsymbol{\hat{\beta}})^{-1}\left\\{\frac{1}{n^{*}}\mathcal{I}(\boldsymbol{\hat{\beta}^{*}})-\frac{1}{n}\mathcal{I}(\boldsymbol{\hat{\beta}})\right\\}\right],$
(4)

where $\mathcal{I}$ is Fisher-information. The determinants $|V|$ and
$|V^{*}|$ are often called the generalized variance (Wilks 1932; Johnson et
al. 2014).

###  2.2 Relaxing Theorem 1 Conditions

Although this paper focuses on prediction error estimation and model selection
for LMM and GPR, Theorem 1 is more general and doesn’t assume the paradigm
applied in LMM and GPR, i.e., predicting using
$E(\boldsymbol{y^{*}}|\boldsymbol{y})$ and estimating the marginal mean
parameters with MLE. Theorem 1 assumes:

  1. 1.

Normality of $\boldsymbol{y^{*}}$ and $\boldsymbol{y}.$

  2. 2.

$\mathbb{E}\boldsymbol{y}=X\boldsymbol{\beta},\,\mathbb{E}\boldsymbol{y^{*}}=X^{*}\boldsymbol{\beta}$

  3. 3.

$H\boldsymbol{\mu}=\boldsymbol{\mu},\,H^{*}\boldsymbol{\mu}=\boldsymbol{\mu^{*}}$

and therefore can be used in other cases satisfying the above conditions.

When the normality assumption cannot be taken, another model selection
criterion, which is based on similar approach as $tAI$ can be implemented. For
more details see Section 5.

In case the normality assumption can be taken, however the fitted model
doesn’t satisfy condition 3 of unbiasedness, the following extended version of
Theorem 1 results can be used instead:

$\displaystyle\mathbb{E}_{\boldsymbol{y}}\left\\{-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})+\frac{1}{n}l(\boldsymbol{y})\right\\}$
$\displaystyle\qquad\qquad=\frac{1}{n}tr\left(R^{-1}HV\right)-\frac{1}{n^{*}}tr\left(R^{*^{-1}}H^{*}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})\right)$
$\displaystyle\qquad\qquad\,\,\,\,\,\,\,+\frac{1}{2}\left\\{\log\left(\frac{|R^{*}|^{\frac{1}{n^{*}}}}{|R|^{\frac{1}{n}}}\right)+\frac{1}{n^{*}}tr\left(R^{*^{-1}}V^{*}\right)-\frac{1}{n}tr\left(R^{-1}V\right)\right\\}$
$\displaystyle\qquad\qquad\,\,\,\,\,\,\,+\frac{1}{2}\left\\{\frac{1}{n^{*}}tr\left(R^{*^{-1}}H^{*}VH^{*^{\prime}}\right)-\frac{1}{n}tr\left(R^{-1}HVH^{{}^{\prime}}\right)\right\\}$
$\displaystyle\qquad\qquad\,\,\,\,\,\,\,+\frac{1}{2n}tr\left(R^{-1}(2H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}-\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}-H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{{}^{\prime}})\right)$
$\displaystyle\qquad\qquad\,\,\,\,\,\,\,\,-\frac{1}{2n^{*}}tr\left(R^{*^{-1}}(2H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{*^{\prime}}}-\boldsymbol{\mu^{*}}\boldsymbol{\mu^{*^{\prime}}}-H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{*^{\prime}})\right).$

The proof can be found in Appendix A as part of the proof of Theorem 1.

Note that this expression is less useful, as it depends on $\boldsymbol{\mu}$
and $\boldsymbol{\mu^{*}}$ which are unknown.

##  3 Use cases

In this section, typical use cases of using $tAI$ and $tAIC$ are presented.

###  3.1 Predicting interpolation and extrapolation in spatial array and
longitudinal temporal data

As was described in the introduction, predicting interpolated and extrapolated
data points using LMM and GPR is common in Biomedicine, health, Climatology
and other research fields, where temporal and spatial datasets are common. The
flexible definition of $X^{*},Z^{*},R^{*}$ and $V^{*}$ in $tAI$ makes it
applicable when the goal is to estimate prediction error at interpolated and
extrapolated data points along time and space dimensions.

In Section 4 we analyze numerically a repeated measures clinical study,
containing child growth measurements (Potthoff and Roy 1964), where
interpolation and extrapolation objectives can be defined and application of
$tAI$ is demonstrated. The following example, built on the application of
Tsanas et al. (2010), demonstrates that appropriate use of $tAIC$ can also
simplify and improve on existing methodology.

######  Example 3.1.

Tsanas et al. introduced a new method for measuring progression of Parkinson’s
disease. Their motivation is that the standard methodology for measuring
Parkinson progression, which uses UPDRS score (Unified Parkinson’s Disease
Rating Scale), is costly and requires a physician visit. Their alternative
methodology is creating a formula that approximates the UPDRS score with
speech signals which are not costly. Six months data was collected for their
study, containing large amount of longitudinal speech signal measurements per
patient, however, UPDRS scores were collected only at a small number of the
time points. In order to select the best covariates with respect to the whole
speech signals sample, they suggested to interpolate the UPDRS scores using
’straightforward linear interpolation’, then to fit several alternative models
and to select one of them using $AIC$ and other model selection criteria. An
alternative paradigm that doesn’t require imputing UPDRS score is by using
$tAIC.$ Since $tAIC$ doesn’t assume $\\{X^{*},Z^{*}\\}=\\{X,Z\\},$ there is no
need in interpolating the UPDRS score in order to select a model minimizing
the estimated prediction error with respect to the whole speech signals
sample.

We note that in Example 3.1, one may think that $\boldsymbol{y^{*}}$ is used
twice, for model building and for prediction error estimation, and therefore
over-fitting can occur. However, since in $tAI$ approach, unlike in cross-
validation approach, $\boldsymbol{y^{*}}$ is used as a conceptual idea in
order to derive $C_{tAI}$ and not as real observations, $\boldsymbol{y^{*}}$
is not used twice.

In the spatial data analysis domain, common application areas include
geographical data (Li and Heap 2014) and neuroimaging data (Salimi-Khorshidi
et al. 2011). Such studies usually use GPR rather than LMM. Although GPR and
LMM reflect different perspectives — while GPR is based on functional data
analysis, LMM is based on multivariate analysis — and use different techniques
for expressing the covariance matrices, both models use conditional
expectation, $\boldsymbol{\hat{f}^{*}},$ for prediction, hence $tAI$ is also
applicable for GPR. In the Introduction we demonstrated this by the use case
of creating high-resolution climate maps (Brown and Comrie 2002). A similar
use case, containing chemical concentration in soil data is analyzed
numerically in Section 4.

###  3.2 Other Transductive Settings

LMM and GPR are also used for modeling data without spatial or temporal
correlation structure, and the prediction problems that arise often involve
prediction outside the training sample.

One interesting example is modeling the effect of SNPs (Single Nucleotide
Polymorphism) on a phenotype as part of a Genome-Wide Association Study
(GWAS). In this case the common practice is to consider the SNPs as random
effects and other explanatory variables (e.g. age, height and gender) as fixed
effects (Zhang et al. 2010). When using LMM for modeling the effect of SNPs on
phenotype, $tAI$ allows estimating the prediction error for an extended
population compared to the training sample. It is directly useful in the
important case when $\\{X^{*},Z^{*}\\}$ can be collected from other studies
which investigate different phenotype, however contain the SNPs and the
explanatory variables that are used in the training sample (Wray et al. 2013).

Missing values of the dependent variable which is a common phenomenon in
statistical analysis and in particular in clinical trial with repeated
measures study design (Wood et al. 2004; O’neill and Temple 2012). There are
many methods for handling missing values in repeated measures studies, some of
the methods involving missing values imputation (Mallinckrodt et al. 2003). In
case of having missing data of the dependent variable at some known points but
the goal is to estimate the prediction error with respect to the original
study design (Hogan et al. 2004), $tAI$ can be used without imputing the
missing values.

##  4 Numerical Results

This section focuses on comparison between $tAI,\,cAI$ and $mAI,$ as well as
between their corresponding model selection criteria, $tAIC,\,cAIC$ and
$mAIC,$ using simulation and real data analyses.

###  4.1 Simulation Analyses

The goal of the following analyses is to investigate the accuracy of $tAI,cAI$
and $mAI$ in estimating
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*},$
for different sample sizes and variance setups. In addition, $tAIC,cAIC$ and
$mAIC$ will also be analyzed and compared with respect to the oracle solution

$\displaystyle h_{best}$
$\displaystyle=\underset{h\in\mathcal{H}}{\mathrm{argmin}}-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l_{h}(\boldsymbol{y^{*}}).$

Additional numerical results with respect to a potentially different oracle
solution

$\displaystyle h_{best}$
$\displaystyle=\underset{h\in\mathcal{H}}{\mathrm{argmin}}-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l_{h}(\boldsymbol{y^{*}})$

are presented in Appendix C.

##### Simulation setup

The simulation demonstrates prediction error estimation and model selection
for the following LMM setting:

$\displaystyle\phi_{i,j}=0.5\times
time_{i,j}+\sum_{k=0}^{k=2}x_{i,j,k}+2\times\sum_{k=3}^{k=6}x_{i,j,k}+b_{i,1}+time_{i,j}\times
b_{i,2}+\epsilon_{i,j},$

where $i\in\\{1,...,S\\}$ is the subject number and $j\in\\{1,...,12\\}$ is
the measurement number. $b_{i,1}$ is distributed $N(0,15),$ $b_{i,2}$ is
distributed $N(0,1)$ and $\epsilon_{i,j}$ is distributed $N(0,\sigma^{2}).$ In
addition, $time_{i,j}=j,\,\forall j\leq 10,$
$time_{i,11}=15,\,time_{i,12}=20,$ $x_{i,j,0}=1,$ $x_{i,j,1}$ was drawn from
$Ber(0.5)$ and $x_{i,j,k},\,\forall k\geq 2$ was drawn form $N(0,1).$

The dependent variable in the training set $,\boldsymbol{y},$ was defined as
$\phi_{i,j}\,\forall j\in\\{1,..,10\\},$ the dependent variable in the
prediction set, $\boldsymbol{y^{*}},$ was defined as $\phi_{i,j}\,\forall
j\in\\{11,12\\}.$ Therefore, this setting demonstrates predicting at
extrapolation time points. This setting was generated nine times, for
different number of subjects, $S\in\\{10,20,30\\},$ and different residual
variance values, $\sigma^{2}\in\\{15,20,25\\}.$ Each simulation was repeated
$200$ times.

Three linear mixed models were fitted given the true covariance matrices, all
the models contain the time covariate, in addition, model number 1 contains
$x_{i,j,k},\,\forall k\leq 2,$ model number 2 contains $x_{i,j,k},\,\forall
k\leq 4$ and model number 3 contains $x_{i,j,k},\,\forall k$ which is also the
model that generates the data.

##### Results

Figure 1 presents the densities of $tAI,\,cAI,\,mAI$ and
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}$
as a function of the sample size and $\sigma^{2}$ for model number $3,$ as
generated from the $200$ simulation runs.

Figure 1: Densities of $tAI,\,cAI,\,mAI$ and
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}$
as a function of the number of subjects and $\sigma^{2}$.

As can be seen from Figure 1, $tAI$ density is concentrated around the mean of
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}.$
$cAI$ and $mAI$ are stochastically smaller than
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}$
since their corrections, $tr(H)/n$ and $p/n,$ are unsuitable for this case of
predicting at extrapolation points.

In addition, since $tAI,\,cAI$ and $mAI$ share the same random part,
$l(\boldsymbol{y})/n,$ but different mean,
$-\mathbb{E}_{\boldsymbol{y}}l(\boldsymbol{y})/n$ plus
$C_{tAI},\,tr(H)/n,\,p/n$ respectively, their densities have the same shape
however shifted with respect to the corrections. In contrast,
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}$
has the same mean as $tAI$ but different variance, since
$\mathrm{Var}\left(-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}\right)$
depends on $H^{*},\,R^{*}$ and $n^{*}$ that do not appear in
$\mathrm{Var}\left(tAI\right)=\mathrm{Var}\left(-l(\boldsymbol{y})/n\right).$
In our case, $H^{*}$ contains large values compared to $H$ and therefore

$\mathrm{Var}\left(\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}-\frac{1}{n^{*}}l(\boldsymbol{y^{*}})\right)>\mathrm{Var}\left(tAI\right).$

Figure 2 presents for each criterion, $tAIC,\,cAIC$ and $mAIC,$ the error

$\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}-\frac{1}{n^{*}}l_{h_{best}}(\boldsymbol{y^{*}}),$

where $h_{best}$ is the selected model by the relevant criterion. This error
reflects the true average error that is obtained when implementing the
different model selection criteria. In addition, the average error of the
oracle criterion,

$\displaystyle h_{best}$
$\displaystyle=\underset{h\in\\{1,2,3\\}}{\mathrm{argmin}}-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l_{h}(\boldsymbol{y^{*}}),$

is presented as well.

Figure 2: For each setup, each symbol refers to the prediction error
$\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}-l_{h_{best}}(\boldsymbol{y^{*}})/n^{*}$
of the relevant criterion, $mAIC,\,cAIC\,tAIC$ and the oracle criterion.

As can be seen from Figure 2, $tAIC$ obtain better results than $cAIC$ and
$mAIC$ in all the nine setups. Similar analysis with respect to the error

$\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}-\frac{1}{n^{*}}l_{h_{best}}(\boldsymbol{y^{*}})$

is presented in in Appendix C.

Figure 3 presents the agreement rate of the criteria, $tAIC,\,cAIC$ and $mAIC$
with the oracle criterion

$\displaystyle h_{best}$
$\displaystyle=\underset{h\in\\{1,2,3\\}}{\mathrm{argmin}}-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l_{h}(\boldsymbol{y^{*}}).$

Figure 3: For each setup, each bar refers to the agreement rate of the
relevant criterion with the oracle criterion

As can be seen from Figure 3, $tAIC$ achieves the best results in this case as
well.

Similar analysis with respect to the oracle criterion

$\displaystyle h_{best}$
$\displaystyle=\underset{h\in\\{1,2,3\\}}{\mathrm{argmin}}-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l_{h}(\boldsymbol{y^{*}})$

is presented in Appendix C.

###  4.2 Real data analyses

The analyses below focus on comparison between $tAI,\,cAI,\,mAI$ and

$-\frac{1}{n^{*}}l(\boldsymbol{y^{*}}).$

Here, $-l(\boldsymbol{y^{*}})/n^{*}$ is used as a ground truth instead of
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}$
since the latter is unknown for the real data sets.

####  4.2.1 Meuse data

##### Data description

Meuse data set was introduced by Rikken and Van Rijn (1993) and is available
in R software. The data was collected in a floodplain area of the river Meuse,
near the village of Stein, Netherlands, and contains 155 measurements of
topsoil concentrations of Zinc, Lead, Copper and Cadmium, along with location
(latitude and longitude) and other covariates. In addition, another data set,
Meuse.grid, is analyzed. Meuse.grid is a higher resolution grid of the same
area, containing 3103 observations of location and some of the covariates that
are available in the Meuse data set, however it doesn’t contain the metal
concentration measurements. The Meuse.grid is available in R software as well.

##### Results

The Meuse data set was partitioned randomly into training and test samples.
Four Gaussian process regression models were fitted to the log of the Lead
concentration.222Only $\log(Lead)$ can be analyzed under the normality
assumption. All the models share the same kernel structure, squared-
exponential kernel,

$K(\boldsymbol{Z_{i}},\boldsymbol{Z_{j}})=\sigma_{f}^{2}\exp\left[-\frac{1}{2}\left\\{\frac{1}{l_{1}^{2}}\left(Z_{i,1}-Z_{j,1}\right)^{2}+\frac{1}{l_{2}^{2}}\left(Z_{i,2}-Z_{j,2}\right)^{2}\right\\}\right],$

where $Z_{i,1}$ refers to the latitude of measurement $i,$ $Z_{i,2}$ refers to
longitude of measurement $i$ and $l_{1},\,l_{2}$ and $\sigma_{f}$ lie in
$\mathbb{R}^{+}.$ Each model has a different marginal mean, see Table 1. The
descriptions of the covariates can be found in R software.

Table 1: Meuse data: Covariates Model | Covariates
---|---
| $Intercept,dist,ffreq,soil$ | $dist\times ffreq$ | $dist\times Soil$
1 | ✓ |  |
2 | ✓ | ✓ |
3 | ✓ |  | ✓
4 | ✓ | ✓ | ✓
Figure 4: For each model, each symbol refers to a prediction error, estimated
by a different method.

As can be seen in Figure 4, $tAI$ estimates $-l(\boldsymbol{y^{*}})/n^{*}$
most accurately. The other prediction error estimators consistently under
estimate $-l(\boldsymbol{y^{*}})/n^{*}.$

Figure 5 is based on Meuse and on Meuse.grid data sets where the whole Meuse
data set is used as training data and the Meuse.grid data set is used as the
prediction set, $\\{X^{*},Z^{*}\\}.$ Since the Lead consternation is not given
in the Meuse.grid data set, then $-l(\boldsymbol{y^{*}})/n^{*}$ is unknown.
Therefore $tAI,\,cAI$ and $mAI$ are compared without having a ground truth.

Figure 5: For each model, each symbol refers to a prediction error, estimated
by a different method.

As can be seen from Figure 5, the differences between the $tAI,\,cAI$ and
$mAI$ are sustained and the results are consistent with the previous figures,
i.e., $cAI,\,mAI$ give lower error estimates, which likely underestimate the
prediction error.

####  4.2.2 Growth data

##### Data description

The Growth data was introduced by Potthoff and Roy (1964) and contains four
skull length measurements for 27 children at ages $8,\,10,\,12$ and $14$
(total of $27\times 4$ measurements) along with the child’s age and gender.

##### Results

Figure 6 presents a scenario where the training sample is defined as the skull
length measurements at ages $8,\,10,\,12$ and the prediction set is defined as
the skull length measurements at age $14.$ Three linear mixed models are
fitted, all have the same variance structure, containing random intercept per
child and random slope for the child’s age, however each model has a different
set of fixed effects (see Table 2).

Table 2: Growth data: Covariates Model | Covariates
---|---
| $Intercept$ | $Age$ | $Gender$ | $Age\times Gender$
1 | ✓ | ✓ |  |
2 | ✓ | ✓ | ✓ |
3 | ✓ | ✓ | ✓ | ✓
Figure 6: For each model, each symbol refers to a prediction error, estimated
by a different method.

As can be seen in Figure 6, in general perspective, $tAI$ estimates
$-l(\boldsymbol{y^{*}})/n^{*}$ most accurately. The other prediction error
estimators under estimate $-l(\boldsymbol{y^{*}})/n^{*}.$

Figure 7 presents three similar analyses as is presented in Figure 6, however
where the other time-points measurements are designated as holdout.

Figure 7: For each model, each symbol refers to a prediction error, estimated
by a different method.

When $age=8,$ the results are similar to the results in Figure 6, however,
when $age=10$ and $age=12,$ $tAI$ and $cAI$ have similar performance. This is
not surprising since in these cases $\\{X^{*},Z^{*},R^{*}\\}$ is similar to
$\\{X,Z,R\\}.$

##  5 Optimism for Prediction at Interpolation and Extrapolation Points

The formulation of $tAI$ and the derivation of $C_{tAI}$ are based on the
normality assumptions of $\boldsymbol{y}$ and $\boldsymbol{y^{*}}$ which is
commonly assumed when LMM and GPR are implemented. However, the approach that
is used for developing $tAI$ can be used for creating other prediction error
estimators that are not based on the normality assumption of
$\boldsymbol{y^{*}}$ and $\boldsymbol{y}.$ For example, in the standard
formulation of the prediction error estimator that is based on expected
Optimism correction (Efron 1986),

$Loss(Opt)=\frac{1}{n}\|\boldsymbol{y}-H\boldsymbol{y}\|^{2}_{2}+w,$

where

$w=\mathbb{E}_{\boldsymbol{y}}\left(\frac{1}{n}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}\|\boldsymbol{y^{*}}-H\boldsymbol{y}\|^{2}_{2}-\frac{1}{n}\|\boldsymbol{y}-H\boldsymbol{y}\|^{2}_{2}\right)$

and it is assumed that $\boldsymbol{y^{*}}$ and $\boldsymbol{y}$ are drawn
from the same distribution and have the same predictor, $H\boldsymbol{y}.$
However, as was already discussed in the previous sections, these conditions
are not satisfied in many use cases. The following prediction error
generalizes $Loss(Opt),$

$Loss(Opt_{t})=\frac{1}{n}\|\boldsymbol{y}-H\boldsymbol{y}\|^{2}_{2}+w_{t},$

where

$w_{t}=\mathbb{E}_{\boldsymbol{y}}\left(\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}\|\boldsymbol{y^{*}}-H^{*}\boldsymbol{y}\|^{2}_{2}-\frac{1}{n}\|\boldsymbol{y}-H\boldsymbol{y}\|^{2}_{2}\right).$

Similarly to $tAIC$ definition, given a set of models $\mathcal{H},$
$Loss(Opt_{t})$ can be used for model selection as follows

$\displaystyle
h_{best}=\underset{h\in\mathcal{H}}{\mathrm{argmin}}\,Loss_{h}(Opt_{t}),$ (5)

where $Loss_{h}(Opt_{t})$ is $Loss(Opt_{t})$ for model $h.$

Lemma 4 introduces a general expression of $w_{t}$ for predictors that are
linear in $\boldsymbol{y}.$

######  Lemma 4.

Let $\boldsymbol{y}\in\mathbb{R}^{n}$ be a random variable with mean
$\boldsymbol{\mu}$ and variance $V.$ Similarly, let
$\boldsymbol{y^{*}}\in\mathbb{R}^{n^{*}}$ be a random variable with mean
$\boldsymbol{\mu^{*}}$ and variance $V^{*}.$ In addition, let
$H\boldsymbol{y}\in R^{n}$ and $H^{*}\boldsymbol{y}\in R^{n^{*}}$ be the
predictors of $\boldsymbol{y}$ and $\boldsymbol{y^{*}}$ respectively when
$H,H^{*}$ don’t contain $\boldsymbol{y}$ and $\boldsymbol{y^{*}}$. Then

$\displaystyle w_{t}=$
$\displaystyle\frac{2}{n}tr\left(HV\right)-\frac{2}{n^{*}}tr\left(H^{*}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})\right)$
$\displaystyle+\frac{1}{n^{*}}tr\left(V^{*}\right)-\frac{1}{n}tr\left(V\right)+\frac{1}{n^{*}}tr\left(H^{*}VH^{*^{\prime}}\right)-\frac{1}{n}tr\left(HVH^{{}^{\prime}}\right)$
$\displaystyle+\frac{1}{n}tr\left(2H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}-\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}-H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{{}^{\prime}}\right)-\frac{1}{n^{*}}tr\left(2H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{*^{\prime}}}-\boldsymbol{\mu^{*}}\boldsymbol{\mu^{*^{\prime}}}-H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{*^{\prime}}\right).$

######  Corollary 5.

Given the definitions in lemma 4, when $H\boldsymbol{\mu}=\boldsymbol{\mu}$
and $H^{*}\boldsymbol{\mu}=\boldsymbol{\mu^{*}}$

$\displaystyle w_{t}=$
$\displaystyle\frac{2}{n}tr\left(HV\right)-\frac{2}{n^{*}}tr\left(H^{*}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})\right)$
$\displaystyle+\frac{1}{n^{*}}tr\left(V^{*}\right)-\frac{1}{n}tr\left(V\right)+\frac{1}{n^{*}}tr\left(H^{*}VH^{*^{\prime}}\right)-\frac{1}{n}tr\left(HVH^{{}^{\prime}}\right).$

In case $H^{*}=H,$ $V^{*}=V$ and
$V-\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})=\sigma^{2}I_{n},$ then

$w_{t}=w=\frac{2}{n}tr\left(H\right),$

which is the same result as was introduced by Hodges and Sargent (2001) for
Linear Hierarchical models.

$Loss(Opt_{t})$ is based on the squared error loss function which reflects
Euclidean distance. Other prediction error estimators which are based on
different distances, such as on Mahalanobis distance (Mahalanobis 1936), might
be suggested as well. Corollary 6 presents a penalty correction for a
prediction error estimator which is based on Mahalanobis distance.

######  Corollary 6.

Given the definitions in lemma 4, when $H\boldsymbol{\mu}=\boldsymbol{\mu}$
and $H^{*}\boldsymbol{\mu}=\boldsymbol{\mu^{*}}$

$\displaystyle\mathbb{E}_{\boldsymbol{y}}\left(\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}\|\boldsymbol{y^{*}}-H^{*}\boldsymbol{y}\|^{2}_{M}-\frac{1}{n}\|\boldsymbol{y}-H\boldsymbol{y}\|^{2}_{M}\right)$
(6)
$\displaystyle\qquad\qquad\qquad\qquad=\frac{2}{n}tr\left(R^{-1}HV\right)-\frac{2}{n^{*}}tr\left(R^{*^{-1}}H^{*}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})\right)$
$\displaystyle\qquad\qquad\qquad\qquad\,\,\,\,\,\,\,+\frac{1}{n^{*}}tr\left(R^{*^{-1}}V^{*}\right)-\frac{1}{n}tr\left(R^{-1}V\right)$
$\displaystyle\qquad\qquad\qquad\qquad\,\,\,\,\,\,\,+\frac{1}{n^{*}}tr\left(R^{*^{-1}}H^{*}VH^{*^{\prime}}\right)-\frac{1}{n}tr\left(R^{-1}HVH^{{}^{\prime}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad=2C_{tAI}-\log\left(\frac{|R^{*}|^{\frac{1}{n^{*}}}}{|R|^{\frac{1}{n}}}\right),$

where

$\displaystyle\|\boldsymbol{y^{*}}-H^{*}\boldsymbol{y}\|^{2}_{M}$
$\displaystyle=(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})^{{}^{\prime}}R^{*{-1}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})$
$\displaystyle\|\boldsymbol{y}-H\boldsymbol{y}\|^{2}_{M}$
$\displaystyle=(\boldsymbol{y}-H\boldsymbol{y})^{{}^{\prime}}R^{-1}(\boldsymbol{y}-H\boldsymbol{y}).$

The relation between eq. (6) and $C_{tAI}$ arises due to the relation between
Mahalanobis distance and the normal likelihood which $tAI$ is based on.

It is natural to use $Loss(Opt_{t})$ instead of $tAI$ for linear predictors
that don’t assume normality, such as the predictors that are used in nearest
neighbors, Nadaraya-Watson kernel regression and smoothing spline models.
Moreover, due to the form of the normal density function, many predictors that
seem to be based on the normality assumption can be alternatively interpreted
as a solution of a least squares problem or complex versions of least squares
problems like weighed least squares and penalized least squares problems. For
example, GLS can be interpreted as the solution of weighted least squares
problem with the weight matrix $V^{-1}$. Similarly, $\boldsymbol{\hat{f}^{*}}$
can be interpreted as the solution of the following problem,

$\underset{\boldsymbol{a}\in\mathbb{R}^{n^{*}},B\in\mathbb{R}^{n^{*}\times
n}}{\mathrm{min}}\mathbb{E}_{\boldsymbol{y^{*}},\boldsymbol{y}}\|\boldsymbol{y^{*}}-(\boldsymbol{a}+B\boldsymbol{y})\|^{2}_{2}.$

The proof is attached in in Appendix A. These alternative interpretations are
free from normality assumption and therefore $Loss(Opt_{t})$ can be suitable
for them. Since many predictors can be interpreted in different ways, then the
assignation of predictors to $tAI$ or to $Loss(Opt_{t})$ should refer to the
possibility to assume normality rather than to the predictor type.

##  6 Discussion and Conclusions

$tAI$ is an extension of the prediction error estimators that are used in
$cAIC$ and $mAIC,$ extending them to estimate prediction error at
interpolation and extrapolation points. As it is demonstrated in Section 3,
these use cases are common in various research fields, and particularly in
Geostatistics and health, when GPR and LMM are used for predicting at
interpolation and extrapolation points. Since GLS, linear regression and
smoothing splines can be expressed as LMM (Brumback et al. 1999), $tAI$ is
applicable for them as well.

The correction in $tAI$ is more complicated than the corrections in $cAIC$ and
$mAIC,$ which are $tr(H)/n$ and $p/n$ respectively. The correction in $tAI,$
is affected by the relations between $\mathrm{Var}(\boldsymbol{y})$ to
$\mathrm{Var}(\boldsymbol{y^{*}}),$ $\mathrm{Var}(\boldsymbol{\hat{f}})$ to
$\mathrm{Var}(\boldsymbol{\hat{f}^{*}})$ and between
$\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})$ to
$\mathrm{Var}(\boldsymbol{y}).$ When interpreting the correction as a measure
of over-fitting, the differences between the corrections gives a new
perspective about how the over-fitting is composed as a function of the
variance structure of the problem.

In many cases the variances parameters are unknown in advance and therefore
are estimated by various procedures prior the model fitting, e.g. REML in LMM
(Verbeke 1997). Estimating the variance parameters implies an extra variation
for $tAI,$ especially when the sample size is small. Estimating the in-sample
prediction error under the LMM setup when the variance parameters are unknown
was addressed by (Liang et al. 2008). Extending this to a transductive setup
is a challenge for a future work.

The numerical analyses emphasize the practical importance in using $tAI$ in
scenarios where $\\{X^{*},Z^{*},R^{*}\\}\neq\\{X,Z,R\\}$ are different. It is
noticeable especially when predicting at extrapolation points, since in this
case the differences between $\mathrm{Var}(\boldsymbol{y})$ to
$\mathrm{Var}(\boldsymbol{y^{*}})$ and between
$\mathrm{Var}(\boldsymbol{\hat{f}})$ to
$\mathrm{Var}(\boldsymbol{\hat{f}^{*}})$ can be large.

$Loss(Opt_{t})$ is another prediction error estimator for cases involving
predicting at interpolation and extrapolation points. Unlike $tAI,$
$Loss(Opt_{t})$ doesn’t assume that the observations are normally distributed
and therefore it is also applicable in various non-parametric applications.
Since many predictors that are apparently based on normal linear model can be
alternatively interpreted as solutions for the generalized least squares
problems, the assignation of predictors to $tAI$ or to $Loss(Opt_{t})$ should
refer to the possibility to assume normality rather than to the predictor
formula.

##  Appendix A Proofs

###### Proof of Theorem 1.

$\displaystyle C_{tAI}=$
$\displaystyle\mathbb{E}_{\boldsymbol{y}}\left[-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})-\left\\{-\frac{1}{n}l(\boldsymbol{y})\right\\}\right]$
$\displaystyle=$
$\displaystyle-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}l(\boldsymbol{y^{*}})+\frac{1}{n}\mathbb{E}_{\boldsymbol{y}}l(\boldsymbol{y})$
$\displaystyle=$
$\displaystyle\frac{1}{2n^{*}}\left\\{\log|R^{*}|+n^{*}\log(2\pi)+\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})^{{}^{\prime}}R^{*^{-1}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})\right\\}$
$\displaystyle-\frac{1}{2n}\left\\{\log|R|+n\log(2\pi)+\mathbb{E}_{\boldsymbol{y}}(\boldsymbol{y}-H\boldsymbol{y})^{{}^{\prime}}R^{-1}(\boldsymbol{y}-H\boldsymbol{y})\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\log\left(\frac{|R^{*}|^{\frac{1}{n^{*}}}}{|R|^{\frac{1}{n}}}\right)$
$\displaystyle+\frac{1}{2}\left\\{\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})^{{}^{\prime}}R^{*^{-1}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})-\frac{1}{n}\mathbb{E}_{\boldsymbol{y}}(\boldsymbol{y}-H\boldsymbol{y})^{{}^{\prime}}R^{-1}(\boldsymbol{y}-H\boldsymbol{y})\right\\}$

Since

$\displaystyle\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})^{{}^{\prime}}R^{*^{-1}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})$
$\displaystyle\qquad\qquad\qquad=tr\left\\{R^{*^{-1}}\mathbb{E}_{\boldsymbol{y^{*}}}(\boldsymbol{y^{*}}\boldsymbol{y^{*^{\prime}}})\right\\}+tr\left\\{H^{*^{\prime}}R^{*^{-1}}H^{*}\mathbb{E}_{\boldsymbol{y}}(\boldsymbol{y}\boldsymbol{y^{{}^{\prime}}})\right\\}$
$\displaystyle\qquad\qquad\qquad\,\,\,\,\,\,\,-2tr\left\\{R^{*^{-1}}H^{*}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}(\boldsymbol{y}\boldsymbol{y^{*^{\prime}}})\right\\}$
$\displaystyle\qquad\qquad\qquad=tr\left\\{R^{*^{-1}}(V^{*}+\boldsymbol{\mu^{*}}\boldsymbol{\mu^{*^{\prime}}})\right\\}+tr\left\\{R^{*^{-1}}H^{*}(V+\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}})H^{*^{\prime}}\right\\}$
$\displaystyle\qquad\qquad\qquad\,\,\,\,\,\,\,-2tr\left\\{R^{*^{-1}}H^{*}(\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})+\boldsymbol{\mu}\boldsymbol{\mu^{*^{\prime}}})\right\\}$

and

$\displaystyle\mathbb{E}_{\boldsymbol{y}}(\boldsymbol{y}-H\boldsymbol{y})^{{}^{\prime}}R^{-1}(\boldsymbol{y}-H\boldsymbol{y})=$
$\displaystyle
tr\left\\{R^{-1}\mathbb{E}_{\boldsymbol{y}}(\boldsymbol{y}\boldsymbol{y^{{}^{\prime}}})\right\\}+tr\left\\{(H^{{}^{\prime}}R^{-1}H\mathbb{E}_{\boldsymbol{y}}(\boldsymbol{y}\boldsymbol{y^{{}^{\prime}}})\right\\}$
$\displaystyle-2tr\left\\{R^{-1}H\mathbb{E}_{\boldsymbol{y}}(\boldsymbol{y}\boldsymbol{y^{{}^{\prime}}})\right\\}$
$\displaystyle=$ $\displaystyle
tr\left\\{R^{-1}(V+\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}})\right\\}+tr\left\\{R^{-1}H(V+\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}})H^{{}^{\prime}}\right\\}$
$\displaystyle-2tr\left\\{R^{-1}H(V+\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}})\right\\}$

then

$\displaystyle\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})^{{}^{\prime}}R^{*^{-1}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})-\frac{1}{n}\mathbb{E}_{\boldsymbol{y}}(\boldsymbol{y}-H\boldsymbol{y})^{{}^{\prime}}R^{-1}(\boldsymbol{y}-H\boldsymbol{y})$
$\displaystyle=\frac{2}{n}tr\left(R^{-1}HV\right)-\frac{2}{n^{*}}tr\left(R^{*^{-1}}H^{*}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})\right)$
$\displaystyle\,\,\,\,\,\,\,+\frac{1}{n^{*}}tr\left(R^{*^{-1}}V^{*}\right)-\frac{1}{n}tr\left(R^{-1}V\right)+\frac{1}{n^{*}}tr\left(R^{*^{-1}}H^{*}VH^{*^{\prime}}\right)-\frac{1}{n}tr\left(R^{-1}HVH^{{}^{\prime}}\right)$
$\displaystyle\,\,\,\,\,\,\,+\frac{1}{n}tr\left(R^{-1}(2H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}-\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}-H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{{}^{\prime}})\right)$
$\displaystyle\,\,\,\,\,\,\,-\frac{1}{n^{*}}tr\left(R^{*^{-1}}(2H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{*^{\prime}}}-\boldsymbol{\mu^{*}}\boldsymbol{\mu^{*^{\prime}}}-H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{*^{\prime}})\right)$

Since $H\boldsymbol{\mu}=\boldsymbol{\mu}$

$\displaystyle H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}$
$\displaystyle=\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}$ $\displaystyle
H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{{}^{\prime}}$
$\displaystyle=\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}$

and similarly since $H^{*}\boldsymbol{\mu}=\boldsymbol{\mu^{*}}$

$\displaystyle H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{*^{\prime}}}$
$\displaystyle=\boldsymbol{\mu^{*}}\boldsymbol{\mu^{*^{\prime}}}$
$\displaystyle
H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{*^{\prime}}$
$\displaystyle=\boldsymbol{\mu^{*}}\boldsymbol{\mu^{*^{\prime}}}$

then

$\displaystyle\frac{1}{n}tr\left(R^{-1}(2H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}-\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}-H\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{{}^{\prime}})\right)$
$\displaystyle=0$
$\displaystyle\frac{1}{n^{*}}tr\left(R^{*^{-1}}(2H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{*^{\prime}}}-\boldsymbol{\mu^{*}}\boldsymbol{\mu^{*^{\prime}}}-H^{*}\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}H^{*^{\prime}})\right)$
$\displaystyle=0$

and therefore

$\displaystyle\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})^{{}^{\prime}}R^{*^{-1}}(\boldsymbol{y^{*}}-H^{*}\boldsymbol{y})-\frac{1}{n}\mathbb{E}_{\boldsymbol{y}}(\boldsymbol{y}-H\boldsymbol{y})^{{}^{\prime}}R^{-1}(\boldsymbol{y}-H\boldsymbol{y}))$
$\displaystyle=\frac{2}{n}tr\left(R^{-1}HV\right)-\frac{2}{n^{*}}tr\left(R^{*^{-1}}H^{*}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})\right)$
$\displaystyle\,\,\,\,\,\,\,+\frac{1}{n^{*}}tr\left(R^{*^{-1}}V^{*}\right)-\frac{1}{n}tr\left(R^{-1}V\right)+\frac{1}{n^{*}}tr\left(R^{*^{-1}}H^{*}VH^{*^{\prime}}\right)-\frac{1}{n}tr\left(R^{-1}HVH^{{}^{\prime}}\right)$

which gives

$\displaystyle C_{tAI}=$
$\displaystyle\frac{1}{n}tr\left(R^{-1}HV\right)-\frac{1}{n^{*}}tr\left(R^{*^{-1}}H^{*}\mathrm{Cov}(\boldsymbol{y},\boldsymbol{y^{*}})\right)$
$\displaystyle+\frac{1}{2}\left\\{\log\left(\frac{|R^{*}|^{\frac{1}{n^{*}}}}{|R|^{\frac{1}{n}}}\right)+\frac{1}{n^{*}}tr\left(R^{*^{-1}}V^{*}\right)-\frac{1}{n}tr\left(R^{-1}V\right)\right\\}$
$\displaystyle+\frac{1}{2}\left\\{\frac{1}{n^{*}}tr\left(R^{*^{-1}}H^{*}VH^{*^{\prime}}\right)-\frac{1}{n}tr\left(R^{-1}HVH^{{}^{\prime}}\right)\right\\}$

∎

###### Proof of Lemma 4.

Under some regularity conditions

$\displaystyle\frac{\partial\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}\lVert\boldsymbol{y^{*}}-(\boldsymbol{a}+B\boldsymbol{y})\rVert_{2}^{2}}{\partial\boldsymbol{a}}$
$\displaystyle=\frac{\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}\partial\lVert\boldsymbol{y^{*}}-(\boldsymbol{a}+B\boldsymbol{y})\rVert_{2}^{2}}{\partial\boldsymbol{a}}$
$\displaystyle\frac{\partial\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}\lVert\boldsymbol{y^{*}}-(\boldsymbol{a}+B\boldsymbol{y})\rVert_{2}^{2}}{\partial
B}$
$\displaystyle=\frac{\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}\partial\lVert\boldsymbol{y^{*}}-(\boldsymbol{a}+B\boldsymbol{y})\rVert_{2}^{2}}{\partial
B}.$

Since

$\displaystyle\frac{\partial\rVert\boldsymbol{y^{*}}-(\boldsymbol{a}+B\boldsymbol{y})\lVert_{2}^{2}}{\partial\boldsymbol{a}}$
$\displaystyle=\frac{\partial\left(-\boldsymbol{y^{*^{\prime}}}\boldsymbol{a}-\boldsymbol{a^{{}^{\prime}}}\boldsymbol{y^{*}}+\boldsymbol{a^{{}^{\prime}}}\boldsymbol{a}+\boldsymbol{a^{{}^{\prime}}}B\boldsymbol{y}+\boldsymbol{y^{{}^{\prime}}}B^{{}^{\prime}}\boldsymbol{a}\right)}{\partial\boldsymbol{a}}$
$\displaystyle=-2\boldsymbol{y^{*}}+2\boldsymbol{a}+2B\boldsymbol{y},$

then

$\frac{\partial\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}\rVert\boldsymbol{y^{*}}-(\boldsymbol{a}+B\boldsymbol{y})\lVert_{2}^{2}}{\partial\boldsymbol{a}}=-2\boldsymbol{\mu^{*}}+2\boldsymbol{a}+2B\mu.$

Similarly

$\displaystyle\frac{\partial\rVert\boldsymbol{y^{*}}-(\boldsymbol{a}+B\boldsymbol{y})\lVert_{2}^{2}}{\partial
B}$
$\displaystyle=\frac{\partial\left(-\boldsymbol{y^{*^{\prime}}}B\boldsymbol{y}+\boldsymbol{a^{{}^{\prime}}}B\boldsymbol{y}-\boldsymbol{y^{{}^{\prime}}}B^{{}^{\prime}}\boldsymbol{y^{*}}+\boldsymbol{y^{{}^{\prime}}}B^{{}^{\prime}}\boldsymbol{a}+\boldsymbol{y^{{}^{\prime}}}B^{{}^{\prime}}B\boldsymbol{y}\right)}{\partial
B}$
$\displaystyle=-2\boldsymbol{y^{*}}\boldsymbol{y^{{}^{\prime}}}+2\boldsymbol{a}\boldsymbol{y^{{}^{\prime}}}+2B\boldsymbol{y}\boldsymbol{y^{{}^{\prime}}}$

and therefore

$\frac{\partial\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}}\rVert\boldsymbol{y^{*}}-(\boldsymbol{a}+B\boldsymbol{y})\lVert_{2}^{2}}{\partial
B}=-2\left(\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})+\boldsymbol{\mu^{*^{\prime}}}\boldsymbol{\mu}\right)+2\boldsymbol{a}\boldsymbol{\mu^{{}^{\prime}}}+2B\left(V+\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}\right).$

Since the optimized function is convex, the solution of the following
equations achieves the global minimum where

$\displaystyle 0=$
$\displaystyle-\boldsymbol{\mu^{*}}+\boldsymbol{a}+B\boldsymbol{\mu}$
$\displaystyle 0=$
$\displaystyle-\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})-\boldsymbol{\mu^{*}}\boldsymbol{\mu^{{}^{\prime}}}+\boldsymbol{a}\boldsymbol{\mu^{{}^{\prime}}}+B\left(V+\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}\right).$

The solution for $B$ is:

$\displaystyle B\left(V+\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}}\right)$
$\displaystyle=\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})+\boldsymbol{\mu^{*}}\boldsymbol{\mu^{{}^{\prime}}}-\boldsymbol{a}\boldsymbol{\mu^{{}^{\prime}}}$
$\displaystyle=\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})+\boldsymbol{\mu^{*}}\boldsymbol{\mu^{{}^{\prime}}}+\left(\boldsymbol{\mu^{*}}-B\mu\right)\boldsymbol{\mu^{{}^{\prime}}}$
$\displaystyle=\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})+B\boldsymbol{\mu}\boldsymbol{\mu^{{}^{\prime}}},$

which gives

$B=\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})V^{-1}.$

The solution for $\boldsymbol{a}$ is

$\displaystyle\boldsymbol{a}$
$\displaystyle=\boldsymbol{\mu^{*}}-B\boldsymbol{\mu}$
$\displaystyle=\boldsymbol{\mu^{*}}-\mathrm{Cov}(\boldsymbol{y^{*}},\boldsymbol{y})V^{-1}\boldsymbol{\mu}.$

Therefore the optimal linear equation is the same as
$\boldsymbol{\hat{f}^{*}}.$ ∎

##  Appendix B Scenarios in mixed model where $R\neq\sigma^{2}_{\epsilon}I$

######  Example 7.

Consider the following model

$\boldsymbol{y}=X\boldsymbol{\beta}+Z_{1}\boldsymbol{b_{1}}+Z_{2}\boldsymbol{b_{2}}+\boldsymbol{\epsilon},$

where $\boldsymbol{\epsilon}\sim
N(0,\sigma^{2}_{\epsilon}I),\,\boldsymbol{b_{1}}\sim
N(0,G_{1}),\,\boldsymbol{b_{2}}\sim N(0,G_{2}),$

$\boldsymbol{y^{*}}=X\boldsymbol{\beta}+Z_{1}\boldsymbol{b_{1}}+Z_{2}^{*}\boldsymbol{b^{*}_{2}}+\boldsymbol{\epsilon^{*}}$

where $\boldsymbol{\epsilon^{*}}\sim
N(0,\sigma^{2}_{\epsilon}I),\,\boldsymbol{b^{*}_{2}}\sim N(0,G_{2}),$

$\boldsymbol{\epsilon}\perp\boldsymbol{\epsilon^{*}}\perp\boldsymbol{b_{1}}\perp\boldsymbol{b_{2}}\perp\boldsymbol{b^{*}_{2}},$

and
$Z_{1}\in\mathbb{R}^{nXq_{1}},\,Z_{2}\in\mathbb{R}^{nXq_{2}},\,Z_{2}^{*}\in\mathbb{R}^{nXq_{2}}.$

Since $\boldsymbol{b_{1}}$ is common for $\boldsymbol{y}$ and
$\boldsymbol{y^{*}},$ its estimate can be utilized for achieving a better
accuracy in predicting $\boldsymbol{y^{*}}.$ Since $\boldsymbol{y^{*}}$
doesn’t contain $\boldsymbol{b_{2}}$ and $\boldsymbol{\epsilon}$, estimating
them doesn’t contribute achieving a better accuracy in predicting
$\boldsymbol{y^{*}}.$ Therefore, in terms of predicting $\boldsymbol{y^{*}}$
using $\boldsymbol{\hat{f}^{*}}$, the following model definition has the same
predicting formula as the previous one,

$\boldsymbol{y}=X\boldsymbol{\beta}+Z_{1}\boldsymbol{b_{1}}+\boldsymbol{\epsilon},$

where $\boldsymbol{\epsilon}\sim
N(0,\sigma^{2}_{\epsilon}I+Z_{2}G_{2}Z_{2}^{{}^{\prime}}),\,\boldsymbol{b_{1}}\sim
N(0,G_{1}),$

$\boldsymbol{y^{*}}=X\boldsymbol{\beta}+Z_{1}\boldsymbol{b_{1}}+\boldsymbol{\epsilon^{*}},$

where $\boldsymbol{\epsilon^{*}}\sim
N(0,\sigma^{2}_{\epsilon}I+Z_{2}^{*}G_{2}Z_{2}^{*^{\prime}}),$

$\boldsymbol{\epsilon}\perp\boldsymbol{\epsilon^{*}}\perp\boldsymbol{b_{1}}.$

Since the second formulation is simpler it can be preferred when the goal is
predicting $\boldsymbol{y^{*}}.$

######  Example 8.

Consider the standard LMM setup when $\boldsymbol{y}\in\mathbb{R}^{n}$ is
drawn from $K$ clusters, however each observation, $y_{i},$ is an average of
$w_{i}$ i.i.d observations, $y_{i}=(\sum_{l=1}^{w_{i}}\phi_{i,l})/w_{i},$
where $\phi_{i,l}\sim N(0,\sigma^{2}_{\epsilon}).$ Assume $\phi_{i,l}$ are
unknown, however $w_{i}$ is known. The variance of the residual, in this case
is

$\left(\begin{array}[]{cccc}\frac{\sigma^{2}_{\epsilon}}{w_{1}}&0&\cdots&0\\\
0&\frac{\sigma^{2}_{\epsilon}}{w_{2}}&&\vdots\\\ \vdots&&\ddots&0\\\
0&\cdots&0&\frac{\sigma^{2}_{\epsilon}}{w_{n}}\end{array}\right).$

Another common use case is when due to poor available data, technical
restrictions or other reasons, part of the correlation of $\boldsymbol{y}$ is
not explained by the random effects. In those cases, this part will be
expressed by the residual, $\boldsymbol{\epsilon},$ and therefore
$\mathrm{Var}(\boldsymbol{\epsilon})$ will be a non-diagonal matrix.

##  Appendix C Additional Numerical Results

Figures 8(a) and 8(b) present the distributions of $tAI,\,cAI,\,mAI$ and
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}$
for models $1$ and $2$. For more details, see Section 4.

(a) Density: Model 1

(b) Density: Model 2

Figure 8: Densities of $tAI,\,cAI,\,mAI$ and
$-\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l(\boldsymbol{y^{*}})/n^{*}$
as a function of the sample size and $\sigma^{2}$.

Figure 9 presents the error

$\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}-\frac{1}{n^{*}}l_{h_{best}}(\boldsymbol{y^{*}})$

for each one of the model selection criteria, $tAIC,\,cAIC\,mAIC$ and the
oracle criterion

$\displaystyle h_{best}$
$\displaystyle=\underset{h\in\\{1,2,3\\}}{\mathrm{argmin}}-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l_{h}(\boldsymbol{y^{*}})$

in the nine setups. For more details see Section 4.

Figure 9: For each setup, each symbol refers to the prediction error
$\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}-l_{h_{best}}(\boldsymbol{y^{*}})/n^{*}$
of the relevant criterion, $mAIC,\,cAIC\,tAIC$ and the oracle criterion.

Figure 10 presents the agreement rate of the criteria, $tAIC,\,cAIC$ and
$mAIC$ with the oracle criterion

$\displaystyle h_{best}$
$\displaystyle=\underset{h\in\\{1,2,3\\}}{\mathrm{argmin}}-\frac{1}{n^{*}}\mathbb{E}_{\boldsymbol{y}}\mathbb{E}_{\boldsymbol{y^{*}}|\boldsymbol{y}}l_{h}(\boldsymbol{y^{*}}).$

For more details see Section 4.

Figure 10: For each setup, each bar refers to the agreement rate of the
relevant criterion with the oracle criterion

## References

  * Akaike (1974) Akaike, H. (1974).  A new look at the statistical model identification.  IEEE transactions on automatic control 19(6), 716–723.
  * Baltagi (2008) Baltagi, B. H. (2008).  Forecasting with panel data.  Journal of Forecasting 27(2), 153–173.
  * Brown and Comrie (2002) Brown, D. P. and A. C. Comrie (2002).  Spatial modeling of winter temperature and precipitation in arizona and new mexico, usa.  Climate Research 22(2), 115–128.
  * Brumback et al. (1999) Brumback, B. A., D. Ruppert, and M. P. Wand (1999).  Comment.  Journal of the American Statistical Association 94(447), 794–797.
  * 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.
  * Harville (1976) Harville, D. (1976).  Extension of the gauss-markov theorem to include the estimation of random effects.  The Annals of Statistics, 384–395.
  * Ho et al. (2011) Ho, B.-C., N. C. Andreasen, S. Ziebell, R. Pierson, and V. Magnotta (2011).  Long-term antipsychotic treatment and brain volumes: a longitudinal study of first-episode schizophrenia.  Archives of general psychiatry 68(2), 128–137.
  * Hodges and Sargent (2001) Hodges, J. S. and D. J. Sargent (2001).  Counting degrees of freedom in hierarchical and other richly-parameterised models.  Biometrika 88(2), 367–379.
  * Hogan et al. (2004) Hogan, J. W., J. Roy, and C. Korkontzelou (2004).  Handling drop-out in longitudinal studies.  Statistics in medicine 23(9), 1455–1497.
  * Joachims (1999) Joachims, T. (1999).  Transductive inference for text classification using support vector machines.  In ICML, Volume 99, pp. 200–209.
  * Johnson et al. (2014) Johnson, R. A., D. W. Wichern, et al. (2014).  Applied multivariate statistical analysis, Volume 4.  Prentice-Hall New Jersey.
  * Kyriakidis and Journel (1999) Kyriakidis, P. C. and A. G. Journel (1999).  Geostatistical space–time models: a review.  Mathematical geology 31(6), 651–684.
  * Le et al. (2006) Le, Q. V., A. J. Smola, T. Gärtner, and Y. Altun (2006).  Transductive gaussian process regression with automatic model selection.  In ECML, pp. 306–317. Springer.
  * Li et al. (2008) Li, G., P. Zhang, J. Wang, E. W. Gregg, W. Yang, Q. Gong, H. Li, H. Li, Y. Jiang, Y. An, et al. (2008).  The long-term effect of lifestyle interventions to prevent diabetes in the china da qing diabetes prevention study: a 20-year follow-up study.  The Lancet 371(9626), 1783–1789.
  * Li and Heap (2011) Li, J. and A. D. Heap (2011).  A review of comparative studies of spatial interpolation methods in environmental sciences: Performance and impact factors.  Ecological Informatics 6(3), 228–241.
  * Li and Heap (2014) Li, J. and A. D. Heap (2014).  Spatial interpolation methods applied in the environmental sciences: A review.  Environmental Modelling & Software 53, 173–189.
  * Liang et al. (2008) Liang, H., H. Wu, and G. Zou (2008).  A note on conditional aic for linear mixed-effects models.  Biometrika 95(3), 773–778.
  * Mahalanobis (1936) Mahalanobis, P. C. (1936).  On the generalised distance in statistics.  Proceedings of the National Institute of Sciences of India, 1936, 49–55.
  * Mallinckrodt et al. (2003) Mallinckrodt, C. H., T. M. Sanger, S. Dubé, D. J. DeBrota, G. Molenberghs, R. J. Carroll, W. Z. Potter, and G. D. Tollefson (2003).  Assessing and interpreting treatment effects in longitudinal clinical trials with missing data.  Biological psychiatry 53(8), 754–760.
  * Manton et al. (2012) Manton, K. G., B. Singer, and R. M. Suzman (2012).  Forecasting the health of elderly populations.  Springer Science & Business Media.
  * Murray et al. (2004) Murray, D. M., S. P. Varnell, and J. L. Blitstein (2004).  Design and analysis of group-randomized trials: a review of recent methodological developments.  American journal of public health 94(3), 423–432.
  * O’neill and Temple (2012) O’neill, R. and R. Temple (2012).  The prevention and treatment of missing data in clinical trials: an fda perspective on the importance of dealing with it.  Clinical Pharmacology & Therapeutics 91(3), 550–554.
  * Pope III et al. (2002) Pope III, C. A., R. T. Burnett, M. J. Thun, E. E. Calle, D. Krewski, K. Ito, and G. D. Thurston (2002).  Lung cancer, cardiopulmonary mortality, and long-term exposure to fine particulate air pollution.  Jama 287(9), 1132–1141.
  * Potthoff and Roy (1964) Potthoff, R. F. and S. Roy (1964).  A generalized multivariate analysis of variance model useful especially for growth curve problems.  Biometrika 51(3-4), 313–326.
  * Rikken and Van Rijn (1993) Rikken, M. and R. Van Rijn (1993).  Soil Pollution with Heavy Metals: In Inquiry Into Spatial Variation, Cost of Mapping and the Risk Evaluation of Copper, Cadmium, Lead and Zinc in the Floodplains of the Meuse West of Stein, The Netherlands: Field Study Report.  University of Utrecht.
  * Salimi-Khorshidi et al. (2011) Salimi-Khorshidi, G., T. E. Nichols, S. M. Smith, and M. W. Woolrich (2011).  Using gaussian-process regression for meta-analytic neuroimaging inference based on sparse observations.  IEEE transactions on medical imaging 30(7), 1401–1416.
  * Stahl et al. (2006) Stahl, K., R. Moore, J. Floyer, M. Asplin, and I. McKendry (2006).  Comparison of approaches for spatial interpolation of daily air temperature in a large region with complex topography and highly variable station density.  Agricultural and Forest Meteorology 139(3), 224–236.
  * Stewart et al. (2009) Stewart, S. T., D. M. Cutler, and A. B. Rosen (2009).  Forecasting the effects of obesity and smoking on us life expectancy.  New England Journal of Medicine 361(23), 2252–2260.
  * Tsanas et al. (2010) Tsanas, A., M. A. Little, P. E. McSharry, and L. O. Ramig (2010).  Accurate telemonitoring of parkinson’s disease progression by noninvasive speech tests.  IEEE transactions on Biomedical Engineering 57(4), 884–893.
  * Vaida and Blanchard (2005) Vaida, F. and S. Blanchard (2005).  Conditional akaike information for mixed-effects models.  Biometrika 92(2), 351–370.
  * Verbeke (1997) Verbeke, G. (1997).  Linear mixed models for longitudinal data.  In Linear mixed models in practice, pp. 63–153. Springer.
  * Vicente-Serrano et al. (2003) Vicente-Serrano, S. M., M. A. Saz-Sánchez, and J. M. Cuadrat (2003).  Comparative analysis of interpolation methods in the middle ebro valley (spain): application to annual precipitation and temperature.  Climate research 24(2), 161–180.
  * Wilks (1932) Wilks, S. S. (1932).  Certain generalizations in the analysis of variance.  Biometrika, 471–494.
  * Wood et al. (2004) Wood, A. M., I. R. White, and S. G. Thompson (2004).  Are missing outcome data adequately handled? a review of published randomized controlled trials in major medical journals.  Clinical trials 1(4), 368–376.
  * Wray et al. (2013) Wray, N. R., J. Yang, B. J. Hayes, A. L. Price, M. E. Goddard, and P. M. Visscher (2013).  Pitfalls of predicting complex traits from snps.  Nature Reviews Genetics 14(7), 507–515.
  * Zhang et al. (2010) Zhang, Z., E. Ersoz, C.-Q. Lai, R. J. Todhunter, H. K. Tiwari, M. A. Gore, P. J. Bradbury, J. Yu, D. K. Arnett, J. M. Ordovas, et al. (2010).  Mixed linear model approach adapted for genome-wide association studies.  Nature genetics 42(4), 355–360.