# MSM_pseudo_CS

Samuel Anyaso-Samuel1 and Somnath Datta2
1Division of Cancer Epidemiology and Genetics, National Cancer Institute,
Rockville, MD, USA
2Department of Biostatistics, University of Florida, Gainesville, FL, USA

# Nonparametric estimation of a future entry time distribution given the
knowledge of a past state occupation in a progressive multistate model with
current status data

Samuel Anyaso-Samuel1 and Somnath Datta2
1Division of Cancer Epidemiology and Genetics, National Cancer Institute,
Rockville, MD, USA
2Department of Biostatistics, University of Florida, Gainesville, FL, USA

###### Abstract

Case-I interval-censored (current status) data from multistate systems are
often encountered in cancer and other epidemiological studies. In this
article, we focus on the problem of estimating state entry distribution and
occupation probabilities, contingent on a preceding state occupation. This
endeavor is particularly complex owing to the inherent challenge of the
unavailability of directly observed counts of individuals at risk of
transitioning from a state, due to the cross-sectional nature of the data. We
propose two nonparametric approaches, one using the fractional at-risk set
approach recently adopted in the right-censoring framework and the other a new
estimator based on the ratio of marginal state occupation probabilities. Both
estimation approaches utilize innovative applications of concepts from the
competing risks paradigm. The finite-sample behavior of the proposed
estimators is studied via extensive simulation studies where we show that the
estimators based on severely censored current status data have good
performance when compared with those based on complete data. We demonstrate
the application of the two methods to analyze data from patients diagnosed
with breast cancer.

##  1 Introduction

Multistate models are a popular type of survival model that allows for
studying the progression of the disease using multivariate time-to-event data.
In biomedical applications, situations arise in which patients progress
reversibly or irreversibly through a succession of disease stages that
describe their current condition. The patients usually begin from an initial
disease state and move through one or more intermediate states before entering
a final absorbing state. The disease stages correspond to the states of the
multistate system. Moreover, in certain applications, the transition of
patients through their multistate system may be subject to a severe form of
(case-I) interval censoring often referred to as current status censoring
(Jewell and van der Laan, 2003; Sun, 2006). Under such an observation scheme,
the researcher is only restricted to information on the state occupied by the
patient at a (single) random inspection time during the entire study period.

Current status data from multistate models often arise in epidemiological and
biomedical research, as well as in observational studies where continuous
monitoring of study participants is not feasible. With event history data from
a multistate system, practitioners are often interested in estimating the
probability of occupying a given state, at a specific time. Datta and Sundaram
(2006) proposed marginal estimators for state occupation probabilities
following the product integration approach of Aalen and Johansen (1978). Their
nonparametric approach is valid without imposing the Markov assumption, and it
has been further studied under different settings (Lan and Datta, 2010; Lan et
al., 2017; Anyaso-Samuel et al., 2023).

In progressive multistate systems, where individuals navigate through
successive states within a disease process, a pertinent question arises: “What
is the probability of occupying a future state $k$, given a prior visit to
state $j$?” We denote this conditional probability as $\Psi_{k|j}$.
Furthermore, interest may extend to studying the probability of reaching a
future state $k$ by time $t$, given a previous visit to state $j$, represented
as $\Psi_{k|j}(t)$. Historically, these conditional probabilities have
garnered less attention, as multistate systems are often assumed to adhere to
Markovity. However, practitioners increasingly recognize the limitations of
imposing a Markov assumption, which may be unverifiable and violated. Given a
progressive multistate system, understanding the proportion of individuals on
a unique path capable of occupying a future state becomes crucial. Our paper
addresses this gap by focusing on the conditional estimation of future state
occupation probabilities given prior occupation of a specific state, within
the context of current status censoring.

In the context of right-censored data, Yang et al. (2023) recently proposed an
estimator for the future state occupational probability conditional on a past
state visit. Constructing such estimators is difficult since information on
the transitions out of the non-absorbing states within the multistate system
is only partially observed. With current status data, the problem becomes even
more challenging since the path that individuals would have taken after they
were observed is unknown. To tackle this problem of estimating $\Psi_{k|j}$
and $\Psi_{k|j}(t)$, we propose two estimation methods. The first is based on
fractional at-risk sets. Here, we formulate the at-risk processes with weights
defined by the fractional observation that a person contributes to the at-risk
set of transitioning out of state $j$. We estimate this fractional weight as
the probability that a patient will eventually reach state $j$ given that they
were occupying a known state at a random inspection time. The idea for
utilizing fractional weighted at-risk sets follows from the seminal work of
Datta and Satten (2000) and the work of Yang et al. (2023). We begin the
construction of the second estimator by providing insight that the conditional
quantity $\Psi_{k|j}$ can be expressed as the ratio of two marginal state
occupation probabilities for multistate models with a tree structure. This
leads to the formulation of an estimator based on product limit estimators of
state occupation probabilities (Datta and Sundaram, 2006). For both estimation
methods, we do not impose the assumption that the multistate systems are
Markov, but the methods require that the systems have a directed tree
structure. Following our construction of estimators for the conditional future
state occupation probabilities, we derive estimators for the distribution
function of the state entry conditional on a prior occupation of a given
state.

The rest of the paper is organized as follows. In Section 2, we introduce the
notation along with marginal nonparametric estimators for current status
multistate data. In Section 3, we discuss the estimation of $\Psi_{k|j}$ based
on the fractional weight approach. The alternate estimator based on the
product-limit method is discussed in Section 4. In Section 5, we describe the
estimation of pointwise confidence intervals for the estimators based on both
approaches. Detailed simulation studies to evaluate the performance of the
proposed estimators are presented in Section 6. In Section 7, we illustrate
the application of both estimation methods in analyzing data from a breast
cancer study. The paper concludes with a discussion in Section 8.

##  2 Nonparametric estimators for a multistate system

###  2.1 Notation and Convention

Consider a multistate system with a directed tracking structure, where
$\mathcal{S}=\\{0,1,\ldots,Q\\}$ represents the finite state space. In this
system, each individual $i=1,\ldots,n$ occupies a state denoted by
$S_{i}(t)\in\mathcal{S}$ at time $t\geq 0$. Additionally, $C_{i}$ stands for
the random inspection time for individual $i$. For an individual, the observed
current-status data is represented as the pair $\\{C_{i},S_{i}(C_{i})\\}$,
where $C_{i}$ is assumed to be independent of the multistate process
$\\{S_{i}(t),t\geq 0\\}$, based on the random censoring assumption. To model
the multistate current-status data, we adopt a framework where each state
$j\in\mathcal{S}$ can be reached from the initial state $0$ through a unique
path: $0=s_{1}\rightarrow s_{2}\rightarrow\ldots\rightarrow s_{(j)+1}=j$,
where $(j)$ denotes the number of steps it takes to reach $j$ from the root
node 0.

Let $U_{jj^{\prime}}$ denote the (unobserved) transition time at which an
individual transitions from $j$ to $j^{\prime}$. To estimate temporal
functions of the multistate system, two essential functionals are the counting
process denoted by $N_{jj^{\prime}}(t)$ and the at-risk process denoted by
$Y_{j}(t)$. The counting process $N_{jj^{\prime}}(t)$ counts the number of
transitions from state $j$ to $j^{\prime}$ in the time interval $[0,t)$. On
the other hand, the at-risk process $Y_{j}(t)$ represents the number of
individuals who are at risk of transitioning out of state $j$ at time $t$.
These processes are expressed as follows

$\displaystyle N_{jj^{\prime}}(t)$
$\displaystyle=\sum^{n}_{i=1}\mathcal{I}(U_{i,jj^{\prime}}\leq t),\quad\text{
and }$ $\displaystyle Y_{j}(t)$
$\displaystyle=\sum^{n}_{i=1}\mathcal{I}(S_{i}(t-)=j),$

where $\mathcal{I}(\cdot)$ denotes the indicator function and $X(t-)$ is the
state occupied just before time $t$. By the law of large numbers,

$\displaystyle
n^{-1}N_{jj^{\prime}}(t)\xrightarrow{P}\text{Pr}(U_{jj^{\prime}}\leq t),$

for any $t\geq 0$. Since the inspection time $C$ is independent of the
multistate process, we can write

