# Gaussian Process Regression and Classification using International
Classification of Disease Codes as Covariates

Sanvesh Srivastava<EMAIL_ADDRESS>Department of Statistics and
Actuarial Science, The University of Iowa  Zongyi Xu<EMAIL_ADDRESS>Department of Statistics and Actuarial Science, The University of Iowa  Yunyi
Li<EMAIL_ADDRESS>McCombs School of Business, University of Texas
Austin  W. Nick Street<EMAIL_ADDRESS>Department of Business Analytics,
The University of Iowa  Stephanie Gilbertson-White stephanie-gilbertson-
<EMAIL_ADDRESS>College of Nursing, The University of Iowa

###### Abstract

International Classification of Disease (ICD) codes are widely used for
encoding diagnoses in electronic health records (EHR). Automated methods have
been developed over the years for predicting biomedical responses using EHR
that borrow information among diagnostically similar patients. Relatively less
attention has been paid on developing patient similarity measures that model
the structure of ICD codes and the presence of multiple chronic conditions,
where a chronic condition is defined as a set of ICD codes. Motivated by this
problem, we first develop a type of string kernel function for measuring
similarity between a pair of subsets of ICD codes, which uses the definition
of chronic conditions. Second, we extend this similarity measure to define a
family of covariance functions on subsets of ICD codes. Using this family, we
develop Gaussian process (GP) priors for Bayesian nonparametric regression and
classification using diagnoses and other demographic information as
covariates. Markov chain Monte Carlo (MCMC) algorithms are used for posterior
inference and predictions. The proposed methods are free of any tuning
parameters and are well-suited for automated prediction of continuous and
categorical biomedical responses that depend on chronic conditions. We
evaluate the practical performance of our method on EHR data collected from
1660 patients at the University of Iowa Hospitals and Clinics (UIHC) with six
different primary cancer sites. Our method has better sensitivity and
specificity than its competitors in classifying different primary cancer sites
and estimates the marginal associations between chronic conditions and primary
cancer sites.

##  1 Introduction

EHR data are useful for organizing information about clinical care and
developing tools for personalized patient care. They are also used for
answering interesting and important research questions that cannot be
practically addressed through traditional prospective research. We develop a
Bayesian model based on GP priors for predicting biomedical responses
depending on the presence of a set of chronic conditions. We assume that the
EHR encode diagnoses using the 10th edition of ICD billing code (hereafter
called ICD code), so a chronic condition is defined using one or more specific
codes from a predefined set of ICD codes [2]. An ICD code is a three to seven
character long string, where the first three characters define the diagnosis
category, the next three define the etiology, anatomic site, and severity, and
the final character identifies encounters. We develop a family of kernel
functions on a pair of ICD code subsets that is based on the definition of ICD
codes and chronic conditions. Using any such kernel as a covariance function
for a GP prior, we perform Bayesian nonparametric regression and
classification using patients’ chronic conditions and demographic information
as covariates.

Encoding diagnoses using ICD coding system has become the most widely adopted
medical diagnostic classification system since it was endorsed by the World
Health Assembly in 1990. Automated extraction of data from the EHR using ICD
codes allows us to obtain large amounts of precise medical information about
chronic conditions. This is a cheap and efficient alternative to the
significantly more expensive and slower manual chart review by a healthcare
provider. Unfortunately, ICD codes are prone to encoding errors, so care must
be exercised in fitting a model with ICD codes as covariates [6]. This has
motivated a significant interest in automating the ICD encoding process to
eliminate any sources of errors [21]. In a properly curated EHR, models that
use ICD codes as covariates in predicting disease risks and related medical
conditions are known to be more accurate than those which do not use ICD codes
[7]. Building on these ideas, we achieve an additional but significant
improvement in predictive accuracy by exploiting the grouping of diagnoses
into chronic conditions and the hierarchical structure of ICD codes. Our
Bayesian methods require minimal tuning and provide results with uncertainty
estimates, so they can be easily used by a healthcare provider for assisting
with patient care.

EHR data have been analyzed using Bayesian nonparametric methods, but the
literature on using ICD codes as covariates is sparsely populated. Henao et
al. [5] have developed a deep hierarchical factor model for identifying latent
subgroups of patients. They summarize the EHR data in terms of counts of
different medications usage, laboratory tests, and diagnoses specified using
ICD codes. If the EHR data only contains information about the presence or
absence of disease, then Bayesian nonparametric extensions of low-rank matrix
factorization perform better in discovering latent disease groups [11].
Recently, deep learning has been widely used for automated regression and
classification using EHR data, but the major focus remains on mining clinical
notes and not on using the information encoded in ICD codes or on providing
uncertainty estimates [15]. A major limitation of all these approaches is that
they ignore the additional structure provided by the membership of diagnoses
in different chronic conditions.

Addressing this gap, we first develop a similarity measure on diagnoses
encoded as ICD codes that accounts for the grouping of diagnoses in chronic
conditions and the hierarchical structure of ICD codes. A widespread practice
is to define the similarity of two patients to be proportional to the number
of common ICD codes in their diagnoses or Jaccard index. Our similarity
measure is superior to this practice in that its value could be positive for a
pair of non-identical ICD codes, and its magnitude depends on the degree of
overlap between the pair. In fact, our similarity measure belongs to the class
of string kernels and extends a restricted version of the _boundrange_ kernel,
which has been widely used in text mining and information retrieval [9, 18].
The similarity between a pair of ICD codes is defined as the total number of
common sub-strings between them, where a sub-string always begins at the first
character of an ICD code. This similarity measure is extended to a pair of
subsets as a weighted average of similarities of all ICD code pairs formed
using the two subsets.

Our second major contribution is to use the similarity between patients for
nonparametric regression and classification using GP priors. This is done in
two stages. First, using the similarity measure, we define the equivalents of
polynomial, exponential, and squared exponential kernels on subsets of ICD
codes, which represent chronic conditions. Second, we use these kernels as
covariance functions for defining GP priors indexed by subsets of ICD codes.
These GP priors are used for fitting Bayesian nonparametric regression and
classification models that are tuned for EHR data analysi [16]. We develop
MCMC algorithms for automated fitting of these models that provide
theoretically guaranteed estimation of posterior uncertainty [4]. Minimal
tuning requirements makes our method extremely attractive for routine
applications in analyzing EHR data and for assisting healthcare providers in
predicting risks related to chronic conditions. We verify our claims through
empirical studies and show that the proposed method provides a more detailed
quantification of dependence of primary cancer sites on chronic conditions.

##  2 Motivating EHR Data

Multimorbidity is the simultaneous occurrence of more than one chronic
condition. Developing operational measures of multimorbidity is an active
research area due to its importance in improving patient care. Calderón-
Larrañaga et al. [2] have operationalized the concept of multimorbidity into a
list of 918 ICD codes that map on to a list of 60 chronic conditions. This
serves as a motivation for developing a classifier that accounts for the
grouping of diagnoses into chronic conditions. We show that this yields better
accuracy in predicting primary cancer sites in UIHC EHR data.

The data for this analysis came from the EHR at UIHC relating to the set of
patients who received cancer treatment beginning in 2017. The cohort included
18 years or older patients who were diagnosed with malignant solid neoplasm of
any type or site. The data included diagnoses of the patients encoded as ICD
codes and the race and marital status of every patient. We filtered patients
who had at least one of the 58 non-cancerous chronic conditions, resulting in
a sample size of 1660. The main question here is to identify chronic
conditions or diagnoses that are predictive of a primary cancer site.
Unfortunately, the Bayesian toolkit provides limited options for answering
such biomedical questions. Using the UIHC EHR data as a representative
example, we develop GP-based regression and classification models that are
tuned for solving such biomedical problems, where the subsets of ICD codes are
covariates.

The main goal is to predict marginal associations between chronic conditions
and cancer primary sites and a minor goal is to predict the primary cancer
site using the ICD codes, which are structured strings; for example, J301
encodes allergy and D50, D51 denote different kinds of anaemia. Clearly, using
dummy variables for representing these ICD codes is inappropriate because D50,
D51 denote the same disease but are assigned different dummy variables.
Motivated from the text mining literature, a useful similarity measure for
strings is defined by counting the number of matching substrings of different
lengths. Accounting for the hierarchical structure of ICD codes, if we modify
this definition by enforcing the substrings to always start at the beginning,
then the similarities of the three pairs (J301, D50), (J301, D51), and (D50,
D51) are 0, 0, 2, respectively, which are more realistic. In the next section,
we develop this idea further to account for the extra structure provided by
the 58 chronic conditions while down-weighting the contributions of commonly
occurring diagnoses.

##  3 Kernel Functions for Diagnoses

###  3.1 Kernel Function on ICD Codes

Consider the structure of ICD codes. Let $t$ be an ICD code, $t^{*}$ be the
longest ICD code, $|t|$ denote the number of characters in $t$, and $t_{1:i}$
be the contiguous substring of length $i$ that starts at the first character
of $t$. The first character of $t$ is $t_{1:1}$ and $t_{1:l}=t$ if $l>|t|$.
Let $\operatorname{\mathcal{T}}$ be the set of all ICD codes,
$\operatorname{\mathcal{A}}$ and $\operatorname{\mathcal{D}}$ be the sets of
letters and numbers that are used to represent an ICD code, and
$\operatorname{\mathcal{B}}=\operatorname{\mathcal{A}}\left(\operatorname{\mathcal{A}}\cup\operatorname{\mathcal{D}}\right)^{|t^{*}|-1}$
is set of all alphanumeric substrings that have an letter as their first
character and have lengths $1,2,\ldots,$ or $|t^{*}|$. Then,
$\operatorname{\mathcal{B}}$ contains all contiguous substrings of ICD codes
in $\operatorname{\mathcal{T}}$ that start at the first character in an ICD
code and
$|\operatorname{\mathcal{T}}|<|\operatorname{\mathcal{B}}|=|\operatorname{\mathcal{A}}|\left(|\operatorname{\mathcal{A}}|+|\operatorname{\mathcal{D}}|\right)^{|t^{*}|-1}$.

The kernel function on two ICD codes is defined as a Euclidean dot product
between their feature maps. Let $\operatorname{{\mathbf{\psi}}}(t)$ be the
feature map of an ICD code $t\in\operatorname{\mathcal{T}}$ and
$\operatorname{\mathcal{F}}=\\{\operatorname{{\mathbf{\psi}}}(t):t\in\operatorname{\mathcal{T}}\\}$
be the feature space. Then, for every $t\in\operatorname{\mathcal{T}}$,
$\operatorname{{\mathbf{\psi}}}(t)\in\\{0,1\\}^{|t^{*}||\operatorname{\mathcal{B}}|}$
and is defined as

$\displaystyle\operatorname{{\mathbf{\psi}}}:\operatorname{\mathcal{T}}\mapsto\operatorname{\mathcal{F}},\quad\operatorname{{\mathbf{\psi}}}(t)=\\{\operatorname{{\mathbf{\psi}}}_{b}(t)\\}_{b\in\operatorname{\mathcal{B}}},\quad\operatorname{{\mathbf{\psi}}}_{b}(t)=\\{1(b=t_{1:1}),\ldots,1(b=t_{1:|t^{*}|})\\}^{\mathrm{\scriptscriptstyle
T}},$ (1)

where $1(\cdot)$ is 1 if $\cdot$ is true and 0 otherwise. The kernel function
on $\operatorname{\mathcal{T}}\times\operatorname{\mathcal{T}}$ is defined
using $\operatorname{{\mathbf{\psi}}}(\cdot)$ and known $|t^{*}|\times|t^{*}|$
symmetric positive semi-definite matrices
$(\operatorname{{\mathbf{\Lambda}}}_{b})_{b\in\operatorname{\mathcal{B}}}$ as

$\displaystyle\kappa_{0}(t,t^{\prime}\mid\operatorname{{\mathbf{\Lambda}}})=\underset{b\in\operatorname{\mathcal{B}}}{\sum}\operatorname{{\mathbf{\psi}}}^{\mathrm{\scriptscriptstyle
T}}_{b}(t)\operatorname{{\mathbf{\Lambda}}}_{b}\operatorname{{\mathbf{\psi}}}_{b}(t^{\prime})\equiv\operatorname{{\mathbf{\psi}}}^{\mathrm{\scriptscriptstyle
T}}(t)\operatorname{{\mathbf{\Lambda}}}\operatorname{{\mathbf{\psi}}}(t^{\prime}),\;t,t^{\prime}\in\operatorname{\mathcal{T}},$
(2)

where
$\operatorname{{\mathbf{\Lambda}}}={\mathrm{diag}}\\{\operatorname{{\mathbf{\Lambda}}}_{b}:b\in\operatorname{\mathcal{B}}\\}$
is a block diagonal matrix. Two useful special cases of (2) are obtained by
setting $\operatorname{{\mathbf{\Lambda}}}_{b}$ to be a diagonal matrix in (2)
such that (i) $\operatorname{{\mathbf{\Lambda}}}_{b}=\lambda_{b}I$ and (ii)
$(\operatorname{{\mathbf{\Lambda}}}_{b})_{aa}=\lambda^{a}$ for every
$b\in\operatorname{\mathcal{B}}$ and $a=1,\ldots,|t^{*}|$, where $I$ is an
identity matrix, $\lambda_{b}>0$ for every $b\in\operatorname{\mathcal{B}}$,
and $\lambda\in(0,1)$. We prove in the next section that
$\kappa_{0}(t,t^{\prime}\mid\operatorname{{\mathbf{\Lambda}}})$ in (2) is a
valid kernel function on
$\operatorname{\mathcal{T}}\times\operatorname{\mathcal{T}}$.

The kernel in (2) belongs to the class of string kernels that are widely used
in text mining and information retrieval [18]. Two popular string kernels are
spectrum and boundrange kernels [9]. The $k$-spectrum kernel defines the
similarity to be the number of matching substrings between two strings of
length $k$ exactly. The $k$-boundrange kernel instead defines the similarity
to be the number of matching substrings between two strings of length less
than or equal to $k$. The substrings can be non-contiguous in both kernels and
the substring matches can be weighted by a factor that decays exponentially
with the substring length. The kernel in (2) is restricted version of
$|t^{*}|$-boundrange kernel where the substrings always start at the beginning
of an ICD code.

###  3.2 Kernel Function on Subsets of ICD Codes

We now define the kernel function on subsets of ICD codes that represent
diagnoses. Let $2^{\operatorname{\mathcal{T}}}$ be the power set of
$\operatorname{\mathcal{T}}$,
$\operatorname{\mathbf{t}}=\\{t_{1},\ldots,t_{r}\\}$ be an element of
$2^{\operatorname{\mathcal{T}}}$,
$\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})$ be the feature map
of $\operatorname{\mathbf{t}}$, and
$\operatorname{\mathcal{F}}_{\operatorname{{\mathbf{\Psi}}}}=\\{\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}):\operatorname{\mathbf{t}}\in
2^{\operatorname{\mathcal{T}}}\\}$ be the feature space of subsets of ICD
codes. The feature vector
$\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})$ and the kernel
function $\kappa_{1}$ on $2^{\operatorname{\mathcal{T}}}\times
2^{\operatorname{\mathcal{T}}}$ are defined using
$\operatorname{{\mathbf{\psi}}}(\cdot)$ in (1) and
$\kappa_{0}(\cdot,\cdot\mid\operatorname{{\mathbf{\Lambda}}})$ in (2),
respectively, as

$\displaystyle\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})=\sum_{{t\in\operatorname{\mathbf{t}}}}\operatorname{{\mathbf{\psi}}}(t),\quad\kappa_{1}(\operatorname{\mathbf{t}},\operatorname{\mathbf{t}}^{\prime}\mid\operatorname{{\mathbf{\Lambda}}})=\sum_{{t\in\operatorname{\mathbf{t}}}}\sum_{{t^{\prime}\in\operatorname{\mathbf{t}}^{\prime}}}\kappa_{0}(t,t^{\prime}\mid\operatorname{{\mathbf{\Lambda}}})=\operatorname{{\mathbf{\Psi}}}^{\mathrm{\scriptscriptstyle
T}}(\operatorname{\mathbf{t}})\operatorname{{\mathbf{\Lambda}}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}^{\prime}),$
(3)

for $\operatorname{\mathbf{t}},\operatorname{\mathbf{t}}^{\prime}\in
2^{\operatorname{\mathcal{T}}}$, where
$\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})$ is a
$(|t^{*}||\operatorname{\mathcal{B}}|)$-dimensional vector with entries in
$\\{0\\}\cup\operatorname{\mathbb{N}}$. The kernel $\kappa_{1}$ is a _derived
subsets_ kernel with the _base_ kernel $\kappa_{0}$; see Definition 9.41 in
Shawe-Taylor and Cristianini [18, Page 317].

