# A Class of Semiparametric Yang and Prentice Frailty Models

Cássius Henrique Xavier Oliveira
Department of Statistics
Federal University of Minas Gerais
Belo Horizonte, MG, Brazil
<EMAIL_ADDRESS>
&Fábio Nogueira Demarqui
Department of Statistics
Federal University of Minas Gerais
Belo Horizonte, MG, Brazil
<EMAIL_ADDRESS>
Vinícius Diniz Mayrink
Department of Statistics
Federal University of Minas Gerais
Belo Horizonte, MG, Brazil
<EMAIL_ADDRESS>

###### Abstract

The Yang and Prentice (YP) regression models have garnered interest from the
scientific community due to their ability to analyze data whose survival
curves exhibit intersection. These models include proportional hazards (PH)
and proportional odds (PO) models as specific cases. However, they encounter
limitations when dealing with multivariate survival data due to potential
dependencies between the times-to-event. A solution is introducing a frailty
term into the hazard functions, making it possible for the times-to-event to
be considered independent, given the frailty term. In this study, we propose a
new class of YP models that incorporate frailty. We use the exponential
distribution, the piecewise exponential distribution (PE), and Bernstein
polynomials (BP) as baseline functions. Our approach adopts a Bayesian
methodology. The proposed models are evaluated through a simulation study,
which shows that the YP frailty models with BP and PE baselines perform
similarly to the generator parametric model of the data. We apply the models
in two real data sets.

_K_ eywords Survival $\cdot$ Frailty $\cdot$ Piecewise exponential
distribution $\cdot$ Bernstein polynomials

##  1 Introduction

In Survival analysis, a common goal is to evaluate potential risk factors on
the occurrence of events. Proportional hazards (PH) regression models, such as
Cox, (1972), and proportional odds (PO) model Bennett, (1983) are possible
approaches. Another alternative is the Yang and Prentice (YP) regression model
Yang and Prentice, (2005), which includes the PH and PO models as particular
cases. In YP models, the survival functions are allowed to intersect and this
provides an advantage over the PH and PO models Demarqui and Mayrink, (2021).
Several works address the YP model in the literature. Diao et al., (2013)
studied the extension of the YP model to accommodate potentially time-
dependent covariates. Wang, (2013) developed methods to calculate the sample
size based on YP models. Demarqui et al., (2019) introduced an approach to fit
the YP model using Bernstein polynomials (BP) for handling the baseline hazard
and odds, applicable within both frequentist and Bayesian frameworks. Demarqui
and Mayrink, (2021) proposed a semiparametric model for survival data using
the YP regression model and the piecewise exponential (PE) distribution as the
baseline hazard function. The fit of this model can be done using R package
YPPE (Demarqui,, 2020). Miranda Filho, (2022) proposed a class of multivariate
survival models based on Archimedean copulas with margins modeled by the YP
model.

Due to the flexibility of BP in approximating continuous functions, some works
in the literature bring applications of these polynomials in survival
analysis. Panaro, (2020) developed an R package named spsurv Panaro et al.,
(2020) model survival times using BP coupled to some regression structures as
PH and PO models. Demarqui et al., (2019) used the BP to model the baseline
functions of the YP model.

When it comes to the construction of survival models, in various research, it
is commonly assumed that the times until the event are mutually independent.
However, this assumption may not hold in certain scenarios. Consider, for
instance, a patient who is infected repeatedly by a virus. We say that the
individual experienced recurrent events of infections. Assuming independence
in these cases may be inappropriate (Klein and Moeschberger,, 2006). Immunity
acquired due to an infection can change the likelihood of future infections.
Therefore, the recurrent events are not mutually independent, as one event can
alter susceptibility to subsequent events.

When survival times exhibit correlations among them, the data are classified
as multivariate. This scenario arises, for example, in cases of recurrent
events experienced by individuals. Conversely, when the assumption of
independence between time-to-events holds true, the data are considered
univariate. Klein and Moeschberger, (2006), Hanagal, (2011) and others argue
that a commonly used approach to deal with some dependence on survival data is
to assume that the time-to-event is conditionally independent on a set of
unobserved variables, called frailties. The concept of frailty, introduced by
Vaupel et al., (1979), defines it as a latent and multiplicative random
variable. The authors used frailties, also called random effects, to explain
the effect of unobserved heterogeneity on the mortality of a population.
Clayton and Cuzick, (1985) used frailties to explain the heterogeneity about
the hazard function in an extension of the PH model for multivariate survival
data. Frailty models can also be used to accommodate the association between
recurrent events, as in (Lawless,, 1987). Huang and Wang, (2004) considered a
joint modeling of recurrent events and a terminal event, utilizing frailties
to model the correlation between the intensity of the recurrent event process
and the hazard of the failure time. Liu et al., (2004) considered frailty PH
models for the recurrent and terminal event processes in which shared frailty
is included in both hazard functions. Mazroui et al., (2012) proposed a joint
frailty model to analyze recurrences and death, using two gamma-distributed
frailties to handle both the recurrences dependence and the dependence between
the recurrences and the death times. Schneider et al., (2020) used the frailty
model to fit survival data subjected to dependent censoring. Some works use
frailty in PO models such as Economou and Caroni, (2007), Lin and Wang,
(2011), and Gupta and Peng, (2014).

The primary objective of this study is to develop a class of YP frailty
models, within the Bayesian framework, with three baseline functions -
exponential, PE, and BP. It allows one to analyze survival data arranged in
two configurations.

  * •

Univariate survival datasets. The frailty term serves to explain unobserved
heterogeneities, that is, variations in survival time that are not explained
by the fixed effects of the models. In this case, we refer to the element of
the frailty as individual frailty;

  * •

Multivariate survival datasets in which individuals present recurrent events.
Thus, frailty is used to accommodate the association between the survival
times of the same individual. In this context, we can understand the
individual as a cluster, and frailty here is referred to as shared frailty.

This work is organized as follows. Section 2 presents the fundamental concepts
of the YP model. Furthermore, it discusses the BP and the PE model, which will
be used to model the baseline hazard functions, and presents the proposed
models. Section 3 analyses the results of the Monte Carlo study. Section 4
illustrates two real applications of our models. We close this text with
discussions of some results and perspectives for future research in Section 5.

##  2 Model formulation

This section presents the YP frailty models whose baseline functions are
modeled using exponential distribution, BP, and PE distribution. Let’s start
by discussing some theoretical aspects of the YP model. Yang and Prentice,
(2005) proposed a model in which survival curves can intersect. Let $T>0$ be
the random variable denoting the time-to-event. The YP model can be
characterized in terms of the survival function

$\displaystyle S(t|\mathbf{x})=\left[1+\frac{\nu}{\xi}R_{0}(t)\right]^{-\xi},$

where $\nu=\exp(\mathbf{x}\boldsymbol{\psi})$ and
$\xi=\exp(\mathbf{x}\boldsymbol{\phi})$,
$\boldsymbol{\psi}=(\psi_{1},\ldots,\psi_{p})^{\prime}$ and
$\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{p})^{\prime}$ are vectors of
regression parameters without intercepts and $\mathbf{x}=(x_{1},\ldots,x_{p})$
is a row vector of explanatory variables. The function $R_{0}(t)$ is the
baseline odds function and is defined as $R_{0}(t)=\frac{F_{0}(t)}{S_{0}(t)}$,
where $S_{0}(t)=1-F_{0}(t)$ is the baseline survival function, and $F_{0}(t)$
is baseline cumulative distribution function. We can express the hazard
function of this model as