$\displaystyle\text{Pr}(U_{jj^{\prime}}\leq
t)=\mathbb{E}(\mathcal{I}(U_{jj^{\prime}}\leq C)\ |\ C=t),$

where $\mathcal{I}(U_{jj^{\prime}}\leq C)$ is the indicator of the event that
the transition from state $j$ to $j^{\prime}$ has taken place by time $C$.

Using the theory of nonparametric regression, Datta and Sundaram (2006)
proposed nonparametric estimators of $N_{jj^{\prime}}(t)$ and $Y_{j}(t)$ under
this context of current status data. For $N_{jj^{\prime}}(t)$, they proposed a
two-step estimation approach where first, isotonic regression via the pooled
adjacent violators (PAV) algorithm (Barlow et al., 1972) is performed based on
the pairs $(C_{i},\ \mathcal{I}(U_{i,jj^{\prime}}\leq C_{i}))$. Let
$\widehat{N}^{P}_{jj^{\prime}}(t)$ denote the fitted responses from the
isotonic regression. Then secondly, kernel smoothing (Mukerjee, 1988;
Nadaraya, 1964; Watson, 1964) with a log-concave density, $K>0$ is performed
to remove long flat parts of $\widehat{N}^{P}_{jj^{\prime}}(t)$. The resulting
estimator of $N_{jj^{\prime}}(t)$ is given by

$\displaystyle\widehat{N}_{jj^{\prime}}(t)=\sum^{n}_{i=1}\left\\{\widehat{N}^{P}_{jj^{\prime}}(t)\times\frac{K_{h}(C_{i}-t)}{\sum^{n}_{i=1}K_{h}(C_{i}-t)}\right\\},$
(1)

where $K_{h}(\cdot)=h^{-1}K(\cdot/h)$ with a bandwidth sequence
$0<h=h(n)\downarrow 0$ chosen by the criteria specified by Wand and Jones
(1995). Moreover, to estimate $\text{Pr}(X(t-)=j)$, which represents the
probability limit of $n^{-1}Y_{j}(t)$, Datta and Sundaram (2006) employed the
kernel smoothing approach. Accordingly, the estimation of $Y_{j}(t)$ is given
by:

$\displaystyle\widehat{Y}_{j}(t)=\sum^{n}_{i=1}\left\\{\mathcal{I}(S_{i}(C_{i})=j)\times\frac{K_{h}(C_{i}-t)}{n^{-1}\sum^{n}_{i=1}K_{h}(C_{i}-t)}\right\\},$
(2)

where it is assumed that $S_{i}(C_{i}-)=S_{i}(C_{i})$ with probability 1. The
estimators outlined in (1) and (2) pave the way for feasible inference in
multistate models with current status data. For example, Datta and Sundaram
(2006) elaborates on the estimation of key marginal functionals, such as the
transition probabilities and state occupation probabilities, employing the
product limit approach described by Andersen et al. (1993).

This paper primarily proposes estimators for future state occupation
probabilities and state entry distributions conditioned on a prior state
visit. We employ the methodology described above to estimate the counting
process and at-risk set, which serves as the foundational framework for
constructing our proposed conditional estimators.

##  3 Estimators based on fractional at-risk sets

###  3.1 An illustration using a five-state model

To provide context for the aim of this study, let’s examine the progressive
five-state illness-death model depicted in Figure 1. Our focus centers on
estimating the probability of ever visiting state 3, given a prior stay in
state 1 over an individual’s lifespan. We denote this conditional future
occupation probability as $\Psi_{3|1}$. Estimating $\Psi_{3|1}$ poses a
challenge due to incomplete observation of transition times. With current
status data, we only have access to information about the state a patient
occupies at a random inspection time.

0\. Disease-free 1\. Ill 2\. Death without illness 3\. Death after illness 4\.
Death from other causes State 0* Figure 1: A five-state illness-death model.

Let $X_{ij}$ denote the (unobserved) indicator of the event that the $i$th
individual would ever enter state $j$. In the typical current status
framework, $X_{ij}$ will likely be unknown, making estimating the future entry
states difficult. For instance, if the $i$th individual was observed to have
$S_{i}(C_{i})=0$, then $X_{i1}$ will be unobserved. Suppose the $i$th
individual transitioned to state 1, they will be eligible for the competing
risks (state 3 and state 4), but we would not know this because $X_{i1}$ is
not observed. To overcome the difficulty in estimating $\Psi_{3|1}$, first, we
define a new state $0^{*}$ which is constructed by combining states 0, 1, and
2 (see Figure 1). Note that not everyone (i.e., those who transitioned to
state 2) in state $0^{*}$ will be subject to the competing risks.

Next, we compute the fractional observation, denoted as
$\phi_{i1}(C_{i},S_{i}(C_{i}))$, representing the contribution of the $i$th
individual to the “at-risk” set for transitioning out of state 1. This
fractional observation estimates the probability of the $i$th individual
eventually reaching state 1 from state 0. Computing
$\phi_{i1}(C_{i},S_{i}(C_{i}))$ is straightforward when $S_{i}(C_{i})\neq 0$.
For $S_{i}(C_{i})=0$, we compute $\phi_{i1}(C_{i},S_{i}(C_{i}))$ as the
transition probability of moving from state 0 to state 1, that is,
$\phi_{i1}(C_{i},S_{i}(C_{i}))={P}_{01}(C_{i},\infty)$. Using the Aalen-
Johansen formula for transition probability, we estimate
$\phi_{i1}(C_{i},S_{i}(C_{i}))$ by the following expression

$\displaystyle\widehat{\phi}_{i1}(C_{i},S_{i}(C_{i}))=\int_{(C_{i},\infty)}\left\\{\prod_{(C_{i},u)}\left[1-\frac{d\widehat{N}_{01}(v)+d\widehat{N}_{02}(v)}{\widehat{Y}_{0}(v)}\right]\right\\}\frac{d\widehat{N}_{01}(u)}{\widehat{Y}_{0}(u)}du,$
(3)

where $\widehat{N}_{\ell\ell^{\prime}}(t)$ and $\widehat{Y}_{\ell}(t)$ are
obtained by (1) and (2), respectively. The calculation of (3) does not require
the Markov assumption since it pertains to the transition out of the root
node. If $S_{i}(C_{i})=2$, then $\widehat{\phi}_{i1}(C_{i},S_{i}(C_{i}))=0$.
For $S_{i}(C_{i})\in\\{1,3,4\\}$, indicating that the individual followed the
$0\rightarrow 1$ path, $\widehat{\phi}_{i1}(C_{i},S_{i}(C_{i}))=1$.

To obtain the estimator of $\Psi_{3|1}$, we introduce a modified multistate
system that includes state 0* as the initial stage, leading to either state 3
or state 4 (see Figure 1). Applying the Aalen-Johansen formula, commonly used
in competing risks calculations, to this modified multistate system yields an
expression that serves as the estimator for eventually reaching state 3:

$\displaystyle\widehat{\Psi}^{[1]}_{3|1}=\int^{\infty}_{0}\left\\{\prod_{(0,u)}\left[1-\frac{d\widehat{N}_{3|1}(v)+d\widehat{N}_{4|1}(v)}{\widehat{Y}_{0^{*}|1}(v)}\right]\right\\}\frac{d\widehat{N}_{3|1}(u)}{\widehat{Y}_{0^{*}|1}(u)}du.$
(4)

Here the superscript $[1]$ is used to distinguish this estimator from the
alternate estimator discussed in Section 4. The estimators of the counting
processes and the at-risk set in (4) are defined by

$\displaystyle\widehat{N}_{3|1}(t)$
$\displaystyle=\sum^{n}_{i=1}\left\\{\frac{\widehat{N}^{P}_{13}(t)K_{h}(C_{i}-t)}{\sum^{n}_{i=1}K_{h}(C_{i}-t)}\right\\},$
$\displaystyle\widehat{N}_{4|1}(t)$
$\displaystyle=\sum^{n}_{i=1}\left\\{\frac{\widehat{N}^{P}_{14}(t)K_{h}(C_{i}-t)}{\sum^{n}_{i=1}K_{h}(C_{i}-t)}\right\\},$
$\displaystyle\widehat{Y}_{0^{*}|1}(t)$
$\displaystyle=\sum^{n}_{i=1}\left\\{\widehat{\phi}_{i1}(C_{i},S_{i}(C_{i}))\times\frac{\mathcal{I}(S_{i}(C_{i})\in\\{0,1,2\\})\times
K_{h}(C_{i}-t)}{n^{-1}\sum^{n}_{i=1}K_{h}(C_{i}-t)}\right\\},$

