Bayesian significance test for discriminating between survival distributions

TL;DR

This paper evaluates the Fully Bayesian Significance Test (FBST) for selecting the best survival distribution among lognormal, gamma, and Weibull by using a mixture model and applying it to both simulated and real patient data.

Contribution

It introduces a Bayesian mixture model with reparametrized distributions for survival analysis and applies FBST to discriminate among competing survival distributions.

Findings

FBST effectively distinguishes between survival distributions in simulations.

The mixture model provides clear posterior weights indicating the most suitable distribution.

Application to real patient data demonstrates practical utility in medical survival analysis.

Abstract

An evaluation of FBST, Fully Bayesian Significance Test, restricted to survival models is the main objective of the present paper. A Survival distribution should be chosen among the tree celebrated ones, lognormal, gamma, and Weibull. For this discrimination, a linear mixture of the three distributions, for which the mixture weights are defined by a Dirichlet distribution of order three, is an important tool: the FBST is used to test the hypotheses defined on the mixture weights space. Another feature of the paper is that all three distributions are reparametrized in that all the six parameters - two for each distribution - are written as functions of the mean and the variance of the population been studied. Note that the three distributions share the same two parameters in the mixture model. The mixture density has then four parameters, the same two for the three discriminating densities and two for the mixture weights. Some numerical results from simulations with some right-censored data are considered. The lognormal-gamma-Weibull model is also applied to a real study with dataset being composed by patient's survival times of patients in the end-stage of chronic kidney failure subjected to hemodialysis procedures; data from Rio de Janeiro hospitals. The posterior density of the weights indicates an order of the mixture weights and the FBST is used for discriminating between the three survival distributions. Keywords: Model choice; Separate Models; Survival distributions; Mixture model; Significance test; FBST

Tables

Table 1. Table 1: Mean of estimates for LGW model parameters and percentages of correct decisions made by FBST on selecting the true distribution used to generate the survival times, using samples with different right-censoring percentages

$\%$ of Rc†	Model	$n$	$\mu$	${\sigma }^2$	$p_L$	$p_G$	$p_W$	$\%$ of Cd‡
			$20$	$50$	-	-	-
$10$	Lognormal	$100$	$19.92$	$48.82$	$0.48$	$0.36$	$0.16$	$84$
		$200$	$20.02$	$48.87$	$0.59$	$0.30$	$0.11$	$84$
		$300$	$19.96$	$48.43$	$0.63$	$0.28$	$0.09$	$84$
		$500$	$19.97$	$48.14$	$0.69$	$0.25$	$0.07$	$93$
	Gamma	$100$	$20.07$	$52.95$	$0.38$	$0.36$	$0.26$	$45$
		$200$	$20.01$	$50.60$	$0.38$	$0.41$	$0.21$	$53$
		$300$	$20.06$	$50.90$	$0.36$	$0.44$	$0.20$	$57$
		$500$	$20.05$	$50.91$	$0.34$	$0.48$	$0.18$	$69$
	Weibull	$100$	$20.17$	$52.06$	$0.19$	$0.26$	$0.55$	$86$
		$200$	$20.04$	$51.27$	$0.15$	$0.22$	$0.63$	$94$
		$300$	$19.76$	$50.24$	$0.12$	$0.18$	$0.70$	$100$
		$500$	$19.97$	$51.26$	$0.09$	$0.13$	$0.78$	$100$
$30$	Lognormal	$100$	$20.03$	$49.20$	$0.47$	$0.35$	$0.19$	$48$
		$200$	$20.01$	$48.71$	$0.55$	$0.33$	$0.13$	$63$
		$300$	$20.02$	$48.20$	$0.59$	$0.31$	$0.10$	$70$
		$500$	$19.97$	$47.06$	$0.64$	$0.28$	$0.08$	$86$
	Gamma	$100$	$20.13$	$53.03$	$0.38$	$0.36$	$0.26$	$35$
		$200$	$19.96$	$50.67$	$0.40$	$0.38$	$0.22$	$47$
		$300$	$20.12$	$55.09$	$0.41$	$0.41$	$0.18$	$51$
		$500$	$20.00$	$50.77$	$0.35$	$0.47$	$0.18$	$70$
	Weibull	$100$	$20.07$	$54.25$	$0.21$	$0.28$	$0.51$	$81$
		$200$	$20.11$	$52.35$	$0.16$	$0.23$	$0.61$	$92$
		$300$	$19.97$	$50.77$	$0.15$	$0.21$	$0.64$	$100$
		$500$	$19.97$	$50.78$	$0.10$	$0.15$	$0.75$	$100$
$50$	Lognormal	$100$	$19.98$	$48.14$	$0.45$	$0.35$	$0.20$	$40$
		$200$	$19.91$	$45.79$	$0.50$	$0.35$	$0.15$	$50$
		$300$	$19.97$	$47.14$	$0.56$	$0.33$	$0.11$	$65$
		$500$	$19.92$	$46.82$	$0.63$	$0.29$	$0.08$	$78$
	Gamma	$100$	$20.06$	$53.77$	$0.36$	$0.36$	$0.28$	$31$
		$200$	$19.99$	$51.17$	$0.37$	$0.38$	$0.25$	$43$
		$300$	$20.10$	$52.55$	$0.37$	$0.41$	$0.22$	$47$
		$500$	$20.08$	$51.65$	$0.40$	$0.43$	$0.17$	$57$
	Weibull	$100$	$20.27$	$58.79$	$0.24$	$0.30$	$0.46$	$74$
		$200$	$20.09$	$54.65$	$0.18$	$0.26$	$0.56$	$90$
		$300$	$20.05$	$53.08$	$0.15$	$0.23$	$0.62$	$95$
		$500$	$19.92$	$51.21$	$0.13$	$0.19$	$0.68$	$95$
$†$ percentage of right-censoring
$‡$ percentage of correct decision

Table 2. Table 2: Measures of evidence provided by HD data

Comparison	Null hypothesis	Evidence in favor of null hypothesis
		e-value	p-value∗
$H_L\times H_W$	$H_L$	$0.874$	$0.404$
	$H_W$	$0.043$	$0.004$
$H_L\times H_G$	$H_L$	$0.901$	$0.446$
	$H_G$	$0.757$	$0.277$
$H_G\times H_W$	$H_G$	$1.000$	$1.000$
	$H_W$	$0.113$	$0.015$
*p-value calculated according to Diniz et al. (2012)

Table 3. Table 3: Summary of the posterior distribution of the LGW parameters

Parameter	Mean	SD	$2.5\%$	Median	$97.5\%$
$p_1\text{-lognormal}$	$0.574$	$0.299$	$0.028$	$0.656$	$0.957$
$p_2\text{-gamma}$	$0.283$	$0.256$	$0.005$	$0.187$	$0.845$
$p_3\text{-Weibull}$	$0.143$	$0.157$	$0.003$	$0.087$	$0.606$
$\mu$	$18.537$	$1.441$	$15.819$	$18.561$	$21.079$
${\sigma }^2$	$204.744$	$78.788$	$84.547$	$197.103$	$339.416$

Table 4. Table 4: Hypothesis testing on the mixture weights of LGW model

Hip tese	e-valor	p-valor∗
$p_1=0$	$0.009$	$0.000$
$p_2=0$	$0.656$	$0.119$
$p_3=0$	$0.878$	$0.273$
*p-value calculated according to Diniz et al. (2012)

Table 5. Table 5: Summary of the posterior distribution of lognormal parameters

Par metro	Mean	SD	$2.5\%$	Median	$97.5\%$
$\mu$	$20.298$	$1.543$	$17.493$	$20.251$	$23.285$
${\sigma }^2$	$355.367$	$101.927$	$191.546$	$343.222$	$589.779$

Equations

$$S(t)=P(T>t)=1-F(t),\text{mbox}fort>0,$$

$$f(t)=\frac{dF(t)}{dt}=\frac{d\{1-S(t)\}}{dt}=-\frac{dS(t)}{dt}.$$

$$h(t)=\underset{\Delta t\to 0}{\lim }\frac{\text{mbox}P(t\le T<t+\Delta t|T\ge t)}{\Delta t}=\frac{f(t)}{S(t)}.$$

$$H(t)={\int }_0^{t}h(u)d(u).$$

$$f_L(t|\alpha )$$

$$S_L(t|\alpha )$$

$$=1-\Phi \left[\frac{(\log t-{\alpha }_1)}{\sqrt{{\alpha }_2}}\right];$$

$$h_L(t|\alpha )$$

$$f_G(t|\gamma )$$

$$\text{vspace}\ast 0.8cm{S}_G(t|\gamma )$$

$$\text{vspace}\ast 0.8cm{h}_G(t|\gamma )$$

$$f_W(t|\beta )$$

$$\text{vspace}\ast 0.8cm{S}_W(t|\beta )$$

$$\text{vspace}\ast 0.8cm{h}_W(t|\beta )$$

$$q^{\ast }=\underset{H}{\sup }q(\theta |y)\text{mbox}andT(y)=\{\theta :q(\theta |y)>q^{\ast }\}.$$

$$\overline{ev}=\text{mbox}Pr(\theta \in T(y)|y)={\int }_{T(y)}q(\theta |y)d\theta ;\text{mbox}consequently,ev=1-\overline{ev}.$$

$$f(y_j|\mathbf{\theta })=p_1{f}_1(y_j|{\psi }_1)+\dots +p_m{f}_m(y_j|{\psi }_m)p_k\ge 0,\sum _{k=1}^m{p}_k=1.$$

$$L(y,\mathbf{\theta })=\prod _{j=1}^n\sum _{k=1}^m{p}_k{f}_k(y_j|\mu ,\sigma ).$$

$$f_L(y|{\alpha }_1,{\alpha }_2)=\frac{1}{y\sqrt{2\pi {\alpha }_2}}\exp \left\{-\frac{(\log y-{\alpha }_1{)}^2}{2{\alpha }_2}\right\}.$$

$$\left\{\begin{array}[]{ll}\vspace*{0.3 cm}\mu=E(y|\alpha_{1},\alpha_{2})=\mbox{e}^{\alpha_{1}+\alpha_{2}/2}\\ \sigma^{2}=Var(y|\alpha_{1},\alpha_{2})=(\mbox{e}^{\alpha_{2}}-1)\mbox{e}^{2\alpha_{1}+\alpha_{2}}\end{array}\right.\Rightarrow\left\{\begin{array}[]{ll}\vspace*{0.3 cm}\alpha_{1}=\log\frac{\mu^{2}}{\sqrt{\mu^{2}+\sigma^{2}}}\\ \alpha_{2}=\sqrt{\log\frac{\mu^{2}+\sigma^{2}}{\mu^{2}}}.\end{array}\right.$$

$$f_G(y|{\gamma }_1,{\gamma }_2)=\frac{1}{\Gamma ({\gamma }_2){\gamma }_1^{{\gamma }_2}}{y}^{{\gamma }_2-1}\exp \left\{-\frac{y}{{\gamma }_1}\right\}.$$

$$\left\{\begin{array}[]{ll}\vspace*{0.3 cm}\mu=E(y|\gamma_{1},\gamma_{2})=\gamma_{1}\gamma_{2}\\ \sigma^{2}=Var(y|\gamma_{1},\gamma_{2})=\gamma_{2}\gamma_{1}^{2}\end{array}\right.\Rightarrow\left\{\begin{array}[]{ll}\vspace*{0.3 cm}\gamma_{1}=\frac{\sigma^{2}}{\mu}\\ \gamma_{2}=\frac{\mu^{2}}{\sigma^{2}}.\end{array}\right.$$

$$f_W(y|{\beta }_1,{\beta }_2)=\frac{{\beta }_2}{{\beta }_1^{{\beta }_2}}{y}^{{\beta }_2-1}\exp \left\{-{\left(\frac{y}{{\beta }_1}\right)}^{{\beta }_2}\right\},$$

$$\displaystyle\left\{\begin{array}[]{ll}\vspace*{0.3 cm}\mu=E(y|\beta_{1},\beta_{2})=\beta_{1}\Gamma(1+1/\beta_{2})\\ \sigma^{2}=Var(y|\beta_{1},\beta_{2})=\beta_{1}^{2}\Gamma(1+2/\beta_{2})-\beta_{1}^{2}\Gamma^{2}(1+1/\beta_{2})\end{array}\right.$$

$$\Rightarrow$$

$$L(y,\mathbf{\theta })$$

$$\propto \prod _j^n[f(y_j|\mathbf{\theta }){]}^{{\delta }_j}[S(y_j|\mathbf{\theta }){]}^{1-{\delta }_j}$$

$$\propto \prod _j^n{\left[\sum _{k=1}^m{p}_k{f}_k(y_j|\mu ,\sigma )\right]}^{{\delta }_j}{\left[\sum _{k=1}^m{p}_k{S}_k(y_j|\mu ,\sigma )\right]}^{1-{\delta }_j},$$

$$f(\mathbf{\theta }|y)\propto L(y,\mathbf{\theta })\pi (\mathbf{\theta }).$$

$$H_k:p_k=1\wedge p_i=0,i\ne k.$$

$$H:p_k=0\wedge \sum _{i\ne k}{p}_i=1.$$

$$p_c$$

$$=\text{mbox}Pr(C\le T\le \infty ,0\le C\le \infty )$$

$$=1-\text{mbox}Pr(0\le T\le C,0\le C\le \infty )$$

$$=1-{\int }_0^{\infty }g(c|\lambda ){\int }_0^c{f}_L(t|\mu ,\sigma )dtdc$$

$$=1-{\int }_0^{\infty }g(c|\lambda )F_L(c|\mu ,\sigma )dc,$$

$$f(y|\mathbf{p},\mu ,\sigma )=p_1{f}_L(y|\mu ,\sigma )+p_2{f}_G(y|\mu ,\sigma )+p_3{f}_W(y|\mu ,\sigma ),$$

Full text

Keywords: Model choice; Separate Models; Survival distributions; Mixture model; Significance test; FBST

1 Introduction

In many scientific disciplines, researchers are constantly faced with the fundamental problem of choosing among alternative statistical models. The Neyman-Pearson theory of hypothesis testing applies only if the models belong to the same family of distributions. Alternatively, special procedures are required if the models belong to families that are separate (or non-nested) in the sense that an arbitrary member of one family cannot be obtained as a limit of members of the other. The set of separate families of probability distributions includes the ones used here: lognormal, gamma, and Weibull models (Pereira, 1981; Araujo and Pereira, 2007; Pereira and Pereira, 2017) which have been used widely to describe survival data (Lawless, 2002; Lee and Wang, 2003).

A considerable amount of research on separate families of hypotheses has been realized since the fundamental work of Cox (1961, 1962), who first dealt with the problem. For reviews and references, see (2005); Araujo and Pereira (2007); and Pereira and Pereira (2017).

The Fully Bayesian Significance Test (FBST) introduced by Pereira and Stern (1999) is an alternative test to the ones that are based on Bayes factor or on the classical p-value; mostly for the case of precise hypotheses. The basis for the FBST is an index known as e-value (e stands for evidence) that measures the inconsistency of the hypothesis. For this, it considers the tangent set, $T$; the set of all parameter values for which their posterior density values are greater than the values of the posterior densities of all points that attend the hypothesis. For reviews and further references on FBST, see Pereira et al. (2008) and Stern and Pereira (2014). For a few interesting applications illustrating the use of e-values and the FBST to practical problems, see Diniz et al. (2012), Lauretto at al. (2003), Lauretto at al. (2007), and Pereira and Stern (1999).

In the present work, we consider the FBST for discriminating between the lognormal, gamma and Weibull distributions. We formulate this problem in the context of linear mixture model, as suggested by Cox (1961). It means that, the models under comparison are considered as components of a finite mixture model. The FBST is used for testing hypotheses defined on the mixture weights space. The e-value is the complementary of the posterior probability of the tangent set $T$; $ev=1-Pr(T|Data)$,

Additionally, the density functions of the mixture components are reparametrized in terms of the mean $\mu$ and the variance $\sigma^{2}$ of the population. Hence, the models under discrimination share common parameters (Kamary et al., 2014; Pereira and Pereira, 2017). A standard Bayesian approach to finite mixture models is to consider different pairs of parameters for each of these models and to adopt independent prior distributions for each pair of parameters and a Dirichlet prior on the mixture weights (Lauretto and Stern, 2005; Lauretto at al., 2007). However, since the comparison between the models is based on the same dataset and on the same sample, we believe that it would be inappropriate to consider different means and variances for these models. Note that this reparameterization reduces the number of the parameters to be estimated: in our case, including the weights, from eight to only four.

To illustrate the procedure, numerical results based on simulated right-censored survival times were considered. Also, a real example is introduced to use the lognormal-gamma-Weibull mixture model to the dataset of patients, from Rio de Janeiro hospitals, with end-stage chronic kidney failure who received hemodialysis.

Section 2 presents a brief review of basic concepts and notation for survival analysis. The parametric distributions used in this paper are also described. Section 3 reviews the basic concepts o FBST. Section 4 discusses the FBST formulation for discriminating between survival distributions in the context of mixture models. Section 5 presents the results of the simulation study. Section 6 is about the use of the lognormal-gamma-Weibull on the real dataset. Final remarks are presented in Section 7.

2 Survival analysis

2.1 Basic concepts and notation

Survival analysis is concerned with the analysis of time to occurrence of a certain event of interest, such as failure, death, relapse or development of a given disease.

Let $T$ be a non-negative random variable representing the time until some event of interest. There are three functions of primary interest used to characterize the distribution of $T$, namely the survival function, the probability density function and the hazard function (Lee and Wang, 2003).

The survival function, denoted by $S(t)$, is defined as the probability that an individual survives beyond time $t$:

[TABLE]

where $F(t)$ is the distribution function of $T$. Note that $S(t)$ is a nonincreasing continuous function of time $t$ with $S(0)=1$ and $S(\infty)=\displaystyle\lim_{\infty}S(t)=0$.

The probability density function, denoted by $f(t)$, is the probability of failure in a small interval per unit time. It can be expressed as

[TABLE]

The hazard function, denoted by $h(t)$, represents the probability of failure during a very small time interval, assuming that the individual has survived to the beginning of the interval:

[TABLE]

This function is also known as the conditional failure rate. The cumulative hazard function is defined as

[TABLE]

Therefore, when $t=0$ then, $S(t)=1$ and $H(t)=0$; and when $t=\infty$ then, $S(t)=0$ and $H(t)=\infty$.

2.2 Parametric survival distributions

In this paper, we consider the the FBST for discriminating between the lognormal, gamma and Weibull distributions which are most frequently used in modeling survival data (Lawless, 2002; Lee and Wang, 2003). The probability density functions, the survival functions and the hazard functions of these distributions are highlighted below.

• i)

Let $T$ be a lognormal random variable with parameters $\alpha=(\alpha_{1},\alpha_{2})$, denoted by $T\sim LN(\alpha_{1},\alpha_{2})$,

[TABLE]

• ii)

If $T$ has a Gamma distribution with parameters $\gamma=(\gamma_{1},\gamma_{2})$, denoted by $T\sim G(\gamma_{1},\gamma_{2})$, then

[TABLE]

• iii)

If $T$ has a Weibull distribution with parameters $\beta=(\beta_{1},\beta_{2})$, denoted by $T\sim W(\beta_{1},\beta_{2})$, then

[TABLE]

3 Fully Bayesian Significance Test (FBST)

The FBST of Pereira and Stern (1999), which is reviewed in Pereira et al. (2008), is a Bayesian version of significance testing, as considered by Cox (1977) and Kempthorne (1976), for precise (or sharp) hypotheses.

First, let us consider a real parameter $\theta$, a point in the parameter space $\Theta\subset\Re$, and an observation $y$ of the random variable $Y$. A frequentist looks for the set $I\in\Re$ of sample points that are at least as inconsistent with the hypothesis as $y$ is. A Bayesian looks for the tangential set $T(y)\subset\Theta$ (Pereira et al., 2008), which is a set of parameter points that are more consistent with the observed $y$ than the hypothesis is. An example of a sharp hypothesis in a parameter space of the real line is of the type $H:\theta=\theta_{0}$. The evidence value in favor of $H$ for a frequentist is the usual p-value, $P(Y\in I|\theta_{0})$, whereas for a Bayesian, the evidence in favor of $H$ is the e-value, $ev=1-\mbox{Pr}(\theta\in T(y)|y)$.

In the general case of multiple parameters, $\Theta\subset\Re^{k}$, let the posterior distribution for $\theta$ given $y$ be denoted by $q(\theta|y)\propto\pi(\theta)L(y,\theta)$, where $\pi(\theta)$ is the prior probability density of $\theta$ and $L(y,\theta)$ is the likelihood function. In this case, a sharp hypothesis is of the type $H:\theta\in\Theta_{H}\subset\Theta$, where $\Theta_{H}$ is a subspace of smaller dimension than $\Theta$. Letting $\displaystyle\sup_{H}$ denote the supremum of $\Theta_{H}$, we define the general Bayesian evidence and the tangential set, $T(y)$, as follows:

[TABLE]

The Bayesian evidence value against $H$ is the posterior probability of $T(y)$,

[TABLE]

It is important to note that evidence that favors $H$ is not evidence against the alternative, $\overline{H}=\Theta\setminus H$, because it is not a sharp hypothesis. This interpretation also holds for p-values in the frequentist paradigm. As in Pereira et al. (2008), we would like to point out that this Bayesian significance index uses only the posterior distribution, with no need for additional artifacts such as the inclusion of positive prior probabilities for the hypotheses or the elimination of nuisance parameters. The computation of the e-values does not require asymptotic methods, and the only technical tools needed are numerical optimization and integration methods.

4 Mixture of survival models

Let us consider a dataset $y=\{y_{1},\ldots,y_{n}\}$ and $m$ alternative parametric survival distributions with densities $f_{1}(y|\psi_{1}),f_{2}(y|\psi_{2}),\ldots,f_{m}(y|\psi_{m})$. Here, $\psi_{k},k=1,\ldots,m$, are unknown (vector) parameters and the families of distributions are separate. The problem of interest is to measure the evidence in favor of each model for fitting the dataset. As suggested by Cox (1961), we can consider a general model including all candidate distributions where the choice of a specific distribution is a special case. In this work, we formulate the FBST for the linear mixture of the survival models as a selection procedure. Denoting $\boldsymbol{\theta}=(\psi_{1},\ldots,\psi_{m},\boldsymbol{p})$, the density function for $m-$component mixture model is

[TABLE]

where $\boldsymbol{p}=(p_{1},\ldots,p_{m})$ is the vector of the mixture weights.

In the presente work, the density functions of the mixture components in (4.1) are reparametrized in terms of the mean $\mu$ and the variance $\sigma^{2}$ of the population. Hence, the models under comparison share common parameters (Kamary et al., 2014; Pereira and Pereira, 2017). The main reason for this reparametrization is that, since the comparison between the models is based on the same dataset and on the same sample, we believe that it would be inappropriate to consider different means and variances for these models as is commonly performed in traditional Bayesian approach to finite mixture model. Therefore, we have $\boldsymbol{\theta}=(\mu,\sigma^{2},\boldsymbol{p})$ denoting all parameters of the mixture model, where $\mu$ and $\sigma^{2}$ are the connecting parameters, with $\boldsymbol{p}$ corresponding to the vector of the mixture weights.

Assuming that the $y_{i}$ are conditionally (on the parameter) independent, the likelihood function is defined as

[TABLE]

The families of distributions considered include the lognormal, gamma and Weibull models. Hence, the relationship between the parameters of these models through the $\mu$ and $\sigma^{2}$ is described as follows.

• (i)

Let $y$ be a $\mbox{lognormal}(\alpha_{1},\alpha_{2}),\alpha_{1}\in\ \mathbb{R}\ \mbox{and}\alpha_{2}>0$, with probability density function

[TABLE]

We then have

[TABLE]

• (ii)

Let $y$ be a $\mbox{gamma}(\gamma_{1},\gamma_{2}),\gamma_{1}>0\ \mbox{and}\ \gamma_{2}>0$, with probability density function

[TABLE]

Therefore

[TABLE]

• (iii)

When $y\sim\mbox{Weibull}(\beta_{1},\beta_{2}),\beta_{1}>0\ \mbox{and}\ \beta_{2}>0$, with probability density function

[TABLE]

then

[TABLE]

In order to find $\beta_{2}$, the Newton-Rapson method can be used to solve the nonlinear equation. Here, we use the nleqslv function in the R package of the same name.

A special feature of survival data is that survival times are frequently censored. The survival time of an individual is said to be censored when the event of interest has not been observed for that individual, but is known only to occur in a certain period of time. There are various categories of censoring, such as right censoring, left censoring and interval censoring (see Klein and Moeschberger (2003) for more details). In this paper, we restrict ourselves to data in which the survival times are subject to right censoring, which is the most common censoring mechanism in medical research.

In the model for right-censored data, it is convenient to consider the following notation. Each individual $j$ is assumed to have an event time $T_{j}$ and a censoring time $C_{j}$. The observations consist of $(y_{1},\delta_{1}),(y_{2},\delta_{2}),\ldots,(y_{n},\delta_{n})$, where $y_{j}=\min\{T_{j},C_{j}\}$ and $\delta_{j}=\mathbb{I}(T_{j}\leq C_{j})$, indicating whether $T_{j}$ was observed ($\delta_{j}=1$) or not ($\delta_{j}=0$).

Note that the likelihood function given by (4.2) is for uncensored (or exact) observations. Assuming noninformative censoring, i.e, independence between $T_{j}$ and $C_{j}$, then, the likelihood function for right-censored observations is

[TABLE]

where, $S_{k}$ is the survival function associated with the mixture component $k$.

Assuming independence, the joint prior density function of $\boldsymbol{\theta}=(\mu,\sigma^{2},\boldsymbol{p})$ is given by $\pi(\boldsymbol{\theta})=\pi_{1}(\boldsymbol{p})\pi_{2}(\mu)\pi_{3}(\sigma^{2})$. Therefore, according to the Bayesian paradigm, the posterior density of $\boldsymbol{\theta}$ is

[TABLE]

In this paper, the prior distributions for the connecting parameters, $\mu$ and $\sigma^{2}$, are assumed to be independent gamma distributions, both with a mean of one and a variance of 100, that is, $\mu,\sigma^{2}\sim gamma(0.01,100)$ (Pereira and Pereira, 2017). For the mixture weights, we use a Dirichlet prior, $\boldsymbol{p}\sim Dir(1,1,1)$ when all families of models are considered ($m=3$) or a Beta prior with parameters (1,1) (uniform$(0,1)$) for any combination of $m=2$.

In order to measure the evidence in favour of each model, the hypotheses on the mixture weights are tested (Kamary et al., 2014; Pereira and Pereira, 2017).

The hypothesis specifying that $y$ has the density function $f_{k}(y|\psi_{k})$ is equivalent to

[TABLE]

On the other hand, the hypothesis that $y$ has not the density $f_{k}(y|\psi_{k})$ is equivalent to

[TABLE]

The alternative hypotheses to (4.13) and (4.14) are $A_{k}:p_{k}<1$ and $A_{k}:p_{k}>0$, respectively, which are not sharp anyway.

The FBST procedure is used to test $H_{k},k=1,\ldots,m$, according to the expressions (3.1) and (3.2). For the optimization step, we used the conjugate gradient method (Fletcher and Reeves, 1964). In order to perform the integration over the posterior measure, we used an Adaptive Metropolis Markov chain Monte Carlo algorithm (MCMC) of Haario et al. (2001).

In this paper, the implementation of the Bayesian models is carried out using LaplacesDemon R package. The LaplacesDemon is an open-source package that provides a complete environment for simulation in Bayesian inference (Statisticat, LCC, 2016).

5 Simulations

In this section we present some numerical results based on simulated right-censored survival times in order to evaluate the performance of the FBST for discriminating between the survival distributions via lognormal-gamma-Weibull mixture model (LGW). The main purpose is to measure the convergence rate of correct decisions, concerning the identification of the true model used to generate the survival times $T$.

The simulations of this paper were performed on a Intel(R) Core(TM) i7-5500U CPU@ 2.40GHz computer.

5.1 Simulation scheme of sample points

Let $H_{L}$, $H_{G}$ and $H_{W}$ be the hypotheses specifying the probability density functions of the lognormal, gamma and Weibull distributions, respectively. From each distribution, we generate $200$ samples of sizes $n=100$, $200$, $300$, and $500$. Each sample contain a desired proportion of right-censored observations.

The steps used to simulate a sample, $y$, of size $n$, in which part of the observations is right-censored, are shown below. For this example, we assume that the true survival times has a lognormal distribution.

Assign values to parameters $\mu$ e $\sigma^{2}$; 2. 2.

Calculate the lognormal parameters $(\alpha_{1},\alpha_{2})$ using the expressions (4.3); 3. 3.

For $j=1,\ldots,n$,

• •

Generate the survival time $T_{j}$ from $lognormal(\alpha_{1},\alpha_{2})$;

• •

Generate the right-censoring time $C_{j}$ from a exponential distribution, i.e, $C_{j}\sim Exp(\lambda)$, where the parameter $\lambda$ is chosen such that approximately a desired percentage of simulated observations are right-censored;

• •

Obtain the observed time $y_{j}=\min\{T_{j},C_{j}\}$

• •

Create an indicator random variable $\delta_{j}=\mathbb{I}(T_{j}\leq C_{j})$

Using this generated sample, we obtain the posterior samples for the mixture parameters from Adaptive Metropolis algorithm and we use the FBST to calculate the evidence measures in favor of each model.

The value for the censoring distribution parameter, $\lambda$, is determined by numerical methods (Wan, 2017). We let $p_{c}$ denote the right-censoring probability. We suppose that the censoring time $C$ has exponencial density function $g(c|\lambda)$ and the independence assumption between $T$ and $C$ holds. In order to simulate a sample with approximately $p_{c}\%$ of right-censored observations, the value of $\lambda$ is obtained by solving the following equation:

[TABLE]

where $f_{L}$ and $F_{L}$ are the lognormal probability density and distribution functions of survival times, respectively.

For generating right-censored survival times from the gamma and Weibull distributions, an analogous procedure to that used for the lognormal distribution is employed.

5.2 Criteria for evaluating the performance of the FBST

In order to evaluate the performance of the FBST on selecting the true distribution used to generate the survival times, we have compared the measures of evidence in favor of the hypotheses $H:p_{k}=0$ and $H:p_{k}=1$, $k=L,G,W$, where $p_{k}$ are respectively the mixture weights associated with the lognormal, gamma and Weibull components in the LGW mixture model.

For instance, suppose again that the true survival time has a lognormal distribution. We consider that the FBST has made a correct choice on the LGW model, if the evidence in favor of $H:p_{L}=0$ is less than that in favor of $H:p_{G}=0$ and $H:p_{W}=0$, and the evidence in favor of $H:p_{L}=1$ is greater than that in favor of $H:p_{G}=1$ e $H:p_{W}=1$.

The calculation of the proportions of correct decisions made by FBST is based on $200$ replicates. In these simulations, we have assigned $\mu=20$ and $\sigma^{2}=50$. The FBST procedure is evaluated considering the samples with different censoring percentages: $10\%$, $30\%$ and $50\%$.

5.3 Simulation results

Table 1 presents the mean of the estimates for the LGW mixture model parameters and the percentages of correct decisions made by FBST on selecting the true distribution used to generate the survival times. It is observed that, regardless of the distribution used for generating the survival times and the sample sizes, the estimates for the mean $\mu$ are very close to each other and to the true value of the parameter. For the estimates of the variance $\sigma^{2}$, we observe a variation between them but, in general, they approach the true value of the parameter as the sample size increases.

It is observed that the FBST presents a high performance on identifying the Weibull distribution as the true data generation process and low performance on identifying the gamma distribution. This happens because, regarding the parameters chosen for these simulations, the gamma and lognormal densities are very similar. The general pattern of the simulation results shows that the FBST achieves good performance even for samples with $50\%$ right-censoring.

6 Application: Choice of a survival model for patients with end-stage kidney disease

6.1 Dataset

The dataset used in this paper refers to a cohort study of 473 patients with end-stage chronic kidney failure who received hemodialysis (HD) in four centers in the State of Rio de Janeiro, Brazil. The patients were followed up $11$ years. The observed time for each patient was the number of months from admission to hemodialysis until death or the end of the observation period (kidney transplant or end of the study) which indicates a right-censored survival time. For a complete description of this dataset, see (2014).

In this paper, our main interest is to apply the LGW model to the survival data for HD patients and use the FBST procedure to examine the mixture parameters in order to choose the parametric distribution that best fits the observed data. But before that, we have performed pairwise comparisons by fitting the lognormal-Weibull, lognormal-gamma, and gamma-Weibull mixture models.

6.2 Results

The measures of evidence provided by HD data in favor of the three models concerning the pairwise comparisons are presented in Table 2. For the comparison between the lognormal and Weibull distributions, the FBST indicates to choose the lognormal model since the e-values $ev(H_{L})=0.874$ and $ev(H_{W})=0.043$. For selecting between the lognormal and gamma distributions, the evidence measures indicate that both models provide good fit to the dataset. Nevertheless, also we would prefer to choose the lognormal model which is the most plausible. The results of the tests for comparison between the gamma and Weibull distributions indicate that the Weibull distribution does not provide reasonable fit to the dataset.

Discrimination based on the LGW mixture model

In order to test simultaneously the three hypotheses, we have applied the the LGW model,

[TABLE]

to the HD data.

The estimates for the parameters of the model (6.1) are presented in Table 3. Here, SD, $2.5\%$ and $97.5\%$ denote the standard deviation, the $2.5$th and the $97.5$th percentiles of the posterior distribution of the LGW parameters, respectively. Both the classical and the Bayesian measures of evidence, presented in Table 4, indicate that neither the gamma and Weibull models should be considered because the null hypotheses $H:p_{2}=0$ e $H:p_{3}=0$ are not rejected. Consequently, among the three models, the lognormal model is the most appropriate for modeling HD data.

Figure 1 displays the survival curves calculated using Bayesian estimates of the lognormal model (Table 5), the LGW mixture model (Table 3) and a procedure called the piecewise exponential estimator (PEXE), introduced by Kim and Proschan (1976), representing the observed data. Unlike the well-known Kaplan-Meier estimator, the PEXE is smooth and continuous estimator of the survival function.

It appears reasonable to disregard both the gamma and the Weibull models; the lognormal model by itself produces a good estimate of survival function.

Note that the preference for the lognormal model is evident in evaluating the LGW mixture model more than in the comparison between the lognormal and gamma distributions, where the evidence measures in favor of both models are very close. It means that the discrimination power provided by LGW model is much higher than the power of the pairwise comparisons. This finding is in agreement with the discussion of Sawyer (1984).

7 Final Remarks

In this paper we considered the FBST for discriminating between survival distributions in the context of linear mixture model. The mixture approach allows us to compare between all alternative models at once by testing the hypotheses on the mixture weights space. The families of survival distributions considered include the lognormal, gamma and Weibull models. In this work, the density functions of the mixture components were reparametrized in terms of the mean $\mu$ and the variance $\sigma^{2}$ of the population so that all models under discrimination share common parameters (Kamary et al., 2014; Pereira and Pereira, 2017).

From the simulation results, we observed that the FBST achieves good performance on identifying the true distribution used to generate the survival times.

The application of the LGW mixture model to the survival data for HD patients allowed us to identify the lognormal distribution as the most appropriate in modeling observed data. Therefore, one can construct a regression model to the HD data considering the lognormal model as the distribution of the response variable.

It would be of interesting to apply the proposed procedure to survival data also considering another censoring mechanisms.

Acknowledgements

The authors are grateful for the support of CNPq, COPPE/UFRJ and IME/USP.