# Joint Dynamic Models and Statistical Inference for Recurrent Competing
Risks, Longitudinal Marker, and Health Status

Lili Tong
Department of Biostatistics
University of Nebraska Medical Center
Omaha, NE 68198
<EMAIL_ADDRESS>
&Piaomu Liu
Department of Mathematical Sciences
Bentley University
Waltham, MA 02452
<EMAIL_ADDRESS>
&Edsel A. Peña
Department of Statistics
University of South Carolina
Columbia, SC 29208
<EMAIL_ADDRESS>

###### Abstract

Consider a subject or unit in a longitudinal biomedical, public health,
engineering, economic or social science study which is being monitored over a
possibly random duration. Over time this unit experiences recurrent events of
several types and a longitudinal marker transitions over a discrete state-
space. In addition, its “health” status also transitions over a discrete
state-space with at least one absorbing state. A vector of covariates will
also be associated with this unit. Of major interest for this unit is the
time-to-absorption of its health status process, which could be viewed as the
unit’s lifetime. Aside from being affected by its covariate vector, there
could be associations among the recurrent competing risks processes, the
longitudinal marker process, and the health status process in the sense that
the time-evolution of each process is associated with the other processes. To
obtain more realistic models and enhance inferential performance, a joint
dynamic stochastic model for these components is proposed and statistical
inference methods are developed. This joint model, formulated via counting
processes and continuous-time Markov chains, has the potential of facilitating
‘personalized’ interventions. This could enhance, for example, the
implementation and adoption of precision medicine in medical settings. Semi-
parametric and likelihood-based inferential methods for the model parameters
are developed when a sample of these units is available. Finite-sample and
asymptotic properties of estimators of model parameters, both finite- and
infinite-dimensional, are obtained analytically or through simulation studies.
The developed procedures are illustrated using a real data set.

_K_ eywords Continuous-time Markov chain $\cdot$ Counting process $\cdot$
Dynamic models $\cdot$ Intensity-based model $\cdot$ Personalized medicine
$\cdot$ Conditional independence $\cdot$ Parametric and semi-parametric
estimation

MSC2020 Subject Classifications: Primary: 62N01, 60J27; Secondary: 62G05,
62M02.

##  1 Introduction and Motivation

Consider a unit — for example, a human subject in a medical or a social
science study, an experimental animal in a biological experiment, a machine in
an engineering or reliability setting, or a company in an economic or business
situation — in a longitudinal study monitored over a period $[0,\tau$], where
$\tau$ is possibly random. Associated with the unit is a covariate row vector
$\mathbf{X}=(X_{1},X_{2},\ldots,X_{p})$. Over time, the unit will experience
occurrences of $Q$ competing types of recurrent events, its recurrent
competing risks (RCR) component; transitions of a longitudinal marker (LM)
process $W(t)$ over a discrete state space $\mathfrak{W}$; and transitions of
a ‘health’ status (HS) process $V(t)$ over a discrete state space
$\mathfrak{V}=\mathfrak{V}_{0}\bigcup\mathfrak{V}_{1}$, with
$\mathfrak{V}_{0}\neq\emptyset$ being absorbing states. If the health status
process transitions into an absorbing state prior to $\tau$, then monitoring
of the unit ceases, so time-to-absorption serves as the lifetime of the unit.
To demonstrate pictorially, the two panels in Figure 1 depict the time-
evolution for two distinct units, where $Q=3$,
$\mathfrak{W}=\\{w_{1},w_{2},w_{3}\\}$, $\mathfrak{V}=\\{v_{0},v_{1},v_{2}\\}$
with $v_{0}$ an absorbing state. In panel 1, the unit did not transition to an
absorbing state prior to reaching $\tau$; whereas in panel 2, the unit reached
an absorbing state prior to $\tau$. Two major questions arise: (a) how do we
specify a dynamic stochastic model that could be a generative model for such a
data, and (b) how do we make statistical inferences for the model parameters
and predictions of a unit’s lifetime from a sample of such data?

|
---|---
Figure 1: Realized data observables from two distinct study units. The first
plot panel is for a unit which was right-censored prior to reaching an
absorbing state, while the second plot panel is for a unit which reached an
absorbing state prior to getting right-censored.

To address these questions, the major goals of this paper are (i) to propose a
joint stochastic model for the random observables consisting of the RCR, LM,
and HS components for such units, and (ii) to develop appropriate statistical
inference methods for the proposed joint model when a sample of units are
observed. Achieving these two goals will enable statistical prediction of the
(remaining) lifetime of a, possibly new, unit; allow for the examination of
the synergistic association among the RCR, LM, and HS components; and provide
a vehicle to compare different groups of units and/or study the effects of
concomitant variables or factors. More importantly, a joint stochastic model
endowed with proper statistical inference methods could potentially enable
unit-based interventions which are performed after a recurrent event
occurrence or a transition in either the LM or HS processes. As such it could
enhance the implementation of precision or personalized decision-making; for
instance, precision medicine in a medical setting.

A specific situation where such a data accrual occurs is in a medical study.
For example, a subject may experience different types of recurring cancer,
with the longitudinal marker being the prostate-specific antigen (PSA) level
categorized into a finite ordinal set, while the health status is categorized
into either a healthy, diseased, or dead state, with the last state absorbing.
A variety of situations in a biomedical, engineering, public health,
sociology, and economics settings, where such data structure arise, are
further described in Section 2. Several previous works dealing with modeling
have either focused in the marginal modeling of each of the three data
components, or in the joint modeling of two of the three data components. In
this paper we tackle the problem of simultaneously modeling all three data
components: RCR, LM, and HS, in order to account for the associations among
these components, which would not be possible using either the marginal
modeling approach or the joint modeling of pairwise combinations of these
three components. A joint full model could also subsume these previous
marginal or joint models – in fact, our proposed class of models subsumes as
special cases models that have been considered in the literature. In contrast,
only by imposing restrictive assumptions, such as the independence of the
three model components, could one obtain a joint full model from marginal or
pairwise joint models. As such, a joint full model will be less likely to be
mis-specified, thereby reducing the potential biases that could accrue from
mis-specified models among estimators of model parameters or when predicting
residual lifetime.

A joint modeling approach has been extensively employed in previous works. For
instance, joint models for an LM process and a survival or time-to-event (TE)
process have been proposed in [29], [32], [27], [16], [30], and [22]. Also,
the joint modeling of an LM process and a recurrent event process has also
been discussed in [15] and [11], while the joint modeling of a recurrent event
and a lifetime has also been done such as in [19]. An important and critical
theoretical aspect that could not be ignored in these settings is that when an
event occurrence is terminal (e.g., death) or when there is a finite
monitoring period, informative censoring naturally occurs in the RCR, LM, or
HS components, since when a terminal event or when the end of monitoring is
reached, then the next recurrent event occurrence, the next LM transition, or
the next HS transition will not be observed.

Another aspect that needs to be taken into account in a dynamic modeling
approach is that of performed interventions, usually upon the occurrence of a
recurrent event. For instance, in engineering or reliability systems, when a
component in the system fails, this component will either be minimally or
partially repaired, or altogether replaced with a new component (called a
perfect repair); while, with human subjects, when a cancer relapses, a
hospital infection transpires, a gout flare-up, or alcoholism recurs, some
form of intervention will be performed or instituted. Such interventions will
impact the time to the next occurrence of the event, hence it is critical that
such intervention effects be reflected in the model; see, for instance, [14]
and [15]. In addition, models should take into consideration the impacts of
the covariates and the effects of accumulating event occurrences on the unit.
Models that take into account these considerations have been studied in [23]
and [25]. Appropriate statistical inference procedures for these dynamic
models of recurrent events and competing risks have been developed in [23] and
[28]. Extensions of these joint dynamic models for both RCR and TE can be
found in [20]. Some other recent works in joint modeling included the modeling
of the three processes: LM, RCR (mostly, a single recurrent event), and TE
simultaneously in [17], [7], [18], [21] and [5]. The joint model that will be
proposed in this paper will take into consideration these important aspects.

We now outline the remainder of this paper. Prior to describing formally the
joint model in Section 3, we first present in Section 2 some concrete
situations in science, medicine, engineering, social, and economic disciplines
where the data accrual described above could arise and where the joint model
will be relevant. Section 3 formally describes the joint model using counting
processes and continuous-time Markov chains (CTMCs), and provide
interpretations of model parameters. In subsection 3.3 we discuss in some
detail a special case of this joint model obtained using independent Poisson
processes and homogeneous CTMCs. Section 4 deals with the estimation of the
parameters. In subsection 4.1 we demonstrate the estimation of the model
parameters in the afore-mentioned special case to pinpoint some intricacies of
joint modeling and its inferential aspects. The general joint model contains
nonparametric (infinite-dimensional) parameters, so in subsection 4.2 we will
describe a semi-parametric estimation procedure for this general model.
Section 5 will present asymptotic properties of the estimators, though we will
not present in this paper the rigorous proofs of these analytical results but
defer them to a more theoretically-focused paper. Section 6 will then
demonstrate the semi-parametric estimation approach through the use of a
synthetic or simulated data using an R [26] program we developed. In Section
7, we perform simulation studies to investigate the finite-sample properties
of the estimators arising from semi-parametric estimation procedure and
compare these finite-sample results to the theoretical asymptotic results. An
illustration of the semi-parametric inference procedure using a real medical
data set is presented in Section 8. Section 9 contains concluding remarks
describing some open research problems.

##  2 Concrete Situations of Relevance

To demonstrate potential applicability of the proposed joint model, we
describe in this section some concrete situations arising in biomedical,
reliability, engineering, and socio-economic settings where the data accrual
described in Section 1 could arise.

  * •

A Medical Example: Gout is a form of arthritis characterized by sudden and
severe attacks of pain, swelling, redness and tenderness in one of more joints
in the toes, ankles, knees, elbows, wrists, and fingers (see, for instance,
Mayo Clinic publications about gout). When a gout flare occurs, it renders the
person incapacitated (personally attested by the senior author) and the
debilitating condition may last for several days. Since the location of the
gout flare could vary, we may consider gout as competing recurrent events —
competing with respect to the location of the flare, and recurrent since it
could keep coming back. Gout occurs when urate crystals accumulate in the
joints, which in turn is associated with high levels of uric acid in the
blood. The level of uric acid is measured by the Serum Urate Level (SUR),
which can be categorized as Hyperuricemia (if SUR $>$ 7.2 mg/dL for males; if
SUR > 6.0 mg/dL for females), or Normal (if 3.5 mg/dL $\leq$ SUR $\leq$ 7.2
mg/dL for males; if 2.6 mg/dL $\leq$ SUR $\leq$ 6.0 mg/dL for females). The
SUR level could be considered a longitudinal marker. Kidneys are associated
with the excretion of uric acid in the body. Thus, chronic kidney disease
(CKD) impacts the level of uric acid in the body, hence the occurrence of
gout. The ordinal stages of CKD, based on the value of the Glomerular
Filtration Rate (GFR), are as follows: Stage 1 (Normal) if GFR $\geq$ 90
mL/min; Stage 2 (Mild CKD) if 60 mL/min $\leq$ GFR $\leq$ 89 mL/min; Stage 3A
(Moderate CKD) if 45 mL/min $\leq$ GFR $\leq$ 59 mL/min; Stage 3B (Moderate
CKD) if 30 mL/min $\leq$ GFR $\leq$ 44; Stage 4 (Severe CKD) if 15 mL/min
$\leq$ GFR $\leq$ 29 mL/min; and Stage 5 (End Stage CKD) if GRF $\leq$ 14
mL/min. The state of Stage 5 (End Stage CKD) can be viewed as an absorbing
state. The CKD status could be viewed as the “health status” of the person.
Other covariates, such as gender, blood pressure, weight, etc., could also
impact the occurrence of gout flares, uric acid level, and CKD. When a gout
flare occurs, lifestyle interventions could be performed such as (i) consuming
skim milk powder enriched with the two dairy products glycomacropeptide (GMP)
and G600 milk fat extract; or (ii) consuming standard milk or lactose powder.
The purpose of such interventions is to lessen gout flare recurrences. Of
major interest is to jointly model the competing gout recurrences, the
categorized SUR process, and the CKD process. A study consisting of $n$
subjects could be performed with each subject monitored over some period, with
the time of gout flare recurrences of each type, SUR levels, and CKD states
recorded over time, aside from relevant covariates. Based on such a data, it
will be of interest to estimate the model parameters and to develop a
prediction procedure for time-to-absorption to End Stage CKD for a person with
gout recurrences.

  * •

A Reliability Example: Observe $n$ independent cars over each of their own
monitoring period until the car is declared inoperable or the monitoring
period ends. Cars are complex systems in that they are constituted by
different components, which could be subsystems or modules, configured
according to some coherent structure function [3]. For each car, the states of
$Q$ components (such as its engine subsystem; transmission subsystem; brake
subsystem; electrical subsystem; etc.) are monitored. Furthermore, its
covariates such as weight; initial mileage; current mileage; years of
operation; and other characteristics (for example, climate in which car is
mostly being driven) are observed. Also, its ‘health status’, which is either
functioning, functioning with some problems, or total inoperability (an
absorbing state), is tracked over the monitoring period. Meanwhile, a
longitudinal marker such as its oil quality indicator (which is either
excellent; good; or poor) and the occurrences of failures of any of the $Q$
components are also recorded over the monitoring period. When a component
failure occurs, a repair or replacement of the component is undertaken. Given
the data for these $n$ cars, an important goal is to predict the time to
inoperability of another car. Note that this type of application could occur
in more complex systems such as space satellites, nuclear power plants,
medical equipments, etc.

  * •

A Social Science Example: Observe $n$ independent newly-married couples over a
period of years (say, 20 years). Over this follow-up period, the marriage
could end in separation or divorce, remain intact, or end due to the death of
at least one of them. Each couple will have certain characteristics: their
ages when they got married; working status of each; income level of each;
education level of each; number of children (this is a time-dependent
covariate); net worth of couple; etc. Possibly through a regularly
administered questionnaire, the couple provides information from which their
“marriage status” could be ascertained (either very satisfied; satisfied;
poor; separated or divorced). Competing recurrent events for the couple could
be changes in job status for either; addition in the family; educational
changes of children; and major disagreements. A longitudinal marker could be
the financial health status of the couple reflected by their categorized FICO
scores. A goal is to infer about the parameters of the joint model based on
the observations from these $n$ couples, and to predict if separation or
divorce will occur for a married couple, who are not in the original sample,
and if so, to obtain a prediction interval of the time of such an occurrence.

  * •

A Financial Example: Track $n$ independent publicly-listed companies over
their monitoring periods. At the end of its monitoring period, a company could
be bankrupt, still in business, or could have been acquired by another
company. Each company has its own characteristics, such as total assets,
number of employees, number of branches, etc. Note that these are all time-
dependent characteristics. The “health status” of a company is rated according
to four categories (A: Exceptional; B: Less than Exceptional; C: Distressed;
D: Bankrupt). The bankrupt status is the absorbing state. The company’s
liability relative to its asset, categorized into High, Medium, Low, Non-
Existent could be an important longitudinal marker. Recurrent competing risks
will be the occurrence of an increase (decrease) of at least 5% during a
trading day in its stock share price. Based on data from a sample of these
companies, it could be of interest to predict the time to bankruptcy of
another company that is not in the sample.

  * •

COVID-19 Example: Consider a vaccine trial where $n$ subjects are randomized
into different vaccine groups, including a no vaccine group. Group membership
could be coded using dummy covariates. Each subject is monitored over an
observation period, until loss to follow-up, or until death. Aside from the
group membership, other covariates (e.g., age or age-group, gender, BMI, race,
pre-existing conditions such whether immuno-compromised or not, political
affiliation, religious affiliation, etc.) for each subject are also observed.
A longitudinal marker for each subject could be the viral load, categorized
into none, low, medium, or high. See [31] for other examples of longitudinal
medical markers observed in Covid-19 studies. The health status for each
subject could be classified into free of Covid-19, moderately infected,
severely infected, or dead. Competing recurrent events could be the occurrence
of abdominal problems, severe coughing, body temperature reaching 103 degrees
Fahrenheit. Possible goals of the study are to compare the different vaccine
groups in terms of preventing Covid-19 infection; with respect to the mean or
median time to absorption; or with respect to mean or median time to
transitioning out of infection state given Covid-19 infection.

##  3 Joint Model of RCR, LM, and HS Processes

###  3.1 Data Observables for One Unit

Denote by $(\Omega,\mathfrak{F},\mathbf{P})$ the basic filtered probability
space with filtration $\mathcal{F}=\\{\mathcal{F}_{s}:s\geq 0\\}$ where all
random entities under consideration are defined. We begin by describing the
joint model for the data observable components for one unit.

Let $\tau$, the end of monitoring period, have a distribution function
$G(\cdot)$, which may be degenerate. The covariate vector will be
${X}=(X_{1},\ldots,X_{p})$, assumed to be time-independent, though the
extension to time-dependent covariates are possible with additional
assumptions. For the RCR component, let
$N^{R}=\\{N_{q}^{R}(s)\equiv(N^{R}(s;q),q\in\mathfrak{I}_{Q}):\ s\geq 0\\}$,
with index set $\mathfrak{I}_{Q}=\\{1,\ldots,Q\\}$, be a $Q$-dimensional
multivariate counting (column) vector process such that, for
$q\in\mathfrak{I}_{Q}$, $N_{q}^{R}(s)$ is the number of observed occurrences
of the recurrent event of type $q$ over $[0,s]$, with $N_{q}^{R}(0)=0$. Thus,
the sample path $s\mapsto N_{q}^{R}(s)$ takes values in
$\mathbb{Z}_{0,+}=\\{0,1,2,\ldots\\}$, is a non-decreasing step-function, and
is right-continuous with left-hand limits. We denote by
$dN_{q}^{R}(s)=N_{q}^{R}(s)-N_{q}^{R}(s-)$, the jump at time $s$ of
$N_{q}^{R}$.

For the LM process, let $W=\\{W(s):s\geq 0\\}$, where $W(s)$ takes values in a
finite state space $\mathfrak{W}$ with cardinality $|\mathfrak{W}|$. $W(s)$
represents the state of the longitudinal marker at time $s$. The sample path
$s\mapsto W(s)$ is a step-function which is right-continuous with left-hand
limits. By introducing the index set
$\mathfrak{I}_{\mathfrak{W}}=\\{(w_{1},w_{2}):w_{1},w_{2}\in\mathfrak{W},w_{1}\neq
w_{2}\\}$, we can convert $W$ into a (column)
$2{\binom{|\mathfrak{W}|}{2}}$-dimensional multivariate counting process
$\\{N^{W}\equiv(N^{W}(s;w_{1},w_{2}),(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}):s\geq
0\\}$, where $N^{W}(s;w_{1},w_{2})$ is the number of observed transitions of
$W$ from state $w_{1}$ into $w_{2}$ over the period $[0,s]$, that is, for
$(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}$,

$N^{W}(s;w_{1},w_{2})=\sum_{t\leq s}I\\{W(t-)=w_{1},W(t)=w_{2}\\},$

with $I(\cdot)$ denoting indicator function. Thus, for each
$(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}$, the sample path $s\mapsto
N^{W}(s;w_{1},w_{2})$ takes values in $\mathbb{Z}_{0,+}$, is a nondecreasing
step-function, right-continuous with left-hand limits, and with
$N^{W}(s;w_{1},w_{2})=0$. Furthermore,
$\sum_{(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}}dN^{W}(s;w_{1},w_{2})\in\\{0,1\\}$
for every $s\geq 0$.