$\displaystyle h(t|\mathbf{x})=\frac{\nu\xi}{\nu F_{0}(t)+\xi
S_{0}(t)}h_{0}(t),$

where $h_{0}(t)=-\frac{d}{dt}\log[S_{0}(t)]$ is baseline hazard function. The
parameter $\nu$ is interpreted as the short-term hazards ratio, because

$\lim_{t\to 0}\frac{h(t|\mathbf{x})}{h(t|\mathbf{0})}=\nu,$

and $\boldsymbol{\psi}$ as the short-term regression coefficients vector. We
can interpret $\xi$ as the long-term hazards ratio, since

$\lim_{t\to\infty}\frac{h(t|\mathbf{x})}{h(t|\mathbf{0})}=\xi,$

and $\boldsymbol{\phi}$ is the long-term regression coefficients vector.

The YP model can be reduced to the PH and PO models. Note that, when
$\boldsymbol{\psi}=\boldsymbol{\phi}$,
$h(t|\mathbf{x})=h_{0}(t)\exp(\mathbf{x}\boldsymbol{\phi})$, and this is the
hazard function of the PH model. On the other hand, if
$\boldsymbol{\phi}=\textbf{0}$, we have

$S(t|\mathbf{x})=\left[1+R_{0}(t)\exp(\mathbf{x}\boldsymbol{\psi})\right]^{-1}\Rightarrow
R(t|\mathbf{x})=\frac{F(t|\mathbf{x})}{S(t|\mathbf{x})},$

and this is the expression of the odds function in the PO model. In YP model,
when $\psi_{j}\phi_{j}<0$, for any pair of coefficients $(\psi_{j},\phi_{j})$,
with $j\in\\{1,\ldots,p\\}$, the survival curves intersect.

We establish the YP frailty model by multiplying expressions of $\nu$ and
$\xi$ by a random effect $z=\exp(w)$. Klein and Moeschberger, (2006) highlight
that $w$ is usually assumed to have a distribution with mean zero and unknown
variance $\sigma^{2}_{w}$. Consider individuals $k$ and $k^{\prime}$ of which
$\mathbf{x}_{k}=(x_{k,1},\ldots,x_{k,p})$ and
$\mathbf{x}_{k}^{\prime}=(x_{k^{\prime},1},\ldots,x_{k^{\prime},p})$ are row
vectors of covariates, respectively. The ratio between the hazard functions
depends not only on the observed characteristics but also on the random
effects of the two individuals, since

$\displaystyle\lim_{t\to
0}\frac{h(t|\mathbf{x}_{k},z_{k})}{h(t|\mathbf{x}_{k^{\prime}},z_{k^{\prime}})}=$
$\displaystyle\frac{z_{k}}{z_{k^{\prime}}}\exp\left[\psi\left(\mathbf{x}_{k}-\mathbf{x}_{k^{\prime}}\right)\right],\text{
and }$
$\displaystyle\lim_{t\to\infty}\frac{h(t|\mathbf{x}_{k},z_{k})}{h(t|\mathbf{x}_{k^{\prime}},z_{k^{\prime}})}=$
$\displaystyle\frac{z_{k}}{z_{k^{\prime}}}\exp\left[\phi\left(\mathbf{x}_{k}-\mathbf{x}_{k^{\prime}}\right)\right].$

We model the baseline function using BP. The BP is a linear combination of
bases introduced by Bernstein, (1912). Consider the continuous function
$C(\cdot)$ defined in a range $(0,\tau]$. It can be approximated arbitrarily
by BP Lorentz, (1986). The BP of degree $m$ evaluated in $t\in(0,\tau]$ with
base $B_{m}=(B_{1,m},B_{2,m},\ldots,B_{m,m})$ and coefficients
$b_{m}=(b_{1,m},b_{2,m},\ldots,b_{m,m})$ to approximate the function
$C(\cdot)$ is

$B^{C}_{m}(t)=\sum_{k=1}^{m}b_{k,m}B_{k,m}(t),$

where $b_{k,m}=C\left(\frac{k\tau}{m}\right)$ and $B_{k,m}={m\choose
k}\left(\frac{t}{\tau}\right)^{k}\left(1-\frac{t}{\tau}\right)^{m-k}$;
$k=1,2,\ldots,m$. The derivative of BP with respect to $t$ can be written as

$b_{m}^{C}(t)=\sum_{k=1}^{m}\left\\{C\left(\frac{k}{m}\tau\right)-C\left(\frac{k-1}{m}\tau\right)\right\\}\frac{f_{\beta}\left(\frac{t}{\tau},m-k+1\right)}{\tau},$

where $f_{\beta}\left(\frac{t}{\tau},k,m-k+1\right)$ is the density of a Beta
distribution with parameters $k$ and $m-k+1$ valued at $\frac{t}{\tau}$.
Lorentz, (2012) shows that $B^{C}_{m}(t)\to_{m\to\infty}C(t)$ uniformly, and
$b^{C}_{m}(t)\to_{m\to\infty}\frac{d}{dt}C(t)$ uniformly on $(0;\tau]$.

In Survival analysis, Osman and Ghosh, (2012) used the expression of
$b^{C}_{m}(t)$ to model the hazard function $h(t|\boldsymbol{\gamma})$ and
$B^{C}_{m}(t)$ to handle the accumulated hazard function
$H(t|\boldsymbol{\gamma})$. Following the reasoning of these authors, assume
$\boldsymbol{\gamma}=(\gamma_{1},...,\gamma_{m})$, with
$\gamma_{k}=C^{*}\left(\frac{k}{m}\tau\right)-C^{*}\left(\frac{k-1}{m}\tau\right);\gamma_{k}\geq
0,k=1,...,m$. Note that $\boldsymbol{\gamma}$ is not time-dependent and its
values are unknown. Also consider
$\mathbf{g}_{m}(t)=\left(g_{1,m}(t),...,g_{m,m}(t)\right)^{\prime}$, where
$g_{k,m}(t)=\frac{1}{\tau}f_{\beta}\left(\frac{t}{\tau},k,m-k+1\right);g_{k,m}(t)\geq
0,k=1,...,m$. Thus,
$b_{m}^{C^{*}}(t)=\boldsymbol{\gamma}~{}\textbf{g}_{m}(t)$. Osman and Ghosh,
(2012) used this expression to model the hazard function. That is,
$h(t|\boldsymbol{\gamma})=\boldsymbol{\gamma}~{}\mathbf{g}_{m}(t);t\in[0,\infty)$.