In the motivating EHR data, covariates include subsets of ICD codes that have
an additional structure specified by their relation to a chronic condition.
Let $C$ be the total number of chronic conditions and
$\operatorname{\mathbf{T}}_{c}\in 2^{\operatorname{\mathcal{T}}}$ be the
subset of ICD codes that defines the chronic condition $c$ $(c=1,\ldots,C)$. A
patient has chronic condition $c$ if the patient’s diagnoses consists of an
ICD code that belongs to $\operatorname{\mathbf{T}}_{c}$. Accounting for the
$C$ chronic conditions, represent a patient’s diagnoses as the set
$\operatorname{\mathbf{T}}=\\{\operatorname{\mathbf{t}}_{1},\ldots,\operatorname{\mathbf{t}}_{c},\ldots,\operatorname{\mathbf{t}}_{C}\\}\in\left(2^{\operatorname{\mathcal{T}}}\right)^{C}$,
where $\operatorname{\mathbf{t}}_{c}\subseteq\operatorname{\mathbf{T}}_{c}$. A
kernel function on
$\left(2^{\operatorname{\mathcal{T}}}\right)^{C}\times\left(2^{\operatorname{\mathcal{T}}}\right)^{C}$
that accounts for the structure imposed by $C$ chronic conditions is

$\displaystyle\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})=\sum_{c=1}^{C}w_{c}\,\kappa_{1}(\operatorname{\mathbf{t}}_{c},\operatorname{\mathbf{t}}^{\prime}_{c}\mid\operatorname{{\mathbf{\Lambda}}}),\quad\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\in\left(2^{\operatorname{\mathcal{T}}}\right)^{C},\quad
w_{c}>0,$ (4)

where $\operatorname{{\mathbf{\Lambda}}}$ is defined in (2),
$\operatorname{\mathbf{w}}=(w_{1},\ldots,w_{C})$, and $w_{c}$ indicates the
importance of chronic condition $c$ in defining the similarity between two
patients with diagnoses $\operatorname{\mathbf{T}}$ and
$\operatorname{\mathbf{T}}^{\prime}$, respectively. The kernel $\kappa_{2}$
assumes the effect of $C$ chronic conditions on the response is additive, but
other valid methods for combining kernels
$\kappa_{1}(\operatorname{\mathbf{t}}_{1},\operatorname{\mathbf{t}}^{\prime}_{1}\mid\operatorname{{\mathbf{\Lambda}}}),\ldots,\kappa_{1}(\operatorname{\mathbf{t}}_{C},\operatorname{\mathbf{t}}^{\prime}_{C}\mid\operatorname{{\mathbf{\Lambda}}})$
can be also used in (4).

The kernel $\kappa_{2}$ in (4) is affected by the number of elements in
$\operatorname{\mathbf{T}}$ and $\operatorname{\mathbf{T}}^{\prime}$. Let
$\|\cdot\|$ be the Euclidean norm. The value of $\kappa_{1}$ in (3) is
impacted by $\operatorname{\mathbf{t}}$ in that
$\|\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})\|$ increases with
the cardinality of $\operatorname{\mathbf{t}}$. This implies that $\kappa_{2}$
in (4) is large if the cardinality of $\operatorname{\mathbf{T}}$ is large;
therefore, the length of feature maps impacts the similarity measure instead
of the chronic conditions. This undesirable effect is removed by normalizing
the feature maps of
$\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}$ before
computing the dot product, which results in the kernel

$\displaystyle\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\frac{\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})}{\left\\{\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{1/2}\left\\{\kappa_{2}(\operatorname{\mathbf{T}}^{\prime},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{1/2}}.$
(5)

All kernels in this paper are derived from $\kappa$ in (5), which depends on
the kernels $\kappa_{0}$, $\kappa_{1}$, and $\kappa_{2}$ defined earlier.

The following theorem proves that $\kappa$ is a valid kernel function on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$
and defines a metric and the equivalents of polynomial, exponential, and
squared exponential (SE) kernels on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$,
respectively.

######  Theorem 3.1.

Let
$\\{\operatorname{{\mathbf{\Lambda}}}_{b}\\}_{b\in\operatorname{\mathcal{B}}}$
and $\operatorname{\mathbf{w}}$ be defined as in (2) and (4),
$\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\in(2^{\operatorname{\mathcal{T}}})^{C}$,
$\sigma^{2}>0$, $\phi>0$, and

$d(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})=\left\\{\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})+\kappa(\operatorname{\mathbf{T}}^{\prime},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})-2\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{1/2}.$

Then, given $\operatorname{{\mathbf{\Lambda}}}$ and
$\operatorname{\mathbf{w}}$,

  1. 1.

$\kappa(\cdot,\cdot\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
and
$d(\cdot,\cdot\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
are valid kernel and distance functions on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$;

  2. 2.

if
$\operatorname{{\mathbf{\Lambda}}}=\operatorname{\mathbf{R}}^{\top}\operatorname{\mathbf{R}}$,
where $\operatorname{\mathbf{R}}$ is an upper triangular matrix, then the
feature map of $\kappa$ is

$\displaystyle\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})=\frac{\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})}{\|\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\|},\quad\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})=\\{w_{1}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}_{1}),\ldots,w_{C}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}_{C})\\}^{\mathrm{\scriptscriptstyle
T}},$ (6)

where
$\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
is a $(C|t^{*}||\operatorname{\mathcal{B}}|)$-dimensional vector; and

  3. 3.

the equivalents of polynomial kernel of order $s$ and $\gamma$-exponential
kernel on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$,
respectively, are defined as

$\displaystyle\kappa_{s}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\sigma^{2},\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\sigma^{2}\left\\{1+\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{s},\quad
s\in\operatorname{\mathbb{N}},$
$\displaystyle\kappa_{\gamma}^{e}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\sigma^{2},\phi,\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\sigma^{2}e^{-\phi\,\\{d(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\\}^{\gamma}},\quad\gamma\in[1,2],$
(7)

where $\sigma^{2}$ and $\phi$ play the role of variance and inverse length-
scale parameters, respectively, and $\kappa_{\gamma}^{e}$ reduces to the
exponential and SE kernels when $\gamma$ equals 1 and 2, respectively.

The proof of this theorem is in the supplementary material along with other
proofs. The polynomial kernel has finite dimensional feature. If $s>1$, then
the polynomial kernel $\kappa_{s}$ is obtained by a non-linear embedding of
the kernel $\kappa$. The dimension of the feature space for $\kappa_{s}$ is
${C|t^{*}||\operatorname{\mathcal{B}}|+s\choose s}$ [Proposition 9.2, 18]. The
feature space of $\kappa_{\gamma}^{e}$ lies in a Hilbert spaces of functions,
where the “smoothness” of a function depends on the chosen kernel; see Section
4.4 for the precise definitions and greater details. We use the SE kernel in
our experiments due to its popularity.

(a) Kernel matrix of $\kappa_{2}$.

(b) Kernel matrix of $\kappa$.

(c) Exponential kernel matrix.

(d) Radial basis function kernel matrix.

(e) Spectrum kernel matrix.

(f) Boundrange kernel matrix.

Figure 1: Heatmaps of kernel matrices for patients in the UIHC EHR data. The
color palette represents 0 and 1 using dark blue and dark red colors,
respectively. The first four kernel matrices in (a)–(d) set
$\operatorname{{\mathbf{\Lambda}}}=\operatorname{\mathbf{I}}$, $\sigma^{2}=1$,
$\phi=1$, and
$w_{c}\propto\kappa_{1}(\operatorname{\mathbf{T}}_{c},\operatorname{\mathbf{T}}_{c})$
for kernels $\kappa_{2},\kappa,\kappa_{1}^{e}$, and $\kappa_{2}^{e}$, where
$\operatorname{\mathbf{T}}_{c}$ is the set of ICD codes defining the $c$th
chronic condition. The last two kernel matrices are obtained using spectrum
and boundrange kernels, which are popular string kernels in text mining, with
substring length equalling 3. The block correlation structures are better
captured in (b), (c), and (d) than (a), (e), and (f), indicating the
superiority of the proposed set of kernel functions.

We illustrate the application of six kernel functions in capturing the
similarity of patients in the UIHC EHR data, where we expect patients with
common primary sites to have similar sets of chronic conditions (Figure 1).
Let $\operatorname{\mathbf{T}}_{i}$ be the ICD codes for the patent $i$, $n$
be the sample size, $k(\cdot,\cdot)$ be a kernel function on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$.
Then, the kernel matrix defined using $k$ is an $n\times n$ matrix with
$(i,j)$-th entry
$k(\operatorname{\mathbf{T}}_{i},\operatorname{\mathbf{T}}_{j})$. To highlight
the similarity of patients, we have grouped the patients into six blocks
depending on their primary cancer sites. The kernel matrix with $k=\kappa_{2}$
in (4) fails to capture the similarities of patients with common primary sites
(Figure 1(a)). On the other hand, the kernel matrix with $k=\kappa$ in (5),
which is the normalized version of $\kappa_{2}$, captures the patient
similarities, where the diagonal blocks denote patient with common primary
sites and the off-diagonal blocks capture similarities of patients with
cancers at two different sites (Figure 1(b)). The same observation is true for
the $\gamma$-exponential kernel in (3.1) with $\gamma$ equals 1 and 2,
respectively (Figures 1(c) and 1(d)). The diagonal and off-diagonal blocks are
clearest for radial basis function kernel. The spectrum and boundrange kernel
matrices fail to capture similarities of patients except those with cancer at
brain and nervous system, where the diagnoses are very similar compared to
cancers at other sites (Figures 1(e) and 1(f)).

##  4 Modeling With Subsets of ICD Codes as Covariates

###  4.1 Setup

Consider the setup for regression and classification problems with ICD codes
as covariates. Let $n$ be the sample size and
$(y_{i},\operatorname{\mathbf{x}}_{i},\operatorname{\mathbf{T}}_{i})$ be the
data for subject $i$ ($i=1,\ldots,n$), where $y_{i}$ is the response,
$\operatorname{\mathbf{x}}_{i}^{\mathrm{\scriptscriptstyle
T}}=(x_{i1},\ldots,x_{ip})$ is the vector of $p$-dimensional covariates
excluding ICD codes with $x_{i1}=1$, and
$\operatorname{\mathbf{T}}_{i}=\\{\operatorname{\mathbf{t}}_{i1},\ldots,\operatorname{\mathbf{t}}_{iC}\\}\in(2^{\operatorname{\mathcal{T}}})^{C}$.
For every $i$, $\operatorname{\mathbf{t}}_{ic}$ denotes the diagnoses of
subject $i$ that belong to chronic condition $c$ and
$\operatorname{\mathbf{t}}_{ic}$ can be an empty set $(c=1,\ldots,C)$. In
regression and classification problems, $y_{i}\in\operatorname{\mathbb{R}}$
and $y_{i}\in\\{0,1\\}$, respectively, and
$(x_{i2},\ldots,x_{ip})^{{\mathrm{\scriptscriptstyle
T}}}\in\operatorname{\mathbb{R}}^{p-1}$ for every $i$. We also predict the
responses and estimate the covariate effects for a given collection of
covariates including ICD codes
$(\operatorname{\mathbf{x}}_{j}^{*},\operatorname{\mathbf{T}}_{j}^{*})$
$(j=1,\ldots,m)$, where
$\operatorname{\mathbf{T}}_{j}^{*}\in(2^{\operatorname{\mathcal{T}}})^{C}$.

We model the effects of $\operatorname{\mathbf{x}}_{i}$ and
$\operatorname{\mathbf{T}}_{i}$ independently. Let
$\operatorname{{\mathbf{\beta}}}\in\operatorname{\mathbb{R}}^{p}$ and
$f:(2^{\operatorname{\mathcal{T}}})^{C}\mapsto\operatorname{\mathbb{R}}$.
Then, the covariate effects of $\operatorname{\mathbf{x}}_{i}$ and
$\operatorname{\mathbf{T}}_{i}$ are
$\operatorname{\mathbf{x}}^{{\mathrm{\scriptscriptstyle
T}}}_{i}\operatorname{{\mathbf{\beta}}}$ and
$f(\operatorname{\mathbf{T}}_{i})$, respectively, and the overall covariate
effect for subject $i$ as
$\operatorname{\mathbf{x}}^{{\mathrm{\scriptscriptstyle
T}}}_{i}\operatorname{{\mathbf{\beta}}}+f(\operatorname{\mathbf{T}}_{i})$. In
practice, $\operatorname{{\mathbf{\beta}}}$ and $f$ are unknown. Choosing a
Bayesian approach, we estimate the posterior distributions of
$\operatorname{{\mathbf{\beta}}}$ and $f$ given
$\\{(y_{i},\operatorname{\mathbf{x}}_{i},\operatorname{\mathbf{T}}_{i})\\}_{i=1}^{n}$.
We use a parametric model for the estimating the effects of
$\operatorname{\mathbf{x}}_{i}$ because our focus in on non-parametric
regression and classification using
$\operatorname{\mathbf{T}}\in(2^{\operatorname{\mathcal{T}}})^{C}$..

###  4.2 Regression Using Subsets of ICD Codes as Covariates

Consider the nonparametric regression model. Let $\epsilon_{i}$ denote the
idiosyncratic error for subject $i$, $\epsilon_{i}$s are independent and
identically distributed (iid) with mean 0, and $\tau^{2}>0$ denotes the error
variance. Then, the model assumes that

$\displaystyle y_{i}=\operatorname{\mathbf{x}}_{i}^{\mathrm{\scriptscriptstyle
T}}\operatorname{{\mathbf{\beta}}}+f(\operatorname{\mathbf{T}}_{i})+\epsilon_{i},\quad\epsilon_{i}\overset{\text{iid}}{\sim}N(0,\tau^{2}),\quad
i=1,\ldots,n,$ (8)

where $N(a,b)$ denotes the Gaussian density with mean $a$ and variance $b$. We
put a GP prior on $f(\cdot)$ with 0 mean function and covariance function
$\kappa_{f}(\cdot,\cdot\mid\operatorname{{\mathbf{\theta}}})$ on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$,
where $\operatorname{{\mathbf{\theta}}}$ are the parameters specifying
$\kappa_{f}$. The model is completed by putting a prior distribution on
$(\operatorname{{\mathbf{\beta}}},\tau^{2},\operatorname{{\mathbf{\theta}}})$
and its form depends on the choice of $\kappa_{f}$.

We develop posterior inference algorithms for $\kappa_{f}$ equalling
$\kappa_{\gamma}^{e}$ defined in Theorem 3.1. We exclude the polynomial kernel
because it yields a finite dimensional regression model, which is inflexible
compared to the non-parametric regression model obtained using
$\kappa_{\gamma}^{e}$. Setting $\operatorname{\mathbf{w}}=(1,\ldots,1)$ and
$\operatorname{{\mathbf{\Lambda}}}_{b}$ to be an identity matrix
$(b\in\operatorname{\mathcal{B}})$ in (6) and (3.1), we define

$\displaystyle\kappa_{\text{exp}}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime})=\kappa_{1}^{e}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}),\quad\kappa_{\text{SE}}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime})=\kappa_{2}^{e}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}),\quad\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\in(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$
(9)

The parameter for $\kappa_{\text{exp}}$ and $\kappa_{\text{SE}}$ is
$\operatorname{{\mathbf{\theta}}}=\\{\sigma^{2},\phi\\}$, where
$\sigma^{2},\phi$ are positive scalars. Denote the kernel matrix as
$\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}$ with
$\\{\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}\\}_{ij}=\kappa_{\text{exp}}(\operatorname{\mathbf{T}}_{i},\operatorname{\mathbf{T}}_{j})$
or
$\\{\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}\\}_{ij}=\kappa_{\text{SE}}(\operatorname{\mathbf{T}}_{i},\operatorname{\mathbf{T}}_{j})$,
where $i,j\in\\{1,\ldots,n\\}$. Let
$\operatorname{\mathbf{y}}=(y_{1},\ldots,y_{n})^{\top}$,
$\operatorname{\mathbf{X}}$ be the $n\times p$ matrix with
$\operatorname{\mathbf{x}}_{i}^{\mathrm{\scriptscriptstyle T}}$ as its $i$th
row,
$\operatorname{\mathbf{f}}=\\{f(\operatorname{\mathbf{T}}_{1}),\ldots,f(\operatorname{\mathbf{T}}_{n})\\}^{\mathrm{\scriptscriptstyle
T}}$, and
$\operatorname{{\mathbf{\epsilon}}}=(\epsilon_{1},\ldots,\epsilon_{n})^{\mathrm{\scriptscriptstyle
T}}$. Then, the hierarchical Bayesian model for $\operatorname{\mathbf{y}}$
based on (8) and parameter-expansion is defined as

$\displaystyle\operatorname{\mathbf{y}}=\operatorname{\mathbf{X}}\operatorname{{\mathbf{\beta}}}+\operatorname{\mathbf{f}}+\operatorname{{\mathbf{\epsilon}}},\quad\operatorname{{\mathbf{\epsilon}}}\mid\sigma^{2},\alpha\sim
N(\operatorname{\mathbf{0}},\sigma^{2}\alpha\operatorname{\mathbf{I}}),\quad(\operatorname{{\mathbf{\beta}}},\sigma^{2})\propto\sigma^{-2},\quad\operatorname{\mathbf{f}}\mid\operatorname{{\mathbf{\theta}}}\sim
N(\operatorname{\mathbf{0}},\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}),$
(10)

