HTML

ePub

PDF

1 Introduction

Statistical lifetime methods are widely used in studies considering the analysis of data sets in the health area, when the variable of interest is related to the time until the occurrence of an event (Klein and Moeschberger, 2005). This event can be, for example, the first recurrence of a disease after treatment, death due to a specific cause or discharge after a hospital admission. While many methods of survival analysis are readily available in statistical computer packages, most of these techniques assume that all individuals are susceptible to the occurrence of the event of interest. However, this assumption can be unrealistic in some specific situations, when there is a subpopulation of individuals that is immune or not susceptible to the occurrence of the event of interest. As an example, in a clinical research, a proportion of patients can respond favourably to the treatment under investigation and be regarded as "cured". A number of statistical tools have been developed to handle such data, especially the cure fraction models (Lambert et al., 2007).

These models include as a special case the mixture model (Farewell, 1982), that explicitly includes a parameter accounting for the cure rate. The mixture model assumes that the probability of the time-to-event to be greater than a specified time t is given by the survival function

where p is a parameter which represents the proportion of "cured patients", regarding the event of interest (0 < p< 1), and S₀(t) is the baseline survival function for the susceptible individuals (Maller and Zhou, 1996). Common choices for S₀(t) are the Gompertz, exponential and Weibull distributions. The probability density function for the lifetime T is given by

where F (t)=1 S(t) and f₀(t) is the baseline probability density function for the susceptible individuals.

Models based on defective distributions are alternatives to the mixture model. A defective distribution is defined as a distribution that is not normalized to one for some values of their parameters, and it can be useful to fit data including both immune and susceptible individuals without explicitly including the parameter p. The use of defective distribution in survival analysis have been considered by a number of authors, including Balka et al. (2011), Cancho and Bolfarine (2001) and Rocha et al. (2014).

In this paper we propose Bayesian and frequentist approaches to a model based on the modified Gompertz distribution, considering survival data in presence of a cure fraction. More recently a number of new families of defective distributions have been proposed in the literature, such as those based on the Kumaraswamy Rocha et al. (2015) and the Marshall–Olkin families Rocha et al. (2017). However, in the present study, the modified Gompertz distribution was preferred because of its computational simplicity and the ease of interpretation of their parameters. Our methodology is motivated by a real data set that arises from the cervical carcinoma study by Brenna et al. (2004) and it extends the model described by Rocha et al. (2014) to the situation in which covariates are present. The computational codes used in this article are provided in the Appendix.

2 Methods

2.1 The modified Gompertz model

The modified Gompertz hazard function was first used by Cantor and Shuster (1992) for estimation of cure rates from paediatric clinical trials and after extended by Gieser et al. (1998) to include covariate effects. This model assumes a survival function given by

where t> 0, α> 0 is a shape parameter, β > 0 is a scale parameter, and the corresponding cure rate η is given by

where 0 <η < 1. In addition, the corresponding probability density function is given by

and the hazard function is given by

2.2 Maximum likelihood estimation

Considering a random sample (t_i,δ_i) of size n, i = 1,...,n, the contribution of the ith subject for the likelihood function is given by

where δ_iis a censoring indicator variable, that is, δ_i = 1 for an observed lifetime and δ_i = 0 for a censored lifetime.

Assuming the modified Gompertz model, the likelihood function for θ = (α,β) is given by

and the corresponding log-likelihood function is

By deriving the log-likelihood function with respect to α and β, we have the following equations:

and

Setting these expressions equal to zero, we get the corresponding score equations whose numerical solution leads to the maximum likelihood estimators (MLE). Although we cannot obtain explicit expressions for the MLEs for the parameters α and β, they can be estimated numerically using iterative algorithms such as the Newton-Raphson method and its variants (Rocha et al., 2014). R (R Development Core Team, 2009) and SAS (Littell et al., 2006) codes for implementing these procedures are presented in the Appendix.

The asymptotic variances of MLEs are given by the elements of the inverse of the Fisher’s information matrix I(α,β). The second partial derivatives of the maximum likelihood function are given as follows:

and

Therefore, the expected Fisher’s information matrix is given by

In practical applications we could use the observed Fisher’s information matrix given the difficulties to get the expected matrix in the determination of the asymptotical normality of the MLEs used to construct confidence intervals and hypotheses tests for the parameters of the model. Wald-type 95% confidence intervals for the parameters can be thus obtained from the respective estimates of the standard errors.

2.3 Model with covariates

In order to include covariates, the parameter α in the likelihood function L(θ) can be replaced by a function α(x_i) such as

