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.

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–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 investigating 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.

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 variableT has a hypertabastic distribution if its cumulative distribution function is

The parameters α and β are both positive andSech [•] andCoth [•] are hyperbolic secant and hyperbolic cotangent respectively. We often read as "T is hypertabastically distributed with parametersα andβ " and write it asH (α,β). The probability density function ofT is given by

whereW(t) =α(1 -t^βCoth(t^β))/β.

Definition 3.1 (The Hypertabastic Survival Function) LetT be a continuous random variable representing the waiting time until the occurrence of an event. Then the hypertabastic survival function is defined asS(t) =P(T >t) =Sech [α(1 -t^βCoth(t^β))/β]

whereP (T >t) is the probability that waiting time exceedst.

Definition 3.2 (The Hypertabastic Hazard Function) The hypertabastic hazard functionh(t) which represents the instantaneous failure rate at timet given survival up to timet is defined ash(t) =α(t^2β-1Csch(t^β)² -t^β-1Coth(t^β))Tanh [W(t)]

The cumulative hazard functionH(t) is defined asH(t) = - ln(Sech [W(t)]).

The hazard function is a conditional failure rate which gives the instantaneous potential for failing at timet per unit time for an individual surviving to timet.

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 can provide helpful medical information. The hypertabastic baseline hazard function shapes are given as follows along with illustrative examples from published literature that could fit a pattern of that particular behavior where applicable:

1. The hypertabastic baseline hazard function is monotone decreasing if 0 <β ≤ 0.25 (Figure1a). 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.

a) Hypertabastic hazard curve for 0 <β ≤ 0.25; b) Hypertabastic hazard curve for 0.25 <β < 1; c) Hypertabastic hazard curve forβ = 1; d) Hypertabastic hazard curve for 1 ≤β < 2; e) Hypertabastic hazard curve forβ = 2; f) Hypertabastic hazard curve forβ > 2.

Full size image

2. The hypertabastic baseline hazard function first increases with time until it reaches its maximum and then decreases (unimodal-shaped) if 0.25 <β < 1 (Figure1b). Demicheli at al. [22] concluded that the hazard rate for node-positive post-menopausal women was unimodal-shaped. 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.

3. The hypertabastic baseline hazard initially increases with time, then it reaches its horizontal asymptoteα provided thatβ = 1 (Figure1c). 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 (Figure1d). 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. (Figure1e).

6. The hypertabastic baseline hazard function is increasing with upward concavity ifβ > 2. (Figure1f).

Definition 3.3 (Disease progression and regression) Let the hazard function be defined on intervalI andh'(t) be its derivative.

1. Disease progresses on the time intervalI ifh'(t) > 0 for allt inI.

2. Disease regresses on the time intervalI ifh'(t) < 0 for allt inI.

3. Disease neither progresses nor regresses on the interval I ifh'(t) = 0 for allt inI.

Definition 3.4 (Speed of progression and regression) Leth(t) be a hazard function defined on time intervalI and leth'(t) andh"(t) be first and second derivative ofh(t) respectively.

1. Ifh'(t) > 0 andh"(t) > 0 for allt inI, then the disease progression speeds up on the time intervalI.

2. Ifh'(t) > 0 andh"(t) < 0 for allt inI, then the disease progression slows down on the time intervalI.

3. Ifh'(t) < 0 andh"(t) > 0 for allt inI, then the disease regression slows down on the time intervalI.

4. Ifh'(t) < 0 andh"(t) < 0 for allt inI, then the disease regression speeds up on the time intervalI.

Definitions 3.3 and 3.4 can assist clinicians, researches and pharmacologists to monitor disease status over time. 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.