For the HS process, let $V=\\{V(s):s\geq 0\\}$, where $V(s)$ takes values in
the finite state space $\mathfrak{V}=\mathfrak{V}_{0}\bigcup\mathfrak{V}_{1}$,
where states in $\mathfrak{V}_{0}$ are absorbing states, and with
$|\mathfrak{V}_{0}|>0$. $V(s)$ is the state occupied by the HS process at time
$s$, so that if $V(s)\in\mathfrak{V}_{0}$ then $V(s^{\prime})=V(s)$ for all
$s^{\prime}>s$. Similar to the LM process, let
$\mathfrak{I}_{\mathfrak{V}}=\\{(v_{1},v):v_{1}\in\mathfrak{V}_{1},v\in\mathfrak{V};v_{1}\neq
v\\}$, whose cardinality is
$|\mathfrak{I}_{\mathfrak{V}}|=|\mathfrak{V}_{1}||\mathfrak{V}|-|\mathfrak{V}_{1}|$.
We convert $V$ into a (column) $|\mathfrak{I}_{\mathfrak{V}}|$-dimensional
multivariate counting process
$\\{N^{V}\equiv(N^{V}(s;v_{1},v),(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}):s\geq
0\\}$, where $N^{V}(s;v_{1},v)$ is the number of observed transitions of $V$
from state $v_{1}$ into state $v$ over the period $[0,s]$, that is, for
$(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}$,

$N^{V}(s;v_{1},v)=\sum_{t\leq s}I\\{V(t-)=v_{1},V(t)=v\\}.$

For each $(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}$, the sample path $s\mapsto
N^{V}(s;v_{1},v)$ takes values in $\mathbb{Z}_{0,+}$, and is a nondecreasing
step-function, right-continuous with left-hand limits, and with
$N^{W}(0;v_{1},v)=0$. In addition,
$\sum_{(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}}dN^{V}(s;v_{1},v)\in\\{0,1\\}$
for every $s\geq 0$. Next, we combine the multivariate counting processes
$N^{R}$, $N^{W}$, and $N^{V}$ into one (column) multivariate counting process
$N=\\{N(s):s\geq 0\\}$ of dimension
$Q+|\mathfrak{I}_{\mathfrak{W}}|+|\mathfrak{I}_{\mathfrak{V}}|$, where, with T
denoting vector/matrix transpose,

$N(s)=\left[(N^{R}(s))^{\tt T},(N^{W}(s))^{\tt T},(N^{V}(s))^{\tt
T}\right]^{\tt T}.$

An important point needs to be stated regarding the observables in the study,
which will have an impact in the interpretation of the parameters of the joint
model. This pertains to the “competing risks” nature of all the possible
events at each time point $s$. The possible $Q$ recurrent event types, as well
as the potential transitions in the LM and HS processes, are all competing
with each other. Thus, suppose that at time $s_{0}$, the event that occurred
is a recurrent event of type $q_{0}$, that is, $dN^{R}(s_{0};q_{0})=1$. This
means that this event has occurred in the presence of the potential recurrent
events from the other $Q-1$ risks, and the potential transitions from either
the LM and HS processes. This will entail the use of crude hazards, instead of
net hazards, in the joint modeling, and this observation will play a critical
role in the dynamic joint model since each of the competing event occurrences
at a given time point $s$ from all the possible event sources (RCR, LM, and
HS) will be affected by the history of all these processes just before time
$s$. This is the aspect that exemplifies the synergistic association among the
three components.

Another observable process for our joint model is a vector of effective (or
virtual) age processes
$\mathcal{E}=\\{(\mathcal{E}_{1}(s),\ldots,\mathcal{E}_{Q}(s)):s\geq 0\\}$,
whose components are $\mathcal{F}$-predictable processes with sample paths
that are non-negative, left-continuous, piecewise nondecreasing and
differentiable. These effective age processes will carry the impact of
interventions performed after each recurrent event occurrence or a transition
in either the LM process or the HS process. For recent articles dealing with
virtual ages, see the philosophically-oriented article [12] and the very
recent review article [4].

Finally, we define the time-to-absorption of the unit to be
$\tau_{A}=\inf\\{s\geq 0:V(s)\in\mathfrak{V}_{0}\\}$ with the convention that
$\inf\emptyset=\infty$. Using this $\tau_{A}$ and $\tau$, we define the unit’s
at-risk process to be $Y=\\{Y(s):s\geq 0\\}$, with
$Y(s)=I\\{\min(\tau,\tau_{A})\geq s\\}.$ In addition, we define LM-specific
and HS-specific at-risk processes as follows: For $w\in\mathfrak{W}$, define
$Y^{W}(s;w)=I\\{W(s-)=w\\}$; and, for $v\in\mathfrak{V}_{1}$, define
$Y^{V}(s;v)=I\\{V(s-)=v\\}$. For a unit, we could then concisely summarize the
random observables in terms of stochastic processes as:

$D=(X,N,\mathcal{E},Y,Y^{W},Y^{V})\equiv\\{X,(N(s),\mathcal{E}(s),Y(s),Y^{W}(s),Y^{V}(s):s\geq
0)\\}.$ (1)

Note that the processes are undefined for
$s>\min(\tau_{A},\tau)\equiv\inf\\{s\geq 0:Y(s+)=0\\}$ since monitoring of the
unit had by then ceased.

###  3.2 Joint Model Specification for One Unit

The joint model specification will be through the specification of the
compensator process vector and the predictable quadratic variation (PQV)
process matrix of the multivariate counting process $N$. The predictable
process vector $A=\\{A(s):s\geq 0\\}$ is of dimension
$Q+|\mathfrak{I}_{\mathfrak{W}}|+|\mathfrak{I}_{\mathfrak{V}}|$ and is such
that the vector process $M=N-A=\\{M(s)=N(s)-A(s):s\geq 0\\}$ is a zero-mean
square-integrable martingale process with PQV matrix process $\langle
M,M\rangle$. The vectors $A$ and $M$ are actually partitioned into three
vector components reflecting the RCR, LM, and HS components, according to

$A=\left[(A^{R})^{\tt T},(A^{W})^{\tt T},(A^{V})^{\tt T}\right]^{\tt
T}\quad\mbox{and}\quad M=\left[(M^{R})^{\tt T},(M^{W})^{\tt T},(M^{V})^{\tt
T}\right],$

where, with $q\in\mathfrak{I}_{Q}$,
$(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}$,
$(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}$, and $s\geq 0$,