BP has some advantages because offers flexibility in modeling different shapes
of hazard functions, facilitating the model’s adaptation to the specific
characteristics of the data. Furthermore, BP has good derivation properties,
and the log-likelihood function has a friendly form. The monotonicity of the
accumulated hazard function is naturally modeled by BP, since that
$\gamma_{k}\geq 0,\forall~{}k\in\\{1,...,m\\}$. This function is expressed by
$H(t|\boldsymbol{\gamma})=\int_{0}^{t}h(u|\boldsymbol{\gamma})du=\boldsymbol{\gamma}~{}\mathbf{G}_{m}(t),$
with $\mathbf{G}_{m}(t)=\left(G_{1,m}(t),...,G_{m,m}(t)\right)^{\prime},$
where
$G_{m,k}(t)=\int_{0}^{t}f_{\beta}\left(\frac{u}{\tau};k,m-k+1\right)d\left(\frac{u}{\tau}\right),\forall~{}k\in\\{1,...,m\\}.$
The function $G_{m,k}(t)$ is the Beta cumulative distribution function with
parameters $k$ and $m-k+1$.

These authors also discuss some aspects of choosing $\tau$. It is necessary
that $\tau<\infty$, such that $\tau=\inf\\{t:S(t)=0\\}$. In practice, in
survival analysis, $\tau$ is chosen as the maximum value among the times
observed until the occurrence of the event of interest or until the follow-up
stops. Here, we will denote it by $\hat{\tau}$. But, using this choice, it is
not possible to satisfy $H(\tau|\boldsymbol{\gamma})=\infty$. Besides, there
is no information about survival times in the region $t>\hat{\tau}$ Demarqui
et al., (2019). Therefore, this choice requires an adjustment in the hazard
and cumulative hazard functions. As a solution, Osman and Ghosh, (2012)
suggest some alterations in these functions, as follows:

$\displaystyle h^{*}(t|\boldsymbol{\gamma})$
$\displaystyle=\begin{cases}h(t|\boldsymbol{\gamma}),\mbox{ if }0\leq
t<\hat{\tau},\\\ m~{}\frac{\boldsymbol{\gamma}}{\hat{\tau}},\mbox{ if
}t\geq\hat{\tau},\end{cases}$ $\displaystyle H^{*}(t|\boldsymbol{\gamma})$
$\displaystyle=\begin{cases}H(t|\boldsymbol{\gamma}),\mbox{ if }0\leq
t<\hat{\tau},\\\
H(t|\boldsymbol{\gamma})+m(t-\hat{\tau})~{}\frac{\boldsymbol{\gamma}}{\hat{\tau}},\mbox{
if }t\geq\hat{\tau}.\end{cases}$

Another way also found in the literature to model baseline function is the use
of the PE distribution. This model was introduced by Kalbfleisch and Prentice,
(1973). Consider a time grid $\rho=\\{\rho_{0},\ldots,\rho_{m}\\}$. Note that
$\rho$ makes a partition of the time axis in $m$ intervals at the points
$\rho_{0},\rho_{1},\ldots,\rho_{m}$, with
$0=\rho_{0}<\rho_{1}<\ldots<\rho_{m}<\infty.$ The intervals generated from
that partition are
$I_{1}=(\rho_{0},\rho_{1}],I_{2}=(\rho_{1},\rho_{2}],\ldots,I_{m}=(\rho_{m-1},\rho_{m}]$.
The set $\rho$ (or $m$, alternatively) can be established in different ways.
Breslow, (1974) and Demarqui and Mayrink, (2021) assume that $\rho$ is a known
set composed by each of different time-to-event observations. Kalbfleisch and
Prentice, (1973) declare that the choice of $\rho$ can be independent of the
data set. On the other hand, Demarqui and Mayrink, (2021) argues that large
$m$ values can provide unstable estimates. In other approaches, $\rho$ is
treated as being random; see Demarqui et al., (2011, 2012).

In truth, the choice of $\rho$ influences the inferential results since we
assume that the hazard function in each interval is constant and given by
$h(t|\gamma_{j})=\gamma_{j},\mbox{ for }t\in I_{j},j=1,\ldots,m\mbox{ and
}\gamma_{j}>0.$ The cumulative hazard function is given by
$H(t|\boldsymbol{\gamma})=\sum_{j=1}^{m}\gamma_{j}(t_{j}-\rho_{j-1}),$ where

$\displaystyle t_{j}=\begin{cases}\rho_{j-1}\mbox{, if }t<\rho_{j-1};\\\
t\mbox{, if }\rho_{j-1}<t\leq\rho_{j};\\\ \rho_{j}\mbox{, if
}t>\rho_{j},\end{cases}$

for $j=1,\ldots,m.$ From the calculation of the cumulative hazard function, we
can use the expression
$S(t|\boldsymbol{\gamma})=\exp\\{-H(t|\boldsymbol{\gamma})\\}$ to find the
survival function.

One of the main advantages of PE is its flexibility in adjusting to different
shapes of the hazard function over time. By dividing the study period into
intervals and assuming a constant hazard function within each interval, but
allowing these rates to vary between the intervals, the model can adapt to a
variety of hazard functions that would not be well captured by an exponential
model. Furthermore, the estimated parameters for each time interval can be
directly interpreted as hazard functions, making the model intuitive for
researchers.

Once we have discussed the necessary elements to construct the class of models
that we propose in this work, we now proceed with the establishment of the
likelihood function. Let $L$ be the number of individuals. Denote by $C_{i}$
the time to the administrative censoring, that is, the time until loss of
follow-up for some reason external to the study, with $i=1,\ldots,L$. Denote
by $R_{i,j}$ the gap-time between the $(j-1)$-th and $j$-th occurrences of the
recurrent event, and let $T_{i,j}=\sum_{j^{\prime}=1}^{j}R_{i,j^{\prime}}$ be
the total observation time until the $j$-th recurrent event. Suppose the
$i$-th subject experiences a total of $n_{i}$ recurrent events. When
$j=n_{i}+1$, $R_{i,n_{i}+1}=C_{i}-\sum_{j=1}^{n_{i}}R_{i,j}$, which can be
interpreted as the gap-time between the $n_{i}$-th recurrent event and the end
of follow-up. Define $\delta_{ij}=I(T_{i,j}<C_{i})$ is the failure state
indicator for the $j$-th recurrent event. When $\delta_{i,j}=0$, it indicates
that the observed time is an administrative censoring time.

Assuming that the survival times $R_{i,j},\ldots,R_{L,n_{L}}$ are mutually
independent conditioned on the frailty term $w_{i}$, we can obtain the
likelihood function as:

$\displaystyle\mathcal{L}(\Theta|\mathcal{D},\mathbf{w})$
$\displaystyle=\prod_{i=1}^{L}\left\\{S(r_{i,n_{i}+1}|w_{i})\prod_{j=1}^{n_{i}}f(r_{i,j}|w_{i})\right\\},$

where $\mathcal{D}$ as the set of observed data, such that
$\mathcal{D}=\left\\{r_{i,j},r_{i,n_{i}+1},\delta_{i,j},\mathbf{x}_{i,j};i=1,\ldots,L;j=1,\ldots,n_{i}\right\\}.$
Let
$\Theta=\left\\{\boldsymbol{\gamma},\boldsymbol{\psi},\boldsymbol{\phi},\sigma^{2}_{w}\right\\}$
denote the set of parameters to be estimated in the models.

##  3 Monte Carlo simulation study