where $n^{-1}\widehat{N}^{P}_{jj^{\prime}}(\cdot)$ is a step function for
taking values $n^{-1}\widehat{N}^{P}_{jj^{\prime}}(C_{(i)})=R_{i}$, that
minimizes the sum of squares
$\sum_{i=1}^{n}\Big{\\{}R_{i}-\mathcal{I}(U_{[i],jj^{\prime}}\leq
C_{(i)})\Big{\\}}^{2}$ subject to $R_{1}\leq...\leq R_{n}$, and $[i]$
corresponds to the index such that $C_{[i]}$ equals the $i$th-order statistic
$C_{(i)}$. Subsequently, to obtain, $\widehat{N}_{j^{\prime}|j}(\cdot)$,
similar to the approach used in constructing (1), we use kernel smoothing to
remove the long flat part from $\widehat{N}^{P}_{jj^{\prime}}(\cdot)$. Note,
the estimator, $\widehat{Y}_{0^{*}|1}(t)$ is formulated by the estimate of the
fractional weight, $\widehat{\phi}_{i1}(C_{i},S_{i}(C_{i}))$. Here,
$\widehat{Y}_{0^{*}|1}(t)$ estimates the number of persons at risk of
transitioning out of the artificial state $0^{*}$ at time $t$.

Additionally, the estimation of the entry time distribution into state 3,
given a prior visit to state 1, is expressed as follows:

$\displaystyle\widehat{F}_{3|1}(t)=\frac{\widehat{\Psi}^{[1]}_{3|1}(t)}{\widehat{\Psi}^{[1]}_{3|1}}.$

Here, the estimation of the subdistribution function
$\widehat{\Psi}^{[1]}_{3|1}(t)$ is achieved by constraining the integral’s
bounds as defined in the expression provided in Equation 4. Therefore,

$\displaystyle\widehat{\Psi}^{[1]}_{3|1}(t)=\int^{t}_{0}\left\\{\prod_{(0,u)}\left[1-\frac{d\widehat{N}_{3|1}(v)+d\widehat{N}_{4|1}(v)}{\widehat{Y}_{0^{*}|1}(v)}\right]\right\\}\frac{d\widehat{N}_{3|1}(u)}{\widehat{Y}_{0^{*}|1}(u)}du.$
(5)

###  3.2 Extension to general multistate model with a tree structure

Next, we expand the explanation of the estimator for
$\Psi_{k|j}=P(X_{k}=1~{}|~{}X_{j}=1)$ to encompass a broader multistate model
context. The estimator (based on fractional weights) of the probability of
occupying state $k$ given a prior visit to state $j$ is denoted as
$\widehat{\Psi}^{[1]}_{k|j}$. The estimation of ${\Psi}_{k|j}$ relies on the
assumption that state $j$ lies on the unique path connecting states $0$ and
$k$; otherwise, ${\Psi}_{k|j}=0$. Let $m$ represent the number of transitions
needed to move from state $0$ to state $k$. Constructing the estimator for
$\Psi_{k|j}$ follows a recursive approach based on $m$.

When $m=1$, the computation of $\widehat{\Psi}^{[1]}_{k|j}$ involves using the
Aalen-Johansen estimator for transition probabilities within a competing risks
framework. On the other hand, when $m=2$, then, one can obtain
$\widehat{\Psi}^{[1]}_{k|j}$ by following the procedure outlined in Section
3.1. More generally, for any state $k$ that is situated at a distance of $m$
from state 0, the estimation of ${\Psi}_{k|j}$ relies on the application of
the chain rule of conditional probability, given that the states of interest
are aligned along a distinct path originating from the root node.

Following the application of the chain rule, for every state $j$ on the path
from state 0 to $k$, where $k$ is a distance $m$ from state 0; define

$\displaystyle\widehat{\Psi}^{[1]}_{k|j}=\widehat{\Psi}^{[1]}_{\tilde{k}|j}\cdot\widehat{\Psi}^{[1]}_{k|\tilde{k}}$
(6)

where $\tilde{k}$ is the state just before state $k$ on the path from $0$ to
$k$. From this recursive formula, $\widehat{\Psi}^{[1]}_{\tilde{k}|j}$ is
available since the distance from state 0 to state $\tilde{k}$ is less than
$m$. To construct the estimator $\widehat{\Psi}^{[1]}_{k|\tilde{k}}$, we
utilize the idea of estimating transition probabilities in the competing risk
framework and the concept of “pooling of states”. Thus, we have

$\displaystyle\widehat{\Psi}^{[1]}_{k|\tilde{k}}=\int^{\infty}_{0}\widehat{S}_{0^{*}|\tilde{k}}(s-)\frac{d\widehat{N}_{k|\tilde{k}}(s)}{\widehat{Y}_{0^{*}|\tilde{k}}(s)}.$
(7)

where

$\displaystyle\widehat{S}_{0^{*}|\tilde{k}}(t)$
$\displaystyle=\prod_{(0,t)}\left[1-\frac{d\widehat{N}_{.|\tilde{k}}(v)}{\widehat{Y}_{0^{*}|\tilde{k}}(v)}\right].$

Here, state $0^{*}$ denotes an artificial state constructed by combining the
states preceding (and including) $\tilde{k}$. Also,

$\displaystyle\widehat{N}_{k|\tilde{k}}(t)$
$\displaystyle=\sum^{n}_{i=1}\left\\{\frac{\widehat{N}^{P}_{\tilde{k}k}(t)K_{h}(C_{i}-t)}{\sum^{n}_{i=1}K_{h}(C_{i}-t)}\right\\}\text{
and }\widehat{N}_{.|\tilde{k}}(t)=\sum_{k}\widehat{N}_{k|\tilde{k}}(t),$
$\displaystyle\widehat{Y}_{0^{*}|\tilde{k}}(t)$
$\displaystyle=\sum^{n}_{i=1}\left\\{\widehat{\phi}_{i\tilde{k}}\times\frac{\mathcal{I}(S_{i}(C_{i})\in
0^{*})K_{h}(C_{i}-t)}{n^{-1}\sum^{n}_{i=1}K_{h}(C_{i}-t)}\right\\},$

where $\widehat{N}^{P}_{\tilde{k}k}(\cdot)$ is obtained by the isotonic
regression of the pairs $(C_{i},\ \mathcal{I}(U_{i,\tilde{k}k}\leq C_{i}))$.
The estimator of the number of persons at risk in state $0^{*}$,
$\widehat{Y}_{0^{*}|\tilde{k}}(t)$, is based on the fractional weight
$\widehat{\phi}_{i\tilde{k}}(C_{i},S_{i}(C_{i}))$ that estimates the
probability of individual $i$ reaching state $\tilde{k}$ based on the
information observed at their inspection time $C_{i}$. We interpret
$\widehat{\phi}_{i\tilde{k}}(C_{i},S_{i}(C_{i}))$ as the contribution of the
$i$th individual to the “at-risk” set of transitioning out of state
$\tilde{k}$. If individual $i$ is observed to occupy state $\tilde{k}$ at
their inspection time, then,
$\widehat{\phi}_{i\tilde{k}}(C_{i},S_{i}(C_{i}))=1$, otherwise if they were
observed to occupy a state not on the path from state 0 to state $\tilde{k}$,
then $\widehat{\phi}_{i\tilde{k}}(C_{i},S_{i}(C_{i}))=0$. The estimation of
$\widehat{\phi}_{i\tilde{k}}(C_{i},S_{i}(C_{i}))=\widehat{\mathbb{P}}(X_{i\tilde{k}}=1\
|\ S_{i}(C_{i})=\ell)$ only becomes intricate when such individual occupies a
state $\ell$ which lies on the path from state 0 to state $\tilde{k}$ at their
inspection time. To calculate
$\widehat{\phi}_{i\tilde{k}}(C_{i},S_{i}(C_{i}))$ for such cases, we utilize a
competing risk calculation where the states prior of $\tilde{k}$ are combined
to form an event-free state.

