# Variable Selection Using Nearest Neighbor Gaussian Processes

Konstantin Posch<EMAIL_ADDRESS>Maximilian Arbeiter  Martin
Pleschberger  Juergen Pilz  Department of Statistics, University Klagenfurt,
Austria  KAI-Kompetenzzentrum Automobil- und Industrieelektronik GmbH,
Villach, Austria

###### Abstract

A novel Bayesian approach to the problem of variable selection using Gaussian
process regression is proposed. The selection of the most relevant variables
for a problem at hand often results in an increased interpretability and in
many cases is an essential step in terms of model regularization. In detail,
the proposed method relies on so-called nearest neighbor Gaussian processes,
that can be considered as highly scalable approximations of classical Gaussian
processes. To perform a variable selection the mean and the covariance
function of the process are conditioned on a random set $\mathcal{A}$. This
set holds the indices of variables that contribute to the model. While the
specification of a priori beliefs regarding $\mathcal{A}$ allows to control
the number of selected variables, so-called reference priors are assigned to
the remaining model parameters. The application of the reference priors
ensures that the process covariance matrix is (numerically) robust. For the
model inference a Metropolis within Gibbs algorithm is proposed. Based on
simulated data, an approximation problem from computer experiments and two
real-world datasets, the performance of the new approach is evaluated.

###### keywords:

Variable Selection, Hierarchical Bayes , Nearest Neighbor Gaussian Process ,
Model Uncertainty , Metropolis-Hastings Algorithm

††journal: Journal of LaTeX Templates

##  1 Introduction

Due to the omnipresent and ongoing digitalization process more and more data
is collected by companies, governments and organizations of all kinds.
Consequently, the urge of analyzing dependencies between observed variables is
also growing. Regression techniques are often the methods of choice for this
task. The application of regression models is straight forward if there is a
large number $n$ of observations of a comparatively small number $d$ of
observed variables. However, often this does not hold true and things get
challenging. Models start to over-fit, since the model complexity cannot be
supported by the available sample size [Buehlmann et al., 2016]. Moreover, it
becomes difficult to identify the most influential variables. For this reason,
approaches that automatically detect the most important variables gain
increasing attention in many fields including genetics [Lee et al., 2003],
astronomy [Zheng and Zhang, 2007], and economics [Foster and Stine, 2004].

In case that the dependencies of interest are purely linear simple linear
regression models are sufficient for the analysis. There are many well-
established approaches to perfrom variable selection in linear regression
models [Tibshirani, 1996, Fan and Li, 2001, Zou and Hastie, 2005, Zou, 2006,
Park and Casella, 2008, Ročková and George, 2014, Alhamzawi and Taha Mohammad
Ali, 2018, Posch et al., 2020]. However, in many practical applications
linearity does not hold and more sophisticated methods are required. Recently,
Gaussian process regression models (GPRMs), originally proposed by O’Hagan
[1978], gained a lot of attention. GPRMs can be considered as an alternative
to neural networks since a large class of neural network based regression
models converges to a Gaussian Process (GP) in the limit of an infinite
network [Neal, 1996]. Alternatively, GPs can also be viewed at from the
perspective of non-parametric Bayesian regression [Rasmussen and Williams,
2005]. Here, Gaussian priors are directly assigned to the space of regression
functions. There are many empirical studies which confirm that GPRMs perform
very well compared to other non-linear models [Chen et al., 2007, Yuan et al.,
2008]. Moroever, GP regression has become a very popular approach for the
approximation of time-expensive computer models [Vollert et al., 2019].

A large class of GPRMs automatically provide some information regarding the
importance of the individual variables the model is fitted with. This property
is known as automatic relevance determination [Rasmussen and Williams, 2005].
However, in case of limited training data such approaches perform poorly,
resulting in a high generalization error. One possibility to overcome this
problem is to add a penalty term, that encourages model fits, where only some
of the variables are considered as important. Respective approaches were
proposed by Yi et al. [2011] and Wu et al. [2017]. A related Bayesian approach
was proposed by Linkletter et al. [2006], who used a spike and slab prior
(mixture prior between a uniform distribution and a point mass) for the
variable selection task. Chen and Wang [2010] proposed a Bayesian approach
that uses a random indicator vector to specify which variables contribute to
the model.

A possible weakness of GPRMs is that the estimation process sometimes becomes
numerically unstable. To address this issue Gu et al. [2018] proposed the
concept of robust estimation. In a Bayesian manner, robust estimates can be
achieved by assigning so-called reference priors to the GP parameters. These
priors are recommended by many authors [Berger et al., 2001, Paulo, 2005,
Kazianka and Pilz, 2012, Vollert et al., 2019]. Another disadvantage of GPRMs
is that they are computationally expensive. Recently, [Datta et al., 2016]
introduced so-called nearest neighbor Gaussian processes (NNGPs). They can be
considered as highly scalable approximations derived from classical GPs.

In this work we propose a variable selection approach combining NNGPs and
reference priors. To the best of our knowledge such a combination has not been
proposed before. The variable selection is based on a random set $\mathcal{A}$
used tho specify variables that contribute to the model. For the model
inference we propose a Metropolis within Gibbs algorithm to sample from the
posterior distribution. In particular, the proposed MCMC algorithm makes use
of Hamiltonian Monte Carlo (HMC) [Duane et al., 1987]. HMC is a family of MCMC
methods based on the concept of dynamic systems in physics. Virtual dynamical
systems are created by augmenting hyperparameters with momentum variables and
in the end a sampling from the distribution of the combined system is
performed. In a Bayesian context HMC uses the derivatives of the posterior
density to generate efficient proposals. In many applications HMC converges
significantly faster than conventional approaches [Duane et al., 1987, Neal,
1997, Rasmussen and Williams, 2005].

##  2 Methodology

###  2.1 Nearest Neighbor Gaussian Process

Formally, a stochastic process $\left\\{y(\mathbf{x});\mathbf{x}\in
X\subseteq\mathbb{R}^{d}\right\\}$ is called a Gaussian Process iff for every
finite set of indices
$\left\\{\mathbf{x}_{1},...,\mathbf{x}_{k}\right\\}\subseteq X$ the random
vector
$y_{\mathbf{x}_{1},...,\mathbf{x}_{k}}=(y(\mathbf{x}_{1}),...,y(\mathbf{x}_{k}))^{T}$
is multivariate normally distributed. A GP is fully determined by its mean
function
$\mu:X\rightarrow\mathbb{R},~{}\mathbf{x}\mapsto\mathbb{E}(y(\mathbf{x}))$ and
its covariance function $C:X\times
X\rightarrow\mathbb{R},~{}(\mathbf{x_{1}},\mathbf{x_{2}})\mapsto\operatorname{Cov}(y(\mathbf{x}_{1}),y(\mathbf{x}_{2}))$.
In the context of nonparametric Bayesian regression GPs are often used as
prior distribution for the unknown regression function. Performing inference
in such a model reduces to the estimation of some unknown parameters used to
define the mean function $\mu(\cdot)$ and the covariance function
$C(\cdot,\cdot)$. In the regression context the indices $\mathbf{x}_{i}$ are
also called covariates or predictors. Let
$D=\\{(\mathbf{x}_{i},y(\mathbf{x}_{i}))|i=1,...,n\\}$ be a finite realization
of a GP, i.e., some training data. Further, let
$S=\\{\mathbf{x}_{1},...,\mathbf{x}_{n}\\}$ denote the set of covariate
constellations for which observations are available. Then,
$\mathbf{C}_{S}=[C(\mathbf{x}_{i},\mathbf{x}_{j})]_{\mathbf{x}_{i},\mathbf{x}_{j}\in
S}$ denotes the corresponding covariance matrix. Due to the form of the
Gaussian density, the model inference is based on the quantities
$\mathbf{C}_{S}^{-1}$ and $|\mathbf{C}_{S}|$. Since the computation of inverse
matrices and also the calculation of determinants is computationally
expensive, the application of GP regression models becomes challenging in case
that $n$ is large. Recently, Datta et al. [2016] introduced so-called nearest
neighbor Gaussian processes as highly scalable approximations derived from
classical GP models. Starting point is a zero mean GP, i.e., $\mu(\cdot)\equiv
0$. Then the main idea behind NNGPs is to write the joint density of
$\mathbf{y}_{S}=(y(\mathbf{x}_{1}),...,y(\mathbf{x}_{n}))^{T}$ as a product of
conditional densities, using the multiplication theorem, and subsequently
replacing the conditioning sets by smaller ones of size at most $m<n$. This is

$\displaystyle p(\mathbf{y}_{S})$
$\displaystyle=p(y(\mathbf{x}_{1}))\prod\limits_{i=2}^{n}p(y(\mathbf{x}_{i})|y(\mathbf{x}_{i-1}),...,y(\mathbf{x}_{1}))$
(2.1) $\displaystyle\tilde{p}(\mathbf{y}_{S})$
$\displaystyle=\prod\limits_{i=1}^{n}p(y(\mathbf{x}_{i})|y_{N(\mathbf{x}_{i})})$
(2.2)

where $y_{N(\mathbf{x}_{i})}$ denotes the vector formed by stacking all the
$y(\mathbf{x})$ for which $\mathbf{x}$ lies in the $N(\mathbf{x}_{i})\subset
S\setminus\\{\mathbf{x}_{i}\\}$, $|N(\mathbf{x}_{i})|\leq m$, the so-called
neighbor set of $\mathbf{x}_{i}$. Defining the neighbor set
$N(\mathbf{x}_{i})$ to be any subset of
$\\{\mathbf{x}_{1},...,\mathbf{x}_{i-1}\\}$ ensures that Equation 2.2
specifies a valid probability density [Datta et al., 2016]. In this work we
focus on the specification of Vecchia [1988], who specified
$N(\mathbf{x}_{i})$ to be the $m$ nearest neighbors of $\mathbf{x}_{i}$ among
$\\{\mathbf{x}_{1},...,\mathbf{x}_{i-1}\\}$ with respect to the Euclidean
distance. Note that the density $\tilde{p}(\mathbf{y}_{S})$ is given by

$\displaystyle\tilde{p}(\mathbf{y}_{S})=\prod\limits_{i=1}^{n}\mathcal{N}(y(\mathbf{x}_{i})|\mathbf{B}_{\mathbf{x_{i}}}y_{N(\mathbf{x}_{i})},F_{\mathbf{x}_{i}})$
(2.3)

where
$\mathbf{B}_{\mathbf{x_{i}}}=\mathbf{C}_{\mathbf{x_{i}},N(\mathbf{x}_{i})}\mathbf{C}_{N(\mathbf{x}_{i})}^{-1}$
and
$F_{\mathbf{x}_{i}}=\mathbf{C}(\mathbf{x}_{i},\mathbf{x}_{i})-\mathbf{C}_{\mathbf{x_{i}},N(\mathbf{x}_{i})}\mathbf{C}_{N(\mathbf{x}_{i})}^{-1}\mathbf{C}_{N(\mathbf{x}_{i}),\mathbf{x_{i}}}$.
Equation 2.3 holds due to Equation 2.2 and the fact that for

$\mathbf{x}=\begin{pmatrix}\mathbf{x}_{1}\\\
\mathbf{x}_{2}\end{pmatrix}\sim\mathcal{N}\left(\begin{pmatrix}\boldsymbol{\mu}_{1}\\\
\boldsymbol{\mu}_{2}\end{pmatrix},\begin{pmatrix}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12}\\\
\boldsymbol{\Sigma}_{21}&\boldsymbol{\Sigma}_{22}\end{pmatrix}\right)$

$(\mathbf{x}_{1}|\mathbf{x}_{2}=\mathbf{a})\sim\mathcal{N}(\bar{\boldsymbol{\mu}},\bar{\boldsymbol{\Sigma}})$
with
$\bar{\boldsymbol{\mu}}=\boldsymbol{\mu}_{1}+\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}(\mathbf{a}-\boldsymbol{\mu}_{2})$
and
$\bar{\boldsymbol{\Sigma}}=\boldsymbol{\Sigma}_{11}-\boldsymbol{\Sigma}_{12}\boldsymbol{\Sigma}_{22}^{-1}\boldsymbol{\Sigma}_{21}$.