The simulation study was conducted to evaluate two scenarios:
($\mathcal{S}_{1}$) individual frailty $\left(n_{i}=1\right)$, and
($\mathcal{S}_{2}$) shared frailty, where individuals experience recurrent
events $\left(n_{i}\geq 1\right)$; $\forall i\in\\{1,\ldots,L\\}$. The data
generation and model fitting were executed in the R programming language (R
Core Team,, 2024). We utilized rstan package Stan Development Team, (2018) to
generate four Markov Chain Monte Carlo (MCMC) chains for each parameter, each
chain comprising 2000 iterations, with 1000 warm-up iterations. This approach
yielded posterior sample sizes of 4000 for each parameter.

We generated $M_{C}=250$ Monte Carlo replicas, each one with $L=300$
individuals using a YP model and an exponential baseline distribution (YPEX)
with rate parameter $\gamma=1.2$. For the individual $i$, we generated the
frailty $w_{i}\sim\mbox{Normal}~{}(0,1)$ and two covariates
$X_{i,1}\sim\mbox{Bernoulli}~{}(0.5)$ and $X_{i,2}\sim\mbox{Normal}~{}(0,1)$.
We generate the administrative censoring $C_{i}\sim U(0,5)$. The gap time of
the $j$-th recurrence was generated by applying the inverse of the survival
function, represented by $r_{i,j}=S^{-1}(u_{i,j}|\mathbf{x}_{i},w_{i})$, where
$U_{i,j}\sim U(0,1)$. See Oliveira, (2024) for more details on how to get
$S^{-1}(u|\mathbf{x},w_{i})$. In $\mathcal{S}_{1}$, we generate $r_{i,j}$ only
once for each $i$. In contrast, in $\mathcal{S}_{2}$, the process of obtaining
$r_{i,j}$ is repeated as long as $t_{i,j}\leq c_{i}$. In both scenarios, the
true values of the parameters established are $\phi_{1}=-1$, $\phi_{2}=2$,
$\psi_{1}=2$, $\psi_{2}=2$.

The following prior distributions were used:
$\boldsymbol{\psi},\boldsymbol{\phi}\sim\mbox{Normal}\left(0,4^{2}\right)$,
$\boldsymbol{\gamma}\sim\mbox{LogNormal}~{}(0,2)$, and
$\sigma_{w}\sim\mbox{Gamma}~{}(0.1,0.1).$ The prior distributions established
for the regression coefficients $\boldsymbol{\psi}$, $\boldsymbol{\phi}$, and
the the standard deviation of the frailty $\sigma_{w}$ are weakly informative
(Stan Development Team,, 2023). For the parameters $\boldsymbol{\gamma}$, we
choose a prior distribution that provides greater stability in the inferential
process, as suggested by (Demarqui et al.,, 2019). Regarding the estimation of
the parameters for the YPBP model, the choice of the BP degree is motivated by
(Osman and Ghosh,, 2012). These authors suggest a range of values, within
which we consider $m=L^{0.4}$ in Scenario $\mathcal{S}_{1}$, and
$m=\left[\sum_{i=1}^{L}n_{i}\right]^{0.4}$ in Scenario $\mathcal{S}_{2}$. We
use the same criteria to define the value of $m$ in models with baseline PE
distribution. All simulation study results are also available in
cassiushenrique.shinyapps.io/appSimulationsFrailty. The statistics used to
evaluate the simulation study are described in Appendix A.

We start by evaluating the Scenario $\mathcal{S}_{1}$. In simulated data,
approximately $68.4\%$ of the individuals experienced the event of interest,
on average. The Monte Carlo summary statistics in Table 1 compare the
performance of the models. The estimates from our models are close to the true
values. The estimated values of $\sigma_{w}$ and $\phi_{2}$ by the YPBP model
are slightly less accurate, deviating more from their true values, compared to
the corresponding estimates obtained from the other models. The generator
model (YPEX) presented smaller biases for the short-term effect coefficients
($\psi_{1}$ and $\psi_{2}$). The YPEX and YPPE models show greater RB for the
parameters of dichotomous variables in comparison to continuous covariates
effects. For all the models, credibility intervals are similar. Furthermore,
the ASE values are close to SDE, and CP values are close to the desired level
of $0.95$, indicating good performance.

Table 1: Scenario $\mathcal{S}_{1}$ \- Monte Carlo summary statistics of the YPEX, YPPE, and YPBP with individual frailties, for $L=300$ and $M_{C}=250$. | 95% CI |
---|---|---
fitted model | par | true | est | RB (%) | ASE | SDE | LW | UP | CP
YPEX | $\phi_{1}$ | -1 | -0.9406 | 5.9395 | 0.2970 | 0.2935 | -1.4732 | -0.3098 | 0.9360
| $\phi_{2}$ | 2 | 2.0795 | 3.9738 | 0.3404 | 0.3388 | 1.4428 | 2.7755 | 0.9440
| $\psi_{1}$ | 2 | 2.0175 | 0.8726 | 0.2941 | 0.2940 | 1.4410 | 2.5919 | 0.9560
| $\psi_{2}$ | 2 | 2.0094 | 0.4704 | 0.1713 | 0.1713 | 1.6731 | 2.3453 | 0.9640
| $\sigma_{w}$ | 1 | 0.9971 | -0.2909 | 0.1763 | 0.1763 | 0.6489 | 1.3418 | 0.9400
YPPE | $\phi_{1}$ | -1 | -0.9399 | 6.0081 | 0.3085 | 0.3049 | -1.4925 | -0.2813 | 0.9520
| $\phi_{2}$ | 2 | 2.0681 | 3.4029 | 0.4606 | 0.4595 | 1.2328 | 3.0187 | 0.9360
| $\psi_{1}$ | 2 | 2.1097 | 5.4830 | 0.3482 | 0.3452 | 1.4506 | 2.8150 | 0.9600
| $\psi_{2}$ | 2 | 2.0802 | 4.0083 | 0.2497 | 0.2481 | 1.6254 | 2.6026 | 0.9360
| $\sigma_{w}$ | 1 | 0.9944 | -0.5598 | 0.3186 | 0.3185 | 0.4260 | 1.6263 | 0.9320
YPBP | $\phi_{1}$ | -1 | -0.9632 | 3.6782 | 0.3076 | 0.3063 | -1.5138 | -0.3090 | 0.9520
| $\phi_{2}$ | 2 | 2.1899 | 9.4949 | 0.5176 | 0.5086 | 1.2725 | 3.2830 | 0.9240
| $\psi_{1}$ | 2 | 2.0925 | 4.6234 | 0.3439 | 0.3418 | 1.4595 | 2.8115 | 0.9640
| $\psi_{2}$ | 2 | 2.0851 | 4.2546 | 0.2574 | 0.2556 | 1.6606 | 2.6741 | 0.9640
| $\sigma_{w}$ | 1 | 1.0628 | 6.2832 | 0.3346 | 0.3306 | 0.4848 | 1.7816 | 0.9440

Now we consider the Scenario $\mathcal{S}_{2}$ whose results are Table 2. The
YPEX shows a high degree of accuracy in estimating the parameters. The
estimates for all the parameters are very close to their true values, with
relative biases (RB). In contrast, when the YPPE is fitted, the relative
biases for the same parameters increase slightly, in magnitude, except for
$\sigma_{w}$. The YPBP demonstrates a further increase in relative biases for
parameters such as $\phi_{2}$. Additionally, the ASE and SDE estimates are
close to each other. The CP values are approximately $0.95$, deviating by no
more than 0.026 from this level in all cases. Notably, across all models, the
estimation of the frailty parameter $\sigma_{w}$ is consistent and accurate,
with relative biases under 2%. This indicates that all three models are
reliable in capturing the shared frailty component, which is an essential
aspect.

