# Ensemble Prediction of Time to Event Outcomes with Competing Risks: A Case
Study of Surgical Complications in Crohn’s Disease

Michael C Sachs  Department of Medicine, Karolinska Institutet Stockholm,
Sweden<EMAIL_ADDRESS>Andrea Discacciati  Institute of Environmental
Medicine, Karolinska Institutet, Stockholm, Sweden.  Åsa H Everhov  Department
of Medicine, Karolinska Institutet, Stockholm, Sweden.  Ola Olén  Department
of Medicine, Karolinska Institutet, Stockholm, Sweden.  Erin E Gabriel
Department of Medical Epidemiology and Biostatistics, Karolinska Institutet,
Stockholm, Sweden.

###### Abstract

We develop a novel algorithm to predict the occurrence of major abdominal
surgery within 5 years following Crohn’s disease diagnosis using a panel of 29
baseline covariates from the Swedish population registers. We model pseudo-
observations based on the Aalen-Johansen estimator of the cause-specific
cumulative incidence with an ensemble of modern machine learning approaches.
Pseudo-observation pre-processing easily extends all existing or new machine
learning procedures to right-censored event history data. We propose pseudo-
observation based estimators for the area under the time varying ROC curve,
for optimizing the ensemble, and the predictiveness curve, for evaluating and
summarizing predictive performance.
Keywords: competing risks; Crohn’s disease; inflammatory bowel disease;
machine learning; pseudo observations

##  1 Introduction

###  1.1 Motivating study and statistical approaches

Crohn’s disease (CD) is a chronic debilitating condition characterized by
periods of inflammatory activity in the bowel that causes symptoms such as
abdominal pain, diarrhea, and weight loss. Pharmacologic treatment for CD
includes medications such as steroids, immunomodulating drugs, and biological
therapy. Despite these available medications, many people with CD are
escalated to surgical interventions from small to extensive resections of the
bowel or colon (Gomollón et al., 2016). Previous studies have estimated that
up to 50% of patients with CD undergo surgery within 10 years after diagnosis;
however, surgical rates have decreased over time, possibly due to the
introduction of modern treatments such as thiopurines and anti-TNF (Lakatos et
al., 2012; Ramadas et al., 2010). The aim of this study is to determine
whether clinical and demographic characteristics observed at the time of
diagnosis can be used to predict the occurrence of major abdominal surgery
within 5 years, with the goal of personalized disease management. Although it
might be of interest to determine which characteristics are associated with
risk, we instead focus on obtaining the most accurate predictions, which, at
least in this registry based setting, are useful even without knowing what is
driving them.

One of the statistical challenges is to deal with the fact that patients may
die or emigrate during follow-up. These competing risks must be accounted for
in development and evaluation of the prediction model. There are several
possible approaches available for modeling a censored time to event outcome
with competing risks. The time to event can be converted into a binary
indicator of having the event within 5 years, and then the remaining
observations weighted by the inverse of the probability of not being censored.
This can be less efficient as subjects censored prior to the time horizon do
not directly contribute to estimation of the model. What’s more, the
estimation method must be able to incorporate weights. A series of Cox models
or direct modeling of the cumulative incidence as a function of covariates is
also an option via the Fine-Gray model (Fine and Gray, 1999).

As it is rarely known which model will give the best predictions, we instead
wish to use an ensemble of models and machine learning (ML) procedures for
prediction. There are existing methods that adapt ML procedures to censored
data. Binder et al. (2009) develops a boosting algorithm for the Fine-Gray
model under competing risks to select a subset of variables and form a
prediction model. Ishwaran et al. (2014) develops a survival random forest
that accounts for competing risks. Hothorn et al. (2006) develop two methods
for ensemble learning that deal with right censored data by considering only
observed event times and weighting by the probability of not being censored.
Scheike et al. (2008) describes an alternative approach that involves direct
binomial regression of the counting process formulation of the event in the
presence of competing risks. Gerds and Ozenne (2018) provides a suite of
methods for estimating the performance of prediction models in the competing
risks setting, using inverse probability of censoring weighted (IPCW)
estimates of the Brier score and AUC, however, one must fit the prediction
models first. Under many of the existing methods based on binary data,
including those in Scheike et al. (2008) and Hothorn et al. (2006), one needs
to model the distribution of time to censoring and then use these as weights
in the ML procedure. Although weight estimation may be easy, particularly in
our motivating setting, there still needs to be a way to use the IPC weights
in the ML procedure, making IPC weighting a useful but procedure specific
solution for dealing with censoring.

Our suggested procedure provides a general method that allows the application
of all existing and any newly developed ML algorithm, without extension, to
censored data with or without competing risks. We propose to generate pseudo-
observations of the cause-specific cumulative incidence of surgery at 5 years,
accounting for the competing risks of death and emigration, and use these
pseudo-observations as the outcome variable in an ensemble learning approach
using standard software for continuous outcomes. Mogensen and Gerds (2013)
takes a similar approach to our proposed method, also using pseudo-
observations, but only for the class of tree-based methods, building an
ensemble of trees to construct a random forest. One may view our proposed
procedure as an extension of concepts in Mogensen and Gerds (2013).

In addition, we demonstrate how pseudo-observations can be used to estimate
the time-varying ROC curve (Zheng et al., 2012; Heagerty et al., 2000) and the
time varying predictiveness curve (Pepe et al., 2007). To our knowledge,
neither of these results has been demonstrated before. Using these results, we
combine a set of ML procedures into an ensemble by optimizing the time-varying
area under the ROC curve.

The use of pseudo-observations has been described and evaluated fairly
extensively for their use in the estimation of association models (Andersen et
al., 2003; Andersen and Pohar Perme, 2010). Relatively little attention has
been given towards their use in the development and evaluation of prediction
models. In addition to Mogensen and Gerds (2013), Nicolaie et al. (2013)
develop methods related to the development of dynamic prediction models using
pseudo-observations and one paper has suggested evaluating existing predictive
models using pseudo-observations to estimate the Brier score (Cortese et al.,
2013).

###  1.2 Description of the dataset and variables

The motivating example comes from the Swedish population registers, in which
we identified all individuals with a first ever diagnosis of CD in the
Inpatient or Outpatient register from January 1, 2003 until December 31, 2014
using the international classification of disease (ICD) code K50 (ICD-10,
1997-2014). We required at least two diagnostic listings with CD as main or
contributory diagnosis in inpatient or outpatient care or CD associated
abdominal surgery at the same time as the first CD diagnosis. Patients
registered with a record of one of the major abdominal surgery events before 1
January 2003 were excluded.

The outcome variable of major abdominal surgery was defined according to a set
of procedure codes based on the NOMESCO Classification of Surgical Procedures.
We used six categories of procedures based on anatomical localization and
expected functional outcome: resection of the small bowel, ileocecal/ileocolic
resection, segmental resection of the colon, colectomy, and/or proctectomy.
The occurrence of any of those surgeries and possibly in combination with each
other defined the major abdominal surgery outcome. Endoscopic, perianal, and
upper gastrointestinal procedures were not included.

Baseline characteristics and demographics are gathered from the patient
register, the multigeneration register, and the prescribed drug register which
are described in more detail elsewhere (Ludvigsson et al., 2016). A detailed
description of the data sources and methods for defining variables is
available in a recent manuscript (Olén et al., 2017), and descriptive
statistics for our study population are given in the Supplementary Materials.

###  1.3 Aims of this paper

Our aims are four-fold: illustrate our proposed general pseudo-observation
based method for ensemble learning optimizing the area under the time varying
ROC curve, compare our approach to existing methods, develop prediction
algorithms for risk of surgery within 5 years in patients with CD, and
evaluate those predictions with the pseudo-observation based predictiveness
curve.

After giving some notation in Section 2.1, we describe pseudo-observations for
the cause-specific cumulative incidence in Section 2.2. We then develop
pseudo-observation based estimators for the time-varying cause-specific area
under the ROC curve and the predictiveness curve in Section 2.3. Then, using
these results, we detail our modification of the SuperLearner algorithm for
this setting in Section 2.4. All of these methods are evaluated and compared
to existing approaches in a simulation study which is described in Section 3.
Finally, we analyze and discuss the results from our motivating study in
Section 4.

##  2 Methods

###  2.1 Notation

  * •

$Y_{i}$: right censored event time and $\Delta_{i}\in\\{0,1,2,3\\}$: event
indicator where 0 indicates censoring, 1 indicates surgery, 2 indicates death,
and 3 indicates emigration, for subjects $i=1,\ldots,n$.

  * •

$T_{i}$: is the true event time, $\delta_{i}\in\\{1,2,3\\}$: event indicator 1
indicates surgery, 2 indicates death, and 3 indicates emigration,

  * •

$X_{ij}$: biomarker measurement/covariate measured at baseline, for
$j=1,\ldots,p$. $\mathbf{X_{i}}$ is the vector of covariates for subject $i$.

  * •

$B_{ik}=f_{k}(\mathbf{X_{i}})$, where
$f_{k}:\mathbb{R}^{p}\rightarrow\mathbb{R}$ is a covariate signature indexed
by $k$. We denote $\hat{f}_{k}$ the estimated signature. The index $k$ is used
to distinguish between signatures estimated by different algorithms or on
different subsets of the data.

We are interested in the cumulative incidence of failure due to surgery at 5
years:

$C_{1}(t_{\star})=E\\{I(T_{i}\leq
t_{\star},\delta_{i}=1)\\}=\int_{0}^{t_{\star}}S(u)\alpha_{1}(u)\,du,$

where $S(u)$ is the overall survival probability, $\alpha_{1}(u)$ is the
cause-specific hazard for failures from abdominal surgery, and $t_{\star}=5$
years.

###  2.2 Pseudo-observations for cumulative incidence

Following Andersen and Pohar Perme (2010), we summarize the survival data as
the counting process:

$N_{1}(t)=\sum_{i}I(Y_{i}\leq t,\Delta_{i}=1),$