When the effect of risk factors is to change the baseline hazard function by a proportionate amount at all timest, we call the model a proportional hazards model. Suppose vector X is a p-dimensional vector of covariates and assumeg(X|θ) is a non negative function of X satisfying the condition thatg(0|θ) = 1. Leth₀(t) which has been defined in section 3 ash(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 functionh(t|X,θ) of the formh(t|X,θ) =h₀(t)g(X|θ)

whereθ is a vector of unknown parameters. The hypertabastic survival functionS(t|X,θ) for the proportional hazards model is defined asS(t|X,θ) = [S₀(t)]^g(X|θ)

where the baseline survival functionS₀(t) is defined asS(t) in section 3. The hypertabastic probability density function for the proportional hazard model is denoted byf(t|X,θ) and is equal tof(t|X,θ) =f₀(t) [S₀(t)]^g(X|θ)-1g(X|θ)

where the baseline probability density functionf₀(t) is defined asf(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 individualst₁,t₂, ...,t_nwith associated p-dimensional covariate vectorsX₁,X₂, ...,X_nand an unknown parameter vectorθ = (θ₁,θ₂, ...,θ_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λ = (α,β,θ₁,θ₂, ...,θ_p) is the vector$\widehat{\lambda }=(\widehat{\alpha },\widehat{\beta },{\widehat{\theta }}_1,{\widehat{\theta }}_2,\mathrm{...},\widehat{\theta }p)$. Asymptotically,$\widehat{\lambda }$ 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 ratioHR(X_i,X_j) for individualsi andj with covariate vectorsX_iandX_jis given byHR(X_i,X_j) =g(X_i|θ)/g(X_j|θ)

The most common form ofg(X|θ) is exp [-θ^TX] [1]. If we use$g\left(X|\theta \right)=\exp \left[-{\theta }^{T}X\right]=\exp \left(-{\displaystyle \sum _{k=1}^p{\theta }_k{x}_k}\right)$, then the hazard ratioHR(X_i|X_j) = exp [-θ^T(X_i-X_j)] whereHR(X_i,X_j) is independent of timet. IfX_i=X andX_j= 0, then the hazard ratio becomesHR(X_i, 0) =g(X|θ) and if we use$g\left(X|\theta \right)=\exp \left(-{\displaystyle \sum _{k=1}^p{\theta }_k{x}_k}\right)$, then$\ln \left[g\left(X|\theta \right)\right]=-{\displaystyle \sum _{k=1}^p{\theta }_k{x}_k}$, where -θ_jx_jis the elasticity of hazard rate with respect to covariatex_j. Symbolically,∂h(t|X,θ)/∂x_j]·[x_j/h(t|X,θ)] = -θ_jx_j.

The elasticity of hazard function with respect to covariatex_jis a measure of the responsiveness of the hazard function to change in covariatex_j. Intuitively, elasticity of the hazard function with respect to change in covariatex_jis percent change in failure rate divided by percent change in covariatex_j.

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 functionh(t|X,θ) of the formh(t|X,θ) =h₀(tg|X,θ))g(X,θ)

whereh₀(•) is the baseline hypertabastic hazard function and the hypertabastic survival function for the accelerated failure time model isS(t|X,θ) =S₀(tg(X|θ))

whereS₀(•) is the baseline hypertabastic survival function. Finally, the hypertabastic probability density function for the accelerated failure time model isf(t|X,θ) =f₀(tg(X|θ))g(X|θ)

wheref₀(•) 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 significance of model covariates. For the right censored data, the hypertabastic accelerated failure time model would have a log-likelihood function of the form

The maximum likelihood estimate of (p+2) dimensional parameter vector (α,β,θ₁,θ₂, ...,θ_p) and testing of hypothesis regarding model parameters are similar to methods of section 4.

If the effect of risk factor is to change the odds of survival beyond timet by a proportionate amount, then the model is called proportional odds model. The odds of surviving beyond timet are expressed asS(t|X,θ)/(1 -S(t|X,θ)) =g(X|θ)S₀(t)/(1 -S₀(t))

whereS₀(t) is the baseline hypertabastic survival function. The hypertabastic survival function for the proportional odds model is given byS(t|X,θ) = 1/[1+((1 -S₀(t))/S₀(t))g(X|θ)]

and the hypertabastic baseline odds function of survival beyond timet is given byS₀(t)/(1 -S₀(t)) =Sech [W(t)]/(1 -Sech [W(t)])

The ratio of the odds of survival of individuali relative to individualj is equal to exp [-θ^T(X_i-X_j)]. The hypertabastic proportional odds model can be fitted by maximizing an appropriate likelihood function.

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 time-to-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 Table1.

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 log-logistic 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.

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 bothtime andlog(time). However, onlylog(time) results are presented for the multiple myeloma proportional hazards and brain cancer accelerated failure time models sincelog(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 withtime andlog(time) for multiple myeloma data are provided [see Additional file1].

8.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 figure2. The Cox regression did identify the same variables as the most significant prognostic factors. Table2 gives results for the hypertabastic and the four comparison models. It shows that the -2log-likelihood and AIC statistics are lowest for the hypertabastic model when fitted to the multiple myeloma data. This indicates that the hypertabastic proportional hazards model fits the multiple myeloma data best. The most significant single variable identified by the hypertabastic model is the logarithm of blood urea nitrogen with an estimated chi-square value of 11.11 (p = 0.0014). The second significant variable is hemoglobin. All other models have identified these two variables as the most significant ones. The mean levels of patients' hemoglobin and the logarithm of blood urea nitrogen are 10.20 and 1.39 respectively. At this level, the median survival time for the hypertabastic, log-normal, log-logistic, Weibull and Cox are 20.04, 19.29, 21.01, 21.92 and 19 months respectively. Figures3a and3b show the graph of hypertabastic hazard and survival functions for multiple myeloma data. Figure3a clearly shows that the hypertabastic hazard function is an increasing function of time. By examining figure3b, 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 (Figure3c). Figures4a and4b 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).

Kaplan-Meier survival curve for multiple myeloma.

Full size image

a) Hypertabastic hazard curve for multiple myeloma; b) Hypertabastic survival curve for multiple myeloma; c) Hypertabastic velocity of hazard curve for multiple myeloma.