Table 2: Scenario $\mathcal{S}_{2}$ \- Monte Carlo summary statistics of the YPEX, YPPE, and YPBP with shared frailties, for $L=300$ and $M_{C}=250$. | 95% CI |
---|---|---
fitted model | par | true | est | RB (%) | ASE | SDE | LW | UP | CP
YPEX | $\phi_{1}$ | -1 | -0.9831 | 1.6948 | 0.2586 | 0.2583 | -1.4664 | -0.4492 | 0.9760
| $\phi_{2}$ | 2 | 2.0664 | 3.3213 | 0.2951 | 0.2940 | 1.5314 | 2.6889 | 0.9560
| $\psi_{1}$ | 2 | 1.9834 | -0.8318 | 0.2090 | 0.2090 | 1.5739 | 2.3930 | 0.9640
| $\psi_{2}$ | 2 | 2.0076 | 0.3822 | 0.1377 | 0.1377 | 1.7418 | 2.2819 | 0.9520
| $\sigma_{w}$ | 1 | 1.0172 | 1.7159 | 0.0911 | 0.0908 | 0.8519 | 1.2085 | 0.9600
YPPE | $\phi_{1}$ | -1 | -0.9692 | 3.0824 | 0.2676 | 0.2666 | -1.4629 | -0.4144 | 0.9760
| $\phi_{2}$ | 2 | 2.0885 | 4.4248 | 0.3237 | 0.3218 | 1.5101 | 2.7747 | 0.9640
| $\psi_{1}$ | 2 | 1.9628 | -1.8625 | 0.2128 | 0.2125 | 1.5461 | 2.3796 | 0.9600
| $\psi_{2}$ | 2 | 1.9958 | -0.2118 | 0.1387 | 0.1386 | 1.7282 | 2.2718 | 0.9440
| $\sigma_{w}$ | 1 | 1.0133 | 1.3316 | 0.0906 | 0.0904 | 0.8491 | 1.2037 | 0.9600
YPBP | $\phi_{1}$ | -1 | -0.9755 | 2.4497 | 0.3115 | 0.3109 | -1.5326 | -0.3207 | 0.9560
| $\phi_{2}$ | 2 | 2.1264 | 6.3217 | 0.3226 | 0.3186 | 1.5465 | 2.8140 | 0.9640
| $\psi_{1}$ | 2 | 1.9396 | -3.0183 | 0.2435 | 0.2426 | 1.4703 | 2.4252 | 0.9360
| $\psi_{2}$ | 2 | 1.9755 | -1.2260 | 0.1447 | 0.1446 | 1.6976 | 2.2649 | 0.9440
| $\sigma_{w}$ | 1 | 1.0099 | 0.9855 | 0.0938 | 0.0937 | 0.8401 | 1.2073 | 0.9720

##  4 Real applications

To illustrate the application of our proposal, we use two databases
readmission and diarrhea. The dataset readmission from the frailtypack package
Rondeau et al., (2012), was previously applied in the study of González et
al., (2005). These data consist of times between hospital readmissions and
time to death for patients with colorectal cancer. Here, we are interested in
modeling only the times until terminal events (death). In this application, we
use a frailty to handle unobserved heterogeneity. The second dataset was first
used by (Barreto et al.,, 1994). It is characterized by the presence of
multiple events because the individuals in the sample had at least one episode
of diarrhea throughout their stay in the study. Our interest is to evaluate
the effects of some covariates on the times between diarrhea recurrences. The
use of frailties, in this case, aims to deal with a possible association
between times until the recurring events. Applications of readmission and
diarrhea will be referred to throughout the text as $\mathcal{A}_{1}$ and
$\mathcal{A}_{2}$, respectively. For the inference procedure, in both
applications, we choose weakly informative prior:
$\boldsymbol{\psi},\boldsymbol{\phi}\sim\mbox{Normal}~{}(0,3);\boldsymbol{\gamma}\sim\mbox{LogNormal}~{}(0,2);\text{
and }\sigma_{w}\sim\mbox{Gamma}~{}(1,1)$. In all models whose baseline is PE
or BP, we use $m=5$.

To compare the quality of the fits of the models in both applications, we
applied the Widely Applicable Information Criterion (WAIC) criterion. The WAIC
is a goodness-of-fit measure for statistical models, especially useful in
Bayesian contexts (Ninomiya,, 2021). It is a generalization of the well-known
Akaike Information Criterion (AIC) and is applicable even when the model is
complex or when the number of parameters is large concerning the number of
observations Akaike, (2011). We applied the function waic available in the
package loo (Vehtari et al.,, 2021). The WAIC calculation steps are in
Appendix B. Among the frailty models fitted in applications $\mathcal{A}_{1}$
and $\mathcal{A}_{2}$ (PHEX, PHPE, PHBP, POEX, POPE, POBP, YPEX, YPPE, and
YPBP), the primary interest is in the model with the best WAIC score, for
which we will report estimates. The WAIC values obtained in applications
$\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ are available in Table 3.

Table 3: WAIC values obtained in applications $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ | WAIC
---|---
Model | $\mathcal{A}_{1}$ | $\mathcal{A}_{2}$
PHEX | 1832.3565 | 2828.3798
PHPE | 1834.1247 | 2663.4610
PHBP | 1819.1796 | 2621.1369
POEX | 1837.6724 | 2657.5942
POPE | 1838.0630 | 2595.1940
POBP | 1829.8205 | 2577.9906
YPEX | 1835.2884 | 2671.5956
YPPE | 1837.6724 | 2642.7094
YPBP | 1818.9457 | 2528.7280

Let’s analyze the application $\mathcal{A}_{1}$. The dataset readmission
originates from Bellvitge’s Public University Hospital in Barcelona, Spain,
capturing medical records from January 1996 to December 1998. The study
focused on 403 individuals who underwent surgery. Some follow-up termination
occurred in cases of patient death, migration, or hospital transfer. Only 103
patients died. Four time-fixed effects recorded in the file will be
considered: (1) sex (Male, when sex = 0, or Female, when sex=1); (2) chemo
which represents whether there was chemotherapy treatment (Treated, when
chemo=1, or nonTreated, when chemo=0) and, (3) dukes which represents the
Dukes’ stage. González et al., (2005) classified the sample pacients as (A-B,
C or D). Table 4 shows how we configure two dummy variables to accommodate the
three levels of Dukes’ stages, A-B, C, and D.

Table 4: Dummy variable for variable dukes. Dukes’ stages | Dukes1 | Dukes2
---|---|---
A-B | 0 | 0
C | 1 | 0
D | 0 | 1

We observed that 164 were women, 239 were men; 217 received chemotherapy
treatment and 186 did not. Colorectal cancer of 180 patients was classified as
Dukes’ stage A-B, of 148 as Dukes’ stage C, and of 75 as Dukes’ stage D.