$\displaystyle A^{R}=\\{(A^{R}(s;q))\\}\quad\mbox{and}\quad
M^{R}=\\{(M^{R}(s;q))\\};$ $\displaystyle
A^{W}=\\{(A^{W}(s;(w_{1},w_{2}))\\}\quad\mbox{and}\quad
M^{W}=\\{(M^{W}(s;(w_{1},w_{2}))\\};$ $\displaystyle
A^{V}=\\{(A^{V}(s;(v_{1},v))\\}\quad\mbox{and}\quad
M^{V}=\\{(M^{V}(s;(v_{1},v))\\},$

with $A^{R}$ and $M^{R}$ of dimensions $Q$; $A^{W}$ and $M^{W}$ of dimensions
$|\mathfrak{I}_{\mathfrak{W}}|$; and $A^{V}$ and $M^{V}$ of dimensions
$|\mathfrak{I}_{\mathfrak{V}}|$. The matrix $\langle M,M\rangle$ could then be
partitioned similarly to reflect these block components.

We can now proceed with the specification of the compensator process vector
and the PQV process matrix. For conciseness, we introduce the generic mapping
$\iota$ defined as follows: For a set $A$ with $m$ elements,
$A=\\{a_{1},\ldots,a_{m}\\}$, let
$\iota_{A}(a)=(I(a_{2}=a),\ldots,I(a_{m}=a))$, a row vector of $m-1$ elements.
Here $\iota_{A}(a)$ is the indicator vector of $a$ excluding the first element
of $A$, so that $\iota_{A}(a_{1})=(\underline{0})$. Excluding $a_{1}$ in the
$\iota$ mapping is for purposes of model identifiability. One may think of the
mapping $\iota$ as a converter to dummy variables. We will also need the
following quantities or functions.

  * •

For each $q\in\mathfrak{I}_{Q}$ there is a baseline (crude) hazard rate
function $\lambda_{0q}(\cdot)$ with associated baseline (crude) cumulative
hazard function $\Lambda_{0q}(\cdot)$. We also denote by
$\bar{F}_{0q}(\cdot)=\prodi_{v=0}^{\cdot}[1-\Lambda_{0q}(dv)]$ the associated
baseline (crude) survivor function, where $\prodi$ is the product-integral
symbol.

  * •

For each $q\in\mathfrak{I}_{Q}$ there is a mapping
$\rho_{q}(\cdot;\cdot):\mathbb{Z}_{0,+}^{Q}\times\Re^{d_{q}}\rightarrow\Re_{0,+}$,
where the $d_{q}$’s are known positive integers. There is an associated vector
$\alpha_{q}\in\Re^{d_{q}}$.

  * •

There is a collection of non-negative real numbers
$\eta=\\{\eta(w_{1},w_{2}):(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}\\}$,
and we define for every $w_{1}\in\mathfrak{W}$,
$\eta(w_{1},w_{1})=-\sum_{w\in\mathfrak{W};w\neq w_{1}}\eta(w_{1},w).$

  * •

There is a collection of non-negative real numbers
$\xi=\\{\xi(v_{1},v):(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}\\}$, and we
define for every $v_{1}\in\mathfrak{V}$,
$\xi(v_{1},v_{1})=-\sum_{v\in\mathfrak{V};v\neq v_{1}}\xi(v_{1},v),$ and with
$\xi(v_{1},v_{2})=0$ for every $v_{1}\in\mathfrak{V}_{0}$ and
$v_{2}\in\mathfrak{V}$.

We then define the observables and associated finite-dimensional parameters,
respectively, for each of the three model components according to

$\displaystyle
B^{R}(s)=[X,\iota_{\mathfrak{V}}(V(s)),\iota_{\mathfrak{W}}(W(s))]\quad\mbox{and}\quad\theta^{R}=[(\beta^{R})^{\tt
T},(\gamma^{R})^{\tt T},(\kappa^{R})^{\tt T}]^{\tt T};$ $\displaystyle
B^{W}(s)=[X,\iota_{\mathfrak{V}}(V(s)),N^{R}(s)]\quad\mbox{and}\quad\theta^{W}=[(\beta^{W})^{\tt
T},(\gamma^{W})^{\tt T},(\nu^{W})^{\tt T}]^{\tt T};$ $\displaystyle
B^{V}(s)=[X,\iota_{\mathfrak{W}}(W(s)),N^{R}(s)]\quad\mbox{and}\quad\theta^{V}=[(\beta^{V})^{\tt
T},(\kappa^{V})^{\tt T},(\nu^{V})^{\tt T}]^{\tt T}.$

In addition to the $\lambda_{0q}$’s, $\alpha_{q}$’s, $\eta$, $\xi$, the
$\theta^{R}$, $\theta^{W}$, and $\theta^{V}$ will constitute all of the model
parameters. For the proposed model, the compensator process components are
given by, for $q\in\mathfrak{I}_{Q}$,
$(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}$, and
$(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}$:

$\displaystyle
A^{R}(s;q)=\int_{0}^{s}Y(t)\lambda_{0q}[\mathcal{E}_{q}(t)]\rho_{q}[N^{R}(t-);\alpha_{q}]\exp\\{B^{R}(t-)\theta^{R}\\}dt;$
$\displaystyle
A^{W}(s;w_{1},w_{2})=\int_{0}^{s}Y(t)Y^{W}(t;w_{1})\eta(w_{1},w_{2})\exp\\{B^{W}(t-)\theta^{W}\\}dt;$
$\displaystyle
A^{V}(s;v_{1},v)=\int_{0}^{s}Y(t)Y^{V}(t;v_{1})\xi(v_{1},v)\exp\\{B^{V}(t-)\theta^{V}\\}dt.$

In the left-hand side of the equations above, we have suppressed writing the
dependence on the model parameters. With $\mbox{Dg}(a)$ denoting the diagonal
matrix formed from vector $a$, the PQV process is specified to be

$\langle M,M\rangle(s)=\mbox{Dg}[(A^{R}(s))^{\tt T},(A^{W}(s))^{\tt
T},(A^{V}(s))^{\tt T}].$

Observe the dynamic nature of this model in that an event occurrence at an
infinitesimal time interval $[s,s+ds)$ depends on the history of the processes
before time $s$. According to the theory of counting processes, we have the
the following probabilistic interpretations (statements are almost surely):

$\displaystyle E\\{dN^{R}(s;q)|\mathcal{F}_{s-}\\}=dA^{R}(s;q);$
$\displaystyle
E\\{dN^{W}(s;w_{1},w_{2})|\mathcal{F}_{s-}\\}=dA^{W}(s;w_{1},w_{2});$
$\displaystyle E\\{dN^{V}(s;v_{1},v)|\mathcal{F}_{s-}\\}=dA^{V}(s;v_{1},v),$

and

$\displaystyle Var\\{dN^{R}(s;q)|\mathcal{F}_{s-}\\}=dA^{R}(s;q);$
$\displaystyle
Var\\{dN^{W}(s;w_{1},w_{2})|\mathcal{F}_{s-}\\}=dA^{W}(s;w_{1},w_{2});$
$\displaystyle Var\\{dN^{V}(s;v_{1},v)|\mathcal{F}_{s-}\\}=dA^{V}(s;v_{1},v),$

together with the conditional covariance, given $\mathcal{F}_{s-}$, between
any pair of elements of $dN(s)$ being equal to zero, e.g.,
$Cov\\{dN^{R}(s),dN^{W}(s;w_{1},w_{2})|\mathcal{F}_{s-}\\}=0$. However, note
that we are not assuming that the components of $N^{R}$, $N^{W}$, and $N^{V}$
are independent, nor even conditionally independent. A way to view this model
is that, given $\mathcal{F}_{s-}$, the history just before time $s$, on the
infinitesimal interval $[s,s+ds)$, $dN(s)=(dN^{R}(s)^{\tt T},dN^{W}(s)^{\tt
T},dN^{V}(s)^{\tt T})^{\tt T}$ has a multinomial distribution with parameters
$1$ and $dA(s)=(dA^{R}(s)^{\tt T},dA^{W}(s)^{\tt T},dA^{V}(s)^{\tt T})^{\tt
T}$. As such, for every $s\geq 0$, the following constraint holds:

$dN_{\bullet}(s)=dN_{\bullet}^{R}(s)+dN_{\bullet}^{W}(s)+dN_{\bullet}^{V}(s)\in\\{0,1\\},$

with the notation that a subscript of $\bullet$ means the sum over all the
appropriate index set, e.g.,
$dN_{\bullet}^{R}(s)=\sum_{q\in\mathfrak{I}_{Q}}dN^{R}(s;q)$ and
$dA_{\bullet}(s)=dA_{\bullet}^{R}(s)+dA_{\bullet}^{W}(s)+dA_{\bullet}^{V}(s)$.
The multinomial distribution above could actually be approximated by
independent Bernoulli distributions. To see this, if we have real numbers
$p_{k},k=1,\ldots,K,$ with $0<p_{k}\approx 0$ for each $k=1,\ldots,K,$ and
with $\sum_{k=1}^{K}p_{k}\approx 0$, then we have the approximation

$\left(1-\sum_{k=1}^{K}p_{k}\right)\approx\prod_{k=1}^{K}(1-p_{k}).$

Consequently, in the equation below, the multinomial probability on the left-
hand side is approximately the product of (independent) Bernoulli
probabilities in the right-hand side.

$\displaystyle\left\\{\prod_{q\in\mathfrak{I}_{Q}}[dA_{q}(s)]^{dN_{q}^{R}(s)}\right\\}\left\\{\prod_{(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}}[dA^{W}(s;w_{1},w_{2})]^{dN^{W}(s;w_{1},w_{2})}\right\\}\times$
$\displaystyle\left\\{\prod_{(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}}[dA^{V}(s;v_{1},v)]^{dN^{V}(s;v_{1},v)}\right\\}\left\\{[1-dA_{\bullet}(s)]^{1-dN_{\bullet}(s)}\right\\}$
$\displaystyle\approx$
$\displaystyle\left\\{\prod_{q\in\mathfrak{I}_{Q}}[dA_{q}(s)]^{dN_{q}^{R}(s)}[1-dA_{q}^{R}(s)]^{1-dN_{q}^{R}(s)}\right\\}$
$\displaystyle\left\\{\prod_{(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}}[dA^{W}(s;w_{1},w_{2})]^{dN^{W}(s;w_{1},w_{2})}[1-dA^{W}(s;w_{1},w_{2})]^{1-dN^{W}(s;w_{1},w_{2})}\right\\}$
$\displaystyle\left\\{\prod_{(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}}[dA^{V}(s;v_{1},v)]^{dN^{V}(s;v_{1},v)}[1-dA^{V}(s;v_{1},v)]^{1-dN^{V}(s;v_{1},v)}\right\\}.$

This approximate equivalence informs the manner in which we generate data from
the model later in the sections dealing with an illustrative data set (Section
6) and the simulation studies (Section 7) where we used this product of
Bernoulli approach.

Consider a unit that is still at risk just before time $s$ whose LM process is
at state $w_{1}$ and HS process at state $v_{1}\notin\mathfrak{V}_{0}$. Two
questions of interest are:

  * (a)

What is the distribution of the next event occurrence?

  * (b)

Given in addition that an event occurred at time $s+t$, what are the
conditional probabilities of each of the possible events?

Denote by $T$ the time to the next event occurrence starting from time $s$.
Then,

$\displaystyle\mathbf{P}\\{T>t|\mathcal{F}_{s-}\\}$ $\displaystyle=$
$\displaystyle\prodi_{u=s}^{s+t}\left[1-\left(\sum_{q=1}^{Q}\lambda_{0q}[\mathcal{E}_{q}(u)]\rho_{q}[N^{R}(u-);\alpha_{q}]\exp\\{B^{R}(u-)\theta^{R}\\}-\right.\right.$
$\displaystyle\left.\left.\eta(w_{1},w_{1})\exp\\{B^{W}(u-)\theta^{W}\\}-\xi(v_{1},v_{1})\exp\\{B^{V}(u-)\theta^{V}\\}\right)du\right]$
$\displaystyle=$
$\displaystyle\exp\left\\{-\int_{s}^{s+t}\left(\sum_{q=1}^{Q}\lambda_{0q}[\mathcal{E}_{q}(u)]\rho_{q}[N^{R}(u-);\alpha_{q}]\exp\\{B^{R}(u-)\theta^{R}\\}-\right.\right.$
$\displaystyle\left.\left.\eta(w_{1},w_{1})\exp\\{B^{W}(u-)\theta^{W}\\}-\xi(v_{1},v_{1})\exp\\{B^{V}(u-)\theta^{V}\\}\right)du\right\\}$
$\displaystyle=$
$\displaystyle\exp\left\\{-\exp\\{B^{R}(s-)\theta^{R}\\}\sum_{q=1}^{Q}\rho_{q}[N^{R}(s-);\alpha_{q}]\int_{s}^{s+t}\lambda_{0q}[\mathcal{E}_{q}(u)]du+\right.$
$\displaystyle\left.\eta(w_{1},w_{1})\exp\\{B^{W}(s-)\theta^{W}\\}t+\xi(v_{1},v_{1})\exp\\{B^{V}(s-)\theta^{V}\\}t\right\\},$

with the second equality obtained by invoking the product-integral identity

$\prodi_{s\in I}[1-dA(s)]^{1-dN(s)}=\exp\left\\{-\int_{s\in I}dA(s)\right\\}$

and since no events in $[s,s+t)$ means
$dN_{\bullet}^{R}(u)+dN_{\bullet}^{W}(u)+dN_{\bullet}^{V}(u)=0$ for
$u\in[s,s+t)$, and the last equality arising since, prior to the next event,
there will be no changes in the values of $N^{R}$, $B^{R}$, $B^{W}$, and
$B^{V}$ from their respective values just before time $s$. In particular, if
the $\lambda_{0q}$s are constants, corresponding to the hazard rates of an
exponential distribution, then the distribution of the time to the next event
occurrence is exponential. Given that the event occurred at time $s+t$, then
the conditional probability that it was an RCR-type $q$ event is
$e^{\\{B^{R}(s-)\theta^{R}\\}}\rho_{q}[N^{R}(s-);\alpha_{q}]\lambda_{0q}[\mathcal{E}_{q}(s+t)]/C(s,t),$
with

$\displaystyle
C(s,t)=e^{\\{B^{R}(s-)\theta^{R}\\}}\sum_{q^{\prime}=1}^{Q}\rho_{q^{\prime}}[N^{R}(s-);\alpha_{q^{\prime}}]\lambda_{0q^{\prime}}[\mathcal{E}_{q^{\prime}}(s+t)]-$
$\displaystyle\eta(w_{1},w_{1})e^{\\{B^{W}(s-)\theta^{W}\\}}-\xi(v_{1},v_{1})e^{\\{B^{V}(s-)\theta^{V}\\}}$

Similarly, the conditional probability that it was a transition to state
$w_{2}$ for the LM process is
$\eta(w_{1},w_{2})e^{\\{B^{W}(s-)\theta^{W}\\}}/C(s,t),$ and the conditional
probability that is was a transition to state $v$, possibly to an absorbing
state, for the HS process is
$\xi(v_{1},v)e^{\\{B^{V}(s-)\theta^{V}\\}}/C(s,t).$ These probabilities
demonstrate the competing risks nature of the different possible events. They
also provide a computational approach to iteratively generate data from the
joint model for use in simulation studies, with the basic idea being to first
generate a time to any type of event, then to mark the type of event or update
each of the counting processes by using the above conditional probabilities.

Denoting by $\Theta$ the set of all parameters of the model, the likelihood
function arising from observing $D$, with $p_{(W,V)}(\cdot,\cdot)$ the initial
joint probability mass function of $(W(0),V(0))$, is given by

$\displaystyle\mathcal{L}(\Theta|D)=p_{(W,V)}(W(0),V(0))\times$
$\displaystyle\prodi_{s=0}^{\infty}\left\\{\left[\prod_{q\in\mathfrak{I}_{Q}}[dA^{R}(s;q)]^{dN^{R}(s;q)}\right]\left[\prod_{(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}}[dA^{W}(s;w_{1},w_{2})]^{dN^{W}(s;w_{1},w_{2})}\right]\right.\times$
$\displaystyle\left.\left[\prod_{(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}}[dA^{V}(s;v_{1},v)]^{dN^{V}(s;v_{1},v)}\right]\left[1-dA_{\bullet}(s)\right]^{1-dN_{\bullet}(s)}\right\\}$

The likelihood in (3.2) could be rewritten in the form

$\displaystyle\mathcal{L}(\Theta|D)=p_{(W,V)}(W(0),V(0))\times$
$\displaystyle\prodi_{s=0}^{\infty}\left\\{\left[\prod_{q\in\mathfrak{I}_{Q}}[dA^{R}(s;q)]^{dN^{R}(s;q)}\right]\left[\prod_{(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}}[dA^{W}(s;w_{1},w_{2})]^{dN^{W}(s;w_{1},w_{2})}\right]\right.$
$\displaystyle\left.\left[\prod_{(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}}[dA^{V}(s;v_{1},v)]^{dN^{V}(s;v_{1},v)}\right]\right\\}\exp\\{-A_{\bullet}(\infty)\\}.$

Let us examine the meaning of $A_{\bullet}(\infty)\equiv
A_{\bullet}^{R}(\infty)+A_{\bullet}^{W}(\infty)+A_{\bullet}^{V}(\infty)$.
Simplifying, we see that this equals
$A_{\bullet}(\infty)=\int_{0}^{\infty}Y(s)T(s)ds=\int_{0}^{\tau\wedge\tau_{A}}T(s)ds$,
where

$\displaystyle T(s)$ $\displaystyle=$
$\displaystyle\sum_{q\in\mathfrak{I}_{Q}}\lambda_{0q}[\mathcal{E}_{q}(s)]\rho_{q}[N^{R}(s-);\alpha_{q}]\exp\\{B^{R}(s-)\theta^{R}\\}-$
$\displaystyle\sum_{w_{1}\in\mathfrak{W}}Y^{W}(s;w_{1})\eta(w_{1},w_{1})\exp\\{B^{W}(s-)\theta^{W}\\}-$
$\displaystyle\sum_{v_{1}\in\mathfrak{V}_{1}}Y^{V}(s;v_{1})\xi(v_{1},v_{1})\exp\\{B^{V}(s-)\theta^{V}\\}.$

Recall that $\eta(w_{1},w_{1})$ and $\xi(v_{1},v_{1})$ are non-positive real
numbers. Thus, $T(s)ds$ could be interpreted as the total risk of an event,
either a recurrent event in the RCR component or a transition in the LM or HS
components, occurring from all possible sources (RCR, LM, HS) that the unit is
exposed to at the infinitesimal time interval $[s,s+ds)$, given the history
$\mathcal{F}_{s-}$ just before time $s$.

We provide further explanations of the elements of the joint model. First,
there is a tacit assumption that no more than one event of any type could
occur at any given time $s$. Second, the event rate at any time point $s$ for
any type of event is in the presence of the other possible risk events. Thus,
consider a specific $q_{0}\in\mathfrak{I}_{Q}$ and assume that the unit is
still at-risk at time $s_{0}$. Then,

$\displaystyle\mathbf{P}\\{dN_{q_{0}}^{R}(s_{0})=1|\mathcal{F}_{s_{0}-}\\}$
$\displaystyle\approx$
$\displaystyle\lambda_{0q_{0}}[\mathcal{E}_{q_{0}}(s_{0})]\rho_{q_{0}}[N^{R}(s_{0}-);\alpha_{q_{0}}]\exp\\{B^{R}(s_{0}-)\theta^{R}\\}d{s_{0}}$
(4)

is the conditional probability, given $\mathcal{F}_{s_{0}-}$, that an event
occurred at $[s_{0},s_{0}+ds_{0})$ and is of RCR type $q_{0}$ and all other
event types did not occur, which is the essence of what is called a crude
hazard rate, instead of a net hazard rate. Third, the effective (or virtual)
age functions $\mathcal{E}_{q}(\cdot)$s, which are assumed to be dynamically
determined and not dependent on any unknown parameters, encodes the impact of
performed interventions that are performed after each event occurrence.
Several possible choices of these functions are:

  * •

$\mathcal{E}_{q}(s)=s$ for all $s\geq 0$ and $q\in\mathfrak{I}_{Q}$. This is
usually referred to as a minimal repair intervention, corresponding to the
situation where an intervention simply puts back the system at the age just
before the event occurrence.

  * •

$\mathcal{E}_{q}(s)=s-S_{N_{\bullet}(s-)}$ where $0=S_{0}<S_{1}<S_{2}<\ldots$
are the successive event times. This corresponds to a perfect intervention,
which has the effect of re-setting the time to zero for each of the RCRs after
each event occurrence. In a reliability setting, this means that all $Q$
components (unless having exponential lifetimes) are replaced by corresponding
identical, but new, components at each event occurrence.

  * •

$\mathcal{E}_{q}(s)=s-S_{N_{\bullet}^{R}(s-)}^{R}$ where
$0=S_{0}^{R}<S_{1}^{R}<S_{2}^{R}<\ldots$ are the successive event times of the
occurrences of the RCR events.

  * •

$\mathcal{E}_{q}(s)=s-S_{qN^{R}(s-;q)}^{R}$ where
$0=S_{q0}^{R}<S_{q1}^{R}<S_{q2}^{R}<\ldots$ are the successive event times of
the occurrences of RCR events of type $q$.

  * •

Other general forms are possible, such as those in [9] and [14], the latter
employing ideas of Kijima. See also the discussion on the ‘reality’ of virtual
or effective ages in the paper by [12], as well as the recent review paper by
[4].

Fourth, the impact of accumulating RCR event occurrences, which could be
adverse, but could also be beneficial as in software engineering applications,
is incorporated in the model through the $\rho_{q}(\cdot;\cdot)$ functions.
One possible choice is an exponential function, such as
$\rho_{q}(N^{R}(s-);\alpha_{q})=\exp\\{N^{R}(s-)^{\tt T}\alpha_{q}\\}$, but
other choices could be made as well. Finally, the modulating exponential link
function in the model is for the impact of the covariates as well as the
values of the RCR, LM, and HS processes just before the time of interest, with
the vector of coefficients quantifying the effects of the covariates. The use
of $N(s-)$ in the model could be viewed as using them as internal covariates,
or that the dynamic model is a self-exciting model.

Similar interpretations hold for the parameters
$\\{\eta(w_{1},w_{2}):(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}\\}$ and
$\\{\xi(v_{1},v):(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}\\}$. Thus, if just
before time $s_{0}$, $W(s_{0}-)=w_{1}$ and $V(s_{0}-)=v_{1}$, indicated by
$\mathcal{F}_{s_{0}-}(w_{1},v_{1})$, then

$\displaystyle\mathbf{P}\\{W(s_{0}+ds_{0})=w_{2}|\mathcal{F}_{s_{0}}(w_{1},v_{1})\\}\approx$
(5) $\displaystyle\eta(w_{1},w_{2})$
$\displaystyle\exp\\{X\beta^{W}+\gamma^{W}_{j(v_{1})}+N^{R}(s_{0}-)\nu^{W}\\}ds_{0};$
$\displaystyle\mathbf{P}\\{V(s_{0}+ds_{0})=v|\mathcal{F}_{s_{0}}(w_{1},v_{1})\\}\approx$
(6) $\displaystyle\xi(v_{1},v)$
$\displaystyle\exp\\{X\beta^{V}+\kappa^{V}_{j(w_{1})}+N^{R}(s_{0}-)\nu^{V}\\}ds_{0},$

where $j(v_{1})$ is the index associated with $v_{1}$ in $\mathfrak{V}_{1}$
and $j(w_{1})$ is the index associated with $w_{1}$ in $\mathfrak{W}$.

###  3.3 Special Case: Independent Poisson Processes and CTMCs for One Unit

There is a special case arising from this general joint model obtained when we
set $\lambda_{0q}(s)=\lambda_{0q}$, $q\in\mathfrak{I}_{Q}$; $\rho_{q}=1$;
$\theta^{R}=0$; $\theta^{W}=0$; and $\theta^{V}=0$. In this situation, we have

$\displaystyle dA^{R}(s;q)$ $\displaystyle=$
$\displaystyle\lambda_{0q}ds,q\in\mathfrak{I}_{Q};$ $\displaystyle
dA^{W}(s;w_{1},w_{2})$ $\displaystyle=$
$\displaystyle\eta(w_{1},w_{2})Y(s)Y^{W}(s;w_{1})ds,(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}};$
$\displaystyle dA^{V}(s;v_{1},v)$ $\displaystyle=$
$\displaystyle\xi(v_{1},v)Y(s)Y^{W}(s,v_{1})ds,(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}.$

It is easy to see that this particular model coincides with the model where we
have the following situations:

  * (i)

$N^{R}(\cdot;q),q\in\mathfrak{I}_{Q},$ are independent homogeneous Poisson
processes with respective rates $\lambda_{0q},q\in\mathfrak{I}_{Q}$;

  * (ii)

$W(\cdot)$ is a continuous-time Markov chain (CTMC) with infinitesimal
generator matrix (IGM) consisting of $\\{\eta(w_{1},w_{2})\\}$;

  * (iii)

$V(\cdot)$ is a CTMC with IGM consisting of $\\{\xi(v_{1},v_{2})\\}$;

  * (iv)

$N^{R}$, $W$, and $V$ are independent; and

  * (v)

Processes are observed over $[0,\min(\tau,\tau_{A})]$, where $\tau$ is the end
of monitoring period, while $\tau_{A}$ is the absorption time of $V$ into
$\mathfrak{V}_{0}$.

In this special situation, the $\lambda_{0q}$’s are both crude and net hazard
rates. Also, due to the memoryless property of the exponential distribution,
interventions performed after each event occurrence will have no impact in the
succeeding event occurrences. This specific joint model further allows us to
provide an operational interpretation of the model parameters. Thus, suppose
that at time $s$, the LM process is at state $w_{1}$ and the HS process is at
state $v_{1}\notin\mathfrak{V}_{0}$. Then, the distribution of the time to the
next event occurrence of any type (the holding or sojourn time at the current
state configuration) has an exponential distribution with parameter
$C=\lambda_{0\bullet}-\eta(w_{1},w_{1})-\xi(v_{1},v_{1})$. When an event
occurs, then the (conditional) probability that it is (i) an RCR event of type
$q$ is $\lambda_{0q}/C$; (ii) a transition to LM state $w_{2}\neq w_{1}$ is
$\eta(w_{1},w_{2})/C$; or (iii) a transition to an HS state $v\neq v_{1}$ is
$\xi(v_{1},v)/C$. This is the essence of the competing risks nature of all the
possible event types: an RCR event, an LM transition, and an HS transition. As
such, the more general model could be viewed as an extension of this basic
model with independent Poisson processes for the RCR component and CTMCs for
the LM and HS components. For this special case, the likelihood function in
(3.2) simplifies to the expression

$\displaystyle\mathcal{L}(\Theta|D)=p_{(W,V)}(W(0),V(0))\left[\prod_{q\in\mathfrak{I}_{Q}}\lambda_{0q}^{N^{R}(\tau\wedge\tau_{A};q)}\right]\times$
$\displaystyle\left[\prod_{(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}}\eta(w_{1},w_{2})^{N^{W}(\tau\wedge\tau_{A};w_{1},w_{2})}\right]\left[\prod_{(v_{1},v)\in\mathfrak{I}_{\mathfrak{W}}}\xi(v_{1},v)^{N^{V}(\tau\wedge\tau_{A};v_{1},v)}\right]\times$
$\displaystyle\exp\left\\{-\int_{0}^{\tau\wedge\tau_{A}}T(s)ds\right\\}.$

where
$T(s)=\lambda_{0\bullet}-\sum_{w_{1}\in\mathfrak{W}}\eta(w_{1},w_{1})Y^{W}(s;w_{1})-\sum_{v_{1}\in\mathfrak{V}_{1}}\xi(v_{1},v_{1})Y^{V}(s;v_{1}).$
Note that $T(s)=T(s;\lambda_{0},\eta,\xi)$, that is, it is a quantity instead
of a statistic.

##  4 Estimation of Model Parameters

###  4.1 Parametric Estimation

Having introduced the joint model, we now address in this section the problem
of making inferences about the model parameters. We assume that we are able to
observe $n$ independent units, with the $i$th unit having data
$D_{i}=(X_{i},N_{i},\mathcal{E}_{i},Y_{i},Y_{i}^{W},Y_{i}^{V})$ as in (1). The
full sample data will then be represented by

$\mathbf{D}=(D_{1},D_{2},\ldots,D_{n}),$ (8)

while the model parameters will be represented by, with the convention that
$q\in\mathfrak{I}_{Q}$, $(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}$, and
$(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}$,

$\displaystyle\Theta$ $\displaystyle\equiv$
$\displaystyle\left[\\{\lambda_{0q}(\cdot),\alpha_{q}\\},\\{\eta(w_{1},w_{2})\\},\\{\xi(v_{1},v)\\},\theta^{R},\theta^{W},\theta^{V}\right].$

The $\lambda_{0q}$s could be parametrically-specified, hence will have finite-
dimensional parameters, so $\Theta$ will also then be finite-dimensional.
Except for the special case mentioned above, our main focus will be the case
where the $\lambda_{0q}$s are nonparametric. The distributions, $G_{i}$s, of
the end of monitoring periods, $\tau_{i}$s, also have model parameters, but
they are not of main interest. To visualize the type of sample data set that
accrues, Figure 2 provides a picture of a simulated sample data with $n=50$
units.

|
---|---
|
Figure 2: Simulated sample data with $n=50$ units. Panel 1: Covariate values
of the $n=50$ units. Values of $X_{2}$ are indicated by rectangles if
$X_{1}=0$, while values of $X_{1}$ are indicated by circles if $X_{1}=1$. Here
$X_{1}\sim\mbox{BER}(0.5)$, $X_{2}\sim N(0,1)$. $X_{1}$ and $X_{2}$ are
generated independently of each other. Panel 2: Recurrent competing risks
occurrences with three types of competing risks. Each unit is either censored
(“+”) or reaches the absorbing status (“$\times$”). Panel 3: Marker processes.
Panel 4: Health status processes, with state “1” absorbing.

The full likelihood function, given $\mathbf{D}$, is
$\mathcal{L}(\Theta|\mathbf{D})=\prod_{i=1}^{n}\mathcal{L}(\Theta|\mathbf{D}_{i})$,
where the $\mathcal{L}(\Theta|D_{i})$ is of the form in (3.2). If the
$\lambda_{0q}(\cdot)$s are parametrically-specified, then estimators of the
finite-dimensional model parameters could be obtained as the maximizers of
this full likelihood function, and their finite and asymptotic properties will
follow from the general results for maximum likelihood estimators based on
counting processes; see, for instance, [6] and [1].

We illustrate this situation for the special case of the model given in
subsection 3.3, so that the parameter is simply
$\Theta=[\\{\lambda_{0q}\\},\\{\eta(w_{1},w_{2})\\},\\{\xi(v_{1},v)\\}]$. In
this situation, from (3.3), the full likelihood reduces to, with
$\tau_{i}^{*}=\tau_{i}\wedge\tau_{iA}$,

$\displaystyle\mathcal{L}(\Theta|\mathbf{D})=\prod_{i=1}^{n}p_{(W,V)}(W_{i}(0),V_{i}(0))\times$
$\displaystyle\left[\prod_{q\in\mathfrak{I}_{Q}}\lambda_{0q}^{\sum_{i=1}^{n}N_{i}^{R}(\tau_{i}^{*};q)}\right]\left[\prod_{(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}}\eta(w_{1},w_{2})^{\sum_{i=1}^{n}N_{i}^{W}(\tau_{i}^{*};w_{1},w_{2})}\right]\times$
$\displaystyle\left[\prod_{(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}}\xi(v_{1},v)^{\sum_{i=1}^{n}N_{i}^{V}(\tau_{i}^{*};v_{1},v)}\right]\exp\left\\{-\sum_{i=1}^{n}\int_{0}^{\tau_{i}^{*}}T_{i}(s;\lambda_{0},\eta,\xi)ds\right\\},$

where

$T_{i}(s;\lambda_{0},\eta,\xi)=\lambda_{0\bullet}-\sum_{w_{1}\in\mathfrak{W}}\eta(w_{1},w_{1})Y_{i}^{W}(s;w_{1})-\sum_{v_{1}\in\mathfrak{V}}\xi(v_{1},v_{1})Y_{i}^{V}(s;v_{1}).$

The score function
$U(\Theta|\mathbf{D})=\nabla_{\Theta}\log\mathcal{L}(\Theta|\mathbf{D})$ has
elements

$\displaystyle U^{R}(\Theta;q)$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}^{R}(\tau_{i}^{*};q)}{\lambda_{0q}}-\sum_{i=1}^{n}\tau_{i}^{*},\quad
q\in\mathfrak{I}_{Q};$ $\displaystyle U^{W}(\Theta;w_{1},w_{2})$
$\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}^{W}(\tau_{i}^{*};w_{1},w_{2})}{\eta(w_{1},w_{2})}-\sum_{i=1}^{n}\int_{0}^{\tau_{i}^{*}}Y_{i}^{W}(s;w_{1})ds,\quad(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}};$
$\displaystyle U^{V}(\Theta;v_{1},v)$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}^{V}(\tau_{i}^{*};v_{1},v)}{\xi(v_{1},v)}-\sum_{i=1}^{n}\int_{0}^{\tau_{i}^{*}}Y_{i}^{V}(s;v_{1})ds,\quad(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}.$

Equating these equations to zeros yield the ML estimators of the parameters,
which are given below and possess the interpretation of being the observed
“occurrence-exposure" rates.