Further, we describe the estimation of the entry distribution for state $k$
given a previous occupancy of state $j$, denoted as $\widehat{F}_{k|j}(t)$.
Let $\mathcal{S}^{k}$ denote the set of states $\ell$ such that state $k$ is
on the path from state 0 to state $\ell$. Then, $\widehat{F}_{k|j}(t)$ can be
expressed as:

$\displaystyle\widehat{F}_{k|j}(t)=\frac{\sum_{\ell\in\mathcal{S}^{k}}\widehat{\Psi}^{[1]}_{\ell|j}(t)}{\sum_{\ell\in\mathcal{S}^{k}}\widehat{\Psi}^{[1]}_{\ell|j}},$
(8)

where
$\widehat{\Psi}^{[1]}_{k|j}(t)=\widehat{\Psi}^{[1]}_{\widetilde{k}|j}(t)\cdot\widehat{\Psi}^{[1]}_{k|\widetilde{k}}(t)$,
as defined in Eqn. (6). Note that, $\widehat{F}_{k|j}(t)$ is the normalized
sum of estimated future occupation probabilities (conditioned on a prior visit
to state $j$) of state $k$ and all other states that come after state $k$ in
the progressive system.

For a practical illustration, consider a seven-state multistate system (refer
to Figure 4). Suppose we aim to estimate the probability of occupying state 5
given a prior visit to state 1, denoted by $\Psi_{5|1}$. Utilizing Eqn. (6),
in this case, we have $k=5$, $j=1$, and $\widetilde{k}=3$. Consequently, our
goal is to compute
$\widehat{\Psi}^{[1]}_{5|1}=\widehat{\Psi}^{[1]}_{3|1}\cdot\widehat{\Psi}^{[1]}_{5|3}$.
Following the recursive calculation, $\widehat{\Psi}^{[1]}_{3|1}$ is already
available. To calculate $\widehat{\Psi}^{[1]}_{5|3}$ using (7), define the
artificial state $0^{*}=\\{0,1,2,3,4\\}$. Moreover, for the non-trivial case,
the fractional weight $\widehat{\phi}_{i3}$ is obtained by estimating the
probability of transitioning from another artificial state defined by pooling
states $\\{0,1,2\\}$, to state 3.

0\. Healthy 1\. Mild COPD 2\. Death without COPD 3\. Moderate COPD 4\. Death
following mild COPD 5\. Severe COPD 6\. Death following moderate COPD Figure
2: A seven-state multistate system.

##  4 Product limit estimators

We propose a novel estimator for estimating the conditional state occupation
probability $\Psi_{k|j}$. This estimator is defined in terms of marginal state
occupation probabilities. Due to the tree structure of the multistate system,
any state will be reached from the initial state 0 by a unique path. Define
$\mathcal{S}^{k}$ as the set of states $\ell$ such that state $k$ is on the
path from 0 to $\ell$. Then, the conditional quantity, $\Psi_{k|j}$, can be
expressed as the ratio of two marginal functions given by

$\displaystyle\Psi_{k|j}(t)=\frac{\mathbb{P}(S(t)\in\mathcal{S}^{k})}{\mathbb{P}(S(\infty)\in\mathcal{S}^{j})}$
(9)

where the numerator is the marginal occupation probability of occupying state
$k$ or any of the other states that come after $k$ in the progressive system
at time $t$, while the denominator is the marginal occupation probability of
occupying state $j$ or any of the other states that come after $j$ in the
system during an individual’s lifetime. We denote the resulting estimator
defined by the ratio of marginal estimates of state occupation probabilities
by $\widehat{\Psi}^{[2]}_{k|j}$. Since Equation (9) is a function of marginal
state occupation probabilities, we begin with estimating these marginal
probabilities by calculating the counting process and at-risk set using the
method of Datta and Sundaram (2006), described in Section 2.1. Afterward, we
obtain the Aalen-Johansen estimates of the marginal state occupation
probabilities (see details below). Finally, $\widehat{\Psi}^{[2]}_{k|j}$ is
obtained by the plug-in method.

The product-limit estimator of the marginal state occupation probability
relies on the Nelson-Aalen estimate of the cumulative transition intensity,
denoted as $\boldsymbol{A}(t)$. Let
$\widehat{\boldsymbol{A}}(t)=\\{\widehat{A}_{jj^{\prime}}(t)\\}$ denote the
corresponding estimator given by the expression:

$\displaystyle\widehat{A}_{jj^{\prime}}(t)=\begin{dcases}\int^{t}_{0}\mathcal{I}\Big{(}\widehat{Y}_{j}(u)>0\Big{)}\widehat{Y}_{j}(u)^{-1}\text{d}\widehat{N}_{jj^{\prime}}(u)&j\neq
j^{\prime},\\\ -\sum_{j^{\prime}\neq
j}\widehat{A}_{jj^{\prime}}(t)&j=j^{\prime},\end{dcases}$ (10)

where $\widehat{N}_{jj^{\prime}}(t)$ and $\widehat{Y}_{j}(t)$ respectively
denote the estimators of the counting process and the at-risk set (ref: Eqns
(1)-(2)). The product limit of $\widehat{\boldsymbol{A}}(t)$ leads to the
(marginal) Aalen-Johansen estimator of the transition probability matrix
${\boldsymbol{P}}(s,t)$,

$\displaystyle\widehat{\boldsymbol{P}}(s,t)=\prod_{(s,t]}\\{\boldsymbol{I}+\text{d}\widehat{\boldsymbol{A}}(u)\\},$
(11)

where $\boldsymbol{I}$ denotes the $Q\times Q$ identity matrix. Consequently,
the marginal estimator of the state occupation probability is given by

$\displaystyle\widehat{\pi}_{j}(t)=\sum^{Q}_{j^{\prime}=1}\widehat{\pi}_{j^{\prime}}(0)\widehat{P}_{j^{\prime}j}(0,t),$
(12)

where
$\widehat{\pi}_{j^{\prime}}(0)=\\{\widehat{Y}_{j^{\prime}}(0+)\\}/\\{\sum^{Q}_{j=1}\widehat{Y}_{j}(0+)\\}$,
and $\widehat{P}_{j^{\prime}j}(0,t)$ is the $j^{\prime}j$th element of
$\widehat{\boldsymbol{P}}(0,t)$. Datta and Satten (2001) demonstrated that
Equation (12) provides a consistent estimator for the state occupation
probability, even in non-Markov models. This finding allows us to estimate the
conditional state occupation probability without the Markov assumption using a
product-limit approach similar to the Aalen-Johansen calculation.

Once we construct $\widehat{\Psi}^{[2]}_{k|j}$ using the product-limit
approach, directly estimating the conditional entry distribution $F_{k|j}(t)$
becomes straightforward by substituting $\widehat{\Psi}^{[2]}_{k|j}$ into
Equation (8).

###  4.1 A toy example

We demonstrate the formulation of the new estimator using an illustrative
multistate system (refer to Figure 4) that can describe the progression of
Chronic Obstructive Pulmonary Disease (COPD). This system comprises seven
distinct states representing various health conditions individuals may
experience over time. Commencing from the healthy state (alive and event-
free), which mirrors normal lung function, individuals might either die
without any COPD events or transit to the mild COPD state upon the onset of
mild airflow limitation typical in early-stage COPD. Subsequently, from the
mild COPD state, individuals could either decease or progress to moderate COPD
as their symptoms exacerbate, and airflow limitation intensifies. Advancing
from the moderate COPD state, individuals could die or further deteriorate to
severe COPD, typified by notable limitations in physical activity, frequent
exacerbations, and a heightened risk of respiratory failure. Our focus lies in
estimating the proportion of individuals that will reach the severe COPD state
given a past occupation of the mild COPD state. This conditional probability
is denoted by $\Psi_{5|1}$. Additionally, we will describe the estimation of
$\Psi_{3|1}$ and $\Psi_{5|3}$.

The expression in (9) leads to the respective estimators for these temporal
functions:

$\displaystyle\widehat{\Psi}^{[2]}_{3|1}(t)$
$\displaystyle=\frac{\widehat{\pi}_{3^{*}}(t)}{\widehat{\pi}_{1^{*}}}=\frac{\widehat{\mathbb{P}}\Big{(}S(t)\in\\{3,5,6\\}\Big{)}}{\widehat{\mathbb{P}}\Big{(}S(\infty)\in\\{1,3,4,5,6\\}\Big{)}}$
$\displaystyle\widehat{\Psi}^{[2]}_{5|3}(t)$
$\displaystyle=\frac{\widehat{\pi}_{5^{*}}(t)}{\widehat{\pi}_{3^{*}}}=\frac{\widehat{\mathbb{P}}\Big{(}S(t)\in\\{5\\}\Big{)}}{\widehat{\mathbb{P}}\Big{(}S(\infty)\in\\{3,5,6\\}\Big{)}}$
$\displaystyle\widehat{\Psi}^{[2]}_{5|1}(t)$
$\displaystyle=\frac{\widehat{\pi}_{5^{*}}(t)}{\widehat{\pi}_{1^{*}}}=\frac{\widehat{\mathbb{P}}\Big{(}S(t)\in\\{5\\}\Big{)}}{\widehat{\mathbb{P}}\Big{(}S(\infty)\in\\{1,3,4,5,6\\}\Big{)}}.$

Let’s consider the estimators for ${\Psi}_{3|1}(t)$ and ${\Psi}_{5|1}(t)$. For
$\widehat{\Psi}_{3|1}(t)$, the numerator estimates the probability of
occupying the artificial state $3^{*}$ that comprises state 3 and other
subsequent states along its path, at time $t$. Whereas, for
$\widehat{\Psi}_{5|1}(t)$, the numerator estimates the occupation probability
of state 5 at time $t$. The denominator of the respective estimators has the
same form. Under the competing risk framework, the denominator simply
estimates the transition probability of moving from state 0 to state 1 at some
time during the observation period. In terms of state occupation
probabilities, the denominator is the probability of occupying the artificial
state $1^{*}$ that comprises state 1 and other subsequent states along its
path at some time during the observation period.

For a simple comparison of the new approach based on the ratio of product
limit estimators (PLE, hereinafter) with the approach based on the fractional
risk sets (FRE, hereinafter), Figure 3 shows the curves from the estimation of
the probability of occupying state 5 conditional on a previous visit of state
1, ${\Psi}_{5|1}(t)$, using simulated data. The 95% pointwise confidence
interval (see Section 5) corresponding to each probability curve is also shown
on the plots. Based on this simulated setting, the plots reveal the validity
of both methods in estimating the true probability function. Moreso, the
estimators are comparably similar. In Section 6, we carry out a thorough
investigation comparing the performance of both estimators.

Figure 3: The estimates and 95% pointwise bootstrap confidence intervals for
$\Psi_{5|1}$ based on simulated data from the seven-state illness death model
using the two proposed methods. The dark solid line represents the true
probability curve, the dark dotted line represents the estimated probability
and the grey dashed lines denote the confidence intervals.

##  5 Bootstrap-based Confidence Intervals

The utilization of nonparametric regression in conjunction with the PAV step
for formulating the estimators proposed in Section 3 and 4 imposes challenges
for the asymptotic analysis for these estimators. To quantify uncertainty, a
bootstrap-based procedure may provide a practical alternative. However,
Faraway and Jhun (1990) pointed out that a direct application of the bootstrap
will generally lead to invalid results for the “smoothed” nonparametric
estimators since it cannot capture the bias component of the smoothed
estimator.

Following guidelines from existing work on the bootstrap for nonparametric
regression estimators (Li and Datta, 2001; Shao and Tu, 1995), we describe a
smoothed bootstrap procedure for estimating the pointwise confidence intervals
for our estimators. First, using random sampling with replacement, we generate
bootstrap censored times $\\{C^{*}_{i};\ 1,...,n\\}$ from the original data
$\\{C_{i};\ 1,...,n\\}$. For the $i$th observation, we then generate the
corresponding state $S^{*}_{i}(C^{*}_{i})$ from the distribution given by

$\displaystyle\text{Prob}(S=j)=\frac{\sum^{n}_{l=1}\mathcal{I}(S_{l}(C_{l})=j)K_{\tilde{h}}(C_{l}-C^{*}_{i})}{\sum^{n}_{l=1}K_{\tilde{h}}(C_{l}-C^{*}_{i})}.$

In other words, we generate $S^{*}_{i}(C^{*}_{i})$ from a smoothed
distribution of the states whose inspection times were close to $C^{*}_{i}$.
Following the theory of the smoothed bootstrap, the bandwidth
$\tilde{h}=h^{\omega}$ needs to be larger than the bandwidth $h$ selected in
the estimation based on the original data. We recommend setting $\omega=0.8$
if $0<h<1$, and $\omega=1.2$ if $h>1$.

For $0<\alpha<1$, let $\widehat{\Delta}_{\alpha}(t)$ be the $(1-\alpha)$
bootstrap percentile of the distribution of

$\displaystyle\big{|}\sin^{-1}\big{(}\sqrt{\widehat{\theta}^{*}_{k|j}(t;h)}\big{)}-\sin^{-1}\big{(}\sqrt{\widehat{\theta}_{k|j}(t;\tilde{h})}\big{)}\big{|}$

where $\widehat{\theta}^{*}_{k|j}(t;h)$ is calculated using the bandwidth $h$
from the original estimation but is based on the bootstrap sample, while
$\widehat{\theta}_{k|j}(t;\tilde{h})$ uses the new bandwidth $\tilde{h}$ but
is based on the original sample. The $(1-\alpha)\times 100\%$ pointwise CI for
the temporal function $\theta_{k|j}(t)$ is given by

$\displaystyle\Big{[}\sin^{2}\Big{\\{}$
$\displaystyle\max\Big{(}0,\sin^{-1}\Big{(}\sqrt{\widehat{\theta}_{k|j}(t;h)}\Big{)}-\widehat{\Delta}_{1-\alpha/2}\Big{)}\Big{\\}},$
$\displaystyle\quad\sin^{2}\Big{\\{}\min\Big{(}\frac{\pi}{2},\sin^{-1}\Big{(}\sqrt{\widehat{\theta}_{k|j}(t;h)}\Big{)}+\widehat{\Delta}_{1-\alpha/2}\Big{)}\Big{\\}}\Big{]}.$

The construction of the confidence interval described above is similar to
Datta and Sundaram (2006). In the subsequent section, we conduct a Monte Carlo
study to examine the behavior of the smoothed bootstrap confidence intervals.

##  6 Simulation study

We performed extensive simulation studies to evaluate the performance of the
proposed methods based on: (i) fractional risk sets (FRE), and (ii) the ratio
of product limit estimators (PLE) in calculating the conditional estimators.
We considered two separate multistate systems - the five-state illness-death
model shown in Figure 1, and the seven-state multistate model shown in Figure
4. For each multistate system, the exact transition times were generated from
a common lognormal hazard. We considered sample sizes of 100, 200, 500, and
1000.

###  6.1 Five-State Illness-Death Model

Earlier, in Section 3.1, we introduced the five-state illness-death model. For
our simulation experiment, our primary interest lies in estimating the
probability of occupying state 3 having visited state 1, denoted by
$\Psi_{3|1}$. We also study the performance of the estimators for the
subdistribution function $\Psi_{3|1}(t)$, and the future entry distribution
for state 3 given the prior occupation of state 1, denoted by $F_{3|1}(t)$.
Our data simulation process initiates with each individual starting in state 0
at time 0, with potential transitions to either state 1 or state 2.
Specifically, we simulate the scenario where 60% of individuals follow the
$0\rightarrow 1$ transition path, while the remaining 40% follow the
$0\rightarrow 2$ transition path. The transition out of state 0 is governed by
an independent Bernoulli distribution. For those individuals transitioning to
state 1, we further simulate their progression, with 60% transitioning to
state 3 and 40% to state 4. This transition from state 1 is also controlled by
a separate independent Bernoulli distribution. Following this simulation
scheme, the true value of $\Psi_{3|1}$ is 0.6.

Let $V_{i,0}$ and $V_{i,1}$ denote the respective waiting time in states 0 and
1, for individual $i$. We simulate: $V_{i,0}\sim\text{LogNorm}(0,0.5)$ and
$V_{i,1}\sim\text{LogNorm}(0,0.5)$. We study the performance of our proposed
estimators under the scenarios where the random inspection time for each
individual is generated from (i) the uniform distribution with a lower bound
of 0 and an upper bound equal to the maximum waiting time observed across all
individuals in the sample, and (ii) the Weibull distribution with (shape)
$\eta=3$ and (scale) $\tau=2.5$.