In $\mathcal{A}_{1}$, we generated four MCMC chains for each parameter via
rstan Stan Development Team, (2018) with 5000 iterations, of which 2500 are
warm-ups, resulting in posterior samples of size 10000. This large volume of
posterior samples aimed to allow a better convergence of $\sigma_{w}$. The
model with the best WAIC is YPBP. By this model, the estimate for the effect
of the variable sex suggests that being female might be associated with a
lower likelihood of death compared to being male, but this result is not
statistically significant, since the credibility interval includes zero. The
variable chemo is significant, only in the short term and its estimate
suggests that undergoing chemotherapy treatment increases the likelihood of
the death of the patient compared those to not undergoing chemotherapy at the
beginning of treatment. The estimate for the effect of the variable Dukes2
indicates that being at Dukes’ stage D in comparison to A-B is associated with
a higher likelihood of death. We also note that the parameter $\sigma_{w}$ is
significant since the credibility interval does not include zero. It suggests
that there is significant variability among study individuals that is not
explained by the model’s fixed effects alone. The estimates of the other
models can be seen by accessing cassiushenrique.shinyapps.io/appRealFrailty.

Table 5: Summary of the YPBP model to the readmission data: posterior mean estimate (est), standard deviation (sd) along with the 95% credibility interval (LW; UP). | 95% CI
---|---
par | description | est | sd | LW | UP
$\psi_{1}$ | Sex | -0.3518 | 0.3286 | -1.0093 | 0.3182
$\psi_{2}$ | Chemo | 1.7940 | 0.5734 | 0.5861 | 2.8310
$\psi_{3}$ | Dukes1 | 0.5240 | 0.5693 | -0.4778 | 1.7884
$\psi_{4}$ | Dukes2 | 3.2558 | 0.6201 | 2.1776 | 4.5962
$\phi_{1}$ | Sex | -0.3694 | 0.8577 | -2.0442 | 1.2857
$\phi_{2}$ | Chemo | -3.2284 | 1.3346 | -4.7150 | 0.6903
$\phi_{3}$ | Dukes1 | 2.1147 | 1.2576 | -0.9656 | 4.5003
$\phi_{4}$ | Dukes2 | 3.8488 | 1.2468 | 0.8991 | 6.3029
$\sigma_{w}$ | sd(Frailty) | 0.8439 | 0.3271 | 0.2258 | 1.4582

From Figure 1, we can notice a similarity between the survival functions
estimated by our model (represented by the continuous curves) and the Kaplan-
Meier estimates. This fact suggests that our model is capturing the pattern of
the observed data. However, is perceived the high volume of censoring in this
dataset. Our models do not include cure fraction models. The inclusion of this
approach could increase the precision of long-term effect estimates. Figure
1-A displays the survival curves for both men and women. This figure suggests
that there is insufficient statistical evidence to conclude that a patient’s
sex significantly influences survival time. Figure 1-B presents that
chemotherapy treatment has a significant short-term effect on patient
mortality, with no apparent long-term impact. Our model estimates that the
survival curves intersect for the first time on the day $1397$. Note that the
survival curves cross more than once, however, the long-term effect tends to
be zero. In Figure 1-C, we see the survival functions for Dukes’ stages A-B,
C, and D present well-defined differences between them. The Dukes’ stage D
shows the lowest survival functions over time, while patients with Dukes’s
stages A-B colorectal cancer appear to survive longer.

Figure 1: Kaplan-Meier (step function) and survival curves estimated by YPBP
model (continuous function) about the terminal event for the levels of
variables (A) sex, (B) chemo, and (C) dukes.

We now discuss the results of the application $\mathcal{A}_{2}$. The diarrhea
dataset originates from a community-based, randomized, placebo-controlled
study spearheaded by Barreto et al., (1994), conducted from December 1990 to
December 1991, and supported by the Institute of Collective Health of the
Federal University of Bahia. This research aimed to evaluate the effects of
vitamin A supplementation on diarrhea incidence among children aged 6 to 48
months in the Northeast of Brazil. The participating children had an average
age of 29.32 months, with a standard deviation of 12.18 months and a median
age of 30 months.

This dataset comprises 860 children, categorized into two groups based on the
covariate treatment: 426 in the placebo group (treatment = 0) and 434 who
received vitamin A supplementation (treatment = 1). In the sample, there were
403 females (sex = 0) and 457 males (sex = 1). Given the occurrence of one or
more diarrhea episodes among various participants during the study period, the
dataset encompasses 5592 records. On average, the children experienced 6.502
episodes of diarrhea, with the median frequency at 5 episodes. The number of
recurrences varied from a minimum of 1 to a maximum of 27 with a standard
deviation of 5.24.

The dataset recorded five covariates, three of which are constant (treatment,
sex, and age) and were documented at the onset of the study. The remaining two
covariates, which vary per episode, include the number of days without
diarrhea before the current episode and the average number of liquid or semi-
liquid stools from the preceding diarrhea episode. However, these time-
dependent covariates are not considered in this particular analysis due to the
limitations of the models discussed, which cannot accommodate such variables.

In $\mathcal{A}_{2}$, four MCMC chains for each parameter via rstan Stan
Development Team, (2018) are also generated, each with 2000 iterations, of
which 1000 are warm-ups, resulting in subsequent samples of size 1000. We
report the estimates from YPBP here because this model has the best WAIC
score. The fit results are in Table 6. The estimates of the other models can
be seen online https://cassiushenrique.shinyapps.io/appRealFrailtyDiarrhea.
The treatment is only significant in the long term because the credibility
interval for $\phi_{1}$ does not include zero. Over the long term, the
treatment tends to reduce the risk of new diarrhea episodes. The sex of the
individual is not significant either in the short or long-term concerning
diarrhea occurrences, because the credibility intervals do not include zero.
The age of the individual seems to be significant both in the short and long
term, in that older children tend to have a lower risk of diarrhea recurrence
both at the start and end of the follow-up. The standard deviation of frailty
is significant, indicating that there is some association between the
recurrent events.

Table 6: Summary of the YPBP model to the diarrhea data: posterior mean estimate (est), standard deviation (sd) along with the 95% credibility interval (LW; UP). | 95% CI
---|---
par | description | est | sd | LW | UP
$\psi_{1}$ | treatment | -0.1741 | 0.0839 | -0.3375 | -0.0070
$\psi_{2}$ | sex | 0.1323 | 0.0824 | -0.0262 | 0.2958
$\psi_{3}$ | Age(beginning) | -0.0177 | 0.0030 | -0.0236 | -0.0120
$\phi_{1}$ | treatment | -0.1290 | 0.0646 | -0.2537 | 0.0016
$\phi_{2}$ | sex | -0.0019 | 0.0640 | -0.1211 | 0.1272
$\phi_{3}$ | Age(beginning) | -0.0368 | 0.0023 | -0.0412 | -0.0324
$\sigma_{w}$ | sd(Frailty) | 0.6350 | 0.0293 | 0.5793 | 0.6934

##  5 Final remarks and future research

This work proposed to develop a class of models within a Bayesian framework,
designed to explain the impact of observed characteristics on survival curves
that may intersect. For this finally, we used the YP regression structure for
its ability to encompass and generalize the PH and PO models. The class of
models embraces YP frailty. The incorporation of frailty in these models
constitutes a contribution of this study, since in the literature the YP
models did not incorporate frailty. We combined exponential, PE, and BP
baseline functions. The selection of these last two baseline functions was
motivated by their versatility because can fit a variety of hazard function
shapes. In that regard, the innovations promoted by this work are the YPEX,
YPPE, YPBP frailty models.