where $\operatorname{{\mathbf{\beta}}}$ and $\operatorname{\mathbf{f}}$ are
assumed to be independent apriori, $\alpha=\tau^{2}/\sigma^{2}$ is the inverse
of signal-to-noise ratio, and $\phi\in(a_{\phi},b_{\phi})$, where
$0<a_{\phi}<b_{\phi}$. We assign a prior on $(\phi,\alpha)$ through
$(u_{1},u_{2})\in\operatorname{\mathbb{R}}^{2}$ as

$\displaystyle\phi=a_{\phi}+(b_{\phi}-a_{\phi})/(1+e^{-u_{1}}),\quad\alpha=e^{u_{2}},\quad(u_{1},u_{2})^{\top}\sim
N\\{\operatorname{\mathbf{0}},{\mathrm{diag}}(b_{1},b_{2})\\},$ (11)

where $b_{1},b_{2}>0$. The parameters $\sigma^{2},\phi$ are non-identified in
$\kappa_{\text{exp}},\kappa_{\text{SE}}$, but this does not affect the
inference on $f$ or prediction of the response [16].

The MCMC algorithm for posterior inference and predictions in (10) is a Gibbs-
type sampling algorithm, with a step that draws $(\phi,\alpha)$ using
Elliptical Slice Sampling (ESS) [12]. Let
$\operatorname{\mathbf{f}}^{*}=\\{f(\operatorname{\mathbf{T}}_{1}^{*}),\ldots,f(\operatorname{\mathbf{T}}^{*}_{m})\\}^{\mathrm{\scriptscriptstyle
T}}$ be the value of $f$ and
$\operatorname{\mathbf{y}}^{*}=(y_{1}^{*},\ldots,y_{m}^{*})^{\mathrm{\scriptscriptstyle
T}}$ be the response at
$(\operatorname{\mathbf{x}}_{j}^{*},\operatorname{\mathbf{T}}_{j}^{*})$
$(j=1,\ldots,m)$. For notational convenience, define
$\operatorname{\mathbf{K}}(\phi)$, $\operatorname{\mathbf{K}}_{*}$, and
$\operatorname{\mathbf{K}}_{**}$ to be the matrices satisfying
$\\{\operatorname{\mathbf{K}}(\phi)\\}_{ii^{\prime}}=\kappa_{f}(\operatorname{\mathbf{T}}_{i},\operatorname{\mathbf{T}}_{i^{\prime}})/\sigma^{2}$,
$\\{\operatorname{\mathbf{K}}_{*}\\}_{ij}=\kappa_{f}(\operatorname{\mathbf{T}}_{i},\operatorname{\mathbf{T}}^{*}_{j})/\sigma^{2}$,
and
$\\{\operatorname{\mathbf{K}}_{**}\\}_{jj^{\prime}}=\kappa_{f}(\operatorname{\mathbf{T}}_{j}^{*},\operatorname{\mathbf{T}}^{*}_{j^{\prime}})/\sigma^{2}$,
where $i,i^{\prime}\in\\{1,\ldots,n\\}$ and $j,j^{\prime}\in\\{1,\ldots,m\\}$,
and
$\operatorname{\mathbf{C}}_{yy}=\operatorname{\mathbf{K}}(\phi)+\alpha\operatorname{\mathbf{I}}$.
If $\operatorname{\mathbf{L}}$ is a lower triangular matrix such that
$\operatorname{\mathbf{C}}_{yy}=\operatorname{\mathbf{L}}\operatorname{\mathbf{L}}^{\top}$,
then define

$\displaystyle\tilde{\operatorname{\mathbf{y}}}=\operatorname{\mathbf{L}}^{-1}\operatorname{\mathbf{y}},\quad\tilde{\operatorname{\mathbf{X}}}=\operatorname{\mathbf{L}}^{-1}\operatorname{\mathbf{X}},\quad\tilde{\operatorname{{\mathbf{\eta}}}}=\operatorname{\mathbf{L}}^{-1}\operatorname{{\mathbf{\eta}}},\quad\hat{\operatorname{{\mathbf{\beta}}}}=(\tilde{\operatorname{\mathbf{X}}}^{\mathrm{\scriptscriptstyle
T}}\tilde{\operatorname{\mathbf{X}}})^{-1}\tilde{\operatorname{\mathbf{X}}}^{\top}\tilde{\operatorname{\mathbf{y}}},\quad\hat{\tilde{\operatorname{\mathbf{y}}}}=\tilde{\operatorname{\mathbf{X}}}\hat{\operatorname{{\mathbf{\beta}}}}.$
(12)

The MCMC algorithm for posterior inference on
$\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{f}}^{*},\operatorname{\mathbf{K}}_{**}$
and for predicting $\operatorname{\mathbf{y}}^{*}$ cycles through the
following four steps until convergence to its stationary distribution.

  1. 1.

Draw $(\sigma^{2},\operatorname{{\mathbf{\beta}}})$ given
$\operatorname{\mathbf{y}},\phi,\alpha$ as

$\displaystyle\sigma^{2}\mid\operatorname{\mathbf{y}},\phi,\alpha\sim\frac{\|\tilde{\operatorname{\mathbf{y}}}-\hat{\tilde{\operatorname{\mathbf{y}}}}\|^{2}}{\chi^{2}_{n-p}},\quad\operatorname{{\mathbf{\beta}}}\mid\sigma^{2},\operatorname{\mathbf{y}},\phi,\alpha\sim
N\\{\hat{\operatorname{{\mathbf{\beta}}}},\sigma^{2}(\tilde{\operatorname{\mathbf{X}}}^{\mathrm{\scriptscriptstyle
T}}\tilde{\operatorname{\mathbf{X}}})^{-1}\\},$ (13)

where $\chi^{2}_{n-p}$ is the $\chi^{2}$ random variable with $n-p$ degrees of
freedom.

  2. 2.

Draw $(u_{1},u_{2})$ given
$\operatorname{\mathbf{y}},\sigma^{2},\operatorname{{\mathbf{\beta}}}$ are
drawn using ESS with the prior $(u_{1},u_{2})\sim
N\\{\operatorname{\mathbf{0}},{\mathrm{diag}}(b_{1},b_{2})\\}$ and the
likelihood

$\displaystyle\ell(u_{1},u_{2})=\frac{1}{2\pi\sigma^{2}[\text{det}\\{\operatorname{\mathbf{K}}(\phi)+\alpha\operatorname{\mathbf{I}}\\}]^{1/2}}e^{-\tfrac{1}{2\sigma^{2}}(\operatorname{\mathbf{y}}-\operatorname{\mathbf{X}}\beta)^{\mathrm{\scriptscriptstyle
T}}\\{\operatorname{\mathbf{K}}(\phi)+\alpha\operatorname{\mathbf{I}}\\}^{-1}(\operatorname{\mathbf{y}}-\operatorname{\mathbf{X}}\beta)}.$
(14)

Define $\phi$ and $\alpha$ given $(u_{1},u_{2})$ as in (11) and set
$\tau^{2}=\sigma^{2}\alpha$.

  3. 3.

Draw $\operatorname{\mathbf{f}}^{*}$ given
$\operatorname{\mathbf{y}},\sigma^{2},\operatorname{{\mathbf{\beta}}},\phi,\alpha$
from $N(\operatorname{\mathbf{m}}_{*},\operatorname{\mathbf{V}}_{*})$, where

$\displaystyle\operatorname{\mathbf{m}}_{*}=\operatorname{\mathbf{K}}_{*}^{\top}\operatorname{\mathbf{C}}_{yy}^{-1}(\operatorname{\mathbf{y}}-\operatorname{\mathbf{X}}\operatorname{{\mathbf{\beta}}}),\quad\operatorname{\mathbf{V}}_{*}$
$\displaystyle=\sigma^{2}\left(\operatorname{\mathbf{K}}_{**}-\operatorname{\mathbf{K}}_{*}^{\top}\operatorname{\mathbf{C}}_{yy}^{-1}\operatorname{\mathbf{K}}_{*}\right).$
(15)

  4. 4.

Draw $\operatorname{\mathbf{y}}^{*}$ given
$\operatorname{\mathbf{X}}^{*},\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{f}}^{*},\tau^{2}$
from
$N(\operatorname{\mathbf{X}}^{*}\operatorname{{\mathbf{\beta}}}+\operatorname{\mathbf{f}}^{*},\tau^{2}\operatorname{\mathbf{I}})$.

We collect MCMC draws of
$\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{f}}^{*},\operatorname{\mathbf{y}}^{*},\sigma^{2},\phi,\tau^{2}$
post convergence. The derivation of this algorithm is in the supplementary
materials along with the others. The algorithm in steps (1)–(4) a variant of
Gibbs sampling algorithm for posterior inference in univariate spatial linear
models in that we replace the Metropolis-Hastings step by an ESS step in (2)
and $\operatorname{\mathbf{T}}$ replaces a spatial location [1]. Unlike the
Metropolis-Hastings step, the ESS step is free of any proposal tuning, which
is preferred in automated applications.

###  4.3 Classification Using Subsets of ICD Codes as Covariates

The classification model is based on logistic regression. It assumes that

$\displaystyle\log\,\frac{\Pr(y_{i}=1)}{\Pr(y_{i}=0)}=\operatorname{\mathbf{x}}_{i}^{\mathrm{\scriptscriptstyle
T}}\operatorname{{\mathbf{\beta}}}+f(\operatorname{\mathbf{T}}_{i}),\quad
y_{i}\in\\{0,1\\},\quad i=1,\ldots,n,$ (16)

where $y_{i}$ is the $i$th response. The choice of priors on
$\operatorname{{\mathbf{\beta}}},f(\cdot)$ and the kernels remain the same as
in (10), but the prior on
$\operatorname{{\mathbf{\theta}}}=(\phi,\sigma^{2})^{\top}$ is assigned
through $(u_{1},u_{2})\in\operatorname{\mathbb{R}}^{2}$ as

$\displaystyle\phi=a_{\phi}+(b_{\phi}-a_{\phi})/(1+e^{-u_{1}}),\quad\sigma^{2}=e^{u_{2}},\quad(u_{1},u_{2})^{\top}\sim
N\\{\operatorname{\mathbf{0}},{\mathrm{diag}}(b_{1},b_{2})\\}.$ (17)

where $b_{1},b_{2}>0$. We use Pólya-Gamma data augmentation (PG-DA) for
posterior inference on
$\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{f}}^{*},\operatorname{\mathbf{K}}_{**}$
and prediction of $\operatorname{\mathbf{y}}^{*}$ [13]. This setup ensures
that the MCMC algorithms for inference and predictions in classification and
regression models are very similar. Specifically, we introduce Pólya-Gamma
random variables $\omega_{1},\ldots,\omega_{n}$ specific to every observation
such that $\omega_{i}$ are marginally distributed as PG(1, 0), where PG is the
Pólya-Gamma distribution with parameters $b=1$ and $c=0$, respectively. Define
$\bar{\operatorname{\mathbf{y}}}=(y_{1}-1/2,\ldots,y_{n}-1/2)^{{\mathrm{\scriptscriptstyle
T}}}$, $n\times(p+n)$ matrix
$\operatorname{\mathbf{A}}=[\operatorname{\mathbf{X}}\;\operatorname{\mathbf{I}}]$,
$(p+n)\times 1$ vector
$\operatorname{\mathbf{b}}=(\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{f}})^{\top}$,
$n\times 1$ vector
$\operatorname{{\mathbf{\omega}}}=(\omega_{1},\ldots,\omega_{n})^{\mathrm{\scriptscriptstyle
T}}$, $n\times n$ diagonal matrix
$\operatorname{{\mathbf{\Omega}}}={\mathrm{diag}}(\operatorname{{\mathbf{\omega}}})$,
pseudo responses
$\operatorname{\mathbf{z}}=\operatorname{{\mathbf{\Omega}}}^{-1}\bar{\operatorname{\mathbf{y}}}$.
Then, the conditional likelihood of $\operatorname{\mathbf{b}}$ given
$\operatorname{\mathbf{y}},\operatorname{{\mathbf{\omega}}}$ and the
associated model are defined as

$\displaystyle\ell(\operatorname{{\mathbf{\eta}}}\mid\operatorname{\mathbf{y}},\operatorname{{\mathbf{\omega}}})\propto
e^{-\frac{1}{2}(\operatorname{\mathbf{A}}\operatorname{\mathbf{b}}-\operatorname{\mathbf{z}})^{\top}\Omega(\operatorname{\mathbf{A}}\operatorname{\mathbf{b}}-\operatorname{\mathbf{z}})},\quad\operatorname{\mathbf{z}}=\operatorname{\mathbf{A}}\operatorname{\mathbf{b}}+\operatorname{{\mathbf{\epsilon}}}=\operatorname{\mathbf{X}}\operatorname{{\mathbf{\beta}}}+\operatorname{\mathbf{f}}+\operatorname{{\mathbf{\epsilon}}},\quad\operatorname{{\mathbf{\epsilon}}}\sim
N(0,\operatorname{{\mathbf{\Omega}}}^{-1}),$ (18)

respectively. Furthermore, Theorem 1 in Polson et al. [13] implies that
$\omega_{i}$ given $\operatorname{\mathbf{b}}$ and the $i$th row of $A$,
denoted as $\operatorname{\mathbf{a}}_{i}^{\mathrm{\scriptscriptstyle T}}$,
follows PG(1, $|\operatorname{\mathbf{a}}_{i}^{\mathrm{\scriptscriptstyle
T}}\operatorname{\mathbf{b}}|$).

The MCMC algorithm for posterior inference and predictions in (16) follows
from arguments similar to those used in Section 4.2. Marginalizing over
$\operatorname{\mathbf{f}}$ in (18) implies that $\operatorname{\mathbf{z}}$
given $\operatorname{\mathbf{y}}$, $\operatorname{{\mathbf{\beta}}}$,
$\operatorname{{\mathbf{\omega}}}$, $\phi,\sigma^{2}$ is distributed as
$N(\operatorname{\mathbf{X}}\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{C}}_{zz})$,
where
$\operatorname{\mathbf{C}}_{zz}=\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}+\operatorname{{\mathbf{\Omega}}}^{-1}$
and it plays the same role as $\operatorname{\mathbf{C}}_{yy}$ in (10).
Following (12), if $\operatorname{\mathbf{L}}$ is a lower triangular matrix
such that
$\operatorname{\mathbf{C}}_{zz}=\operatorname{\mathbf{L}}\operatorname{\mathbf{L}}^{\top}$,
then define

$\displaystyle\tilde{\operatorname{\mathbf{z}}}=\operatorname{\mathbf{L}}^{-1}\operatorname{\mathbf{z}},\quad\tilde{\operatorname{\mathbf{X}}}=\operatorname{\mathbf{L}}^{-1}\operatorname{\mathbf{X}},\quad\hat{\operatorname{{\mathbf{\beta}}}}=(\tilde{\operatorname{\mathbf{X}}}^{\mathrm{\scriptscriptstyle
T}}\tilde{\operatorname{\mathbf{X}}})^{-1}\tilde{\operatorname{\mathbf{X}}}^{\top}\tilde{\operatorname{\mathbf{z}}}.$
(19)

The MCMC algorithm for posterior inference on
$\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{f}}^{*},\operatorname{\mathbf{K}}_{**}$
and for predicting $\operatorname{\mathbf{y}}^{*}$ cycles through the
following six steps until convergence to its stationary distribution.

  1. 1.

Draw $\operatorname{{\mathbf{\beta}}}$ given
$\operatorname{\mathbf{y}},\operatorname{{\mathbf{\omega}}},\phi,\sigma^{2}$
as
$N\\{\hat{\operatorname{{\mathbf{\beta}}}},(\tilde{\operatorname{\mathbf{X}}}^{\mathrm{\scriptscriptstyle
T}}\tilde{\operatorname{\mathbf{X}}})^{-1}\\}$.

  2. 2.

Draw $(u_{1},u_{2})$ given
$\operatorname{\mathbf{y}},\operatorname{{\mathbf{\omega}}},\operatorname{{\mathbf{\beta}}}$
using ESS following (14), where the likelihood and prior for $(u_{1},u_{2})$
are defined in (18) and (17). Define $\phi$ and $\sigma^{2}$ given
$(u_{1},u_{2})$ using (17).

  3. 3.

Draw $\operatorname{\mathbf{f}}$ given
$\operatorname{\mathbf{y}},\operatorname{{\mathbf{\omega}}},\operatorname{{\mathbf{\beta}}},\phi,\sigma^{2}$
as $N(\operatorname{\mathbf{m}},\operatorname{\mathbf{V}})$, where