Full size image

a) Hypertabastic 3D survival curve with variables time and Log (BUN) for multiple myeloma data; b) Hypertabastic 3D survival curve with variables time and HGB for multiple myeloma data.

Full size image

8.2 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 misonidazole to radiotherapy without misonidazole. Figure5 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. Table3 gives a statistical summary for the glioma data. For instance, the hypertabastic estimated value for the coefficient of the variable radiosensitizer is 0.4387 with a standard error of 0.3437. The Wald and the likelihood ratio statistics associated with this variable are 1.6294 (p = 0.2018) and 1.6254 (p = 0.2023) respectively. Both tests indicate that the effect of individual variable radiosensitizer is non-significant. The estimated accelerator factor is 1.5507 (exponentiated value of parameter 0.4387). This means that after controlling for the age of the patient, the probability of a patient treated with "therapy with radiosensitizer misonidazole" survivingt days equals to the probability of a patient treated with "therapy without radiosensitizer misonidazole" surviving 1.5507t days. For instance, the hypertabastic accelerated failure time model suggests that for 49 year old patients (median age for all patients under study), the probability of a patient treated with " therapy with radiosensitizer misonidazole " surviving 293 days equals to the probability of a patient treated with" therapy without radiosensitizer misonidazole " surviving 454 days. The Wald statistic for the age covariate is 22.18 (p < 0.0001).

Kaplan-Meier survival curves for glioma brain cancer.

Full size image

Using AIC, we came to the conclusion that both Weibull and hypertabastic models fit the data very well. However, the Weibull AIC value was 75.54 which indicates the best fit. For the hypertabastic model, the AIC value was 76.72 which was slightly less than Weibull. The remaining AIC values for log-normal, log-logistic and Cox are 79.15, 78.45, and 107.43 respectively. At the median age of 49, the median survival time for those patients who underwent radiotherapy with radiosensitizer misonidazole using the hypertabastic, Weibul, log-normal, log-logistic, and the Cox models are 289, 316, 270, 277, and 244 days respectively. The corresponding median survival times for the group without the radiosensitizer misonidazole are 449, 453, 423, 456, and 488 days respectively. By examining the Kaplan-Meier survival functions for the two groups we observe the crossing pattern which is a clear indication of violation of proportional hazards assumption. Figures6a and6b show the graphs of hypertabastic hazard and survival functions for both treatment groups at the age level of 49. The hazards for both groups are increasing function of time. For those patients who received radiotherapy treatment with radiosensitizer misonidazole, hazard function reached its maximum velocity in about 200 days. For those who did receive radiotherapy without misonidazole, hazard function reached its maximum velocity in about 311 days. These are the points in time where the speeds of failure rates (hazards) are highest (Figure6c). Figures7a and7b represent 3-dimensional graphs of survival by time and age for each treatment group separately.

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.

Full size image

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.

Full size image

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 incidences 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.

Authors and Affiliations

1. Department of Mathematical Sciences, Cameron University, Lawton, OK, USA

   Mohammad A Tabatabai

2. Department of Biostatistics, University of Arkansas for Medical Sciences, Little Rock, AR, USA

   Zoran Bursac & David K Williams

3. Department of Biostatistics, University of North Texas Health Science Center, Ft.Worth, TX, USA

   Karan P Singh

Corresponding author

Correspondence toZoran Bursac.

Competing interests

The author(s) declare that they have no competing interests.

Authors' contributions

The survival model equations were developed by MAT and ZB based on their previous work on hyperbolastic growth and survival models. Applications of the models to data were done by MAT, ZB, DKW and KPS. Writing of the manuscript was done by MAT and ZB, and final edits by DKW and KPS.

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions