# Learning and Concentration for High Dimensional Linear Gaussians: an
Invariant Subspace Approach

Muhammad Abdullah Naeem Author would like to acknowledge, helpful insights and
feedback from Dr. Yuansi ChenMuhammad Abdullah Naeem is with the Department of
Electrical and Computer Engineering, Duke University, Durham, NC 27708, USA,
Email<EMAIL_ADDRESS>
###### Abstract

In this work, we study non-asymptotic bounds on correlation between two time
realizations of stable linear systems with isotropic Gaussian noise.
Consequently, via sampling from a sub-trajectory and using _Talagrands’_
inequality, we show that empirical averages of reward concentrate around
steady state (dynamical system mixes to when closed loop system is stable
under linear feedback policy ) reward , with high-probability. As opposed to
common belief of larger the spectral radius stronger the correlation between
samples, _large discrepancy between algebraic and geometric multiplicity of
system eigenvalues leads to large invariant subspaces related to system-
transition matrix_ ; once the system enters the large invariant subspace it
will travel away from origin for a while before coming close to a unit ball
centered at origin where an isotropic Gaussian noise can with high probability
allow it to escape the current invariant subspace it resides in, leading to
_bottlenecks_ between different invariant subspaces that span
$\mathbb{R}^{n}$, to be precise : system initiated in a large invariant
subspace will be stuck there for a long-time: log-linear in dimension of the
invariant subspace and inversely to log of inverse of magnitude of the
eigenvalue. In the problem of Ordinary Least Squares estimate of system
transition matrix via a single trajectory, this phenomenon is even more
evident if spectrum of transition matrix associated to large invariant
subspace is explosive and small invariant subspaces correspond to stable
eigenvalues. Our analysis provide first interpretable and geometric
explanation into intricacies of learning and concentration for random
dynamical systems on continuous, high dimensional state space; exposing us to
surprises in high dimensions and suggesting, whether it is a problem of system
identification or policy evaluation, practitioner should avoid naive isotropic
Gaussian excitations. Either pre-processing or a different choice of
excitations should be chosen in accordance with the worst case possible size
of invariant subspaces they can encounter, ensuring ‘thorough’ exploration of
state space in minimum number of time steps.

##  I Introduction

Over the last decade, we have seen a tremendous surge in sample complexity
analysis for learning in control tasks. Whether it is the problem of learning
value function corresponding to a control policy (see e.g. [1]) or system
identification as in ([2]), analysis heavily relies on tedious probabilistic
and analytic methods offering very less interpretation or geometric insights.
As a result, uncertainty looms over our current understanding for learning of
dynamical systems via single trajectory, and as we will discuss shortly
afterwards; even a thorough understanding of stable Linear Gaussians(LGs) is
absent. In this paper we conclude for good, sampling complexity, role of
spectral radius and address more recent speculations about large noise being
beneficial [3], [4] and [5],[6]. In fact, it was recently pointed out by [7],
that decay of correlation between two distant samples of trajectory of a
dynamical system is associated to spectral gaps(Functional analytic
phenomenon), not the size of spectral radius of finite dimensional matrix. We
provide a finite-dimensional interpretation of their result and throughout
this paper we will assume working with _high dimensional_ underlying state
space.

Recognizing these ambiguities, we study two simple problems in this paper.
Assume that an unknown LG system is simulated under some stable policy $\pi$
and is assumed to have reached stationarity. In an ideal case, one would like
to know expected reward w.r.t stationary distribution, but happens to only
have access to time averages of reward. How good of an approximation are
temporally-dependent time averages of LG for expected reward compared to i.i.d
draws of reward from stationary distribution. It turns out that answer to this
simple problem, explains all the queries in introduction. Secondly, we study
the problem of Ordinary Least Squares (OLS) estimates for system
identification via single trajectory of dynamical system.

To the best of authors’ knowledge we give first analysis of concentration and
system identification by direct sum decomposition of original state transition
matrix onto its’ invariant subspaces. Along with Talagrand’s concentration
inequality and Gaussian projections on subspaces we are able to conclude, as
opposed to standard beliefs, main issue in system identification and
concentration is invariant subspaces of large dimensions with large magnitude
of associated eigenvalue. If the size of an invariant subspace is large and
excitations are isotropic Gaussians’, with overwhelming probability majority
of the excitation signal will lie in the large invariant subspace(Gaussian
concentration of measure phenomenon). Adding to the complications, if the
spectrum associated to the large invariant subspace is unstable or at a
periphery between being stable or unstable, action of projected state-
transition operator on initial excitation will move the new realization away
from origin while keeping it inside the large invariant subspace. As a result
signal never explores smaller invariant subspaces leading to inaccurate
learned behavior. Large invariant subspaces are a consequence of discrepancy
between algebraic and geometric multiplicity of eigenvalues associated to
state transition matrix. In contrast to our approach, most of the system
identification work focuses on crafting bounds on spectrum as a function of
length of simulated trajectory, or initializing the system at origin and
proving their results on low dimensional simulations. As we will show in
simulation results on OLS in high dimensions, by projecting random initial
excitations into large invariant subspaces, leads to incorrect estimates.

This phenomenon is also apparent when we study the problem of concentration of
time averages around their spatial average(defined by the distribution LG
mixes to, which only happens when spectral radius of system matrix is strictly
less than 1). Even after the associated Markov chain has mixed to its’
stationary distribution, only time-averages of a carefully chosen sub-
trajectory concentrates around spatial average , because once the trajectory
enters a large dimensional invariant subspace it will keep moving away from
the origin (while being inside the aforementioned invariant subspace) till
its’ _first contractive hitting time_ when it is already en route origin and
an isotropic Gaussian excitation will let it hop out to a different invariant
subspace and the trajectory seems to _regenerate_. Although, sampling from a
sub-trajectory of correlated samples give similar concentration results as
i.i.d samples from stationary distribution, but the sub-trajectory is
generated by selecting consecutive realizations, with a gap of : maximum over
all linearly independent invariant subspaces of ‘log-linear in dimension of
the subspace and inversely to log of inverse of magnitude of the associated
eigenvalue ’

The paper is organized as follows. In Section II, we introduce notation and
preliminaries. Section III develops two tensorizartion procedures for
Talagrands’ inequality for general dynamical system which will be at heart of
understanding concentration of stable linear random dynamical systems and
inconsistency of Ordinary Least Squares for explosive systems. In Section IV,
we lay down some facts about stable LGs and propose sampling from sub-
trajectory by leveraging upon Gelfands’ formula. Section V gives a concise
introduction to invariant subspaces associated to a state-transition matrix
and non-asymptotic bounds for concentration via sampling from sub-trajectory
are provided. We begin subsection A of Section VI with concise introduction to
OLS problem via single trajectory and conclude with inconsistency of OLS (for
explosive systems) under isotropic Gaussian exciations via a trivial
application of tensorized Talagrands’ inequality. In subsection B we present
an isoperimetric approach to give an intuition of what may go wrong in high
dimensions which leads to incorrect OLS estimates. Simulation results are
presented in the Section VII, where we show as opposed to existing beliefs on
consistency of OLS for regular systems, OLS in high dimensional regular system
fails when large block of an invariants subspace corrresponds to an explosive
eigenvalue. We conclude with a summary and direction on future work in VIII.

##  II Notation and Preliminaries

####  II-1 Notation

We use ${I}_{n}\in\mathbb{R}^{n\times n}$ to denote the $n$ dimensional
identity matrix. For random variables $x$ and $y$, $Cov(x,y)$ denote the
covariance.
$\mathcal{B}_{\alpha}^{n}:=\\{x\in\mathbb{R}^{n}:\|x\|:=\|x\|_{2}\leq\alpha\\}$
is the $\alpha$-ball in $\mathbb{R}^{n}$. Similarly,
$\mathcal{S}^{p-1}:=\\{x\in\mathbb{R}^{p}:\|x\|_{2}=1\\}$, is the unit sphere
in $\mathbb{R}^{p}$. $\chi_{\\{\\}}()$ is the indicator function, whereas
$\rho(A)$, $\|A\|_{2}$, $\|A\|_{F}$, $det(A)$, $tr(A)$ and $\sigma(A)$
represent the spectral radius, matrix 2-norm , Frobenius norm, determinant,
trace and set of eigenvalues(spectrum) of $A$ respectively. For a positive
definite matrix $A$, largest and smallest eigenvaues are denoted by
$\lambda_{max}(A)$ and $\lambda_{min}(A)$, respectively. Associated with every
rectangular matrix $\mathbb{X}\in\mathbb{R}^{n\times N}$ are its’ singular
values
$\sigma_{1}(\mathbb{X})\geq\sigma_{2}(\mathbb{X}),\ldots,\sigma_{n}(\mathbb{X})\geq
0$, where without loss of generality we assume that $N>n$. Of utmost
importance is largest singular value,
$\sigma_{1}(\mathbb{X}):=\sup_{a\in\mathcal{S}^{N-1}}\|\mathbb{X}a\|$ and the
least singular value
$\sigma_{n}(\mathbb{X}):=\inf_{a\in\mathcal{S}^{N-1}}\|\mathbb{X}a\|$.
Condition number of a matrix $\mathbb{X}$ is the ratio of the largest and
least singular value, denoted by
$\kappa(\mathbb{X})=\frac{\sigma_{1}(\mathbb{X})}{\sigma_{n}(\mathbb{X})}$. If
the span of image space of $\mathbb{X}$ is $\mathbb{R}^{n}$, more compactly
written as $Im(\mathbb{X})=\mathbb{R}^{n}$, least singular value equals the
inverse of the norm of inverse of matrix $\mathbb{X}$ i.e.,
$\sigma_{n}(\mathbb{X})=\frac{1}{\|\mathbb{X}^{-1}\|}$. A function
$g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p}$ is Lipschitz with constant $L$ if
for every $x,y\in\mathbb{R}^{n}$, $\|g(x)-g(y)\|\leq L\|x-y\|$.

A sequence ${\\{a(N)\\}_{N\in\mathbb{N}}\in\mathcal{O}(N)}$, if it increases
at most linearly in $N$ (this is not limited to asymptotic results).
$\mathcal{O}(1)$ will be used to denote quantities independent of the size of
the underlying state space or number of the iterations. Space of probability
measure on $\mathcal{X}$(continuous space) is denoted by
$\mathcal{P(\mathcal{X})}$ and space of its Borel subsets is represented by
$\mathbb{B}\big{(}\mathcal{P(\mathcal{X})}\big{)}$. For a function $r$ and
$\mu\in\mathcal{P(X)}$, we use $<r>_{\mu}$ to denote expectation of $r$ w.r.t
$\mu$. Finally, for a set $\mathcal{K}\subseteq\\{1,...,M\\}$, its complement
is $\mathcal{K}^{\complement}:=\\{1,...,M\\}\setminus\mathcal{K}$.

On a metric space $(\mathcal{X},d)$, for $\mu,\nu\in\mathcal{P(\mathcal{X})}$,
we define Wasserstein metric of order $p\in[1,\infty)$ as

$\mathcal{W}_{d}^{p}(\nu,\mu)=\bigg{(}\inf_{(X,Y)\in\Gamma(\nu,\mu)}\mathbb{E}~{}d^{p}(X,Y)\bigg{)}^{\frac{1}{p}};$
(1)

here, $\Gamma(\nu,\mu)\in P(\mathcal{X}^{2})$, and $(X,Y)\in\Gamma(\nu,\mu)$
implies that random variables $(X,Y)$ follow some probability distributions on
$P(\mathcal{X}^{2})$ with marginals $\nu$ and $\mu$. Another way of comparing
two probability distributions on $\mathcal{X}$ is via relative entropy, which
is defined as

$Ent(v||u)=\left\\{\begin{array}[]{lr}\int\log\bigg{(}\frac{d\nu}{d\mu}\bigg{)}d\nu,&\text{if}~{}\nu<<\mu,\\\
+\infty,&\text{otherwise}.\end{array}\right.$ (2)

Before we introduce the mathematical framework to derive concentration for
dependent random variables, we introduce the following results utilized later
in this work.

##### Talagrands’ inequality or Transport-Entropy Inequality

Consider metric space $(\mathcal{X},d)$ and reference probability measure
$\mu\in P(\mathcal{X})$. Then we say that $\mu$ satisfies
$\mathcal{T}_{1}^{d}(C)$ or to be concise $\mu\in\mathcal{T}_{1}^{d}(C)$ for
some $C>0$ if for all $\nu\in P(\mathcal{X})$ it holds that

$\mathcal{W}_{d}(\mu,\nu)\leq\sqrt{2CEnt(\nu||\mu)}.$ (3)

######  Theorem 1 ([8]).

$\mu$ satisfies $\mathcal{T}_{1}^{d}(C)$ if and only if for any Lipschitz
function $f$ with $<f>_{\mu}:=\mathbb{E}_{\mu}f$, it holds that

$\displaystyle\int
e^{\lambda(f-<f>_{\mu})}d\mu\leq\exp(\frac{\lambda^{2}}{2}C\|f\|_{L(d)}^{2}),\hskip
15.0pt$ (4) $\displaystyle\text{where~{}~{}~{}~{}}\|f\|_{L(d)}:=\sup_{x\neq
y}\frac{|f(x)-f(y)|}{d(x,y)}.$

######  Remark 2.

(4) along with the Markov inequality implies that if we sample $x$ from
$\mu\in\mathcal{T}_{1}^{d}(C)$, then

$\mathbb{P}\Bigg{[}\bigg{|}r(x)-<r>_{\mu}\bigg{|}>\epsilon\Bigg{]}\leq
2\exp\bigg{(}-\frac{\epsilon^{2}}{2C\|r\|_{L(d)}^{2}}\bigg{)}.$ (5)

##  III Extending Concentration to Dependent Random Variables via
Tensorization

Under the action of some state dependent policy $\pi$, we consider a closed-
loop random dynamical system of the form

$x_{k+1}=F\big{(}x_{k},\pi(x_{k}),\epsilon_{k}\big{)},\hskip
10.0pt\text{with~{}}\epsilon_{k}\hskip 10.0pti.i.d,$ (6)

where $x_{k}\in\mathbb{R}^{n}$ for all $k\in\mathbb{N}$ and
$F:\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\longrightarrow\mathbb{R}^{n}$.
For the problem of concentration of ergodic averages, we will assume that the
transition kernel converges to some stationary distribution $\mu_{\infty}$
under Wasserstein metric $W_{d}$ equipped with some distance function $d$.
This random dynamical system can be viewed as a Markov chain
$x^{N}:=(x_{i})_{i=1}^{N}$ with distribution $\mu^{N}\in P$
$(\mathcal{X}^{N})$ and
$P^{m}(x,\mathcal{B}):=\mathbb{P}(x_{m}\in\mathcal{B}|x_{0}=x)$, for all Borel
subsets $\mathcal{B}$ of $\mathcal{X}$. We can extend the metric $d$ to
$\mathcal{X}^{N}$ as

$d_{(N)}(x^{N},y^{N})=\sum_{i=1}^{N}d(x_{i},y_{i}).$ (7)

If $\mu^{N}\in\mathcal{T}_{1}^{d_{(N)}}\big{(}\mathcal{O}(N)\big{)}$ and $r$
is one Lipschitz, i.e., $\|r\|_{L(d)}\leq 1$, then
$\Phi(x^{N}):=\frac{1}{N}\sum_{i=1}^{N}r(x_{i})$ satisfies
$\|\Phi\|_{L(d_{(N)})}\leq\frac{1}{N}$; plugging these results into (4), we
obtain that

$\displaystyle\mu^{N}\Bigg{[}\bigg{|}\frac{1}{N}\sum_{i=1}^{N}r(x_{i})-\mathbb{E}\Bigg{(}\frac{1}{N}\sum_{i=1}^{N}r(x_{i})\Bigg{)}\bigg{|}>\epsilon\Bigg{]}$
(8) $\displaystyle\leq 2\exp\bigg{(}-\frac{N\epsilon^{2}}{2C}\bigg{)}.$

###  III-A Contractivity and Uniform Transport Constants

As one would wonder from (7), when does the T-E for process level law of
Markov chain, increases at worse linearly with dimension (in sample term)?
Sufficient conditions (see e.g., [9, 10]) are

$\displaystyle(i)$ $\displaystyle\hskip 10.0ptP(x,\cdot)\in
T_{1}^{d}(C),\hskip 5.0pt\text{for all}~{}{x\in\mathcal{X}},\text{and
some}~{}{C>0},$ (9) $\displaystyle(ii)$ $\displaystyle\hskip
10.0pt\mathcal{W}_{d}(P(x,\cdot),P(y,\cdot))\leq\hat{\lambda}d(x,y),\hskip
5.0pt\text{for all }~{}{(x,y)\in\mathcal{X}^{2}}$ (10)

and some $\hat{\lambda}\in[0,1)$.

Property (9) is often referred to as existence of a uniform transportation
constant and (10) represents contractivity of the Markov Chain in the
Wasserstein metric / _spectral gap in the Wasserstein sense_. Now, the
following result holds.

######  Theorem 3.

If (9) and (10) hold, process level distribution of samples from a Markov
chain $(x_{1},\ldots,x_{N})$, which we will denote as
$Law(x_{1},\ldots,x_{N})$, denoted by $\mathcal{\mu}^{N}$ satisfies
$T_{1}^{d_{(N)}}\bigg{(}\frac{CN}{(1-\hat{\lambda})^{2}}\bigg{)}$, for all
$N\in\mathbb{N}$.

###### Proof.

See Theorem 2.5 of [9] for a detailed proof. ∎

##### Decay of correlation.

By combining conditions from (9) and (10), with Taylor’s expansion for small
$\lambda$ (terms of order up to $\lambda^{2}$) appearing on both sides in
Bobkov-Gotze dual form (4), for all $x\in\mathcal{X}$ it holds that

$\displaystyle|Cov_{P_{x}}[f(x_{n}),f(x_{n+k})]|\leq\frac{\hat{\lambda}^{k}}{1-\hat{\lambda}^{2}}C\|f\|_{L(d)}^{2}.$
(11)

Another tensorization procedure that we will heavily rely on in Section VI, to
better understand deviation inequalities for singular values of data matrix,
is an an $\ell^{2}$ inspired metric on $\mathcal{X}^{N}$ as in [11], i.e.,

$d_{(N)}^{2}(x^{N},y^{N}):=\sqrt{\sum\nolimits_{i=1}^{N}d^{2}(x_{i},y_{i})}.$
(12)

##  IV STABLE DYNAMICAL SYSTEMS AND ALMOST INDEPENDENT SUBTRAJECTORY

###  IV-A Independent sampling from invariant measure

Markov chain under consideration is $n$ dimensional LG with isotropic noise:

$x_{t+1}=Ax_{t}+w_{t},\hskip 10.0pt\rho(A)<1\hskip 10.0pt\text{and i.i.d
}~{}w_{t}\thicksim\mathcal{N}(0,\mathcal{I}_{n}).$ (13)

It mixes to stationary distribution
$\mu_{\infty}\thicksim\mathcal{N}(0,P_{\infty})$, where the controllability
grammian $P_{\infty}$ is the unique positive definite solution of the
following Lyapunov equation:

$A^{T}P_{\infty}A-P_{\infty}+I_{n}=0.$ (14)

$\displaystyle\mathbb{P}\Big{[}\big{|}\frac{1}{N}\sum_{i=1}^{N}r(x_{i})-\mu_{\infty}(r)\big{|}>\epsilon\Big{]}\leq
2\exp\bigg{(}-\frac{N\epsilon^{2}}{2\lambda_{max}(P_{\infty})}\bigg{)}.$ (15)

###  IV-B Sampling from a Sub-trajectory of LGs

Stability in controls community for Linear systems correspond to $\rho(A)<1$
(marginally stable corresponds $\rho(A)\geq 1$ and explosive system when
$\rho(A)>1$) and it is an established result in real analysis (Gelfands
formula) that for all $\rho\in(\rho(A),1)$ and $k\in\mathbb{N}$, there exists
a finite positive constant $L_{\rho}$ such that $\|A^{k}\|\leq
L_{\rho}\rho^{k}$. Naive Wasserstien contractivity condition (10) does not
hold as $L_{\rho}\rho$ can be very large. In order to extend preceding result
to spectral radius case, we define _first contractive hitting time_ as:

$\hat{k}:=\min\\{k\in\mathbb{N}:\|A^{k}\|<1\\},$ (16)

that is smallest natural number such that $\hat{k}$-th step transition kernel
$P(x_{\hat{k}}=\cdot|x_{0})$ contractive in Wasserstein sense.

So, instead of considering the original trajectory
$(x_{0},x_{1},x_{2},\ldots)$, we will consider the sub-trajectory
$(x_{0},x_{\hat{k}},x_{2\hat{k}},\ldots)$ with modified LG dynamics:

$\displaystyle x_{\hat{k}(i+1)}=A^{\hat{k}}x_{\hat{k}(i)}+s_{\hat{k}(i)}$ (17)

where $s_{\hat{k}(i)}\thicksim\mathcal{N}\big{(}0,\Sigma_{\hat{k}}\big{)}$,
i.i.d with $\Sigma_{\hat{k}}:=\sum_{l=0}^{\hat{k}-1}[A^{l}][A^{l}]^{T}$
$\forall i\in\mathbb{N}$. Notice that $\Sigma_{\hat{k}}$ is positive definite:
Consequently,

$\displaystyle(i)$ $\displaystyle\hskip 10.0ptP^{\hat{k}}(x,\cdot)\in
T_{1}^{d}\bigg{(}\|\Sigma_{\hat{k}}^{\frac{1}{2}}\|^{2}\bigg{)},\hskip
65.0pt\text{for all}~{}{x\in\mathcal{X}}.$ (18) $\displaystyle(ii)$
$\displaystyle\hskip
10.0pt\mathcal{W}_{d}(P^{\hat{k}}(x,\cdot),P^{\hat{k}}(y,\cdot))\leq\lambda_{\hat{k}}d(x,y),$
(19)

for all ${(x,y)\in\mathcal{X}^{2}}$ and
${\lambda_{\hat{k}}:=\|A^{\hat{k}}\|\in[0,1)}$.
$\mu_{\hat{k}}^{N}:=Law(x_{\hat{k}(1)},\ldots,x_{\hat{k}(N)})\in
T_{1}^{d_{(N)}}\bigg{(}\frac{\|\Sigma_{\hat{k}}^{\frac{1}{2}}\|^{2}N}{(1-\|A^{\hat{k}}\|)^{2}}\bigg{)}$
and from 3 if we start $x_{0}\thicksim\mu_{\infty}$, we have the following
concentration:

$\displaystyle\mathbb{P}\Big{[}\big{|}\frac{1}{N}\sum_{i=1}^{N}r(x_{\hat{k}(i)})-\mu_{\infty}(r)\big{|}>\epsilon\Big{]}\leq
2\exp\bigg{(}-\frac{N\epsilon^{2}[1-\|A^{\hat{k}}\|]^{2}}{2\lambda_{\max}(P_{\infty})}\bigg{)}$
(20)
$\displaystyle=2\exp\bigg{(}-\frac{N\epsilon^{2}[1-\|A^{\hat{k}}\|]^{2}\lambda_{\min}(P_{\infty}^{-1})}{2}\bigg{)}.$
(21)

Therefore, compared to i.i.d samples from $\mu_{\infty}$, temporally dependent
although identically distributed samples(each individually distributed as
$\mu_{\infty}$) spaced $\hat{k}$ times apart can give us sharp concentration
but we need a trajectory of length $N\hat{k}$. Notice that as
$\hat{k}\rightarrow\infty$ stationary chain concentration becomes i.i.d
concentration. For a detail analysis of non-stationary case we refer to [12].
In order to bound first contractive hitting time we need to understand:

##  V Structure of invariant sub spaces associated to eigenvalue problem of a
non-symmetric operator

Position or magnitude of eigenvalues associated to a linear operator $A$ only
provides partial information about its’ properties (for the ease of
exposition, throughout this paper we will assume that $A$ does not have any
non-trivial null space). In fact knowing $A$ is equivalent to knowing its’
invariant subspaces (see e.g., [13]). Roughly speaking, algebraic multiplicity
of eigenvalues follow from determinant of the matirx.

$det(zI-A)=\prod_{i=1}^{K}(z-\lambda_{i})^{m_{i}},$ (22)

where $\lambda_{i}$ are distinct with multiplicity $m_{i}$. Complication
happens when $dim[N(A-\lambda_{i}I)]<m_{i}$, which leads to invariants
subspace (spanned by more that one linearly independent vector). Consequently,
states space can be written as direct sum decomposition of $A-$ invariant
subspaces.

$\mathbb{R}^{n}=M_{\phi(1)}\oplus M_{\phi(2)}\oplus\ldots\oplus M_{\phi(L)}$
(23)

and respective orthogonal projetcions $[E_{\phi(i)}]_{i=1}^{L}$ such that
identity matrix can be written as:

$I_{n}=E_{\phi(1)}\oplus E_{\phi(2)}\oplus\ldots\oplus E_{\phi(L)}$ (24)

where $\phi$ is a surjective map from $\\{1,\ldots,L\\}$ to $\sigma(A)$.
$\phi$ is bijective iff eigenvectors span $\mathbb{R}^{n}$. _In the case of
gap between between algebraic and geometric multiplicity related to some
element of $\sigma(A)$_. Consider the invariant subspace $M_{\lambda}$, for
some $\lambda\in\sigma(A)$, with algebraic multiplicity of $\lambda$ is
$|B_{\lambda}|$ but only one linearly independent eigenvector $v_{1}$ such
that $Av_{1}=\lambda v_{1}$. So we generate generalized eigenvector
$v_{2},v_{3},\ldots,v_{|B_{\lambda}|}$ recursively as $(A-\lambda
I)v_{2}=v_{1}$ and $(A-\lambda I)v_{3}=v_{2}$ and so on. We have the following
$k$ -th step iteration:

$\displaystyle A^{k}v_{1}=\lambda^{k}v_{1}$ $\displaystyle
A^{k}v_{2}=\lambda^{k}v_{2}+\binom{k}{1}\lambda^{k-1}v_{1}$ $\displaystyle
A^{k}v_{3}=\lambda^{k}v_{3}+\binom{k}{1}\lambda^{k-1}v_{2}+\binom{k}{2}\lambda^{k-2}v_{1}$
$\displaystyle\ldots=\ldots$ $\displaystyle
A^{k}v_{|B_{\lambda}|}=\lambda^{k}v_{|B_{\lambda}|}+\binom{k}{1}\lambda^{k-1}v_{|B_{\lambda}|-1}+\ldots$
$\displaystyle+\binom{k}{|B_{\lambda}|-1}v_{1}$ (25)

######  Theorem 4.

Although evident from the preceding iterations, we can rigorously upper and
lower bound norm of the $k-th$ iteration associated to action of matrix $A$ on
invariant subspace $M_{\lambda}$, precisely given as:

$\displaystyle|\lambda|^{k}\sum_{m=0}^{|B_{\lambda}|-1}\frac{1}{|\lambda|^{m}}\leq\|A^{k}_{M_{\lambda}}\|_{2}\leq|\lambda|^{k}k^{|B_{\lambda}|}\sum_{m=0}^{|B_{\lambda}|-1}\frac{1}{|\lambda|^{m}},$
(26)

where $A^{k}_{M_{\lambda}}:=A^{k}E_{\lambda}$

###### Proof.

The result follows by a simple variation of bounds provided in [14].

$\displaystyle|\lambda|^{k}\sum_{m=0}^{|B_{\lambda}|-1}\frac{1}{|\lambda|^{m}}\leq\|A^{k}_{M_{\lambda}}\|_{2}$
$\displaystyle=\sum_{m=0}^{|B_{\lambda}|-1}\binom{k}{m}|\lambda|^{k-m}$ (27)
$\displaystyle\leq|\lambda|^{k}k^{|B_{\lambda}|}\sum_{m=0}^{|B_{\lambda}|-1}\frac{1}{|\lambda|^{m}}$

In the first and second inequality we have used the fact
$1\leq\binom{k}{m}\leq k^{|B_{\lambda}|}$ for $m\in[0,\ldots,B_{|\lambda|}-1]$

∎

######  Remark 5.

Moreover, if $|\lambda|\in(0,1)$ then:

$\|A^{k}_{M_{\lambda}}\|_{2}\leq
k^{|B_{\lambda}|}|B_{\lambda}||\lambda|^{k+1-|B_{\lambda}|}$ (28)

and

$k\geq\frac{\ln(|B_{\lambda}|)}{\ln(\frac{1}{|\lambda|})}+\frac{|B_{\lambda}|\ln(k)}{\ln(\frac{1}{|\lambda|})}+(|B_{\lambda}|-1)$
(29)

suffices for $\|A^{k}_{M_{\lambda}}\|_{2}<1$.

Now we are in a position to give a conclusive analytic remark on the how to
pick a sub-trajectory to get sharp concentration for time averages around
their spatial average as raised in subsection IV-B.

######  Theorem 6.

First contractive hitting time for operator $A$ restricted to invariant
subspace $M_{\phi(i)}$ is
$\mathcal{O}\bigg{(}\frac{|B_{\phi(i)}|\ln|B_{\phi(i)}|}{\ln(\frac{1}{|\phi(i)|})}\bigg{)}$.
Therefore, contractive hitting time for linear operator $A$ is the worst
contractive hitting time over all invariant subspaces.

$\hat{k}=\min\big{[}k\in\mathbb{N}:k\geq\max_{i\in
1,\ldots,L}\bigg{(}\frac{4|B_{\phi(i)}|\ln|B_{\phi(i)}|}{\ln\frac{1}{|\lambda_{\phi(i)}|}}\bigg{)}\big{]},$
(30)

which is also verified via simulations shown in Fig 1.

Figure 1: First contractive hitting time of $n\times n$ Jordan forms. case 1:
single Jordan block with eigenvalue of 0.86. case 2: single Jordan block with
eigenvalue of 0.9. case 3: $(n-1)\times(n-1)$ Jordan block with eigenvalue 0.9
and a single block with eigenvalue 0.9.

######  Remark 7.

Since a very lose upper-bound on first contractive hitting time is
$\bigg{(}\frac{n\ln n}{\ln\frac{1}{\rho(A)}}\bigg{)}$, most of the literature
on system identification for stable system uses this as an overhead for
choosing block sizes (see e.g., [3] and [2]).

As $\|A^{\hat{k}}\|$ is contractive in Wasserstein sense, we have the
following exponential convergence of the sub-trajectory (17) to stationary
distribution:

######  Theorem 8.

$\lambda_{\hat{k}}:=\|A^{\hat{k}}\|_{2}<1$ and for all $m\in\mathbb{N}$

$\displaystyle W_{d}(P_{x}^{\hat{k}(m)},\mu_{\infty})$
$\displaystyle\leq\lambda_{\hat{k}}^{m}W_{d}(P_{x}^{\hat{k}},\mu_{\infty})$
$\displaystyle\leq\lambda_{\hat{k}}^{m}\sqrt{\lambda_{\hat{k}}\|x\|+Tr\big{(}[\sqrt{\Sigma_{\hat{k}}}-\sqrt{P_{\infty}}]^{2}\big{)}}.$
(31)

###### Proof.

It suffices to show that for all $x,y$ in $\mathbb{R}^{n}$,
$W_{d}(P_{x}^{\hat{k}},P_{y}^{\hat{k}})\leq\lambda_{\hat{k}}\|x-y\|$, see
e.g., [15]. Given any $x\in\mathbb{R}^{n}$, we can write it as a direct sum
$x=\sum_{i=1}^{L}E_{\phi(i)}x$. Leveraging on orthogonality:
$E_{\phi(i)}E_{\phi(j)}=0$ for $i\neq j$, we also have
$\|x\|_{2}=\sum_{i=1}^{L}\|E_{\phi(i)}x\|_{2}$. Now,
$\|A^{\hat{k}}x\|_{2}=\|A^{\hat{k}}\sum_{i=1}^{L}E_{\phi(i)}x\|_{2}\leq\sum_{i=1}^{L}\|A_{M_{\phi(i)}}^{\hat{k}}\|_{2}\|E_{\phi(i)}x\|_{2}\leq\lambda_{\hat{k}}\|x\|_{2}$,
where the last inequality follows from hypothesis and first contractive
hitting time for individual block given in (29). Therefore
$W_{d}(P_{x}^{\hat{k}},P_{y}^{\hat{k}})\leq\lambda_{\hat{k}}\|x-y\|$ for all
$x,y$ in $\mathbb{R}^{n}$ and the result follows by realizing general
expression for Wasserstein distance between two Gaussians in terms of their
mean and covariance see e.g., [16]. ∎

##  VI SYSTEM IDENTIFICATION VIA SINGLE TRAJECTORY

###  VI-A Ordinary Least Square

In this section we analyse the problem of OLS estimation for system transition
matrix $A$ from single observed (as in [17], [3], [2]) trajectory of
$(x_{0},x_{1},\ldots,x_{N})$ satisyfing:

$x_{t+1}=Ax_{t}+\eta_{t},\hskip 10.0pt\text{ where
}\eta_{t}\thicksim\mathcal{N}(0,\mathcal{I}_{n}).$ (32)

OLS solution is:

$\hat{A}=\arg\min_{B\in\mathbb{R}^{n\times
n}}\sum_{t=0}^{N-1}\|x_{t+1}-Bx_{t}\|.$ (33)

Let $\mathbb{X}_{+}=[x_{1},x_{2},\ldots,x_{N}]$ and
$\mathbb{X}_{-}=[x_{0},x_{1},\ldots,x_{(N-1)}]$, and noise covariates
$E=[\eta_{0},\eta_{1},\ldots,\eta_{N-1}]$ then the closed form expression for
Least squares solution and error are:

$\displaystyle\hat{A}=\mathbb{X}_{+}\mathbb{X}_{-}^{T}(\mathbb{X}_{-}\mathbb{X}_{-}^{T})^{-1}$
(34)
$\displaystyle\|A-\hat{A}\|=\|E\mathbb{X}_{-}^{T}(\mathbb{X}_{-}\mathbb{X}_{-}^{T})^{-1}\|$
(35)

Error can be upper bounded:

$\displaystyle\|A-\hat{A}\|=\|E\mathbb{X}_{-}^{T}(\mathbb{X}_{-}\mathbb{X}_{-}^{T})^{-1}\|$
$\displaystyle\leq\|E\|\sigma_{1}(\mathbb{X}_{-})\frac{1}{\sigma_{n}^{2}(\mathbb{X}_{-})}$
$\displaystyle=\frac{\sigma_{1}(E)\kappa(\mathbb{X}_{-})}{\sigma_{n}(\mathbb{X}_{-})},$
(36)

where, recall $\kappa(\mathbb{X}_{-})$ is the condition number of
$\mathbb{X}_{-}$. It is a well known result in Random Matrix Theory(see e.g.,
[18]) if $n$ and $N$ are increased while maintaining their ratio
$\frac{n}{N}=\gamma\in(0,1)$, then:

$\displaystyle\sigma_{1}(E)\thicksim\sqrt{N}+\sqrt{n},$ (37)

where $\thicksim$ here denotes typical behavior explained in discussion below
Theorem 12. So we are left with task of bounding singular values of data
matrix $\mathbb{X}_{-}$ (which contains dependent random variables, but here
we will see Talagrands’ inequality in all of its’ glory )

######  Theorem 9.

We have the following conentration bounds on all singular values
$[\sigma_{k}(\mathbb{X}_{-})]_{k=1}^{n}$ of the data matrix

  1. 1.

if $\|A\|_{2}<1$:

$\displaystyle\mathbb{P}\Big{[}\big{|}\sigma_{k}(\mathbb{X}_{-})-\mathbb{E}\sigma_{k}(\mathbb{X}_{-})\big{|}>\epsilon\Big{]}$
$\displaystyle\leq 2\exp\bigg{(}-\epsilon^{2}[1-\|A\|]^{2}\bigg{)}$

  2. 2.

if $\|A\|_{2}=1$ :

$\displaystyle\mathbb{P}\Big{[}\big{|}\sigma_{k}(\mathbb{X}_{-})-\mathbb{E}\sigma_{k}(\mathbb{X}_{-})\big{|}>\epsilon\Big{]}$
$\displaystyle\leq 2\exp\bigg{(}-\frac{\epsilon^{2}[e-1]}{N(N+1)}\bigg{)}$

  3. 3.

$\|A\|_{2}>1$ :

$\displaystyle\mathbb{P}\Big{[}\big{|}\sigma_{k}(\mathbb{X}_{-})-\mathbb{E}\sigma_{k}(\mathbb{X}_{-})\big{|}>\epsilon\Big{]}$
$\displaystyle\leq
2\exp\bigg{(}-\frac{\epsilon^{2}[N-1]}{\|A\|^{N}e(N+1)}\bigg{)}$

###### Proof.

The idea of the proof follows from tensorization of Talagrands’ inequality for
dependent covariates. Notice that for any $a\in S^{N-1}$:

$(x_{0},x_{1},\ldots,x_{N-1})\mapsto\mathbb{X}_{-}a\mapsto\|\mathbb{X}_{-}a\|_{2}$
(38)

is a 1-Lipschitz map from $(\mathbb{R}^{n})^{N}$ to $\mathbb{R}$ under metric
$d_{(N)}^{2}$. Since $S^{N-1}$ is compact, $\inf$ and $\sup$ are attained,
least and largest singular values are 1-Lipschitz (so are intermediate
singular values via their min-max characterization/ Courant-Fischer theorem).
Let $\mu^{N}:=Law(x_{0},x_{1},\ldots,x_{N-1})$

  * •

if $\|A\|_{2}<1$:

$\mu^{N}\in T_{1}^{d_{(N)}^{2}}\bigg{(}\frac{1}{[1-\|A\|_{2}]^{2}}\bigg{)}$
(39)

  * •

if $\|A\|_{2}=1$:

$\mu^{N}\in T_{1}^{d_{(N)}^{2}}\bigg{(}\frac{N(N+1)}{e-1}\bigg{)}$ (40)

  * •

if $\|A\|_{2}>1$:

$\mu^{N}\in T_{1}^{d_{(N)}^{2}}\bigg{(}\frac{\|A\|^{N}e(N+1)}{N-1}\bigg{)}$
(41)

see Proposition 4.1 in [19] and the result follows. ∎

######  Remark 10.

Since the case $\|A\|_{2}>1$ includes explosive systems and variance seems to
deteriorate with number of samples; therefore, least singular value does not
concentrate and OLS is inconsistent. Also notice that how tensorization of a
sub-trajectory for stable systems give better concentration estimates as
$\mu_{\hat{k}}^{N}:=Law(x_{\hat{k}(1)},\ldots,x_{\hat{k}(N)})\in
T_{1}^{d_{(N)}^{2}}\bigg{(}\frac{\|\Sigma_{\hat{k}}^{\frac{1}{2}}\|^{2}}{(1-\|A^{\hat{k}}\|)^{2}}\bigg{)}$
and consequently, dimension free concentration inequalities for its singular
values.

Now we provide an intuitive explanation of why OLS is inconsistent for
explosive systems via isoperimetric reasonings.

###  VI-B Projection of Isortopic Gaussian

######  Remark 11.

This subsection is only for instructive purposes, employing notations like
approximately, which we do not justify as the results here are isoperimetric
in nature which we plan on considering for future work but even now
intuitively explains what may go wrong in learning for high dimensional
dynamical systems

######  Theorem 12.

Let $\gamma_{n}$ be isotropic Gaussian in $\mathbb{R}^{n}$ and
$S\subset\mathbb{R}^{n}$ be a $k-$ dimensional subspace. Given
$x\in\mathbb{R}^{n}$, let $x_{S}$ denote the projection of $x$ onto $S$. Then
for any $\delta\in(0,1)$

$\displaystyle\gamma_{n}\bigg{(}x\in\mathbb{R}^{n}:\frac{\|x_{S}\|}{\|x\|}\geq(1-\delta)^{-1}\sqrt{\frac{k}{n}}\bigg{)}\leq
e^{-\frac{\delta^{2}k}{4}}+e^{-\frac{\delta^{2}n}{4}}$
$\displaystyle\gamma_{n}\bigg{(}x\in\mathbb{R}^{n}:\frac{\|x_{S}\|}{\|x\|}\leq(1-\delta)\sqrt{\frac{k}{n}}\bigg{)}\leq
e^{-\frac{\delta^{2}k}{4}}+e^{-\frac{\delta^{2}n}{4}},$ (42)

see Lemma 3.2 in [20].

That is ratio of the norm of projection onto a $k-$ dimensional subspace is
_typically_ ($\thicksim$) $\sqrt{\frac{k}{n}}$ and a remarkable advantage of
this observation is that we can apply it on $A$\- invariant sub-spaces to get
an intuitive understanding of bottlenecks between invariant sub-spaces.
Formally speaking if $E_{\lambda}$ is a projection onto $A-$invariant subspace
$M_{\lambda}$ than:

$\frac{\|E_{\lambda}\eta_{0}\|}{\|\eta_{0}\|}\thicksim\sqrt{\frac{|B_{\lambda}|}{n}}$
(43)

and as the size of underlying state space increases and so does the size of
invariant subspace while maintaining $n-|B_{\lambda}|=\mathcal{O}(1)$ we can
almost certainly conclude that almost all the excitation signal $\eta_{0}$
lies inside $M_{\lambda}$ (follows from (42))

######  Remark 13.

Estimate in (42) is at the heart of Johnson-Lindenstrauss (JL) lemma which is
a very powerful tool for dimensionality reduction in high dimensional
euclidean space to circumvent curse of dimensionality, see e.g., Lemma 2.1 of
chapter 1 in [21]

Conditioned on $x_{0}=0$, we can express $N-th$ realization of the signal as:

$x_{N}=\sum_{t=1}^{N}A^{N-t}\eta_{t-1}.$ (44)

For the explosive case i.e., $\lambda\in\sigma(A)$ such that $|\lambda|>1$, if
$\eta_{0}\thicksim\mathcal{N}(0,I_{n})$, with overwhelming probability we have
the following norm bound

$|\lambda|^{k-|B_{\lambda}|+1}\bigg{(}\frac{|B_{\lambda}|^{\frac{3}{2}}}{n^{\frac{1}{2}}}\bigg{)}\lesssim\frac{\|A^{k}E_{\lambda}\eta_{0}\|}{\|\eta_{0}\|}\lesssim
k^{|B_{\lambda}|}|\lambda|^{k}\bigg{(}\frac{|B_{\lambda}|^{\frac{3}{2}}}{n^{\frac{1}{2}}}\bigg{)}$
(45)

Since $A^{k}E_{\lambda}=A^{k}_{M_{\lambda}}E_{\lambda}$, signal
$A^{k}E_{\lambda}\eta_{0}$ lies entirely inside subspace $M_{\lambda}$ and
(45) can be interpreted as: if the block size is huge and corresponding
eigenvalue is explosive, predominantly the realizations of our dynamical
system will lie inside $M_{\lambda}$ and consequently any realistic algorithm
would fail to learn the behavior of the system in other parts of the state
space.

Assume that $\mathbb{R}^{n}=M_{\lambda}\oplus M_{\lambda}^{\perp}$, with
$\sigma(A_{M_{\lambda}^{\perp}})$ being stable with largest eigenvalue
$\lambda^{\perp}\in(0,1)$. Since $n-|B_{\lambda}|=\mathcal{O}(1)$, implies
that $|B_{\lambda}^{\perp}|=\mathcal{O}(1)$ it trivially follows that for even
small values of $N$

$\displaystyle\frac{\|A^{N}E_{\lambda}^{\perp}\eta_{0}\|}{\|A^{N}E_{\lambda}\eta_{0}\|}\approx
0.$ (46)

##  VII Simulation Result on OLS

Figure 2: Ordinary Least Squares on regular explosive systems with isotropic
excitations

$A\in\mathbb{R}^{50\times 50}$ with only two distinct eigenvalues
$\lambda_{1}=1.5$ with Block size $|B_{\lambda_{1}}|=47$ and
$\lambda_{2}=-0.5$ with Block size of $|B_{\lambda_{2}}|=3$, using direct sum
decomposition:

$A=A_{\lambda_{1}}\oplus A_{\lambda_{2}}.$ (47)

Similarly the state space can be represented as direct sum decomposition of
two-A-invariant subspaces:

$\mathbb{R}^{n}=M_{\lambda_{1}}\oplus M_{\lambda_{2}}.$ (48)

We can get span of two invariant subsapces related to the two Jordan blocks by
computing $[M,D]=eig(A)$, where $M=[M_{\lambda_{1}}\hskip
3.0ptM_{\lambda_{2}}]$ and orthogonal projections onto $M_{\lambda_{1}}$ and
$M_{\lambda_{2}}$ by applying projection operator
$E_{1}=M_{\lambda_{1}}(M_{\lambda_{1}}^{T}M_{\lambda_{1}})^{\dagger}M_{\lambda_{1}}^{T}$
and
$E_{2}=M_{\lambda_{2}}(M_{\lambda_{2}}^{T}M_{\lambda_{2}})^{\dagger}M_{\lambda_{2}}^{T}$,
respectively and $\dagger$ is used to represent pseudo-inverse. In Fig 2 we
perform four different simulations and record error in operator norm as the
length of trajectory increase from $50$ to $100$. Curves in blue correspond to
initial condition being an orthogonal projection defined by $E_{1}$ on
randomly sampled isotropic Gaussian and similarly curves in red correspond to
initial condition being an orthogonal projection defined by $E_{2}$ on
randomly sampled isotropic Gaussian. As suggested by Talagrands-inequality for
variance of least singular value of explosive systems in Theorem (9), OLS is
inconsistent.

######  Remark 14.

System considered in the case study is completely regular (as defined in [22])
but still least square estimates are incorrect (Gaussian concentration of
measure phenomenon is at work)

##  VIII Conclusion and Future Work

In this paper, we began with the study of the correlation between two distinct
time realization of stable linear systems in high dimensions with excitations
of isotropic Gaussian. We employ a novel approach, where rather than just
basing our analysis on the magnitude of eigenvalue, we also took into
consideration the _geometric content_ related to the operator via information
on the size of its invariant sub-spaces w.r.t state transition matrix. Which
not only provided us with a geometric insight but also improved concentration
results via sampling from sub-trajectory with smaller gaps (compared to what
was previously believed ) between two almost uncorrelated samples. Leveraging
on these geometric insights along with Talagrands’ inequality in the later
half of the paper we analyse inconsistency issues with OLS for explosive
systems with isotropic Gaussian as excitation signal. It turns out that
problem of system identification for high dimensional dynamical systems is
inseparable from concentration of measure phenomenon and isoperimetric
inequalities in high dimensions, where we were able to show issues in OLS with
isotropic excitations by employing results associated to concentration of
projections of Gaussian measures on large sub-spaces. In the future we intend
on formalizing these results via tools in high dimensional geometry.

## References

  * [1] S. Tu and B. Recht, “Least-squares temporal difference learning for the linear quadratic regulator,” in _International Conference on Machine Learning_ , 2018, pp. 5005–5014.
  * [2] A. Tsiamis and G. J. Pappas, “Linear systems can be hard to learn,” in _2021 60th IEEE Conference on Decision and Control (CDC)_. IEEE, 2021, pp. 2903–2910.
  * [3] M. Simchowitz, H. Mania, S. Tu, M. I. Jordan, and B. Recht, “Learning without mixing: Towards a sharp analysis of linear system identification,” in _Conference On Learning Theory_ , 2018, pp. 439–473.
  * [4] S. Oymak and N. Ozay, “Revisiting ho–kalman-based system identification: Robustness and finite-sample analysis,” _IEEE Transactions on Automatic Control_ , vol. 67, no. 4, pp. 1914–1928, 2021.
  * [5] A. Tsiamis and G. J. Pappas, “Online learning of the kalman filter with logarithmic regret,” _IEEE Transactions on Automatic Control_ , 2022.
  * [6] A. Tsiamis, I. Ziemann, N. Matni, and G. J. Pappas, “Statistical learning theory for control: A finite sample perspective,” _arXiv preprint arXiv:2209.05423_ , 2022.
  * [7] M. A. Naeem, “Concentration phenomenon for random dynamical systems: An operator theoretic approach,” _arXiv preprint arXiv:2212.03670_ , 2022.
  * [8] S. G. Bobkov and F. Götze, “Exponential integrability and transportation cost related to logarithmic sobolev inequalities,” _Journal of Functional Analysis_ , vol. 163, no. 1, pp. 1–28, 1999.
  * [9] H. Djellout, A. Guillin, L. Wu _et al._ , “Transportation cost-information inequalities and applications to random dynamical systems and diffusions,” _Annals of Probability_ , vol. 32, no. 3B, pp. 2702–2732, 2004.
  * [10] F. Bolley and C. Villani, “Weighted csiszár-kullback-pinsker inequalities and applications to transportation inequalities,” in _Annales de la Faculté des sciences de Toulouse: Mathématiques_ , vol. 14, no. 3, 2005, pp. 331–352.
  * [11] F. Malrieu, “Logarithmic sobolev inequalities for some nonlinear pde’s,” _Stochastic processes and their applications_ , vol. 95, no. 1, pp. 109–132, 2001.
  * [12] M. A. Naeem and M. Pajic, “Transportation-inequalities, lyapunov stability and sampling for dynamical systems on continuous state space,” _arXiv preprint arXiv:2205.12448_ , 2022.
  * [13] T. Kato, _Perturbation theory for linear operators_. Springer Science & Business Media, 2013, vol. 132.
  * [14] J. A. Tropp, “An elementary proof of the spectral radius formula for matrices,” 2001.
  * [15] M. Hairer and J. C. Mattingly, “Yet another look at harris’ ergodic theorem for markov chains,” in _Seminar on Stochastic Analysis, Random Fields and Applications VI_. Springer, 2011, pp. 109–117.
  * [16] C. R. Givens and R. M. Shortt, “A class of wasserstein metrics for probability distributions.” _Michigan Mathematical Journal_ , vol. 31, no. 2, pp. 231–240, 1984.
  * [17] T. Sarkar and A. Rakhlin, “Near optimal finite time identification of arbitrary linear dynamical systems,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 5610–5618.
  * [18] M. Rudelson and R. Vershynin, “Smallest singular value of a random rectangular matrix,” _Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences_ , vol. 62, no. 12, pp. 1707–1739, 2009.
  * [19] G. Blower and F. Bolley, “Concentration inequalities on product spaces with applications to markov processes,” _arXiv preprint math/0505536_ , 2005.
  * [20] A. Barvinok, “Math 710: Measure concentration,” _Lecture notes_ , 2005.
  * [21] L. Ambrosio, “Optimal transportation, geometry and functional inequalities,” _(No Title)_.
  * [22] M. K. S. Faradonbeh, A. Tewari, and G. Michailidis, “Finite time identification in unstable linear systems,” _Automatica_ , vol. 96, pp. 342–353, 2018.


# Assessing survival models by interval testing

Ben Lee

###### Abstract

When considering many survival models, decisions become more challenging in
health economic evaluation. In this paper, we present a set of methods to
assist with selecting the most appropriate survival models. The methods
highlight areas of particularly poor fit. Furthermore, plots and overall
p-values provide guidance on whether a survival model should be rejected or
not.

Keywords: Survival analysis, calibration, health economics, decision making

An R-package for these methods, including examples is available here.

##  1 Introduction

The primary results of a clinical trial are often right-censored survival
data. To understand the efficacy of a treatment, survival models must be
fitted to this data. A choice between fitted survival models must then be
made, and this is often the primary driving factor in the economic evaluation
of a new medical product. During this process, a large number of models are
usually considered. However, with a growing emphasis on flexible survival
mdoels[1] and a requirement for standard parametric models[2], the number of
models considered can easily exceed 20, making decisions challenging. Our
proposed methods can assist this process by indicating which survival models
provide implausible estimates.

When we plot these survival models, we will refer to them as survival curves.
Recommendations on selecting the most appropriate survival models are provided
in NICE TSD 14[2]. These include visual curve inspection, AIC/BIC tests,
clinical validity, and the use of external data. Except for AIC/BIC, these
approaches are subjective and often time-consuming to implement, especially
when considering many curves. It would be useful if curves were quickly shown
to be implausible. This would reduce the number of feasible curves, making the
decision process easier. One approach that may facilitate this is calibration,
but this is rarely used in health technology assessment.

In survival analysis, calibration refers to the agreement between predicted
probabilities and observed events within a specific time interval. The most
common calibration methods are extensions of the Hosmer-Lemeshow goodness-of-
fit test[3],[4]. These include the Grønnesby and Borgan test (specific to the
Cox proportional hazard model)[5], the Nam and D’Agostino test (applicable to
more general models)[6], as well as smoothed calibration curves and their
associated metrics[7]. In all of these approaches, patients are assigned a
risk score based on the fitted model and are grouped (often into deciles)
based on these scores. The groups are then evaluated to determine how well-
calibrated a survival model is. Our approach is different from this because
patients do not need to be assigned risk scores. This is especially useful for
models where the risk scores are equal among all patients, which is the case
for our worked examples.

The paper is structured as follows: methods, application, and discussion. The
methods section describes: a new calibration approach with intervals based on
censoring, an extension to this where intervals can be specified, and finally
methods for generating an overall p-value. The application section considers
two examples. The first example evaluates a simple exponential model,
providing an introduction to the developed methods. The second example is more
practical and considers survival data which is challenging to model. In this
example, seven parametric survival models are fitted. Subsequently, the full
range of developed methods evaluate which curves are most appropriate.
Finally, before the discussion is given, the results of a simulation
experiment are summarized.

##  2 Methods

###  2.1 Interval testing

We consider right-censored survival data, where for each subject we have time
and event given by $(t,E)\in(\mathbb{R}_{>0},\\{1,0\\})$. We split this data
into two types:

  * •

The observed data (event=1) with times $t_{i,1}$ such that $0<t_{i,1}\leq
t_{i+1,1}$ and $i=1,\ldots,E.$

  * •

The censored data (event=0) with times $t_{i,0}$ such that $0<t_{i,0}\leq
t_{i+1,0}$ and $i=1,\ldots,C.$

Subsequently, we define the set of unique censor times $t_{j,0}^{\prime}$ for
$j=1,\ldots,J$. That is, the values $t_{j,0}^{\prime}$ are the unique values
from the set of $t_{i,0}$ values, and thereby the $t_{j,0}^{\prime}$ values
satisfy $0<t_{j,0}^{\prime}<t_{j+1,0}^{\prime}$. We now define the intervals
of sequential censor times such that for $i=1,...,J-1$

$\begin{split}I_{0}&=(0,\ t_{1,0}^{\prime}]\\\ I_{i}&=(t_{i,0}^{\prime},\
t_{i+1,0}^{\prime}]\\\ I_{J}&=(t_{J,0}^{\prime},\ \infty).\end{split}$

We are interested in the hypothesis that the event generating process is
described by the model $M(\boldsymbol{\beta})$, with probability density
function $f_{M}$ and fixed parameters $\boldsymbol{\beta}$.

Thereby, imagine we have not seen the data. However, we still have the
intervals $I_{j}$ and the model $M(\boldsymbol{\beta})$. We now make the
assumption that the intervals themselves, $I_{j}$, are uninformative of the
event generating process. It follows that each patient that enters the
interval $I_{j}$:

  * •

either experiences an event or does not experience an event in interval
$I_{j}$.

  * •

has the same probability of experiencing an event as any other patient
entering $I_{j}$.

  * •

has a probability of experiencing an event defined by the event generating
process $M(\boldsymbol{\beta})$ (when assuming the null hypothesis holds).

It follows that under the null, the following process occurs within each
interval:

$events(I_{j})\sim binomial(N_{I_{j}},p_{I_{j}}),$ (1)

where $events(I_{j})$ are the number of events observed within the interval
$I_{j}$.

$p_{I_{j}}=$“the probability of a subject experiencing an event in interval
$I_{j}$ given the subject makes it to interval $I_{j}$”, which is given by

$\begin{split}p_{I_{j}}&=\mathbb{P}(T\in I_{j}\ |\ T>min(I_{j})\ )\\\
&=\mathbb{P}(T\in I_{j})\ /\ \mathbb{P}(T>min(I_{j}))\\\
&=\int_{I_{j}}f_{M}(t)dt\Big{/}\int_{min(I_{j})}^{\infty}f_{M}(t)dt.\end{split}$
(2)

Finally, $N_{I_{j}}=$“the number of subjects at the start of interval
$I_{j}$.” That is, the number of starting patients $(C+E)$, minus the patients
lost to events in the previous $j-1$ intervals, minus the patients lost to
censors in the previous $j-1$ intervals. That is,

$N_{I_{j}}=\ C+E\ -\sum_{k=0}^{j-1}events(I_{k})\
-\sum_{k=0}^{j-1}censors(I_{k}),$ (3)

where we define the number of censors observed within the interval $I_{j}$ as
$censors(I_{j})$. Note that this quantity usually takes a value of 1.
(However, there are exceptions to this. For example, the final interval $J$
has $censors(I_{J})=0$ and intervals with multiple censors of equal time at
the upper boundary of the interval have $censors(I_{j})>1$.)

For each interval $I_{j}$, we have specified a binomial model with known
$N_{I_{j}}$, $p_{I_{j}}$, and observed $events(I_{j})$. Each interval $I_{j}$
thereby has a midpoint p-value for observing $events(I_{j})$ in that interval.

Extension to specified intervals

The methodology above describes how events can be modelled within each
interval. However, the intervals themselves are entirely determined by the
censor times. It is often of interest to understand events in specified
intervals (for example 10 evenly-spaced intervals up to the last censor.)
Below, we extend the methodology to allow for such an analysis.

Define the intervals of interest for $k=1,\dots,K$ as

$\begin{split}V_{1}&=(S_{0},\ S_{1}]\\\ ...\\\ V_{k}&=(S_{k-1},\
S_{k}].\end{split}$

We have unique censor times $t_{j,0}^{\prime}$ for $j=1,\ldots,J$. Define the
set of unique times $\tau_{i}$ as the unique set of the union of the
$t_{j,0}^{\prime}$ values and the $S_{k}$ values, with $\tau_{i}<\tau_{i+1}$
for $i=1,...,L_{\tau}-1$, where $L_{\tau}$ denotes the number of unique $\tau$
values. We now construct the intervals

$\begin{split}I_{0}^{\prime}&=(0,\tau_{1}]\\\ ...\\\
I_{i}^{\prime}&=(\tau_{i},\tau_{i+1}]\\\ ...\\\
I_{L_{\tau}}^{\prime}&=(\tau_{L_{\tau}},\infty].\end{split}$

As in Equation 1, we now have a known distribution for each interval
$I_{i}^{\prime}$. The observations in said intervals are described by

$events(I_{i}^{\prime})\sim binomial(N_{I_{i}^{\prime}},p_{I_{i}^{\prime}}),$

with $p_{I_{i}^{\prime}}$ and $N_{I_{i}^{\prime}}$ defined analogously to
Equations 2 and 3, with the sum in Equation 3 now being over the $I^{\prime}$
intervals.

We now consider the $I_{i}^{\prime}$ intervals as the subintervals that make
up each $V_{k}$ interval. Define the collection of subintervals
$l=1,...,L_{k}$ as

$\\{I_{l,V_{k}}^{\prime}\\}=\\{I_{i}^{\prime}\ \ \ such\ that\ \
I_{i}^{\prime}\subseteq V_{k}\\}.$

It follows that for each interval $V_{k}$ its subintervals follow

$\begin{split}events(I_{1,V_{k}}^{\prime})&\sim
binomial(N_{I_{1,V_{k}}^{\prime}},p_{I_{1,V_{k}}^{\prime}})\\\ ...\\\
events(I_{L_{k},V_{k}}^{\prime})&\sim
binomial(N_{I_{L_{k},V_{k}}^{\prime}},p_{I_{L_{k},V_{k}}^{\prime}}).\\\
\end{split}$

Thereby $\sum^{L_{k}}_{l=1}events({I_{l,V_{k}}^{\prime}})$ is distributed as
the realization of a sum of binomials with known parameters. That is, the
number of observed events in each interval $V_{k}$ follows a known
distribution. We can thus obtain a midpoint p-value for each specified
interval $V_{k}$.

###  2.2 Test statistics

The primary recommendation of this paper is the use of interval plots, such as
Figure 1b. To support these plots, a single overall p-value can be a useful
tool. Such a p-value can be obtained by condensing the plot’s constituent
p-values into a single overall p-value.

When deriving these test statistics, we assume that the midpoint p-value for
each interval follows a continuous uniform[0,1] distribution under H0. When it
comes to the examples, this assumption does not necessarily hold. This is
justified in detail in Appendix A and C.

We derive two test statistics for obtaining an overall p-value. From Birnbaum
(1954)[8], we know that selecting a best test statistic for H0 is not trivial,
and a more modern discussion is given in [9]. Our first test statistic is
continuous and allows for extreme p-values to have a high influence on the
overall p-value. The second test statistic is discrete and aims to mitigate
the influence of extreme p-values.

####  2.2.1 Transformed Fisher test (TFT) statistic

To motivate the test, we consider a p-value for a given interval. A high
p-value indicates more events have occurred in the interval than would
generally be expected under H0. A low p-value indicates fewer events than
expected. It follows that both high and low p-values indicate poor model fit
within a given interval. As such, we will construct a test statistic that can
detect an uncommon amount of high and low p-values.

For each interval, we have a midpoint p-value, $p_{j}$. Under H0, we make the
assumption these p-values are $uniform[0,1]$. The p-values are transformed to
the form

$U_{j}=\begin{cases}2p_{j}\ \ \ \ \ \ \ \ \ \ \ for\ \ p_{j}\leq 0.5\\\
2(1-p_{j})\ \ \ \ for\ \ p_{j}>0.5\end{cases}$

It follows that $U_{j}\sim uniform[0,1]$. By the classical result from Fisher,
we obtain the test statistic

$T_{cont}=-2\sum_{j=1}^{I}log(U_{j})\sim\chi^{2}_{{}_{2I}},$

where we reject for large $T_{cont}$ and $I$ is the total number of intervals.
The intuition for this test is that unfavourable p-values (close to 0 or 1)
are transformed such that they correspond to $U_{j}$ values close to 0. This
makes the overall test able to detect unfavourable p-values.

####  2.2.2 Protection Against Very Small Intervals (PAVSI) test statistic

The transformed Fisher test is useful as it enables us to leverage the full
range of reported p-values. However, in some situations, this may not be
ideal. For example, a single p-value that is very close to 0 can heavily
influence the test. We argue that a test that is not dominated by extreme
p-values is a useful tool.

For each interval, we have a midpoint p-value, $p_{j}$. Under H0, we make the
assumption these p-values are $uniform[0,1]$. The p-values are transformed to
the form

$V_{j}=\begin{cases}1\ \ \ for\ \ p_{j}\geq 0.975\ \ or\ \ p_{j}\leq 0.025\\\
0\ \ \ else\end{cases}$

It follows that under H0 we have $V_{j}\sim bernoulli(0.05)$, and thereby

$T_{pavsi}=\sum_{j=1}^{I}V_{j}\sim binom(I,0.05),$

where we reject for large $T_{pavsi}$ values. It is noted that $T_{pavsi}$ is
a discrete test statistic, and we recommend reporting the midpoint p-value for
this, see Appendix A.

Compared to the transformed Fisher test, the PAVSI is more pragmatic in terms
of extreme p-values. However, for PAVSI, a p-value of say 0.026 would
contribute equivalently to a p-value of 0.5, which is arguably a limitation of
this approach.

##  3 Application

The methods above described two general approaches to interval testing:

  * •

Use the censor times to specify the intervals

  * •

Specify the intervals yourself (“extension to specified intervals”)

In example 1, we focus on the first approach only. In the second example, both
approaches are described. Note that the two test statistics (2.2) can be used
in either approach.

For our examples, we use real world overall survival data for melanoma as is
described in [10]. Example 1, uses the ‘Dabrafenib’ treatment arm of the
BREAK-3 trial[11]. Example 2, uses the ‘Dabrafenib+Trametinib’ treatment arm
of the COMBI-d trial[12]. All data and code is available on Github here.

###  3.1 Example 1: A fitted exponential model

In this example, we fit the simplest survival model (exponential) to our data.
Following this, we want to determine how appropriate the exponential model is.
First, we do a visual check of the curve fit, Figure 1a. Visually, the Kaplan-
Meier curve and the exponential curve seem to be misaligned, suggesting poor
fit. However, an interval test provides a more rigorous assessment of this,
Figure 1b. For this example, we are only interested in an interval test where
the intervals are selected using the censor times.

To perform the test, we obtain the values $I_{j}$, $N_{I_{j}}$, and
$events(I_{j})$ directly from the data. From Equation 1, it remains to obtain
$p_{I_{j}}$ for the fitted exponential model. As the probability density
function for the exponential distribution is given by $f_{M}(t)=\lambda
e^{-\lambda t}$, using Equation 2, it follows that

$p_{I_{j}}=1-exp(\lambda[min(I_{j})-max(I_{j})]\ ).$

As the value $\lambda$ is known (from fitting the exponential model to the
data), it follows that we have fully specified binomial realisations for each
interval $I_{j}$. This information for the first five intervals is presented
in Table 1.

$I_{j}$ | N.risk | $p_{I_{j}}$ | $events(I_{j})$ |  $\mathbb{E}$(events) | p-mid | p-rand
---|---|---|---|---|---|---
(0, 1] | 187 | 0.026 | 1 | 4.8 | 0.027 | 0.034
(1, 1.6] | 185 | 0.016 | 1 | 3 | 0.127 | 0.111
(1.6, 2.5] | 183 | 0.024 | 1 | 4.4 | 0.037 | 0.059
(2.5, 3.4] | 181 | 0.023 | 4 | 4.2 | 0.491 | 0.538
(3.4, 4.2] | 176 | 0.021 | 5 | 3.6 | 0.771 | 0.762
… | … | … | … | … | … | …
Table 1: Summary of interval data for Example 1

The p-value for each interval is discrete as it comes from a binomial model.
It follows that each discrete p-value does not follow a continuous
uniform[0,1] distribution under H0. This can be rectified by using randomized
p-values[13]. However, in the interest of reproducibility and clarity, we will
focus our discussion on the results of the midpoint p-values. (This is
discussed in more detail in Appendix A.)

Figure 1: Part a: The KM curve for the data in purple, censors are given by
dashed vertical black lines, and the fitted exponential curve is the overlaid
curve in black. Part b: the binomial interval test for this data. Intervals
are defined by the censors, where for each interval, red indicates a
Bonferroni rejection, and grey indicates an individual flag. Part c: Numbers
at risk for the data

For each interval, we consider the midpoint p-value with respect to a 2-sided
test. That is, if $p_{I_{j}}\leq 0.025$ or $p_{I_{j}}\geq 0.975$, we have
reason to suspect the interval generating said value. We call this the
individual test.

Of greater concern, are the intervals that fail the Bonferroni test[14]. The
Bonferroni test is chosen to strongly control the familywise error rate at a
level of 0.05. In particular, we reject if $p_{I_{j}}\leq 0.025/I$ or
$p_{I_{j}}\geq 0.975/I$, where $I=$“the total number of intervals”. This means
that under H0 (that the model holds for all intervals) a single Bonferroni
rejection of an interval is only expected to occur for 1 in 20 times the data
is generated under H0.

In Figure 1, we present the results for the Bonferroni and individual tests
alongside the Kaplan-Meier (KM) curve, the fitted exponential model, and the
numbers at risk. A single Bonferroni rejection of an interval is good
reasoning for rejecting a model. In Figure 1, the extremely poor fit
(Bonferroni rejection) is at the end of the curve. This is very unfavourable
in health economic evaluation, where plausible extrapolations beyond the
observed data are paramount to model selection. That is, the Bonferroni
rejection at the end of the curve indicates that the extrapolation will be
poor, which is more detrimental compared with a Bonferroni rejection at the
start of the curve.

An overall p-value for the fitted exponential model is also of interest. For
the transformed Fisher test statistic, there are $I=43$ intervals, and our
test statistic was computed as $-2\sum_{j=1}^{I}log(U_{j})=81.84$, giving a
p-value of 0.607. For the PAVSI test statistic, there are $I=43$ intervals.
From Figure 1, we can read off that $T_{pavsi}=4$ (as there are 3 individual
rejections and 1 Bonferroni rejection). This gives us an overall midpoint
p-value of 0.114. (These test statistics should be used with caution. In
Appendix C, the test statistics are shown to be highly conservative when using
the censor interval approach.)

As we can see, the two test statistics yield very different results. It is our
recommendation that a model should not be accepted or rejected based entirely
on the outcome of such tests. Instead, they are tools that provide an initial
overview of model fit. The greatest detail and insight can be obtained by
assessing the binomial interval plot alongside specific interval p-values and
an overall p-value. For this particular example, we would recommend that the
model is rejected. This is primarily due to the Bonferroni rejection at the
end of the curve.

###  3.2 Example 2: Selecting a parametric survival model for challenging data

In this example, we fit seven standard parametric survival models. The
survival data has a challenging shape, making it difficult for the models to
provide good visual fit. Our aim is to use the interval test methods to
determine whether any of the survival models are appropriate or whether
flexible survival models are required.

Intervals selected based on censors

First, we perform the interval test with intervals chosen using censors as in
Example 1. The 7 parametric models were fitted using flexsurv [15]. The final
plots and tables were generated by the BITsurv package which is available
here.

The test statistics for each parametric model are presented in Table 2. The
PAVSI and TFT overall p-values are well above the 0.05 threshold for all of
the survival models. This might lead us to believe that the models are
providing good fit. However, as discussed in Appendix C, these test statistics
are highly conservative when using the censor interval approach, meaning the
overall p-values are not informative unless they indicate a rejection (less
than 0.05.)

For each model, there are between 1 and 4 individual flags. However, these
flags are primarily descriptive. Of greater interest are the Bonferroni
rejections, of which there are none. This seems good, and from Table 2 alone
we might believe that our curves provide good fit. However, this is an
unfortunate result of using intervals based on censors, see poor fit in Figure
2. Due to the large number of censors, the majority of intervals are very
small, which makes it difficult to detect poor fit within these intervals. For
example, there are many small intervals for the last 5 months of the curve. It
would be more informative to consider a single interval for those 5 months.
Because of these limitations, a specified interval approach is required.

Figure 2: Part a: The KM curve (purple), censors (vertical dashed lines), and the fitted generalized gamma curve (black). Part b: The binomial interval test results for 10 evenly-spaced intervals. Part c: Numbers at risk. Parametric survival model | Bonferroni rejections | Individual flags | PAVSI | TFT
---|---|---|---|---
Exponential | 0 | 3 | 0.468 | 0.854
Gamma | 0 | 2 | 0.686 | 0.887
Generalised gamma | 0 | 1 | 0.881 | 0.955
Gompertz | 0 | 4 | 0.267 | 0.884
Log-logistic | 0 | 1 | 0.881 | 0.977
Log-normal | 0 | 2 | 0.696 | 0.968
Weibull | 0 | 3 | 0.468 | 0.871
Table 2: Example 2: Interval test results using censor defined intervals.
(Note that Bonferroni rejected intervals are also counted as individual
flags.)

Specified intervals (10 evenly-spaced)

There are many ways these intervals could be chosen, and several approaches
are briefly mentioned in the discussion. For the sake of brevity, we will
outline one of the simpler approaches. The data is split into 10 evenly-spaced
intervals (owing to the classical convention of the Hosmer-Lemeshow test.)
These intervals start at time 0 and end at the largest censor time. That is,
for $k=1,...,10$

$V_{k}=\frac{1}{10}\cdot t.max\cdot(k-1,\ k].$

To simplify code implementation, p-values were calculated using the sinib R
package [16], which uses a saddle point approximation for the distribution of
the sum of binomials. The results are presented in Table 3, and the interval
test plot for generalised gamma is presented in Figure 2.

Parametric survival model | Bonferroni rejections | Individual flags | PAVSI | TFT
---|---|---|---|---
Exponential |  [HTML]fc86721 | 4 |  [HTML]fc86720.0005 |  [HTML]fc86720.0040
Gamma |  [HTML]fc86721 | 2 |  [HTML]fc86720.0488 |  [HTML]fc86720.0052
Generalised gamma |  [HTML]fc86721 | 2 |  [HTML]fc86720.0488 |  [HTML]fc86720.0345
Gompertz |  [HTML]fc86721 | 4 |  [HTML]fc86720.0005 |  [HTML]fc86720.0056
Log-logistic |  [HTML]fc86721 | 1 | 0.2437 |  [HTML]f3aa7f0.0631
Log-normal | 0 | 1 | 0.2437 |  [HTML]f3aa7f0.0699
Weibull |  [HTML]fc86721 | 3 |  [HTML]fc86720.0063 [HTML]fc8672 |  [HTML]fc86720.0041
Table 3: Example 2: Interval test results using 10 evenly-spaced intervals.

As we can see, Table 3 tells a very different story to the previous interval
selection approach. The overall PAVSI p-value is below 0.05 for all models
except log-logistic and log-normal, indicating poor fit in 5 out of 7 models.
Similarly, the overall TFT p-value rejects 5 out of 7 models and is below 0.07
in the final two models, which again raises concerns about poor fit. Finally,
all of the parametric models have a Bonferroni rejection except log-normal,
suggesting only the log-normal is a feasible model. We investigate the log-
normal model further by presenting the results for each of its 10 intervals in
Table 4.

$V_{j}$ |  $\mathbb{E}$(Events) | Obs Events | p-mid | Individual test | Bonferroni test
---|---|---|---|---|---
(0, 4.34] | 18.76 | 11 | 0.0289 | Accept | Accept
(4.34, 8.67] | 22.52 | 23 | 0.5534 | Accept | Accept
(8.67, 13.01] | 17.90 | 26 | 0.9677 | Accept | Accept
(13.01, 17.34] | 13.67 | 21 | 0.9710 | Accept | Accept
(17.34, 21.68] | 10.70 | 11 | 0.5549 | Accept | Accept
(21.68, 26.01] | 8.99 | 7 | 0.2550 | Accept | Accept
(26.01, 30.35] | 7.71 | 9 | 0.6947 | Accept | Accept
(30.35, 34.69] | 6.12 | 1 | 0.0087 |  [HTML]dedcdcFlag | Accept
(34.69, 39.02] | 5.02 | 5 | 0.5247 | Accept | Accept
(39.02, 43.36] | 1.72 | 3 | 0.8279 | Accept | Accept
Table 4: Example 2: Interval test for the log-normal model with 10 evenly-
spaced intervals.

From Table 4, we see that the interval resulting in an individual flag for the
log-normal model was between times 30.35 and 34.69 with a p-value of 0.0087,
which is above the 0.0025 threshold of a Bonferroni rejection. Furthermore, we
note that the 4th interval had a p-value of 0.971. If this value was only
0.004 higher, then the log-normal would have received a second individual flag
and thereby a PAVSI overall p-value rejection. Furthermore, the log-normal
performed the worst of all 7 models for both AIC and BIC (not presented). For
these reasons, as well as the overall TFT p-value of 0.0699 (close to 0.05) we
can be justifiably suspicious of the log-normal model. However, we have no
strong evidence to outright reject it.

In conclusion, due to the Bonferroni rejections, we would suggest that all
presented models except for the log-normal are unfeasible. This leaves us with
only the log-normal model remaining, which is somewhat restrictive and has
also been flagged for concern as described above. For this reason, it is
recommended that exploring flexible models such as splines and piecewise would
be appropriate.

##  4 Simulation summary

A simulation experiment was performed to verify results under a true null
hypothesis. Additional details are provided in Appendix C. The simulation
experiment aimed to confirm our methodological assumptions with a simple but
realistic example. For the approach where intervals were chosen using
censoring, the test statistics had a very low type I error rate, suggesting
that they are implausibly conservative in practice. For the 10 evenly-spaced
interval approach, the test statistics performed as expected with a type I
error rate close to 0.05. The Bonferroni test performed well in both
approaches but was somewhat conservative for the censor approach. In
conclusion, the two primary assumptions hold fairly well for the 10 evenly-
spaced intervals approach. (The two primary assumptions are that censor
intervals $I_{j}$ are uninformative of the event generating process and that
the mid-point p-values are approximately uniform[0,1].) We also note that
varying the maturity of the data had only a minor impact on the evenly-spaced
interval approach.

##  5 Discussion

In terms of applying the interval test, two approaches were presented. In the
first, intervals of interest are defined entirely by the censors. In the
second, intervals of interest were specified as 10 evenly-spaced intervals up
to the final censor. There are strengths and weaknesses to both approaches. We
address this in the discussion, where the following topics are covered: the
censor interval approach, how to specify intervals, the benefits of the 10
evenly-spaced interval approach, testing with limited patient numbers, and a
conclusion.

In our first approach, intervals are selected by the data using the censor
times. Unfortunately, this approach is generally unhelpful in describing
regions with heavy censoring. This is because very small intervals are
unlikely to be rejected due to low power. Additionally, with this approach,
our test statistics are very conservative in standard cases (see Appendix C).
On the other hand, we argue that in some cases (not presented) this approach
is indicative of informative or unusual censoring. That is, it is very
apparent if a lot of censoring has occurred, and furthermore, we can see if an
unusual number of events occurred within these intervals. As a final point,
this approach underlies the specified interval approach, so it will always
provide helpful supporting evidence.

In the second approach, intervals can be specified. There are several ways
this could be done. In our example, we provide one approach (10 evenly-spaced
intervals). If submitting to a health technology assessment body, we recommend
using this approach. This prevents intervals from being systematically chosen
such that model selection favours a treatment’s efficacy. However, in terms of
further research, there are certainly additional options for how the intervals
are chosen. One option is automatic interval selection based on equal
statistical power. Another topic of consideration is restricting chosen
intervals to have boundaries that are censor times.

There are several benefits to the 10 evenly-spaced intervals approach compared
to the censor interval approach. The 10 intervals are generally easier to
interpret with plots. Furthermore, the approach is informative for areas with
heavy censoring (such as the tail of the KM curve), which is often the area of
greatest importance for long-term extrapolation. Additionally, overall
p-values perform much better for this method. It is for these reasons that in
general we recommend the fixed interval approach. However, we acknowledge that
the censor interval approach is still valid and has its own benefits.

In situations with limited data, say under 20 patients, the tests can still be
used. However, they will likely be under-powered and accept the fitted model.
This is not necessarily a bad thing. When we look at data with only 20
patients, visual inspection of curve fit can be challenging, and after a
subjective evaluation, we may think that flexible survival models are
appropriate. That is, we have potentially over-fitted the data. The interval
tests can prevent this by providing an objective evaluation. The tests will
likely confirm that the observed data is feasible for simple models,
indicating that flexible models are not required.

In conclusion, the interval tests provide additional information that can
support the model selection process. The tests are able to highlight areas of
particularly poor fit. It is recommended that they are especially useful in
aiding decisions when the curve fit is poor.

##  6 Appendix

##  Appendix A In defense of midpoint p-values

When deriving test statistics and results, we will often make the assumption
that midpoint p-values are uniform[0,1] under H0. If the validity of this
approach is in question, then randomized p-values can be used, which are known
to be uniform[0,1] under H0[13]. For this reason, the randomized p-values are
superior from this purely statistical standpoint. However, we argue below that
the midpoint p-values are more practical and that for our applications the
uniform[0,1] assumption is often reasonable.

In the case of reducing the number of intervals to say 10, as is described in
example 2, the generating process is the sum of binomials. If the number of
patients is reasonable, one would expect these intervals to contain a
reasonable number of feasible discrete outcomes. This means the uniform[0,1]
assumption is more valid, compared to when there are only a small number of
possible outcomes.

We also acknowledge that in practice, randomized p-values are more difficult
to deal with. When an analysis is rerun, the p-values are rerandomized. This
would then involve rerunning the analysis to determine the level of variation
in the randomised p-values and their resulting test statistic. This is
helpful. However, this approach leads you down the path of using the expected
value of these reruns as your final quoted value. That is, the midpoint
p-value, which is the expected value of randomized p-value.

##  Appendix B Extension to individual patient models

In the introduction, we mentioned that the most commonly used calibration
methods rely on risk scores [4]. However, in our approach, risk scores were
not considered. Furthermore, in our examples, all of our models inherently
assume equal risk scores for each patient. It is, however, possible to extend
our methodology to models where the risk scores are not assumed constant. This
would involve amending Equation 1 such that each individual is considered
separately for each interval. That is, every patient is assigned a specific
probability of experiencing an event (as defined by the model) for that
interval. Subsequently, Equation 1 would be replaced by a collection of
$N_{I_{j}}$ Bernoulli random variables with different probabilities.

##  Appendix C Additional details on simulation

In practice, the parametric model never truly underlies the data, meaning type
I errors don’t make logical sense. However, it does provide an interesting
thought-experiment; if the observed data was generated by its fitted model, is
the observed data behaving as you would expect? If the answer is yes, then the
model is reasonable. If the answer is no, then the model is likely
unreasonable.

The primary purpose of this investigation by simulation was to confirm that
the methodology behaves as expected with respect to type I errors. The
realistic model was described by an event generating process $T\sim
exp(\lambda)$ and a censor generating process for $C$ given by

$\begin{split}&C1\sim uniform[0,100],\\\ &C2\sim uniform[18,22],\\\
&C=min(C1,C2).\end{split}$

Finally, $Time=min(C,T)$; and if $T\leq C$, $Event=1$, else $Event=0$. Taking
the time to be in months, we are saying some patients are heavily censored in
the last 4 months of the trial, patients have an expected event time of
$1/\lambda$, and some patients are censored throughout the trial.

For the censored interval approach, test statistics were shown to be
unrealistically conservative. This is due to an inflated number of interval
p-values that are close to 0.5, which occurs due to using the midpoint p-value
approach. That is, most of the small intervals will yield a midpoint p-value
of 0.5. These are then transformed for the TFT to give effective values of 1,
which inflates the overall TFT p-value, pushing it closer to 1. In general,
more censor intervals result in the tests being even more conservative. In our
simulation, with 200 patients, $\lambda=1/{10}$, and 10,000 iterations the
PAVSI type I error rate was 0.0041 and for TFT this was 0.0001. This is not
ideal as it indicates that in practice these test statistics would be overly
conservative. Thereby, both test statistics are not recommended in practice
when using the censor interval approach.

For the 10 evenly-spaced intervals approach, the test statistics behaved as
expected. When the number of patients increased, the test statistics appeared
to converge to the expected values. Even with small numbers of patients, the
test statistics remained reasonable. For example, with 50 patients,
$\lambda=1/{10}$, and 10,000 iterations the PAVSI type I error rate was 0.0331
and for TFT this was 0.0146.

For the 10 evenly-spaced interval approach with 200 patients, the PAVSI type I
error rate was 0.071 and for TFT this was 0.04. Note that here the PAVSI type
I error rate is above 0.05. This is as expected. In particular, when assuming
the interval p-values are uniform[0,1], the type I error is 0.08614. (This
result can be shown by considering our approach, which uses 10 intervals, a
0.05 rejection level, and a test statistic obtained by using the mid-point
p-value.) As the number of patients increases, our PAVSI test statistic
appears to converge towards this value, 0.08614. That is, the midpoint
p-values are behaving as uniform[0,1] random variables. Similarly, the TFT
appears to converge towards 0.05 as the number of patients increases.

We make a quick digression here to discuss the PAVSI type I error converging
towards 0.08614. If we were to select a different number of intervals, say 7,
then the PAVSI type I error would converge towards 0.0444 instead. However,
this is not ideal. With lower patient numbers, this value decreases even
further due to the midpoint p-values. Thereby the 7 interval approach would
always be conservative, which we argue is a bad property for model selection.
The 10 evenly-spaced interval method approaches a type 1 error of 0.08614 as
the patients increase. However, a type I error slightly above 0.05 isn’t
necessarily a bad thing, especially in practice. Also, note that when using 10
intervals with a simulation of 100 patients, the PAVSI type I error rate was
0.0556, which is approximately what we want as it is close to 0.05.

The Bonferroni type I error rates were approximately 0.02 for the censored
interval approach and 0.04 for the 10 evenly-spaced interval approach. These
values appear to remain approximately stable as the number of patients
changes. That is, approximately within $(0.03,0.05)$ for the 10 evenly-spaced
interval approach and within (0.017,0.022) for the censor approach, when
considering 50 to 2000 patients with $\lambda=1/{10}$.

Finally, we experimented with different levels of maturity in the data. In
particular, we used $\lambda=1/{10}$ (mature data), $\lambda=1/{30}$ (somewhat
mature), and $\lambda=1/{70}$ (immature). For these simulated data, the KM
curves finished at approximate survival levels of 0.11, 0.48, and 0.73
respectively. In the censored interval approach, when the data was less
mature, the test statistics performed worse. (They became even more
conservative.) This is due to an increased amount of censoring in the less
mature data. For 10 evenly-spaced interval approach, the different levels of
maturity in the data had a somewhat minimal impact on the type I error rates.
In particular, when weakening the data from mature to immature
($\lambda=1/{10}$ to $\lambda=1/{70}$), this would approximately result in a
decrease of around 0.01 in the type 1 error rate for the PAVSI, TFT and
Bonferroni tests.

## References

  * [1] Mark J Rutherford et al. (2020) _NICE DSU Technical Support Document 21: Flexible Methods for Survival Analysis_.
  * [2] Nicholas Latimer. (2011) _NICE DSU Technical Support Document 14: Survival Analysis for Economic Evaluations alongside Clinical Trials - Extrapolation with Patient-Level Data_.
  * [3] David W. Hosmer and Stanley Lemeshow. (2000) _Applied Logistic Regression 2nd edition_ , Wiley series in probability and statistics: Section 5.2.2
  * [4] Olga V. Demler, Nina P. Paynter, and Nancy R. Cook. (2015) _Tests of Calibration and Goodness of Fit in the Survival Setting_ , Statistics in Medicine, Vol.34: 1659–1680
  * [5] Grønnesby JK, Borgan Ø. (1996) _A method for checking regression models in survival analysis based on the risk score_ , Lifetime Data Analysis. 2:315–320. 9
  * [6] D’Agostino RB and Byung-Ho N. (2004) _Evaluation of the performance of survival analysis models: discrimination and calibration measures_ , Handbook of statistics. 2004; 23:1–25
  * [7] Peter C. Austin, Frank E. Harrell, and David van Klaveren. (2020) _Graphical calibration curves and the integrated calibration index (ICI) for survival models_ , Statistics in Medicine, Vol.39: 2714–2742
  * [8] Allan Birnbaum (1954) _Combining Independent Tests of Significance,_ Journal of the American Statistical Association, Vol 49: No 267 559-574
  * [9] Nicholas A. Heard and Patrick Rubin-Delanchy (2018) _Choosing Between Methods of Combining p-values,_ Biometrika, Vol 105: 239-246
  * [10] Suzanne C Freeman et al. (2022) _Challenges of modelling approaches for network meta-analysis of time-to-event outcomes in the presence of non-proportional hazards to aid decision making: Application to a melanoma network_ , Statistical Methods in Medical Research, Vol. 31(5) 839–861
  * [11] Axel Hauschild et al. (2020) _Long-term outcomes in patients with BRAF V600-mutant metastatic melanoma receiving dabrafenib monotherapy: Analysis from phase 2 and 3 clinical trials_ , European Journal of Cancer, Vol. 125:114-120.
  * [12] Georgina Long et al. (2017) _Dabrafenib plus trametinib versus dabrafenib monotherapy in patients with metastatic BRAF V600E/K-mutant melanoma: long-term survival and safety analysis of a phase 3 study_ , Annal of Oncology, Vol. 28(7):1631-1639.
  * [13] Patrick Rubin-Delanchy, Nicholas A. Heard, and Daniel J. Lawson (2019) _Meta-Analysis of Mid-p-Values: Some New Results based on the Convex Order,_ Journal of the American Statistical Association, Vol 114, No. 527: 1105–1112: Theory and Methods
  * [14] Richard A. Armstrong. (2014) _When to use the Bonferroni correction,_ Ophthalmic and Physiological Optics, Vol 34: 502-508
  * [15] Christopher Jackson. (2016) _flexsurv: A Platform for Parametric Survival Modeling in R_ , Journal of Statistical Software, Vol. 70 Issue 8
  * [16] Boxiang Liu and Thomas Quertermous. (2017) _Sinib: Approximating the Sum of Independent Non-Identical Binomial Random Variables,_ arXiv