Following Datta et al. [2016] the density $\tilde{p}(\mathbf{y}_{S})$ can be
extended to a legitimate stochastic process conform to Komogorov’s consistency
criteria. Let $\mathbf{u}\in X\setminus S$. Further, let $N(\mathbf{u})$ be
the set of $m$-nearest neighbors of $\mathbf{u}$ in $S$. For any finite set
$U=\\{\mathbf{u}_{1},...,\mathbf{u}_{r}\\}$ with $S\cap U=\emptyset$ we define
the nearest neighbor density of $\mathbf{y}_{U}$ conditional on
$\mathbf{y}_{S}$ as

$\displaystyle\tilde{p}(\mathbf{y}_{U}|\mathbf{y}_{S})=\prod\limits_{i=1}^{r}p(y(\mathbf{u}_{i})|\mathbf{y}_{N(\mathbf{u}_{i})}).$
(2.4)

This conditional density is similar to 2.2 except that all the neighbor sets
are subsets of $S$. This ensures a proper conditional density. Indeed Equation
2.2 and Equation 2.4 are sufficient to describe the joint (nearest neighbor)
density of any finite set $V=\\{\mathbf{v}_{1},...,\mathbf{v}_{n}\\}\subseteq
X$

$\displaystyle\tilde{p}(\mathbf{y}_{V})=\int\tilde{p}(\mathbf{y}_{U}|\mathbf{y}_{S})\tilde{p}(\mathbf{y}_{S})~{}\prod\limits_{\\{\mathbf{s}_{i}\in
S\setminus V\\}}d\mathbf{y}(\mathbf{s}_{i})$ (2.5)

where $U=V\setminus S$. If $U=\emptyset$, then Equation 2.4 implies that
$\tilde{p}(\mathbf{y}_{U}|\mathbf{y}_{S})=1$ in Equation 2.5. If $S\setminus
V=\emptyset$ the integration in Equation 2.5 is not needed. Moreover, Datta et
al. [2016] could show that $\tilde{p}(\mathbf{y}_{V})$ is the density of a
zero mean GP with covariance function

$\displaystyle\tilde{C}(\mathbf{v}_{1},\mathbf{v}_{2})=\begin{cases}\tilde{C}_{\mathbf{x}_{i},\mathbf{x}_{j}}&\text{if
}\mathbf{v}_{1}=\mathbf{x}_{i},\mathbf{v}_{2}=\mathbf{x}_{j}\text{ are both in
}S\\\
\mathbf{B}_{\mathbf{v}_{1}}\tilde{\mathbf{C}}_{N(\mathbf{v}_{1}),\mathbf{x}_{j}}&\text{if
}\mathbf{v}_{1}\notin S,\mathbf{v}_{2}=\mathbf{x}_{j}\in S\\\
\mathbf{B}_{\mathbf{v}_{1}}\tilde{\mathbf{C}}_{N(\mathbf{v}_{1}),N(\mathbf{v}_{2})}\mathbf{B}_{\mathbf{v}_{2}}^{T}+\delta_{(\mathbf{v}_{1}=\mathbf{v}_{2})}F_{\mathbf{v}_{1}}&\text{if
}\mathbf{v}_{1},\mathbf{v}_{2}\notin S\end{cases}$ (2.6)

where

$\displaystyle\tilde{\mathbf{C}}_{S}=\mathbf{B}_{S}^{-1}\mathbf{F}_{S}(\mathbf{B}_{S}^{T})^{-1},$
(2.7)

with

$\displaystyle\mathbf{B}_{S}$ $\displaystyle=\begin{pmatrix}1\\\
&1&&\mathcal{O}\\\ &&\ddots\\\ &B_{\mathbf{x}_{i},j}&&1\\\
&&&&&1\end{pmatrix},~{}~{}~{}B_{\mathbf{x}_{i},j}=\begin{cases}1&i=j\\\
-\mathbf{B}_{\mathbf{x}_{i}}[l]&\mathbf{x}_{j}=\mathbf{x}_{i}^{l}\text{ for
some $l$},\\\ 0&\text{else}\end{cases}$ (2.8)
$\displaystyle\mathbf{B}_{\mathbf{v}_{i}}$
$\displaystyle=\mathbf{C}_{\mathbf{v}_{i},N(\mathbf{v}_{i})}\mathbf{C}_{N(\mathbf{v}_{i})}^{-1},$
(2.9) $\displaystyle\mathbf{F}_{S}$
$\displaystyle=\operatorname{diag}(F_{\mathbf{x}_{1}},...,F_{\mathbf{x}_{n}}),$
(2.10) $\displaystyle F_{\mathbf{v}_{i}}$
$\displaystyle=C_{\mathbf{v}_{i}}-\mathbf{B}_{\mathbf{v}_{i}}\mathbf{C}_{N(\mathbf{v}_{i}),\mathbf{v}_{i}},$
(2.11)

$\mathbf{x}_{j}=\mathbf{x}_{i}^{l}$ means that $\mathbf{x}_{j}$ is the $l$-th
neighbor of $\mathbf{x}_{i}$
($N(\mathbf{x}_{i})=\left\\{\mathbf{x}_{i}^{1},\mathbf{x}_{i}^{2},...,\mathbf{x}_{i}^{l},...,\mathbf{x}_{i}^{|N(\mathbf{x}_{i})|}\right\\}$)
, $\mathbf{v}_{1},\mathbf{v}_{2}\in X$, $\tilde{\mathbf{C}}_{A,B}$ denotes
submatrices of $\tilde{\mathbf{C}}$ indexed by the elements of the sets $A$
and $B$, and $\delta_{(\mathbf{v}_{1}=\mathbf{v}_{2})}$ is the Kronecker
delta.

###  2.2 Reference Prior

In terms of robust estimation of the learnable parameters of a GP the
application of so-called reference priors has become popular [Gu et al., 2018,
Vollert et al., 2019]. Paulo [2005] has computed the general form of the
reference prior for a broad class of Gaussian processes. Assume that
$y(\cdot)$ is a GP on $X\subseteq\mathbb{R}^{d}$ with mean and covariance
given by

$\displaystyle\mu(\mathbf{x})=\mathbf{x}^{T}\boldsymbol{\beta},~{}C(\mathbf{x}_{1},\mathbf{x}_{2})=\sigma^{2}K(\mathbf{x}_{1},\mathbf{x}_{2}|\boldsymbol{\xi}),$
(2.12)

where $\boldsymbol{\beta}\in\mathbb{R}^{d}$,
$\boldsymbol{\xi}\in\mathbb{R}_{+}^{r}$, and $\sigma^{2}>0$ is the process
variance. In order to define the likelihood function corresponding to the
observed data $D$, consider the following additional notations:

$\displaystyle\mathbf{X}$
$\displaystyle=(\mathbf{x}_{1}\cdots\mathbf{x}_{n})^{T},$ (2.13)
$\displaystyle\boldsymbol{\eta}$
$\displaystyle=(\sigma^{2},\boldsymbol{\beta}^{T},\boldsymbol{\xi}^{T})^{T},$
(2.14) $\displaystyle\mathbf{K}_{S}$
$\displaystyle=[K(\mathbf{x}_{i},\mathbf{x}_{j}|\boldsymbol{\xi})]_{\mathbf{x}_{i},\mathbf{x}_{j}\in
S}.$ (2.15)

Then,
$\mathbf{y}_{S}|\boldsymbol{\eta}\sim\mathcal{N}(\mathbf{X}\boldsymbol{\beta},\sigma^{2}\mathbf{K}_{S})$,
which implies that the likelihood is given by:

$\displaystyle
L(\boldsymbol{\eta}|\mathbf{y}_{S})\propto(\sigma^{2})^{-\frac{n}{2}}|\mathbf{K}_{S}|^{-\frac{1}{2}}\operatorname{exp}\left\\{-\frac{1}{2\sigma^{2}}(\mathbf{y}_{S}-\mathbf{X}\boldsymbol{\beta})^{T}\mathbf{K}_{S}^{-1}(\mathbf{y}_{S}-\mathbf{X}\boldsymbol{\beta})\right\\}$
(2.16)

In order to define the reference prior $(\sigma^{2},\boldsymbol{\xi}^{T})^{T}$
is considered as parameter of interest and $\boldsymbol{\beta}$ is considered
as nuisance parameter. Moreover, the reference prior distribution
$\pi^{R}(\boldsymbol{\eta})$ is factored as
$\pi^{R}(\boldsymbol{\eta})=\pi^{R}(\boldsymbol{\beta}|\sigma^{2},\boldsymbol{\xi})\pi^{R}(\sigma^{2},\boldsymbol{\xi})$.
The conditional prior of $\boldsymbol{\beta}$ is specified proportional to
$1$, $\pi^{R}(\boldsymbol{\beta}|\sigma^{2},\boldsymbol{\xi})\propto 1$,
because this is the Jeffreys-rule prior for the model at hand when
$\sigma^{2}$ and $\boldsymbol{\xi}$ are considered to be known. Next
$\pi^{R}(\sigma^{2},\boldsymbol{\xi})$ is calculated as the Jeffreys-rule
prior for the marginal experiment defined via the integrated likelihood with
respect to $\pi^{R}(\boldsymbol{\beta}|\sigma^{2},\boldsymbol{\xi})$. The
integrated likelihood is given by:

$\displaystyle L^{I}(\sigma^{2},\boldsymbol{\xi}|\mathbf{y}_{S})$
$\displaystyle=\int\limits_{\mathbb{R}^{d}}L(\boldsymbol{\eta}|\mathbf{y}_{S})\pi^{R}(\boldsymbol{\beta}|\sigma^{2},\boldsymbol{\xi})~{}d\boldsymbol{\beta}$
(2.17)
$\displaystyle=\int\limits_{\mathbb{R}^{d}}L(\boldsymbol{\eta}|\mathbf{y}_{S})~{}d\boldsymbol{\beta}$
(2.18)
$\displaystyle\propto(\sigma^{2})^{-\frac{n-d}{2}}|\mathbf{K}_{S}|^{-\frac{1}{2}}||\mathbf{X}^{T}\mathbf{K}_{S}^{-1}\mathbf{X}|^{-\frac{1}{2}}\operatorname{exp}\left\\{-\frac{S_{\boldsymbol{\xi}}^{2}}{2\sigma^{2}}\right\\}$
(2.19)