giving the number of observed failures due to surgery on or before time $t$,
and $R(t)=\sum_{i}I(Y_{i}\geq t)$ gives the number of subjects still at risk
just before time $t$. Our estimator of the cumulative incidence function is
the Aalen-Johansen estimator (Aalen and Johansen, 1978)

$\hat{C}_{1}(t)=\int_{0}^{t}\hat{S}(u)\,d\hat{A}_{1}(u),$

where $\hat{A}_{1}(u)=\int_{0}^{t}\,dN_{1}(u)/R(u)$ is the Nelson-Aalen
estimator for the cumulative cause-specific hazard for surgery failures, and
$\hat{S}$ is the Kaplan-Meier estimator of the overall survival from any
cause.

Then the $i$th pseudo-observation for cause 1 at time $t$ is

$\hat{C}^{i}_{1}(t)=n\hat{C}_{1}(t)-(n-1)\hat{C}_{1}^{-i}(t),$

where $\hat{C}^{-i}_{1}(t)$ is the cumulative incidence function estimator
computed by using the sample excluding the $i$th observation. By construction
and the unbiasedness of the Aalen-Johansen estimator, the pseudo-observations
are unbiased for the cumulative cause-specific incidence:
$E\\{\hat{C}^{i}_{1}(t)\\}=C_{1}(t).$ Moreover, observe that the survival from
any cause can be estimated
$E\\{1-\sum_{j=1}^{3}\hat{C}^{i}_{j}(t)\\}=P(Y_{i}>t)$.

Graw et al. (2009) proved that the mean of pseudo-observations conditional on
covariates is asymptotically unbiased for the conditional cause-specific
cumulative incidence. The validity of this proof and asymptotic properties of
pseudo-observations in regression settings was further studied in Overgaard et
al. (2017) and Jacobsen and Martinussen (2016). Thus we have,

$E(\hat{C}^{i}_{1}(t_{\star})|\mathbf{X_{i}})=P(T_{i}\leq
t_{\star},\delta_{i}=1|\mathbf{X_{i}})+o_{P}(1).$

This property of the pseudo-observations for the cumulative incidence allows
us to use them as the outcome in prediction models in which we aim to estimate
a mapping $f(\mathbf{X_{i}})$. Then, as we show in the next section, we can
evaluate the mapping $f$ for predicting surgery by estimating the time-varying
ROC curve, again using pseudo-observations.

Conditional asymptotic unbiasedness of the pseudo-observations relies on the
assumption that censoring time is independent of event time, event type and
all other covariates; this was called ‘completely independent censoring’ in
Overgaard et al. (2017). In our dataset, due to the quality of the Swedish
population register, there is no loss to follow-up that is not due to death or
emigration. Thus it is critical to model death and emigration as competing
risks rather than simply censoring those subjects because the association
between certain covariates and death or emigration may violate this
assumption. All other censoring before 5 years is determined by the time of
entry into the study (the date of CD diagnosis). For this reason, we use
stratified pseudo-observation procedures. This was suggested in Andersen and
Pohar Perme (2010), as a way to relax the deal with censoring that is
dependent on a known categorical covariate. In other studies with loss to
follow-up, censoring dependence maybe more complex. Procedures to extend
pseudo-observations by modeling the censoring distribution and reweighting
have been developed (Binder et al., 2014).

###  2.3 Time varying measures of prediction accuracy estimated using pseudo
observations

The time-varying ROC curve (Heagerty et al., 2000; Saha and Heagerty, 2010)
for biomarker signature $B_{i}$ at time $t$ and cutoff $c$ is defined by the
time varying, cause-specific true-positive fraction

$TP(t,c)=P(B_{i}>c|T_{i}\leq t,\delta_{i}=1)$

and the false positive fraction that is defined across all event types:

$FP(t,c)=P(B_{i}>c|T_{i}>t).$

While the $TP$ can be defined for all causes, we omit the cause-type subscript
for the $TP$ because we are only interested in the failures due to cause 1
(abdominal surgery). We can estimate these using pseudo-observations as
follows. Bayes’ rule gives us

$P(B_{i}>c|T_{i}\leq t,\delta_{i}=1)=\frac{P(T_{i}\leq
t,\delta_{i}=1|B_{i}>c)P(B_{i}>c)}{P(T_{i}\leq t,\delta_{i}=1)},$

which is the ratio of a conditional cumulative incidence to the marginal
cumulative incidence. This suggests the following estimators:

$\widehat{TP}(t,c)=\frac{\sum_{i}\hat{C}^{i}_{1}(t)I(B_{i}>c)}{\sum_{i}\hat{C}^{i}_{1}(t)}$

and

$\widehat{FP}(t,c)=\frac{\sum_{i}(1-\sum_{j=1}^{3}\hat{C}^{i}_{j}(t))I(B_{i}>c)}{\sum_{i}1-\sum_{j=1}^{3}\hat{C}^{i}_{j}(t)}.$

We conjecture that these estimators are asymptotically unbiased, based on the
properties of the pseudo-observations shown in Overgaard et al. (2017) and
application of the continuous mapping theorem.

The ROC curve at time $t$ is a plot of the pairs $\\{FP(t,c),TP(t,c)\\}$ as
the cutoff $c$ varies. The estimated area under the time varying ROC curve at
a fixed time $t$ for signature $B$ is computed numerically using the
trapezoidal rule, and is denoted $\widehat{AUC}(B,t)$. To our knowledge, this
is a novel estimator of the time-varying ROC curve and, as a result, the time-
varying AUC. The time-varying AUC has recently been shown to be a proper
scoring function in the context of risk prediction of events at a fixed time,
unlike the survival concordance index (Blanche et al., 2018), making it a
sensible target for prediction. The Brier score is a proper scoring rule also,
but is difficult to interpret in this case because the pseudo-observations are
unbounded and because the true binary outcome is not completely observed. The
time-varying predictiveness curve is then used to comprehensively evaluate the
performance of the resulting prediction models.

We can also estimate a time varying version of the predictiveness curve (Pepe
et al., 2007) using pseudo-observations. The cause-specific predictiveness
curve can be defined $P(T_{i}\leq t,\delta_{i}=1|B_{i}=c).$ We assess the
usefulness of the model that creates $B_{i}$ by plotting the percentiles of
$B_{i}$ versus the estimated cause specific cumulative incidence at those
percentiles. As described by Pepe et al. (2007), this curve displays whether
there are meaningful variations in the risk of the event in question over the
range of predicted quantities. To our knowledge, pseudo-observations have not
been used to estimate predictiveness curves previously.

This curve could be estimated by binning $B_{i}$, into deciles and taking the
average of the pseudo-observations within each of the bins. Further, we can
fit a locally weighted smooth average with the pseudo-observations as the
outcome and the quantiles of $B_{i}$ as the predictor. The fitted values from
this model provide an estimate of the predictiveness curve, from which we can
display the estimated population cumulative incidence values as a smooth
function of the predicted probabilities.

The estimated true positive fractions and false positive fractions plotted as
functions of the predicted risk values, combined with ROC and predictiveness
curves provide a comprehensive summary of the performance and operating
characteristics of the predictive models. We illustrate this in the data
analysis below in Figure 3.

###  2.4 Super-Learner predictive model development

It is unknown which class of prediction procedures will perform best on a
given real dataset. Instead of selecting one specific procedure or machine
learning method, we use an ensemble of methods, and then form predictions by
stacking the predictions from the particular methods into an ”ensemble
predictor”. The stacking is performed by selecting the linear combination of
individual predictions that maximizes the estimated $AUC(t)$ in a cross-
validation procedure. This stacked modeling approach is called the
Superlearner (Van der Laan et al., 2007), and is implemented in the
SuperLearner R package (Polley et al., 2017). Van der Laan et al. (2007)
demonstrate that the resulting prediction algorithm performs at least as well
as the best predictor in the library of algorithms. The general algorithm is
described in that paper, and here we present the specific steps adapted to our
setting.

Let $\mathcal{L}$ denote the library of $K$ prediction algorithms. In general,
a prediction algorithm with index $k$ for $k=1,\ldots,K$, could be any method
to develop $\hat{f}_{k}$. The specific algorithms that we use in our
simulations and data example are given in Table 3. The steps are as follows:

  1. 1.

For a fixed, finite grid of time points
$\\{t_{1},\ldots,t_{\star},\ldots,t_{m}\\}$, where $t_{\star}$ is the time
point of interest for predicting incidence, calculate the pseudo-observations
$\hat{C}^{i}_{j}(t_{l})$, for $j=1,2,3;l=1,\ldots,m$; $i=1,\ldots,n$.

  2. 2.

Split the dataset $\mathcal{X}$ into a training and validation sample,
according to a V-fold cross validation scheme: split the ordered $n$
observations into $V$ equal size groups. Let the $\nu$-th group be the
validation sample, and the remaining group the training sample, $\nu=1,...,V$.
Define $T(\nu)$ to be the $\nu$th training data split and $V(\nu)$ to be the
corresponding validation data split where $T(\nu)=\mathcal{X}\setminus
V(\nu),\nu=1,...,V$.

  3. 3.

Stack the datasets so that there are $m$ rows per subject $i$, with the same
covariates $\mathbf{X}_{i}$, additional variables for the pseudo-observations,
and an additional variable for the times
$t_{1},\ldots,t_{\star},\ldots,t_{m}$.

  4. 4.

Fit each algorithm on the training set for each fold $T(\nu)$, to obtain
$\hat{f}_{k\nu}$ for $k=1,\ldots,K$ and $\nu=1,...,V$.

  5. 5.

For each $\hat{f}_{k\nu}$, obtain predicted pseudo-observations at $t_{\star}$
for the validation set for each fold, i.e.,
$\tilde{C}^{i}_{1}(t)_{k}=\hat{f}_{k\nu}(X_{i},t_{\star})$, for $i\in V(\nu)$
and $k=1,\ldots,K$. This yields $K$ predicted pseudo-observations for each
subject $i$. Denote the vector of these for subject $i$ and time $t_{\star}$
as $\mathbf{\tilde{C}}^{i}_{1}(t_{\star})$.

  6. 6.