$\displaystyle\hat{\lambda}_{0q}$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}^{R}(\tau_{i}^{*};q)}{\sum_{i=1}^{n}\tau_{i}^{*}}=\frac{\sum_{i=1}^{n}\int_{0}^{\infty}dN_{i}^{R}(s;q)}{\sum_{i=1}^{n}\int_{0}^{\infty}Y_{i}(s)ds},q\in\mathfrak{I}_{Q};$
$\displaystyle\hat{\eta}(w_{1},w_{2})$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}^{W}(\tau_{i}^{*};w_{1},w_{2})}{\sum_{i=1}^{n}\int_{0}^{\tau_{i}^{*}}Y_{i}^{W}(s;w_{1})ds}=\frac{\int_{0}^{\infty}dN_{i}^{W}(s;w_{1},w_{2})}{\sum_{i=1}^{n}\int_{0}^{\infty}Y_{i}(s)Y_{i}^{W}(s;w_{1})ds},(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}};$
$\displaystyle\hat{\xi}(v_{1},v)$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}^{V}(\tau_{i}^{*};v_{1},v)}{\sum_{i=1}^{n}\int_{0}^{\tau_{i}^{*}}Y_{i}^{V}(s;v_{1})ds}=\frac{\int_{0}^{\infty}dN_{i}^{V}(s;v_{1},v)}{\sum_{i=1}^{n}\int_{0}^{\infty}Y_{i}(s)Y_{i}^{V}(s;v_{1})ds},(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}.$

In these estimators, the numerators are total event counts, e.g.,
$\sum_{i=1}^{n}N_{i}^{R}(\tau_{i}^{*};q)$ is the total number of observed RCR
type $q$ events; $\sum_{i=1}^{n}N_{i}^{W}(\tau_{i}^{*};w_{1},w_{2})$ is the
total number of observed transitions in the LM process from state $w_{1}$ into
$w_{2}$; and $\sum_{i=1}^{n}N_{i}^{V}(\tau_{i}^{*};v_{1},v)$ is the total
number of observed transitions in the HS process from state $v_{1}$ into $v$.
On the other hand, the denominators are total observed exposure times, e.g.,
$\sum_{i=1}^{n}\int_{0}^{\infty}Y_{i}(s)ds$ is the total time at-risk for all
the units; $\sum_{i=1}^{n}\int_{0}^{\infty}Y_{i}(s)Y_{i}^{W}(s;w_{1})dw$ is
the total observed time of all units that they were at-risk for a transition
in the LM process from state $w_{1}$; and
$\sum_{i=1}^{n}\int_{0}^{\infty}Y_{i}(s)Y_{i}^{V}(s;v_{1})ds$ is the total
observed time of all units that they were at-risk for a transition in the HS
process from state $v_{1}$. An important and crucial point to emphasize here
is that in these exposure times, they all take into account the time after the
last observed events in each component process until the end of monitoring,
whether it is a censoring (reaching $\tau_{i}$) or an absorption (reaching
$\tau_{iA}$). If one ignores these right-censored times, then the estimators
could be severely biased. This is a critical aspect we mentioned in the
introductory section and re-iterate at this point that this should not be
glossed over when dealing with recurrent event models.

The elements of the observed information matrix,
$I(\Theta;\mathbf{D})=-\nabla_{\Theta^{\tt T}}U(\Theta|\mathbf{D})$, which is
a diagonal matrix, have diagonal elements given by:

$\displaystyle I^{R}(\Theta;q)$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}^{R}(\tau_{i}^{*};q)}{\lambda_{0q}^{2}},q\in\mathfrak{I}_{Q};$
$\displaystyle I^{W}(\Theta;w_{1},w_{2})$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}^{W}(\tau_{i}^{*};w_{1},w_{2})}{\eta(w_{1},w_{2})^{2}},(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}};$
$\displaystyle I^{V}(\Theta;v_{1},v)$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}N_{i}^{V}(\tau_{i}^{*};v_{1},v)}{\xi(v_{1},v)^{2}},(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}.$

Abbreviating the estimators into
$\hat{\Theta}=(\hat{\lambda}_{0},\hat{\eta},\hat{\xi})$, we obtain the
asymptotic result, that as $n\rightarrow\infty$,

$\hat{\Theta}\sim\mbox{AsyMVN}(\Theta,I(\hat{\Theta};\mathbf{D})^{-1}),$ (9)

with AsyMVN meaning asymptotically multivariate normal. Thus, this result
seems to indicate that the RCR, LM, and HS components or the estimators of
their respective parameters do not have anything to do with each other, which
appears intuitive since the RCR, LM, and HS processes were assumed to be
independent processes to begin with. But, let us examine this issue further.
The result in (9) is an approximation to the theoretical result that

$\hat{\Theta}\sim\mbox{AsyMVN}\left(\Theta,\frac{1}{n}\mathfrak{I}(\Theta)^{-1}\right),$
(10)

where $\frac{1}{n}I(\hat{\Theta};\mathbf{D})\stackrel{{\scriptstyle
pr}}{{\rightarrow}}\mathfrak{I}(\Theta)$. Evidently, $\mathfrak{I}$ is a
diagonal matrix, so let us examine its diagonal elements. Let
$q\in\mathfrak{I}_{Q}$. Then we have, with ‘pr-lim’ denoting in-probability
limit,

$\displaystyle\lambda_{0q}^{2}\mathfrak{I}^{R}(\Theta;q)$ $\displaystyle=$
$\displaystyle\mbox{pr-
lim}_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}dN_{i}^{R}(s;q)$
$\displaystyle=$ $\displaystyle\mbox{pr-
lim}_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}dM_{i}^{R}(s;q)+\lambda_{0q}\left[\mbox{pr-
lim}_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}Y_{i}(s)ds\right].$

The first term on the right-hand side (RHS) converges in probability to zero
by the weak law of large numbers and the zero-mean martingale property. The
second term in the RHS converges in probability to its expectation, hence

$\mathfrak{I}^{R}(\Theta;q)=\frac{1}{\lambda_{0q}}\left[\int_{0}^{\infty}\left\\{\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i}^{n}E[Y_{i}(s)]\right\\}ds\right].$

But, now,

$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i}^{n}E[Y_{i}(s)]$
$\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i}^{n}\mathbf{P}\\{\tau_{i}\geq
s\\}\mathbf{P}\\{\tau_{i}^{A}\geq s\\}$ $\displaystyle=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i}^{n}\bar{G}_{i}(s-)\mathbf{P}\\{V_{i}(u)\notin\mathfrak{V}_{0},u\leq
s\\}.$

with $\bar{G}_{i}=1-G_{i}$. The last probability term above will depend on the
generators $\\{\xi(v_{1},v):(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}\\}$ of the
CTMC $\\{V(s):s\geq 0\\}$, so that the theoretical Fisher information or the
asymptotic variance associated with the estimator $\hat{\lambda}_{0q}$ depends
after all on the HS process, as well as on the $G_{i}$s, contrary to the
seemingly intuitive expectation that it should not depend on the LM and HS
processes. This result is a subtle one which arise because of the structure of
the observation processes. If $G_{i}=G,i=1,\ldots,n$, then we have that

$\mathfrak{I}^{R}(\Theta;q)=\frac{1}{\lambda_{0q}}\int_{0}^{\infty}\bar{G}(s-)\mathbf{P}\\{V(u)\notin\mathfrak{V}_{0},u\leq
s\\}ds,$

since $\mathbf{P}\\{V_{i}(u)\notin\mathfrak{V}_{0},u\leq
s\\}=\mathbf{P}\\{V(u)\notin\mathfrak{V}_{0},u\leq s\\},i=1,\ldots,n$. If we
denote by $\Gamma$ the generator matrix of $\\{V(s):s\geq 0\\}$ and let
$\Gamma_{1}$ be the sub-matrix associated with the $\mathfrak{V}_{1}$ states,
then if $p_{0}^{V}\equiv(p_{V}(v_{1}),v_{1}\in\mathfrak{V}_{1})^{\tt T}$ is
the initial probability mass function of $V(0)$, we have

$\mathbf{P}\\{V(u)\notin\mathfrak{V}_{0},u\leq
s\\}=\mathbf{P}\\{V(s)\in\mathfrak{V}_{1}\\}=(p_{0}^{V})^{\tt
T}\left[e^{s\Gamma_{11}}\right]1_{|\mathfrak{V}_{1}|}$

where the matrix exponential is
$e^{s\Gamma_{11}}\equiv\sum_{k=0}^{\infty}\frac{s^{k}\Gamma_{11}^{k}}{k!}$ and
$1_{K}$ is a column vector of $1$s of dimension $K$. Thus, we obtain

$\mathfrak{I}^{R}(\Theta;q)=\frac{1}{\lambda_{0q}}(p_{0}^{V})^{\tt
T}\sum_{k=0}^{\infty}\left[\int_{0}^{\infty}\bar{G}(s-)\frac{s^{k}}{k!}ds\right]\Gamma_{11}^{k}1_{|\mathfrak{V}_{1}|}.$

For example, if $\bar{G}(s)=\exp(-\nu s)$, that is, $\tau_{i}$’s are
exponentially-distributed with mean $1/\nu$, then the above expression
simplifies to

$\mathfrak{I}^{R}(\Theta;q)=\frac{1}{\lambda_{0q}}\frac{1}{\nu}(p_{0}^{V})^{\tt
T}\left[\sum_{k=0}^{\infty}\left(\frac{\Gamma_{11}}{\nu}\right)^{k}\right]1_{|\mathfrak{V}_{1}|}.$

For computational purposes, one may use an eigenvalue decomposition of
$\Gamma_{11}$: $\Gamma_{11}=U\mbox{Dg}(d)U^{-1}$ where $d$ consists of the
eigenvalues of $\Gamma_{11}$ and $U$ is the matrix of eigenvectors associated
with the eigenvalues $d$. The main point of this example though is the
demonstration that estimators of the parameters associated with the RCR, LM,
or HS process will depend on features of the other processes, even when one
starts with independent processes.

We remark that the estimators $\hat{\lambda}_{0q}$s,
$\hat{\eta}(w_{1},w_{2})$s, and $\hat{\xi}(v_{1},v)$s could also be derived as
method-of-moments estimators using the martingale structure. The inverse of
the observed Fisher information matrix coincides with an estimator using the
optional variation (OV) matrix process, while the inverse of the Fisher
information matrix coincides with the limit-in-probability of the predictable
quadratic variation matrix. To demonstrate for $\lambda_{0q}$, we have that

$\left\\{\sum_{i=1}^{n}M_{i}^{R}(s;q)=\sum_{i=1}^{n}\left[N_{i}^{R}(s;q)-\int_{0}^{s}Y_{i}(t)\lambda_{0q}dt\right]:s\geq
0\right\\}$

is a zero-mean square-integrable martingale. Letting $s\rightarrow\infty$,
setting $\sum_{i=1}^{n}M_{i}^{R}(\infty;q)=0$, and solving for $\lambda_{0q}$
yields $\hat{\lambda}_{0q}$. Next, we have

$\displaystyle\hat{\lambda}_{0q}$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}\int_{0}^{\infty}dN_{i}^{R}(s;q)}{\sum_{i=1}^{n}\int_{0}^{\infty}Y_{i}(s)ds}=\lambda_{0q}+\frac{\sum_{i=1}^{n}\int_{0}^{\infty}dM_{i}^{R}(s;q)}{\sum_{i=1}^{n}\int_{0}^{\infty}Y_{i}(s)ds}$

so that

$\displaystyle\sqrt{n}[\hat{\lambda}_{0q}-\lambda_{0q}]=\left(\int_{0}^{\infty}\frac{1}{n}\sum_{i=1}^{n}Y_{i}(s)ds\right)^{-1}\frac{1}{\sqrt{n}}\sum_{I=1}^{n}\int_{0}^{\infty}dM_{i}^{R}(s;q).$

We have already seen where
$\int_{0}^{\infty}\frac{1}{n}\sum_{i=1}^{n}Y_{i}(s)ds$ converges in
probability, whereas by the Martingale Central Limit Theorem, we have that

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\int_{0}^{\infty}dM_{i}^{R}(s;q)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\sigma^{2}_{R}(q))$

with

$\displaystyle\sigma_{R}^{2}(q)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\left\\{\mbox{pr-
lim}_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}d\langle
M_{i}^{R}(\cdot;q),M_{i}^{R}(\cdot;q)\rangle(s)\right\\}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\left\\{\mbox{pr-
lim}_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}Y_{i}(s)\lambda_{0q}\right\\}ds.$

Therefore, we have

$\displaystyle\sqrt{n}[\hat{\lambda}_{0q}-\lambda_{0q}]\stackrel{{\scriptstyle
d}}{{\rightarrow}}$ $\displaystyle
N\left(0,\frac{\lambda_{0q}}{\int_{0}^{\infty}\left[\mbox{pr-
lim}_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}Y_{i}(s)\right]ds}=\frac{\lambda_{0q}}{\int_{0}^{\infty}\left[\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}E[Y_{i}(s)]\right]ds}\right),$

which is the same result stated above using ML theory. Analogous asymptotic
derivations can be done for $\hat{\eta}(w_{1},w_{2})$ and
$\hat{\xi}(v_{1},v)$, though the resulting limiting variances will involve
expected occupation times for their respective states of the $W_{i}$-processes
coming from the $Y_{i}^{W}(\cdot;w_{1})$ terms and the $V_{i}$-processes from
the $Y_{i}^{V}(\cdot;v_{1})$ terms. Note that, asymptotically, these
estimators are independent, but their limiting variances depend on the
parameters from the other processes.

###  4.2 Semi-Parametric Estimation

In this section we consider the estimation of model parameters when the hazard
rate functions $\lambda_{0q}(\cdot)$s are specified nonparametrically. We
shall denote by
$\Lambda_{0q}(t)=\int_{0}^{t}\lambda_{0q}(u)du,q\in\mathfrak{I}_{Q}$, the
associated cumulative hazard functions. To simplify notation, we let

$\displaystyle\psi_{R}(s;\theta^{R})=\exp\\{B^{R}(s)\theta^{R}\\};\psi^{W}(s;\theta^{W})=\exp\\{B^{W}(s)\theta^{W}\\};\psi^{V}(s;\theta^{V})=\exp\\{B^{V}(s)\theta^{V}\\}.$

Using these functions, for $q\in\mathfrak{I}_{Q}$,
$(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}$, and
$(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}$, we then have

$\displaystyle dA^{R}(s;q)$ $\displaystyle=$ $\displaystyle
Y(s)\lambda_{0q}[\mathcal{E}_{q}(s)]\rho_{q}(N^{R}(s-);\alpha_{q})\psi^{R}(s-;\theta^{R})ds;$
$\displaystyle dA^{W}(s;w_{1},w_{2})$ $\displaystyle=$ $\displaystyle
Y(s)Y^{W}(s;w_{1})\eta(w_{1},w_{2})\psi^{W}(s-;\theta^{W})ds;$ $\displaystyle
dA^{V}(s;v_{1},v)$ $\displaystyle=$ $\displaystyle
Y(s)Y^{V}(s;v_{1})\xi(v_{1},v)\psi^{V}(s-;\theta^{V})ds.$

We also abbreviate the vector of model parameters into
$\Theta\equiv(\Lambda_{0},\alpha,\eta,\xi,\theta^{R},\theta^{W},\theta^{V})$.
Our goal is to obtain estimators for these parameters based on the sample data
$\mathbf{D}=(D_{1},D_{2},\ldots,D_{n})$. In a nutshell, the basic approach to
obtaining our estimators is to first assume that
$(\alpha,\theta^{R},\theta^{W},\theta^{V})$ are known, then obtain
‘estimators’ of $(\Lambda_{0},\eta,\xi)$. Having obtained these ‘estimators’,
in quotes since they are not yet estimators when
$(\alpha,\theta^{R},\theta^{W},\theta^{V})$ are unknown, we plug them into the
likelihood function to obtain a profile likelihood function. From the
resulting profile likelihood function, which depends on
$(\alpha,\theta^{R},\theta^{W},\theta^{V})$, we obtain its maximizers with
respect to these finite-dimensional parameters to obtain their estimators.
These estimators are then plugged into the ‘estimators’ of
$(\Lambda_{0},\eta,\xi)$ to obtain their estimators.

The full likelihood function based on the sample data
$\mathbf{D}=(D_{1},\ldots,D_{n})$ could be written as a product of three
“major” likelihood functions corresponding to the three model components:

$\displaystyle\mathcal{L}(\Theta|\mathbf{D})=\left[\prod_{q\in\mathfrak{I}_{Q}}\mathcal{L}^{R}(\Lambda_{0q},\alpha,\theta^{R};q|\mathbf{D})\right]\times$
$\displaystyle\left[\prod_{(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}}\mathcal{L}^{W}(\eta,\theta^{V};w_{1},w_{2}|\mathbf{D})\right]\times\left[\prod_{(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}}}\mathcal{L}^{V}(\xi,\theta^{W};v_{1},v|\mathbf{D})\right],$

where, suppressing writing of the parameters in the functions,

$\displaystyle\mathcal{L}^{R}(\Lambda_{0q},\alpha,\theta^{R};q|\mathbf{D})$
$\displaystyle=$
$\displaystyle\left\\{\prod_{i=1}^{n}\prodi_{s=0}^{\infty}\left[dA_{i}^{R}(s;q)\right]^{dN_{i}^{R}(s;q)}\right\\}\times$
$\displaystyle\exp\left\\{-\sum_{i=1}^{n}\int_{0}^{\infty}dA_{i}^{R}(s;q)\right\\};$
$\displaystyle\mathcal{L}^{W}(\eta,\theta^{V};w_{1},w_{2}|\mathbf{D})$
$\displaystyle=$
$\displaystyle\left\\{\prod_{i=1}^{n}\prodi_{s=0}^{\infty}\left[dA_{i}^{W}(s;w_{1},w_{2})\right]^{dN_{i}^{W}(s;w_{1},w_{2})}\right\\}\times$
$\displaystyle\exp\left\\{-\sum_{i=1}^{n}\int_{0}^{\infty}dA_{i}^{W}(s;w_{1},w_{2})\right\\};$
$\displaystyle\mathcal{L}^{V}(\eta,\theta^{V};v_{1},v|\mathbf{D})$
$\displaystyle=$
$\displaystyle\left\\{\prod_{i=1}^{n}\prodi_{s=0}^{\infty}\left[dA_{i}^{V}(s;v_{1},v)\right]^{dN_{i}^{V}(s;v_{1},v)}\right\\}\times$
$\displaystyle\exp\left\\{-\sum_{i=1}^{n}\int_{0}^{\infty}dA_{i}^{V}(s;v_{1},v)\right\\}.$

Let $0=S_{0}<S_{1}<S_{2}<\ldots<S_{K}<S_{K+1}=\infty$ be the ordered distinct
times of any type of event occurrence for all the $n$ sample units. Also, let
$0=T_{0}<T_{1}<T_{2}<\ldots<T_{L}<T_{L+1}=\infty$ be the ordered distinct
values of the set
$\\{\mathcal{E}_{iq}(S_{j}):i=1,\ldots,n;q\in\mathfrak{I}_{Q};j=0,1,\ldots,S_{K}\\}$.
Recall that $\tau_{i}^{*}=\tau_{i}\wedge\tau_{iA}$. Observe that both
$\\{S_{k}:k=0,1,\ldots,K,K+1\\}$ and $\\{T_{l}:l=0,1,\ldots,L,L+1\\}$
partition $[0,\infty)$. For each $i=1,\ldots,n$, and $q\in\mathfrak{I}_{Q}$,
$\mathcal{E}_{iq}(\cdot)$ is observed, hence defined, only on
$[0,\tau_{i}^{*})$. However, for notational convenience, we define
$\mathcal{E}_{iq}(s)=0$ for $s>\tau_{i}^{*}$. In addition, on each non-empty
interval $(S_{j-1}\wedge\tau_{i}^{*},S_{j}\wedge\tau_{i}^{*}]$,
$\mathcal{E}_{iq}(\cdot)$ has an inverse which will be denoted by
$\mathcal{E}_{iqj}^{-1}(\cdot)$. Henceforth, for brevity of notation, we adopt
the mathematically imprecise convention that $0/0=0$.

######  Proposition 4.1.

