Bayesian index to compare assumed survival functions of two hazards following exponential distributions

In clinical studies, the period before an event occurs is frequently considered as an outcome. The general term used for this analytical method in clinical studies is “survival time analysis”. The method employs Cox regression to compare hazards between groups. In this paper, we propose an index defined in the Bayesian framework to compare hazards. This index compares the sizes of the hazards directly. Further, we propose two methods to calculate the sizes of the hazards. One is an approximate method, by which we can easily calculate the size, and the other is an exact method. In addition, we apply the proposed index to two actual clinical studies to assess its usability. Subjects: Mathematics & Statistics; Medical Statistics; Medical Statistics & Computing; Medicine; Medicine, Dentistry, Nursing & Allied Health; Science; Statistics; Statistics & Probability


PUBLIC INTEREST STATEMENT
In clinical studies, the period before an event occurs is frequently considered as an outcome. The general term used for this analytical method in clinical studies is "survival time analysis". The method employs Cox regression to compare hazards between groups. In this paper, we propose an index defined in the Bayesian framework to compare hazards. We apply the proposed index to two actual clinical studies to assess its usability.
It is not an exaggeration to say that in the frequentist framework, the basic methodology of survival time analysis has almost been settled (Colett, 2003;Hosmer, Lemeshow, & May, 2008). Besides, recently the method of direct comparison of hazards was proposed (Heller & Mo, 2015).
In the Bayesian framework, on the other hand, a methodology has not yet been established. In this paper, Bayes introduced the idea of inverse probability, which, following the work of Laplace and others, came to be known as Bayes' Theorem. The Bayesian approach to statistical inference is useful for several reasons. Perhaps foremost among them is the formal manner in which Bayesian methods incorporate prior information, whether it be expert opinion or historical data. The sequential learning nature of the method, with prior knowledge being updated to yield posterior distributions, forms the basis for adaptive methods, in which models are continually updated with new information. Other advantages of the method include the ease of predictive model building and the avoidance of large-sample approximations in some cases. Furthermore, Bayesian analysis provides straightforward probabilistic statements about parameters of interest and inferential questions. The Bayesian influence can be found in a broad array of fields such as statistics, medicine, epidemiology, computer science, economics, engineering, physics, and philosophy of science. A number of papers on biomedical applications have suggested the use of Bayesian methods in healthcare evaluation (Spiegelhalter, 2004;Zhao, Feng, Chen, & Taylor, 2015), the pharmaceutical industry (Grieve, 2007;Racine, Grieve, Fluhler, & Smith, 1986), and clinical trials (Berry, 2006;Cornfield, 1976;Lee & Chu, 2012;Spiegelhalter, Freedman, & Parmar, 1993, 1994. As mentioned above, Bayesian statistics are employed in clinical trials to evaluate outcomes, but no one has so far proposed an index for the easy interpretation of results. In Bayesian model selection, the most widely used tools are Bayes factors (BFs) (Kass & Raftery, 1995). However, these are unscaled when improper priors are used. In order to overcome this problem and approximate the underlying BFs, fractional and intrinsic BFs (Berger & Pericchi, 1996;O'Hagan, 1995) have been proposed in the literature. Thus, BFs are useful in Bayesian model selection. However, we cannot directly compare parameters of posterior distributions using BFs. Moreover, it is often difficult to interpret BFs for several reasons; one reason is that the coefficients of BFs do not have intuitive meaning. Zaslavsky (2010Zaslavsky ( , 2013 and Kawasaki and Miyaoka (2012b) were the first to propose a new type of index, called the "Bayesian index", for the comparison of binomial proportions and to apply it to clinical trial results. These indexes can compare parameters directly, and since the results are evaluated as a probability, it is easy to understand them intuitively. In this paper, we propose the index η = P(h 1 < h 2 |T 1 , T 2 ), which is developed from the index proposed by Miyaoka (2012a, 2012b), for comparison of hazards, where h i indicates the hazard of group i, with the exponential distribution assumed for the survival function. We further assume that the prior distribution is a gamma distribution with parameters α i and β i . T i denotes the total observation duration of group i.
This index indicates the probability that the hazard of group 1 is less than that of group 2.
The remainder of the paper is organized as follows. We propose a method for the exact calculation of the index η as well as a method for an approximate calculation, with which we can easily calculate the probability of the index, in Section 2. We compare the approximate probability with the exact probability and present the results of actual clinical trials to show the utility of η in Section 3. Finally, we conclude with a brief summary in Section 4.

Notation
Let C ij and Y ij be independent. For patient number j in group i, let Y ij be the random variable representing the time to event occurrence, C ij be the random variable representing the time to censoring (where i = 1, 2 and j = 1, 2, …, n i ), and T ij = min(Y ij ,C ij ) be the variable representing the observation duration. Additionally, to judge the executing or censoring event for patient number j in group i, we define the event index as where the realized value of T ij is t ij . Furthermore, let the hazard function of patient number j in group i be constant as follows The probability density function and survival time function can be denoted, respectively, as

Posterior distribution
We give the likelihood function for group i as Assuming gamma distributions with parameters α i and β i as prior distributions, and i = 1, 2, we have

Exact expression
Theorem The exact expression for the posterior hazard of group 1 to be smaller than that of group 2 is where 2 F 1 a 1 , a 2 ; b; x denotes the hypergeometric series, I x a, b denotes the regularized incomplete beta function and Beta(a, b) denotes the beta function.
Proof The posterior probability is given by Hence, its joint probability density function is From the definition, η can be given as From Theorem 1 in Kawasaki and Miyaoka (2012a), Furthermore, from Formulas 26.5.23 and 26.5.2 of Abramowitz and Stegun (2010), i.e. and I z a, b = 1 − I 1−z b, a , η can be given as We can calculate the exact probability by Equation (1) and using expressions containing various hypergeometric series. In this paper, we use one such expression. Kawasaki and Miyaoka (2012b) give a more detailed proof in a similar problem. For the SAS programming code used for calculating the probability of η, see Appendix 1.

Approximate expression
We can easily formulate an approximate expression using the expected value of the posterior distribution and the variance. The expected value and the variance of the difference between the posterior distributions are, respectively, where h i,post denote the posterior hazard, and i = 1, 2. Therefore, we can rewrite Equation (2) as where Φ(·) is the cumulative distribution function of the standard normal distribution. Thus, η can easily be calculated if we know the cumulative distribution function of the standard normal distribution.

Simulation
Our aim of simulation in Section 3.1 is to verify the condition under which the differences of the probability, indicated by the approximated method and that by the exact method, become greater. Hence, we compare the exact method and the approximation method using simulated data in this section. The simulation method is as follows. Assume the data occur randomly from the exponential distribution when the hazard ratio is between 0.40 and 1.50 in steps of 0.10, and that they occur randomly with binomial proportions when the censoring ratio is 10 or 40%. Further, let the sample sizes for the groups be either balanced or unbalanced. For the balanced cases, we set n 1 = n 2 = 10, n 1 = n 2 = 20, n 1 = n 2 = 50, or n 1 = n 2 = 100, and for the unbalanced cases, we set n 1 = 15 and n 2 = 5, or n 1 = 75 and n 2 = 25. In this simulation, we use the non-informative prior, that is, α i = β i = 0.001 and i = 1, 2.
At the beginning, we verify the difference between the probability by the approximated method and that by the exact method with small sample size. Figure 1 compares the results under the censoring ratios 10 and 40% when n 1 = n 2 = 10 and n 1 = n 2 = 20. The box plots show the difference between the exact probability and the approximate probability, with the vertical axis representing the difference between the exact probability and the approximate probability, and the horizontal axis representing the hazard ratio. We make the following observations from Figure 1: • When the censoring ratio is 10%, the difference between the exact probability and the approximate probability is about 4% at maximum.
• The difference is about 12% at maximum when the censoring ratio is 40% and the hazard ratio is 0.40.
• The difference shows a variance of around 0% when the hazard ratio is 1.
Next, we verify that with large sample size. Figure 2 compares the results under the censoring ratios 10 and 40% when n 1 = n 2 = 50 and n 1 = n 2 = 100. Figure 2 shows the following: Notes: The vertical axis shows the difference between the exact probability and the approximate probability, and the horizontal axis indicates the hazard ratio. HR n 1 =n 2 =10, censoring ratio=10% n 1 =n 2 =10, censoring ratio=40% n 1 =n 2 =20, censoring ratio = 10% n 1 =n 2 =20, censoring ratio = 40% • The results are not affected by censoring, and the difference between the exact probability and the approximate probability goes up to about 2.5%.
• The interquartile range is very narrow.
• The difference is even larger relative to the case where both groups have the same sample size in total.

Chronic active hepatitis cases
In this section, we apply the proposed index to clinical trial data used in Kirk, Jain, Pocock, Thomas, and Sherlock (1980). In this trial, 44 chronic active hepatitis cases were assigned randomly to either the prednisolone group (treatment group) or the control group with no treatment (control group). From the treatment group's 22 cases, 11 were observed and the remaining 11 were censored, and from the control group's 22 cases, 16 were observed and the remaining 6 were censored. The most interesting outcome in the trial was the survival time of patients after the trial commenced.
The p-value was 0.0309 by the log-rank test and 0.0105 by the Wilcoxon test. This indicates a longer survival time for the treatment group at the significance level of 5%. The proposed index indicates 98.35% as the approximate probability and 99.01% as the exact probability, with the hazard of the treatment group lower than that of the control group. We assume the prior distribution to be non-informative for both groups.

Lung cancer data analysis
In this section, we apply the proposed index to clinical trial data used in Chen et al. (2007).
We obtained the non-small-cell lung cancer data used in Chen et al. (2007) from http://www.ncbi. nlm.nih.gov/projects/geo/ with accession number GSE4882. The data contained 672 gene profiles of 125 lung cancer patients. From among them, 38 patients died, and the remainder were censored.
The median follow-up of 101 patients was 20 months. The patients with a high-risk gene signature showed a shorter median overall survival than those with a low-risk gene signature. The p-value was 0.001 by the log-rank test; the proposed index gave 99.99% probability as calculated approximately as well as exactly. A high-risk gene signature was associated with a median relapse-free survival period of 13 months, whereas a low-risk gene signature was associated with a median relapse-free survival period of 29 months. Here, the p-value was 0.002 by the log-rank test, and the proposed index gave 99.53% as the approximate probability and 99.72% as the exact probability.

Discussion and conclusions
In this paper, we have proposed an index to compare the hazards defined in the Bayesian framework. Further, we have presented two probability functions: the approximation function, which provides for easy calculation, and the exact function.
Our simulation experiment shows a difference between the exact probability and the approximate probability. Specifically, when we censor more, the difference becomes larger. By contrast, we found that the difference becomes smaller with larger sample sizes; the effect of censoring too becomes less. With unbalanced sample sizes for the groups, the difference becomes larger. Thus, we recommend using the exact probability calculation for cases with small samples, unbalanced sample sizes, or numerous censorings.
We applied the proposed index to actual clinical trial data and confirmed that it shows the difference between groups very clearly. Because, we indicated it in example 1 that the hazard of the treatment group becomes lower than that of the control group with about 99% probability. Thus, as this index provides comparison of the hazards with probability, the result can be indicated within the range of 0-100%. So, the index can be very easy to understand intuitively.
Thus, the index proposed in this paper to compare hazards is easy to understand and intuitive and should be beneficial for use in research studies. http://dx.doi.org/10.1080/23311835.2016.1193927

Appendix 1 SAS programming code
The following SAS macro can be used to calculate η = P(h 1 < h 2 |T 1 , T 2 ).
The parameters are defined as follows: N1 = sample size of group 1; N2 = sample size of group 2; X1 = total survival time of group 1; X2 = total survival time of group 2; e1 = total event count for group 1; e2 = total event count for group 2; alpha1 and beta1 = hyperparameters for group 1; alpha2 and beta2 = hyperparameters for group 2.