where x_i = (1,x_1i,x_2i,...,x_pi) is a vector containing observations on p independent variables and α∗ = (α₀,α₁,...,α_p) is a vector of unknown parameters. Analogously, the parameter β in L(θ) can be replaced by an function β(w_i) such as

where w_i = (1,w_1i,w_2i,...,w_qi) is a vector which may or may not be equal to x_i and β* = (β₀,β₁,...,β_q) is a vector of unknown parameters.

Comparisons between different model formulations can be based on the Akaike’s information criterion (AIC) (Akaike, 1973). The model with the lowest AIC value suggest a better fit to the data.

2.4 Bayesian analysis

The Bayesian method considers that prior distributions of the parameters of the model are used to derive their corresponding posterior densities (Gelman et al., 2013). According to the Bayes theorem, we can write the joint posterior density by combining the joint prior distribution with the likelihood function for α and β. In this case, Markov chain Monte Carlo (MCMC) is a useful method to sample from posterior probability distributions by means of simulation.

For a Bayesian analysis of the model without covariates, we assume inverse gamma (IG) prior distributions for the parameters α and β, considering that these parameters are real and positive numbers. Thus,

and

where a₁, b₁, a₂ and b₂ are known hyperparameters. The mean and variance of a random variable following the inverse gamma distribution with parameters a and b are given by

respectively. Therefore, an improper prior distribution is considered if ≥ 2. An improper prior distribution is not a true probability distribution in that it does not integrate to 1. In Bayesian inference, an improper prior distribution is acceptable, in the sense that the corresponding posterior distribution is proper.

In the case of the model with covariates, consider the following prior distributions for the parameters:

and

where N (c,d) denotes a normal distribution with mean c and variance d, and c_j, d_j, e_kand f_k (j = 0,1,...,p and k = 0,1,...,q) are known hyperparameters.

Posterior summaries of interest will be obtained in the examples section from simulated samples for the joint posterior distribution using standard MCMC procedures, as the Gibbs sampling. We will generate 1,005,000 samples for each parameter of interest. The 5,000 first simulated samples will be discarded as a burn-in period, which is usually used to minimize the effect of the initial values. The posterior summaries of interest will be based on 10,000 samples, taking every 100th sample to have approximately uncorrelated values. We assume prior independence between all model parameters. The Bayes estimates of the parameters will be obtained as the mean of samples drawn from the joint posterior distribution. In addition, 95% credible intervals for the parameters are given by the 0.025th and 0.975th percentiles of the respective posterior distributions. Convergence of the MCMC algorithm will be monitored by usual time series plots for the simulated samples and also using some existing Bayesian convergence methods.

Posterior summaries of interest will be obtained using the OpenBUGS software (Lunn et al., 2000), that only requires the specification of the distribution for the data and the prior distributions for the parameters. See Appendix for details about the OpenBUGS codes used.

The Deviance Information Criterion (DIC) will be used to compare the fit of different model formulations (Spiegelhalter et al., 2014). Smaller values correspond to models that provide better fit to the data.

2.5 A simulation method for a random variable with a modified Gompertz distribution

The algorithm used to simulate a sample of size n from the modified Gompertz distribution with right-censored data follows the steps:

1. Step 1. Fix values of α and β.

2. Step 2. Generate n random samples from

where η is the corresponding cure rate given by

1. Step 3. Consider

where

and

1. Step 4. Generate n random samples from

considering only the finite t´_i.

1. Step 5. Calculate t_i = min(t´_i,u´_i).

2. Step 6. Pairs of values (t₁,δ₁), (t₂,δ₂),...,(t_n,δ_n) are thus obtained, where δ_i = 1 if t_i ≥ u´_iand δ_i= 0 if t_i ≥ u´_i, i = 1,...,n.

This algorithm is similar to that presented by Rocha et al. (2017). An R function for generating random samples based on these steps is presented in the Appendix.

3 Results

3.1 Simulation study

A brief simulation study of the performance of the maximum likelihood method of estimation was carried on based on simulated samples of sizes n = 25, 50, 75, 100, 150 and 200. Each sample was replicated 5000 times. The variance of the cure rate η was estimated using the delta method. The results from the simulation study are presented in Table 1, considering a nominal confidence coefficient of 95 per cent and two sets of arbitrary values for the parameters, given by (α,β) = (0.4,0.2) and (α,β) = (0.3,0.7).

The results in Table 1 show that the coverage probability of the Wald-type confidence intervals for the parameters α and β is quite close to the nominal confidence level of 95%. Furthermore, the coverage probability of the confidence intervals for η approaches 95% as the sample size is increased. In all simulations, the biases and the mean squared errors (MSE) always approached zero as the sample size increased.

3.2 An application to real data