$\displaystyle\operatorname{\mathbf{V}}=\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}-\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}\operatorname{\mathbf{C}}_{zz}^{-1}\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}},\quad\operatorname{\mathbf{m}}=\operatorname{\mathbf{V}}\operatorname{{\mathbf{\Omega}}}(\operatorname{\mathbf{z}}-\operatorname{\mathbf{X}}\operatorname{{\mathbf{\beta}}}).$
(20)

  4. 4.

Draw $\operatorname{\mathbf{f}}^{*}$ given
$\operatorname{\mathbf{y}},\operatorname{{\mathbf{\beta}}},\sigma^{2},\phi$
from $N(\operatorname{\mathbf{m}}_{*},\operatorname{\mathbf{V}}_{*})$, where

$\displaystyle\operatorname{\mathbf{m}}_{*}=\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}*}^{\top}\operatorname{\mathbf{C}}_{zz}^{-1}(\operatorname{\mathbf{z}}-\operatorname{\mathbf{X}}\operatorname{{\mathbf{\beta}}}),\quad\operatorname{\mathbf{V}}_{*}$
$\displaystyle=\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}**}-\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}*}^{\top}\operatorname{\mathbf{C}}_{zz}^{-1}\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}*},$
(21)

where
$\\{\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}*}\\}_{ij}=\kappa_{f}(\operatorname{\mathbf{T}}_{i},\operatorname{\mathbf{T}}^{*}_{j})$,
and
$\\{\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}**}\\}_{jj^{\prime}}=\kappa_{f}(\operatorname{\mathbf{T}}_{j}^{*},\operatorname{\mathbf{T}}^{*}_{j^{\prime}})$,
where $i\in\\{1,\ldots,n\\}$ and $j,j^{\prime}\in\\{1,\ldots,m\\}$.

  5. 5.

Draw $y^{*}_{i}$ given
$\operatorname{\mathbf{X}}^{*},\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{f}}^{*}$
from Bernoulli($p_{i}^{*}$) for $i=1,\ldots,m$, where
$p^{*}_{i}=e^{v_{i}}/(1+e^{v_{i}})$ and
$v_{i}=\operatorname{\mathbf{x}}_{i}^{*\top}\operatorname{{\mathbf{\beta}}}+(\operatorname{\mathbf{f}}^{*})_{i}$.

  6. 6.

Draw $\omega_{i}$ given
$\operatorname{\mathbf{X}},\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{f}}$
from
PG$\\{1,|\operatorname{\mathbf{x}}_{i}^{\top}\operatorname{{\mathbf{\beta}}}+(\operatorname{\mathbf{f}})_{i}|\\}$
for $i=1,\ldots,n$.

We collect MCMC draws of
$\operatorname{{\mathbf{\beta}}},\operatorname{\mathbf{f}}^{*},\operatorname{\mathbf{y}}^{*},\sigma^{2},\phi$
post convergence.

Our MCMC algorithm is similar to other algorithms based on the PG-DA strategy.
First, Polson et al. [13, Section 5, Table 3] present an application of the
PG-DA strategy for nonparametric regression with a negative binomial response
and a GP prior, but a linear predictor similar to
$\operatorname{\mathbf{x}}_{i}^{\top}\operatorname{{\mathbf{\beta}}}$ in (8)
is absent in their model. Second, Wang and Roy [19] develop a sampling
algorithm based on the PG-DA strategy for posterior inference in a Bayesian
logistic linear mixed model with independent and Gaussian random effects and
prove its geometric ergodicity. Our MCMC algorithm is similar to theirs in
that $f(\cdot)$ plays the role of random effects and the GP prior on
$f(\cdot)$ induces dependence among patients that have similar diagnoses.
Overall, the MCMC algorithm above broadens the range of application of the PG-
DA strategy.

###  4.4 Theoretical Properties

The convergence rates of the posterior distributions of $f(\cdot)$ in Sections
4.2 and 4.3 follow from the general theoretical setup for GP regression and
classification Ghosal and van der Vaart [4, Section 11.4.4]. Our focus is on
the stationary GPs with the Matérn and SE covariance kernels because
$\kappa_{\gamma}^{e}$ is based on them. For stationary GP priors, the
posterior convergence rates are described using tail-decay properties of the
_spectral measure_ associated with the covariance kernel and _regularity_ of
the true $f(\cdot)$. In this section, we adapt existing results for our
feature maps
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}$
assuming that $\operatorname{\mathbf{w}}$ and $\operatorname{\mathbf{R}}$ are
known.

Consider the Matérn GP prior and its spectral measure. Following Rasmussen and
Williams [16, Page 84], if $d=C|t^{*}||\operatorname{\mathcal{B}}|$, then the
Matérn covariance kernel with smoothness and range parameters $\nu$ and
$\phi$, respectively, is defined through its spectral measure $\lambda$ on
$\operatorname{\mathbb{R}}^{d}$ as

$\displaystyle\kappa_{\text{Mat}}(\operatorname{\mathbf{h}}\mid\nu,\phi)$
$\displaystyle=\int_{\operatorname{\mathbb{R}}^{d}}e^{-\iota
2\pi\operatorname{\mathbf{h}}^{\top}\operatorname{\mathbf{s}}}d\lambda(\operatorname{\mathbf{s}}),\quad\operatorname{\mathbf{h}}\in\operatorname{\mathbb{R}}^{d},\quad\nu>0,\quad\phi>0,\quad\iota=\sqrt{-1},$
$\displaystyle d\lambda(\operatorname{\mathbf{s}})$
$\displaystyle=\frac{2^{d}\pi^{d/2}\Gamma(\nu+d/2)}{\Gamma(\nu)}(2\nu\phi^{2})^{\nu}\left(2\nu\phi^{2}+4\pi^{2}\|\operatorname{\mathbf{s}}\|^{2}\right)^{-\nu-d/2}d\operatorname{\mathbf{s}},\quad\operatorname{\mathbf{s}}\in\operatorname{\mathbb{R}}^{d},$
$\displaystyle\kappa_{\text{Mat}}(\operatorname{\mathbf{h}}\mid\nu,\phi)$
$\displaystyle=\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\phi\|\operatorname{\mathbf{h}}\|\right)^{\nu}\operatorname{\mathcal{K}}_{\nu}\left(\sqrt{2\nu}\phi\|\operatorname{\mathbf{h}}\|\right),$
(22)

where $\operatorname{\mathcal{K}}_{\nu}$ is the modified Bessel function of
second kind. If $\nu\in\operatorname{\mathbb{N}}$, then sample paths of the
Matérn GP have partial derivatives up to order $\nu$ and all of them are
square integrable.

The GP prior in Theorem 3.1 is indexed by the feature vectors
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\in[0,1]^{d}$
for $\operatorname{\mathbf{T}}\in(2^{\operatorname{\mathcal{T}}})^{C}$. This
is equivalent to embedding the string $\operatorname{\mathbf{T}}$ in
$[0,1]^{d}$ using the feature map. The GP covariance kernels are defined using
this embedding. The equivalent of Matérn covariance kernel on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$
with variance $\sigma^{2}$ is defined through
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}},\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}^{\prime}}$
as
$\sigma^{2}\kappa_{\text{Mat}}(\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}}-\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}^{\prime}}\mid\nu,\phi)$,
where $\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}}$
denotes
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$.
If $\nu=1/2$, then this kernel equals $\kappa_{1}^{e}$. Taking point-wise
limit of
$\sigma^{2}\kappa_{\text{Mat}}(\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}}-\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}^{\prime}}\mid\nu,\phi)$
as $\nu\rightarrow\infty$, we obtain the equivalent of SE covariance kernel on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$
as
$\sigma^{2}\exp({-\phi^{2}\|\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}}-\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}^{\prime}}\|^{2}/2})$.

We now define the regularity of functions using three function spaces. Let
$\lfloor\nu\rfloor$ denote the largest integer less than or equal to $\nu$ and
$\check{g}(\operatorname{\mathbf{s}})$ denote the Fourier transform of
$g:\operatorname{\mathbb{R}}^{d}\rightarrow\operatorname{\mathbb{R}}$. First,
the Hölder space $\mathfrak{C}^{\nu}[0,1]^{d}$ is the space of all functions
on $[0,1]^{d}$ whose partial derivatives of order $(k_{1},\ldots,k_{d})$ exist
for all nonnegative integers $k_{1},\ldots,k_{d}$ such that
$k_{1}+\ldots+k_{d}\leq\lfloor\nu\rfloor$ and whose $\lfloor\nu\rfloor$-th
derivative is Lipschitz of order $\nu-\lfloor\nu\rfloor$. Second, the Sobolev
space $\mathfrak{H}^{\nu}[0,1]^{d}$ is space of all functions on $[0,1]^{d}$
that are the restrictions of a function $g$ for which
$\int(1+\|\operatorname{\mathbf{s}}\|^{2})^{\nu}|\check{g}(\operatorname{\mathbf{s}})|^{2}d\operatorname{\mathbf{s}}$
is finite. Finally, the space of infinitely smooth functions on $[0,1]^{d}$
are restrictions of a function $g$ on $[0,1]^{d}$ for which
$\int\exp(\gamma\|\operatorname{\mathbf{s}}\|^{r})|\check{g}(\operatorname{\mathbf{s}})|^{2}d\operatorname{\mathbf{s}}$
is finite and is denoted as $\mathfrak{A}^{r,\gamma}[0,1]^{d}$ with $r\geq 1$
and $\gamma>0$. A function is $\nu$-regular and $\infty$-regular on
$[0,1]^{d}$ if it belongs to
$\mathfrak{C}^{\nu}[0,1]^{d}\cap\mathfrak{H}^{\nu}[0,1]^{d}$ and
$\mathfrak{A}^{r,\gamma}[0,1]^{d}$ for some $r\geq 1$ and $\gamma>0$,
respectively.

We introduce some notations for stating the theoretical results. Reformulating
the models in Sections 4.2 and 4.3 in terms of feature maps
$\Psi\in[0,1]^{d}$, the population versions of (8) and (16) after setting
$\operatorname{{\mathbf{\beta}}}=\operatorname{\mathbf{0}}$ are

$\displaystyle y=\mu_{0}(\Psi)+\epsilon,\;\epsilon\sim
N(0,\tau^{2}_{0}),\;\tau^{2}_{0}>0,\quad\log\,\frac{\Pr(y=1)}{\Pr(y=0)}=\mu_{0}(\Psi),\quad$
(23)

respectively, where $\tau_{0}^{2}$ is the error variance, any
$\operatorname{\mathbf{T}}\in(2^{\operatorname{\mathcal{T}}})^{C}$ is embedded
in $[0,1]^{d}$ as
$\Psi=\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$,
and the effect of $\operatorname{\mathbf{T}}$ on $y$ is modeled through the
feature map
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$,
which is fixed for a given
$\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}$. The model is (23) is
well-defined because $\Psi$ is defined in Theorem 3.1 and is fixed if
$\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}$ are known. The true
relation between the feature map and response is denoted as $\mu_{0}(\cdot)$,
which maps $[0,1]^{d}$ to $\operatorname{\mathbb{R}}$. The training data
$y_{1},\ldots,y_{n}$ are independently generated following (23) in the two
models for fixed
$\operatorname{\mathbf{T}}_{1},\ldots,\operatorname{\mathbf{T}}_{n}$ that have
been determined apriori. Every $\operatorname{\mathbf{T}}_{i}$ is embedded in
$[0,1]^{d}$ as
$\Psi_{i}=\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}_{i})$
and (23) is satisfied for all the $n$ embeddings. The true distribution of
$y_{i}$ is $N\\{\mu_{0}(\Psi_{i}),\tau_{0}^{2}\\}$ and
Bernoulli[$1/\\{1+e^{-\mu_{0}(\Psi_{i})}\\}$] in (8) and (16), respectively.
The $n$-fold product of these distributions is denoted as
$\operatorname{\mathbb{P}}_{0}^{n}$ in both models and
$\operatorname{\mathbf{y}}\sim\operatorname{\mathbb{P}}_{0}^{n}$. The GP prior
on $f(\cdot)$ in (8) and (16) is induced using the Matérn and SE covariance
kernels if $\mu_{0}(\cdot)$ in (23) has a finite or infinite regularity,
respectively. The posterior distributions of $f$ given
$\operatorname{\mathbf{y}}$ are denoted as $\Pi_{\text{Mat},n}^{\text{Reg}}$,
$\Pi_{\text{SE},n}^{\text{Reg}}$ and $\Pi_{\text{Mat},n}^{\text{Clas}}$,
$\Pi_{\text{SE},n}^{\text{Clas}}$ in (8) and (16) depending on the GP prior.

The convergence rates are described in terms of posterior probabilities
assigned to small neighborhoods of $\mu_{0}$ as $n\rightarrow\infty$. Using
(23) and the embedding of $(2^{\operatorname{\mathcal{T}}})^{C}$ in
$[0,1]^{d}$ based on the feature map
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}$,
the distance between $f(\cdot)$ and $\mu_{0}(\cdot)$ is defined as the
empirical norm

$\displaystyle\|f-\mu_{0}\|_{n}=\left(\frac{1}{n}\sum_{i=1}^{n}\left[f\left\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}_{i})\right\\}-\mu_{0}\left\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}_{i})\right\\}\right]^{2}\right)^{1/2},$
(24)

where $\|\cdot\|_{n}$ conditions on the covariates
$\operatorname{\mathbf{T}}_{1},\ldots,\operatorname{\mathbf{T}}_{n}$ and its
dependence on $\operatorname{\mathbf{w}}$ and $\operatorname{\mathbf{R}}$ is
suppressed. Given $\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}$, the
true and estimated effects of $\operatorname{\mathbf{T}}_{i}$ on $y_{i}$ are
$\mu_{0}\left\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}_{i})\right\\}$
and
$f\left\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}_{i})\right\\}$,
respectively. The next theorem quantifies the rate of convergence of the
posterior distribution in terms of $\|\cdot\|_{n}$ norm.

######  Theorem 4.1.

Assume that the regression and classification models defined in (23) are true,
$\sigma$ and $\phi$ are known, and the Matérn covariance kernel has smoothness
$\nu_{1}$.

  1. 1.

If $\mu_{0}$ is $\nu_{2}$-regular in the sense that
$\mu_{0}\in\mathfrak{C}^{\nu_{2}}[0,1]^{d}\cap\mathfrak{H}^{\nu_{2}}[0,1]^{d}$
with $\nu_{2}\leq\nu_{1}$, then as $n\rightarrow\infty$ there are positive
constants $C_{1},C_{2}$ such that

$\displaystyle\Pi_{\text{Mat},n}^{\text{Reg}}\left\\{f:\|f-\mu_{0}\|_{n}>C_{1}\epsilon_{n}\mid\operatorname{\mathbf{y}}\right\\}\rightarrow
0,\quad\Pi_{\text{Mat},n}^{\text{Clas}}\left\\{f:\|f-\mu_{0}\|_{n}>C_{2}\epsilon_{n}\mid\operatorname{\mathbf{y}}\right\\}\rightarrow
0$ (25)

in $\operatorname{\mathbb{P}}_{0}^{n}$-probability for every
$\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}$ with
$\epsilon_{n}=n^{-\nu_{2}/(2\nu_{1}+d)}$.

  2. 2.

If $\mu_{0}$ is infinitely smooth functions in the sense that
$\mu_{0}\in\mathfrak{A}^{r,\gamma}[0,1]^{d}$ for some $r\geq 1$ and
$\gamma>0$, then as $n\rightarrow\infty$ there are positive constants
$C_{1},C_{2}$ such that

$\displaystyle\Pi_{\text{SE},n}^{\text{Reg}}\left\\{f:\|f-\mu_{0}\|_{n}>C_{1}\epsilon_{n}\mid\operatorname{\mathbf{y}}\right\\}\rightarrow
0,\quad\Pi_{\text{SE},n}^{\text{Clas}}\left\\{f:\|f-\mu_{0}\|_{n}>C_{2}\epsilon_{n}\mid\operatorname{\mathbf{y}}\right\\}\rightarrow
0$ (26)

in $\operatorname{\mathbb{P}}_{0}^{n}$-probability for every
$\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}$ with
$\epsilon_{n}=n^{-1/2}(\log n)^{\max(1/r,1/2+d/4)}$.

The proof of Theorem 4.1 is based on the standard setup [4, Theorems 11.23 and
11.22]. The theorem is also applicable to the case when
$\operatorname{{\mathbf{\beta}}}$ is non-zero if we absorb
$\operatorname{\mathbf{x}}$ as a part of the
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}}$
vector. Specifically, we scale the covariates
$\operatorname{\mathbf{x}},\operatorname{\mathbf{x}}^{\prime}$ so that they
lie in $[0,1]^{p}$, redefine $d$ to be
$p+C|t^{*}||\operatorname{\mathcal{B}}|$, and set the Matérn and SE kernels to
be
$\sigma^{2}\kappa_{\text{Mat}}(\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}}-\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}^{\prime}}\mid\nu,\phi)\cdot\kappa_{\text{Mat}}(\operatorname{\mathbf{x}}-\operatorname{\mathbf{x}}^{\prime}\mid\nu,\phi)$
and
$\sigma^{2}\exp({-\phi^{2}\|\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}}-\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{T}}^{\prime}}\|^{2}/2})\cdot\exp({-\phi^{2}\|\operatorname{\mathbf{x}}-\operatorname{\mathbf{x}}^{\prime}\|^{2}/2})$,
respectively, for any
$(\operatorname{\mathbf{x}},\operatorname{\mathbf{T}}),(\operatorname{\mathbf{x}}^{\prime},\operatorname{\mathbf{T}}^{\prime})$
in $[0,1]^{p}\times(2^{\operatorname{\mathcal{T}}})^{C}$.