where $S_{\boldsymbol{\xi}}^{2}=\mathbf{y}_{S}^{T}\mathbf{Q}\mathbf{y}_{S}$,
$\mathbf{Q}=\mathbf{K}_{S}^{-1}\mathbf{P}$, and
$\mathbf{P}=\mathbf{I}-\mathbf{X}(\mathbf{X}^{T}\mathbf{K}_{S}^{-1}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{K}_{S}^{-1}$.
Finally, the reference prior $\pi^{R}(\boldsymbol{\eta})$ is of the form
$\pi^{R}(\boldsymbol{\eta})\propto\frac{\pi^{R}(\boldsymbol{\xi})}{\sigma^{2}}$
with
$\pi^{R}(\boldsymbol{\xi})\propto|\tilde{I}_{R}(\boldsymbol{\xi})|^{\frac{1}{2}}$
where

$\displaystyle\mathbf{W}_{k}$
$\displaystyle=\frac{\partial\mathbf{K}_{S}}{\partial\xi_{k}}\mathbf{Q}~{}~{}~{}k=1,...,r,$
(2.20) $\displaystyle\tilde{I}_{R}(\boldsymbol{\xi})$
$\displaystyle=\begin{pmatrix}n-d&\operatorname{tr}(\mathbf{W}_{1})&\operatorname{tr}(\mathbf{W}_{2})&\cdots&\operatorname{tr}(\mathbf{W}_{r})\\\
\operatorname{tr}(\mathbf{W}_{1})&\operatorname{tr}(\mathbf{W}_{1}^{2})&\operatorname{tr}(\mathbf{W}_{1}\mathbf{W}_{2})&\cdots&\operatorname{tr}(\mathbf{W}_{1}\mathbf{W}_{r})\\\
\vdots&\vdots&\vdots&&\vdots\\\
\operatorname{tr}(\mathbf{W}_{r})&\operatorname{tr}(\mathbf{W}_{1}\mathbf{W}_{r})&\operatorname{tr}(\mathbf{W}_{2}\mathbf{W}_{r})&\cdots&\operatorname{tr}(\mathbf{W}_{r}^{2})\end{pmatrix},$
(2.21)

see [Paulo, 2005].

###  2.3 Proposed Variable Selection Model

In this section, the proposed variable selection model is described. To this
aim at first the likelihood function is defined. Afterwards, in a Bayesian
manner, diverse prior distributions are assigned to the model parameters.

####  2.3.1 Likelihood Function

Let $y\in\mathbb{R}$ denote a target variable of interest. Further, let
$\mathbf{x}\in\mathbb{R}^{d}$ denote predictors. It is assumed that the
relationship

$\displaystyle y=\mathbf{x}^{T}\boldsymbol{\beta}+f(\mathbf{x})+\varepsilon$
(2.22)

holds, where $f$ denotes an unknown function and $\varepsilon$ additive noise.
Moreover, let $\mathcal{A}$ denote a (random) subset of the index set
$\left\\{1,...,d\right\\}$, with $|\mathcal{A}|=k>0$. This set is used to
specify the covariates that contribute to the regression model. For given
$\mathcal{A}$ we assign a GP prior with zero mean and covariance function
$\tilde{C}$ to the sum $z:=f+\varepsilon$:

$\displaystyle z\sim
GP(0,\tilde{C}(\cdot,\cdot|\sigma^{2},\gamma,\rho,\mathcal{A})),$ (2.23)

with $\sigma^{2}>0,~{}\rho>0$ and $\gamma\in(0,1)$. In particular,
$GP(0,\tilde{C}(\cdot,\cdot|\sigma^{2},\gamma,\rho,\mathcal{A}))$ is defined
as the NNGP derived from the parent
$GP(0,C(\cdot,\cdot|\sigma^{2},\gamma,\rho,\mathcal{A}))$, with

$\displaystyle
C(\mathbf{x},\mathbf{x}^{\prime}|\sigma^{2},\gamma,\rho,\mathcal{A})=\sigma^{2}\underbrace{\left[\delta_{(\mathbf{x}=\mathbf{x}^{\prime})}(1-\gamma)+\gamma
K(\mathbf{x},\mathbf{x}^{\prime}|\rho,\mathcal{A})\right]}_{=\mathcal{K}(\mathbf{x},\mathbf{x}^{\prime}|\gamma,\rho,\mathcal{A})},$
(2.24)

where

$\displaystyle\mathcal{K}(\mathbf{x},\mathbf{x}^{\prime}|\rho,\mathcal{A})=\left(1+\frac{\sqrt{5}d_{\mathcal{A}}}{\rho}+\frac{5d_{\mathcal{A}}^{2}}{3\rho^{2}}\right)\operatorname{exp}\left(-\frac{\sqrt{5}d_{\mathcal{A}}}{\rho}\right),$
(2.25)

with

$\displaystyle
d_{\mathcal{A}}=\sqrt{\sum\limits_{i\in\mathcal{A}}(x_{i}-x_{i}^{\prime})^{2}}.$
(2.26)

Therefore, the parent GP has a Matern-$\frac{5}{2}$ covariance function.
Moreover, the Euclidean distance function $d_{\mathcal{A}}$ is based on
$\mathcal{A}$ and considers only those predictors with indices in
$\mathcal{A}$. Please note that $d_{\mathcal{A}}$ is also used to determine
the neighbor sets required for the computation of the NNGP.

We want to point out that the parent GP results from assigning a zero mean GP
with covariance function
$\sigma_{f}^{2}K(\mathbf{x},\mathbf{x}^{\prime}|\rho,\mathcal{A})$ to $f$ and,
moreover, assuming the noise $\varepsilon$ to be $i.i.d.$
$\mathcal{N}(0,\sigma_{\varepsilon}^{2})$. The reparameterization
$\sigma^{2}=\sigma_{f}^{2}+\sigma_{\varepsilon}^{2},\gamma=1-\frac{\sigma_{\varepsilon}^{2}}{\sigma^{2}}$
in (2.24) leads to:

$\displaystyle
C(\mathbf{x},\mathbf{x}^{\prime}|\sigma_{f}^{2},\sigma_{\varepsilon}^{2},\rho,\mathcal{A})$
$\displaystyle=\sigma_{\varepsilon}^{2}\delta_{(\mathbf{x}=\mathbf{x}^{\prime})}+\sigma_{f}^{2}K(\mathbf{x},\mathbf{x}^{\prime}|\rho,\mathcal{A})$
(2.27)
$\displaystyle=\sigma^{2}\left[\frac{\sigma_{\varepsilon}^{2}}{\sigma^{2}}\delta_{(\mathbf{x}=\mathbf{x}^{\prime})}+\frac{\sigma^{2}-\sigma_{\varepsilon}^{2}}{\sigma^{2}}K(\mathbf{x},\mathbf{x}^{\prime}|\rho,\mathcal{A})\right]$
(2.28)
$\displaystyle=\sigma^{2}\left[(1-\gamma)\delta_{(\mathbf{x}=\mathbf{x}^{\prime})}+\gamma
K(\mathbf{x},\mathbf{x}^{\prime}|\rho,\mathcal{A})\right].$ (2.29)

Above specifications directly imply that $y$ follows a Gaussian process with
mean function $\mu(\mathbf{x})=\mathbf{x}^{T}\boldsymbol{\beta}$ and
covariance function
$\tilde{C}(\mathbf{x},\mathbf{x}^{\prime}|\sigma^{2},\gamma,\rho,\mathcal{A})$.
Let $D=\\{(\mathbf{x}_{i},y_{i})|i=1,...,n\\}$ denote some training data and,
further, let $S=\\{\mathbf{x}_{1},...,\mathbf{x}_{n}\\}$ denote the
corresponding set of covariate constellations. Consider the following
additional notations:

$\displaystyle\tilde{\mathbf{C}}_{S}$
$\displaystyle=[\tilde{C}(\mathbf{x}_{i},\mathbf{x}_{j}|\sigma^{2},\gamma,\rho,\mathcal{A})]_{\mathbf{x}_{i},\mathbf{x}_{j}\in
S}.$ (2.30)

According to Section 2.1 $\tilde{\mathbf{C}}_{S}$ computes as:

$\displaystyle\tilde{\mathbf{C}}_{S}=\mathbf{B}_{S}^{-1}\mathbf{F}_{S}(\mathbf{B}_{S}^{T})^{-1}$
(2.31)

with

$\displaystyle\mathbf{B}_{S}$ $\displaystyle=\begin{pmatrix}1\\\
&1&&\mathcal{O}\\\ &&\ddots\\\ &B_{\mathbf{x}_{i},j}&&1\\\
&&&&&1\end{pmatrix},~{}~{}~{}B_{\mathbf{x}_{i},j}=\begin{cases}1&i=j\\\
-\mathbf{B}_{\mathbf{x}_{i}}[l]&\mathbf{x}_{j}=\mathbf{x}_{i}^{l}\text{ for
some $l$},\\\ 0&\text{else}\end{cases}$ (2.32)
$\displaystyle\mathbf{B}_{\mathbf{x}_{i}}$
$\displaystyle=\mathbf{C}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}\mathbf{C}_{N(\mathbf{x}_{i})}^{-1}=\sigma^{2}\boldsymbol{\mathcal{K}}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}(\sigma^{2})^{-1}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}$
(2.33)
$\displaystyle=\boldsymbol{\mathcal{K}}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}~{}~{}~{}\text{(independent
of $\sigma^{2}$)},$ (2.34) $\displaystyle\mathbf{F}_{S}$
$\displaystyle=\operatorname{diag}(F_{\mathbf{x}_{1}},...,F_{\mathbf{x}_{n}}),$
(2.35) $\displaystyle F_{\mathbf{x}_{i}}$
$\displaystyle=C_{\mathbf{x}_{i}}-\mathbf{B}_{\mathbf{x}_{i}}\mathbf{C}_{N(\mathbf{x}_{i}),\mathbf{x}_{i}}$
(2.36)
$\displaystyle=\sigma^{2}\mathcal{K}_{\mathbf{x}_{i}}-\mathbf{B}_{\mathbf{x}_{i}}\sigma^{2}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i}),\mathbf{x}_{i}}$
(2.37)
$\displaystyle=\sigma^{2}\underbrace{\left[\mathcal{K}_{\mathbf{x}_{i}}-\mathbf{B}_{\mathbf{x}_{i}}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i}),\mathbf{x}_{i}}\right]}_{=\mathcal{F}_{\mathbf{x}_{i}}~{}~{}~{}\text{(independent
of $\sigma^{2}$)}}.$ (2.38)

Overall one obtains:

$\displaystyle\tilde{\mathbf{C}}_{S}$
$\displaystyle=\sigma^{2}\underbrace{\mathbf{B}_{S}^{-1}\operatorname{diag}(\mathcal{F}_{\mathbf{x}_{1}},...,\mathcal{F}_{\mathbf{x}_{n}})(\mathbf{B}_{S}^{T})^{-1}}_{=:\boldsymbol{\tilde{\mathcal{K}}}_{S}\text{~{}(independent
of $\sigma^{2}$, dependent of $\rho,\gamma,\mathcal{A}$)}}$ (2.39)
$\displaystyle=\sigma^{2}\boldsymbol{\tilde{\mathcal{K}}}_{S}.$ (2.40)

Thus,
$\mathbf{y}_{S}|\sigma^{2},\gamma,\rho,\mathcal{A}\sim\mathcal{N}(\mathbf{X}\boldsymbol{\beta},\sigma^{2}\boldsymbol{\tilde{\mathcal{K}}}_{S})$,
which implies that the likelihood function is given by:

$\displaystyle
L(\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A}|\mathbf{y}_{S})\propto(\sigma^{2})^{-\frac{n}{2}}|\boldsymbol{\tilde{\mathcal{K}}}_{S}|^{-\frac{1}{2}}\operatorname{exp}\left\\{-\frac{1}{2\sigma^{2}}(\mathbf{y}_{S}-\mathbf{X}\boldsymbol{\beta})^{T}\boldsymbol{\tilde{\mathcal{K}}}_{S}^{-1}(\mathbf{y}_{S}-\mathbf{X}\boldsymbol{\beta})\right\\}.$
(2.41)

####  2.3.2 Reference Prior

For given $\mathcal{A}$ we denote the components of $\boldsymbol{\beta}$ with
indices in $\mathcal{A}$ as $\beta_{\mathcal{A}}$ an the remaining ones as
$\boldsymbol{\beta}_{\mathcal{A}^{C}}$. We assume that
$\boldsymbol{\beta}_{\mathcal{A}^{C}}$ is equal to zero and, moreover, a
uniform prior is assigned to $\beta_{\mathcal{A}}$. The joint prior density
with respect to the measure $(\lambda+\delta_{0})^{d}$ is given by

$\displaystyle
p(\boldsymbol{\beta}|\mathcal{A})=p(\boldsymbol{\beta}_{\mathcal{A}}|\mathcal{A})p(\boldsymbol{\beta}_{\mathcal{A}^{C}}|\mathcal{A})\propto
I_{\\{\boldsymbol{\beta}_{\mathcal{A}^{C}}=\mathbf{0}\\}}(\boldsymbol{\beta}_{\mathcal{A}^{C}})I_{\\{\beta_{i}\neq
0|\forall i\in\mathcal{A}\\}}(\boldsymbol{\beta}_{\mathcal{A}}),$ (2.42)