Let us consider the data from a study that included a total of 148 women diagnosed and treated for invasive cervical carcinoma between 1992 and 2002 (Brenna et al., 2004). For the purposes of the present study, it was considered a subsample of 118 women who received the standard treatment recommended by the International Federation of Gynecology and Obstetrics (FIGO). Let us define the disease-free survival (DFS) as the time in complete months from the date of surgery to the first event of disease recurrence. Nearly 48% of the data are censored observations.

Figure 1 is a contour plot of the log of the likelihood function for these data, considering the parametric model based on the modified Gompertz distribution. The plot indicates that the maximum of the log-likelihood function is at (α,β) = (0.04178, 0.03995), as showed in the Table 2. Table 2 shows the maximum likelihood estimates of α, β and the cure rate η, as well as the estimates of their standard errors and confidence intervals. Figure 2 shows plots of the Kaplan-Meier estimates for the survival function and the predicted values obtained from the parametric model against times (months). From Figure 2, it is possible to note that the predicted values obtained from the model based on the modified Gompertz distribution are closest to the empirical values, suggesting that this model is well fitted for the data.

The cure rate is estimated by

and the respective standard error was obtained by the Delta method as a first-order approximation of a Taylor-series expansion (Table 2).

Table 3 shows the Bayesian estimates of α, β and the cure rate η, with their respective 95% credible intervals. In this analysis, it was considered the hyperparameters a₁ = b₁ = a₂ = b₂= 0.001 in order to obtain noninformative priors, that is, α~IG (0.001,0.001) and β~ IG (0.001,0.001). We can note that the Bayesian estimates for the parameters are relatively close to those obtained by the maximum likelihood method (Table 2).

3.3 Model with covariates

In order to illustrate an application of the regression model based on the modified Gompertz distribution, in this analysis we consider the following variables:

• Age of the patients at start of follow-up (x₁), classified as less than 50 years (x₁ = 0) versus greater or equal to 50 years (x₁= 1).

• Clinical stage of the disease at the start of treatment, classified as I, II, or III. This variable is represented in the regression model by using two dummy variables (x₂ and x₃), where x₂ = 0 and x₃ = 0 if stage I, x₂= 1 and x₃ = 0 if stage II, and x₂ = 0 and x₃= 1 if stage III.

Thus, the regression model for the data considers that

and

where α₁₂, α₁₃, β₁₂ and β₁₃ are interaction terms. Let us consider the following model formulations:

• Model 1: Model 1 considers only the effect of age (x₁). In this case, α₂, α₃, α₁₂, α₁₃, β₂, β₃, β₁₂ and β₁₃ are considered equal to zero.

• Model 2: Model 2 considers only the effect of the stage of the disease (dummy variables x₂ and x₃). In this case, α₁, β₁ and the interaction terms α₁₂, α₁₃, β₁₂ and β₁₃, are considered equal to zero.

• Model 3: Model 3 considers the effects of age (x₁) and the stage of the disease (x₂ and x₃) but it does not includes the interaction terms α_12, α₁₃, β₁₂ and β₁₃.

• Model 4: Model 4 considers the effects of age (x₁) and the stage of the disease (x₂ and x₃) and it includes the interaction terms α₁₂, α₁₃, β₁₂ and β₁₃.

For all regression coefficients, we assumed normal prior distributions with mean 0 and large variance. Both Bayesian and maximum likelihood estimates are shown in Table 4. The fit of the Models 1 and 2 also considered the estimation of the cure fraction (η) for each level of the respective independent variable. In the frequentist estimation, standard errors for the cure fraction were calculated using the Delta method. As an illustration, Figure 3 shows the survival curves estimated by the Kaplan-Meier method and by the Model 2 using the maximum likelihood approach.

Table 4 shows that the maximum likelihood and the Bayesian estimates are fairly close to each other. Model 1 suggests that the age do not have a significant effect on the disease-free survival, since the confidence and credible intervals for the parameters α₁ and β₁include the value 0. Model 2 suggests a significant effect of the stage of the disease on the disease-free survival, given that the confidence and credible intervals for the shape parameters α₂ and α₃do not include the zero value (Figure 3).

Although the Model 1 does not evidences a significant effect of the age on the disease-free survival, Model 4 suggest that the interaction between age and the clinical stages is important to understand the role of the age on the time to progression of the disease. The confidence and credible intervals for the interaction terms β₁₂and β₁₃do not include the zero value, suggesting a significant effect of interaction. In fact, the survival curves showed in the Figure 4 help us to understand this effect. The parametric curves in Figure 4 were obtained based on the Model 4, and shown that the age has a strong effect on the disease-free survival of patients in clinical stage II, but the age is not an important factor influencing the survival of patients in clinical stage III given that the survival curves showed in the panel (c) are close to each other.