For a fixed $\lambda$ and $t_{\star}$, find the $K$ dimensional parameter
vector $\hat{\alpha}$ that satisfies

$\max_{\alpha}\widehat{AUC}(\alpha^{T}\mathbf{\tilde{C}}^{i}_{1}(t_{\star}),t_{\star})+\lambda\sum_{k}|\alpha_{k}|.$

To ensure proper calibration, normalize the coefficients by their sum,
$\hat{\alpha}^{*}=\hat{\alpha}/\mathbf{1}^{T}\hat{\alpha}$. Although multiple
time points may be used in the prediction procedures, note that only the time
point $t_{\star}$ is used in the loss function optimization procedure.

  7. 7.

For $k=1,\dots,K$, fit each algorithm on the full dataset to obtain
$\hat{f}_{k}(\mathbf{X}_{i},t_{\star})$, for $i=1,\dots,n$. This forms an $n$
by $k$ matrix where the rows are predicted pseudo-observations for each
subject, and the columns correspond to different prediction methods. The final
predictor is the linear combination $\hat{\alpha}^{*T}\hat{f}_{k}$ according
to the cross-validated weights obtained in the previous step.

The penalty term in step 6 above is introduced to solve the indentifiability
issues that arise in maximization of the AUC, similar to what was suggested by
Lin et al. (2011) and further discussed by Fong et al. (2016). Briefly, the
issues in maximizing the empirical AUC arise because the ROC curve is
invariant to monotone transformations of the signature. Penalization of the
objective function is one way to solve this issue and yield a unique solution
for $\alpha$.

The calculation of the pseudo-observations at the fixed grid of time points in
step (a) and the stacking thereof allows us to borrow information across the
time points, with the idea that this should perform as well or better than the
counting process approach, or using an single time point. The computational
burden of both calculating the pseudo-observations and fitting the models in
the library increases with the number of time points. The selection of the
density and location of the time points in the grid is not something we
investigate in this paper, however, in principle, one can include all observed
event times. We suggest possible approaches to this issue in the discussion.

##  3 Simulation study

###  3.1 Data generation

For each replicate of the simulation scenario, we generate a dataset with 500
subjects, 20 normally distributed variables that are correlated with each
other to varying degrees, and event times that are distributed Weibull with
about 20% cumulative incidence of the surgery event and 7% cumulative
incidence of the competing event of death at time $t_{\star}$. Censoring is
completely independent (except for scenario D) and is either 20% or 50% on
average random censoring, corresponding to 33% or 65% at the time $t_{\star}$.
In scenario D, we allow censoring to depend on covariates. In all scenarios,
the competing event was generated according to a Weibull distribution with a
scale parameter that depends on a single covariate that was not related to the
covariates that are associated with the time to surgery.

For each simulated dataset, we implement our proposed method by computing
pseudo observations for a grid of 4 time points, 2 before and 1 after the time
point of interest, $t_{\star}$ and using the AUC based SuperLearner procedure
desribed in detail above. For comparison to our proposed method, we fit a
pseudo-observation model using a single time point at $t_{\star}$, a IPCW
SuperLearner model on the binary outcome created by dichotomizing the time to
surgery at $t_{\star}$ using the non-negative log likelihood loss, the model
proposed by Binder et al. (2009) as implemented in the CoxBoost R package, and
the random forests for competing risks survival proposed by (Ishwaran et al.,
2014), as implemented in the R package randomForestsSRC. We estimate and use
inverse probability of censoring weighting in the binary SuperLearner approach
via Kaplan-Meier for censoring event. For the binary SuperLearner only methods
that natively allowed for weights were considered, while more libraries were
used by the pseudo-observations. Table 3 in the application section describes
the libraries used for each SuperLearner type. We generate an independent
dataset on which to estimate the out-of-sample prediction error as estimated
by the AUC.

We consider 2 different censoring rates and 4 different scenarios:

  * •

Scenario A: exactly one variable is associated with the scale parameter of the
surgery outcome in a log-linear way.

  * •

Scenario B: exactly 5 variables are associated with the scale parameter of the
surgery outcome according to a b-spline relationship of degree 3, and with a
linear interaction effect.

  * •

Scenario C: exactly 5 variables are associated with the scale parameter of the
surgery outcome in a complex way, including non-linearities, interactions, and
threshold effects.

  * •

Scenario D: Same as B, but censoring depends on a subset of the covariates
from those associated with the surgery outcome. The covariate effect on
censoring is log-linear on the scale parameter of a Weibull distribution.

  * •

Scenario E: Bootstrap of 500 samples and the first 20 covariates from the CD
data analysis. In this case the true event times and probabilities are
unknown.

Detailed specification of and code for the simulation scenarios is given in
the supplementary materials. In addition, we released an R package called
sachsmc/pseupersims on github that contains the code used to perform this
simulation study.

The purpose of the simulation study is to evaluate the ability of the
SuperLearner model based on pseudo-observations to develop a prediction model
for the cumulative incidence and to compare to alternative approaches, and to
the true model.

###  3.2 Using pseudo-observations to develop prediction models

Our proposed approach performs comparably to competing methods in terms of
accurately predicting the true cumulative incidence, as shown in Tables 1 and
2. The average AUC over the simulation replicates is superior to all other
considered approaches in scenarios B, C, and D, while the other approaches
perform well in scenario A. Compared to the binary SuperLearner approach,
differences in the performance of the pseudo-observation based models were
small, while the CoxBoost model and random forests occasionally showed
markedly inferior performance. Still the classification accuracy of the pseudo
observation model is competitive with other approaches and potential for
performance gains are clear. In the simulation E, 500 person bootstraps from
the CD analysis data, we present the average estimated AUC for each of the
methods. In that case, even the model based pseudo observations estimated at a
single time point outperforms the binary, but the CoxBoost and random forests
models outperform all other methods.

The models based on pseudo observations computed at a single time point
occasionally had convergence issues and therefore have extreme outliers, which
explains the sometimes dramatically high MSE. The pseudo-observation model
based on four time points had smaller variance in the predictions relative to
the binary model in scenarios D and E. However, in some cases such as scenario
B, we observe a higher variance, because, unlike all of the other methods, the
predictions from the pseudo-observation models are not bounded. Finally, our
proposed estimator for the AUC using the pseudo observations appears to be
unbiased, as the estimated values were nearly identical to the AUC estimated
using the true, and typically unobserved binary outcome of surgery at time
$t_{\star}$.