##  5 Experiments

###  5.1 Setup

We evaluate the performance of GP-based classification in (16) using three
simulation studies and an analysis of the UIHC EHR data. In this section,
covariance kernel of the GP is SE and the focus is on classification because
our main goal is to address questions relevant to the UIHC EHR data;
supplementary materials contain the results for GP-based regression using (8).
In the simulated and real data analyses, the cutoff for predicting the
response is estimated using Receiver Operating Characteristics (ROC) curve. We
compare the predictive performance of the all the methods using accuracy, area
under the ROC curve (AUC), sensitivity, and specificity on the test data.

We compare GP-based classification with parametric and nonparametric models.
The parametric competitors are logistic regression and its regularized
versions that penalize the regression coefficients using the lasso and ridge
penalties. The covariates in these models are dummy variables that indicate
the presence of ICD codes in a patient’s diagnosis. Because the SE kernel in
(5) is a restriction of the boundrange kernel, we use support vector machine
(SVM) and kernel ridge regression (KRR) with the spectrum and boundrange
kernels as the nonparametric competitors. We also use the kernel function in
(5) for defining feature vectors. Set $\phi=1,\sigma^{2}=1$ in the kernel
matrix $\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}$. If
$\operatorname{\mathbf{U}}\operatorname{\mathbf{D}}\operatorname{\mathbf{U}}^{\top}$
is the singular value decomposition of $\operatorname{\mathbf{K}}$, then the
first $r$ columns of
$\operatorname{\mathbf{U}}\operatorname{\mathbf{D}}^{1/2}$ are the feature
vectors that are used as covariates in random forest, logistic regression, and
its penalized extensions, where $r$ is selected using the scree plot.

The experiments are performed in R. We used glmnet [3] for regularized
logistic regression, ranger [20] for random forest, and kernlab [8] for SVM
and KRR. The performance metrics are evaluated using the pROC package [17].
The tuning parameters in all the methods are selected using the recommended
settings. We use the spectrum and boundrange kernels in kernlab with 3 and 4
as the length of the matching substring, respectively. This means that the
matches in the two kernels include the first three and four contiguous
characters in a pair of ICD codes, which are also the most informative [2];
therefore, we expect the results for kernlab and our method to be very
similar. We use ESS by setting $a_{\phi}$, $b_{\phi}$, $b_{1}$, and $b_{2}$ in
(17) to be 0,5, 1, and 2, respectively. The MCMC algorithms run for 10,000
iterations. The draws from the first 5,000 iterations are discarded as burn-
ins, and every fifth draw in the remaining chain is chosen for posterior
inference and prediction.

###  5.2 Simulated data analysis

We present three simulation studies. The first two are based on simple
parametric models, whereas the third is based on a nonparametric model based
on the UIHC EHR data. We expect all methods to have comparable performances in
the first two simulations and differ significantly in the third one. The three
simulations are replicated ten times.

_First two simulations._ Both simulations have four chronic conditions. The
first chronic condition is denoted using the codes in $\\{A,B,AA,BB,BA,AB\\}$,
where every code denotes a “diagnosis”. Replacing $(A,B)$ with $(C,D)$,
$(E,F)$, and $(G,H)$ in this set gives the set of codes for the second, third,
and fourth chronic conditions, respectively. The chronic conditions are
further structured into two groups: the first group includes the first and
third chronic conditions and the second group includes the remaining two. We
assign nine codes to every patient, where six out of the nine codes are from
the 12 codes defining first or second group of chronic conditions and the
remaining three come from the other chronic condition group. The probability
of $y=1$ is higher than that of $y=0$ if codes from the second chronic
condition group are in majority and vice-versa.

The first simulation study is based on logistic regression. Let $z_{ij}^{(c)}$
be a dummy variable indicating the presence of the $j$th code in the $c$th
chronic condition in the $i$th patient ($c=1,\ldots,4$; $j=1,\ldots,6$;
$i=1,\ldots,n+m$), where codes follow the dictionary order, $n$ is the
training data size, and $m$ testing data size. The first and second group
dummy variables for the $i$th patient are
$\\{z_{ij}^{(c)}:j=1,\ldots,6;c=1,3\\}$ and
$\\{z_{ij}^{(c)}:j=1,\ldots,6;c=2,4\\}$. The response $y_{i}$ is simulated
independently from Bernoulli ($p_{i}$), where

$\displaystyle\log\frac{p_{i}}{1-p_{i}}=0.1+0.2x_{i}+\sum_{c=1}^{4}\sum_{j=1}^{6}(-1)^{c}z_{ij}^{(c)},\quad
x_{i}\overset{\text{ind.}}{\sim}N(0,1)$ (27)

for $i=1,\ldots,n+m$. The second simulation is a slight modification of the
first. Dropping the
$\operatorname{\mathbf{x}}_{i}^{\top}\operatorname{{\mathbf{\beta}}}$ term in
(27), we modify it to only include the interaction of $\\{AB,BA\\}$,
$\\{CD,DC\\}$, $\\{EF,FE\\}$, and $\\{GH,HG\\}$, respectively, as

$\displaystyle\log\frac{p_{i}}{1-p_{i}}=2\left(z_{i\text{EF}}^{(2)}\cdot
z_{i\text{FE}}^{(2)}+z_{i\text{GH}}^{(4)}\cdot
z_{i\text{HG}}^{(4)}\right)-\left(z_{i\text{AB}}^{(1)}\cdot
z_{i\text{BA}}^{(1)}+z_{i\text{CD}}^{(3)}\cdot z_{i\text{DC}}^{(3)}\right),$
(28)

and $y_{i}=1(p_{i}>0.5)$ for $i=1,\ldots,n+m$.

In first two simulations, GP with SE covariance kernel is among the best
performers (Table 1). The first simulation has a linear decision boundary in
terms of the 24 codes. On the other hand, the second simulation has a
quadratic decision boundary depending on the interactions between
$\\{AB,BA\\}$, $\\{CD,DC\\}$, $\\{EF,FE\\}$, and $\\{GH,HG\\}$. Logistic
regression and its regularized versions perform the best in the first
simulation, but their performance deteriorates slightly in the second
simulation due to the presence of interactions between dummy variables. Random
forest using the codes as covariates has the worst performance in both
simulations. After including SE kernel-based covariates, random forest
achieves the same level of performance as other method in both simulations.
The remaining methods are based on non-linear kernels, so they are more
flexible in adapting to both linear and non-linear decision boundaries.
Compared to the spectrum and boundrange kernels, the SE kernel in (5) is
better tuned for modeling the extra structure provided by the chronic
conditions. This implies that kernel-based methods perform well in both
simulations and that the GP with SE kernel performs better than SVM and KRR
with spectrum and boundrange kernels.

_Third simulation._ For this simulation, we used ICD codes of patients in the
UIHC EHR data with brain and other nervous system or breast as the cancer
primary sites. For the $i$th (pseudo) patient in this subset, let
$z_{i}=\|\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{T}}_{i})\|^{2}$
and $\delta_{i}=1$ if the patient with $\operatorname{\mathbf{T}}_{i}$ code
has brain and other nervous system cancer and $\delta_{i}=-1$ otherwise, where
$\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{T}}_{i})$ is defined in
(3). The response $y_{i}$ is simulated independently from Bernoulli ($p_{i}$),
where

$\displaystyle\log\frac{p_{i}}{1-p_{i}}=0.1+0.2x_{i}+3\tan(z_{i})+3\delta_{i},\quad
x_{i}\overset{\text{ind.}}{\sim}N(0,1)$ (29)

for $i=1,\ldots,n+m$. We simulate the data after setting $n=1000$ and $m=100$.

GP with SE covariance kernel is still among the top performers (Table 1). This
simulation has a nonlinear periodic decision boundary; therefore, the
performance of random forest and parametric models, including logistic
regression and its regularized, deteriorate significantly. While all the
kernel-based methods are suited for modeling the non-linear periodic decision
boundary in (29), the SE kernel is better tuned than spectrum and boundrange
kernels for modeling the hierarchical structure of ICD codes. This implies
that all kernel-based methods perform well, but the GP with SE kernel performs
better than SVM and KRR with spectrum and boundrange kernels.

_Summary._ Our simulation studies suggest that the kernel-based methods are
better suited for applications in practice where we expect interactions among
the diagnoses. Furthermore, GP with the SE kernel is easily extended to
account for any biomedical information, is tuned for modeling the structure of
ICD codes, and produces similar results as SVM or KRR with spectrum and
boundrange kernels in the absence of any additional structure. The
distinguishing feature of the proposed method is that the posterior MCMC draws
can be used for quantifying uncertainty in predictions and parameter
estimates, which is key in biomedical applications; therefore, we conclude
that the GP with SE covariance kernel is better than other kernel based
methods in quantifying uncertainty and in modeling the effect of chronic
conditions on the response.

Table 1: Performance comparisons for the three simulation studies. Every entry
in the table is the average of its values across ten replications. SE-GP is
the proposed method and SE kernel-based features are obtained from the kernel
matrix estimated using the proposed method. If a method fails to produce
results, then we indicate this using ‘-’. FIRST SIMULATION
---
| AUC | Accuracy | Sensitivity | Specificity
SE-GP | 0.96 | 0.96 | 0.96 | 0.96
| Logistic Regression with SE Kernel-Based Covariates
No Penalty | 0.95 | 0.95 | 0.96 | 0.95
Ridge Penalty | 0.95 | 0.96 | 0.96 | 0.95
Lasso Penalty | 0.95 | 0.95 | 0.96 | 0.95
| Logistic Regression
No Penalty | 0.96 | 0.96 | 0.96 | 0.96
Ridge Penalty | 0.96 | 0.96 | 0.96 | 0.96
Lasso Penalty | 0.96 | 0.96 | 0.96 | 0.96
| Kernel Ridge Regression
Spectrum | 0.95 | 0.96 | 0.96 | 0.95
Boundrange | 0.95 | 0.95 | 0.96 | 0.95
| Support Vector Machine
Spectrum | 0.95 | 0.96 | 0.96 | 0.96
Boundrange | 0.95 | 0.95 | 0.96 | 0.95
| Random Forest
Default | 0.50 | 0.50 | 0.99 | 0.01
SE Kernel-Based Covariates | 0.95 | 0.95 | 0.96 | 0.95
SECOND SIMULATION
| AUC | Accuracy | Sensitivity | Specificity
SE-GP | 1.00 | 0.99 | 1.00 | 0.98
| Logistic Regression with SE Kernel-Based Covariates
No Penalty | 0.90 | 0.82 | 0.77 | 0.86
Ridge Penalty | 0.90 | 0.82 | 0.81 | 0.84
Lasso Penalty | 0.90 | 0.82 | 0.75 | 0.89
| Logistic Regression
No Penalty | 0.90 | 0.82 | 0.78 | 0.86
Ridge Penalty | 0.90 | 0.82 | 0.80 | 0.84
Lasso Penalty | 0.90 | 0.83 | 0.77 | 0.87
| Kernel Ridge Regression
Spectrum | 0.92 | 0.86 | 0.86 | 0.86
Boundrange | 0.91 | 0.84 | 0.83 | 0.85
| Support Vector Machine
Spectrum | 0.92 | 0.86 | 0.86 | 0.86
Boundrange | 0.91 | 0.84 | 0.83 | 0.85
| Random Forest
Default | 0.50 | 0.45 | 1.00 | 0.00
SE Kernel-Based Covariates | 0.88 | 0.88 | 0.85 | 0.92
THIRD SIMULATION
| AUC | Accuracy | Sensitivity | Specificity
SE-GP | 0.94 | 0.91 | 0.90 | 0.92
| Logistic Regression with SE Kernel-Based Covariates
No Penalty | - | - | - | -
Ridge Penalty | 0.76 | 0.73 | 0.66 | 0.82
Lasso Penalty | 0.75 | 0.73 | 0.67 | 0.81
| Logistic Regression
No Penalty | - | - | - | -
Ridge Penalty | 0.89 | 0.84 | 0.85 | 0.84
Lasso Penalty | 0.90 | 0.86 | 0.85 | 0.88
| Kernel Ridge Regression
Spectrum | 0.84 | 0.80 | 0.82 | 0.78
Boundrange | 0.81 | 0.78 | 0.81 | 0.74
| Support Vector Machine
Spectrum | 0.85 | 0.81 | 0.82 | 0.81
Boundrange | 0.83 | 0.79 | 0.77 | 0.81
| Random Forest
Default | 0.85 | 0.86 | 0.95 | 0.76
SE Kernel-Based Covariates | 0.89 | 0.89 | 0.91 | 0.88

###  5.3 Real data analysis

The UIHC EHR data contains information about 1660 patients, including their
diagnoses, marital status, and the primary cancer site. There are six types of
cancer sites: brain and other nervous system (brain), breast, urinary system,
respiratory system, female genital system, and digestive system. Calderón-
Larrañaga et al. [2] define 58 chronic conditions using ICD codes that do not
include any diagnosis related to neoplasms. For the feature map in (6) that
defines the SE kernel, we set $C=58$,
$\operatorname{\mathbf{R}}=\operatorname{\mathbf{I}}$, and
$w_{c}\propto\kappa_{1}(\operatorname{\mathbf{T}}_{c},\operatorname{\mathbf{T}}_{c})$,
where $\operatorname{\mathbf{T}}_{c}$ is the set of ICD codes definining $c$th
chronic condition. We include marital status as the only demographic
predictor, so $p=2$ in (16). The major and minor goals of the analysis are to
estimate the marginal associations of the 58 chronic conditions with the
cancer sites and to predict the cancer primary site using the patient
diagnoses and marital status, respectively.

We use the methods from the previous section to achieve both goals. Five
classification models are used for predicting the six primary sites using
$y=1$ to denote the presence of cancer in brain, breast, urinary system,
respiratory system, and female genital system, respectively, in the five
models. For every model, all the methods are trained on 80% of the full data,
their performance is evaluated on the excluded data, and the setup is
replicated ten times. There is no principled way of using chronic conditions
as covariates in logistic regression and its penalized extension, so we only
use KRR, SVM, and SE-GP for estimating the marginal associations between the
58 chronic conditions and primary cancer site. To this end, we include
$\operatorname{\mathbf{T}}_{1},\ldots,\operatorname{\mathbf{T}}_{58}$ as the
diagnoses of 58 pseudo patients in the test data and summarize the marginal
associations using the estimate of Pr$(y=1\mid\operatorname{\mathbf{T}}_{c})$,
where $\operatorname{\mathbf{T}}_{c}$ is the diagnosis for the $c$th pseudo
patient ($c=1,\ldots,58$). Unlike SVM and KRR, SE-GP also provides 95%
credible intervals for the marginal associations using the MCMC draws of
Pr$(y=1)$.

SE-GP is among the top performers in all the five classification models
(Tables 2–6). We do not provide results for KRR and logistic regression with
no penalty because the former fails due to a line search error and the latter
cannot be used because the number of dummy variables is larger than the sample
size. Calderón-Larrañaga et al. [2] have developed a comprehensive list of
chronic conditions after a careful study, so we expect a subset of them to be
related to primary cancer sites. This is confirmed by a relatively high
accuracy of logistic regression with the ridge and lasso penalties for
predicting the cancer in brain and urinary system; however, these methods do
not use the information encoded in ICD codes, so their sensitivity and
specificity vary a lot; for example, logistic regression with the lasso
penalty has a high accuracy of 0.90 for predicting the cancer in urinary
system, but its specificity in the same model is 0.36.

On the other hand, SE-GP is among the top performers in terms of accuracy,
sensitivity, and specificity in all the five classification models. This also
implies that the AUC of SE-GP is the highest among all methods in all the five
classification models. The features obtained from the SE kernel matrix are
also promising in that if we use them as covariates in logistic regression and
its penalized extensions, then the AUCs are relatively large in all the five
classification models. The performance of default random forest is the worst
in all five classification models. AUC and accuracy increase after including
SE kernel-based covariates, but the results are still poorer than all other
methods.

