Theoretical Biology and Medical Modelling Open Access Hypertabastic Survival Model

A new two-parameter probability distribution called hypertabastic is introduced to model the survival or time-to-event data. A simulation study was carried out to evaluate the performance of the hypertabastic distribution in comparison with popular distributions. We then demonstrate the application of the hypertabastic survival model by applying it to data from two motivating studies. The first one demonstrates the proportional hazards version of the model by applying it to a data set from multiple myeloma study. The second one demonstrates an accelerated failure time version of the model by applying it to data from a randomized study of glioma patients who underwent radiotherapy treatment with and without radiosensitizer misonidazole. Based on the results from the simulation study and two applications, the proposed model shows to be a flexible and promising alternative to practitioners in this field.


Introduction
Time to event models, commonly known as survival or reliability models, have been studied and applied in a variety of scientific disciplines such as medicine, engineering and business. The Hosmer and Lemeshow [1], Lee and Wang [2], Kleinbaum and Klein [3], and Collet [4] books give a detailed overview of survival data modeling techniques. Non-parametric and semi-parametric survival models such as the Cox regression analysis have been the most widely used models in the analysis of time to event survival data [5]. On the other hand, if the assumption for parametric probability distribution is met for the data set under consideration, it will result in more efficient and easier to interpret estimates than non-parametric or semi parametric models. A comprehensive review was given by Efron [6] and Lee and Go [7].
Parametric hazard functions can enable clinicians and researchers to model various disease scenarios, assess disease prognosis and progression, give valuable insights on the pattern of failure, and understand the pathogenesis of a chronic disease and how they are affected by different treatment effects [8][9][10]. Estimation of hazard function is also useful in the analysis of change-point hazard rate models. It helps policy makers with cost effective health care policy decisions [11]. Lundin et al. [12] estimated the survival probabilities in breast cancer patients and concluded that parametric survival estimates may be more precise than Kaplan-Meier estimates when there are few patients in a particular stratum. Royston and Parmar [13] modeled the baseline distribution function by restricted cubic regression spline. Kay et al. [14] discussed the use of hazard functions in breast cancer studies. They believe that the hazard function is an important tool in investigat-ing disease curability and can help the clinician to express his ideas regarding disease progression and the biology of treatment effect. Foulkes et al. [15] used parametric modeling to assess the prognostic factors in the recurrence of ischemic strokes. Sama et al. [16] used five parametric models to analyze the survival time data of infections and they found that the best fit could be obtained using parametric models. They also indicated that parametric models can be used to model the duration of malaria infections. Kannan et al. [17] used log-logistic probability distribution to model altitude decompression sickness (DCS) risk and symptom onset time. They concluded that the log-logistic model could provide good estimates of the probability of DCS over time. Nardi and Schemper [18,19] emphasized the role of residuals for the selection of survival models. They argued that when empirical data is sufficient, parametric models provide some insight into the shape of the baseline hazard function.
In the Cox model, the baseline hazard function is regarded as a nuisance parameter, while in parametric models, the hazard function reflects the time course of the process under study. In this paper, we introduce a new two-parameter continuous probability distribution called hypertabastic probability distribution. The hypertabastic hazard function can assume a different variety of shapes. It can be used to analyze biomedical data such as cancer recurrence time. Based on the hypertabastic distribution, we introduce the hypertabastic survival model which includes the hypertabastic proportional hazards model with parametric baseline hazard function, the hypertabastic accelerated failure model and the hypertabastic proportional odds model. The hypertabastic distribution can be used to analyze the accelerated hazards regression model of Chen [20]. It can also be used to monitor the disease progression and regression and provide the clinicians with the time interval(s) on which the disease progresses or regresses and the time interval(s) on which the disease progression or regression speeds up or slows down. This vital information will make it easier for the physicians to take appropriate action regarding their patients.

Hypertabastic distribution
In this section we introduce a new probability distribution which can be used in many scientific disciplines. Such disciplines may include, but are not limited to, biomedical, engineering, and business fields. Definition 2.1 (Hypertabastic Distribution) We say a continuous random variable T has a hypertabastic distribution if its cumulative distribution function is The parameters α and β are both positive and Sech [•] and Coth [•] are hyperbolic secant and hyperbolic cotangent respectively. We often read as "T is hypertabastically distributed with parameters α and β "and write it as H (α, β).
The probability density function of T is given by

Hypertabastic survival and hazard Functions
Definition 3.1 (The Hypertabastic Survival Function) Let T be a continuous random variable representing the waiting time until the occurrence of an event. Then the hypertabastic survival function is defined as where P (T > t) is the probability that waiting time exceeds t.
Definition 3.2 (The Hypertabastic Hazard Function) The hypertabastic hazard function h(t) which represents the instantaneous failure rate at time t given survival up to time t is defined as The cumulative hazard function H(t) is defined as The hazard function is a conditional failure rate which gives the instantaneous potential for failing at time t per unit time for an individual surviving to time t.
The most commonly used baseline hazard functions are Weibull, log-normal and log-logistic. The Weibull baseline hazard function has monotone increasing, monotone decreasing and constant forms. The log-normal baseline hazard function is non-monotonic. It follows a behavior that increases to a maximum, and then decreases (unimodal shaped). The log-logistic baseline hazard function can be monotone decreasing or have unimodal shape. These different hazard shapes may reveal different biological mechanisms of disease progression and regression and 1. The hypertabastic baseline hazard function is monotone decreasing if 0 <β ≤ 0.25 ( Figure 1a). Clark et al. [21] analyzed data for 825 patients diagnosed with primary epithelial ovarian carcinoma and observed that the hazard rate was initially high after the diagnosis and gradually decreased afterwards.
2. The hypertabastic baseline hazard function first increases with time until it reaches its maximum and then decreases (unimodal-shaped) if 0.25 <β < 1 (Figure 1b). Demicheli at al. [22] concluded that the hazard rate for node-positive post-menopausal women was unimodalshaped. Schulman et al. [23] studied the influence of donor and recipient HLA locus mismatching on the development of obliterative bronchiolitis (OB) after lung transplantation and estimated the risk of OB after the transplant. Their estimated hazard function of developing OB indicated a unimodal-shaped curve as well.

The hypertabastic baseline hazard initially increases
with time, then it reaches its horizontal asymptote α provided that β = 1 (Figure 1c). Weitz and Fraser [24] concluded that hazard rate plateaus are explained as a generic consequence of considering death in terms of first passage time for processes undergoing a random walk with drift. They analyzed the hazard rate plateau in populations of fruit flies, yeast and other organisms.
4. The hypertabastic baseline hazard function is increasing with upward concavity until it reaches its inflection point and then it continues to increase with downward concavity thereafter if 1 <β < 2 ( Figure 1d). Such a hazard can be used in many applied sciences when failure rate increases with respect to increase in time. For instance, it can be used to model leukemia survival times for patients not responding to treatment [3].
5. The hypertabastic baseline hazard function is increasing with upward concavity for a while and then it becomes a linear function with slope α if β = 2. (Figure 1e).

The hypertabastic baseline hazard function is increasing
with upward concavity if β > 2. (Figure 1f).   If for instance the goal of the drug treatment study is to slow down the disease progression and to understand the time course and management of disease, the above mentioned definitions may be very useful. See for example, the paper by Chan and Holford [25] who studied the rate of deterioration of degenerative diseases over time.

Hypertabastic proportional hazards model
When the effect of risk factors is to change the baseline hazard function by a proportionate amount at all times t, we call the model a proportional hazards model. Suppose vector X is a p-dimensional vector of covariates and assume g(X|θ) is a non negative function of X satisfying the condition that g(0|θ) = 1. Let h 0 (t) which has been defined in section 3 as h(t), be called baseline hazard function. This is the hazard function when there are no covariates in the model. The hypertabastic proportional hazard model assumes a hazard function h(t|X, θ) of the where θ is a vector of unknown parameters. The hypertabastic survival function S(t|X, θ) for the proportional hazards model is defined as  where the baseline survival function S 0 (t) is defined as S(t) in section 3. The hypertabastic probability density function for the proportional hazard model is denoted by f(t|X, θ) and is equal to where the baseline probability density function f 0 (t) is defined as f(t) in section 2. Depending on the type of censoring, the maximum likelihood function technique along with an appropriate log-likelihood function may be used to estimate the unknown parameters in this model. Most data used in survival analysis have only right censoring, therefore we will focus on right censoring. Consider a sample of right censored survival time's data of n individuals t 1 , t 2 , ..., t n with associated p-dimensional covariate vectors X 1 , X 2 , ..., X n and an unknown parameter vector θ = (θ 1 , θ 2 , ..., θ p ).
Then, the hypertabastic proportional hazards log-likelihood function can be written as where Some statistical software packages use logarithm of survival time as their survival time variable in their model fittings. If this is the case, then the following alternative formula can be used as the proportional hazards log-likelihood function: The maximum likelihood estimate of (p+2) dimensional parameter vector λ = (α, β, θ 1 , θ 2 , ..., θ p ) is the vector . Asymptotically, is normally distributed with mean vector λ and variance-covariance matrix V. An estimate of V can be obtained by calculating the inverse of the observed information matrix.
To assess the statistical significance of model parameters, one can use the well known statistical tests such as likelihood ratio test, Wald test, or the score test.
The hazard ratio HR(X i , X j ) for individuals i and j with covariate vectors X i and X j is given by is independent of time t. If X i = X and X j = 0, then the hazard ratio becomes HR(X i , 0) = g(X|θ) and if we use , then , where -θ j x j is the elasticity of hazard rate with respect to covariate x j . Symbolically, The elasticity of hazard function with respect to covariate x j is a measure of the responsiveness of the hazard function to change in covariate x j . Intuitively, elasticity of the hazard function with respect to change in covariate x j is percent change in failure rate divided by percent change in covariate x j .

Hypertabastic accelerated failure time model
When the covariates act multiplicatively on the time-scale, the model is called accelerated failure time model [3,4]. The hypertabastic accelerated failure time model assumes a hazard function h(t|X, θ) of the form where h 0 (•) is the baseline hypertabastic hazard function and the hypertabastic survival function for the accelerated failure time model is where S 0 (•) is the baseline hypertabastic survival function. Finally, the hypertabastic probability density function for the accelerated failure time model is where f 0 (•) is the baseline hypertabastic probability density function. The maximum likelihood technique can be used to estimate the parameters of this model and statistical tests similar to section 4 can be used to assess the sig-  where S 0 (t) is the baseline hypertabastic survival function. The hypertabastic survival function for the proportional odds model is given by

Hypertabastic proportional odds model
and the hypertabastic baseline odds function of survival beyond time t is given by The ratio of the odds of survival of individual i relative to individual j is equal to exp [-θ T (X i -X j )]. The hypertabastic proportional odds model can be fitted by maximizing an appropriate likelihood function.

Simulation study
To evaluate the performance of the hypertabastic model we conduct a simulation study in which we compare the overall fit of the proposed model with Weibull, log-logistic and log-normal models. Since all distributions under our consideration have exactly two parameters, we will use the negative of the log-likelihood as a measure of goodness-of-fit. This measure would result in the same conclusion as the Akaike's Information Criterion (AIC) [26]. Thus the smallest value suggests a better fit. We conduct 1000 simulations with sample size of 200 and random censoring of approximately 40%. We sample timeto-event from 11 different parameter combinations of two parameter Weibull, log-normal and gamma distributions for a total of 33 combinations or 33000 simulations. We fit four mentioned models and average the -log likelihood over 1000 runs to determine which model fits the simulated data with the overall most precision on the average. Simulation results are presented in Table 1.
In simulations where we sample from a two-parameter Weibull distribution, obviously, the Weibull model fits the data with highest precision in all instances. Hypertabastic model is a close second, outperforming log-normal and log-logistic models for all combinations of parameters, with log-normal being the worst. Similarly when sampling is from a two-parameter log-normal distribution, the log-normal model outperforms all other models. The hypertabastic model and the log-logistic show similar results, with the log-logistic being slightly better in eight combinations of parameters and the Weibull model performing the worst. Finally, when sampling from a two parameter gamma distribution, the Weibull model fits with the most precision in seven out of eleven combinations. That is because the Weibull distribution and the gamma distribution have hazard functions which are similar in shapes. In another four instances, the hypertabastic model slightly outperforms the Weibull model. The loglogistic model comes in third, however it is close to the hypertabastic and the Weibull for several combinations, while the log-normal does the worst in all instances.

Application
All models presented in the next two examples were fitted using Mathematica 6.0 (Wolfram Research, Inc.) and SASv9 (SAS Institute Inc., Cary, NC). They were fitted using both time and log(time). However, only log(time) results are presented for the multiple myeloma proportional hazards and brain cancer accelerated failure time models since log(time) is commonly used as a default approach in many packages and procedures. Besides the hypertabastic survival model we fit Weibull, log-normal, log-logistic and Cox models. We present all parameter estimates along with their standard errors and compare -2log-likelihood estimates, as well as AIC, to assess the goodness-of-fit of the models. SAS programs used to fit the hypertabastic proportional hazards models with time and log(time) for multiple myeloma data are provided [see Additional file 1].

Analysis of multiple myeloma data
Multiple myeloma is a cancer formed by abnormal white blood cells, called malignant plasma cells. A malignant monoclonal plasma cell is called plasmacytoma. If the disease spreads throughout multiple bone marrow sites in the body, it is called multiple myeloma. This disease weakens a patient's immune system and is usually difficult to cure. To investigate the performance of the hypertabastic proportional hazards model and to compare it with Cox, Weibull, log-logistic and log-normal models, we analyze a cancer data set obtained from a study conducted by Krall et al. [27]. The data contains information on 65 patients with multiple myeloma in which the patients were treated with alkylating agents. This drug is designed to interfere with the cell's DNA and inhibits cancer cell growth. Of these 65 patients, only 17 survived the duration of the study. The data is right-censored and the survival time from the date of diagnosis is measured in months. The covariates have been measured at diagnosis and the covariate list consists of a logarithm of white blood cell count, serum calcium, presence or absence of Bence Jones protein, proteinuria, gender, age, percent myeloid cells in peripheral blood, percent lymphocytes in peripheral blood, logarithm of percent plasma cells in bone marrow, total serum protein, presence or absence of infection, serum globin, logarithm of blood urea nitrogen, fractures, platelets, and hemoglobin. Our first task is to select risk factors that are statistically significant. To accomplish this task we use the hypertabastic log-likelihood function for proportional hazards model and follow the general variable selection strategy outlined by Collet [4]. We also examine the possibility of interactions. The two most significant risk factors that we found are the logarithm of blood urea nitrogen and hemoglobin. Using the stepwise regression, we fit the Cox proportional hazards model. The Kaplan-Meier survival curve for multiple myeloma data is shown in figure 2. The Cox regression did identify the same variables as the most significant prognostic factors. Table 2 gives results for the hypertabastic and the four comparison models. It shows that the -2log-  Kaplan-Meier survival curve for multiple myeloma Figure 2 Kaplan-Meier survival curve for multiple myeloma.  Figures 3a and 3b show the graph of hypertabastic hazard and survival functions for multiple myeloma data. Figure 3a clearly shows that the hypertabastic hazard function is an increasing function of time. By examining figure 3b, we realize that for patients with a Log (BUN) reading of 1.39 and a hemoglobin level of 10.20 there is about a10% chance of survival beyond 65 months. At the above mentioned mean levels for the hemoglobin and logarithm of blood urea nitrogen, the hypertabastic hazard function shows that the failure rate (hazard) reaches its maximum velocity in about 5.15 months. At this point the disease progression has its highest speed (Figure 3c). Figures 4a and 4b show the 3-dimensional graphs for survival of multiple myeloma patients as functions of time and Log(BUN), as well as time and HGB. Other models under consideration had monotone increasing hazard functions (graphs not shown).

Analysis of glioma brain cancer data
Glioma is a cancer of the brain which begins in glial cells. These cells support the neurons. These cells have a very high rate of growth which can quickly destroy the normal cells. The primary types of glioma cancers are astrocytomas, ependymomas and oligodendrogliomas. Shin et al. [28] studied ways to improve the effectiveness of radiation in the treatment of cerebral malignant astrocytoma. Their study focused on the assessment of the effect of multiple daily fractionated radiation therapy with and without misonidazole. They concluded that the addition of misonidazole did not significantly improve the patients' survival. In this section we apply the hypertabastic accelerated failure time technique to model the survival time of a sample of 30 patients from the randomized trials of radiotherapy with and without the radiosensitizer misinidazole. The data was obtained from the Medical Research Council Working Party (MRC) on misonidazole in gliomas. This data is right censored and has been previously analyzed for the selection of variables by MRC [29]. Survival time was measured in days and the longest survival time was 1098 days. We compared radiotherapy treatment of brain cancer patients with radiosensitizer misoni-dazole to radiotherapy without misonidazole. Figure 5 shows a Kaplan-Meier plot of the estimated survival curves for both groups. The log-rank, Wilcoxon, and likelihood ratio tests are all non-significant, suggesting no significant difference between survival curves for the two groups. The Kaplan-Meier estimates of median survival time for radiotherapy with misonidazole group is 258.5 days and for the radiotherapy without misonidazole group is 488 days. The overall median survival time for both groups combined is 361 days. Using the Kaplan- Meier estimates of survival function, a plot of log-cumulative hazard function against the logarithm of the survival time for individuals in two groups indicates that the data is coming from a Weibul distribution. Knowing this information led us to evaluate the performance of the hypertabastic model. First, we fit a hypertabastic accelerated failure time model to analyze the brain cancer data. Then the hypertabastic accelerated failure time model is compared with Weibul, log-logistic, log-normal accelerated failure time models and the Cox proportional hazards model. This model incorporates a binary covariate coded as treatment = 1 for the type of radiotherapy with misonidazole, and as treatment = 0 for radiotherapy without misonidazole. The second covariate is the age of the patient. Thus, this model contains two covariates-treatment and age. We use the method of maximum likelihood to maximize hypertabastic accelerated failure time log-likelihood function for right censored data discussed in section 5.  Kaplan-Meier survival curves for glioma brain cancer Figure 5 Kaplan-Meier survival curves for glioma brain cancer.

Survival Distribution Function
Days Without Misonidazole With Misonidazole

Discussion and conclusion
In this paper we have introduced a new survival model called hypertabastic survival model. The overall results of our simulation indicate that the hypertabastic model performs best when the data is generated from the Weibull distribution. When the data comes from the gamma distribution, both the Weibull and hypertabastic distributions perform well. In the case when we generate survival data from the log-normal distribution, the log-logistic and hypertabastic distributions perform well but the gamma and Weibull perform poorly. In the last case, the Weibull distribution performance is the worst. We believe that the hypertabastic survival model is to some extent robust with respect to variations in the distribution. We believe that the hypertabastic hazard function can play a role in modeling failure rates in medical, biological and engineering fields. The hypertabastic model can also assist physicians and clinicians with their treatment planning through its ability to predict patients' outcome as well as the risk of disease recurrence. Gilbert et al. [30] argued that the pattern of instantaneous risk over time is more interesting than the cumulative risk. For instance, in a cancer study where recurrence time of malignant tumor is of interest, a bimodal hazard curve may represent the elevated inci-a) Hypertabastic 3D survival curve with variables time and age for group using radiotherapy with radiosentisizer misoni-dazole for glioma brain cancer; b) Hypertabastic 3D survival curve with variables time and age for group using radiother-apy without radiosentisizer misonidazole for glioma brain cancer Figure 7 a) Hypertabastic 3D survival curve with variables time and age for group using radiotherapy with radiosentisizer misonidazole for glioma brain cancer; b) Hypertabastic 3D survival curve with variables time and age for group using radiotherapy without radiosentisizer misonidazole for glioma brain cancer.
a) Hypertabastic hazard curves for glioma brain cancer; b) Hypertabastic survival curves for glioma brain cancer; c) Hypertabastic velocity of hazard curves for glioma brain can-cer Figure 6 a) Hypertabastic hazard curves for glioma brain cancer; b) Hypertabastic survival curves for glioma brain cancer; c) Hypertabastic velocity of hazard curves for glioma brain cancer.
dences of early and late recurrences and the magnitude of the hazard rates at the peaks may reveal the intensity of the failure rates. Therefore we recommend that clinicians, practitioners and analysts consider using this model along with other models, and compare it to the models they ordinarily use before making any decision as to which model would provide the best fit and prediction.