Models 2, 3 and 4 have similar AIC and DIC values (Table 4), and these values are lower than the respective values calculated for the Model 1. However, under a clinical point of view, the Model 4 seems to be the most appropriate to the data, given that the Figure 3 shows the importance of considering the interaction terms for the interpretation of the pattern of association between age, clinical stages and the time to progression of the disease.

4 Conclusions

Time-to-an-event data including a proportion of individuals who are immune to the event of interest are common, especially in medical studies. Basic tools for survival analysis usually consider that the survival function S(t) tends to zero as the time t tends to infinity, and this assumption is unrealistic if immune individuals are present. The use of models based on defective distributions is a suitable way to analyse data in this situation. In this way, the parametric model based on the modified Gompertz distribution allows for the estimation of the cure fraction and allows the insertion of a vector of covariates. Moreover, the model can be easily implemented in computational programs as SAS, R and OpenBUGS, as showed in the Appendix. To illustrate the use of the model we used a real data set from a cervical cancer study. We can note that the model satisfactorily fitted the data.

Acknowledgements

This research was supported by grants from CNPq (Brazil). The authors are grateful to the anonymous reviewers for their valuable comments and suggestions.

Referências

Akaike, H. (1973). Information theory and an extension of the maximum likelihood principle. In: Proceedings of the 2nd International Symposium on Information Theory, pp. 267–281.

Balka, J., Desmond, A. F., McNicholas, P. D. (2011). Bayesian and likelihood inference for cure rates based on defective inverse Gaussian regression models. Journal of Applied Statistics, 38(1), 127–144.

Brenna, S. M., Silva, I. D., Zeferino, L. C., Pereira, J. S., Martinez, E. Z., Syrjänen, K. J. (2004). Prognostic value of P53 codon 72 polymorphism in invasive cervical cancer in Brazil. Gynecologic Oncology, 93(2), 374–380.

Cancho, V. G., Bolfarine, H. (2001). Modeling the presence of immunes by using the exponentiated-Weibull model. Journal of Applied Statistics, 28(6), 659–671.

Cantor, A. B., Shuster, J. J. (1992) Parametric versus non-parametric methods for estimating cure rates based on censored survival data. Statistics in Medicine, 11(7), 931–937.

Farewell, V. T. (1982). The use of mixture models for the analysis of survival data with long-term survivors. Biometrics, 38(4), 1041–1046.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., Rubin, D. B. (2013). Bayesian Data Analysis, 3o edn. Chapman and Hall/CRC.

Gieser, P. W., Chang, M. N., Rao, P. V, Shuster, J. J., Pullen, J. (2014). Modelling cure rates using the Gompertz model with covariate information. Statistics in Medicine, 17(8), 831–839.

Henningsen, A., Toomet, O. (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics, 26(3), 443–458.

Klein, J. P., Moeschberger, M. L. (2005). Survival analysis: techniques for censored and truncated data, 2o edn. Springer Science & Business Media.

Lambert, P. C., Thompson, J. R., Weston, C. L., Dickman, P. W. (2007). Estimating and modeling the cure fraction in population- based cancer survival analysis. Biostatistics, 8(3), 576–594.

Littell, R. C., Milliken, G. A., Stroup, W. W., Wolfinger, R. D., Schabenberger, O. (2006). SAS for Mixed Models, 2o edn. SAS Institute.

Lunn, D. J., Thomas, A., Best, N., Spiegelhalter, D. (2000). WinBUGS-a Bayesian modelling framework: concepts, structure, and extensibility. Statistics and Computing, 10(4), 325–337.

Maller, R. A., Zhou, X. (1996). Survival Analysis with Long-Term Survivors, Wiley.

R Development Core Team (2009). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria, URL http://www.R-project.org.

Rocha, R., Nadarajah, S., Tomazella, V., Louzada, F., Eudes, A. (2015). New defective models based on the Kumaraswamy family of distributions with application to cancer data sets. Statistical Methods in Medical Research, 1–23.

Rocha, R., Nadarajah, S., Tomazella, V., Louzada, F. (2017). A new class of defective models based on the Marshall–Olkin family of distributions for cure rate modeling. Computational Statistics & Data Analysis, 107, 48–63.

Rocha, R. F., Tomazella, V. L. D., Louzada, F. (2014). Bayesian and classic inference for the Defective Gompertz Cure Rate Model. Revista Brasileira de Biometria, 32(1), 104–114.

Spiegelhalter, D. J., Best, N. G., Carlin, B. P., Linde, A. (2014). The deviance information criterion: 12 years on. Journal of the Royal Statistical Society: Series B Statistical Methodology, 76(3), 485–493.

IR