SE-GP also performs better than SVM with spectrum and boundrange kernels
because SE kernel accounts for the additional structure provided by the 58
chronic conditions. The SE-GP is a restriction of the boundrange kernel that
is tuned for EHR data analysis, so its sensitivity and specificity are much
higher than those of the SVM in all the five classification models. Most
importantly, the estimates of marginal associations between chronic conditions
and primary cancer sites agree closely with those obtained using the SVM with
spectrum and boundrange kernels (Table 7). A key feature of the SE-GP approach
is that it provides 95% credible intervals in addition to the point estimates,
which are important for characterizing uncertainty in biomedical applications.
Additionally, these credible intervals cover the corresponding estimates
obtained using SVM. Based on our results in simulations and this section, we
conclude that SE-GP outperforms its competitors in estimating the marginal
associations among chronic conditions and primary cancer sites and in
predicting the cancer primary site using diagnoses and demographic information
as predictors.

Table 2: Performance in the classification model with $y=1$ denoting the presence of cancer in the brain and other nervous system. Every entry in the table is the average of its values across ten replications.  | AUC | Accuracy | Sensitivity | Specificity
---|---|---|---|---
SE-GP | 0.93 | 0.92 | 0.86 | 0.93
| Logistic Regression with SE Kernel-Based Covariates
No Penalty | 0.91 | 0.88 | 0.82 | 0.89
Lasso Penalty | 0.91 | 0.89 | 0.81 | 0.91
Ridge Penalty | 0.90 | 0.90 | 0.81 | 0.91
| Logistic Regression
Lasso Penalty | 0.90 | 0.92 | 0.81 | 0.94
Ridge Penalty | 0.93 | 0.90 | 0.89 | 0.90
| Support Vector Machine
Boundrange | 0.94 | 0.91 | 0.88 | 0.92
Spectrum | 0.94 | 0.90 | 0.90 | 0.90
| Random Forest
Default | 0.50 | 0.11 | 1.00 | 0.00
SE Kernel-Based Covariates | 0.72 | 0.92 | 0.45 | 0.98
Table 3: Performance in the classification model with $y=1$ denoting the presence of cancer in the breast. Every entry in the table is the average of its values across ten replications.  | AUC | Accuracy | Sensitivity | Specificity
---|---|---|---|---
SE-GP | 0.77 | 0.72 | 0.75 | 0.71
| Logistic Regression with SE Kernel-Based Covariates
No Penalty | 0.69 | 0.54 | 0.80 | 0.51
Lasso Penalty | 0.68 | 0.59 | 0.74 | 0.56
Ridge Penalty | 0.68 | 0.58 | 0.74 | 0.56
| Logistic Regression
Lasso Penalty | 0.74 | 0.66 | 0.74 | 0.66
Ridge Penalty | 0.73 | 0.72 | 0.69 | 0.73
| Support Vector Machine
Boundrange | 0.68 | 0.67 | 0.63 | 0.68
Spectrum | 0.74 | 0.67 | 0.75 | 0.66
| Random Forest
Default | 0.50 | 0.11 | 1.00 | 0.00
SE Kernel-Based Covariates | 0.53 | 0.80 | 0.18 | 0.88
Table 4: Performance in the classification model with $y=1$ denoting the presence of cancer in the urinary system. Every entry in the table is the average of its values across ten replications.  | AUC | Accuracy | Sensitivity | Specificity
---|---|---|---|---
SE-GP | 0.80 | 0.78 | 0.71 | 0.79
| Logistic Regression with SE Kernel-Based Covariates
No Penalty | 0.75 | 0.68 | 0.71 | 0.68
Lasso Penalty | 0.75 | 0.70 | 0.70 | 0.70
Ridge Penalty | 0.75 | 0.75 | 0.63 | 0.77
| Logistic Regression
Lasso Penalty | 0.69 | 0.89 | 0.36 | 0.97
Ridge Penalty | 0.74 | 0.76 | 0.68 | 0.77
| Support Vector Machine
Boundrange | 0.60 | 0.71 | 0.43 | 0.75
Spectrum | 0.68 | 0.76 | 0.53 | 0.79
| Random Forest
Default | 0.50 | 0.12 | 1.00 | 0.00
SE Kernel-Based Covariates | 0.58 | 0.89 | 0.17 | 0.98
Table 5: Performance in the classification model with $y=1$ denoting the presence of cancer in the respiratory system. Every entry in the table is the average of its values across ten replications.  | AUC | Accuracy | Sensitivity | Specificity
---|---|---|---|---
SE-GP | 0.83 | 0.75 | 0.81 | 0.73
| Logistic Regression with SE Kernel-Based Covariates
No Penalty | 0.67 | 0.56 | 0.78 | 0.52
Lasso Penalty | 0.67 | 0.57 | 0.76 | 0.53
Ridge Penalty | 0.67 | 0.56 | 0.78 | 0.52
| Logistic Regression
Lasso Penalty | 0.81 | 0.74 | 0.73 | 0.74
Ridge Penalty | 0.83 | 0.75 | 0.83 | 0.74
| Support Vector Machine
Boundrange | 0.77 | 0.70 | 0.73 | 0.70
Spectrum | 0.82 | 0.75 | 0.78 | 0.74
| Random Forest
Default | 0.50 | 0.16 | 1.00 | 0.00
SE Kernel-Based Covariates | 0.55 | 0.84 | 0.13 | 0.97
Table 6: Performance in the classification model with $y=1$ denoting the presence of cancer in female gential system. Every entry in the table is the average of its values across ten replications.  | AUC | Accuracy | Sensitivity | Specificity
---|---|---|---|---
SE-GP | 0.86 | 0.83 | 0.72 | 0.86
| Logistic Regression with SE Kernel-Based Covariates
No Penalty | 0.78 | 0.73 | 0.66 | 0.75
Lasso Penalty | 0.78 | 0.75 | 0.65 | 0.78
Ridge Penalty | 0.78 | 0.75 | 0.66 | 0.78
| Logistic Regression
Lasso Penalty | 0.84 | 0.77 | 0.75 | 0.77
Ridge Penalty | 0.85 | 0.81 | 0.77 | 0.82
| Support Vector Machine
Boundrange | 0.77 | 0.75 | 0.68 | 0.77
Spectrum | 0.83 | 0.80 | 0.70 | 0.83
| Random Forest
Default | 0.50 | 0.21 | 1.00 | 0.00
SE Kernel-Based Covariates | 0.70 | 0.84 | 0.45 | 0.95
Table 7: Summary of the marginal associations between 58 chronic conditions and the six primary cancer sites. The marginal association between a chronic condition and a primary cancer site denoted using $y=1$ is defined as the point estimate of Pr$(y=1\mid\operatorname{\mathbf{T}}_{c})$, where $\operatorname{\mathbf{T}}_{c}$ is the diagnosis of the pseudo patient representing the chronic condition $c$. For SE-GP, we define the point estimate of Pr$(y=1\mid x^{*})$ to be its posterior mean and use the 95% credible interval to provide an estimate of uncertainty. For the 58 pseudo patients and six primary cancer sites, we obtain the point estimate and 95% credible interval using the MCMC draws. The table includes only those chronic conditions for which the point estimates of the marginal associations are greater than 0.5 and a ‘-’ indicates the absence of such chronic conditions. Kernel | Chronic Condition | Estimate
---|---|---
Brain and Other Nervous System
SE-GP | OTHER NEUROLOGICAL DISEASES | 0.8075, (0.5199, 0.9649)
Boundrange | OTHER NEUROLOGICAL DISEASES | 0.8057
Spectrum | OTHER NEUROLOGICAL DISEASES | 0.6570
Breast
SE-GP | DEPRESSION AND MOOD DISEASES | 0.6909, (0.2960, 0.9451)
Boundrange | - | -
Spectrum | OTHER CARDIOVASCULAR DISEASES | 0.5799
Spectrum | OSTEOPOROSIS | 0.5884
Spectrum | ALLERGY | 0.6005
Spectrum | DEPRESSION AND MOOD DISEASES | 0.6454
Spectrum | VENOUS AND LYMPHATIC DISEASES | 0.7290
Urinary System
SE-GP | - | -
Boundrange | - | -
Spectrum | - | -
Respiratory System
SE-GP | OTHER RESPIRATORY DISEASES | 0.7973, (0.5421, 0.9483)
SE-GP | COPD, EMPHYSEMA, CHRONIC BRONCHITIS | 0.8815, (0.6181, 0.9880)
Boundrange | OTHER RESPIRATORY DISEASES | 0.6901
Boundrange | COPD, EMPHYSEMA, CHRONIC BRONCHITIS | 0.7575
Spectrum | OTHER RESPIRATORY DISEASES | 0.5024
Spectrum | OTHER PSYCHIATRIC AND BEHAVIORAL DISEASES | 0.5871
Spectrum | COPD, EMPHYSEMA, CHRONIC BRONCHITIS | 0.8992
Female Genital System
SE-GP | OBESITY | 0.6616, (0.3507, 0.8915)
Boundrange | OBESITY | 0.5276
Spectrum | - | -
Digestive System
SE-GP | COLITIS AND RELATED DISEASES | 0.5114, (0.1576, 0.8618)
SE-GP | OTHER METABOLIC DISEASES | 0.5382, (0.2198, 0.8286)
SE-GP | INFLAMMATORY BOWEL DISEASES | 0.5874, (0.2478, 0.8788)
SE-GP | CHRONIC PANCREAS, BILIARY TRACT AND GALLBLADDER DISEASES | 0.7002, (0.3003, 0.9464)
SE-GP | CHRONIC LIVER DISEASES | 0.7893, (0.4196, 0.9686)
SE-GP | ESOPHAGUS, STOMACH AND DUODENUM DISEASES | 0.8413, (0.5631, 0.9754)
Boundrange | OTHER METABOLIC DISEASES | 0.5736
Boundrange | INFLAMMATORY BOWEL DISEASES | 0.8326
Boundrange | OTHER DIGESTIVE DISEASES | 0.9170
Boundrange | CHRONIC PANCREAS, BILIARY TRACT AND GALLBLADDER DISEASES | 0.9281
Boundrange | COLITIS AND RELATED DISEASES | 0.9389
Boundrange | ESOPHAGUS, STOMACH AND DUODENUM DISEASES | 0.9949
Boundrange | CHRONIC LIVER DISEASES | 0.9958
Spectrum | OTHER METABOLIC DISEASES | 0.5052
Spectrum | OTHER DIGESTIVE DISEASES | 0.6525
Spectrum | COLITIS AND RELATED DISEASES | 0.6580
Spectrum | CHRONIC PANCREAS, BILIARY TRACT AND GALLBLADDER DISEASES | 0.7498
Spectrum | CHRONIC LIVER DISEASES | 0.8930
Spectrum | ESOPHAGUS, STOMACH AND DUODENUM DISEASES | 0.9323

##  6 Discussion

Our approach can be extended in several ways. First, the sample size in the
motivating application is small so that repeated computations and evaluation
of the kernel function is relatively inexpensive. This, however, becomes
problematic as the sample size becomes moderately large. The kernels developed
in this work can be immediately extended for nonparametric regression and
classification in large sample settings using low-rank kernels based on
inducing points [14]. Second, Every element of the kernel matrix is a
similarity measure between a pair of patients. This can be used an input in an
algorithm for unsupervised learning using ICD codes as inputs, such as
clustering of subsets, data visualization, and dimension reduction [18,
Chapters 6, 8]. The PG-DA strategy outlined in Section 4.3 is easily extended
to model zero-inflated negative binomial responses following Neelon [10].
Finally, We are exploring application of product partition model with
regression on diagnoses, where the similarity measure on diagnoses is defined
using the kernel matrix.

## Acknowledgement

We thank UIHC for sharing the EHR data. Sanvesh Srivastava’s research is
supported by grants from the Office of Naval Research (ONR-BAA N000141812741)
and the National Science Foundation (DMS-1854667/1854662). Stephanie
Gilbertson-White’s research is supported by the University of Iowa, College of
Nursing Center on Advancing Multimorbidity Science (CAMS) funded be the
National Institute for Nursing Research (P20NR018081). The R code for
simulations is available at https://github.com/blayes/ehr-data-analysis.

## References

  * Banerjee et al. [2014] Banerjee, S., B. P. Carlin, and A. E. Gelfand (2014).  Hierarchical Modeling and Analysis for Spatial Data.  CRC Press, Boca Raton, FL.
  * Calderón-Larrañaga et al. [2017] Calderón-Larrañaga, A., D. L. Vetrano, G. Onder, L. A. Gimeno-Feliu, C. Coscollar-Santaliestra, A. Carfí, M. S. Pisciotta, S. Angleman, R. J. Melis, and G. Santoni (2017).  Assessing and measuring chronic multimorbidity in the older population: a proposal for its operationalization.  Journals of Gerontology Series A: Biomedical Sciences and Medical Sciences 72(10), 1417–1423.
  * Friedman et al. [2009] Friedman, J., T. Hastie, and R. Tibshirani (2009).  glmnet: Lasso and elastic-net regularized generalized linear models.  R package version 1(4).
  * Ghosal and van der Vaart [2017] Ghosal, S. and A. van der Vaart (2017).  Fundamentals of nonparametric Bayesian inference, Volume 44.  Cambridge University Press.
  * Henao et al. [2016] Henao, R., J. T. Lu, J. E. Lucas, J. Ferranti, and L. Carin (2016).  Electronic health record analysis via deep poisson factor models.  The Journal of Machine Learning Research 17(1), 6422–6453.
  * Hong et al. [2019] Hong, C., K. P. Liao, and T. Cai (2019).  Semi-supervised validation of multiple surrogate outcomes with application to electronic medical records phenotyping.  Biometrics 75(1), 78–89.
  * Horng et al. [2017] Horng, S., D. A. Sontag, Y. Halpern, Y. Jernite, N. I. Shapiro, and L. A. Nathanson (2017).  Creating an automated trigger for sepsis clinical decision support at emergency department triage using machine learning.  PloS one 12(4), e0174708.
  * Karatzoglou et al. [2004] Karatzoglou, A., A. Smola, K. Hornik, and A. Zeileis (2004).  kernlab – an S4 package for kernel methods in R.  Journal of Statistical Software 11(9), 1–20.
  * Lodhi et al. [2002] Lodhi, H., C. Saunders, J. Shawe-Taylor, N. Cristianini, and C. Watkins (2002).  Text classification using string kernels.  Journal of Machine Learning Research 2(Feb), 419–444.
  * Neelon [2019] Neelon, B. (2019).  Bayesian zero-inflated negative binomial regression based on pólya-gamma mixtures.  Bayesian Analysis 14(3), 829–855.
  * Ni et al. [2019] Ni, Y., P. Müller, and Y. Ji (2019).  Bayesian double feature allocation for phenotyping with electronic health records.  Journal of the American Statistical Association, 1–15.
  * Nishihara et al. [2014] Nishihara, R., I. Murray, and R. P. Adams (2014).  Parallel MCMC with generalized elliptical slice sampling.  The Journal of Machine Learning Research 15(1), 2087–2112.
  * Polson et al. [2013] Polson, N. G., J. G. Scott, and J. Windle (2013).  Bayesian inference for logistic models using pólya–gamma latent variables.  Journal of the American statistical Association 108(504), 1339–1349.
  * Quiñonero-Candela and Rasmussen [2005] Quiñonero-Candela, J. and C. E. Rasmussen (2005).  A unifying view of sparse approximate gaussian process regression.  Journal of Machine Learning Research 6(Dec), 1939–1959.
  * Rajkomar et al. [2018] Rajkomar, A., E. Oren, K. Chen, A. M. Dai, N. Hajaj, M. Hardt, P. J. Liu, X. Liu, J. Marcus, M. Sun, et al. (2018).  Scalable and accurate deep learning with electronic health records.  NPJ Digital Medicine 1(1), 18.
  * Rasmussen and Williams [2006] Rasmussen, C. E. and C. K. Williams (2006).  Gaussian processes for machine learning.  MIT Press.
  * Robin et al. [2011] Robin, X., N. Turck, A. Hainard, N. Tiberti, F. Lisacek, J.-C. Sanchez, and M. Müller (2011).  proc: an open-source package to analyze and compare roc curves.  BMC Bioinformatics 12, 77.
  * Shawe-Taylor and Cristianini [2004] Shawe-Taylor, J. and N. Cristianini (2004).  Kernel methods for pattern analysis.  Cambridge University Press.
  * Wang and Roy [2018] Wang, X. and V. Roy (2018).  Analysis of the pólya-gamma block gibbs sampler for bayesian logistic linear mixed models.  Statistics & Probability Letters 137, 251–256.
  * Wright and Ziegler [2017] Wright, M. N. and A. Ziegler (2017).  ranger: A fast implementation of random forests for high dimensional data in C++ and R.  Journal of Statistical Software 77(1), 1–17.
  * Xu et al. [2019] Xu, K., M. Lam, J. Pang, X. Gao, C. Band, P. Mathur, F. Papay, A. K. Khanna, J. B. Cywinski, K. Maheshwari, et al. (2019).  Multimodal machine learning for automated icd coding.  In Machine Learning for Healthcare Conference, pp. 197–215. PMLR.

## Supplemental Document for Gaussian Process Regression and Classification
using International Classification of Disease Codes as Covariates

##  A Proof of Theorem 1

We recall the definitions of kernel functions on ICD codes. Let
$\operatorname{{\mathbf{\psi}}}(t)$ be the feature map of an ICD code
$t\in\operatorname{\mathcal{T}}$ and
$\operatorname{\mathcal{F}}=\\{\operatorname{{\mathbf{\psi}}}(t):t\in\operatorname{\mathcal{T}}\\}$
be the feature space. For every $t\in\operatorname{\mathcal{T}}$,
$\operatorname{{\mathbf{\psi}}}(t)\in\\{0,1\\}^{|t^{*}||\operatorname{\mathcal{B}}|}$
and is defined as

$\displaystyle\operatorname{{\mathbf{\psi}}}:\operatorname{\mathcal{T}}\mapsto\operatorname{\mathcal{F}},\quad\operatorname{{\mathbf{\psi}}}(t)=\\{\operatorname{{\mathbf{\psi}}}_{b}(t)\\}_{b\in\operatorname{\mathcal{B}}},\quad\operatorname{{\mathbf{\psi}}}_{b}(t)=\\{1(b=t_{1:1}),\ldots,1(b=t_{1:|t^{*}|})\\}^{\mathrm{\scriptscriptstyle
T}},$ (30)

where $1(\cdot)$ is 1 if $\cdot$ is true and 0 otherwise. The kernel function
on $\operatorname{\mathcal{T}}\times\operatorname{\mathcal{T}}$ is defined
using $\operatorname{{\mathbf{\psi}}}(\cdot)$ and known $|t^{*}|\times|t^{*}|$
symmetric positive semi-definite matrices
$(\operatorname{{\mathbf{\Lambda}}}_{b})_{b\in\operatorname{\mathcal{B}}}$ as

$\displaystyle\kappa_{0}(t,t^{\prime}\mid\operatorname{{\mathbf{\Lambda}}})=\underset{b\in\operatorname{\mathcal{B}}}{\sum}\operatorname{{\mathbf{\psi}}}^{\mathrm{\scriptscriptstyle
T}}_{b}(t)\operatorname{{\mathbf{\Lambda}}}_{b}\operatorname{{\mathbf{\psi}}}_{b}(t^{\prime})\equiv\operatorname{{\mathbf{\psi}}}^{\mathrm{\scriptscriptstyle
T}}(t)\operatorname{{\mathbf{\Lambda}}}\operatorname{{\mathbf{\psi}}}(t^{\prime}),\;t,t^{\prime}\in\operatorname{\mathcal{T}},\quad\operatorname{{\mathbf{\Lambda}}}={\mathrm{diag}}\\{\operatorname{{\mathbf{\Lambda}}}_{b}:b\in\operatorname{\mathcal{B}}\\},$
(31)

where $\operatorname{{\mathbf{\Lambda}}}$ is a block diagonal matrix.

We use $\kappa_{0}$ in (31) to define the kernel function on subsets of ICD
codes that represent diagnoses. Let $2^{\operatorname{\mathcal{T}}}$ be the
power set of $\operatorname{\mathcal{T}}$,
$\operatorname{\mathbf{t}}=\\{t_{1},\ldots,t_{r}\\}$ be an element of
$2^{\operatorname{\mathcal{T}}}$, $\Psi(\operatorname{\mathbf{t}})$ be the
feature map of $\operatorname{\mathbf{t}}$, and
$\operatorname{\mathcal{F}}_{\Psi}=\\{\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}):\operatorname{\mathbf{t}}\in
2^{\operatorname{\mathcal{T}}}\\}$ be the feature space of subsets of ICD
codes. The feature vector
$\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})$ and the kernel
function $\kappa_{1}$ on $2^{\operatorname{\mathcal{T}}}\times
2^{\operatorname{\mathcal{T}}}$ are defined using
$\operatorname{{\mathbf{\psi}}}(\cdot)$ in (30) and
$\kappa_{0}(\cdot,\cdot\mid\operatorname{{\mathbf{\Lambda}}})$ in (31),
respectively, as

$\displaystyle\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})=\sum_{{t\in\operatorname{\mathbf{t}}}}\operatorname{{\mathbf{\psi}}}(t),\quad\kappa_{1}(\operatorname{\mathbf{t}},\operatorname{\mathbf{t}}^{\prime}\mid\operatorname{{\mathbf{\Lambda}}})=\sum_{{t\in\operatorname{\mathbf{t}}}}\sum_{{t^{\prime}\in\operatorname{\mathbf{t}}^{\prime}}}\kappa_{0}(t,t^{\prime}\mid\operatorname{{\mathbf{\Lambda}}})=\operatorname{{\mathbf{\Psi}}}^{\mathrm{\scriptscriptstyle
T}}(\operatorname{\mathbf{t}})\operatorname{{\mathbf{\Lambda}}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}^{\prime}),\quad\operatorname{\mathbf{t}},\operatorname{\mathbf{t}}^{\prime}\in
2^{\operatorname{\mathcal{T}}},$ (32)

where $\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})$ is a
$(|t^{*}||\operatorname{\mathcal{B}}|)$-dimensional vector with entries in
$\\{0\\}\cup\operatorname{\mathbb{N}}$ and $\operatorname{{\mathbf{\Lambda}}}$
is defined as in (31).

We account the additional structure provided by the $C$ chronic conditions. A
patient’s diagnoses are represented using the set
$\operatorname{\mathbf{T}}=\\{\operatorname{\mathbf{t}}_{1},\ldots,\operatorname{\mathbf{t}}_{c},\ldots,\operatorname{\mathbf{t}}_{C}\\}\in\left(2^{\operatorname{\mathcal{T}}}\right)^{C}$,
where $\operatorname{\mathbf{t}}_{c}\subseteq\operatorname{\mathbf{T}}_{c}$. A
kernel function on
$\left(2^{\operatorname{\mathcal{T}}}\right)^{C}\times\left(2^{\operatorname{\mathcal{T}}}\right)^{C}$
that accounts for the structure imposed by $C$ chronic conditions is

$\displaystyle\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})=\sum_{c=1}^{C}w_{c}\,\kappa_{1}(\operatorname{\mathbf{t}}_{c},\operatorname{\mathbf{t}}^{\prime}_{c}\mid\operatorname{{\mathbf{\Lambda}}}),\quad\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\in\left(2^{\operatorname{\mathcal{T}}}\right)^{C},\quad
w_{1},\ldots,w_{C}>0,$ (33)

where $\operatorname{{\mathbf{\Lambda}}}$ is defined in (31),
$\operatorname{\mathbf{w}}=(w_{1},\ldots,w_{C})$, and $w_{c}$ indicates the
importance of chronic condition $c$ in defining the similarity between two
patients with diagnoses $\operatorname{\mathbf{T}}$ and
$\operatorname{\mathbf{T}}^{\prime}$, respectively. Finally, we normalize
$\kappa_{2}$ to obtain our kernel function

$\displaystyle\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\frac{\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})}{\left\\{\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{1/2}\left\\{\kappa_{2}(\operatorname{\mathbf{T}}^{\prime},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{1/2}},\quad\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\in\left(2^{\operatorname{\mathcal{T}}}\right)^{C}.$
(34)

We restate Theorem 1.

######  Theorem A.1.

Let
$\\{\operatorname{{\mathbf{\Lambda}}}_{b}\\}_{b\in\operatorname{\mathcal{B}}}$
and $\operatorname{\mathbf{w}}$ be defined as in (31) and (33),
$\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\in(2^{\operatorname{\mathcal{T}}})^{C}$,
$\sigma^{2}>0$, $\phi>0$, and

$\displaystyle
d(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\left\\{\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})+\kappa(\operatorname{\mathbf{T}}^{\prime},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})-2\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{1/2}.$
(35)

Then,

  1. 1.

given $\operatorname{{\mathbf{\Lambda}}}$ and $\operatorname{\mathbf{w}}$,
$\kappa(\cdot,\cdot\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
in (34) and
$d(\cdot,\cdot\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
in (35) are valid kernel and distance functions on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$;

  2. 2.

if
$\operatorname{{\mathbf{\Lambda}}}=\operatorname{\mathbf{R}}^{\top}\operatorname{\mathbf{R}}$,
where $\operatorname{\mathbf{R}}$ is an upper triangular matrix, then the
feature map of $\kappa$ is

$\displaystyle\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})=\frac{\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})}{\|\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\|},\quad\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})=\\{w_{1}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}_{1}),\ldots,w_{C}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}_{C})\\}^{\mathrm{\scriptscriptstyle
T}},$ (36)

where
$\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
is a $(C|t^{*}||\operatorname{\mathcal{B}}|)$-dimensional vector; and

  3. 3.

the equivalents of polynomial kernel of order $s$ and $\gamma$-exponential
kernel on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$,
respectively, are defined as

$\displaystyle\kappa_{s}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\sigma^{2},\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\sigma^{2}\left\\{1+\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{s},\quad
s\in\operatorname{\mathbb{N}},$
$\displaystyle\kappa_{\gamma}^{e}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\sigma^{2},\phi,\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\sigma^{2}e^{-\phi\,\\{d(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\\}^{\gamma}},\quad\gamma\in[1,2],$
(37)

where $\sigma^{2}$ and $\phi$ play the role of variance and inverse length-
scale parameters, respectively, and $\kappa_{\gamma}^{e}$ reduces to the
exponential and SE kernels when $\gamma$ equals 1 and 2, respectively.

###### Proof.

Proof of (1): Because $\operatorname{{\mathbf{\Lambda}}}_{b}$ is a symmetric
positive semi-definite matrix, then Cholesky decomposition implies that there
exists an upper triangular matrix with a non-negative diagonal elements
$\operatorname{\mathbf{R}}_{b}$ such that
$\operatorname{\mathbf{R}}_{b}^{\top}\operatorname{\mathbf{R}}_{b}=\operatorname{{\mathbf{\Lambda}}}_{b}$.
Let
$\operatorname{\mathbf{R}}={\mathrm{diag}}\\{\operatorname{\mathbf{R}}_{b}:b\in\operatorname{\mathcal{B}}\\}$
and
$\operatorname{{\mathbf{\Lambda}}}={\mathrm{diag}}\\{\operatorname{{\mathbf{\Lambda}}}_{b}:b\in\operatorname{\mathcal{B}}\\}$.
Then,
$\operatorname{\mathbf{R}}^{\top}\operatorname{\mathbf{R}}=\operatorname{{\mathbf{\Lambda}}}$
and (31) implies that

$\displaystyle\kappa_{0}(t,t^{\prime}\mid\operatorname{{\mathbf{\Lambda}}})=\operatorname{{\mathbf{\psi}}}(t)^{\top}\operatorname{\mathbf{R}}^{\top}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\psi}}}(t^{\prime})=\left\\{\operatorname{\mathbf{R}}\operatorname{{\mathbf{\psi}}}(t)\right\\}^{\top}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\psi}}}(t^{\prime})=\left\\{\operatorname{\mathbf{R}}\operatorname{{\mathbf{\psi}}}(t^{\prime})\right\\}^{\top}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\psi}}}(t)=\kappa_{0}(t^{\prime},t\mid\operatorname{{\mathbf{\Lambda}}}).$
(38)

For any $t_{1},\ldots,t_{n}\in\operatorname{\mathcal{T}}$, define
$\operatorname{\mathbf{G}}$ to be the kernel matrix with
$\operatorname{\mathbf{G}}_{ij}=\kappa_{0}(t_{i},t_{j})$
($i,j\in\\{1,\ldots,n\\}$). Then, $\operatorname{\mathbf{G}}$ is symmetric
because
$\operatorname{\mathbf{G}}_{ij}=\kappa_{0}(t_{i},t_{j})=\kappa_{0}(t_{j},t_{i})=\operatorname{\mathbf{G}}_{ji}$
using (38) and, for any
$\operatorname{\mathbf{v}}\in\operatorname{\mathbb{R}}^{n}$,

$\displaystyle\operatorname{\mathbf{v}}^{\top}\operatorname{\mathbf{G}}\operatorname{\mathbf{v}}$
$\displaystyle=\sum_{i,j=1}^{n}v_{i}\operatorname{\mathbf{G}}_{ij}v_{j}=\sum_{i=1}^{n}\sum_{j=1}^{n}v_{i}\kappa_{0}(t_{i},t_{j})v_{j}\overset{(i)}{=}\sum_{i=1}^{n}\sum_{j=1}^{n}v_{i}\left\\{\operatorname{\mathbf{R}}\operatorname{{\mathbf{\psi}}}(t_{i})\right\\}^{\top}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\psi}}}(t_{j})v_{j}$
$\displaystyle=\left\\{\operatorname{\mathbf{R}}\sum_{i=1}^{n}v_{i}\operatorname{{\mathbf{\psi}}}(t_{i})\right\\}^{\top}\operatorname{\mathbf{R}}\sum_{j=1}^{n}v_{j}\operatorname{{\mathbf{\psi}}}(t_{j})=\|\operatorname{\mathbf{R}}\sum_{i=1}^{n}v_{i}\operatorname{{\mathbf{\psi}}}(t_{i})\|^{2}\geq
0,$ (39)

where $(i)$ follows from (38); therefore,
$\kappa_{0}(\cdot,\cdot\mid\operatorname{{\mathbf{\Lambda}}})$ is a kernel
function with feature map
$\operatorname{\mathbf{R}}\operatorname{{\mathbf{\psi}}}(t)$ for any
$t\in\operatorname{\mathcal{T}}$. Similarly, (32) implies that for any
$\operatorname{\mathbf{t}},\operatorname{\mathbf{t}}^{\prime}\in
2^{\operatorname{\mathcal{T}}}$,

$\displaystyle\kappa_{1}(\operatorname{\mathbf{t}},\operatorname{\mathbf{t}}^{\prime}\mid\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\left\\{\operatorname{\mathbf{R}}\sum_{t\in\operatorname{\mathbf{t}}}\operatorname{{\mathbf{\psi}}}(t)\right\\}^{\top}\operatorname{\mathbf{R}}\sum_{t^{\prime}\in\operatorname{\mathbf{t}}^{\prime}}\operatorname{{\mathbf{\psi}}}(t^{\prime})=\left\\{\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})\right\\}^{\top}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}^{\prime})$
$\displaystyle=\left\\{\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}^{\prime})\right\\}^{\top}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})=\left\\{\operatorname{\mathbf{R}}\sum_{t^{\prime}\in\operatorname{\mathbf{t}}^{\prime}}\operatorname{{\mathbf{\psi}}}(t^{\prime})\right\\}^{\top}\operatorname{\mathbf{R}}\sum_{t\in\operatorname{\mathbf{t}}}\operatorname{{\mathbf{\psi}}}(t).$
(40)

Following the same arguments used in (A), we have that
$\kappa_{1}(\cdot,\cdot\mid\operatorname{{\mathbf{\Lambda}}})$ is a kernel
function with feature map
$\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}})$
for any $\operatorname{\mathbf{t}}\in 2^{\operatorname{\mathcal{T}}}$.

The closure properties of kernel functions implies that $\kappa_{2}$ in (33)
is a valid kernel function. Proposition 3.22 on Page 75 in Shawe-Taylor and
Cristianini [18] implies that
$w_{c}\kappa_{1}(\cdot,\cdot\mid\operatorname{{\mathbf{\Lambda}}})$ is a
kernel function because
$\kappa_{1}(\cdot,\cdot\mid\operatorname{{\mathbf{\Lambda}}})$ is a kernel
function and $w_{c}>0$ ($c=1,\ldots,C$). The same proposition also implies
that
$\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})=\sum_{c=1}^{C}w_{c}\kappa_{1}(\operatorname{\mathbf{t}}_{c},\operatorname{\mathbf{t}}_{c}\mid\operatorname{{\mathbf{\Lambda}}})$
is a valid kernel function for any
$\operatorname{\mathbf{T}}=\\{\operatorname{\mathbf{t}}_{1},\ldots,\operatorname{\mathbf{t}}_{C}\\},\operatorname{\mathbf{T}}^{\prime}=\\{\operatorname{\mathbf{t}}_{1}^{\prime},\ldots,\operatorname{\mathbf{t}}_{C}^{\prime}\\}\in(2^{\operatorname{\mathcal{T}}})^{C}$.
The feature map of
$\kappa_{2}(\cdot,\cdot\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
is
$\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
for any $\operatorname{\mathbf{T}}\in(2^{\operatorname{\mathcal{T}}})^{C}$,
and it is obtained following the same steps in (A) and (A) as

$\displaystyle\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\sum_{c=1}^{C}w_{c}\kappa_{1}(\operatorname{\mathbf{t}}_{c},\operatorname{\mathbf{t}}^{\prime}_{c}\mid\operatorname{{\mathbf{\Lambda}}})=\sum_{c=1}^{C}\left\\{w_{c}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}_{c})\right\\}^{\top}w_{c}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}^{\prime}_{c})$
$\displaystyle=\\{w_{1}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}_{1}),\ldots,w_{C}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}_{C})\\}^{\top}\\{w_{1}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}_{1}^{\prime}),\ldots,w_{C}^{1/2}\operatorname{\mathbf{R}}\operatorname{{\mathbf{\Psi}}}(\operatorname{\mathbf{t}}_{C}^{\prime})\\}$
$\displaystyle\equiv\\{\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\\}^{\top}\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime}).$
(41)

The dimension of
$\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
is $C$ times the dimension of $\operatorname{{\mathbf{\psi}}}(t)$ for any
$t\in\operatorname{\mathcal{T}}$, which equals
$C|t^{*}||\operatorname{\mathcal{B}}|$.

Finally, the kernel $\kappa$ in (34) is obtained by the normalization of
$\kappa_{2}$. First, normalize the feature map
$\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
for any $\operatorname{\mathbf{T}}\in(2^{\operatorname{\mathcal{T}}})^{C}$ by
its Euclidean norm as

$\displaystyle\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
$\displaystyle=\frac{\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})}{\|\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\|}\overset{(i)}{\in}[0,1]^{C|t^{*}||\operatorname{\mathcal{B}}|},$
$\displaystyle\|\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\|^{2}$
$\displaystyle=\\{\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\\}^{\top}\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})=\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}}),$
(42)

where $(i)$ follows from the Cauchy-Schwartz inequality. Second, noting that
$\|\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\|=\\{\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\\}^{1/2}$,
define the normalized $\kappa_{2}$ kernel with feature map
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
as

$\displaystyle\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\\}^{\top}\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})=\frac{\\{\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\\}^{\top}\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})}{\|\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\|\|\operatorname{{\mathbf{\Psi}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})\|}$
$\displaystyle=\frac{\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})}{\\{\kappa_{2}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\\}^{1/2}\\{\kappa_{2}(\operatorname{\mathbf{T}}^{\prime},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\\}^{1/2}}.$
(43)

Using the feature map
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
in the definition of distance function in (35), we have that

$\displaystyle\\{d(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\\}^{2}$
$\displaystyle=\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})+\kappa(\operatorname{\mathbf{T}}^{\prime},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})-2\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
$\displaystyle=\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\\}^{\top}\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})+\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})\\}^{\top}\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})-2\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\\}^{\top}\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})$
$\displaystyle=\|\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\|^{2}+\|\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})\|^{2}-2\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})\\}^{\top}\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})$
$\displaystyle=\|\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})-\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})\|^{2};$
(44)

therefore,
$d(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
is a valid distance function because it equals
$\|\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})-\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})\|$,
the Euclidean distance between
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
and
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}^{\prime})$.

Proof of (2): Using (A), we have that
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
is the feature map of $\kappa$.

Proof of (3): Proposition 3.24 on Page 76 in Shawe-Taylor and Cristianini [18]
implies that if $\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}})$
is a kernel over
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$,
then
$\left\\{1+\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{s}$,
where $s\in\operatorname{\mathbb{N}}$, is also a valid kernel over
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$
because it is a degree-$s$ polynomial in
$\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
with positive coefficients. Furthermore, scaling by $\sigma^{2}>0$ implies
that
$\kappa_{s}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\sigma^{2},\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})=\sigma^{2}\left\\{1+\kappa(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})\right\\}^{s}$
in (A.1) is also a valid kernel over
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$
for any $s\in\operatorname{\mathbb{N}}$.

Finally,
$\kappa_{\gamma}^{e}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\sigma^{2},\phi,\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
in (A.1) is the $\gamma$-exponential kernel defined over
$[0,1]^{(C|t^{*}||\operatorname{\mathcal{B}}|)}\times[0,1]^{(C|t^{*}||\operatorname{\mathcal{B}}|)}$
because
$d(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})=\|\overline{\operatorname{{\mathbf{\Psi}}}}(\operatorname{\mathbf{T}})-\overline{\operatorname{{\mathbf{\Psi}}}}(\operatorname{\mathbf{T}})\|$
and
$\overline{\operatorname{{\mathbf{\Psi}}}}(\operatorname{\mathbf{T}})\in[0,1]^{(C|t^{*}||\operatorname{\mathcal{B}}|)}$
using its definition in (A). Section 4.2.1 in Rasmussen and Williams [16]
shows that this is a valid kernel function, so
$\kappa_{\gamma}^{e}(\operatorname{\mathbf{T}},\operatorname{\mathbf{T}}^{\prime}\mid\sigma^{2},\phi,\operatorname{\mathbf{w}},\operatorname{{\mathbf{\Lambda}}})$
is a valid kernel function on
$(2^{\operatorname{\mathcal{T}}})^{C}\times(2^{\operatorname{\mathcal{T}}})^{C}$;
see Equation 4.18 on Page 86 in Rasmussen and Williams [16]. ∎

##  B Proof of Theorem 2

Recall that the population regression and classification models are

$\displaystyle y=\mu_{0}(\Psi)+\epsilon,\;\epsilon\sim
N(0,\tau^{2}_{0}),\;\tau^{2}_{0}>0,\quad\log\,\frac{\Pr(y=1)}{\Pr(y=0)}=\mu_{0}(\Psi),\quad\Psi\in[0,1]^{d},$
(45)

where $\tau_{0}^{2}$ is the error variance,
$\operatorname{\mathbf{T}}\in(2^{\operatorname{\mathcal{T}}})^{C}$ is embedded
in $[0,1]^{d}$ as
$\Psi=\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$,
$d=C|t^{*}||\operatorname{\mathcal{B}}|$ using the definition of
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$
in (A), and the effect of $\operatorname{\mathbf{T}}$ on $y$ is modeled
through the feature map
$\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}})$,
which is fixed for a given
$\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}$. The distance between
$f(\cdot)$ and $\mu_{0}(\cdot)$ is defined as the empirical norm

$\displaystyle\|f-\mu_{0}\|_{n}$
$\displaystyle=\left(\frac{1}{n}\sum_{i=1}^{n}\left[f\left\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}_{i})\right\\}-\mu_{0}\left\\{\overline{\operatorname{{\mathbf{\Psi}}}}_{\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}}(\operatorname{\mathbf{T}}_{i})\right\\}\right]^{2}\right)^{1/2}$
$\displaystyle\equiv\left[\frac{1}{n}\sum_{i=1}^{n}\\{f(\Psi_{i})-\mu_{0}(\Psi_{i})\\}^{2}\right]^{1/2},$
(46)

where $\|\cdot\|_{n}$ conditions on the covariates
$\operatorname{\mathbf{T}}_{1},\ldots,\operatorname{\mathbf{T}}_{n}$ and its
dependence on $\operatorname{\mathbf{w}}$ and $\operatorname{\mathbf{R}}$ is
suppressed. We now restate Theorem 2.

######  Theorem B.1.

Assume that the regression and classification models in (45) are true, the
Matérn covariance kernel has smoothness $\nu_{1}$, and $\sigma$ and $\phi$
parameters of the covariance function are known.

  1. 1.

If $\mu_{0}$ is $\nu_{2}$-regular in the sense that
$\mu_{0}\in\mathfrak{C}^{\nu_{2}}[0,1]^{d}\cap\mathfrak{H}^{\nu_{2}}[0,1]^{d}$
with $\nu_{2}\leq\nu_{1}$, then as $n\rightarrow\infty$,

$\displaystyle\Pi_{\text{Mat},n}^{\text{Reg}}\left\\{f:\|f-\mu_{0}\|_{n}>C_{1}\epsilon_{n}\mid\operatorname{\mathbf{y}}\right\\}\rightarrow
0\text{ in }\operatorname{\mathbb{P}}_{0}^{n}\text{-probability},$
$\displaystyle\Pi_{\text{Mat},n}^{\text{Clas}}\left\\{f:\|f-\mu_{0}\|_{n}>C_{2}\epsilon_{n}\mid\operatorname{\mathbf{y}}\right\\}\rightarrow
0\text{ in }\operatorname{\mathbb{P}}_{0}^{n}\text{-probability},$ (47)

are satisfied for every $\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}$
with $\epsilon_{n}=n^{-\nu_{2}/(2\nu_{1}+d)}$, where $C_{1},C_{2}$ are
positive constants.

  2. 2.

If $\mu_{0}$ is infinitely smooth functions in the sense that
$\mu_{0}\in\mathfrak{A}^{r,\gamma}[0,1]^{d}$ for some $r\geq 1$ and
$\gamma>0$, then as $n\rightarrow\infty$

$\displaystyle\Pi_{\text{SE},n}^{\text{Reg}}\left\\{f:\|f-\mu_{0}\|_{n}>C_{1}\epsilon_{n}\mid\operatorname{\mathbf{y}}\right\\}\rightarrow
0\text{ in }\operatorname{\mathbb{P}}^{n}_{0}\text{-probability},$
$\displaystyle\Pi_{\text{SE},n}^{\text{Clas}}\left\\{f:\|f-\mu_{0}\|_{n}>C_{2}\epsilon_{n}\mid\operatorname{\mathbf{y}}\right\\}\rightarrow
0\text{ in }\operatorname{\mathbb{P}}_{0}^{n}\text{-probability},$ (48)

are satisfied for every $\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}$
with $\epsilon_{n}=n^{-1/2}(\log n)^{\max(1/r,1/2+d/4)}$, where $C_{1},C_{2}$
are positive constants.

###### Proof.

We note that $\Psi$ indexes the GP prior and
$\Psi\in[0,1]^{C|t^{*}||\operatorname{\mathcal{B}}|}$, where
$C,|t^{*}|,|\operatorname{\mathcal{B}}|$ are defined in Theorem A.1;
therefore, we can directly apply Theorem 11.22, Theorem 11.23, and the results
in Section 11.4.4 of Ghosal and van der Vaart [4] for Matern and SE GP priors
indexed by a covariate lying in $[0,1]^{d}$ with
$d=C|t^{*}||\operatorname{\mathcal{B}}|$. Specifically, Lemma 11.36 and Lemma
11.37 imply that (1) is satisfied with
$\epsilon_{n}=n^{-\nu_{2}/(2\nu_{1}+d)}$ for a given
$\operatorname{\mathbf{w}},\operatorname{\mathbf{R}}$. Similarly, Lemma 11.38
and Lemma 11.41 imply that (2) is satisfied with $\epsilon_{n}=n^{-1/2}(\log
n)^{\max(1/r,1/2+d/4)}$. ∎

##  C Sampling Algorithm for Posterior Inference

###  C.1 GP Regression

Recall Equation (11) in the main manuscript. Denote the kernel matrix as
$\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}$ with
$\\{\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}\\}_{ij}=\kappa_{f}(\operatorname{\mathbf{T}}_{i},\operatorname{\mathbf{T}}_{j})$,
where $f\in\\{\text{exp},\text{SE}\\}$ and $i,j\in\\{1,\ldots,n\\}$. Let
$\operatorname{\mathbf{y}}=(y_{1},\ldots,y_{n})^{\top}$,
$\operatorname{\mathbf{X}}$ be the $n\times p$ matrix with
$\operatorname{\mathbf{x}}_{i}^{\mathrm{\scriptscriptstyle T}}$ as its $i$th
row,
$\operatorname{\mathbf{f}}=\\{f(\operatorname{\mathbf{T}}_{1}),\ldots,f(\operatorname{\mathbf{T}}_{n})\\}^{\mathrm{\scriptscriptstyle
T}}$, and
$\operatorname{{\mathbf{\epsilon}}}=(\epsilon_{1},\ldots,\epsilon_{n})^{\mathrm{\scriptscriptstyle
T}}$. Then, the hierarchical Bayesian model for $\operatorname{\mathbf{y}}$ is

$\displaystyle\operatorname{\mathbf{y}}=\operatorname{\mathbf{X}}\operatorname{{\mathbf{\beta}}}+\operatorname{\mathbf{f}}+\operatorname{{\mathbf{\epsilon}}},\quad\operatorname{{\mathbf{\epsilon}}}\mid\sigma^{2},\alpha\sim
N(\operatorname{\mathbf{0}},\sigma^{2}\alpha\operatorname{\mathbf{I}}),\quad(\operatorname{{\mathbf{\beta}}},\sigma^{2})\propto\sigma^{-2},\quad\operatorname{\mathbf{f}}\mid\operatorname{{\mathbf{\theta}}}\sim
N(\operatorname{\mathbf{0}},\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}),$
(49)

where $\operatorname{{\mathbf{\beta}}}$ and $\operatorname{\mathbf{f}}$ are
assumed to be independent apriori, $\alpha=\tau^{2}/\sigma^{2}$ is the inverse
of signal-to-noise ratio, and $\phi\in(a_{\phi},b_{\phi})$, where
$0<a_{\phi}<b_{\phi}$. We assign a prior on $(\phi,\alpha)$ through
$(u_{1},u_{2})\in\operatorname{\mathbb{R}}^{2}$ as

$\displaystyle\phi=a_{\phi}+(b_{\phi}-a_{\phi})/(1+e^{-u_{1}}),\quad\alpha=e^{u_{2}},\quad(u_{1},u_{2})^{\top}\sim
N\\{\operatorname{\mathbf{0}},{\mathrm{diag}}(b_{1},b_{2})\\},$ (50)

where $b_{1},b_{2}>0$. The parameters $\sigma^{2},\phi$ are non-identified in
$\kappa_{\text{exp}},\kappa_{\text{SE}}$, but this does not affect the
inference on $f$ or prediction of the response [16].

First, we derive the full conditional of
$(\sigma^{2},\operatorname{{\mathbf{\beta}}})$ and the likelihood for drawing
$\phi$ and $\alpha$. Marginalizing over $\operatorname{\mathbf{f}}$ in (49)
gives

$\displaystyle\operatorname{\mathbf{y}}=\operatorname{\mathbf{X}}\operatorname{{\mathbf{\beta}}}+\operatorname{{\mathbf{\eta}}},\quad\operatorname{{\mathbf{\eta}}}\sim
N(\operatorname{\mathbf{0}},\sigma^{2}\operatorname{\mathbf{C}}_{yy}),\quad\operatorname{\mathbf{C}}_{yy}=\operatorname{\mathbf{K}}_{\operatorname{{\mathbf{\theta}}}}/\sigma^{2}+\alpha\operatorname{\mathbf{I}}.$
(51)

If $\operatorname{\mathbf{L}}$ is a lower triangular matrix such that
$\operatorname{\mathbf{C}}_{yy}=\operatorname{\mathbf{L}}\operatorname{\mathbf{L}}^{\top}$,
then define

$\displaystyle\tilde{\operatorname{\mathbf{y}}}=\operatorname{\mathbf{L}}^{-1}\operatorname{\mathbf{y}},\quad\tilde{\operatorname{\mathbf{X}}}=\operatorname{\mathbf{L}}^{-1}\operatorname{\mathbf{X}},\quad\hat{\operatorname{{\mathbf{\beta}}}}=(\tilde{\operatorname{\mathbf{X}}}^{\mathrm{\scriptscriptstyle
T}}\tilde{\operatorname{\mathbf{X}}})^{-1}\tilde{\operatorname{\mathbf{X}}}^{\top}\tilde{\operatorname{\mathbf{y}}},\quad\hat{\tilde{\operatorname{\mathbf{y}}}}=\tilde{\operatorname{\mathbf{X}}}\hat{\operatorname{{\mathbf{\beta}}}},\quad\tilde{\operatorname{{\mathbf{\eta}}}}=\operatorname{\mathbf{L}}^{-1}\operatorname{{\mathbf{\eta}}},$
(52)

which reduces (51) to

$\displaystyle\tilde{\operatorname{\mathbf{y}}}=\tilde{\operatorname{\mathbf{X}}}\operatorname{{\mathbf{\beta}}}+\tilde{\operatorname{{\mathbf{\eta}}}},\quad
N(\operatorname{\mathbf{0}},\sigma^{2}\operatorname{\mathbf{I}}).$ (53)

The full conditional of $(\sigma^{2},\operatorname{{\mathbf{\beta}}})$ given
the rest is

$\displaystyle\sigma^{2},\operatorname{{\mathbf{\beta}}}\mid\operatorname{\mathbf{y}},\alpha,\phi\propto(2\pi)^{-n/2}(\sigma^{2})^{-n/2-1}e^{-\frac{1}{2\sigma^{2}}\|\tilde{\operatorname{\mathbf{y}}}-\hat{\tilde{\operatorname{\mathbf{y}}}}\|^{2}}e^{-\frac{1}{2\sigma^{2}}(\operatorname{{\mathbf{\beta}}}-\hat{\operatorname{{\mathbf{\beta}}}})^{\top}(\tilde{\operatorname{\mathbf{X}}}\tilde{\operatorname{\mathbf{X}}})(\operatorname{{\mathbf{\beta}}}-\hat{\operatorname{{\mathbf{\beta}}}})};$
(54)

therefore, we generate $\sigma^{2}\mid\operatorname{\mathbf{y}},\alpha,\phi$
and
$\operatorname{{\mathbf{\beta}}}\mid\sigma^{2},\operatorname{\mathbf{y}},\alpha,\phi$
as