For $q\in\mathfrak{I}_{Q}$, if $(\alpha_{q},\theta^{R})$ is known, then the
nonparametric maximum likelihood estimator (NPMLE) of $\Lambda_{0q}(\cdot)$ is
given by

$\displaystyle\hat{\Lambda}_{0q}(t;\alpha_{q},\theta^{R})=\sum_{l:T_{l}\leq
t}\left[\frac{\sum_{i=1}^{n}\sum_{j=1}^{K}I\\{\mathcal{E}_{iq}(S_{j})=T_{l}\\}dN_{i}^{R}(S_{j};q)}{S_{q}^{0R}(T_{l}|\alpha_{q},\theta^{R})}\right]$
(12)

where

$\displaystyle S_{q}^{0R}(u|\alpha_{q},\theta^{R})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{K}\left\\{\left[\frac{\rho_{q}[N_{i}^{R}(\mathcal{E}_{iqj}^{-1}(u)-);\alpha_{q}]\psi_{i}^{R}(\mathcal{E}_{iqj}^{-1}(u)-;\theta^{R})}{\mathcal{E}_{iq}^{\prime}[\mathcal{E}_{iqj}^{-1}(u)]}\right]\times\right.$
$\displaystyle\left.I\\{\mathcal{E}_{iq}(S_{j-1}\wedge\tau_{i}^{*})<u\leq\mathcal{E}_{iq}(S_{j}\wedge\tau_{i}^{*})\\}\right\\}.$

###### Proof.

The likelihood
$\mathcal{L}_{q}^{R}(\Lambda_{0q},\alpha_{q},\theta^{R}|\mathbf{D})$ could be
written as follows:

$\displaystyle\mathcal{L}_{q}^{R}$ $\displaystyle=$
$\displaystyle\left[\prod_{i=1}^{n}\prod_{j=1}^{K}[Y_{i}(S_{j})\lambda_{0q}[\mathcal{E}_{iq}(S_{j})\rho_{q}[N_{i}^{R}(S_{j}-);\alpha_{q}]\psi_{i}(S_{j}-;\theta^{R})]^{dN_{i}^{R}(S_{j};q)}\right]\times$
$\displaystyle\exp\left\\{-\sum_{i=1}^{n}\sum_{j=1}^{K+1}\int_{S_{j-1}}^{S_{j}}Y_{i}(s)\lambda_{0q}[\mathcal{E}_{iq}(s)]\rho_{q}[N_{i}^{R}(s-);\alpha_{q}]\psi_{i}(s-;\theta^{R})ds\right\\}.$

Focusing on the nonparametric parameter $\Lambda_{0q}(\cdot)$, the first term
of $\mathcal{L}_{q}^{R}$ could be written as

$\displaystyle\prod_{i=1}^{n}\prod_{j=1}^{K}[\lambda_{0q}[\mathcal{E}_{iq}(S_{j})]^{dN_{i}^{R}(S_{j};q)}$
$\displaystyle=$
$\displaystyle\prod_{l=1}^{L}[\lambda_{0q}(T_{l})]^{dN_{\bullet}^{R}(\infty,T_{l};q)}=\prod_{l=1}^{L}[d\Lambda_{0q}(T_{l})]^{N_{\bullet}^{R}(\infty,T_{l};q)}$

where

$dN_{\bullet}^{R}(\infty,T_{l};q)=\sum_{i=1}^{n}\sum_{j=1}^{K}I\\{\mathcal{E}_{iq}(S_{j})=T_{l}\\}dN_{i}^{R}(S_{j};q).$

The exponent in the second term of $\mathcal{L}_{q}^{R}$ could be written as
follows, the second equality obtained after an obvious change-of-variable and
using the definition of $S_{q}^{0R}(\cdot;\cdot,\cdot)$ in the proposition:

$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{K+1}\int_{S_{j-1}}^{S_{j}}Y_{i}(s)\lambda_{0q}[\mathcal{E}_{iq}(s)]\rho_{q}[N_{i}^{R}(s-);\alpha_{q}]\psi_{i}(s-;\theta^{R})ds$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{K+1}\int_{S_{j-1}\wedge\tau_{i}^{*}}^{S_{j}\wedge\tau_{i}^{*}}\lambda_{0q}[\mathcal{E}_{iq}(s)]\rho_{q}[N_{i}^{R}(s-);\alpha_{q}]\psi_{i}(s-;\theta^{R})ds$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}S_{q}^{0R}(u|\alpha_{q},\theta^{R})d\Lambda_{0q}(u)=\sum_{l=1}^{L+1}\int_{T_{l-1}}^{T_{l}}S_{q}^{0R}(u;\alpha_{q},\theta^{R})d\Lambda_{0q}(u)$
$\displaystyle=$
$\displaystyle\sum_{l=1}^{L}S_{q}^{0R}(T_{l}|\alpha_{q},\theta^{R})d\Lambda_{0q}(T_{l})+\sum_{l=1}^{L+1}\int_{u\in(T_{l-1},T_{l})}S_{q}^{0R}(u|\alpha_{q},\theta^{R})d\Lambda_{0q}(u)$

Therefore, $\mathcal{L}_{q}^{R}$, when viewed solely in terms of the parameter
$\Lambda_{0q}$ equals

$\displaystyle\mathcal{L}_{q}^{R}$ $\displaystyle=$
$\displaystyle\prod_{l=}^{L}[d\Lambda_{0q}(T_{l})]^{dN_{\bullet}^{R}(\infty,T_{l};q)}\times$
$\displaystyle\exp\left\\{-\left[\sum_{l=1}^{L}S_{q}^{0R}(T_{l}|\alpha_{q},\theta^{R})d\Lambda_{0q}(T_{l})+\sum_{l=1}^{L+1}\int_{u\in(T_{l-1},T_{l})}S_{q}^{0R}(u|\alpha_{q},\theta^{R})d\Lambda_{0q}(u)\right]\right\\}$

Since
$\sum_{l=1}^{L+1}\int_{u\in(T_{l-1},T_{l})}S_{q}^{0R}(u|\alpha_{q},\theta^{R})d\Lambda_{0q}(u)\geq
0$, then we obtain the upper bound for $\mathcal{L}_{q}^{R}$ by setting this
term to be equal to zero:

$\mathcal{L}_{q}^{R}\leq\left[\prod_{l=1}^{L}[d\Lambda_{0q}(T_{l})]^{dN_{\bullet}^{R}(\infty,T_{l};q)}\right]\exp\left\\{-\sum_{l=1}^{L}S_{q}^{0R}(T_{l}|\alpha_{q},\theta^{R})d\Lambda_{0q}(T_{l})\right\\}.$

The upper bound is maximized by setting

$d\hat{\Lambda}_{0q}(T_{l}|\alpha_{q},\theta^{R})=\frac{dN_{\bullet}^{R}(\infty,T_{l};q)}{S_{q}^{0R}(T_{l}|\alpha_{q},\theta^{R})},l=1,2,\ldots,L.$

For $u\in(T_{l-1},T_{l})$, we then take
$\hat{\Lambda}_{0q}(u|\alpha_{q},\theta^{R})=\hat{\Lambda}_{0q}(T_{l-1}|\alpha_{q},\theta^{R})$
which will satisfy the condition
$\int_{u\in(T_{l-1},T_{l})}S_{q}^{0R}(u|\alpha_{q},\theta^{R})d\hat{\Lambda}_{0q}(u)=0$
for all $l=1,2,\ldots,L+1$. Thus,

$\hat{\Lambda}_{0q}(t|\alpha_{q},\theta^{R})=\sum_{l:T_{l}\leq
t}d\hat{\Lambda}_{0q}(T_{l};\alpha_{q},\theta^{R})=\sum_{l:T_{l}\leq
t}\frac{dN_{\bullet}^{R}(\infty,T_{l};q)}{S_{q}^{0R}(T_{l}|\alpha_{q},\theta^{R})},$

which is a step-function with jumps only on $T_{l}$s with
$dN_{\bullet}^{R}(\infty,T_{l};q)>0$, maximizes
$\Lambda_{0q}(\cdot)\mapsto\mathcal{L}_{q}^{R}(\Lambda_{0q}(\cdot)|\alpha_{q},\theta^{R}|\mathbf{D})$
for given $(\alpha_{q},\theta^{R})$, completing the proof of the proposition.
∎

A more elegant representation of the Aalen-Breslow-Nelson type estimator
$\hat{\Lambda}_{0q}(\cdot;\alpha_{q},\theta^{R})$ in Proposition 4.1, which
shows that the estimator is also moment-based, aside from being useful in
obtaining finite and asymptotic properties, is through the use doubly-indexed
processes as in [25] for a setting with only one recurrent event type and
without LM and HS processes. Define the doubly-indexed processes
$\\{(N_{i}^{R}(s,t;q),A_{i}^{R}(s,t;q|\Lambda_{0q},\alpha_{1},\theta^{R}):(s,t)\in\Re_{+}^{2}\\}$
where

$\displaystyle N_{i}^{R}(s,t;q)=\int_{0}^{s}I\\{\mathcal{E}_{iq}(v)\leq
t\\}dN_{i}^{R}(v;q);$ $\displaystyle
A_{i}^{R}(s,t;q|\Lambda_{0q},\alpha_{q},\theta^{R})=\int_{0}^{s}I\\{\mathcal{E}_{iq}(v)\leq
t\\}dA_{i}^{R}(v;q|\Lambda_{0q},\alpha_{q},\theta^{R}).$

Also, let

$\displaystyle N_{\bullet}^{R}(s,t;q)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}N_{i}^{R}(s,t;q)\ \mbox{and}\
A_{\bullet}^{R}(s,t;q|\Lambda_{0q},\alpha_{q},\theta^{R})=\sum_{i=1}^{n}A_{i}^{R}(s,t;q|\Lambda_{0q},\alpha_{q},\theta^{R}).$

Then, for fixed $t$,
$\\{M_{\bullet}^{R}(s,t;q|\Lambda_{0q},\alpha_{q},\theta^{R})=N_{\bullet}^{R}(s,t;q)-A_{\bullet}^{R}(s,t;q|\Lambda_{0q},\alpha_{q},\theta^{R}):s\geq
0\\}$ is a zero-mean square-integrable martingale. Thus,
$E[N_{\bullet}^{R}(s,t;q)]=E[A_{\bullet}^{R}(s,t;q|\Lambda_{0q},\alpha_{q},\theta^{R})]$.

######  Proposition 4.2.

For $q\in\mathfrak{I}_{Q}$,
$A_{\bullet}^{R}(s,t;q|\Lambda_{0q},\alpha_{q},\theta^{R})=\int_{0}^{t}S_{q}^{0R}(s,u|\alpha_{q},\theta^{R})d\Lambda_{0q}(u)$,
where

$\displaystyle S_{q}^{0R}(s,u|\alpha_{q},\theta^{R})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{K}\left\\{\left[\frac{\rho_{q}[N_{i}^{R}(\mathcal{E}_{iqj}^{-1}(u)-);\alpha_{q}]\psi_{i}^{R}(\mathcal{E}_{iqj}^{-1}(u)-;\theta^{R})}{\mathcal{E}_{iq}^{\prime}[\mathcal{E}_{iqj}^{-1}(u)]}\right]\times\right.$
$\displaystyle\left.I\\{\mathcal{E}_{iq}(S_{j-1}\wedge\tau_{i}^{*}\wedge
s)<u\leq\mathcal{E}_{iq}(S_{j}\wedge\tau_{i}^{*}\wedge s)\\}\right\\};$

and $dN_{\bullet}^{R}(s,T_{l};q)=\sum_{i=1}^{n}\sum_{j=1}^{K}I\\{S_{j}\leq
s;\mathcal{E}_{iq}(S_{j})=T_{l}\\}dN_{i}^{R}(S_{j};q)$.

###### Proof.

Similar to steps in the proof of Proposition 4.1. ∎

By first assuming that $(\alpha_{q},\theta^{R})$ is known, then from the
identities in Proposition 4.2, a method-of-moments type estimator of
$\Lambda_{0q}(t)$ is given by

$\hat{\Lambda}_{0q}(s,t|\alpha_{q},\theta^{R})=\sum_{l:T_{l}\leq
t}\frac{dN_{\bullet}^{R}(s,T_{l};q)}{S_{q}^{0R}(s,T_{l}|\alpha_{q},\theta^{R})}=\int_{0}^{t}\frac{N_{\bullet}^{R}(s,du;q)}{S_{q}^{0R}(s,u|\alpha_{q},\theta^{R})}.$
(14)

When $s\rightarrow\infty$, this
$\hat{\Lambda}_{0q}(s,t|\alpha_{q},\theta^{R})$ converges to the estimator
$\hat{\Lambda}_{0q}(t|\alpha_{q},\theta^{R})$ in Proposition 4.1. Next, we
obtain estimators of the $\eta(w_{1},w_{2})$s and $\xi(v_{1},v)$s, again by
first assuming first that $\theta^{W}$ and $\theta^{V}$ are known.

######  Proposition 4.3.

If $(\theta^{W},\theta^{V})$ are known, the ML estimators of the
$\eta(w_{1},w_{2})$s and $\xi(v_{1},v)$s are the “occurrence-exposure” rates

$\displaystyle\hat{\eta}(w_{1},w_{2}|\theta^{W})$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}\sum_{j=1}^{K}dN_{i}^{W}(S_{j};w_{1},w_{2})}{S^{0W}(w_{1};\theta^{W})},\forall(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}};$
(15) $\displaystyle\hat{\xi}(v_{1},v|\theta^{V})$ $\displaystyle=$
$\displaystyle\frac{\sum_{i=1}^{n}\sum_{j=1}^{K}dN_{i}^{V}(S_{j};v_{1},v)}{S^{0V}(v_{1};\theta^{V})},\forall(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}},$
(16)

where

$S^{0W}(w_{1};\theta^{W})=\int_{0}^{\infty}\sum_{i=1}^{n}Y_{i}(s)Y_{i}^{W}(s;w_{1})\psi_{i}^{W}(s-;\theta^{W})ds;$
$S^{0V}(v_{1};\theta^{V})=\int_{0}^{\infty}\sum_{i=1}^{n}Y_{i}(s)Y_{i}^{V}(s;v_{1})\psi_{i}^{V}(s-;\theta^{V})ds.$

###### Proof.

Follows immediately by maximizing the likelihood functions $\mathcal{L}^{W}$
and $\mathcal{L}^{V}$ with respect to the $\eta(w_{1},w_{2})$s and
$\xi(v_{1},v)$s, respectively. ∎

We can now form the profile likelihoods for the parameters
$(\\{\alpha_{q},q\in\mathfrak{I}_{Q}\\},\theta^{R},\theta^{W},\theta^{V})$.
These are the likelihoods that are obtained after plugging-in the ‘estimators’
$\hat{\Lambda}_{0q}(\cdot;\alpha_{q},\theta^{R})$s,
$\hat{\eta}(w_{1},w_{2})$s, and $\hat{\xi}(v_{1},v)$s in the full likelihoods.
The resulting profile likelihoods are reminiscent of the partial likelihood
function in Cox’s proportional hazards model [8, 2].

######  Proposition 4.4.

The three profile likelihood functions $\mathcal{L}_{pl}^{R}$,
$\mathcal{L}_{pl}^{W}$ and $\mathcal{L}_{pl}^{V}$ are given by

$\displaystyle\mathcal{L}^{R}_{pl}(\alpha_{q},q\in\mathfrak{I}_{Q},\theta^{R}|\mathbf{D})=$
$\displaystyle\prod_{q\in\mathfrak{I}_{Q}}\prod_{i=1}^{n}\prod_{j=1}^{K}\prod_{l=1}^{L}\left[\frac{\rho_{q}[N_{i}^{R}(S_{j}-);\alpha_{q}]\psi_{i}^{R}(S_{j}-;\theta^{R})}{S_{q}^{0R}(T_{l}|\alpha_{q},\theta^{R})}\right]^{I\\{\mathcal{E}_{iq}(S_{j})=T_{l}\\}dN_{i}^{R}(S_{j};q)};$
$\mathcal{L}_{pl}^{W}(\theta^{W}|\mathbf{D})=\prod_{w_{1}\in\mathfrak{W}}\prod_{i=1}^{n}\prod_{j=1}^{K}\left[\frac{\psi_{i}^{W}(S_{j}-;\theta^{W})}{S^{0W}(w_{1};\theta^{W})}\right]^{dN_{i}^{W}(S_{j};w_{1},\bullet)};$
$\mathcal{L}_{pl}^{V}(\theta^{V}|\mathbf{D})=\prod_{v_{1}\in\mathfrak{V}_{1}}\prod_{i=1}^{n}\prod_{j=1}^{K}\left[\frac{\psi_{i}^{V}(S_{j}-;\theta^{V})}{S^{0V}(v_{1};\theta^{V})}\right]^{dN_{i}^{V}(S_{j};v_{1},\bullet)}.$

with $dN_{i}^{W}(S_{j};w_{1},\bullet)=\sum_{w_{2}\in\mathfrak{W};\ w_{2}\neq
w_{1}}dN_{i}^{W}(S_{j};w_{1},w_{2})$, the number of transitions from state
$w_{1}$ at time $S_{j}$ for unit $i$, and
$dN_{i}^{V}(S_{j};v_{1},\bullet)=\sum_{v\in\mathfrak{V};\ v\neq
v_{1}}dN_{i}^{V}(S_{j};v_{1},v)$, the number of transitions from state $v_{1}$
at time $S_{j}$ for unit $i$.

###### Proof.

These follow immediately by plugging-in the ‘estimators’ into the three main
likelihoods in (4.2) and then simplifying. ∎

From these three profile likelihoods, we could obtain estimators of the
parameters $\alpha_{q}$s, $\theta^{R}$, $\theta^{W}$, and $\theta^{V}$ as
follows:

$\displaystyle(\hat{\alpha}_{q},q\in\mathfrak{I}_{Q},\hat{\theta}^{R})={\arg\max}_{(\alpha_{q},\theta^{R})}\mathcal{L}^{R}_{pl}(\alpha_{q},q\in\mathfrak{I}_{Q},\theta^{R}|\mathbf{D});$
$\displaystyle\hat{\theta}^{W}=\arg\max_{\theta^{W}}\mathcal{L}_{pl}^{W}(\theta^{W}|\mathbf{D})\
\mbox{and}\
\hat{\theta}^{V}=\arg\max_{\theta^{V}}\mathcal{L}_{pl}^{V}(\theta^{V}|\mathbf{D}).$

Equivalently, these estimators are maximizers of the logarithm of the profile
likelihoods. These log-profile likelihoods are more conveniently expressed in
terms of stochastic integrals as follows:

$\displaystyle l_{pl}^{R}$ $\displaystyle=$
$\displaystyle\sum_{q\in\mathfrak{I}_{Q}}\sum_{i=1}^{n}\int_{0}^{\infty}\int_{0}^{\infty}\left[\log\rho_{q}[N_{i}^{R}(s-);\alpha_{q}]+\log\psi_{i}^{R}(s-;\theta^{R})-\right.$
$\displaystyle\left.\log
S_{q}^{0R}(t;\alpha_{q},\theta^{R})\right]N_{i}^{R}(ds,dt;q);$ $\displaystyle
l_{pl}^{W}$ $\displaystyle=$
$\displaystyle\sum_{w_{1}\in\mathfrak{W}}\sum_{i=1}^{n}\int_{0}^{\infty}\left[\log\psi_{i}^{W}(s-;\theta^{W})-\log
S^{0W}(s;w_{1}|\theta^{W})\right]N_{i}^{W}(ds;w_{1},\bullet);$ $\displaystyle
l_{pl}^{V}$ $\displaystyle=$
$\displaystyle\sum_{v_{1}\in\mathfrak{V}_{1}}\sum_{i=1}^{n}\int_{0}^{\infty}\left[\log\psi_{i}^{V}(s-;\theta^{V})-\log
S^{0V}(s;v_{1}|\theta^{V})\right]N_{i}^{V}(ds;v_{1},\bullet).$