Table 1: Results of simulation study using a setting with approximately 20% random censoring (approximately 33% missingness in the binary outcome). We compare the estimation of the area under the time varying ROC curve using the true binary outcome that is unobserved (labelled tbauc), the estimation of the AUC using pseudo observations as described in Section 2 (labelled pauc), and the bias, standard deviation, and mean squared error (mse) of the predicted outcomes compared to the true probabilities of the outcome.  model | mean.tbauc | sd.tbauc | mean.pauc | sd.pauc | mean.bias | sd.prob | mse
---|---|---|---|---|---|---|---
| Scenario A
binary | 0.735 | 0.031 | 0.736 | 0.033 | 0.003 | 0.146 | 0.043
CoxBoost | 0.734 | 0.026 | 0.731 | 0.029 | 0.002 | 0.116 | 0.040
pseudo | 0.729 | 0.028 | 0.729 | 0.031 | 0.002 | 0.147 | 0.086
pseudo.single | 0.717 | 0.033 | 0.718 | 0.035 | 0.000 | 0.171 | 0.132
rfsrc | 0.716 | 0.027 | 0.717 | 0.029 | -0.008 | 0.123 | 0.049
true | 0.747 | 0.025 | 0.747 | 0.028 | 0.000 | 0.228 | 0.000
| Scenario B
binary | 0.569 | 0.042 | 0.570 | 0.045 | -0.001 | 0.064 | 0.068
CoxBoost | 0.522 | 0.034 | 0.514 | 0.038 | 0.009 | 0.041 | 0.044
pseudo | 0.546 | 0.073 | 0.548 | 0.075 | -0.056 | 0.505 | 103.615
pseudo.single | 0.533 | 0.051 | 0.535 | 0.054 | -0.009 | 0.527 | 5.799
rfsrc | 0.579 | 0.036 | 0.580 | 0.040 | -0.011 | 0.072 | 0.055
true | 0.666 | 0.031 | 0.665 | 0.035 | 0.000 | 0.122 | 0.000
| Scenario C
binary | 0.839 | 0.039 | 0.840 | 0.041 | 0.010 | 0.256 | 0.062
CoxBoost | 0.630 | 0.032 | 0.622 | 0.036 | 0.002 | 0.081 | 0.025
pseudo | 0.848 | 0.062 | 0.848 | 0.062 | 0.010 | 0.268 | 0.257
pseudo.single | 0.833 | 0.040 | 0.834 | 0.042 | -0.007 | 0.300 | 0.117
rfsrc | 0.780 | 0.034 | 0.780 | 0.036 | -0.010 | 0.109 | 0.034
true | 0.927 | 0.013 | 0.927 | 0.013 | 0.000 | 0.521 | 0.000
| Scenario D
binary | 0.839 | 0.037 | 0.843 | 0.037 | 0.009 | 0.258 | 0.054
CoxBoost | 0.625 | 0.032 | 0.621 | 0.034 | -0.002 | 0.081 | 0.027
pseudo | 0.853 | 0.035 | 0.858 | 0.036 | 0.011 | 0.252 | 0.155
pseudo.single | 0.830 | 0.075 | 0.834 | 0.075 | 0.010 | 0.380 | 4.797
rfsrc | 0.782 | 0.032 | 0.785 | 0.033 | -0.013 | 0.107 | 0.042
true | 0.926 | 0.013 | 0.929 | 0.014 | 0.000 | 0.520 | 0.000
| Scenario E: bootstrap of CD data, 36% missingness
binary | - | - | 0.525 | 0.035 | - | 0.324 | -
CoxBoost | - | - | 0.585 | 0.043 | - | 0.062 | -
pseudo | - | - | 0.560 | 0.059 | - | 0.297 | -
pseudo.single | - | - | 0.552 | 0.061 | - | 0.184 | -
rfsrc | - | - | 0.598 | 0.042 | - | 0.089 | -
Table 2: Results of simulation study using a setting with approximately 50% random censoring (approximately 65% missingness in the binary outcome). We compare the estimation of the area under the time varying ROC curve using the true binary outcome that is unobserved (labeled tbauc), the estimation of the AUC using pseudo observations as described in Section 2 (labeled pauc), and the bias, standard deviation, and mean squared error (mse) of the predicted outcomes compared to the true probabilities of the outcome.  model | mean.tbauc | sd.tbauc | mean.pauc | sd.pauc | mean.bias | sd.pred |  $100*$mse
---|---|---|---|---|---|---|---
| Scenario A
binary | 0.730 | 0.034 | 0.728 | 0.039 | 0.024 | 0.140 | 0.114
CoxBoost | 0.734 | 0.028 | 0.729 | 0.036 | 0.026 | 0.110 | 0.117
pseudo | 0.722 | 0.041 | 0.719 | 0.048 | 0.037 | 0.188 | 8.272
pseudo.single | 0.704 | 0.041 | 0.701 | 0.045 | 0.031 | 0.161 | 0.168
rfsrc | 0.711 | 0.031 | 0.708 | 0.038 | 0.017 | 0.110 | 0.077
true | 0.747 | 0.028 | 0.745 | 0.035 | 0.000 | 0.230 | 0.000
| Scenario B
binary | 0.552 | 0.043 | 0.554 | 0.051 | 0.020 | 0.061 | 0.123
CoxBoost | 0.519 | 0.033 | 0.513 | 0.042 | 0.032 | 0.042 | 0.144
pseudo | 0.532 | 0.065 | 0.534 | 0.069 | 0.013 | 0.382 | 6.204
pseudo.single | 0.517 | 0.044 | 0.517 | 0.051 | 0.122 | 0.708 | 96.630
rfsrc | 0.557 | 0.035 | 0.558 | 0.045 | 0.013 | 0.067 | 0.068
true | 0.661 | 0.029 | 0.659 | 0.037 | 0.000 | 0.121 | 0.000
| Scenario C
binary | 0.814 | 0.039 | 0.802 | 0.043 | 0.033 | 0.231 | 0.182
CoxBoost | 0.620 | 0.034 | 0.607 | 0.044 | 0.044 | 0.073 | 0.232
pseudo | 0.820 | 0.057 | 0.811 | 0.057 | 0.075 | 0.319 | 15.444
pseudo.single | 0.792 | 0.074 | 0.780 | 0.075 | 0.045 | 0.320 | 3.285
rfsrc | 0.738 | 0.036 | 0.726 | 0.046 | 0.035 | 0.089 | 0.161
true | 0.926 | 0.013 | 0.911 | 0.017 | 0.000 | 0.521 | 0.000
| Scenario D
binary | 0.822 | 0.042 | 0.812 | 0.047 | 0.026 | 0.234 | 0.139
CoxBoost | 0.616 | 0.033 | 0.602 | 0.040 | 0.031 | 0.075 | 0.139
pseudo | 0.824 | 0.040 | 0.814 | 0.043 | 0.032 | 0.322 | 0.403
pseudo.single | 0.788 | 0.086 | 0.781 | 0.087 | 0.016 | 0.432 | 6.134
rfsrc | 0.738 | 0.035 | 0.727 | 0.039 | 0.023 | 0.090 | 0.090
true | 0.925 | 0.013 | 0.911 | 0.017 | 0.000 | 0.522 | 0.000

##  4 CD complications results

The total sample size included 9605 subjects diagnosed with adult-onset CD
after January 1, 2003. The demographic characteristics we considered for
prediction include gender, age at diagnosis, calendar year of diagnosis,
highest parental education, whether or not parents were born outside of
Sweden, whether the subject was born outside of Sweden, and whether there is
any family history of CD or ulcerative colitis. Clinical characteristics
considered at the time of diagnosis are whether it was an inpatient or
outpatient diagnosis, Montréal classification of disease, presence of primary
sclerosing cholangitis, presence of extra-intestinal manifestations, and
exposure to CD related medications. Comorbidities considered were diabetes,
hypertension, ischemic heart disease, cerebrovascular disease, congestive
heart failure, chronic obstructive pulmonary disorder, and kidney failure. The
covariates are described in Supplemental Tables 1 - 3.

Figure 1: Cumulative incidence curves of surgery, death or emigration since
onset of Crohn’s disease in a sample of 9605 patients with onset between 2003
and 2014.

The cumulative incidence curves for the two competing risks are shown in
Figure 1. We are interested in developing a model to predict the incidence of
abdominal surgery within 5 years, while accounting for the competing risks of
death and emigration.

###  4.1 Predictive performance

We randomly split the dataset into two, equal-size subsamples, a training and
a validation set. On the training set, we performed the algorithm as described
above, including cross-validation within the training sample. Time of entry
into the study (CD diagnosis) was categorized into the sets January
2003-December 2008, and then in 6 month periods from January 2009 to December
2015. Pseudo observations were calculated conditionally on the categorized
time of entry to account for censoring that is determined by time of entry and
for subsets with no subjects with more than 5 years of follow up, the last
estimated value of the cumulative incidence was carried forward. We used the
grid of time points $\\{2.5,4,5,5.5\\}$ years, and calculated pseudo
observations at these times. We ran one SuperLearner model using the full
grid, and another one including all nested subsets of the grid (meaning no
non-adjacent times were included) combined with each ML procedure in the
ensemble. For comparison, we ran all methods considered in the simulations.
The binary outcome had 36% missingness on the training sample. The censoring
weights used in the binary method were estimated using Kaplan-Meier for
censoring, stratified by the times of entry in the same way as for the pseudo
observations.

The performance of the candidate learning algorithms are described in Table 3
by the cross-validated AUCs and the estimated SuperLearner coefficients, when
applicable. Precise specification of the various settings and tuning
parameters for each procedure are supplied in the supplementary materials. The
performance of each of the models is modest, with AUCs ranging from 0.57 to
0.64. The non-negative log-likelihood risk is also presented in Table 3 for
the binary methods; smaller values of this indicate better performance. The
weights estimated using the SuperLearner ensemble optimization are given in
the weight columns. The performance of the ensemble over all nested subsets of
the times in the grid was similar, but slightly superior, to that using only
the four point time grid.

The out-of-sample predictive performance results are shown in Figure 3. The
SuperLearner model using the pseudo-observations was superior to the weighted
binary SuperLearner, with estimated AUCs of 0.64 and 0.50, respectively. The
random forests model and the CoxBoost model performed slightly worse than the
pseudo-observation approach, as summarized by the AUC, with an AUC of 0.632
and 0.60, respectively. The characteristics of the models differed in other
ways, such as the shape of the ROC curve and predictiveness curves.

The predictiveness curves in Figure 3 show that the pseudo-observation model
is informative regarding risk of abdominal surgery. CD patients in the top
percentile of the predicted risk have an estimated 21% probability of surgery
within 5 years, compared to the marginal risk of 13%. Likewise, the patients
with the lowest predicted risk have less than a 5% risk of surgery within 5
years. The smoothed locally weighted average of the pseudo observations allows
us to estimate the predictiveness curve as a smooth function of the risk
percentile, as shown in the solid line. This also suggests that pseudo-
observation based predictions have value, as it varies monotonically over the
percentiles of the predictions.

Table 3:  Description of the prediction algorithms used and in which settings
they were used for the example. For the binary and pseudo observations
SuperLearner approaches, we summarize the non-negative log-likelihood risk for
binary (smaller is better), and the AUC (larger is better) for pseudo for each
method, along with the estimated coefficients for the method in the ensemble.
The nested pseudo column presents the average (over the subsets) AUC and sum
of the weights for each model in which all nested subsets of the grid of time
points are allowed in the ensemble. For the survival based methods, CoxBoost
and random forests, we compute the out-of-sample AUC based on pseudo-
observations for comparison to the other methods.

Model R package Binary All Pseudo Nested Pseudo Risk weight AUC weight AUC
weight Used for both Binary and Pseudo Observations General Linear model, with
screening stats 25.2 0.17 0.63 0.15 0.63 0.04 Generalized additive model, with
screening gam 24.7 0.00 0.63 0.19 0.63 0.03 Random forest with bagging
randomForest 22.1 0.40 0.59 0.20 0.60 0.04 Recursive partitioning and
regression trees rpart Inf 0.08 0.57 0.04 0.55 0.14 LASSO glmnet 22.6 0.00
0.64 0.20 0.63 0.03 Multivariate adaptive polynomial spline, with screening
polspline 27.6 0.06 0.61 0.10 0.62 0.21 Extreme gradient boosting xgboost 22.0
0.28 0.63 0.46 0.63 0.25 Used only for Pseudo Observations Kernel support
vector machine kernlab 0.59 $-0.07$ 0.58 0.08 Stepwise selection stats 0.63
0.13 0.63 0.05 Used only for censored survival outcomes Cox Boost CoxBoost
0.63 Random Forests randomForestSRC 0.60

Figure 2: Estimated out of sample ROC curves for the different models used in
the example.  Figure 3: (Left panel) Estimated out-of-sample true and false
positive fractions as functions of the risk percentiles for the cumulative
incidence of surgery at 5 years, accounting for the competing risk of death,
since onset of Crohn’s disease in the validation sample. (Right panel)
Estimated predictiveness curves for the cumulative incidence of surgery at 5
years, accounting for the competing risk of death, since onset of Crohn’s
disease in the validation sample. Points are binned estimates, solid line is
the locally weighted smoothed average of the individual pseudo-observations
(not the binned estimates), and dotted line is the marginal risk.