The models of the class enable the analysis of survival data under distinct
scenarios: ($\mathcal{S}_{1}$) individuals with a unique survival time where
individual frailty explains unobserved heterogeneities; ($\mathcal{S}_{2}$)
individuals who experience recurring events for which the shared frailty
accommodates the association between the survival times of the same
individual.

As for numerical results, we executed a Monte Carlo simulation study for the
class aimed at evaluating the influence of model selection on parameter
estimation. This assessment focused on some criteria such as estimation biases
(RB), average standard error (ASE), standard deviation of estimates (SDE),
credibility intervals, and coverage probability (CP). A total of $M_{C}=250$
Monte Carlo replicas were generated, each comprising $L=300$ individuals. In
the two scenarios, our estimates mean and median are generally close to the
true values, indicating a good level of accuracy. In addition, the ASE and SDE
values are close and the CP values are not very far from 95%. These facts
signal a good performance of our models. In the real application, we fitted
our models on the readmission and diarrhea databases and we report the
estimates and interpretations of our models with the best WAIC.

This research has some limitations. In our simulation studies, we did not
apply WAIC to Monte Carlo samples. We acknowledge that the assessment of these
values could enhance the depth of comparative analysis of our models. We did
not consider Weibull baseline distributions in the fit of our models. In the
real application of readmission, we believe that the high volume of
administrative censoring somewhat reduced the predictive ability of our
models. Our models are not yet capable of accommodating time-dependent
variables. In future research, we want to apply the following approaches: (A)
Evaluate the WAIC of our model fits in a simulation study. (B) Conduct a more
extensive simulation study incorporating the baseline Weibull distribution.
(C) Incorporate a cure fraction model into our model classes. (D) Adapt our
regression frameworks to allow us to model time-dependent covariates. (E)
Extend our models to a frequentist approach. (F) Deploy a residual analysis
that provides an additional way of evaluating the quality of our fits. (G)
Publish an R package that provides the functions used in this work,
facilitating the replication of its results.

## Acknowledgments

To the Federal University of Minas Gerais and CAPES, for the financial support
granted during the period of my doctorate.

## Appendix A: Simulation study comparison criteria

Consider a generic parameter, whose true value is $\Phi$, and $\hat{\Phi}_{K}$
the posterior estimate obtained from the $k$-th Monte Carlo replica, which
$k\in\left\\{1,\ldots,M_{C}\right\\}$. The average estimate (est) is given by

$\text{est}(\Phi)=\frac{1}{M_{C}}\sum_{k=1}^{M_{C}}\hat{\Phi}_{k}.$

It is possible to compute bias by

$\text{RB}(\%)=\frac{100}{M_{C}}\sum_{k=1}^{M_{C}}\frac{\hat{\Phi}_{k}-\Phi}{|\Phi|}.$

Additionally, we can compute the average standard error (ASE) of the estimates
by

$\text{ASE}=\frac{1}{M_{C}}\sum_{k=1}^{M_{C}}\text{se}(\hat{\Phi}_{k}),$

where $\text{se}(\hat{\Phi}_{k})$ represents the mean of standard error
estimates of $\Phi$. We are also interested in evaluating the standard
deviation estimate (SDE) de $\Phi$ by

$\text{SDE}=\left\\{\frac{1}{M_{C}-1}\sum_{k=1}^{M_{C}}\left[\hat{\Phi}_{k}-\text{est}(\Phi)\right]^{2}\right\\}^{1/2}.$

In a well-fitted model, we note these characteristics: est$(\Phi)$ should be
close to the true value $\Phi$; SDE and ASE should be similar; RB(%) should
approximate zero; and CP should be close to the pre-defined confidence level
($1-\alpha$). When the ASE $<$ SDE, is expected CP $<1-\alpha$. On the other
hand, if ASE $>$ SDE, it is expected that CP $>1-\alpha$.

## Appendix B: Widely Applicable Information Criterion

The Widely Applicable Information Criterion (WAIC) calculates the log-
likelihood of the data given a model, denoted as $\widehat{\mbox{lppd}}$, and
penalizes for the complexity of this model, considering the effective number
of parameters. Lower WAIC values indicate a model with a better predictive
fit. Unlike AIC, WAIC is based on a weighted average of all parameter
posterior distributions rather than just a point estimate, which can provide a
more robust assessment of model quality (Vehtari et al.,, 2023).
Mathematically, it is an alternative approach to estimating the expected log
pointwise predictive density and is calculated by

$\displaystyle\widehat{\mbox{elppd}}_{\text{waic}}=\widehat{\mbox{lppd}}-\widehat{p}_{\text{waic}},$

where $\widehat{p}_{\text{waic}}$ is the estimated effective number of
parameters and is computed based on the sum of the posterior variance of the
log-likelihood function. In practical terms, we can calculate using the
posterior variance of the log predictive density for each data point $y_{i}$,
i.e.,

$\widehat{p}_{\text{waic}}=V_{s=1}^{S}\left[\log
p\left(y_{i}|\Theta^{(s)}\right)\right],$

in which

$V_{s=1}^{S}a_{s}=\frac{1}{S}\sum_{s=1}^{S}(a_{s}-\bar{a})^{2}.$

Then, we define

$\displaystyle\text{WAIC}=-2~{}\widehat{\mbox{elppd}}_{\text{waic}}.$

To compare the quality of the fits of our models, we applied the function waic
available in the package loo (Vehtari et al.,, 2021).