Associated with each of these log-profile likelihood functions are their
profile score vector (the gradient or vector of partial derivatives) and
profile observed information matrix (negative of the matrix of second partial
derivatives): $U^{R}(\alpha,\theta^{R})$; $U^{W}(\theta^{W})$ and
$U^{V}(\theta^{V})$ for the score vectors, and $I^{R}(\alpha,\theta^{R})$,
$I^{W}(\theta^{W})$ and $I^{V}(\theta^{V})$ for the observed information
matrices. The estimators could then be obtained as the solutions of the set of
equations

$U_{pl}^{R}(\hat{\alpha},\hat{\theta}^{R})=0;U_{pl}^{W}(\hat{\theta}^{W})=0;U_{pl}^{V}(\hat{\theta}^{V})=0.$

A possible computational approach to obtaining the estimates is via the
Newton-Raphson iteration procedure:

$\displaystyle(\hat{\alpha},\hat{\theta}^{R})\leftarrow(\hat{\alpha},\hat{\theta}^{R})+I^{R}(\hat{\alpha},\hat{\theta}^{R})^{-1}U^{R}(\hat{\alpha},\hat{\theta}^{R});$
$\displaystyle\hat{\theta}^{W}\leftarrow\hat{\theta}^{W}+I^{W}(\hat{\theta}^{W})^{-1}U^{W}(\hat{\theta}^{W});\
\hat{\theta}^{V}\leftarrow\hat{\theta}^{V}+I^{V}(\hat{\theta}^{V})^{-1}U^{V}(\hat{\theta}^{V}).$

Having obtained the estimates of $\alpha$, $\theta^{W}$ and $\theta^{V}$, we
plug them into $\hat{\Lambda}_{0q}(t|\alpha_{q},\theta^{R})$s,
$\eta(w_{1},w_{2}|\theta^{W})$s and $\xi(v_{1},v|\theta^{V})$s to obtain the
estimators $\hat{\Lambda}_{0q}(t)$s, $\hat{\eta}(w_{1},w_{2})$s and
$\hat{\xi}(v_{1},v)$s:

$\hat{\Lambda}_{0q}(t)=\hat{\Lambda}_{0q}(t|\hat{\alpha}_{q},\hat{\theta}^{R});\quad\hat{\eta}(w_{1},w_{2})=\hat{\eta}(w_{1},w_{2}|\hat{\theta}^{W});\quad\hat{\xi}(v_{1},v)=\hat{\xi}(v_{1},v|\hat{\theta}^{V}).$
(17)

##  5 Asymptotic Properties of Estimators

In this section we provide some asymptotic properties of the estimators,
though we do not present the rigorous proofs of the results due to space
constraints and instead defer them to a separate paper. To make our exposition
more concise and compact, we introduce additional notation. Consider a real-
valued function $h$ defined on $\mathfrak{I}_{\mathfrak{W}}$. Then,
$h(w_{1},w_{2})$ will represent the value at $(w_{1},w_{2})$, but when we
simply write $h$ it means the $|\mathfrak{I}_{\mathfrak{W}}|\times 1$ vector
consisting of $h(w_{1},w_{2}),(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}$.
Thus, $\eta$ is an $|\mathfrak{I}_{\mathfrak{W}}|\times 1$ (column) vector
with elements $\eta(w_{1},w_{2})$s; $N_{i}^{W}(s)$ is an
$|\mathfrak{I}_{\mathfrak{W}}|\times 1$ vector consisting of
$N_{i}^{W}(s;w_{1},w_{2})$s; $S^{0W}(s|\theta^{W})$ is an
$|\mathfrak{I}_{\mathfrak{W}}|\times 1$ vector with elements
$S^{0W}(s;w_{1}|\theta^{W}),(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}$; and
$\psi_{i}^{W}(s|\theta^{W})$ is an $|\mathfrak{I}_{\mathfrak{W}}|\times 1$
vector consisting of the same elements. Similarly for those associated with
the HS process where functions are defined over $\mathfrak{I}_{\mathfrak{V}}$;
e.g., $\eta=(\eta(v_{1},v),(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}})$,
$N_{i}^{V}(s)=(N_{i}^{V}(s;v_{1},v):(v_{1},v)\in\mathfrak{I}_{\mathfrak{V}})$,
etc.

Let us first consider the profile likelihood function $\mathcal{L}_{pl}^{W}$
and the estimators $\hat{\theta}^{W}$ and $\hat{\eta}$. Using the above
notation, the associated profile log-likelihood function, score function, and
observed information matrix could be written as follows:

$\displaystyle
l_{pl}^{W}(\theta^{W})=\sum_{i=1}^{n}\int_{0}^{\infty}\left[\log\psi_{i}^{W}(s-|\theta^{W})-\log
S^{0W}(s|\theta^{W})\right]^{\tt T}dN_{i}^{W}(s);$ $\displaystyle
U_{pl}^{W}(\theta^{W})=\sum_{i=1}^{n}\int_{0}^{\infty}H_{1i}^{W}(s|\theta^{W})^{\tt
T}dN_{i}^{W}(s);$ $\displaystyle
I_{pl}^{W}(\theta^{W})=\sum_{i=1}^{n}\int_{0}^{\infty}V_{1}^{W}(s|\theta^{W})^{\tt
T}dN_{i}^{W}(s),$

where

$\displaystyle H_{1i}^{W}(s|\theta^{W})$ $\displaystyle=$
$\displaystyle\frac{\stackrel{{\scriptstyle\cdot}}{{\psi}}_{i}^{W}(s-|\theta^{W})}{\psi_{i}^{W}(s-|\theta^{W})}-\frac{\stackrel{{\scriptstyle\cdot}}{{S}}^{0W}(s|\theta^{W})}{S^{0W}(s|\theta^{W})}$
$\displaystyle\equiv$
$\displaystyle\left(\frac{\stackrel{{\scriptstyle\cdot}}{{\psi}}_{i}^{W}(s-|\theta^{W})}{\psi_{i}^{W}(s-|\theta^{W})}-\frac{\stackrel{{\scriptstyle\cdot}}{{S}}^{0W}(s;w_{1}|\theta^{W})}{S^{0W}(s;w_{1}|\theta^{W})}:(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}\right);$

and

$\displaystyle V_{1}^{W}(s|\theta^{W})$ $\displaystyle=$
$\displaystyle\frac{\stackrel{{\scriptstyle\cdot\cdot}}{{S}}^{0W}(s|\theta^{W})}{S^{0W}(s|\theta^{W})}-\left(\frac{\stackrel{{\scriptstyle\cdot}}{{S}}^{0W}(s|\theta^{W})}{S^{0W}(s|\theta^{W})}\right)^{\otimes
2}$ $\displaystyle\equiv$
$\displaystyle\left(\frac{\stackrel{{\scriptstyle\cdot\cdot}}{{S}}^{0W}(s;w_{1}|\theta^{W})}{S^{0W}(s;w_{1}|\theta^{W})}-\left(\frac{\stackrel{{\scriptstyle\cdot}}{{S}}^{0W}(s;w_{1}|\theta^{W})}{S^{0W}(s;w_{1}|\theta^{W})}\right)^{\otimes
2}:(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}\right).$

A ‘$\cdot$’ over a function means gradient with respect to the parameter
vector, while a ‘$\cdot\cdot$’ over a function means the matrix of second-
partial derivatives with respect to the element of the parameter vector. In
obtaining $I_{pl}^{W}$, we also used the fact that $\psi_{i}^{W}$ is an
exponential function. We now present some results about $\hat{\theta}^{W}$. We
also let

$\mathfrak{I}_{pl}^{W}=\mbox{pr-
lim}\left[\frac{1}{n}I_{pl}^{W}(\theta^{W})\right].$

This will be a function of all the model parameters, since the limiting
behavior of the $Y_{i}$s and $Y_{i}^{W}$s will depend on all the model
parameters, owing to the interplay among the RCR, LM, and HS processes, as
could be seen in the special case of Poisson processes and CTMCs.

######  Proposition 5.1.

Under certain regularity conditions, we have, as $n\rightarrow\infty$,

  * (i)

(Consistency) $\hat{\theta}^{W}\stackrel{{\scriptstyle
p}}{{\rightarrow}}\theta^{W}$;

  * (ii)

(Asymptotic Representation)

$\displaystyle\sqrt{n}[\hat{\theta}^{W}-\theta^{W}]=\left[\frac{1}{n}I_{pl}^{W}(\theta^{W})\right]^{-1}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\int_{0}^{\infty}H_{1i}^{W}(s|\theta^{W})^{\tt
T}dM_{i}^{W}(s|\eta,\theta^{W})+o_{p}(1);$

  * (iii)

(Asymptotic Normality)
$\sqrt{n}[\hat{\theta}^{W}-\theta^{W}]\stackrel{{\scriptstyle
d}}{{\rightarrow}}N\left[0,\left(\mathfrak{I}_{pl}^{W}\right)^{-1}\right].$

We point out an important result needed for the proof of Proposition 5.1(iii),
which is that

$\displaystyle\frac{1}{n}I_{pl}^{W}(\theta^{W})=$
$\displaystyle\sum_{(w_{1},w_{2})\in\mathfrak{I}_{\mathfrak{W}}}\left[\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}\left\\{H_{1i}^{W}(s;w_{1}|\theta^{W})\right\\}^{\otimes
2}Y_{i}(s)Y_{i}^{W}(s;w_{1})\eta(w_{1},w_{2})\psi_{i}^{W}(s-|\theta^{W})ds\right]$
$\displaystyle+o_{p}(1)\stackrel{{\scriptstyle
p}}{{\rightarrow}}\mathfrak{I}_{pl}^{W}.$

This asymptotic equivalence indicates where the involvement of the at-risk
processes come into play in the limiting profile information matrix, hence the
dependence on all the model parameters. This also shows that a natural
consistent estimator of $\mathfrak{I}_{pl}^{W}$ is
$I_{pl}^{W}(\hat{\theta}^{W})/n$.

Next, we are now in position to present asymptotic results about $\hat{\eta}$.
First, let $s^{0W}$ and $\stackrel{{\scriptstyle\cdot}}{{s}}^{0W}$ be
deterministic functions satisfying

$\displaystyle\left|\frac{1}{n}\int_{0}^{\infty}S^{0W}(z|\theta^{W})dz-\int_{0}^{\infty}s^{0W}(z)dz\right|\rightarrow
0;$
$\displaystyle\left|\frac{1}{n}\int_{0}^{\infty}\stackrel{{\scriptstyle\cdot}}{{S}}^{0W}(z|\theta^{W})dz-\int_{0}^{\infty}\stackrel{{\scriptstyle\cdot}}{{s}}^{0W}(z)dz\right|\rightarrow
0.$

######  Proposition 5.2.

Under certain regularity conditions, we have as $n\rightarrow\infty$,

  * (i)

(Consistency) $\hat{\eta}\stackrel{{\scriptstyle p}}{{\rightarrow}}\eta$;

  * (ii)

(Asymptotic Normality) $\sqrt{n}(\hat{\eta}-\eta)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\Sigma)$, where

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\mbox{Dg}\left(\int_{0}^{\infty}s^{0W}(z)dz\right)^{-1}\mbox{Dg}(\eta)+$
$\displaystyle\left[\mbox{Dg}\left(\int_{0}^{\infty}s^{0W}(z)dz\right)^{-1}\mbox{Dg}(\eta)\left(\int_{0}^{\infty}\stackrel{{\scriptstyle\cdot}}{{s}}^{0W}(z)dz\right)^{\tt
T}(\mathfrak{I}_{pl}^{W})^{-1/2}\right]^{\otimes 2}.$

  * (iii)

(Joint Asymptotic Normality) $\sqrt{n}(\hat{\eta}-\eta)$ and
$\sqrt{n}(\hat{\theta}^{W}-\theta^{W})$ are jointly asymptotically normal with
means zeros and asymptotic covariance matrix

$\displaystyle\mbox{Acov}(\sqrt{n}(\hat{\theta}^{W}-\theta^{W}),\sqrt{n}(\hat{\eta}-\eta))=$
$\displaystyle-(\mathfrak{I}_{pl}^{W})^{-1}\left(\int_{0}^{\infty}\stackrel{{\scriptstyle\cdot}}{{s}}^{0W}(z)dz\right)\mbox{Dg}(\eta)\mbox{Dg}\left(\int_{0}^{\infty}s^{0W}(z)dz\right)^{-1}.$

We remark that in result (ii) for the asymptotic covariance matrix $\Sigma$ in
Proposition 5.2, the additional variance term in the right-hand side is the
effect of plugging-in the estimator $\hat{\theta}^{W}$ for $\theta^{W}$ in the
‘estimators’ $\hat{\eta}(w_{1},w_{2}|\theta^{W})$s to obtain
$\hat{\eta}(w_{,}w_{2})$s. Without having to write them down explicitly,
similar results are obtainable for the estimators $\hat{\theta}^{V}$ and
$\hat{\xi}$.

Next, we present results concerning the asymptotic properties of the
estimators of $\alpha_{q}$s, $\theta^{R}$, and $\Lambda_{0q}(\cdot)$s. Define
the restricted profile likelihood for $(\alpha_{q},\theta^{R})$ based on data
$\mathbf{D}$ observed over $[0,s^{*}]$ via

$\displaystyle\mathcal{L}_{pl}^{R}(\alpha,\theta^{R})$ $\displaystyle\equiv$
$\displaystyle\mathcal{L}_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})\equiv\mathcal{L}_{pl}^{R}(\alpha,\theta^{R}|\mathbf{D}(s^{*},t^{*}))$
$\displaystyle=$
$\displaystyle\prod_{q=1}^{Q}\prod_{i=1}^{n}\prodi_{s=0}^{s^{*}}\left[\frac{\rho_{q}[N_{i}^{R}(s-);\alpha_{q}]\psi_{i}^{R}(s-;\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}\right]^{I\\{\mathcal{E}_{iq}(s)\leq
t^{*}\\}dN_{i}^{R}(s;q)}$

with $S_{q}^{0R}(\cdot,\cdot|\cdot,\cdot)$ as defined in Proposition 4.2. This
is restricted in the sense that we only consider events that happened when the
effective age are no more than $t^{*}$ and happened over $[0,s^{*}]$. Note
that as we let $t^{*}\rightarrow\infty$ and $s^{*}\rightarrow\infty$, we
obtain the profile likelihood in Proposition 4.4, The log-profile likelihood
function is

$\displaystyle{l}_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})\equiv\log\mathcal{L}_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})$
$\displaystyle=$
$\displaystyle\sum_{q=1}^{Q}\sum_{i=1}^{n}\int_{0}^{s^{*}}\left\\{\log\rho_{q}[N_{i}^{R}(s-);\alpha_{q}]+\log\psi_{i}^{R}(s-|\theta^{R})-\log
S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})\right\\}\times$
$\displaystyle I\\{\mathcal{E}_{iq}(s)\leq t^{*}\\}dN_{i}^{R}(ds;q).$

The associated profile score function is

$U_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})=U_{pl}^{R}(\alpha,\theta^{R})=\sum_{i=1}^{n}\sum_{q=1}^{Q}\int_{0}^{s^{*}}H_{iq}^{R}(s|s^{*},\alpha,\theta^{R})N_{i}^{R}(ds,t^{*};q)$

where, for $q\in\mathfrak{I}_{Q}$,

$H_{iq}^{R}(s|s^{*},\alpha,\theta^{R})^{\tt T}=\left[0^{\tt T},\ldots,0^{\tt
T},H_{i1q}(s|s^{*},\alpha,\theta^{R})^{\tt T},0^{\tt T},\ldots,0^{\tt
T},H_{i2q}(s|s^{*},\alpha,\theta^{R})^{\tt T}\right],$

a $(Q+1)\times 1$ block matrix and with

$H_{i1q}^{R}(s|s^{*},\alpha_{q},\theta^{R})=\frac{\nabla_{\alpha_{q}}\rho_{q}[N_{i}^{R}(s-);\alpha_{q})}{\rho_{q}[N_{i}^{R}(s-);\alpha_{q})}-\frac{\nabla_{\alpha_{q}}S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})};$
$H_{i2q}^{R}(s|s^{*},\alpha_{q},\theta^{R})=\frac{\nabla_{\theta^{R}}\psi_{i}^{R}(s-|\theta^{R})}{\psi_{i}^{R}(s-|\theta^{R})}-\frac{\nabla_{\theta^{R}}S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}.$

Note the dimensions of these block vectors: the $q$th block is of dimension
$\dim(\alpha_{q})\times 1$, while the $(Q+1)$th block is of dimension
$\dim(\theta^{R})\times 1$. An important martingale representation of
$U_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})$, which is straight-forward to
establish, is

$U_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})=\sum_{i=1}^{n}\sum_{q=1}^{Q}\int_{0}^{s^{*}}H_{iq}^{R}(s|s^{*},\alpha,\theta^{R})M_{i}^{R}(ds,t^{*};q).$

The predictable variation associated with this score function is

$\langle
U_{pl}^{R}(\alpha,\theta^{R}|\cdot,t^{*})\rangle(s^{*})=\sum_{i=1}^{n}\sum_{q=1}^{n}\int_{0}^{s^{*}}[H_{iq}^{R}(s|s^{*},\alpha,\theta^{R})]^{\otimes
2}A_{i}^{R}(ds,t^{*};q)$

with $A_{i}^{R}(ds,t;q)=I\\{\mathcal{E}_{iq}(s)\leq t^{*}\\}dA_{i}^{R}(s;q)$.
The estimators $\hat{\alpha}_{q}(s^{*},t^{*}),q=1,\ldots,Q,$ and
$\hat{\theta}^{R}(s^{*},t^{*})$ satisfy the equation
$U_{pl}^{R}(\hat{\alpha},\hat{\theta}^{R}|s^{*},t^{*})=0.$

The observed profile information matrix
$I_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})$ is a $(Q+1)\times(Q+1)$ symmetric
block matrix with the following block elements, for $q,q^{\prime}=1,\ldots,Q$:

$I_{pl,qq^{\prime}}^{R}(\alpha,\theta^{R}|s^{*},t^{*})=\left\\{\begin{array}[]{ccc}\sum_{i=1}^{n}\int_{0}^{s^{*}}V_{i11qq}^{R}(s|s^{*},\alpha,\theta^{R})N_{i}^{R}(ds,t^{*};q)&\mbox{for}&q=q^{\prime}\\\
0&\mbox{for}&q\neq q^{\prime}\end{array}\right.;$
$I_{pl,(Q+1)q}^{R}(\alpha,\theta^{R}|s^{*},t^{*})=\sum_{i=1}^{n}\int_{0}^{s^{*}}V_{i21q}^{R}(s|s^{*},\alpha,\theta^{R})N_{i}^{R}(ds,t^{*};q);$
$I_{pl,(Q+1)(Q+1)}^{R}(\alpha,\theta^{R}|s^{*},t^{*})=\sum_{i=1}^{n}\sum_{q=1}^{Q}\int_{0}^{s^{*}}V_{i22q}^{R}(s|s^{*},\alpha,\theta^{R})N_{i}^{R}(ds,t^{*};q),$

where $V_{i11qq^{\prime}}^{R}$ is of dimension
$\dim(\alpha_{q})\times\dim(\alpha_{q^{\prime}})$ for
$q,q^{\prime}=1,2,\ldots,Q$; $V_{i21q}^{R}$ and $(V_{i12q}^{R})^{\tt T}$ have
dimensions $\dim(\theta^{R})\times\dim(\alpha_{q})$; and $V_{i22q}^{R}$ has
dimension $\dim(\theta^{R})\times\dim(\theta^{R})$. These are given by the
following expressions, for $q,q^{\prime}=1,\ldots,Q$:

$\displaystyle
V_{i11qq}^{R}(s|s^{*},\alpha,\theta^{R})=-\left\\{\frac{\nabla_{\alpha_{q}\alpha_{q}^{\tt
T}}\rho_{q}[N_{i}(s-);\alpha_{q}]}{\rho_{q}[N_{i}(s-);\alpha_{q}]}-\left(\frac{\nabla_{\alpha_{q}}\rho_{q}[N_{i}(s-);\alpha_{q}]}{\rho_{q}[N_{i}(s-);\alpha_{q}]}\right)^{\otimes
2}\right\\}+$ $\displaystyle\left\\{\frac{\nabla_{\alpha_{q}\alpha_{q}^{\tt
T}}S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}-\left(\frac{\nabla_{\alpha_{q}}S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}\right)^{\otimes
2}\right\\};$ $V_{i11qq^{\prime}}^{R}(s|s^{*},\alpha,\theta^{R})=0,q\neq
q^{\prime};$ $\displaystyle
V_{i21q}^{R}(s|s^{*},\alpha,\theta^{R})=-\left\\{\frac{\nabla_{\alpha_{q}(\theta^{R})^{\tt
T}}S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}-\right.$
$\displaystyle\left.\left(\frac{\nabla_{\alpha_{q}}S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}\right)\left(\frac{\nabla_{(\theta^{R})^{\tt
T}}S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}\right)\right\\};$
$\displaystyle V_{i22q}^{R}(s|s^{*},\alpha,\theta^{R})$ $\displaystyle=$
$\displaystyle-\left\\{\frac{\nabla_{\theta^{R}(\theta^{R})^{\tt
T}}S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}-\right.$
$\displaystyle\left.\left(\frac{\nabla_{\theta^{R}}S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}{S_{q}^{0R}(s^{*},\mathcal{E}_{iq}(s)|\alpha_{q},\theta^{R})}\right)^{\otimes
2}\right\\},$

where for the last expression we used the fact that $\psi_{i}(z)=\exp(z)$.

A condition needed for the asymptotic results is that there is an invertible
matrix function $\mathfrak{I}_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})$ which
equals the in-probability limit of
$\frac{1}{n}I_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})$ as
$n\rightarrow\infty$. Under this condition, we have the following consistency
and asymptotic normality results for $(\hat{\alpha},\hat{\theta}^{R})$:

######  Proposition 5.3.

Under regularity conditions and as $n\rightarrow\infty$,

  * (i)

$(\hat{\alpha}(s^{*},t^{*}),\hat{\theta}^{R}(s^{*},t^{*}))\stackrel{{\scriptstyle
p}}{{\rightarrow}}(\alpha,\theta^{R})$;

  * (ii)

$\sqrt{n}\left[\begin{array}[]{c}\hat{\alpha}(s^{*},t^{*})-\alpha\\\
\hat{\theta}^{R}(s^{*},t^{*})-\theta^{R}\end{array}\right]\stackrel{{\scriptstyle
d}}{{\rightarrow}}N\left(0,[\mathfrak{I}_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})]^{-1}\right).$

Important equivalences to note are, for $q,q^{\prime}=1,\ldots,Q$:

$\displaystyle\frac{1}{n}I_{pl,qq^{\prime}}^{R}(\alpha,\theta^{R}|s^{*},t^{*})$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{s^{*}}V_{i11qq^{\prime}}^{R}(s|s^{*},\alpha,\theta^{R})I\\{\mathcal{E}_{iq}(s)\leq
t^{*}\\}\times$ $\displaystyle
Y_{i}(s)\rho_{q}[N_{i}^{R}(s-);\alpha_{q}]\psi_{i}^{R}(s-;\theta^{R})\lambda_{0q}[\mathcal{E}_{iq}(s)]ds+o_{p}(1);$
$\displaystyle\frac{1}{n}I_{pl,(Q+1)q}^{R}(\alpha,\theta^{R}|s^{*},t^{*})$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{s^{*}}V_{i21q}^{R}(s|s^{*},\alpha,\theta^{R})I\\{\mathcal{E}_{iq}(s)\leq
t^{*}\\}\times$ $\displaystyle
Y_{i}(s)\rho_{q}[N_{i}^{R}(s-);\alpha_{q}]\psi_{i}^{R}(s-;\theta^{R})\lambda_{0q}[\mathcal{E}_{iq}(s)]ds+o_{p}(1);$
$\displaystyle\frac{1}{n}I_{pl,(Q+1)(Q+1)}^{R}(\alpha,\theta^{R}|s^{*},t^{*})$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\sum_{q=1}^{Q}\int_{0}^{s^{*}}V_{i22q}^{R}(s|s^{*},\alpha,\theta^{R})I\\{\mathcal{E}_{iq}(s)\leq
t^{*}\\}\times$ $\displaystyle
Y_{i}(s)\rho_{q}[N_{i}^{R}(s-);\alpha_{q}]\psi_{i}^{R}(s-;\theta^{R})\lambda_{0q}[\mathcal{E}_{iq}(s)]ds+o_{p}(1).$

In addition, the regularity conditions must imply that, we have

$\displaystyle\left|\frac{1}{n}I_{pl,qq}^{R}(\alpha,\theta^{R}|s^{*},t^{*})-\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{s^{*}}[H_{i1q}^{R}(s|s^{*},\alpha,\theta^{R})]^{\otimes
2}A_{i}^{R}(ds,t^{*};q)\right|\stackrel{{\scriptstyle p}}{{\rightarrow}}0;$
$\displaystyle\left|\frac{1}{n}I_{pl,(Q+1)(Q+1)}^{R}(\alpha,\theta^{R}|s^{*},t^{*})-\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{s^{*}}[H_{i2q}^{R}(s|s^{*},\alpha,\theta^{R})]^{\otimes
2}A_{i}^{R}(ds,t^{*};q)\right|\stackrel{{\scriptstyle p}}{{\rightarrow}}0;$
$\displaystyle\left|\frac{1}{n}I_{pl,(Q+1)q}^{R}(\alpha,\theta^{R}|s^{*},t^{*})-\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{s^{*}}H_{i1q}^{R}(s|s^{*},\alpha,\theta^{R})[H_{i2q}^{R}(s|s^{*},\alpha,\theta^{R})]^{\tt
T}A_{i}^{R}(ds,t^{*};q)\right|$ $\displaystyle\stackrel{{\scriptstyle
p}}{{\rightarrow}}0.$

These conditions imply that $\left|\frac{1}{n}\langle
U_{pl}^{R}\rangle-\frac{1}{n}I_{pl}^{R}\right|\stackrel{{\scriptstyle
p}}{{\rightarrow}}0$ as $n\rightarrow\infty$. These are the analogous results
to those in the usual asymptotic theory of maximum likelihood estimators,
where the Fisher information is equal to the expected value of the squared
partial derivative with respect to the parameter of the log-likelihood
function and also the negative of the expected value of the second partial
derivative of the log-likelihood function. They are usually satisfied by
imposing a set of conditions that allows for the interchange of the order of
integration with respect to $s$ and the partial differentiation with respect
to the parameters; see, for instance, [6], [2], and [24] in similar but
simpler settings.

The proof of Proposition 5.3 then relies on the asymptotic representation

$\sqrt{n}\left[\begin{array}[]{c}\hat{\alpha}(s^{*},t^{*})-\alpha\\\
\hat{\theta}^{R}(s^{*},t^{*})-\theta^{R}\end{array}\right]=\left[\frac{1}{n}I_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})\right]^{-1}\left[\frac{1}{\sqrt{n}}U_{pl}^{R}(\alpha,\theta^{R}|s^{*},t^{*})\right]+o_{p}(1).$

In fact, this representation is also crucial for finding the asymptotic
properties of the estimators of the $\Lambda_{0q}(\cdot)$s, which we will now
present. By first-order Taylor expansion and under regularity conditions, we
have the representations, for each $q=1,\ldots,Q$, given by

$\displaystyle\hat{\Lambda}_{0q}(s^{*},t)$ $\displaystyle=$
$\displaystyle\int_{0}^{t}\frac{I\\{S_{q}^{0R}(s^{*},w|\hat{\alpha}(s^{*},t^{*}),\hat{\theta}^{R}(s^{*},t^{*}))>0\\}}{S_{q}^{0R}(s^{*},w|\hat{\alpha}(s^{*},t^{*}),\hat{\theta}^{R}(s^{*},t^{*})}\sum_{i=1}^{n}N_{i}^{R}(s^{*},dw;q)$
$\displaystyle=$
$\displaystyle\int_{0}^{t}\frac{I\\{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})>0\\}}{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})}\sum_{i=1}^{n}N_{i}^{R}(s^{*},dw;q)+$
$\displaystyle\left[B_{1q}(s^{*},t|\alpha,\theta^{R}),\
B_{2q}(s^{*},t|\alpha,\theta^{R})\right]\left[\begin{array}[]{c}\hat{\alpha}(s^{*},t^{*})-\alpha\\\
\hat{\theta}^{R}(s^{*},t^{*})-\theta^{R}\end{array}\right]+o_{p}(1/\sqrt{n})$

where

$\displaystyle
B_{1q}(s^{*},t|\alpha,\theta^{R})=-\int_{0}^{t}I\\{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})>0\\}\times$
$\displaystyle\left\\{\frac{\nabla_{\alpha}S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})}{[S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})]^{2}}\right\\}\sum_{i=1}^{n}N_{i}^{R}(s^{*},dw;q);$
$\displaystyle
B_{2q}(s^{*},t|\alpha,\theta^{R})=-\int_{0}^{t}I\\{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})>0\\}\times$
$\displaystyle\left\\{\frac{\nabla_{\theta}^{R}S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})}{[S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})]^{2}}\right\\}\sum_{i=1}^{n}N_{i}^{R}(s^{*},dw;q).$

For $q=1,\ldots,Q,$ let

$\Lambda_{0q}^{*}(s^{*},t)=\int_{0}^{t}I\\{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})>0\\}\lambda_{0q}(w)dw.$

Observe that

$\int_{0}^{t}\frac{I\\{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})>0\\}}{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})}\sum_{i=1}^{n}A_{i}^{R}(s^{*},dw|\alpha,\theta^{R})=\Lambda_{0q}^{*}(s^{*},t),$

implying that

$\displaystyle\int_{0}^{t}\frac{I\\{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})>0\\}}{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})}\sum_{i=1}^{n}N_{i}^{R}(s^{*},dw;q)=$
$\displaystyle\int_{0}^{t}\frac{I\\{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})>0\\}}{S_{q}^{0R}(s^{*},w|\alpha,\theta^{R})}\sum_{i=1}^{n}M_{i}^{R}(s^{*},dw;q)+\Lambda_{0q}^{*}(s^{*},t).$

Let us also define

$\displaystyle\hat{Z}_{q}^{R}(s^{*},t)=[\hat{\Lambda}_{0q}(s^{*},t)-\Lambda_{0q}^{*}(s^{*},t)]-$
$\displaystyle\left[B_{1q}(s^{*},t|\alpha,\theta^{R}),\
B_{2q}(s^{*},t|\alpha,\theta^{R})\right]\left[\begin{array}[]{c}\hat{\alpha}(s^{*},t^{*})-\alpha\\\
\hat{\theta}^{R}(s^{*},t^{*})-\theta^{R}\end{array}\right];$
$\tilde{H}_{iq}^{R}(w|s^{*})=\sum_{j:S_{j}\leq
s^{*}}H_{iq}^{R}[\mathcal{E}_{iqj}^{-1}(w)|s^{*}]I\\{\mathcal{E}_{iq}(S_{j-1})<w\leq\mathcal{E}_{iq}(S_{j})\\};$

and form the vectors of functions $\hat{Z}=(\hat{Z}_{q},q=1,\ldots,Q)$ and
$\tilde{H}_{i}^{R}=(\tilde{H}_{iq}^{R},q=1,\ldots,Q)$. Then, the main
asymptotic representation leading to the asymptotic results for the RCR
component parameters is given by, for $t\in[0,t^{*}]$,

$\displaystyle\sqrt{n}\left[\begin{array}[]{c}\hat{\alpha}(s^{*},t^{*})-\alpha\\\
\hat{\theta}^{R}(s^{*},t^{*})-\theta^{R}\\\
\hat{Z}(s^{*},t)\end{array}\right]$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cc}\left(\frac{1}{n}I_{pl}^{R}(s^{*},t^{*})\right)^{-1}&0\\\
0&I\end{array}\right]\times$ (28)
$\displaystyle\frac{1}{\sqrt{n}}\int_{0}^{t}\left[\begin{array}[]{c}\tilde{H}_{i}^{R}(w|s^{*})\\\
\mbox{Dg}\left(\frac{I\\{S^{0R}(s^{*},w)/n>0\\}}{S^{0R}(s^{*},w)/n}\right)\end{array}\right]M_{i}^{R}(s^{*},dw)+o_{p}(1).$

Using this representation, we obtain the following asymptotic properties,
though not stated in the most general form.

######  Proposition 5.4.

Under certain regularity conditions, we have

  * (i)

$\sqrt{n}\left[\begin{array}[]{c}\hat{\alpha}(s^{*},t^{*})-\alpha\\\
\hat{\theta}^{R}(s^{*},t^{*})-\theta^{R}\end{array}\right]$ and
$\sqrt{n}\hat{Z}(s^{*},t)$ are asymptotically independent;

  * (ii)

For each $q=1,\ldots,Q$, $\hat{\Lambda}_{0q}(s^{*},t)$ converges uniformly in
probability to $\Lambda_{0q}(t)$ for $t\in[0,t^{*}]$;

  * (iii)

For each $q=1,\ldots,Q$, and for $t\in[0,t^{*}]$,
$\sqrt{n}[\hat{\Lambda}_{0q}(s^{*},t)-\Lambda_{0q}(t)]$ converges in
distribution to a normal distribution with mean zero and variance

$\Gamma_{q}(s^{*},t)=\int_{0}^{t}\frac{d\Lambda_{0q}(w)}{s_{q}^{0R}(s^{*},w)}+(b_{1q}(s^{*},t),\
b_{2q}(s^{*},t))[\mathfrak{I}_{pl}^{R}(s^{*},t^{*})]^{-1}(b_{1q}(s^{*},t),\
b_{2q}(s^{*},t))^{\tt T},$

where $(b_{1q}(s^{*},t),\ b_{2q}(s^{*},t))=\mbox{pr-
lim}\left\\{\frac{1}{n}(B_{1q}(s^{*},t),\ B_{2q}(s^{*},t))\right\\}$.

  * (iv)

More, generally, for $q=1,\ldots,Q$, the stochastic process
$\\{\sqrt{n}[\hat{\Lambda}_{0q}(s^{*},t)-\Lambda_{0q}(t)]:\ t\in[0,t^{*}]\\}$
converges weakly in Skorokhod’s space
$(\mathfrak{D}[0,t^{*}],\mathcal{D}[0,t^{*}])$ to a zero-mean Gaussian process
with covariance function

$\displaystyle\Gamma_{q}(s^{*},t_{1},t_{2})$ $\displaystyle=$
$\displaystyle\int_{0}^{\min(t_{1},t_{2})}\frac{d\Lambda_{0q}(w)}{s_{q}^{0R}(s^{*},w)}+$
$\displaystyle(b_{1q}(s^{*},t_{1}),\
b_{2q}(s^{*},t_{1}))[\mathfrak{I}_{pl}^{R}(s^{*},t^{*})]^{-1}(b_{1q}(s^{*},t_{2}),\
b_{2q}(s^{*},t_{2}))^{\tt T}.$

This covariance function is consistently estimated by replacing the unknown
functions by their empirical counterparts and replacing the finite-dimensional
parameters by their estimators.

We point out that the $\hat{\Lambda}_{0q},q=1,\ldots,Q,$ are asymptotically
dependent, and are also not independent of
$(\hat{\alpha}(s^{*},t^{*}),\hat{\theta}^{R}(s*,t^{*}))$. The results stated
above generalize those in [2] and [24] to a more complex and general
situation.

##  6 Illustration of Estimation Approach on Simulated Data

In this section we provide a numerical illustration of the estimation
procedure when given a sample data. The illustrative sample data set with
$n=50$ units is depicted in Figure 2. This is generated from the proposed
model with the following characteristics. For the $i$th sample unit, the
covariate values are generated according to $X_{i1}\ {\sim}\ \mbox{BER}(0.5)$,
$X_{i2}\ {\sim}\ N(0,1)$ with $X_{i1}$ and $X_{i2}$ independent, where
$\mbox{BER}(p)$ is the Bernoulli distribution with success probability of $p$.
The end of monitoring time is $\tau_{i}\ {\sim}\ \mbox{EXP}(5)$, where
$\mbox{EXP}(\lambda)$ is the exponential distribution with mean $\lambda$. For
the RCR component with $Q=3$, the baseline (crude) hazard rate function for
risk $q\in\\{1,2,3\\}$, is a two-parameter Weibull given by

$\lambda_{0q}(t|\kappa_{q}^{*},\theta_{q}^{*})=\frac{\kappa_{q}^{*}}{\theta_{q}^{*}}\left(\frac{t}{\theta_{q}^{*}}\right)^{\kappa_{q}^{*}-1},\quad
t\geq 0,$

with $\kappa_{q}^{*}\in\\{2,2,3\\}$ and $\theta_{q}^{*}\in\\{0.9,1.1,1\\}$.
The associated (crude) survivor function for risk $q$ is

$\bar{F}_{0q}(t|\kappa_{q}^{*},\theta_{q}^{*})=\frac{\kappa_{q}^{*}}{\theta_{q}^{*}}\left(\frac{t}{\theta_{q}^{*}}\right)^{\kappa_{q}^{*}-1}\exp\left\\{-\left(\frac{t}{\theta_{q}^{*}}\right)^{\kappa_{q}^{*}}\right\\},\quad
t\geq 0.$

For risk $q$, the effective age process function is
$\mathcal{E}_{iq}(s)=s-S_{iqN_{iq}^{R}(s-)}^{R}$, the backward recurrence time
for this risk. For the effects of the accumulating event occurrences,
$\rho_{q}(N_{i}^{R}(s-);\alpha_{q})=(\alpha_{q})^{\log(1+N^{R}_{iq}(s-))}$.
For the HS component, $\mathfrak{V}=\\{1,2,3\\}$ with state ‘$1$’ an absorbing
state, so $\gamma$ is a scalar. For the LM component,
$\mathfrak{W}=\\{1=\mbox{High},2=\mbox{Normal},3=\mbox{Low}\\}$, so $\kappa$
is a two-dimensional vector. The infinitesimal generator matrices $\eta$ for
the LM process and $\xi$ for the HS process are, respectively,

$\displaystyle\eta=\begin{bmatrix}-0.3&0.2&0.1\\\ 0.1&-0.2&0.1\\\
0.1&0.2&-0.3\end{bmatrix}\quad\mbox{and}\quad\xi=\begin{bmatrix}0&0&0\\\
0.2&-0.7&0.5\\\ 0.05&0.5&-0.55\end{bmatrix}.$

The values in the first row for the $\xi$-matrix are all zeros because state
$1$ in HS is absorbing. The true values of the remaining model parameters are
given in the second column of Table 1.

Parameter | True | Estimate | Est. Standard Error
---|---|---|---
$\alpha_{1}$ | 1.50 | 1.58 | 0.13
$\alpha_{2}$ | 1.20 | 1.16 | 0.15
$\alpha_{3}$ | 2.00 | 2.10 | 0.21
$\beta^{R}_{1}$ | 1.00 | 1.17 | 0.09
$\beta^{R}_{2}$ | -1.00 | -0.94 | 0.10
$\gamma^{R}_{1}$ | 1.00 | 1.17 | 0.09
$\kappa^{R}_{1}$ | 1.00 | 1.10 | 0.09
$\kappa^{R}_{2}$ | -1.00 | -0.68 | 0.10
$\beta^{W}_{1}$ | 1.00 | 1.02 | 0.28
$\beta^{W}_{2}$ | -1.00 | -0.78 | 0.28
$\gamma^{W}_{1}$ | 1.00 | 0.97 | 0.27
$\nu^{W}_{1}$ | 1.00 | 0.82 | 0.06
$\nu^{W}_{2}$ | 1.00 | 1.11 | 0.09
$\nu^{W}_{3}$ | -2.00 | -2.01 | 0.09
$\beta^{V}_{1}$ | 1.00 | 0.93 | 0.17
$\beta^{V}_{2}$ | -1.00 | -1.00 | 0.16
$\kappa^{V}_{1}$ | 1.00 | 1.00 | 0.20
$\kappa^{V}_{2}$ | -1.00 | -1.59 | 0.40
$\nu^{V}_{1}$ | 1.00 | 1.24 | 0.04
$\nu^{V}_{2}$ | 1.00 | 1.04 | 0.06
$\nu^{V}_{3}$ | -2.00 | -2.30 | 0.06
Table 1: The second column contains the true values of the parameters in the
first column of the model that generated the illustrative sample data in
Figure 2 and also used in the simulation study. The third column are the
estimates of these parameters arising from the illustrative sample data, while
the fourth column contains the information-based estimates of the standard
errors of the estimators. The RCR component model parameters are the $\alpha$s
and those with superscript $R$. The LM component parameters are those with
superscript of $W$. Finally, the HS component parameters are those with
superscript $V$.