###  6.2 Seven-State Multistate Model

Here, we consider the more complex multistate system shown in Figure 4.
Primarily, we are interested in estimating the probability of occupying State
5 having visited State 1, denoted by $\Psi_{5|1}$. From Eqn 6, we have
$\Psi_{5|1}=\Psi_{3|1}\times\Psi_{5|3}$. In addition to studying the
performance of the estimators for $\Psi_{5|1}$, we will also study the
performance of the estimators for $\Psi_{5|3}$, along with the subdistribution
functions; $\Psi_{5|1}(t)$ and $\Psi_{5|3}(t)$, and the future entry
distributions; $F_{5|1}(t)$, and $F_{5|3}(t)$. For the data simulation, all
individuals start in state 0 at time 0, and they have two possible paths for
transition: either to state 1 or state 2. We generate the data as follows: 80%
of the individuals follow the $0\rightarrow 1$ path, while the remaining 20%
follow the $0\rightarrow 2$ path. Among those who transitioned to state 1,
there are two possible paths: 70% transition to state 3, and the remaining 30%
transition to state 4. Moreover, among those that visited state 3, 60%
proceeded to state 5, while the remaining 40% transitioned to state 6. Note
that, the respective transitions out of state 0, state 1, and state 3 are
controlled by three separate Bernoulli distributions. Following this
simulation approach, the true value of $\Psi_{5|1}$ is 0.6.

To explore the performance of our proposed estimators, we simulate individual
waiting times in states 0 and 1 denoted by $V_{i,0}$ and $V_{i,1}$
respectively. We simulate: $V_{i,0}\sim\text{LogNorm}(0,0.5)$,
$V_{i,1}\sim\text{LogNorm}(0,0.5)$, and $V_{i,3}\sim\text{LogNorm}(0,0.7)$. We
consider two scenarios where we generate the random inspection time for each
individual from (i) the uniform distribution with limits 0 and the maximum
waiting time across the entire sample and, (ii) the Weibull distribution with
$\eta=2.5$ and $\tau=4.5$.

###  6.3 Simulation results

We evaluated the performance of the estimator of the conditional state
occupation probability $\Psi_{k|j}$ using absolute bias. For the global
assessment of the estimators of the subdistribution functions,
$\Psi_{k|j}(t)$, and the future state entry time distribution, $F_{k|j}(t)$,
we employed the mean absolute distance (MAD), defined as follows:

$\displaystyle\text{MAD}=\mathbb{E}\int|\widehat{\theta}^{\text{Emp}}(x)-\widehat{\theta}(x)|\
d\theta_{n}(x),$ (13)

where $\widehat{\theta}$ represents the estimator based on the current status
data, and $\widehat{\theta}^{\text{Emp}}$ corresponds to the empirical
estimator based on complete data - this served as the benchmark. To compute
the MAD, we used the empirical distribution function of the inspection times,
denoted by $\theta_{n}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathcal{I}(C_{i}\leq x)$.
When the MAD value is equal to 0, it indicates that the estimators are in
complete agreement on the support of the observed inspection times. We
calculated the absolute bias and MAD via Monte Carlo averaging with a
replication size of 1000.

For the five-state illness death model, Table 1 shows the results for the
estimators constructed with the two proposed methods under the simulation
setting where the true transition times follow a lognormal hazard. The results
are shown for the separate scenarios where the inspection time for each
individual is generated from the uniform and Weibull distribution. For each
estimator, the bias and MAD values decrease with increasing sample size, thus
suggesting the consistency of the estimators. Although the performance of both
methods in estimating the $\Psi_{3|1}$ is comparable in most cases, the
estimator based on the FRE approach tends to have a smaller bias. This finding
is consistent when we consider the other temporal functions of interest.
Generally, the performance of the estimators is reasonable considering the
severe censoring of the transition times.

Table 2 shows the simulation results for the seven-state multistate system
where the true transition times are generated from the lognormal hazard. For
each quantity of interest estimated by either the FRE or PLE approach, the
performance of the estimators appears to improve with increasing sample size.
Under this system, we observe that the estimators based on the FRE method
slightly outperform the estimators based on the PLE method. As pointed out in
Lan and Datta (2010), a possible explanation of the superiority of the FRE
method over the PLE method especially for later states in the system is that
the PLE-based estimation of the functions for the later states depends on the
estimation done in prior states, therefore, the estimation error may propagate
along the path.

Table 1: Absolute bias ($\Psi_{3|1}$) and MAD ($\Psi_{3|1}(t)$ & $F_{3|1}(t)$) for the scenario where the true transition times for the five-state illness-death model were generated from the lognormal distribution. Inspection time for each individual was generated from (i) uniform and (ii) Weibull distribution. The results are shown by the methods proposed: (a) FRE, and (b) PLE. The results are based on a Monte Carlo sample size of 1000. | $\Psi_{3|1}$ | $\Psi_{3|1}(t)$ | $F_{3|1}(t)$
---|---|---|---
Inspection time | $N$ | FRE | PLE | FRE | PLE | FRE | PLE
| 100 | 0.057 | 0.065 | 0.055 | 0.059 | 0.065 | 0.061
| 200 | 0.039 | 0.041 | 0.040 | 0.042 | 0.049 | 0.046
| 500 | 0.026 | 0.025 | 0.029 | 0.029 | 0.035 | 0.032
Uniform | 1000 | 0.019 | 0.018 | 0.023 | 0.022 | 0.027 | 0.025
| 100 | 0.047 | 0.060 | 0.053 | 0.065 | 0.082 | 0.087
| 200 | 0.036 | 0.045 | 0.041 | 0.047 | 0.070 | 0.072
| 500 | 0.028 | 0.030 | 0.029 | 0.032 | 0.054 | 0.054
Weibull | 1000 | 0.025 | 0.025 | 0.023 | 0.024 | 0.046 | 0.045
Table 2: Absolute bias ($\Psi_{5|1}\ \&\ \Psi_{5|3}$) and MAD
($\Psi_{5|1}(t),\ \Psi_{5|3}(t),\ F_{5|1}(t)\ \&\ F_{5|3}(t)$) for the
scenario where the true transition times for the seven-state illness-death
model were generated from the lognormal distribution. Inspection time for each
individual was generated from (i) uniform and (ii) Weibull distribution. The
results are shown by the methods proposed: (a) FRE, and (b) PLE. The results
are based on a Monte Carlo sample size of 1000.

| $\Psi_{5|1}$ | $\Psi_{5|3}$ | $\Psi_{5|1}(t)$ | $\Psi_{5|3}(t)$ | $F_{5|1}(t)$ | $F_{5|3}(t)$
---|---|---|---|---|---|---
Inspection time | $N$ | FRE | PLE | FRE | PLE | FRE | PLE | FRE | PLE | FRE | PLE | FRE | PLE
| 100 | 0.052 | 0.060 | 0.059 | 0.068 | 0.044 | 0.049 | 0.049 | 0.056 | 0.060 | 0.059 | 0.057 | 0.059
| 200 | 0.037 | 0.043 | 0.042 | 0.050 | 0.035 | 0.039 | 0.039 | 0.045 | 0.049 | 0.048 | 0.045 | 0.048
| 500 | 0.024 | 0.027 | 0.027 | 0.031 | 0.026 | 0.029 | 0.028 | 0.034 | 0.038 | 0.040 | 0.033 | 0.040
Uniform | 1000 | 0.018 | 0.019 | 0.021 | 0.022 | 0.021 | 0.025 | 0.022 | 0.029 | 0.030 | 0.035 | 0.025 | 0.035
| 100 | 0.054 | 0.056 | 0.068 | 0.073 | 0.045 | 0.046 | 0.054 | 0.058 | 0.074 | 0.072 | 0.072 | 0.072
| 200 | 0.036 | 0.044 | 0.040 | 0.048 | 0.034 | 0.034 | 0.038 | 0.040 | 0.060 | 0.060 | 0.058 | 0.060
| 500 | 0.027 | 0.031 | 0.030 | 0.034 | 0.026 | 0.023 | 0.028 | 0.028 | 0.049 | 0.048 | 0.047 | 0.048
Weibull | 1000 | 0.022 | 0.023 | 0.025 | 0.027 | 0.023 | 0.018 | 0.022 | 0.022 | 0.044 | 0.041 | 0.038 | 0.041

####  6.3.1 Coverage of smoothed bootstrap confidence intervals

For each multistate system considered in this simulation study, we evaluate
the performance of the smoothed bootstrap confidence intervals. To obtain the
confidence intervals, we used 1000 bootstrap resamples from the original data
to calculate the bootstrap percentiles. The reported empirical coverage
probabilities are based on 1000 Monte Carlo iterations at a nominal level of
0.05. Table 3 shows the empirical coverage probabilities for the estimators of
$\Psi_{3|1}(t)$ based on the five-state illness death model. The coverage
probabilities for the estimators of $\Psi_{5|1}(t)$ based on the seven-state
illness death model are shown in Table 4. The confidence intervals for the
temporal functions of interest are calculated at five fixed time points. The
results are shown for the scenario where the random inspection time for each
individual is generated from a uniform distribution. From Tables 3 and 4, the
coverage probabilities for the estimator based on both FRE and PLE approach
are reasonable.

Table 3: Coverage probabilities for the estimator of $\Psi_{3|1}(t)$ for the two proposed methods; (i) FRE, and (ii) PLE. The confidence intervals are based on 1000 bootstrap samples at a nominal level of 0.05, and the coverage results are based on 1000 Monte Carlo iterations. | FRE | PLE
---|---|---
| Sample size | Sample size
Time | 100 | 200 | 500 | 1000 | 100 | 200 | 500 | 1000
0.75 | 0.985 | 0.953 | 0.928 | 0.905 | 0.940 | 0.971 | 0.987 | 0.993
1.25 | 0.945 | 0.956 | 0.977 | 0.976 | 0.906 | 0.937 | 0.968 | 0.984
1.75 | 0.975 | 0.970 | 0.977 | 0.976 | 0.981 | 0.982 | 0.992 | 0.995
2.25 | 0.989 | 0.989 | 0.978 | 0.989 | 0.995 | 0.996 | 0.997 | 0.999
2.75 | 0.996 | 0.995 | 0.997 | 0.985 | 0.999 | 0.999 | 0.997 | 0.999
Table 4: Coverage probabilities for the estimator of $\Psi_{5|1}(t)$ for the two proposed methods; (i) FRE, and (ii) PLE. The confidence intervals are based on 1000 bootstrap samples at a nominal level of 0.05, and the coverage results are based on 1000 Monte Carlo iterations. | FRE | PLE
---|---|---
| Sample size | Sample size
Time | 100 | 200 | 500 | 1000 | 100 | 200 | 500 | 1000
2.00 | 0.858 | 0.873 | 0.840 | 0.813 | 0.929 | 0.955 | 0.979 | 0.985
3.00 | 0.853 | 0.851 | 0.797 | 0.764 | 0.976 | 0.983 | 0.988 | 0.968
4.00 | 0.970 | 0.957 | 0.945 | 0.935 | 0.994 | 0.992 | 0.996 | 0.978
5.00 | 0.994 | 0.994 | 0.986 | 0.984 | 0.998 | 0.999 | 0.997 | 0.994
6.00 | 0.999 | 1.000 | 0.997 | 0.998 | 0.999 | 0.999 | 0.999 | 1.000

##  7 Breast cancer example (EORTC-trial 10854)

To further illustrate the two proposed estimation methods, we re-analyzed data
on 2793 breast cancer patients who participated in a trial conducted by the
European Organization for Research and Treatment of Cancer (EORTC-trial
10854). The trial which involved patients with early breast cancer, aimed to
compare whether surgery followed by a short intensive course of perioperative
chemotherapy yields better treatment than surgery alone. A comprehensive
description of the study is provided in Clahsen et al. (1996), and further
findings are reported in Van der Hage et al. (2001).

Figure 4 describes the multistate system employed to illustrate our proposed
methods. This system comprises nine states, including event-free and alive
after surgery (state 0), loco-regional relapse following surgery (state 1),
distant metastasis following surgery (state 2), simultaneous occurrence of
loco-regional relapse and distinct metastasis following surgery (state 3),
death without the occurrence of loco-regional relapse or distinct metastasis
following surgery (state 4), distant metastasis following occurrence of loco-
regional relapse after surgery (state 5), death following occurrence of loco-
regional relapse after surgery (state 6), loco-regional relapse following
occurrence of distant metastasis after surgery (state 7), and death following
occurrence of distant metastasis after surgery (state 8). The time origin for
the analysis was set as the time of surgery, meaning all individuals commenced
at state 0. Transition times in the original dataset were right-censored, with
follow-up times ranging from 12 to 6,756 days, and a median of 3,817 days.

To obtain current status data, the inspection time for each patient was
randomly generated from a uniform distribution with a lower bound of 0 and an
upper bound of 6,756 days, which is the maximum follow-up time observed across
all patients in the original data. The inspection times were independent of
the progression between the event states. We defined a patient’s inspection
state as the state occupied at their inspection time. If the patient was
right-censored, the inspection state was the most recent state they occupied
before censoring. To highlight the structure of the current status data from
the multistate system, we provide information for five randomly selected
individuals in the supplementary material. This cross-sectional version of the
data poses a severe type of censoring since we would only know the state
occupied by each patient at a random inspection time rather than their
transition history up to a censoring time.

0\. Surgery 1\. Local recurrence only 2\. Distant metastasis only 3\. Local
recurrence and distant metastasis 4\. Death without local recurrence or
distant metastasis 5\. Distant metastasis following local recurrence 6\. Death
following Local recurrence 7\. Local recurrence following distant metastasis
8\. Death following distant metastasis Figure 4: Multistate system for the
breast cancer study from EORTC-trial 10854.

Table 5 summarizes the distribution of the state occupation at the inspection
times. As expected, most patients were inspected while occupying state 0
(alive after surgery). We are interested in estimating the probability of a
patient ever reaching state 5 given a previous occupation of state 1.
Inference around this quantity is crucial in breast cancer studies (Vicini et
al., 2003; Murata et al., 2023) since identifying patients at higher risk of
distinct metastasis following local-regional relapse informs better treatment
decisions to prevent or manage metastatic disease. This conditional quantity,
denoted by $\Psi_{5|1}$, was estimated using both the FRE and PLE approaches.

The estimates and corresponding 95% bootstrap confidence intervals based on
the FRE and PLE approach are $0.400\ (0.241,0.570)$ and $0.433\
(0.222,0.658)$, respectively. The estimation results from both proposed
methods are comparable. To provide more context to the results, we conducted a
purely marginal analysis (without conditioning on a previous visit to state 1)
where we estimated the probability of ever reaching state 5. The resulting
estimate is $0.050$. As expected, conditioning on the previous occupation of
the relapse state significantly increases the probability of ever reaching
state 5.

We further estimated $\Psi_{5|1}$ using the PLE approach on the original
right-censored data. This yielded an estimate of $0.344$. Notably, the
estimations based on the heavily censored cross-sectional data approximate the
result obtained from the right-censored data, suggesting their potential as an
alternative analysis approach.

Table 5: Distribution of state occupation at inspection based on the current status data from the EORTC-trial 10854. State | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
Frequency | 1941 | 119 | 122 | 66 | 93 | 57 | 29 | 30 | 336
Percent | 0.695 | 0.043 | 0.044 | 0.024 | 0.033 | 0.020 | 0.010 | 0.011 | 0.120

##  8 Discussion

This paper carefully employs techniques from the framework of competing risks
to construct nonparametric estimators of future state entry probabilities and
future state entry distributions, predicated on the prior occupation of a
state, utilizing current status multistate data. The cross-sectional nature of
the censored event history data brings about additional complexities when
attempting to construct nonparametric estimates of temporal functions within a
multistate model. We tackle this problem by harnessing ideas from
nonparametric regression modeling to formulate estimators of the counting
process and at-risk sets that serve as foundational elements in constructing
estimators of the probability functions. To estimate the conditional
quantities in the context of current status data, we employed a method (FRE)
that was only available for right-censored data, and we also introduced a
novel method (based on PLE). Using data from a real-world application and
simulation studies for illustrative purposes, we show that both estimation
methods yield valid and comparable results.

In practice, each estimation method has its benefits and challenges.
Constructing the estimator based on the FRE approach could appear challenging
since it requires the calculation of the fractional weights. Although the PLE
method may be easier to implement once the artificial states have been
specified correctly, the FRE method may be preferred for calculations
involving distant states since it produces more reliable estimates for such
states.

The intuitive approaches we propose provide a more direct route for analyzing
the severely limited event data since they utilize the averaging principle
rather than maximizing a likelihood function. Suppose we consider a less
complicated setting (survival model with two states); the nonparametric
maximum likelihood estimator (NPMLE) of the survival function has undesirable
properties since it is undefined over a set of regions. For the general
multistate model with current status data, an NPMLE approach would be more
problematic.

We operate under the assumption of a progressive multistate system, structured
as a directed tree, where each state is reachable via a unique path
originating from the root node. This implies knowledge of path information
from past states at a random inspection time, although the exact transition
times remain unknown. While this assumption often holds due to the inherent
structure of the multistate system, there are instances where distinctions
among paths can be made by incorporating auxiliary states into an expanded
multistate system. Additionally, a considerable number of inspections in a
specific state are necessary to calculate the occupation probability of that
state. In general, a relatively large sample size is required for valid
nonparametric estimation of occupation probabilities, especially when dealing
with current status data, where information is significantly limited compared
to complete or right-censored data.

We could extend the proposed methods to multistate models with case-II
interval-censored data by utilizing similar formulas for estimating the
counting process and the at-risk sets. Since the estimators of these
functionals are based on averaging, we could sum over the multiple inspection
times for each individual. Ignoring the correlation among the inspection times
will not inhibit the unbiasedness and consistency properties of the resulting
estimators. However, the efficiency of the estimators could be improved by
including a weighted variance-covariance matrix in the PAV step of the
isotonic regression to account for the events recorded multiple times for each
individual.

Further, researchers analyzing current status multistate data may also be
interested in assessing the impact of measured risk factors on the probability
of occupying a future state given a prior visit to a particular state.
Regression modeling based on jackknife pseudo-values (Andersen et al., 2003)
presents a flexible approach for estimating covariate effects in event history
analysis. Recently, Anyaso-Samuel et al. (2023) extended the pseudo-value
approach to model the effect of covariates on state occupation probabilities
with multistate current status data. By jackknifing any of the estimators
proposed here and utilizing the guidelines provided in Anyaso-Samuel et al.
(2023), researchers can easily conduct covariate inference for the future
state entry probabilities.

## Data Availability Statement

The data from the breast cancer study used in Section 7 are available from the
European Organization for Research and Treatment of Cancer (EORTC).
Restrictions apply to the availability of these data, which were used under
license in this paper.

## References

  * Aalen and Johansen (1978) Aalen, O. O. and S. Johansen (1978).  An empirical transition matrix for non-homogeneous Markov chains based on censored observations.  Scandinavian Journal of Statistics 5, 141–150.
  * Andersen et al. (1993) Andersen, P. K., O. Borgan, R. D. Gill, and N. Keiding (1993).  Statistical models based on counting processes.  New York: Springer-Verlag.
  * Andersen et al. (2003) Andersen, P. K., J. P. Klein, and S. Rosthøj (2003).  Generalised linear models for correlated pseudo-observations, with applications to multi-state models.  Biometrika 90(1), 15–27.
  * Anyaso-Samuel et al. (2023) Anyaso-Samuel, S., D. Bandyopadhyay, and S. Datta (2023).  Pseudo-value regression of clustered multistate current status data with informative cluster sizes.  Statistical Methods in Medical Research 32(8), 1494–1510.
  * Barlow et al. (1972) Barlow, R., D. Bartholomew, J. Bremner, and H. Brunk (1972).  Statistical Inference Under Order Restrictions: The Theory and Application of Isotonic Regression.  New York, N.Y.: J. Wiley.
  * Clahsen et al. (1996) Clahsen, P. C., C. Van de Velde, J.-P. Julien, J.-L. Floiras, T. Delozier, F. Y. Mignolet, and T. M. Sahmoud (1996).  Improved local control and disease-free survival after perioperative chemotherapy for early-stage breast cancer. a european organization for research and treatment of cancer breast cancer cooperative group study.  Journal of clinical oncology 14(3), 745–753.
  * Datta and Satten (2000) Datta, S. and G. A. Satten (2000).  Estimating future stage entry and occupation probabilities in a multistage model based on randomly right-censored data.  Statistics & probability letters 50(1), 89–95.
  * Datta and Satten (2001) Datta, S. and G. A. Satten (2001).  Validity of the Aalen–Johansen estimators of stage occupation probabilities and Nelson–Aalen estimators of integrated transition hazards for non-Markov models.  Statistics & probability letters 55(4), 403–411.
  * Datta and Sundaram (2006) Datta, S. and R. Sundaram (2006).  Nonparametric estimation of stage occupation probabilities in a multistage model with current status data.  Biometrics 62(3), 829–837.
  * Faraway and Jhun (1990) Faraway, J. J. and M. Jhun (1990).  Bootstrap choice of bandwidth for density estimation.  Journal of the American Statistical Association 85(412), 1119–1122.
  * Jewell and van der Laan (2003) Jewell, N. P. and M. van der Laan (2003).  Current status data: Review, recent developments and open problems.  Handbook of statistics 23, 625–642.
  * Lan et al. (2017) Lan, L., D. Bandyopadhyay, and S. Datta (2017).  Non-parametric regression in clustered multistate current status data with informative cluster size.  Statistica Neerlandica 71(1), 31–57.
  * Lan and Datta (2010) Lan, L. and S. Datta (2010).  Non-parametric estimation of state occupation, entry and exit times with multistate current status data.  Statistical Methods in Medical Research 19(2), 147–165.
  * Li and Datta (2001) Li, G. and S. Datta (2001).  A bootstrap approach to nonparametric regression for right censored data.  Annals of the Institute of Statistical Mathematics 53, 708–729.
  * Mukerjee (1988) Mukerjee, H. (1988).  Monotone nonparametric regression.  The Annals of Statistics 16, 741–750.
  * Murata et al. (2023) Murata, T., M. Yoshida, S. Shiino, A. Ogawa, C. Watase, K. Satomi, K. Jimbo, A. Maeshima, E. Iwamoto, S. Takayama, et al. (2023).  A prediction model for distant metastasis after isolated locoregional recurrence of breast cancer.  Breast Cancer Research and Treatment 199(1), 57–66.
  * Nadaraya (1964) Nadaraya, E. A. (1964).  On estimating regression.  Theory of Probability & Its Applications 9(1), 141–142.
  * Shao and Tu (1995) Shao, J. and D. Tu (1995).  The Jackknife and Bootstrap.  Springer Series in Statistics. New York: Springer.
  * Sun (2006) Sun, J. (2006).  The Statistical Analysis of Interval-Censored Failure Time Data.  New York: Springer.
  * Van der Hage et al. (2001) Van der Hage, J., C. van De Velde, J.-P. Julien, J.-L. Floiras, T. Delozier, C. Vandervelden, L. Duchateau, et al. (2001).  Improved survival after one course of perioperative chemotherapy in early breast cancer patients: long-term results from the european organization for research and treatment of cancer (eortc) trial 10854.  European journal of cancer 37(17), 2184–2193.
  * Vicini et al. (2003) Vicini, F. A., L. Kestin, R. Huang, and A. Martinez (2003).  Does local recurrence affect the rate of distant metastases and survival in patients with early-stage breast carcinoma treated with breast-conserving therapy?  Cancer: Interdisciplinary International Journal of the American Cancer Society 97(4), 910–919.
  * Wand and Jones (1995) Wand, M. P. and M. C. Jones (1995).  Kernel smoothing.  London: Chapman & Hall.
  * Watson (1964) Watson, G. S. (1964).  Smooth regression analysis.  Sankhyā: The Indian Journal of Statistics, Series A 26(4), 359–372.
  * Yang et al. (2023) Yang, Y., S. Wu, and S. Datta (2023).  Regression analysis of a future state entry time distribution conditional on a past state occupation in a progressive multistate model.  Statistical Methods in Medical Research 32(12), 2285–2298.