##  5 Discussion

In our motivating study, we developed a predictive model for abdominal surgery
within 5 years following CD diagnosis based on routinely collected clinical
characteristics observed at the time of diagnosis. The proposed prediction
procedure has an estimated out of sample AUC of 0.64, which is the highest
among the considered methods. Medical management of CD is becoming
increasingly personalized and individualized risk estimations are needed. For
some patients at high risk of surgery, a more aggressive medication or
monitoring strategy may be warranted. In some patients, where a surgical
intervention is inevitable, it should be performed sooner to improve quality
of life.

Although it may be helpful to further understanding, it is not necessary to
’unpack the black box’ for the predictive model to be useful in this setting
because all variables contributing to the predictions will be collected
regardless of their usefulness in predicting risk. For this reason, the model
for the predictions could be used directly, via electronic health records, for
real-time predictions in clinical settings, without the need to know which
variables are driving the models. Implementation of such real-time predictions
and evaluating its impact on patient outcomes is a future area of clinical
research for the authors.

We have shown that the use of pseudo-observations based on the cumulative
incidence can be used with the all libraries for continuous data that are
compatible with the SuperLearner package, and by extension can be used with
any existing or newly developed ML procedures for continuous data in the same
way. We also show that our proposed loss function, the pseudo-observation
based time-dependent AUC, is a useful measure of predictive value for our
setting and that the pseudo-observation based predictiveness curve is a useful
depiction of those predictions. We provide freely available software which
contains this modification to the SuperLearner package in addition to the
other methods described herein.

Although our proposed pseudo-observation based method outperforms the IPCW
binary method in terms of predictive accuracy in the data analysis and when
using multiple times points in the most settings in the simulations, it
performs more similarly when a single time point is used. In addition, in
most, but not all, of the scenarios, our method outperforms the individual
class survival methods. The pseudo-observation based method, however, we found
easier to apply, as the package randomForestSRC caused R to crash repeatedly
during development; we believe this may be due to the large sample size
leading to memory issues in the data analysis. Although our proposed methods
performs similarly to, sometimes better and sometimes worse than, the other
existing methods, a new ML procedure may be developed that will improve
performance of our methods in the future. As creating pseudo-observations is a
general pre-processing step, a pseudo-observation based ensemble could include
this new method without extension, where direct adaptations of the method for
application to censored data will take time and effort.

The idea of stacking pseudo-observations in regression models is not new, and
was an approach suggested for dynamic prediction models of competing risks
(Nicolaie et al., 2013). One advantage of the SuperLearner approach is that it
is possible to include multiple time-point sets, as we did in the data
analysis, to determine what time points from a pre-selected grid are most
useful. The question of how to choose the grid of time points at which to
calculate the pseudo-observations is one that we did not investigate in great
detail. It is of note that we selected the grid somewhat arbitrarily, thus one
might expect that performance can only be improved via some selection or
optimization process. The pre-selection of the grid is an avenue for future
research.

##  6 Bibliography

## References

  * Aalen and Johansen (1978) Aalen, O. O. and Johansen, S. (1978) An empirical transition matrix for non-homogeneous markov chains based on censored observations.  Scandinavian Journal of Statistics, 141–150.
  * Andersen et al. (2003) Andersen, P. K., Klein, J. P. and Rosthøj, S. (2003) Generalised linear models for correlated pseudo-observations, with applications to multi-state models.  Biometrika, 90, 15–27.
  * Andersen and Pohar Perme (2010) Andersen, P. K. and Pohar Perme, M. (2010) Pseudo-observations in survival analysis.  Statistical methods in medical research, 19, 71–99.
  * Binder et al. (2009) Binder, H., Allignol, A., Schumacher, M. and Beyersmann, J. (2009) Boosting for high-dimensional time-to-event data with competing risks.  Bioinformatics, 25, 890–896.
  * Binder et al. (2014) Binder, N., Gerds, T. A. and Andersen, P. K. (2014) Pseudo-observations for competing risks with covariate dependent censoring.  Lifetime data analysis, 20, 303–315.
  * Blanche et al. (2018) Blanche, P., Kattan, M. W. and Gerds, T. A. (2018) The c-index is not proper for the evaluation of $t$-year predicted risks.  Biostatistics.
  * Cortese et al. (2013) Cortese, G., Gerds, T. A. and Andersen, P. K. (2013) Comparing predictions among competing risks models with time-dependent covariates.  Statistics in medicine, 32, 3089–3101.
  * Fine and Gray (1999) Fine, J. P. and Gray, R. J. (1999) A proportional hazards model for the subdistribution of a competing risk.  Journal of the American statistical association, 94, 496–509.
  * Fong et al. (2016) Fong, Y., Yin, S. and Huang, Y. (2016) Combining biomarkers linearly and nonlinearly for classification using the area under the roc curve.  Statistics in medicine, 35, 3792–3809.
  * Gerds and Ozenne (2018) Gerds, T. A. and Ozenne, B. (2018) riskRegression: Risk Regression Models and Prediction Scores for Survival Analysis with Competing Risks.  URL: https://CRAN.R-project.org/package=riskRegression.  R package version 2018.10.03.
  * Gomollón et al. (2016) Gomollón, F., Dignass, A., Annese, V., Tilg, H., Van Assche, G., Lindsay, J. O., Peyrin-Biroulet, L., Cullen, G. J., Daperno, M., Kucharzik, T. et al. (2016) 3rd european evidence-based consensus on the diagnosis and management of crohn’s disease 2016: part 1: diagnosis and medical management.  Journal of Crohn’s and Colitis, 11, 3–25.
  * Graw et al. (2009) Graw, F., Gerds, T. A. and Schumacher, M. (2009) On pseudo-values for regression analysis in competing risks models.  Lifetime Data Analysis, 15, 241–255.
  * Heagerty et al. (2000) Heagerty, P. J., Lumley, T. and Pepe, M. S. (2000) Time-dependent roc curves for censored survival data and a diagnostic marker.  Biometrics, 56, 337–344.
  * Hothorn et al. (2006) Hothorn, T., Bühlmann, P., Dudoit, S., Molinaro, A. and Van Der Laan, M. J. (2006) Survival ensembles.  Biostatistics, 7, 355–373.
  * Ishwaran et al. (2014) Ishwaran, H., Gerds, T. A., Kogalur, U. B., Moore, R. D., Gange, S. J. and Lau, B. M. (2014) Random survival forests for competing risks.  Biostatistics, 15, 757–773.
  * Jacobsen and Martinussen (2016) Jacobsen, M. and Martinussen, T. (2016) A note on the large sample properties of estimators based on generalized linear models for correlated pseudo-observations.  Scandinavian Journal of Statistics, 43, 845–862.
  * Van der Laan et al. (2007) Van der Laan, M. J., Polley, E. C. and Hubbard, A. E. (2007) Super learner.  Statistical applications in genetics and molecular biology, 6.
  * Lakatos et al. (2012) Lakatos, P. L., Golovics, P. A., David, G., Pandur, T., Erdelyi, Z., Horvath, A., Mester, G., Balogh, M., Szipocs, I., Molnar, C. et al. (2012) Has there been a change in the natural history of crohn’s disease? surgical rates and medical management in a population-based inception cohort from western hungary between 1977–2009.  The American journal of gastroenterology, 107, 579.
  * Lin et al. (2011) Lin, H., Zhou, L., Peng, H. and Zhou, X.-H. (2011) Selection and combination of biomarkers using roc method for disease classification and prediction.  Canadian Journal of Statistics, 39, 324–343.
  * Ludvigsson et al. (2016) Ludvigsson, J. F., Almqvist, C., Bonamy, A.-K. E., Ljung, R., Michaëlsson, K., Neovius, M., Stephansson, O. and Ye, W. (2016) Registers of the swedish total population and their use in medical research.  European journal of epidemiology, 31, 125–136.
  * Mogensen and Gerds (2013) Mogensen, U. B. and Gerds, T. A. (2013) A random forest approach for competing risks based on pseudo-values.  Statistics in medicine, 32, 3102–3114.
  * Nicolaie et al. (2013) Nicolaie, M., van Houwelingen, J., de Witte, T. and Putter, H. (2013) Dynamic pseudo-observations: A robust approach to dynamic prediction in competing risks.  Biometrics, 69, 1043–1052.
  * Olén et al. (2017) Olén, O., Askling, J., Sachs, M., Frumento, P., Neovius, M., Smedby, K., Ekbom, A., Malmborg, P. and Ludvigsson, J. F. (2017) Childhood onset inflammatory bowel disease and risk of cancer: a swedish nationwide cohort study 1964-2014.  Bmj, 358, j3951.
  * Overgaard et al. (2017) Overgaard, M., Parner, E. T., Pedersen, J. et al. (2017) Asymptotic theory of generalized estimating equations based on jack-knife pseudo-observations.  The Annals of Statistics, 45, 1988–2015.
  * Pepe et al. (2007) Pepe, M. S., Feng, Z., Huang, Y., Longton, G., Prentice, R., Thompson, I. M. and Zheng, Y. (2007) Integrating the predictiveness of a marker with its performance as a classifier.  American journal of epidemiology, 167, 362–368.
  * Polley et al. (2017) Polley, E., LeDell, E., Kennedy, C. and van der Laan, M. (2017) SuperLearner: Super Learner Prediction.  URL: https://CRAN.R-project.org/package=SuperLearner.  R package version 2.0-22.
  * Ramadas et al. (2010) Ramadas, A., Gunesh, S., Thomas, G., Williams, G. and Hawthorne, A. (2010) Natural history of crohn’s disease in a population-based cohort from cardiff (1986–2003): a study of changes in medical treatment and surgical resection rates.  Gut, 59, 1200–1206.
  * Saha and Heagerty (2010) Saha, P. and Heagerty, P. (2010) Time-dependent predictive accuracy in the presence of competing risks.  Biometrics, 66, 999–1011.
  * Scheike et al. (2008) Scheike, T. H., Zhang, M.-J. and Gerds, T. A. (2008) Predicting cumulative incidence probability by direct binomial regression.  Biometrika, 95, 205–220.
  * Zheng et al. (2012) Zheng, Y., Cai, T., Jin, Y. and Feng, Z. (2012) Evaluating prognostic accuracy of biomarkers under competing risk.  Biometrics, 68, 388–396.


# Markov Chain Concentration with an Application in Reinforcement Learning

Debangshu Banerjee

(February 2021)

###### Abstract

Given $X_{1},\cdot,X_{N}$ random variables whose joint distribution is given
as $\mu$ we will use the Martingale Method to show any Lipshitz Function $f$
over these random variables is subgaussian. The Variance parameter however can
have a simple expression under certain conditions. For example under the
assumption that the random variables follow a Markov Chain and that the
function is Lipschitz under a Weighted Hamming Metric. We shall conclude with
certain well known techniques from concentration of suprema of random
processes with applications in Reinforcement Learning

##  1 Introduction

In the class we have been introduced to the concentration of measure for
functions of independent random variables. A natural extension of this would
be to consider the concentration of measure phenomenon for dependent random
variables. In this report we shall investigate this. We shall see how Markov
Chain concentration form an example in this setting. We shall highlight the
following result which has already been highlighted in a form for independent
random variables within the class lectures

Lipschitz Functions defined on Markov Chains are Sub-Gaussian under certain
Ergodic Properties of the Markov Chain

The way this report is organized is as follows:

  * •

Theoretical Platform and Definitions

  * •

Martingale Method for Concentration of Measure

  * •

Wasserstein Matrix

  * •

General Result and Extensions

We shall also present a hypothesis of how concentration of Markov Chains may
be utilized in the field of Reinforcement Learning. These results are still in
their infancy but we shall nevertheless include them here with hopes of future
work.

##  2 Theoretical Platform and Definitions

We shall assume the following : Let $T$ be a countable index set,
$(\Omega_{t},\mathcal{B}_{t})$ be a measurable space where $\Omega_{t}$ is a
Polish Space for each $t\in T$ and $d_{t}$ be a metric function for each
$\Omega_{t}$.

We shall define the product space $\Omega=\otimes_{t}\Omega_{t}$, the product
sigma algebra $\mathcal{B}=\otimes_{t}\mathcal{B}_{t}$ and define a borel
probability measure $\mu$ on $\Omega$ (Note that that is not a product
measure). We define the metric $d$ on the product space as $d=\sum_{t}d_{t}$

Local Oscillation of $f$ Let $f:\Omega\rightarrow\mathrm{R}$ be a measurable
function. We will define the Local oscillation of $f$ at $i\in T$ as
$\Delta_{i}(f)=\sup_{\begin{subarray}{c}x,y\\\
x^{-i}=y^{-i}\end{subarray}}\frac{f(x)-f(y)}{d_{i}(x_{i},y_{i})}$
That is we keep all except the $i^{th}$ coordinate fixed and vary the $i$
coordinate in forming the ratio.

Note: This definition is siimilar to the bounded difference assumption that
was introduced in the class.
Note: Under this definition $f(x)-f(y)\leq\sum_{i\in
T}\Delta_{i}(f)d_{i}(x_{i},y_{i})$.
Note: Under the assumtption that $f$ is Lipshitz under the weighted Hamming
Distance $|f(x)-f(y)|\leq\sum_{i\in T}c_{i}\mathrm{1}_{x_{i}\neq y_{i}}$, we
have $\Delta_{i}(f)\leq c_{i}$.

Markov Kernel $K:\Omega\times\mathcal{B}\rightarrow[0,1]$ is a Markov Kernel
if $K$ is a Borel Probability measure for each $x\in\Omega$.

Note: If $f$ is any Bounded Measurable function we have $Kf(x)=\int
f(y)K(x,dy)$ a measurable function for almost every $x$ in $\Omega$
Note: If $\mu$ is any borel probability measure we can define $\mu
K(\mathcal{A})=\int K(x,\mathcal{A})d\mu(x)$ a borel probability measure
$\forall\mathcal{A}\in\mathcal{B}$

The following result has been done in class and is presented here without
proof:
Azuma-Hoeffding: Let $\\{\mathcal{F}_{k}\\}_{k\leq n}$ be any filtration and
$M_{k}$ be a martingale difference sequence such that $A_{k}\leq M_{k}\leq
B_{k}$ a.s., then $\sum_{k=1}^{N}M_{k}$ is subgaussian with variance paramter
$1/4\sum_{k=1}^{N}||B_{k}-A_{k}||_{\infty}^{2}$

##  3 Martingale Method for Concentration of Measure

Here we see how to use the Martingale Method introduced in class to get the
concentration of measure phenomenon. The proof will will take the following
course.

  1. 1.

Construct the Martingales
$\mathbf{E}[f(X_{1},\cdots,X_{n}|X_{1},\cdots,X_{i}]-\mathbf{E}[f(X_{1},\cdots,X_{n}|X_{1},\cdots,X_{i-1}]$

  2. 2.

Bound each term above as
$f(x_{1},\cdots,x_{i},X_{i+1},\cdots,X_{N})-f(x_{1},\cdots,X_{i},\cdots,X_{N})$
under the assumption of existence of an Wasserstein Matrix

  3. 3.

Finally use Azuma-Hoeffding to get bounds in term of a matrix norm of this
Wasserstein Matrix

  4. 4.

Motivate how Wassertein Matrices can be represented as Upper Triangular
Matrices and get simple interpretations under stronger assumptions.

Thus let us begin

###  3.1 Construct the Martingale Differences $M_{k}$

Given any borel probability $\mu$ measure on $\Omega$. (Recall $\Omega$ is the
product space) let us define $\\{X_{i}\\}_{i\in T}$ random variables defined
on each $\Omega_{i}$ with joint distribution $\mu$. Let us define the natural
filtration $\\{\mathcal{F}_{k}\\}_{i\in T}$. For each $i\in T$ we definte the
Markov Kernel

$K_{i}(x,dy)=\delta_{x^{[i-1]}}(dy^{[i-1]})\otimes\mu^{[i,n]}(dy^{[i,n]}|x^{[i-1]})$

. That is $K_{i}$ is nothing but the conditional measure of $\mu$ given
$\mu_{[i-1]}$. Note that $[n]=\\{1,\cdots,n\\}$.

Under this definition of the Markov Kernel we have for any bounded measurable
function $f:\Omega\rightarrow\mathrm{R}$,

$K_{i}f(x)=\mathbf{E}_{\mu}[f(X)|X^{[i-1]}=x^{[i-1]}]$

Now let us denote the Martingale Difference under the natural filtration

$M_{i}=K_{i+1}f(x)-K_{i}(f)\\\
=\mathbf{E}_{\mu}[f(X)|X^{[i]}=x^{[i]}]-\mathbf{E}_{\mu}[f(X)|X^{[i-1]}=x^{[i-1]}]$
(1)

Therefore Note that
$f(X_{1},\cdots,X_{N})-\mathbf{E}[f(X_{1},\cdots,X_{N})]=\sum_{i=1}^{N}M_{i}$

###  3.2 Bounding the Martingale differences

Using the tower property of conditional expectations

$M_{i}=\mathbf{E}_{\mu}[f(X)|X^{[i]}=x^{[i]}]-\mathbf{E}_{\mu}[f(X)|X^{[i-1]}=x^{[i-1]}]$
$=\mathbf{E}_{\mu}[f(X)|X^{[i]}=x^{[i]}]-\mathbf{E}_{\mu}[\mathbf{E}_{\mu}[f(X)|X^{[i-1]}=x^{[i-1]},X_{i}]|X^{[i-1]}=x^{[i-1]}]$
$=\int_{\Omega_{[i,n]}}\big{(}\int_{\Omega_{(i,n]}}f(x_{[i-1]}x_{i}y_{(i,n]})\mu^{(i,n]}(dy^{(i,n]}|x^{[i]})-\int_{\Omega_{(i,n]}}f(x_{[i-1]}y_{i}y_{(i,n]})\mu^{(i,n]}(dy^{(i,n]}|x^{[i]},y_{i})\big{)}\\\
\mu^{[i,n]}(dy^{[i,n]}|x^{[i-1]})$ (2)

$=\int_{\Omega_{[i,n]}}\big{(}K_{i+1}f(x_{[i-1]}x_{i}y_{(i,n]})-K_{i+1}f(x_{[i-1]}y_{i}y_{(i,n]})\big{)}\mu^{[i,n]}(dy^{[i,n]}|x^{[i-1]})$
(3)

This implies that each $M_{i}$ can be bounded by

$A_{i}=\int_{X_{[i,n]}}\inf_{x_{i}\in
X_{i}}\big{(}K_{i+1}f(x_{[i-1]}x_{i}y_{(i,n]})-K_{i+1}f(x_{[i-1]}y_{i}y_{(i,n]})\big{)}\mu^{[i,n]}(dy^{[i,n]}|x^{[i-1]})$
$B_{i}=\int_{X_{[i,n]}}\sup_{x_{i}\in
X_{i}}\big{(}K_{i+1}f(x_{[i-1]}x_{i}y_{(i,n]})-K_{i+1}f(x_{[i-1]}y_{i}y_{(i,n]})\big{)}\mu^{[i,n]}(dy^{[i,n]}|x^{[i-1]})$

Now using the bounded differnces definition of $\Delta_{i}(K_{i+1}f)$, we have

$||B_{i}-A_{i}||_{\infty}<||d_{i}||\Delta_{i}(K_{i+1}f)$

where $||d_{i}||=\sup_{x_{i},z_{i}}d_{i}(x_{i},z_{i})$

So now we are in the vicinity of using Azuma-Hoeffding. The only remaining
step remains is how to bound $\Delta_{i}(K_{i+1}f)$ using some known function
of $\Delta_{i}(f)$.

##  4 Wasserstein Matrices

From the theory of general contractive Markov Kernels $K$ with Dobrushin
Coefficient $\theta$ defined as:

$\sup_{x,y}||K(x,.)-K(y,.)||_{TV}<\theta$

, under the Weighted Hamming Distance $\forall f\in Lip(X,d)$

$\Delta_{i}(Kf)<\frac{\theta}{\alpha_{i}}\sum_{j\in T}\alpha_{j}\Delta_{j}(f)$

We refer to [1] for a proof of the above.

As displayed above it is suggested that all $\Delta_{j}(f)$ $\forall j\in T$
influence $\Delta_{i}(Kf)$. With such a relation as a motivation a Wasserstein
Matrix is defined:

Wasserstein MatrixA Markov Kernel $K$ is said to have a Wasserstrien Matrix
$\mathbf{V}=(V_{ij})_{i,j\in T}$ $V_{ij}>0$ if $\forall f\in Lip(X,d)$

$\Delta_{i}(Kf)<\sum_{j\in T}V_{ij}\Delta_{j}(f)\forall i\in T$

or in vector form

$\mathbf{\Delta(Kf)}<\mathbf{V\Delta(f)}$

##  5 General Result and Extensions

###  5.1 SubGaussianity of Markov Kernels under assumption of existence of
Wassertein Matrix

Let us make the following assumption

For each Markov Kernel $K_{i}$, $i\in T$ defined as before assume there exists
Wasserstein Matrices $\mathbf{V^{i}}=(V^{i}_{lm})_{l,m\in T}$ s.t for each
$i\in T$

$\mathbf{\Delta(K_{i}f)}<\mathbf{V^{i}\Delta(f)}$

for each $f\in Lip(X,d)$

What this assumption is essentially saying is a rather strong statement. It
states that we assume that for each markov kernel $K_{i}$ as defined before,
there exists a Wasserstein Matrix $\mathbf{V^{i}}$. We shall soon see how we
can back up our assumption by explicit construction of Wasserstein Matrices
for general probability measures $\mu$, and later give simple expressions of
Wasserstein Matrices under further assumptions about the structure of $\mu$
and the distance metric $d$.
For the time being let us see how this assumption immediately gives us a
SubGaussianity result

As before, from Azuma-Hoeffidfing,
$f(X_{1},....X_{N})-\mathbf{E}f(X_{1},...X_{N})$ is subgaussian with variance
paramter $\sum_{k=1}^{N}||B_{k}-A_{k}||^{2}_{\infty}$ where

$\sum_{k=1}^{N}||B_{k}-A_{k}||^{2}_{\infty}<=\sum_{k=1}^{N}||d_{k}||^{2}\Delta_{k}(K_{k+1}f)^{2}$
$=\sum_{k=1}^{N}(\sum_{j\in T}||d_{k}||V^{k+1}_{kj}\Delta_{j}(f))^{2}$
$=\sum_{k=1}^{N}(\sum_{j\in T}\Gamma_{kj}\Delta_{j}(f))^{2}$
$=\sum_{k=1}^{N}\mathbf{\Gamma\Delta(f)}_{k}^{2}=||\mathbf{\Gamma\Delta(f)}||^{2}_{l_{2}(T)}$

where we define $\Gamma_{ij}=||d_{i}||V^{i+1}_{ij}$

The remaining of the document is sectioned as follows:

  * •

Construction of Wassersterin Matrices Using Couplings

  * •

Wassersterin Matrices under the discrete Metric and Goldstein Coupling

  * •

Contractive Markov Chains

  * •

Uniformly Ergodic Markov Chains

###  5.2 Construction of Wassersterin Matrices Using Couplings

This section illustrates the most general form of a Wasserstein Matrix.

Let $\mathbf{P}^{[i]}_{xy}$ be any coupling of conditional probabilities
$\mu^{(i,n]}(.|x^{[i]})$ and $\mu^{(i,n]}(.|y^{[i]})$ where $x$ and $y$ only
differ in the $i$ coordinate. Then recall that we are trying to relate
$\Delta_{i}(K_{i+1}f)$, we see

$K_{i+1}f(x)-K_{i+1}f(y)=\\\
\int_{X_{(i,n]}}\int_{X_{(i,n]}}(f(x^{[i]},u^{(i,n]})-f(y^{[i]},v^{(i,n]})\mu^{(i,n]}(du^{(i,n]}|x^{[i]})\mu^{(i,n]}(dv^{(i,n]}|y^{[i]})\\\
=\int_{X_{(i,n]}\times
X_{(i,n]}}\mathbf{P}^{[i]}_{xy}(du^{(i,n]},dv^{(i,n]})(f(x^{[i]},u^{(i,n]})-f(y^{[i]},v^{(i,n]})$
(4)

Further note that by our construction of the Markov Kernels $K_{i}$, we have
$\Delta_{i}(K_{j}f)=0$ for $i>j$ Thus we can express the above as

$\leq\Delta_{i}(f)d_{i}(x_{i},y_{i})+\\\
\sum_{j>i}\Delta_{j}(f)\int_{X_{(i,n]}\times
X_{(i,n]}}\mathbf{P}^{[i]}_{xy}(du^{(i,n]},dv^{(i,n]})d_{j}(u_{j},v_{j})$ (5)

or,

$\frac{K_{i+1}f(x)-K_{i+1}f(y)}{d_{i}(x_{i},y_{i})}\leq\Delta_{i}(f)+\sum_{j>i}\frac{\int\mathbf{P}^{[i]}_{xy}d_{j}}{d_{i}(x_{i},y_{i})}\Delta_{j}(f)$

and therefore,

$\Delta_{i}(K_{i+1}f)\leq\Delta_{i}(f)+\sum_{j>i}\sup_{\begin{subarray}{c}x,y\\\
x^{-i}=y^{-i}\end{subarray}}\frac{\int\mathbf{P}^{[i]}_{xy}d_{j}}{d_{i}(x_{i},y_{i})}\Delta_{j}(f)$

Thus from our definition of Wasserstien Matrices

$V^{i+1}_{ij}=\begin{cases}0&\quad\text{if }i>j\\\ 1&\quad\text{if }i=j\\\
\sup_{\begin{subarray}{c}x,y\\\
x^{-i}=y^{-i}\end{subarray}}\frac{\int\mathbf{P}^{[i]}_{xy}d_{j}}{d_{i}(x_{i},y_{i})}&\quad\text{if
}i<j\end{cases}$

###  5.3 Wassersterin Matrices under the discrete Metric and Goldstein
Coupling

Let us introduce a further assumption, that we are working with a discrete
metric.

Under the discrete metric,

$\Gamma_{ij}=\sup_{\begin{subarray}{c}x,y\\\
x^{-i}=y^{-i}\end{subarray}}\mathbf{P}^{[i]}_{xy}[Y^{(0)}_{j}\neq
Y^{(1)}_{j}]$

where
$(Y^{(0)},Y^{(1)})=((Y^{(0)}_{i+1},\cdots,Y^{(0)}_{n}),(Y^{(1)}_{i+1},\cdots,Y^{(1)}_{n})$
is a random variable taking value in $X_{(i,n]}\times X_{(i,n]}$ with
$Y^{(0)}$ distributed as $\mu^{(i,n]}(.|x^{[i]})$ and $Y^{(1)}$ distributed as
$\mu^{(i,n]}(.|y^{[i]})$

$\mathbf{P}^{[i]}_{xy}[Y^{(0)}_{j}\neq
Y^{(1)}_{j}]\leq\mathbf{P}^{[i]}_{xy}[(Y^{(0)}_{j},\cdots,Y^{(0)}_{n})\neq(Y^{(1)}_{j},\cdots,Y^{(1)}_{n})]$

Goldstein Coupling[2]
There exists a coupling $\mathbf{P}^{[i]}_{xy}$ known as the Goldstein Maximal
Coupling s.t.
$\mathbf{P}^{[i]}_{xy}[(Y^{(0)}_{j},\cdots,Y^{(0)}_{n})\neq(Y^{(1)}_{j},\cdots,Y^{(1)}_{n})]=||\mu^{[j,n]}(.|x^{[i]})-\mu^{[j,n]}(.|y^{[i]})||_{TV}$

Thus using the Goldstein Coupling

$\Gamma_{ij}=\begin{cases}0&\quad\text{if }i>j\\\ 1&\quad\text{if }i=j\\\
\sup_{\begin{subarray}{c}x,y\\\
x^{-i}=y^{-i}\end{subarray}}||\mu^{[j,n]}(.|x^{[i]})-\mu^{[j,n]}(.|y^{[i]})||_{TV}&\quad\text{if
}i<j\end{cases}$

###  5.4 Contractive Markov Chains

Let us revisit the example of Contractive Markov Chain

Here we have a directed Markov Model, where $\mu=\mu_{0}K_{1}....K_{n}$ where
$K_{i}$ is a Markov Kernel from $X_{i}$ to $X_{i+1}$ under the Doeblin
Contraction Coefficient

$\sup_{x,y}||K_{i}(x,.)-K_{i}(y,.)||_{TV}<\theta_{i}$

Under this assumption it can be shown using a recursion argument that,

$||\mu^{[j,n]}(.|x^{[i-1]},x_{i})-\mu^{[j,n]}(.|x^{[i]},y_{i})||_{TV}\leq\theta^{j-i}$

where $\theta=\max\theta_{i}$.
We refer to [3] or [4] for a full proof.
Thus from the previous sections, we have

$\Gamma=\begin{pmatrix}1&\theta&\theta^{2}&\dots&\theta^{n-1}\\\
0&1&\theta&\dots&\theta^{n-2}\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\dots&1\end{pmatrix}$

Now from our general result, and noting that if $f$ is Lipshitz under the
weighted Hamming Distance $|f(x)-f(y)|\leq\sum_{i\in
T}c_{i}\mathrm{1}_{x_{i}\neq y_{i}}$, we have $\Delta_{i}(f)\leq c_{i}$. we
have subgaussianity with $||\Gamma||^{2}||c||^{2}$ where $c=(c_{i})_{i\leq n}$

###  5.5 Uniformly Ergodic Markov Chains

Be recalling a general result that uniform ergodic Markov Chains have finite
mixing times, we have the following variant of the general result.

Let us first define what the mixing time for a Markov Chain is:

$\tau(\epsilon)=\min_{t}[\\{\max_{1<i<N-t}\max_{x,y}||\mathbf{P}(X_{i+t}|X_{i}=x)-\mathbf{P}(X_{i+t}|X_{i}=y)||_{TV}\\}\leq\epsilon]$

Let $n=N/\tau(\epsilon)$ and let us form $n$ partitions for the Markov Process
of size $\tau(\epsilon)$. The main theorem that was presented also holds for
such partitions by redefining the $\sigma-algebra$ defined in the proof. Also
we note for this partitioned Markov Chain we have

$||\mu^{[j,n]}(.|\hat{x}^{[i-1]},\hat{x}_{i})-\mu^{[j,n]}(.|\hat{x}^{[i]},\hat{x}_{i})||_{TV}\leq\theta^{j-i-1}$

. We refer to [5] for a full proof.
This gives the matrix $\Gamma$ as

$\Gamma=\begin{pmatrix}1&1&\epsilon&\dots&\epsilon^{n-2}\\\
0&1&1&\dots&\epsilon^{n-3}\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\dots&1\end{pmatrix}$

Therefore for any Lipschitz function w.r.t the weighted discrete metric would
be subgaussian with Variance Parameter $\tau(\epsilon)||c||^{2}||\Gamma||^{2}$

##  6 Reinforcement Learning: A Hypothesis

We conclude this section by perhaps an interesting phenomenon of how such
Markov Chain concentrates can be used in the field of reinforcement Learning.

Recall that for any policy $\pi$ the evolution of states $X_{i}$ and rewards
$r_{i}$ follow a Markov Process.

In reinforcement Learning we are interested in the cumulative reward over a
fixed time which is a function of the trajectory $X_{1},\cdots,X_{N}$.

Let us have a Lipschitz function
$V_{\pi}(X_{1},\cdots,X_{N})\rightarrow\mathbf{R}$. Let us assume the
underlying Markov Process $X_{1},\cdots,X_{N}$ is ergodic under the policy
$\pi$ then we have the rather nice result that the
$V_{\pi}(X_{1},\cdots,X_{N})$ is subgaussian. Let us denote the variance proxy
by $\sigma_{\pi}^{2}$.

In reinforcement Learning we are interested in the
$\mathbf{E}sup_{\pi\in\Pi}V_{\pi}$ where $\Pi$ is a class of policies.

Under a very trivial case if we assume that $\Pi$ contains only finite number
of policies then under our hypotheses of $V_{\pi}$ is subgaussian with
variance parameter $\sigma_{\pi}^{2}$ for any policy $\pi$.

Then using a Maximal Inequality for subgassians

$\mathbf{E}sup_{\pi\in\Pi}V_{\pi}\leq\sqrt{2\sigma^{2}\log|\Pi|}$

where $\sigma=\max\sigma_{\pi}$

This bound is sharp under assumption of independence. We refer to an exercise
problem in [6] for a more general result.

Under the condition $\Pi$ is not finite we have to introduce a sense of
continuity for the map $\pi\rightarrow V_{\pi}$. The simplest such sense would
be to consider a Lipschitz Process

Let $\Pi$ be equipped with a metric $d$ and assume $(V_{\pi})_{\pi\in\Pi}$ is
a Lipschitz Process i.e.

$|V_{\pi}-V_{\pi^{\prime}}|<Cd(\pi,\pi^{\prime})a.s.$

, where $C$ is a random variable Then

Under the assumption of $(V_{\pi})_{\pi\in\Pi}$ is a Lipschitz Process on
$(\Pi,d)$ with $V_{\pi}$ subgaussian with variance proxy $\sigma^{2}$ for
every $\pi\in\Pi$, then
$\mathbf{E}sup_{\pi\in\Pi}V_{\pi}\leq\inf_{\epsilon}\\{\epsilon\mathbf{E}[C]+\sqrt{2\sigma^{2}\log
N(\Pi,d,\epsilon)}\\}$ where $N(\Pi,d,\epsilon)$ is the covering number for
$\epsilon$-nets.

This result is not sharp.

The sharpest known inequality gives a rather interesting result

Dudley’s Inequality and SubGaussian Processes Let $(V_{\pi})_{\pi\in\Pi}$ on
$(\Pi,d)$ be a subgaussian process:

$\mathbf{E}[\exp{\lambda(V_{\pi}-V_{\pi^{\prime}})}]\leq\exp{\lambda^{2}d(\pi,\pi^{\prime})^{2}/2}$

That we indeed have subgaussian processes for Markov Chains can perhaps be
noted from the following observation

Note that under our hypothesis of $V_{\pi}$ subgaussian
$\mathbf{E}[\exp{\lambda V_{\pi}}]<\exp{\lambda^{2}\tau_{\pi}(\epsilon)}$.
Thus it seems probable that our natural metric induced on the class of
policies could be mixing times under various policies.

Dudley’s Inequality For separable subgaussian processes
$\mathbf{E}sup_{\pi\in\Pi}V_{\pi}\leq 12\int_{0}^{\infty}\sqrt{\log
N(\Pi,d,\epsilon)}d\epsilon$

###  6.1 Further hypotheticals

Value Functions

For Value Functions we define

$V_{\pi}(X_{1},\cdots,X_{N})=r(X_{1},\pi(X_{1}))+\cdots+r(X_{N},\pi(X_{N}))$

. Here $r(X_{i},\pi(X_{i}))$ is a random variable and denotes the reward
incurred under the $i^{th}$ stage at state $X_{i}$ under policy $\pi$. Let us
assume that rewards at any stage $i$ is nonnegative and bounded by $r_{i}$.
That is
$0\leq r(X_{i},\pi(X_{i}))\leq r_{i}$ $\forall X_{i}\in\Omega_{i}$ and for any
$\pi$.

Then

$V_{\pi}(x_{1},\cdots,x_{N})-V_{\pi}(y_{1},\cdots,y_{N})=\sum_{i=1}^{N}(r(x_{i},\pi(x_{i}))-r(y_{i},\pi(y_{i})))$
$\leq\sum_{i=1}^{N}r_{i}\mathrm{1}_{\\{x_{i}\neq y_{i}\\}}$

for any policy $\pi$ This gives that the Value function is Lipshitz under the
weighted Hamming Metric and our previous hypothesis would apply verbatim.

Finite State-Action Space

Let State Space be of cardinality $S$, action-space be of cardinality $A$ and
the time horizon be $H$. This would suggest that the number of policies
available would be no more than $SA$. Since this is a finite set, we can apply
the maximal inequality for subgaussians and conclude

Finite State Space Hypothesis Under the assumption that rewards are bounded
between $[0,1]$ $\mathbf{E}sup_{\pi\in\Pi}V_{\pi}\leq\sqrt{H\tau_{mix}\log
SA}$

Proof Sketch: Note that $V_{\pi}$ is subgaussian with variance paramater
$H\tau_{\pi}$. Let where $\tau_{mix}$ be the maximum mixing time of all such
policies. Then we use Maximal Inequality formula.

A trivial upper bound can be established also where instead of finding the
maximum of $\tau_{\pi}$ we just perform an union bound of $SA\tau_{\pi}$
This gives the following bound

Finite State Space Hypothesis Under the assumption that rewards are bounded
between $[0,1]$
$\mathbf{E}sup_{\pi\in\Pi}V_{\pi}\leq\tilde{\mathcal{O}}(\sqrt{HSA\tau_{\pi}})$
for any policy $\pi$ where $\tilde{\mathcal{O}}$ supreses the $\log$ terms

One can see [7] where an algorithm is given which suffers a similar regret.

##  7 Acknowledgements

I would like to express my thanks to Professor Aditya Gopalan and Prof
Himangshu Tyagi. The references are included herein [8], [9]

## References

  * [1] Aryeh Kontorovich and Maxim Raginsky.  Concentration of measure without independence: a unified approach via the martingale method.  In Convexity and Concentration, pages 183–210. Springer, 2017.
  * [2] Katalin Marton et al.  Measure concentration for euclidean distance in the case of dependent random variables.  Annals of probability, 32(3B):2526–2544, 2004.
  * [3] Aryeh Kontorovich.  Obtaining measure concentration from markov contraction.  Markov Processes and Related Fields, 18(4):613–638, 2012.
  * [4] Paul-Marie Samson et al.  Concentration of measure inequalities for markov chains and $\phi$-mixing processes.  The Annals of Probability, 28(1):416–461, 2000.
  * [5] Daniel Paulin et al.  Concentration inequalities for markov chains by marton couplings and spectral methods.  Electronic Journal of Probability, 20, 2015.
  * [6] Ramon van Handel.  Probability in high dimension.  Technical report, PRINCETON UNIV NJ, 2014.
  * [7] Ronald Ortner.  Regret bounds for reinforcement learning via markov chain concentration.  Journal of Artificial Intelligence Research, 67:115–128, 2020.
  * [8] David A Levin and Yuval Peres.  Markov chains and mixing times, volume 107.  American Mathematical Soc., 2017.
  * [9] Maxim Raginsky and Igal Sason.  Concentration of measure inequalities in information theory, communications and coding.  arXiv preprint arXiv:1212.4663, 2012.