For each replication, the realized data for the $i$th unit among the $n$ units
are generated in the following manner. At time $s=0$, we first randomly assign
an initial LM state (either 1, 2, or 3) and HS state (either 2 or 3) uniformly
among the allowable states and independently for the LM and HS processes. We
specify a fixed length of the intervals partitioning $[0,\infty)$, which in
the simulation runs was set to $ds=0.001$, so we have intervals
$I_{k}=[s_{k},s_{k+1})=[k(ds),(k+1)(ds)),k=0,1,2,\ldots$. A smaller value of
$ds$ will make the data generation coincide more closely to the model, but at
the same time will also lead to higher computational costs, especially in a
simulation study with many replications. The data generation proceeds
sequentially over $k=0,1,2,\ldots$. Suppose that we have reached interval
$I_{k}=[s_{k},s_{k+1})$ with $W(s_{k})=w_{1}$ and $V(s_{k})=v_{1}$. For
$q\in\mathfrak{I}_{Q}$, generate a realization $e_{q}^{R}$ of $E_{q}^{R}$
according to a $\mbox{BER}(p_{q}^{R})$ with $p_{q}^{R}$ given by (4). Also,
for $w_{2}\neq w_{1}$, generate a realization $e_{w_{2}}^{W}$ of
$E_{w_{2}}^{W}$ according to a $\mbox{BER}(p_{w_{2}}^{W})$ with
$p_{w_{2}}^{W}$ given by (5). Finally, for $v\neq v_{1}$, generate a
realization $e_{v}^{V}$ of $E_{v}^{V}$ according to a $\mbox{BER}(p_{v}^{V})$
with $p_{v}^{V}$ given by (6).

  * •

If all the realizations $e_{q}^{R},q\in\mathfrak{I}_{Q}$,
$e_{w_{2}}^{W},w_{2}\in\mathfrak{I}_{\mathfrak{W}},w_{2}\neq w_{1}$, and
$e_{v}^{V},v\in\mathfrak{I}_{\mathfrak{V}},v\neq v_{1}$ are zeros, which means
no events occurred in the interval $I_{k}$, then we proceed to the next
interval $I_{k+1}$, provided that $\tau\notin I_{k}$, otherwise we stop.

  * •

If exactly one of the realizations $e_{q}^{R},q\in\mathfrak{I}_{Q}$,
$e_{w_{2}}^{W},w_{2}\in\mathfrak{I}_{\mathfrak{W}},w_{2}\neq w_{1}$, and
$e_{v}^{V},v\in\mathfrak{I}_{\mathfrak{V}},v\neq v_{1}$ equals one, so that an
event occurred, then we update the values of $N^{R},N^{W},N^{V}$ and proceed
to the next interval $I_{k+1}$, unless $e_{v}^{V}=1$ for a
$v\in\mathfrak{V}_{0}$ (i.e., there is a transition to an absorbing state) or
$\tau\in I_{k}$, in which case we stop.

In our implementation, since $0<ds\approx 0$, the success probabilities in the
Bernoulli distributions above will all be close to zeros, hence the
probability of more than one of the $e_{q}^{R},q\in\mathfrak{I}_{Q}$,
$e_{w_{2}}^{W},w_{2}\in\mathfrak{I}_{\mathfrak{W}},w_{2}\neq w_{1}$, and
$e_{v}^{V},v\in\mathfrak{I}_{\mathfrak{V}},v\neq v_{1}$ taking values of one
is very small. But, in case there were at least two of them with values of
one, then we randomly choose one to take the value of one and the others are
set to zeros. Thus, we always have

$\sum_{q\in\mathfrak{I}_{Q}}e_{q}^{R}+\sum_{w_{2}\in\mathfrak{I}_{\mathfrak{W}};\
w_{2}\neq w_{1}}e_{w_{2}}^{W}+\sum_{v\in\mathfrak{I}_{\mathfrak{V}};\ v\neq
v_{1}}e_{v}^{V}\in\\{0,1\\},$

which means there is at most one event in any interval. Note that whether we
reach an absorbing state or not, there will always be right-censored event
times or sojourn times.

We remark that the event time generation could also have been implemented by
first generating a sojourn time and then deciding which event occurred
according to a multinomial distribution at the realized sojourn time as
indicated in subsection 3.2. In addition, depending on the form of the
baseline hazard rate functions $\lambda_{0q}(\cdot)$s (e.g., Weibulls) and the
effective age processes $\mathcal{E}_{iq}(\cdot)$s (e.g., backward recurrence
times), a more direct and efficient manner of generating the event times using
a variation of the Probability Integral Transformation is possible, without
having to do the partitioning of the time axis as done above. However, the
method above is more general, though approximate, and applicable even if the
$\lambda_{0q}$’s or the $\mathcal{E}_{iq}$’s are of more complex forms.

Graphical plots associated with the generated illustrative sample data are
provided in Figure 2 of Section 4.1. For this sample data set there were $36$
units that reached the absorbing state, with a mean time to absorption of
about $t=1$; while $14$ units did not reach the absorbing state before hitting
their respective end of their monitoring periods, with the mean monitoring
time at about $t=2$. Recall that $\tau_{i}$’s were distributed as
$\mbox{EXP}(5)$ so the mean of $\tau_{i}$ is 5. One may be curious why the
mean monitoring time for those that did not get absorbed is about 2 and not
close to 5. The reason for this is because of an induced selection bias. Those
units who got absorbed will tend to have longer monitoring times, hence those
that were not absorbed will tend to have shorter monitoring times, explaining
a reduction in the mean monitoring times for the subset of units that were not
absorbed.

We have developed programs in R [26] for implementing the semi-parametric
estimation procedure for the joint model described above. We used these
programs on this illustrative sample data set to obtain estimates and their
estimated standard errors of the model parameters. The third and fourth
columns of Table 1 contain the parameter estimates and estimates of their
standard errors, respectively, of the finite-dimensional parameters in the
first column, whose true values are given in the second column. Figure 3
depicts the true baseline survivor functions, which are Weibull survivor
functions, together with their semi-parametric estimates for each risk of the
three risks in the RCR component. The estimates of the baseline survivor
functions are obtained using the product-integral representation of cumulative
hazard functions. For this generated sample data, the estimates obtained and
the function plots demonstrate that there is reasonable agreement between the
true parameter values and true functions and their associated estimates,
indicating that the semi-parametric estimation procedure described earlier
appears to be viable.

Figure 3: The true (red, dashed) and estimated (blue, solid) baseline (crude)
survivor functions for each of the three risks based on the simulated
illustrative sample data set with $n=50$ units depicted in Figure 2.

##  7 Finite-Sample Properties via Simulation Studies

###  7.1 Simulation Design

We have provided asymptotic results of the estimators in Section 5. In this
section we present the results of simulation studies to assess the finite-
sample properties of the estimators of model parameters. This will provide
some evidence whether the semi-parametric estimation procedure, which appears
to perform satisfactorily for the single illustrative sample data set in
Section 6, performs satisfactorily over many sample data sets. These
simulation studies were implemented using R programs we developed, in
particular, the programs utilized in estimating parameters in the preceding
section. In these simulation studies, as in the preceding section, when we
analyze each of the sample data, the baseline hazard rate functions are
estimated semi-parametrically, even though in the generation of each of the
sample data sets, two-parameter Weibull models were used in the RCR
components. Aside from the set of model parameters described in Section 6, the
simulation study have the additional inputs which are the sample size $n$ and
the number of simulation replications Mreps, the latter set to 1000. The
sample sizes used in the two simulation studies are $n\in\\{50,100\\}$. For
fixed $n$, for each of the Mreps replications, the sample data generation is
as described in Section 6.

Table 2 presents some summary results pertaining to the three processes based
on the Mreps replications. The first three rows indicate the means and
standard deviations of the number of event occurrences per unit for each risk,
and the mean time for the first event occurrence of each risk. For example,
risk $1$ occurs about $2.6$ times per unit with a standard deviation $3.57$.
The mean time for risk $1$ to occur for the first time is about $0.48$. We
notice that occurrence frequencies for three risks are ordered according to
$\mbox{Risk 1}\succ\mbox{Risk 2}\succ\mbox{Risk 3}$, and consequently risk $1$
tends to have the shortest mean time to the first event. Also note that the
mean number of event occurrences per unit for each risk is around $2$, which
implies that there are not too few RCR events or too many RCR events (see the
property of “explosion” as discussed in [13]) per unit. This indicates that
the choice of the effective age function $\mathcal{E}_{iq}(\cdot)$ and the
accumulating event function $\rho_{q}(\cdot)$, together with the parameter
values we chose for the data generation, were reasonable.

The fourth to ninth rows show the mean and standard deviation of the number of
transitions to specific states per unit, the mean and standard deviation of
occupation times per unit for specific states, the mean and standard deviation
of sojourn times for specific states. For example, column 4 tells us that (i)
the mean number of transitions to state $2$ of the HS process (HS $2$ for
short) per unit is $2.34$; (ii) a unit would stay in HS $2$ for an approximate
time of $0.8$ on average; (iii) the mean sojourn time for HS $2$ is about
$0.34$. We do not include information for HS $1$ since it is an absorbing
state. Comparing the $V=2$ and $V=3$ columns, we find that units tend to
transit to HS state $2$ more often than to HS state $3$. The mean occupation
time for state $2$ per unit is longer compared to state 3. For the last three
columns, there are more transitions to state $2$ than to other states. Thus, a
unit tends to stay in LM state $1$ more often than in the other two LM states.

| RCR1 | RCR2 | RCR3 | V=2 | V=3 | W=1 | W=2 | W=3
---|---|---|---|---|---|---|---|---
Mean: Count | 2.60 | 2.06 | 1.62 |  |  |  |  |
SD: Count | 3.57 | 2.62 | 3.37 |  |  |  |  |
Mean: TimePerEvent | 0.48 | 0.61 | 0.78 |  |  |  |  |
Mean: NumTransition |  |  |  | 2.34 | 2.10 | 1.16 | 1.36 | 1.10
SD: NumTransition |  |  |  | 2.17 | 2.02 | 1.31 | 1.24 | 1.33
Mean: OccupationTime |  |  |  | 0.80 | 0.46 | 0.59 | 0.31 | 0.36
SD: OccupationTime |  |  |  | 1.83 | 0.78 | 1.51 | 0.58 | 0.71
Mean: SojournTime |  |  |  | 0.34 | 0.22 | 0.51 | 0.23 | 0.32
SD: SojournTime |  |  |  | 1.09 | 0.54 | 1.39 | 0.49 | 0.64
Table 2: Summary statistics for three processes for Mreps replications of data
set. The first three rows are for RCR events. The first three rows indicate
the mean/standard deviation of the number of recurrent event occurrences per
unit for each risk, the mean time for one recurrent event for each risk. The
fourth to ninth rows are for HS and LM events. They indicate the mean/standard
deviation of the number of transitions to the specific state per unit, the
mean/standard deviation of the occupation time for the specific state per
unit, the mean/standard deviation of the sojourn time for the specific state.
State $1$ of the HS process $V$ is absorbing, hence not included in the table.

Also, to obtain some insights into the model-induced dependencies among the
components, we also obtained the correlations among RCR, LM, and HS processes
over time from the simulated data. We first constructed a vector of six
variables over a finite number of time points $\mathcal{S}\subset[0,\tau]$
given by

$\\{Z(s)\equiv[I(V_{i}(s)=2),I(W_{i}(s)=2),I(W_{i}(s)=3),N_{i1}^{R}(s),N_{i2}^{R}(s),N_{i3}^{R}(s)]^{T}:s\in\mathcal{S}\\}.$

For each $s\in\mathcal{S}$, we then obtained the sample correlation matrix
$C(s)$ from $\\{Z_{i}(s),i=1,2,\ldots,n\\}$. Each of the Mreps replications
then yielded a $C(s)$, so we took the mean, element-wise, of these Mreps
correlation matrices. The matrix of scatterplots in Figure 4 provides the
plots of these mean correlation coefficients over time points
$s\in\mathcal{S}$. The point we are making here is that the joint model does
induce non-trivial patterns of dependencies over time among the three model
components.

Figure 4: Plots of sample (Pearson) correlations over a finite set of time
points in $[0,3]$. The random vector, at each time $s$ and for each sample
unit, is
$C(s)=[I(V_{i}(s)=2),I(W_{i}(s)=2),I(W_{i}(s)=3),N_{i1}^{R}(s),N_{i2}^{R}(s),N_{i3}^{R}(s)]$.
The sample correlation matrix is computed based on the $n=50$ sample units.
The element-wise means of the Mreps correlation matrices were then computed.
The plots depict these mean sample correlations for the pairs of variables in
$C(s)$ over time $s$.

###  7.2 Finite-Sample Properties of Estimators

The set of estimates obtained for one sample data set in the last section is
insufficient to assess the performance of the semi-parametric estimation
procedure. To get a sense of its performance we performed simulation studies
with $\mbox{\tt Mreps}=1000$ replications and sample sizes $n\in\\{50,100\\}$.
For each replication, we generated a sample data set according to the same
joint model, then obtained the set of estimates, via the semi-parametric
procedure, for this data set. Summary statistics, such as the means, standard
deviations, asymptotic standard errors, percentiles, and boxplots for all
Mreps estimates were then obtained or constructed. Table 3 shows these summary
statistics of the Mreps estimates for each parameter. The asymptotic standard
errors reported in Table 3 are the means of the asymptotic standard errors
over the Mreps replications. Also included are the percentile $95\%$
confidence intervals for each of the unknown parameters based on the Mreps
replications stratified according to $n=\\{50,100\\}$. Since the
$\Lambda_{0q}$s are functional nonparametric parameters, we only provide the
estimates at four selected time points. Due to space limitation, we also only
provide the results for a small subset of the $\eta$ and $\xi$ parameters in
Table 3. From these simulation results, we observe that the estimates are
close to the true values of the model parameters, and the sample standard
deviations are close to the asymptotic standard errors, providing some
validation to the semi-parametric estimation procedure for the joint model and
empirically lending support to the asymptotic results. By comparing the
results for $n=50$ and $n=100$ in Table 3, we find that when the sample size
$n$ increases, the performance of the estimators of the parameters improves
with biases and standard errors decreasing.

Sample Size | $n=50$ | $n=100$
---|---|---
Parameter | True | Mean | SD | ASE | PL | PU | Mean | SD | ASE | PL | PU
$\Lambda_{01}(0.3)$ | 0.11 | 0.09 | 0.03 | 0.04 | 0.05 | 0.16 | 0.1 | 0.02 | 0.03 | 0.07 | 0.15
$\Lambda_{01}(0.6)$ | 0.44 | 0.43 | 0.09 | 0.07 | 0.22 | 0.66 | 0.42 | 0.06 | 0.05 | 0.3 | 0.6
$\Lambda_{01}(0.9)$ | 1 | 0.97 | 0.19 | 0.18 | 0.57 | 1.61 | 0.99 | 0.12 | 0.12 | 0.71 | 1.46
$\Lambda_{01}(1.2)$ | 1.78 | 1.83 | 0.43 | 0.42 | 1.02 | 2.7 | 1.79 | 0.27 | 0.27 | 1.25 | 2.65
$\Lambda_{02}(0.3)$ | 0.09 | 0.07 | 0.03 | 0.04 | 0.04 | 0.13 | 0.08 | 0.02 | 0.03 | 0.05 | 0.12
$\Lambda_{02}(0.6)$ | 0.36 | 0.32 | 0.08 | 0.07 | 0.18 | 0.59 | 0.34 | 0.06 | 0.05 | 0.23 | 0.48
$\Lambda_{02}(0.9)$ | 0.81 | 0.78 | 0.16 | 0.15 | 0.41 | 1.36 | 0.79 | 0.12 | 0.1 | 0.54 | 1.14
$\Lambda_{02}(1.2)$ | 1.44 | 1.53 | 0.33 | 0.32 | 0.81 | 2.33 | 1.46 | 0.23 | 0.23 | 1.03 | 2.24
$\alpha_{1}$ | 1.5 | 1.52 | 0.15 | 0.13 | 1.13 | 2.06 | 1.5 | 0.1 | 0.09 | 1.25 | 1.83
$\alpha_{2}$ | 1.2 | 1.21 | 0.15 | 0.14 | 0.8 | 1.76 | 1.2 | 0.11 | 0.09 | 0.95 | 1.58
$\alpha_{3}$ | 2 | 2.02 | 0.23 | 0.22 | 1.43 | 2.7 | 2.01 | 0.14 | 0.14 | 1.56 | 2.6
$\beta^{R}_{1}$ | 1 | 1.03 | 0.11 | 0.1 | 0.78 | 1.45 | 1.02 | 0.08 | 0.07 | 0.82 | 1.25
$\beta^{R}_{2}$ | -1 | -1.03 | 0.11 | 0.11 | -1.29 | -0.85 | -1.02 | 0.07 | 0.08 | -1.17 | -0.89
$\gamma^{R}_{1}$ | 1 | 1.03 | 0.1 | 0.08 | 0.8 | 1.47 | 1.02 | 0.06 | 0.06 | 0.86 | 1.25
$\kappa^{R}_{1}$ | 1 | 1.03 | 0.1 | 0.09 | 0.77 | 1.52 | 1.02 | 0.08 | 0.06 | 0.83 | 1.29
$\kappa^{R}_{2}$ | -1 | -1.02 | 0.15 | 0.14 | -1.51 | -0.52 | -1.01 | 0.1 | 0.1 | -1.33 | -0.74
$\eta(2,1)$ | 0.1 | 0.1 | 0.04 | 0.03 | 0.04 | 0.18 | 0.1 | 0.02 | 0.02 | 0.05 | 0.16
$\eta(3,1)$ | 0.1 | 0.1 | 0.04 | 0.03 | 0.04 | 0.19 | 0.1 | 0.02 | 0.02 | 0.06 | 0.17
$\eta(1,2)$ | 0.2 | 0.2 | 0.05 | 0.06 | 0.11 | 0.29 | 0.2 | 0.04 | 0.04 | 0.12 | 0.29
$\eta(3,2)$ | 0.2 | 0.21 | 0.05 | 0.05 | 0.11 | 0.29 | 0.2 | 0.03 | 0.03 | 0.13 | 0.29
$\beta^{W}_{1}$ | 1 | 0.98 | 0.23 | 0.21 | 0.51 | 1.45 | 0.99 | 0.16 | 0.15 | 0.68 | 1.32
$\beta^{W}_{2}$ | -1 | -0.99 | 0.18 | 0.2 | -1.24 | -0.75 | -0.99 | 0.14 | 0.14 | -1.16 | -0.81
$\gamma^{W}_{1}$ | 1 | 0.99 | 0.23 | 0.23 | 0.56 | 1.48 | 0.99 | 0.15 | 0.15 | 0.71 | 1.33
$\nu^{W}_{1}$ | 1 | 0.99 | 0.07 | 0.06 | 0.74 | 1.25 | 0.99 | 0.05 | 0.04 | 0.82 | 1.15
$\nu^{W}_{2}$ | 1 | 0.99 | 0.09 | 0.07 | 0.71 | 1.29 | 0.99 | 0.06 | 0.06 | 0.8 | 1.17