where $\lambda$ denotes the Lebesgue measure and $\delta_{0}$ denotes the
Dirac measure concentrated at zero. Thus, the integrated likelihood (compare
to Section 2.2) is given by:

$\displaystyle L^{I}(\sigma^{2},\rho,\gamma,\mathcal{A}|\mathbf{y}_{S})$
$\displaystyle=\int\limits_{\mathbb{R}^{d}}L(\sigma^{2},\rho,\gamma,\mathcal{A},\boldsymbol{\beta}|\mathbf{y}_{S})p(\boldsymbol{\beta}|\mathcal{A})~{}d\lambda(\boldsymbol{\beta}_{\mathcal{A}})d(\delta_{0}(\boldsymbol{\beta}_{\mathcal{A}}^{C}))$
(2.43)
$\displaystyle\propto(\sigma^{2})^{-\frac{n-|\mathcal{A}|}{2}}|\boldsymbol{\tilde{\mathcal{K}}}_{S}|^{-\frac{1}{2}}||\mathbf{X}_{\mathcal{A}}^{T}\boldsymbol{\tilde{\mathcal{K}}}_{S}^{-1}\mathbf{X}_{\mathcal{A}}|^{-\frac{1}{2}}\operatorname{exp}\left\\{-\frac{S_{\boldsymbol{\lambda}}^{2}}{2\sigma^{2}}\right\\}$
(2.44)

where $\boldsymbol{\lambda}=(\rho,\gamma,\mathcal{A})^{T}$,
$S_{\boldsymbol{\lambda}}^{2}=\mathbf{y}_{S}^{T}\mathbf{Q}\mathbf{y}_{S}$,
$\mathbf{Q}=\boldsymbol{\tilde{\mathcal{K}}}_{S}^{-1}\mathbf{P}$,
$\mathbf{P}=\mathbf{I}-\mathbf{X}_{\mathcal{A}}(\mathbf{X}_{\mathcal{A}}^{T}\boldsymbol{\tilde{\mathcal{K}}}_{S}^{-1}\mathbf{X}_{\mathcal{A}})^{-1}\mathbf{X}_{\mathcal{A}}^{T}\boldsymbol{\tilde{\mathcal{K}}}_{S}^{-1}$,
and $\mathbf{X}_{\mathcal{A}}$ denotes the corresponding to the set
$\mathcal{A}$ reduced design matrix. In accordance with Section 2.2 we use the
prior

$\displaystyle
p(\sigma^{2},\rho,\gamma)\propto\frac{|\tilde{I}_{R}(\rho,\gamma)|^{\frac{1}{2}}}{\sigma^{2}}$
(2.45)

where

$\displaystyle\mathbf{W}_{\rho}$
$\displaystyle=\frac{\partial\boldsymbol{\tilde{\mathcal{K}}}_{S}}{\partial\xi_{k}}\mathbf{Q}$
(2.46) $\displaystyle\mathbf{W}_{\gamma}$
$\displaystyle=\frac{\partial\boldsymbol{\tilde{\mathcal{K}}}_{S}}{\partial\gamma}\mathbf{Q}$
(2.47) $\displaystyle\tilde{I}_{R}(\boldsymbol{\lambda})$
$\displaystyle=\begin{pmatrix}n-|\mathcal{A}|&\operatorname{tr}(\mathbf{W}_{\rho})&\operatorname{tr}(\mathbf{W}_{\gamma}))\\\
\operatorname{tr}(\mathbf{W}_{\rho})&\operatorname{tr}(\mathbf{W}_{\rho}^{2})&\operatorname{tr}(\mathbf{W}_{\rho}\mathbf{W}_{\gamma})\\\
\operatorname{tr}(\mathbf{W}_{\gamma})&\operatorname{tr}(\mathbf{W}_{\rho}\mathbf{W}_{\gamma})&\operatorname{tr}(\mathbf{W}_{\rho}^{2})\end{pmatrix}.$
(2.48)

The computation of the derivatives presented in Equation 2.46 and Equation
2.47 is given in A.

####  2.3.3 Prior for Random Set

Inspired by Posch et al. [2020] we propose the following prior for
$\mathcal{A}$:

$p(\mathcal{A}=\\{\alpha_{1},...,\alpha_{k}\\})\propto(p_{\alpha_{1}}+...+p_{\alpha_{k}})\frac{1}{k}\tilde{p}(k),$
(2.49)

where

  1. (i)

$\\{p_{\alpha_{1}},...,p_{\alpha_{k}}\\}\subseteq\\{p_{1},...,p_{d}\\}$ with
$\sum_{i=1}^{d}p_{i}=1$ and $p_{i}\geq 0$ for $i=1,...,d$ and

  2. (ii)

$\tilde{p}:\\{1,...,d\\}\rightarrow\mathbb{R}_{0}^{+}$.

The mapping $\tilde{p}$ can be arbitrarily chosen and should represent the a
priori belief in the number of predictors contributing to the model. On the
other hand the parameters $p_{1},...,p_{d}$ are used to represent the a priori
belief in the importance of the predictors $x_{1},...,x_{d}$ to the model.

###  2.4 Model Inference

Clearly, the joint posterior
$p(\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A}|\mathbf{y}_{S},\mathbf{X})$
is analytically intractable. Therefore, we propose a Metropolis within Gibbs
algorithm to simulate from this distribution. The algorithm produces a Markov
Chain with the true posterior as limiting distribution. In particular the
algorithm includes 3 main steps:

  1. 1.

Sample from the conditional posterior
$p(\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1}|\mathbf{y}_{S},\mathbf{X},\sigma_{t}^{2},\gamma_{t},\rho_{t})$
using a classical Metropolis Hastings (MH) algorithm.

  2. 2.

Sample from the conditional posterior
$p(\sigma_{t+1}^{2}|\mathbf{y}_{S},\mathbf{X},\gamma_{t},\rho_{t},\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1})$.

  3. 3.

Sample from the conditional posterior
$p(\gamma_{t+1},\rho_{t+1}|\mathbf{y}_{S},\mathbf{X},\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1},\sigma_{t+1}^{2})$
using Hamiltonian MCMC (HMC).

In order to run the procedure one has to define an initial sample
$\boldsymbol{\beta}_{0},\sigma_{0}^{2},\gamma_{0},\rho_{0},\mathcal{A}_{0}$
that lies in the support of
$p(\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A}|\mathbf{y}_{S},\mathbf{X})$.
Additionally, some tuning parameters for the proposal distribution of
$\mathcal{A}$ are required and also the step size $\varepsilon$ and the number
of steps $L$ for the leap frog method, that is used within the HMC step. In
the following sections details regarding the individual steps are presented.

####  2.4.1 Details Step 1:

Following Posch et al. [2020] to specify the proposal for the random set
$\mathcal{A}$ at first a Bernoulli distributed random variable $c_{h}$ is
introduced

$c_{h}\sim\text{Bernoulli}(p_{h}),$ (2.50)

where $p_{h}\in[0,1]$ can be considered as a tuning parameter. The event
$c_{h}=1$ means that the model size (number of predictors used) changes, i.e.
$k_{t+1}\neq k_{t}$. Moreover, a random variable $\alpha$ with support on the
index set $\\{1,...,d\\}$ is introduced in order to describe the model
transition probabilities in case of changing model size. A realization of the
conditional random variable $\alpha|\mathcal{A}_{t}$ corresponds to the index
of a predictor which is going to be added to or removed from the model.
Transitions between models, which differ by two or more parameters, are not
allowed. The probability mass function of $\alpha|\mathcal{A}_{t}$ is defined
as

$q(\alpha|\mathcal{A}_{t})=\text{I}_{\\{1,...,d\\}}(\alpha)\begin{cases}\tilde{p}_{\alpha}&\text{if
}k_{t}>1,\alpha\notin\mathcal{A}_{t}\\\
\frac{\sum\limits_{i\in\mathcal{A}_{t}}\tilde{p}_{i}}{\tilde{p}_{\alpha}\sum\limits_{i\in\mathcal{A}_{t}}1/\tilde{p}_{i}}&\text{if
}k_{t}>1,\alpha\in\mathcal{A}_{t}\\\ 0&\text{if
}k_{t}=1,\alpha\in\mathcal{A}_{t}\\\
\frac{\tilde{p}_{\alpha}}{\sum\limits_{i\in\\{1,...,p\\}\setminus\mathcal{A}_{t}}\tilde{p}_{i}}&\text{else}\end{cases},$
(2.51)

where the parameters $\tilde{p}_{1},...,\tilde{p}_{d}$ are greater or equal to
zero, sum up to one, and can but must not be identical to the parameters
$p_{1},...,p_{d}$ already used in the prior $p(\mathcal{A})$. Using Equation
2.50 and Equation 2.51 the proposal distribution
$q(\mathcal{A}_{t+1}|\mathcal{A}_{t})$ is finally defined by:

$\displaystyle
q(\mathcal{A}_{t+1}|\mathcal{A}_{t},c_{h}=1)=q([\mathcal{A}_{t+1}\setminus\mathcal{A}_{t}]\cup[\mathcal{A}_{t}\setminus\mathcal{A}_{t+1}]|\mathcal{A}_{t})\text{I}_{\\{1\\}}(|k_{t+1}-k_{t}|),$
(2.52) $\displaystyle
q(\mathcal{A}_{t+1}|\mathcal{A}_{t},c_{h}=0)=\begin{cases}1&\text{if
}\mathcal{A}_{t+1}=\mathcal{A}_{t}\\\ 0&\text{else}\end{cases},$ (2.53)
$\displaystyle
q(\mathcal{A}_{t+1}|\mathcal{A}_{t})=\sum\limits_{c_{h}\in\\{0,1\\}}q(\mathcal{A}_{t+1},c_{h}|\mathcal{A}_{t})$
(2.54) $\displaystyle\hskip
48.36958pt=\sum\limits_{c_{h}\in\\{0,1\\}}q(\mathcal{A}_{t+1}|c_{h},\mathcal{A}_{t})q(c_{h}).$
(2.55)

Simple algebra shows that:

$\displaystyle\boldsymbol{\beta}_{t+1,\mathcal{A}_{t+1}}|\mathbf{y}_{S},\mathbf{X},\sigma_{t}^{2},\gamma_{t},\rho_{t},\mathcal{A}_{t+1}$
$\displaystyle\sim\mathcal{N}(\boldsymbol{\mu}=\mathbf{\Sigma}\mathbf{X}_{\mathcal{A}_{t+1}}^{T}\mathbf{\tilde{C}}_{S}^{-1}\mathbf{y}_{S},\mathbf{\Sigma}=(\mathbf{X}_{\mathcal{A}_{t+1}}^{T}\mathbf{\tilde{C}}_{S}^{-1}\mathbf{X}_{\mathcal{A}_{t+1}})^{-1}),$
(2.56)
$\displaystyle\boldsymbol{\beta}_{t+1,\mathcal{A}_{t+1}^{C}}|\mathbf{y}_{S},\mathbf{X},\sigma_{t}^{2},\gamma_{t},\rho_{t},\mathcal{A}_{t+1}$
$\displaystyle\sim\delta_{0}.$ (2.57)

We use the conditional posterior 2.56-2.57 as proposal for
$\boldsymbol{\beta}_{t+1}$. The Hastings ratio is then given by:

$\displaystyle
r=\frac{L(\boldsymbol{\beta}_{t+1},\sigma_{t}^{2},\gamma_{t},\rho_{t},\mathcal{A}_{t+1}|\mathbf{y}_{S})p(\boldsymbol{\beta}_{t+1}|\mathcal{A}_{t+1})p(\mathcal{A}_{t+1})q(\boldsymbol{\beta}_{t}|\sigma_{t}^{2},\gamma_{t},\rho_{t},\mathcal{A}_{t})q(\mathcal{A}_{t}|\mathcal{A}_{t+1})}{L(\boldsymbol{\beta}_{t},\sigma_{t}^{2},\gamma_{t},\rho_{t},\mathcal{A}_{t}|\mathbf{y}_{S})p(\boldsymbol{\beta}_{t}|\mathcal{A}_{t})p(\mathcal{A}_{t})q(\boldsymbol{\beta}_{t+1}|\sigma_{t}^{2},\gamma_{t},\rho_{t},\mathcal{A}_{t+1})q(\mathcal{A}_{t+1}|\mathcal{A}_{t})}.$
(2.58)

####  2.4.2 Details Step 2:

The conditional posterior
$p(\sigma_{t+1}^{2}|\mathbf{y}_{S},\mathbf{X},\gamma_{t},\rho_{t},\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1})$
can be computed in closed analytical form. Simple algebra shows that it is
given by:

$\displaystyle\sigma_{t+1}^{2}|\mathbf{y}_{S},X,\gamma_{t},\rho_{t},\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1}\sim
IG(\alpha=\frac{n}{2},\beta=\frac{1}{2}(\mathbf{y}_{S}-\mathbf{X}\boldsymbol{\beta})^{T}\boldsymbol{\tilde{\mathcal{K}}}_{S}^{-1}(\mathbf{y}_{S}-\mathbf{X}\boldsymbol{\beta})).$
(2.59)

####  2.4.3 Details Step 3:

For the HMC algorithm we reparameterize the posterior
$p(\gamma_{t+1},\rho_{t+1}|\mathbf{y}_{S},\mathbf{X},\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1},\sigma_{t+1}^{2})$,
such that the support of the reparameterized distribution is unbounded. This
can be easily done using Jacobi transformation:

$\displaystyle\gamma_{t+1}$
$\displaystyle=\frac{1}{1+\operatorname{exp}(-\tilde{\gamma}_{t+1})}\text{
with }\tilde{\gamma}_{t+1}\in\mathbb{R},$ (2.60) $\displaystyle\rho_{t+1}$
$\displaystyle=\operatorname{log}(1+\operatorname{exp}(\tilde{\rho}_{t+1}))\text{
with }\tilde{\rho}_{t+1}\in\mathbb{R},$ (2.61) $\displaystyle
p(\tilde{\gamma}_{t+1},\tilde{\rho}_{t+1}|\mathbf{y}_{S},\mathbf{X},\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1},\sigma_{t+1}^{2})$
$\displaystyle\propto
L(\boldsymbol{\beta}_{t+1},\sigma_{t+1}^{2},\gamma_{t+1},\rho_{t+1},\mathcal{A}_{t+1}|\mathbf{y}_{S})p(\rho_{t+1},\gamma_{t+1})$
$\displaystyle\frac{\operatorname{exp}(\tilde{\rho}_{t+1}-\tilde{\gamma}_{t+1})}{(1+\operatorname{exp}(-\tilde{\gamma}_{t+1})^{2})(1+\operatorname{exp}(\tilde{\rho}_{t+1}))}$
(2.62)

Next we specify the so-called potential energy function:

$\displaystyle E(\tilde{\rho}_{t+1},\tilde{\gamma}_{t+1})$
$\displaystyle=-\operatorname{l}og(p(\tilde{\gamma}_{t+1},\tilde{\rho}_{t+1}|\mathbf{y}_{S},\mathbf{X},\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1},\sigma_{t+1}^{2}))$
(2.63)

Due to Equation 2.62 $E(\cdot|\cdot)$ is known up to an additive constant,
that vanishes if derivatives are computed. Moreover, the kinetic energy is
defined as the negative logarithm of the density of a zero mean normal
distribution with a diagonal covariance matrix:

$\displaystyle
K(\boldsymbol{v}=(v_{\tilde{\rho}},v_{\tilde{\gamma}})^{T})=c+\frac{1}{2}\left(\frac{v_{\tilde{\rho}}^{2}}{m_{\tilde{\rho}}}+\frac{v_{\tilde{\gamma}}^{2}}{m_{\tilde{\gamma}}}\right).$
(2.64)

The parameters $m_{\tilde{\rho}}$ and $m_{\tilde{\gamma}}$ are tuning
parameters of the algorithm and, further, $c$ denotes an additive constant.
Consequently, the partial derivatives of the kinetic energy function are given
by:

$\displaystyle\frac{\partial K(\boldsymbol{v})}{\partial
v_{\tilde{\rho}}}=\frac{v_{\tilde{\rho}}}{m_{\tilde{\rho}}},~{}~{}~{}\frac{\partial
K(\boldsymbol{v})}{\partial
v_{\tilde{\gamma}}}=\frac{v_{\tilde{\gamma}}}{m_{\tilde{\gamma}}}.$ (2.65)

With above considerations the HMC algorithm proposes a sample from the
conditional posterior by applying the following procedure:

  1. 1.

Draw a sample $\boldsymbol{v}=(v_{\tilde{\rho}},v_{\tilde{\gamma}})^{T}$ from
$\mathcal{N}(\mathbf{0},\operatorname{diag}(m_{\tilde{\rho}},m_{\tilde{\gamma}}))$.

  2. 2.

Use the leap-frog method to simulate
$\tilde{\rho}^{*},\tilde{\gamma}^{*},v_{\tilde{\rho}}^{*},v_{\tilde{\gamma}}^{*}$.

    1. (i)

Set
$\boldsymbol{v}^{*}=\boldsymbol{v},\tilde{\rho}^{*}=\tilde{\rho}_{t},\tilde{\gamma}^{*}=\tilde{\gamma}_{t}$.

    2. (ii)

For $j=1,...,L:$

      1. (a)

$v_{\tilde{\rho}}^{*}=v_{\tilde{\rho}}^{*}-\frac{\varepsilon}{2}\frac{\partial
E}{\partial\tilde{\rho}^{*}}(\tilde{\rho}^{*},\tilde{\gamma}^{*}),~{}v_{\tilde{\gamma}}^{*}=v_{\tilde{\gamma}}^{*}-\frac{\varepsilon}{2}\frac{\partial
E}{\partial\tilde{\gamma}^{*}}(\tilde{\rho}^{*},\tilde{\gamma}^{*})$

      2. (b)

$\tilde{\rho}^{*}=\tilde{\rho}^{*}+\varepsilon\frac{v_{\tilde{\rho}}^{*}}{m_{\tilde{\rho}}},\tilde{\gamma}^{*}=\tilde{\gamma}^{*}+\varepsilon\frac{v_{\tilde{\gamma}}^{*}}{m_{\tilde{\gamma}}}$

      3. (c)

$v_{\tilde{\rho}}^{*}=v_{\tilde{\rho}}^{*}-\frac{\varepsilon}{2}\frac{\partial
E}{\partial\tilde{\rho}^{*}}(\tilde{\rho}^{*},\tilde{\gamma}^{*}),~{}v_{\tilde{\gamma}}^{*}=v_{\tilde{\gamma}}^{*}-\frac{\varepsilon}{2}\frac{\partial
E}{\partial\tilde{\gamma}^{*}}(\tilde{\rho}^{*},\tilde{\gamma}^{*})$

  3. 3.

Calculate the Hastings ratio
$r=\operatorname{exp}[K(\boldsymbol{v})+E(\tilde{\rho}_{t},\tilde{\gamma}_{t})-K(\boldsymbol{v}^{*})-E(\tilde{\rho}^{*},\tilde{\gamma}^{*})]$.

  4. 4.

With probability $\operatorname{min}(r,1)$ set
$\tilde{\rho}_{t+1}=\tilde{\rho}^{*},\tilde{\gamma}_{t+1}=\tilde{\gamma}^{*}$.
Otherwise set
$\tilde{\rho}_{t+1}=\tilde{\rho}_{t},\tilde{\gamma}_{t+1}=\tilde{\gamma}_{t}$.

####  2.4.4 Tuning Parameters

This section provides some insights regarding the tuning parameter
specifications used for the evaluation of our approach (Section 3). The tuning
parameter $p_{h}$ is set to the value $0.6$. Thus, the proposal for
$\mathcal{A}$ slightly prefers to modify $\mathcal{A}$ over keeping it
unchanged. Moreover, the parameters $\tilde{p}_{1},...,\tilde{p}_{d}$ are all
specified as $1/d$. Consequently, the proposal does not favor special
predictors. The $\tilde{p}_{i}$ are specified like this, since there is no a
priori knowledge available and a specification using empirical investigations
regarding the importance of the individual predictors is not necessary. The
MCMC chains converge without informative values for the $\tilde{p}_{i}$ after
a reasonable amount of MCMC iterations (a few thousand). For simplicity,
nearly all tuning parameters of the HMC step are fixed and set to the values
$m_{\tilde{\rho}}=m_{\tilde{\gamma}}=1,L=2$. The only exception is the step
size $\varepsilon$ of the leap-frog method. This parameter turns out to be the
most influential one and values between $0.15$ and $0.5$ are chosen. The
actual specification is selected in such a way (per trial and error) that the
MCMC chains show a good mixing.

###  2.5 Prediction

Let $\mathbf{X}^{*}=(\mathbf{x}_{1}^{*}\cdots\mathbf{x}_{v}^{*})^{T}$ be the
design matrix corresponding to some predictor specifications for which
predictions $\widehat{\mathbf{y}}^{*}$ are required. Note that
$\mathbf{y}_{S}$ and $\mathbf{y}^{*}$ are joint normally distributed:

$\displaystyle(\mathbf{y}^{*},\mathbf{y}_{S})^{T}|\mathbf{X}^{*},\mathbf{X},\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A}\sim\mathcal{N}\begin{pmatrix}\begin{pmatrix}\mathbf{X}^{*}\boldsymbol{\beta}\\\
\mathbf{X}\boldsymbol{\beta}\end{pmatrix},\begin{pmatrix}\tilde{\mathbf{C}}_{\mathbf{X}^{*}}&\tilde{\mathbf{C}}_{\mathbf{X}^{*},\mathbf{S}}\\\
\tilde{\mathbf{C}}_{\mathbf{X}^{*},\mathbf{S}}&\tilde{\mathbf{C}}_{\mathbf{S}}\end{pmatrix}\end{pmatrix}.$
(2.66)

This directly implies:

$\displaystyle\mathbf{y}^{*}|\mathbf{y}_{S},\mathbf{X}^{*},\mathbf{X},\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A}\sim\mathcal{N}(\mathbf{X}^{*}\boldsymbol{\beta}+\tilde{\mathbf{C}}_{\mathbf{X}^{*},\mathbf{S}}\tilde{\mathbf{C}}_{\mathbf{S}}^{-1}(\mathbf{y}_{S}-\mathbf{X}\boldsymbol{\beta}),\tilde{\mathbf{C}}_{\mathbf{X}^{*}}-\tilde{\mathbf{C}}_{\mathbf{X}^{*},\mathbf{S}}\tilde{\mathbf{C}}_{\mathbf{S}}^{-1}\tilde{\mathbf{C}}_{\mathbf{X}^{*},\mathbf{S}})$
(2.67)

We use the expected value of the posterior predictive distribution as
prediction. Using Monte Carlo Integration this quantity can easily be
approximated:

$\displaystyle\widehat{\mathbf{y}}^{*}$
$\displaystyle=\int\mathbf{y}^{*}p(\mathbf{y}^{*}|\mathbf{y}_{S},\mathbf{X}^{*},\mathbf{X})~{}d\mathbf{y}^{*}$
(2.68) $\displaystyle=\int\mathbf{y}^{*}\left\\{\int\int\int\int\int
p(\mathbf{y}^{*}|\mathbf{y}_{S},\mathbf{X}^{*},\mathbf{X},\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A})p(\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A}|\mathbf{y}_{S},\mathbf{X})~{}d\sigma^{2}d\gamma
d\rho\\#(d\mathcal{A})(\lambda+\delta_{0})^{d}(d\boldsymbol{\beta})\right\\}~{}d\mathbf{y}^{*}$
(2.69)
$\displaystyle=\int\int\int\int\int\left\\{\int\mathbf{y}^{*}p(\mathbf{y}^{*}|\mathbf{y}_{S},\mathbf{X}^{*},\mathbf{X},\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A})~{}d\mathbf{y}^{*}\right\\}p(\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A}|\mathbf{y}_{S},\mathbf{X})~{}d\sigma^{2}d\gamma
d\rho\\#(d\mathcal{A})(\lambda+\delta_{0})^{d}(d\boldsymbol{\beta})$ (2.70)
$\displaystyle\underset{\ref{pred1}}{=}\int\int\int\int\int\left\\{\mathbf{X}^{*}\boldsymbol{\beta}+\tilde{\mathbf{C}}_{\mathbf{X}^{*},\mathbf{S}}\tilde{\mathbf{C}}_{\mathbf{S}}^{-1}(\mathbf{y}_{S}-\mathbf{X}\boldsymbol{\beta})\right\\}p(\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A}|\mathbf{y}_{S},\mathbf{X})~{}d\sigma^{2}d\gamma
d\rho\\#(d\mathcal{A})(\lambda+\delta_{0})^{d}(d\boldsymbol{\beta})$ (2.71)
$\displaystyle\approx\frac{1}{N}\sum\limits_{i=1}^{N}\mathbf{X}^{*}\boldsymbol{\beta}^{(i)}+\tilde{\mathbf{C}}_{\mathbf{X}^{*},\mathbf{S}}^{(i)}\tilde{\mathbf{C}}_{\mathbf{S}}^{{}^{(i)}-1}(\mathbf{y}_{S}-\mathbf{X}\boldsymbol{\beta}^{(i)}),$
(2.72)

where $\\#$ denotes the counting measure. The $(i)$ in Equation 2.72 is used
to outline that each of the individual summands is computed for an $iid$
sample from the posterior
$p(\boldsymbol{\beta},\sigma^{2},\gamma,\rho,\mathcal{A}|\mathbf{y}_{S},\mathbf{X})$.

##  3 Evaluation

We have implemented the proposed model using the package Rcpp [Eddelbuettel
and François, 2011], which allows for a seamless R and C++ integration.
Additionally, we have used the package RcppArmadillo [Eddelbuettel and
Sanderson, 2014] in order to accelerate R with high-performance linear
algebra. In this section, we evaluate and compare the predictive performance
of our approach (VRNNGP) with the accuracy of some other well-established
feature selection methods. The methods under comparison are:

  * 1.

Lasso. The least absolute shrinkage and selection operator (Lasso)
[Tibshirani, 1996] penalizes the absolute values of the regression
coefficients in the classical linear model in order to obtain a sparse
solution. The R package glmnet [Friedman et al., 2010] implements the lasso.
Cross-validation is used to determine the penalization strength.

  * 2.

Adaptive lasso (ALasso). The adaptive Lasso [Zou, 2006] is an extension of the
classical lasso that allows for individual penalization strengths regarding
the regression coefficients. In this comparison study the individual
penalization factors are defined by multiplying a global factor with the
inverse numbers of the absolute values of the optimized regression
coefficients from a ridge-penalized model. Again, the package glmnet is used
for this task.

  * 3.

Bayes linear selection (BLS). Recently, Posch et al. [2020] proposed a
Bayesian variable selection approach for linear regression models. Since our
approach shares some common features (e.g. the prior for the random set
$\mathcal{A}$) we also include their method in our evaluation study.

  * 4.

Random forest (RF). A combination of tree predictors proposed by Breiman
[2001]. The method is implemented in the R package randomForest [Liaw and
Wiener, 2002].

  * 5.

Robust Gaussian process (RGP). The R package RobustGaSP [Gu et al., 2020]
implements robust Gaussian process regression models using reference priors
and separable (anisotropic) covariance functions. We use this package with the
default specifications (i.e. a Matern $5/2$ kernel etc.), except that we force
the method to estimate a nugget term. Note, that GP models with anisotropic
kernels allow for an automatic relevance determination [Rasmussen and
Williams, 2005]. To this aim, the predictors have to be scaled prior to model
inference.

  * 6.

Maximum likelihood Gaussian process (MLGP). The R package mlegp [Dancik and
Dorman, 2008] implements Gaussian Process regression using maximum likelihood
estimates and an anisotropic covariance function. We use this package with the
default settings for the model inference, except that we specify the mean
function to be non-constant.

The accuracy measures used for the performance comparison are the mean squared
error (MSE) and the mean absolute deviation (MAD). Accuracies are always
computed based on testing datasets that do not share common observations with
the training data. Finally, it should be mentioned that the number of nearest
neighbors is specified as $m=10$ in all subsequent experiments. It turned out
that $10$ neighbors are sufficient to obtain promising results.

###  3.1 Pepelyshev Function

In this section, the Pepelyshev function [Dette and Pepelyshev, 2010], a
commonly used test function in computer experiments, is considered. In this
field it has become common to use a GP in order to model the output of a
complex system [Vollert et al., 2019]. The highly curved Pepelyshev function
is given by

$\displaystyle
f(\mathbf{x})=4(x_{1}-2+8x_{2}-8x_{2}^{2})^{2}+(3-4x_{2})^{2}+16\sqrt{x_{3}+1}(2x_{3}-1)^{2},$
(3.1)

for $x_{1},x_{2},x_{3}\in[0,1]$. To make the approximation problem more
challenging and also in order to obtain a true variable selection task, we
extend Function 3.1 by $17$ inert variables $x_{4},...,x_{20}\in[0,1]$.

The evaluation of our approach is based on a training data set with sample
size $31$ and a testing dataset of size $100$. Both datasets are constructed
by generating a maxmin latin hypercube design (LHD), respectively. Maxmin LHDs
are standard designs for computer experiments [Santner et al., 2003]. They can
be easily constructed by using the R package lhs [Carnell, 2020]. Prior to the
analysis the training and the testing dataset are scaled. In Table 1 the
accuracy of the methods considered in this study is shown. It can be seen that
VRNNGP outperforms the other ones. Reason for the superior performance is that
the other methods fail to identify the actual variables $x_{1},x_{2},x_{3}$,
i.e. they also use the inert variables. However, it should be mentioned that
also VRNNGP is unable to identify $x_{1},x_{2},x_{3}$. In this approach
$x_{2},x_{3}$ have a posterior probability of $1$ to contribute to the models,
while all other predictors have a probability near zero.

Finally, we want to illustrate that the proposed Metropolis within Gibbs
algorithm (Section 2.4) converges fast. To this aim the Figures 4-4 show the
Markov chains of some representative model parameters. It can bee seen that
the algorithm converges after about $1000$ iterations.

Table 1: Accuracy of individual approaches for the Pepelyshev function with inert variables. The proposed approach performs best. Method | MSE | MAD
---|---|---
VRNNGP | 0.1949179 | 0.3148653
MLGP | 1.999522 | 1.176881
RGP | 1.483502 | 1.01951
RF | 0.7351364 | 0.7156615
Lasso | 1.003959 | 0.8617405
ALasso | 1.63588 | 1.065925
BLS | 1.30754 | 0.972762

Figure 1: Markov chain of the cardinality of the random set $\mathcal{A}$.

Figure 2: Markov chain of the covariance parameter $\gamma$.

Figure 3: Markov chain of the covariance parameter $\rho$.

Figure 4: Markov chain of the covariance parameter $\sigma^{2}$.

###  3.2 Sine Function

In this simulation study data are generated from the model

$\displaystyle y$
$\displaystyle=\operatorname{sin}(x_{3})+\operatorname{sin}(5x_{4})+\varepsilon,$
(3.2) $\displaystyle\varepsilon$ $\displaystyle\sim\mathcal{N}(0,0.0025),$
(3.3) $\displaystyle(x_{1},...,x_{20})^{T}$
$\displaystyle\sim\mathcal{N}(\mathbf{0},\mathbf{\Sigma}),$ (3.4)
$\displaystyle\operatorname{diag}(\mathbf{\Sigma})$
$\displaystyle=\mathbf{1},$ (3.5) $\displaystyle\Sigma_{ij}$
$\displaystyle=\begin{cases}0.4&\mbox{ if }(i,j)\in\\{(3,13),(13,3)\\},\\\
0.3&\mbox{ if }(i,j)\in\\{(3,12),(12,3)\\},\\\ 0.4&\mbox{ if
}(i,j)\in\\{(4,14),(14,4)\\},\\\ 0.3&\mbox{ if
}(i,j)\in\\{(4,15),(15,4)\\},\\\ 0&\mbox{ else },\end{cases}$ (3.6)

where $x_{1},x_{2},x_{5},...,x_{20}$ are inert variables. In particular, 100
samples are simulated from model 3.2-3.6 to obtain a dataset, that is
subsequently zero centered and scaled. The performance evaluation is based on
a $5$-fold cross validation. In Table 2 the mean MSEs and the mean MADs of the
comparison methods is given. VRNNGP performs significantly better than the
other methods. Indeed, this method identifies the active variables
$x_{3},x_{4}$ in four of the five considered train/tests splits, while the
other methods are unable to do so.

Table 2: Mean accuracy of individual approaches for the simulation model 3.2-3.6. The proposed approach performs best. Method | MSE | MAD
---|---|---
VRNNGP | 0.3580486 | 0.4316482
MLGP | 1.219693 | 0.9032923
RGP | 0.7396772 | 0.7100051
RF | 0.7706562 | 0.7031102
Lasso | 0.8173408 | 0.7418092
ALasso | 0.8286106 | 0.7483805
BLS | 0.8528274 | 0.7460358

###  3.3 APC data

In this section, real-world data originating from the so-called semiconductor
frontend production is considered [Pleschberger, 2021]. In this production
chain, a wafer, i.e. a slice based on semiconductor material, goes through
hundreds of process steps resulting in thousands of semiconductor devices on
the wafer. During each process step, several sensors mounted on the process
equipment record physical parameters like gas flow, voltage, temperature, etc.
resulting in so-called time-dependent trace data. From this trace data, so-
called key numbers are extracted, which are simple statistical quantities like
mean, standard deviation or maximum used to monitor the process stability.
After frontend production, each device on the wafer is fastidiously tested in
order to meet the customer’s requirements resulting in the so-called wafer
test data. These are analog test procedures measuring quantities like current,
voltage or resistance. Now, if certain wafer test parameters are conspicuous
and lead to potentially high failures, the root cause must be found in
frontend production. The following analysis is intended to identify important
influencing variables in order to identify processes in frontend production
causing the variation at the end. The selected and prepared dataset includes
$228$ observations from $77$ possible predictors (key numbers) extracted from
trace data of two crucial process steps. The target variable is an analog
measurement from the wafer test data describing the wafer quality. The higher
the measured value the worse the quality. Prior to the analysis we have scaled
the predictors and zero centered the target variable. Moreover, the complete
dataset was randomly partitioned into a training and a testing dataset. Both
datasets are equal sized.

In Table 3 the accuracy of the methods considered in this study is shown.
VRNNGP returns the best results followed by RF. MLGP fails completely, while
the robust version RGP performs slightly worse than the linear approaches. In
the Figures 7-7 the importance of the individual predictors obtained from the
methods VRNNGP, RF, and BLS is visualized. For VRNNGP and BLS the importance
of a predictor is the relative frequency how often its index lies in the
posterior samples from $\mathcal{A}$. This can be interpreted as the posterior
probability of the predictor being included in the model. For RF the
visualized importance is the mean decrease in accuracy if a particular
variable is removed. This quantity can be easily computed using the
importance() function included in the R package randomForest. It can be seen
that the importance values differ from method to method. However, the two most
important variables in NNGP (variable $34$, variable $41$) are also important
in RF. In BLS at least the variable $41$ is considered as important. We assume
that the influence of the variable $34$ is non-linear, such that a linear
model fails to detect this influence. Since the variables $34$ and $41$ appear
to be important across models we train the two models RGP and RF again, but
with a reduced dataset including only these two variables. In Table 4 the
respective accuracies are presented. Both models perform comparable, but
significantly better than before. In particular, the reduction to the
variables $34$ and $41$ led to the best results overall. The Figures show the
functions learned by RGP and RF on basis of the reduced dataset. Both models
agree on parameter specifications with a poor production quality. Finally, in
Figure 10 the true production quality of the testing dataset(yTest) is plotted
vs the predicted quality from the RGP model with the reduced dataset. It can
be observed, that all of the really poor production results (yTest $>$ $0.5$)
can be identified. This shows, that the variables selected by NNGP are useful
for this application.

Table 3: Accuracy of individual approaches on the APC data. The proposed approach performs best. Method | MSE | MAD
---|---|---
VRNNGP | 0.2389607 | 0.4197644
MLGP | 32.49251 | 1.44142
RGP | 0.3108199 | 0.4695633
RF | 0.2687667 | 0.4441394
Lasso | 0.2843513 | 0.4568827
ALasso | 0.2903017 | 0.4529153
BLS | 0.2846291 | 0.4471566
Table 4: Accuracy of RGP and RF on reduced dataset. Only the variables $34$ and $41$ are considered. Method | MSE | MAD
---|---|---
RGP | 0.2053222 | 0.3692166
RF | 0.2022327 | 0.3539577

Figure 5: Variable Importance obtained from VRNNGP.

Figure 6: Variable Importance obtained from RF.

Figure 7: Variable Importance obtained from BLS.

Figure 8: Model fit obtained from RGP.

Figure 9: Model fit obtained from RF. Figure 10: True production quality of
the testing dataset (yTest) vs predicted quality. Used model: RGP on the
reduced dataset.

###  3.4 Body Fat Data

In this section, the popular body fat dataset [Johnson, 1996] available in the
R package mplot [Tarr et al., 2018] is considered. The dataset includes 128
observations on 14 variables. The target variable is the body fat percentage
of a man and the predictors are age, weight, height, and ten body
circumference measurements. A 5-fold cross validation is performed to evaluate
and compare the performance of VRNNGP with other methods. Prior to the
analysis that the whole dataset was zero centered and scaled. In Table 5 the
mean MSEs and the mean MADs of the comparison methods is given. VRNNGP
achieves the best results. In the Figures 13-13 the importance of the
individual predictors (compare to Section 3.3) obtained from the methods
VRNNGP, RF, and BLS is visualized. While all of these methods agree that the
sixth predictor (abdomen circumference) is the most influential, there are
some disagreements regarding the order of importance of the other predictors.

Table 5: Mean accuracy of individual approaches for the body fat dataset. The proposed approach performs best. Method | MSE | MAD
---|---|---
VRNNGP | 0.2762897 | 0.4211418
MLGP | 0.3657517 | 0.4865936
RGP | 0.308341 | 0.4434407
RF | 0.3035638 | 0.4458874
Lasso | 0.3113187 | 0.4484568
ALasso | 0.314493 | 0.4521927
BLS | 0.3014364 | 0.4356817

Figure 11: Variable Importance obtained from VRNNGP.

Figure 12: Variable Importance obtained from RF.

Figure 13: Variable Importance obtained from BLS.

###  3.5 Prior Specifications

In the Sections 3.1-3.4 the proposed model is evaluated based on different
datasets. The good performance is not least attributable to the fact that
valuable prior information regarding the random set $\mathcal{A}$ can be
formulated (see Section 2.3.3). Now we provide the individual prior
specifications used in the preceding evaluations. All evaluations share the
definition of the a priori importance $p_{1}=...=p_{d}=1/d$ of the predictors.
No predictor is preferred over another. However, the specification of the
mapping $\tilde{p}(\cdot)$, that represents the a priori belief regarding the
number of predictors contributing to the model, differs among the evaluation
studies. In Section 3.1, in Section 3.2, and in Section 3.4 the mapping is
defined as $\tilde{p}(k)=1/k~{}\forall k\in\\{1,...,d\\}$. Thus, smaller
models are quite intuitively preferred. In Section 3.3 the number of potential
predictors is larger than in the other evaluation studies. Therefore,
$\tilde{p}(\cdot)$ is specified using the recommendations from Posch et al.
[2020], who proposed the prior for $\mathcal{A}$. In particular,
$\tilde{p}(\cdot)$ is defined to be proportional to the third power of the
probability mass function of a zero truncated binomial distribution with size
parameter $d$ and second parameter $1/d$.

##  4 Conclusion

In this work, a novel Bayesian approach to the problem of variable selection
was proposed. To this aim nearest neighbor Gaussian processes, reference
priors, and model selection strategies based on random sets were combined. To
the best of our knowledge, such a combination has not been investigated in
previous works. While the application of the nearest neighbor processes
reduces the computational cost compared to classical processes, the usage of
reference priors guarantees a (numerically) robust model inference. The random
set approach allows for a selection of the most predictive variables, whilst
incorporating appropriate a priori knowledge. For the model inference a
Metropolis within Gibbs algorithm was proposed. The algorithm allows for
transitions between different sets of active predictors (predictors that truly
contribute to the model) and, further, uses Hamiltonian dynamics for the
inference of the Gaussian process correlation parameters. The performance of
the proposed method was extensively evaluated. In evaluations including
simulated data and an approximation problem corresponding to computer
experiments, the new approach could identify truly significant numbers. Purely
linear models, anisotropic Gaussian processes (automatic relevance
determination), and random forest failed at these variable selection tasks.
Finally, a real-world application related to the semiconductor manufacturing
industry showed that the proposed method is of real value for companies and
organizations of many kinds.

## Acknowledgments

The publication is one of the results of the project iDev40 (www.idev40.eu).
This project has received funding from the ECSEL Joint Undertaking (JU) under
grant agreement No 783163. The JU receives support from the European Unions
Horizon 2020 research and innovation programme and Austria (BMVIT-IKT der
Zukunft, FFG project no. 865348), Germany, Belgium, Italy, Spain, Romania. The
information and results set out in this publication are those of the authors
and do not necessarily reflect the opinion of the ECSEL Joint Undertaking.

##  Appendix A Derivatives for Reference Prior

Calculation rules for matrices:

$\displaystyle\frac{\partial\mathbf{AB}}{\partial t}$
$\displaystyle=\frac{\partial\mathbf{A}}{\partial
t}\mathbf{B}+\mathbf{A}\frac{\partial\mathbf{B}}{\partial t}$ (A.1)
$\displaystyle\frac{\partial\mathbf{A}^{-1}}{\partial t}$
$\displaystyle=-\mathbf{A}^{-1}\frac{\partial\mathbf{A}}{\partial
t}\mathbf{A}^{-1}$ (A.2)

$\displaystyle\frac{\partial\boldsymbol{\tilde{\mathcal{K}}}_{S}}{\partial\rho}$
$\displaystyle=\frac{\partial}{\partial\rho}\mathbf{B}_{S}^{-1}\overbrace{\operatorname{diag}(\mathcal{F}_{\mathbf{x}_{1}},...,\mathcal{F}_{\mathbf{x}_{n}})}^{=:\boldsymbol{\mathcal{F}}_{S}}(\mathbf{B}_{S}^{T})^{-1}$
(A.3)
$\displaystyle=\frac{\partial\mathbf{B}_{S}^{-1}\boldsymbol{\mathcal{F}}_{S}}{\partial\rho}(\mathbf{B}_{S}^{T})^{-1}+\mathbf{B}_{S}^{-1}\boldsymbol{\mathcal{F}}_{S}\frac{\partial(\mathbf{B}_{S}^{T})^{-1}}{\partial\rho}$
(A.4)
$\displaystyle=\left(\frac{\partial\mathbf{B}_{S}^{-1}}{\partial\rho}\boldsymbol{\mathcal{F}}_{S}+\mathbf{B}_{S}^{-1}\frac{\partial\boldsymbol{\mathcal{F}}_{S}}{\partial\rho}\right)(\mathbf{B}_{S}^{T})^{-1}+\mathbf{B}_{S}^{-1}\boldsymbol{\mathcal{F}}_{S}\frac{\partial(\mathbf{B}_{S}^{T})^{-1}}{\partial\rho}$
(A.5)
$\displaystyle=\left(-\mathbf{B}_{S}^{-1}\frac{\partial\mathbf{B}_{S}}{\partial\rho}\mathbf{B}_{S}^{-1}\boldsymbol{\mathcal{F}}_{S}+\mathbf{B}_{S}^{-1}\frac{\partial\boldsymbol{\mathcal{F}}_{S}}{\partial\rho}\right)(\mathbf{B}_{S}^{T})^{-1}-\mathbf{B}_{S}^{-1}\boldsymbol{\mathcal{F}}_{S}(\mathbf{B}_{S}^{T})^{-1}\frac{\partial(\mathbf{B}_{S}^{T})}{\partial\rho}(\mathbf{B}_{S}^{T})^{-1}$
(A.6)
$\displaystyle=-\mathbf{B}_{S}^{-1}\left\\{\underbrace{\frac{\partial\mathbf{B}_{S}}{\partial\rho}\mathbf{B}_{S}^{-1}\boldsymbol{\mathcal{F}}_{S}}_{=:\mathbf{A}}-\frac{\partial\boldsymbol{\mathcal{F}}_{S}}{\partial\rho}+\underbrace{\boldsymbol{\mathcal{F}}_{S}(\mathbf{B}_{S}^{T})^{-1}\frac{\partial(\mathbf{B}_{S}^{T})}{\partial\rho}}_{\mathbf{A}^{T}}\right\\}(\mathbf{B}_{S}^{T})^{-1}$
(A.7)
$\displaystyle=-\mathbf{B}_{S}^{-1}\left\\{\mathbf{A}+\mathbf{A}^{T}-\frac{\partial\boldsymbol{\mathcal{F}}_{S}}{\partial\rho}\right\\}(\mathbf{B}_{S}^{T})^{-1}$
(A.8)

$\displaystyle\frac{\partial\mathbf{B}_{S}}{\partial\rho}=\frac{\partial[B_{\mathbf{x}_{i},j}]_{i,j=1,...,n}}{\partial\rho}~{}~{}~{}B_{\mathbf{x}_{i},j}=\begin{cases}1&i=j\\\
-\mathbf{B}_{\mathbf{x}_{i}}[l]&\mathbf{x}_{j}=\mathbf{x}_{i}^{l}\text{ for
some $l$}\\\ 0&\text{else}\end{cases}$ (A.9)

$\displaystyle\frac{\partial\mathbf{B}_{\mathbf{x}_{i}}}{\partial\rho}$
$\displaystyle=\frac{\partial\boldsymbol{\mathcal{K}}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}}{\partial\rho}$
(A.11)
$\displaystyle=\frac{\partial\boldsymbol{\mathcal{K}}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}}{\partial\rho}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}+\boldsymbol{\mathcal{K}}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}\frac{\partial\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}}{\partial\rho}$
(A.12)
$\displaystyle=\frac{\partial\boldsymbol{\mathcal{K}}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}}{\partial\rho}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}-\boldsymbol{\mathcal{K}}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}\frac{\partial\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}}{\partial\rho}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}$
(A.13)
$\displaystyle=\left(\frac{\partial\boldsymbol{\mathcal{K}}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}}{\partial\rho}-\mathbf{B}_{\mathbf{x}_{i}}\frac{\partial\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}}{\partial\rho}\right)\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}$
(A.14)

Thus:

$\displaystyle\frac{\partial\mathbf{B}_{S}}{\partial\rho}$
$\displaystyle=\begin{cases}0&i\geq j\\\
\left\\{\left(\mathbf{B}_{\mathbf{x}_{i}}\frac{\partial\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}}{\partial\rho}-\frac{\partial\boldsymbol{\mathcal{K}}_{\mathbf{x}_{i},N(\mathbf{x}_{i})}}{\partial\rho}\right)\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i})}^{-1}\right\\}[l]&i<j\text{
and }\mathbf{x}_{j}=\mathbf{x}_{i}^{l}\text{ for some }l\\\
0&\text{else}\end{cases}$ (A.15)

$\displaystyle\frac{\partial\boldsymbol{\mathcal{F}}_{S}}{\partial\rho}=\operatorname{diag}\left(\frac{\partial\mathcal{F}_{\mathbf{x}_{1}}}{\rho},...,\frac{\partial\mathcal{F}_{\mathbf{x}_{n}}}{\rho}\right)$
(A.16)

$\displaystyle\frac{\partial\mathcal{F}_{\mathbf{x}_{i}}}{\rho}$
$\displaystyle=\frac{\partial\mathcal{K}_{\mathbf{x}_{i}}-\mathbf{B}_{\mathbf{x}_{i}}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i}),\mathbf{x}_{i}}}{\partial\rho}$
(A.17)
$\displaystyle=\frac{\partial\mathcal{K}_{\mathbf{x}_{i}}}{\partial\rho}-\left(\frac{\partial\mathbf{B}_{\mathbf{x}_{i}}}{\partial\rho}\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i}),\mathbf{x}_{i}}+\mathbf{B}_{\mathbf{x}_{i}}\frac{\partial\boldsymbol{\mathcal{K}}_{N(\mathbf{x}_{i}),\mathbf{x}_{i}}}{\partial\rho}\right)$
(A.18)

All the calculations above hold also for $\mathbf{W}_{\gamma}$ when replacing
$\rho$ with $\gamma$.

$\displaystyle\frac{\partial\mathcal{K}_{\mathbf{x}_{i},\mathbf{x}_{j}}}{\partial\rho}$
$\displaystyle=\gamma\frac{\partial
K_{\mathbf{x}_{i},\mathbf{x}_{j}}}{\partial\rho}$ (A.19)
$\displaystyle\frac{\partial\mathcal{K}_{\mathbf{x}_{i},\mathbf{x}_{j}}}{\partial\gamma}$
$\displaystyle=-\delta_{(\mathbf{x}_{i}=\mathbf{x}_{j})}+K_{\mathbf{x}_{i},\mathbf{x}_{j}}$
(A.20)

$\displaystyle\frac{\partial
K_{\mathbf{x}_{i},\mathbf{x}_{j}}}{\partial\rho}=\left[-\frac{\sqrt{5}d_{\mathcal{A}}}{\rho^{2}}-\frac{10d_{\mathcal{A}}^{2}}{3\rho^{3}}\right]\operatorname{exp}\left(-\frac{\sqrt{5}d_{\mathcal{A}}}{\rho}\right)+K_{\mathbf{x}_{i},\mathbf{x}_{j}}\frac{\sqrt{5}d_{\mathcal{A}}}{\rho^{2}}$
(A.21)

## References

  * Alhamzawi and Taha Mohammad Ali [2018] Alhamzawi, R., Taha Mohammad Ali, H., 2018\.  The bayesian adaptive lasso regression.  Mathematical Biosciences 303, 75 – 82.
  * Berger et al. [2001] Berger, J.O., de Oliveira, V., Sanso, B., 2001.  Objective bayesian analysis of spatially correlated data.  Journal of the American Statistical Association 96, 1361–1374.  URL: http://www.jstor.org/stable/3085905.
  * Breiman [2001] Breiman, L., 2001.  Machine learning, volume 45, number 1 - springerlink.  Machine Learning 45, 5–32.  doi:10.1023/A:1010933404324.
  * Buehlmann et al. [2016] Buehlmann, P., Drineas, P., Kane, M., van der Laan, M., 2016\.  Handbook of Big Data.  Chapman & Hall/CRC.
  * Carnell [2020] Carnell, R., 2020.  lhs: Latin Hypercube Samples.  URL: https://CRAN.R-project.org/package=lhs. r package version 1.1.1.
  * Chen et al. [2007] Chen, T., Morris, J., Martin, E., 2007.  Gaussian process regression for multivariate spectroscopic calibration.  Chemometrics and Intelligent Laboratory Systems 87, 59 – 71.  doi:10.1016/j.chemolab.2006.09.004.
  * Chen and Wang [2010] Chen, T., Wang, B., 2010.  Bayesian variable selection for gaussian process regression: Application to chemometric calibration of spectrometers.  Neurocomputing 73, 2718 – 2726.  doi:10.1016/j.neucom.2010.04.014. pattern Recognition in Bioinformatics Advances in Neural Control.
  * Dancik and Dorman [2008] Dancik, G.M., Dorman, K.S., 2008\.  mlegp: statistical analysis for computer models of biological systems using r.  Bioinformatics 24, 1966–1967.
  * Datta et al. [2016] Datta, A., Banerjee, S., Finley, A.O., Gelfand, A.E., 2016\.  Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets.  Journal of the American Statistical Association 111, 800–812.
  * Dette and Pepelyshev [2010] Dette, H., Pepelyshev, A., 2010\.  Generalized latin hypercube design for computer experiments.  Technometrics 52, 421–429.  doi:10.1198/TECH.2010.09157.
  * Duane et al. [1987] Duane, S., Kennedy, A., Pendleton, B.J., Roweth, D., 1987\.  Hybrid monte carlo.  Physics Letters B 195, 216 – 222.  doi:10.1016/0370-2693(87)91197-X.
  * Eddelbuettel and François [2011] Eddelbuettel, D., François, R., 2011\.  Rcpp: Seamless R and C++ integration.  Journal of Statistical Software 40, 1–18.  doi:10.18637/jss.v040.i08.
  * Eddelbuettel and Sanderson [2014] Eddelbuettel, D., Sanderson, C., 2014\.  Rcpparmadillo: Accelerating r with high-performance c++ linear algebra.  Computational Statistics and Data Analysis 71, 1054–1063.  URL: http://dx.doi.org/10.1016/j.csda.2013.02.005.
  * Fan and Li [2001] Fan, J., Li, R., 2001.  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association 96, 1348–1360.
  * Foster and Stine [2004] Foster, D.P., Stine, R.A., 2004\.  Variable selection in data mining: Building a predictive model for bankruptcy.  Journal of the American Statistical Association 99, 303–313.
  * Friedman et al. [2010] Friedman, J., Hastie, T., Tibshirani, R., 2010.  Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software 33, 1–22.  URL: http://www.jstatsoft.org/v33/i01/.
  * Gu et al. [2020] Gu, M., Palomo, J., Berger, J., 2020.  RobustGaSP: Robust Gaussian Stochastic Process Emulation.  URL: https://CRAN.R-project.org/package=RobustGaSP. r package version 0.6.0.
  * Gu et al. [2018] Gu, M., Wang, X., Berger, J.O., 2018.  Robust gaussian stochastic process emulation.  Ann. Statist. 46, 3038–3066.  URL: https://doi.org/10.1214/17-AOS1648, doi:10.1214/17-AOS1648.
  * Johnson [1996] Johnson, R.W., 1996.  Fitting percentage of body fat to simple body measurements.  Journal of Statistics Education 4\.  doi:10.1080/10691898.1996.11910505.
  * Kazianka and Pilz [2012] Kazianka, H., Pilz, J., 2012\.  Objective bayesian analysis of spatial data with uncertain nugget and range parameters.  The Canadian Journal of Statistics / La Revue Canadienne de Statistique 40, 304–327.  URL: http://www.jstor.org/stable/41724530.
  * Lee et al. [2003] Lee, K.E., Sha, N., Dougherty, E.R., Vannucci, M., Mallick, B.K., 2003.  Gene selection: a bayesian variable selection approach.  Bioinformatics 19 1, 90–7.
  * Liaw and Wiener [2002] Liaw, A., Wiener, M., 2002.  Classification and regression by randomforest.  R News 2, 18–22.  URL: https://CRAN.R-project.org/doc/Rnews/.
  * Linkletter et al. [2006] Linkletter, C., Bingham, D., Hengartner, N., Higdon, D., Ye, K.Q., 2006.  Variable selection for gaussian process models in computer experiments.  Technometrics 48, 478–490.  doi:10.1198/004017006000000228.
  * Neal [1997] Neal, R., 1997.  Monte carlo implementation of gaussian process models for bayesian regression and classification.  arXiv: Data Analysis, Statistics and Probability .
  * Neal [1996] Neal, R.M., 1996.  Bayesian Learning for Neural Networks.  Springer-Verlag, Berlin, Heidelberg.
  * O’Hagan [1978] O’Hagan, A., 1978.  Curve fitting and optimal design for prediction.  Journal of the Royal Statistical Society. Series B (Methodological) 40, 1–42.  URL: http://www.jstor.org/stable/2984861.
  * Park and Casella [2008] Park, T., Casella, G., 2008\.  The bayesian lasso.  Journal of the American Statistical Association 103, 681–686.
  * Paulo [2005] Paulo, R., 2005.  Default priors for gaussian processes.  The Annals of Statistics 33, 556–582.
  * Pleschberger [2021] Pleschberger, M., 2021.  Data set of extracted summary statistics from equipment sensor data.  doi:10.5281/zenodo.4462777. dataset on Zenodo, https://doi.org/10.5281/zenodo.4462777.
  * Posch et al. [2020] Posch, K., Arbeiter, M., Pilz, J., 2020.  A novel bayesian approach for variable selection in linear regression models.  Computational Statistics & Data Analysis 144, 106881.  doi:10.1016/j.csda.2019.106881.
  * Rasmussen and Williams [2005] Rasmussen, C.E., Williams, C.K.I., 2005\.  Gaussian Processes for Machine Learning.  MIT Press.
  * Ročková and George [2014] Ročková, V., George, E.I., 2014\.  Emvs: The em approach to bayesian variable selection.  Journal of the American Statistical Association 109, 828–846.
  * Santner et al. [2003] Santner, T.J., Williams, B.J., Notz, W.I., 2003.  The Design and Analysis of Computer Experiments.  Springer.
  * Tarr et al. [2018] Tarr, G., Müller, S., Welsh, A.H., 2018.  mplot: An R package for graphical model stability and variable selection procedures.  Journal of Statistical Software 83, 1–28.  doi:10.18637/jss.v083.i09.
  * Tibshirani [1996] Tibshirani, R., 1996.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological) 58, 267–288.
  * Vecchia [1988] Vecchia, A.V., 1988.  Estimation and model identification for continuous spatial processes.  Journal of the Royal Statistical Society. Series B (Methodological) 50, 297–312.
  * Vollert et al. [2019] Vollert, N., Ortner, M., Pilz, J., 2019.  Robust additive gaussian process models using reference priors and cut-off-designs.  Applied Mathematical Modelling 65, 586–596.  doi:10.1016/j.apm.2018.07.050.
  * Wu et al. [2017] Wu, X., Chan, L.L.T., Chen, J., Xie, L., 2017\.  Application of gaussian processes with variable shrinkage method and just-in-time modeling in the semiconductor industry, in: 2017 6th International Symposium on Advanced Control of Industrial Processes (AdCONIP), pp. 78–83.  doi:10.1109/ADCONIP.2017.7983759.
  * Yi et al. [2011] Yi, G., Shi, J.Q., Choi, T., 2011.  Penalized gaussian process regression and classification for high-dimensional nonlinear data.  Biometrics 67, 1285–1294.  URL: http://www.jstor.org/stable/41434433.
  * Yuan et al. [2008] Yuan, J., Wang, K., yu, z., Fang, M., 2008.  Reliable multi-objective optimization of high-speed wedm process based on gaussian process regression.  International Journal of Machine Tools and Manufacture 48, 47–60.  doi:10.1016/j.ijmachtools.2007.07.011.
  * Zheng and Zhang [2007] Zheng, H., Zhang, Y., 2007.  Feature selection for high dimensional data in astronomy.  Advances in Space Research 41, 1960 – 1964.
  * Zou [2006] Zou, H., 2006.  The adaptive lasso and its oracle properties.  Journal of the American Statistical Association 101, 1418–1429.
  * Zou and Hastie [2005] Zou, H., Hastie, T., 2005.  Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67, 301–320.