## References

  * Akaike, (2011) Akaike, H. (2011).  Akaike’s Information Criterion.  In: Lovric, M. International Encyclopedia of Statistical Science. pages 25–25. Springer, Berlin, Heidelberg.
  * Barreto et al., (1994) Barreto, M. L., Farenzena, G., Fiaccone, R., Santos, L., Assis, A. M. d. O., Araújo, M., and Santos, P. (1994).  Effect of vitamin a supplementation on diarrhoea and acute lower-respiratory-tract infections in young children in brazil.  The Lancet, 344(8917):228–231.
  * Bennett, (1983) Bennett, S. (1983).  Analysis of survival data by the proportional odds model.  Statistics in Medicine, 2(2):273–277.
  * Bernstein, (1912) Bernstein, S. (1912).  On the best approximation of continuous functions by polynomials of a given degree.  Kharkov Mathematical Society, 2(13):49–194.
  * Breslow, (1974) Breslow, N. (1974).  Covariance analysis of censored survival data.  Biometrics, 30(1):89–99.
  * Clayton and Cuzick, (1985) Clayton, D. and Cuzick, J. (1985).  Multivariate generalizations of the proportional hazards model.  Journal of the Royal Statistical Society: Series A, 148(2):82–108.
  * Cox, (1972) Cox, D. R. (1972).  Regression models and life-tables.  Journal of the Royal Statistical Society: Series B, 34(2):187–202.
  * Demarqui, (2020) Demarqui, F. N. (2020).  YPPE: Yang and Prentice Model with Piecewise Exponential Baseline Distribution.  R package version 1.0.1.
  * Demarqui et al., (2011) Demarqui, F. N., Dey, D. K., Loschi, R. H., and Colosimo, E. A. (2011).  Modeling survival data using the piecewise exponential model with random time grid.  Recent Advances in Biostatistics: False Discovery Rates, Survival Analysis, and Related Topics, 4(1):109–122.
  * Demarqui et al., (2012) Demarqui, F. N., Loschi, R. H., Dey, D. K., and Colosimo, E. A. (2012).  A class of dynamic piecewise exponential models with random time grid.  Journal of Statistical Planning and Inference, 142(3):728–742.
  * Demarqui and Mayrink, (2021) Demarqui, F. N. and Mayrink, V. D. (2021).  Yang and prentice model with piecewise exponential baseline distribution for modeling lifetime data with crossing survival curves.  Brazilian Journal of Probability and Statistics, 35(1):172–186.
  * Demarqui et al., (2019) Demarqui, F. N., Mayrink, V. D., and Ghosh, S. K. (2019).  An unified semiparametric approach to model lifetime data with crossing survival curves.  arXiv preprint arXiv:1910.04475.
  * Diao et al., (2013) Diao, G., Zeng, D., and Yang, S. (2013).  Efficient semiparametric estimation of short-term and long-term hazard ratios with right-censored data.  Biometrics, 69(4):840–849.
  * Economou and Caroni, (2007) Economou, P. and Caroni, C. (2007).  Parametric proportional odds frailty models.  Communications in Statistics—Simulation and Computation, 36(6):1295–1307.
  * González et al., (2005) González, J. R., Fernandez, E., Moreno, V., Ribes, J., Peris, M., Navarro, M., Cambray, M., and Borràs, J. M. (2005).  Sex differences in hospital readmission among colorectal cancer patients.  Journal of Epidemiology and Community Health, 59(6):506.
  * Gupta and Peng, (2014) Gupta, R. C. and Peng, C. (2014).  Proportional odds frailty model and stochastic comparisons.  Annals of the Institute of Statistical Mathematics, 66(5):897–912.
  * Hanagal, (2011) Hanagal, D. D. (2011).  Modeling survival data using frailty models.  Springer, New York, USA, first edition.
  * Huang and Wang, (2004) Huang, C.-Y. and Wang, M.-C. (2004).  Joint modeling and estimation for recurrent event processes and failure time data.  Journal of the American Statistical Association, 99(468):1153–1165.
  * Kalbfleisch and Prentice, (1973) Kalbfleisch, J. D. and Prentice, R. L. (1973).  Marginal likelihoods based on cox’s regression and life model.  Biometrika, 60(2):267–278.
  * Klein and Moeschberger, (2006) Klein, J. P. and Moeschberger, M. L. (2006).  Survival analysis: techniques for censored and truncated data.  Springer, New York, USA.
  * Lawless, (1987) Lawless, J. F. (1987).  Regression methods for poisson process data.  Journal of the American Statistical Association, 82(399):808–815.
  * Lin and Wang, (2011) Lin, X. and Wang, L. (2011).  Bayesian proportional odds models for analyzing current status data: univariate, clustered, and multivariate.  Communications in Statistics-Simulation and Computation, 40(8):1171–1181.
  * Liu et al., (2004) Liu, L., Wolfe, R. A., and Huang, X. (2004).  Shared frailty models for recurrent events and a terminal event.  Biometrics, 60(3):747–756.
  * Lorentz, (1986) Lorentz, G. G. (1986).  Bernstein polynomials.  American Mathematical Soc., New York, USA.
  * Lorentz, (2012) Lorentz, G. G. (2012).  Bernstein polynomials.  American Mathematical Society, New York, USA.
  * Mazroui et al., (2012) Mazroui, Y., Mathoulin-Pelissier, S., Soubeyran, P., and Rondeau, V. (2012).  General joint frailty model for recurrent event data with a dependent terminal event: application to follicular lymphoma data.  Statistics in Medicine, 31(11-12):1162–1176.
  * Miranda Filho, (2022) Miranda Filho, W. R. (2022).  Semiparametric modeling for multivariate survival data via copulas.  PhD thesis, Universidade Federal de Minas Gerais.  Available online: https://repositorio.ufmg.br/handle/1843/41685 (accessed on 10 November 2023).
  * Ninomiya, (2021) Ninomiya, Y. (2021).  Prior intensified information criterion.  arXiv preprint arXiv:2110.12145.
  * Oliveira, (2024) Oliveira, C. H. X. (2024).  A Class of Semiparametric Joint Frailty-Copula Models for Recurrent Events Subjected to a Terminal Event.  PhD thesis, Universidade Federal de Minas Gerais, Statistics Department.
  * Osman and Ghosh, (2012) Osman, M. and Ghosh, S. K. (2012).  Nonparametric regression models for right-censored data using bernstein polynomials.  Computational Statistics and Data Analysis, 56(3):559–573.
  * Panaro et al., (2020) Panaro, R., Demarqui, F., and Mayrink, V. (2020).  spsurv: Bernstein Polynomial Based Semiparametric Survival Analysis.  R package version 1.0.0.
  * Panaro, (2020) Panaro, R. V. (2020).  spsurv: An R package for semi-parametric survival analysis.  arXiv preprint arXiv:2003.10548.
  * R Core Team, (2024) R Core Team (2024).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rondeau et al., (2012) Rondeau, V., Marzroui, Y., and Gonzalez, J. R. (2012).  frailtypack: an r package for the analysis of correlated survival data with frailty models using penalized likelihood estimation or parametrical estimation.  Journal of Statistical Software, 47:1–28.
  * Schneider et al., (2020) Schneider, S., Demarqui, F. N., Colosimo, E. A., and Mayrink, V. D. (2020).  An approach to model clustered survival data with dependent censoring.  Biometrical Journal, 62(1):157–174.
  * Stan Development Team, (2018) Stan Development Team (2018).  RStan: the R interface to Stan.  R package version 2.18.1.
  * Stan Development Team, (2023) Stan Development Team (2023).  Stan modeling language users guide and reference manual: Prior choice recommendations.  Available online: https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations (accessed on 10 May 2023).
  * Vaupel et al., (1979) Vaupel, J. W., Manton, K. G., and Stallard, E. (1979).  The impact of heterogeneity in individual frailty on the dynamics of mortality.  Demography, 16(3):439–454.
  * Vehtari et al., (2023) Vehtari, A., Gabry, J., Magnusson, M., Yao, Y., Bürkner, P.-C., Paananen, T., and Gelman, A. (2023).  loo: Efficient leave-one-out cross-validation and waic for bayesian models.  R package version 2.6.0.
  * Vehtari et al., (2021) Vehtari, A., Gelman, A., Gabry, J., and Yao, Y. (2021).  Package ‘loo’.  Efficient Leave-One-Out Cross-Validation and WAIC for Bayesian Models.  R package version 2.6.0.
  * Wang, (2013) Wang, Y. (2013).  Sample Size Calculation Based on the Semiparametric Analysis of Short-term and Long-term Hazard Ratios.  PhD thesis, Columbia University, Statistics Department.  Available online: https://academiccommons.columbia.edu/doi/10.7916/D8ST7X25 (accessed on 13 November 2022).
  * Yang and Prentice, (2005) Yang, S. and Prentice, R. (2005).  Semiparametric analysis of short-term and long-term hazard ratios with two-sample survival data.  Biometrika, 92(1):1–17.