The Poisson-exponential regression model under different latent activation schemes

The Poisson-exponential regression model under different latent activation schemes

Francisco Louzada^I,*; Vicente G. Cancho^I; Gladys D.C. Barriga^II

^ISME, ICMC - Universidade de São Paulo, São Carlos, São Paulo, Brasil. E-mails: [e-mail] / [e-mail]

^IIFEB - Universidade Estadual Paulista. Email: [e-mail]

ABSTRACT

In this paper, a new family of survival distributions is presented. It is derived by considering that the latent number of failure causes follows a Poisson distribution and the time for these causes to be activated follows an exponential distribution. Three different activationschemes are also considered. Moreover, we propose the inclusion of covariates in the model formulation in order to study their effect on the expected value of the number of causes and on the failure rate function. Inferential procedure based on the maximum likelihood method is discussed and evaluated via simulation. The developed methodology is illustrated on a real data set on ovarian cancer.

Mathematical subject classification: 62N01, 62N99.

Key words: activation schemes, exponential distribution, poisson distribution, survival analysis.

1 Introduction

As well known, although the exponential distribution provides a simple, elegant and close form solution to many problems in lifetime testing and reliability studies, it does not provide a reasonable parametric fit for some practical applications where the underlying failure rates are nonconstant, presenting monotone shapes. In recent years, in order to overcame such problem, new classes of models were introduced based on modifications of the exponential distribution. Gupta and Kundu (1999) proposed a generalized exponential distribution, this family can accommodate data with increasing and decreasing failure rate function. Adamidis and Loukas (1998) introduced the exponential-geometric distribution with decreasing failure rate. Kus (2007) proposed a two-parameter distribution known as exponential-Poisson distribution, which has decreasing failure rate. Tahmasbi and Rezaei (2008) proposed another modification of the exponential distribution with decreasing failure rate function. This distribution is known as logarithmic-exponential distribution. Cancho et al. (2011) introduced the Poisson-exponential distribution with increasing failure rate. This model is derived in a complementary risks scenario (Louzada-Neto, 1999), where the lifetime associated with a particular risk is not observable, rather we observe only the maximum lifetime value among all risks.

In this paper we proposed a class of hierarchical models with latent competing risks and different activation schemes, where the lifetime associated with a particular risk is not observable, rather we observe only the minimum, maximum or a randomly lifetime value among all risks, with the lifetimes following an exponential distribution. We consider, that the number of latent competing risks (or causes) is modeled by a zero truncated Poisson distribution, with our formulation has several particular cases, including the models proposed by Kus (2007) and Cancho et al. (2011).

The paper is organized as follows. In Section 2, we introduce the new distribution and present its properties. Section 3 we carry out inference for the Poisson-exponential regression model, discuss some measures of model selection. Section 4 we present the results of a simulation study and in Section 5 the methodology is illustrated on a real ovarian cancer data set. Some final comments are presented in Section 6.

2 The model

Let M denote the unobservable number of causes related to the occurrence of an event of interest for an individual in the population. Assume that M has a zero truncated Poisson distribution with probability mass function given by

The time for the j^th cause to produce the event of interest is denoted by T_j, j = 1,..., M. We assume that, conditional on M, the T_j are independent and identically distributed with cumulative distribution function G(t). Also, we assume that T₁, T₂, ... are independent of M. Further, consider the order statistics of the T_j's, T₍₁₎< T₍₂₎,..., T_(R)< ... < T_(M). The observed lifetime can be defined by the random variable Y = T_(R), where R = 1, 2,..., M is dependent of M. In many biological process, R indicates the resistance factorof the immune system of the individual. In other words, as in Cooner et al. (2007), R out of N causes are required to produce the event of interest. The resistance factor can be a fixed constant, a function of M or a random variable specified through a conditional distribution on M.

In this paper we deal with three specifications for R . Firstly, we assume that given M, the conditional distribution of R is uniform on {1,2,..., M} (random activation scheme). Under this setup, the distribution function for Y given M = m and R = r, is given by

Also, we can demonstrate that the marginal distribution function of Y is given by

where B(x, m, G(y)) is probability mass function (pmg) of binomial distribution with parameters, m and G(y) and p_m is given in (1).

Note that, in (3), the distribution function of Y is the same as the distribution of the latent random variables T_j's. If we consider that R = r (fixed), then the marginal distribution of Y is given

where I B(x; a, b) is incomplete beta function and p_m is given in (1).

As a second setup, the so-called first activation scheme (Cooner et al., 2007), we suppose that the event of interest happens due to any one of the possible causes, but for R = 1. Then, the observed failure is Y = Z₍₁₎ = min{T₁,...,T_M}. In this case the marginal distribution of Y in (4) is given by

In a third scenario, also known as the last activation scheme (Cooner et al., 2007), the event of interest only takes place after all the M causes have been occurred, so that R = M and the observed failure time is Y = Z_(M) = max{T₁,...,T_M}. The marginal distribution of Y in (4)is given by

The relationship between the distribution functions in (3), (5) and (6) isdescribed in next proposition.

Proposition 1 Under some conditions on the models in (3), (5) and (6), for any distribution function G(y), F_last(y) < F_random(y) < F_first(y) for y > 0.

Proof. We know that A(y) = is a decreasing function in y, indeed, let

and lim_y→0B(y) = 0 thus B(y) < 0, ∀y. So

and as

then,

Similarly, we can prove the F_last(y) < F_random(y) ∀y.

Considering different choices for the distribution of latent random variables T_j's, new families of distribution can be obtained. For instance, the exponential-Poisson distribution proposed by Kus (2007) is obtained by considering that the variable T_j's follows an exponential distribution with failure rate function, λ > 0, under the first-activation scheme, members of this family of distribution have decreasing failure rate functions. For properties of this distribution interested readers can refer to Kus (2007). The counterpart distribution is proposed by Cancho et al. (2011), with observed time failure given by Y = max(Z₁,..., Z_M) (last-activation), which has increasing failure ratefunction.

We may consider an extension of the exponential-Poisson (EP) models under latent competitive risks and different activation schemes, by including covariates in the modeling. Particularly, in order to study the effect of covariates on the value expected of number of causes (E(M) = η), we change the parametrization of the models in (5) and (6) in order to include E(M) = η in the expressions. The solution in the parameter θ of the equation, θ/(1 – e^-θ) = η is given by W₀ = η+ W(– ηe^-η ), where W(·) stands for theLambert W function (Corless et al., 1996). The Lambert W(·) function has been used in the solution of many problems in mathematics and engineering, and then the mathematical softwares Mathematica, Maple and R have introduced it in their packages.

In this context, we obtain the survival, failure rate and density functionspresented in Table 1 for the unobserved failure time T following exponentialdistribution with failure rate λ > 0. Figure 1 portrays distinct behaviors of the survival functions (top panel) and the failure rate functions (bottom panel) in Table 1. These plots illustrate the flexibility afforded by our proposal.

3 Inference

Let us consider the situation when the failure time Y in Section 2 is not completely observed and it is subjected to right no-informative censoring. Let C_idenote the censoring time. In a sample of size n, we then observe

t_i = min{Y_i, C_i} and δ_i = I(Y_i< C_i),

where δ_i = 1 if T_i is a failure time and δ_i = 0 if it is right censored, for i = 1,..., n.

Let x_i and z_i denote the covariate vector. We relate the model parameters η and λ to covariates x_i and z_i by adopting the following link functions,

i = 1,..., n, where β and γ denote vectors of coefficients, and x_iand z_i can be equals. This model will be referred to as the Poisson-exponential (PE) regression model.

With the expression (7) we can write the likelihood of = (β^T,γ^T)^T under right non-informative censoring as

where D = (t,δ, x, z), t = (t₁,..., t_n)^T and δ = (δ₁,..., δ_n)^T, whereas f(·;) and S(·;) are the density and surviving functions in Table 1.

From the likelihood function in (10), the maximum likelihood estimation of the parameter is carried out. Numerical maximization of the log-likelihood function () = (; D) is accomplished by using existing software. In this paper, the software R (see, R Development Core Team, 2009) was used to compute maximum likelihood estimates (MLE). The Lambert W function in Table 1 can be found in the R package emdbook. Covariance estimates for the maximum likelihood estimators may also be obtained using the Hessian matrix. For large samples, confidence intervals may be conducted by using the large sample distribution of the MLE which is a normal distribution with the covariance matrix as the inverse of the Fisher information. More specifically, the asymptotic covariance matrix is given by I^-1() with I() = –E[J()] such that J() = ∂²()/∂∂^T. Since it is not possible to compute the Fisher information matrix I() due to the censored observations (censoring is random and non-informative), it is possible to use the matrix of second derivatives of the log likelihood, -J(), evaluated at the MLE = , which is a consistent estimator. The required second derivatives are computed numerically.

Besides estimation, hypothesis testing is another key issue. Let

₁ and ₂ be two proper disjoint subsets of . We aim to test H₀:₁ = ₀₁ against H₁: ₁ ≠ ₀₁, ₂ unspecified. Let ₀ maximize L(; D) constrained to H₀ and define the log-likelihood ratio statistic as w_n = 2(()-(₀)), where (·) is the log-likelihood. Under H₀ and some regularity conditions, w_n converges in distribution to a chi-square distribution with dim(₁) degrees of freedom.

Different models can be compared penalizing over-fitting by considering the Akaike information criterion given by AIC = -2() + 2 #() and the Schwartz information criterion defined by SIC = -2() + #() log(n), where #() is the number of model parameters. The model with the smallest value of any of these criteria (among all models considered) is commonly taken as the preferred model for describing the dataset.

4 Simulation

In this section we present the results of a simulation study performed in order to evaluate the performance of the EMVs of parameters of the Poisson-exponential regression model in terms of their biases, standard errors and root mean squared errors. Without loss of generality, in this study we considered the Poisson-exponential regression model with increasing failure rate function (first-activation) with two covariates, x_i and z_i, being that, x_i generated from a Bernoulli distribution with parameter 0.5 and z_i generated from a uniform distribution on the range [0,1]. The failure times were simulated from the quantile function of model (via the inversion method) with parameter η_i = exp(β₁₀+β₁₁x_i) and λ_i = exp(β₂₀+β₂₁z_i), where β = (1, 1) and γ = (–1, –1). Censoring time were generated from an uniform distribution [0, τ], where τ controlled the proportion of censoring, which is considered here equals to 0% and 10%. We considered sample sizes n equal to 30, 50, 100 and 200. For each configuration, we conducted 1,000 replicates and then we averaged the estimates of the parameters, the standard errors and the square root of themean square errors.

The Table 2 show that the bias of the MLEs β and θ, and their standard errors and root mean squared errors become smaller when the sample size increases and the percentage of censored observations is smaller.

5 Application

In this section we reanalyze the data extracted from Collett (2003). The data describe a study of 26 women diagnosed with ovarian cancer submitted to chemotherapy. The data consist of the survival times y_i (in months) of patients and several prognostic covariates such as: x_1i: treatment (0 = single, 1 = combined), x_2i: age of patient in years, x_3i: extent of residual disease (0 = incomplete, 1 = complete) and x_4i: performance status (0 = good, 1 = poor).

Firstly, in order to identify the shape of a lifetime data failure rate functionwe shall consider, as a crude indicative, a graphical method based on the TTT plot (Aarset, 1985). In its empirical version the TTT plot is given by

where r = 1,..., n and Y _i:n represent the order statistics of the sample. It has been shown that the failure rate function is increasing (decreasing) if the TTT plot is concave (convex). Figure 2 (top panel) shows the TTT plot for the considered data, which is concave indicating an increasing failure rate function, which can be properly accommodated by a Poisson-exponential regression model with increasing failure rate. This proposal indicates that all causes were responsible for activating the event of interest (last-activation scheme).

Firstly, we consider the following EP regression model with all covariates, i.e,

and

Table 3 presents the MLEs of the coefficients. The QQ plot of the normalized randomized quantile residuals (Dunn and Smyth, 1996; Rigby and Stasinopoulos, 2005) in Figure 2 (bottom panel) suggests that the Poisson-exponential regression model yields an adequate fit. Each point in Figure 2 (bottom panel) corresponds to the median of ten sets of ordered residuals.

Considering the LR statistics in Section 3 we tested the following hypothesis, H₀:β₂ = β₃ = β₄ = γ₁ = γ₃ = γ₄ = 0, resulting, a LR statistics w_n equals to 2.18, with g.l = 6 and p-value = 0.902, indicating that the effect coefficients are no significant in the model. Thus, the final Poisson-exponential regression model becomes the one given by

Table 4 shows the MLEs and their standard errors of the parameters for final Poisson-exponential regression model. We note the covariate treatment is significant (at 5%) for the expected value of number of causes, while that the covariate age is significant (at 1%) in the failure rate.

From the standpoint of practical researchers there is interest in inferencesabout the expected number of competing causes that operate in the occurrence of the event of interest. The answer for this question is directly obtained by considering the MLEs of the final Poisson-exponential regression model parameters in Table 4, after an application of the delta method. We then obtain the estimates of the expected number of competing causes related to the occurrence of the event of interest, stratified by treatment, which are equal to 3.605 for single and 9.325 for combining.

The benefit of the combining treatment can be observed from Figure 3, which displays the surviving (top panel) and failure rate (bottom panel) functions for 56 years old (age mean) patients stratified by treatment, implying in an increasing in survival and a decreasing in failure rate when the combining treatment is used instead of considering the control treatment.

Although the Poisson-exponential regression model under the last activation scheme was firstly suggested by the TTT plot graphical procedure presented above, for sake of comparison we fitted the Poisson-exponential regression model under the other two activation schemes, namely, first and random activation schemes. We fitted these models with the same covariates as to our final model (see Table 4). The estimated AIC and SIC for the models are reported in the Table 5. Notice that the Poisson-exponential regression model under the last activation scheme outperforms the other two schemes irrespective of the criterion considered. Finally, we fitted the usual Weibull regression model, with AIC and SIC values equals to 103.90 and 108.93, respectively, which also was outperformed by our Poisson-exponential regression model under the last activation scheme.

6 Concluding remarks

In this paper we proposed the Poisson-exponential regression model, as a possible extension of the well known exponential regression model. The model was built based on a latent competing risk structure with three different activation schemes, minimum, maximum or random. We discussed parameter maximum likelihood estimation and a straightforwardly modeling comparison procedure. The flexibility of our modeling was illustrated in on a real data set on ovarian cancer.

We considered the number of latent competing risks modeled by a zero truncated Poisson distribution, directly generalizing the models proposed by Kus (2007) (minimum activation scheme) and Cancho et al. (2011) (maximum activation scheme). However, at least in principle, other distributions may be considered to represent the random behavior of the number of latent competing risks, such as the geometric, negative binomial and logarithmic distributions, or even a more general one such as the power series distribution. Chahkandi and Ganjali (2009) has unified the models proposed by Adamidis and Loukas (1998), Kus (2007) and Tahmasbi and Rezaei (2008) by compounding an exponential distribution and power series distribution, but they only considered a minimum activation scheme. A possible study involving the power series distribution in the context of our modeling with three different activation schemes shall be considered elsewhere.

Models for survival data with a surviving fraction take a prominent position in survival studies, covering situations where there are sampling units insusceptible to the occurrence of the event of interest (Rodrigues et al., 2009; Louzada et al., 2011; Perdona and Louzada, 2011). The proposed Poisson-exponential regression model for survival data in presence of a survival fraction should be considered in a future study, as well as modeling considering the shape parameter depending on covariates (Louzada-Neto, 1997). Diagnostic methods have been an important tool in survival regression analysis. Influential diagnostics should be investigated further in the context of the proposed Poisson-exponential regression model (Fachini et al., 2008).

Acknowledgments. The authors are grafeful the Editorial Boarding as well as to the referees for their comments and suggestions. The research of Francisco Louzada and Vicente G. Cancho are funded by the Brazilian organization CNPq.

Received: 22/IX/12.

Accepted: 28/IX/12.

#CAM-704/12.

References

• [1] M.V. Aarset, The null distribution for a test of constant versus "bathtub" failure rate. Scandinavian Journal of Statistics, 12(1) (1985), 55-68.
• [2] K. Adamidis and S. Loukas, A lifetime distribution with decreasing failure rate. Statistics & Probability Letters, 39(1) (1998), 35-42.
• [3] V. Cancho, F. Louzada-Neto and G. Barriga, The poisson-exponential lifetime distribution. Computational Statistics & Data Analysis, 55 (2011), 677-686.
• [4] M. Chahkandi and M. Ganjali, On some lifetime distributions with decreasing failure rate. Comput. Stat. Data Anal., 53(12) (2009), 4433-4440.
• [5] D. Collett, Modelling survival data in medical research. CRC press. ISBN 1584883251, (2003).
• [6] F. Cooner, S. Banerjee, B.P. Carlin and D. Sinha, Flexible cure rate modeling under latent activation schemes. Journal of the American Statistical Association, 102 (2007), 560-572.
• [7] R.M. Corless, G.H. Gonnet, D.E.G. Hare, D.J. Jeffrey and D.E. Knuth, On the Lambert W function. Advances in Computational Mathematics, 5 (1996), 329-359.
• [8] P.K. Dunn and G.K. Smyth, Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5(3) (1996), 236-244.
• [9] J.B. Fachini, E.M.M. Ortega and F. Louzada-Neto, Influential diagnostics forpolyhazard models in the presence of covariates. Statistical Methods and Applications, 17 (2008), 413-433.
• [10] R. Gupta and D. Kundu, Generalized exponential distributions. Australian and New Zealand Journal of Statistics, 41(2) (1999), 173-188.
• [11] D. Kus, A new lifetime distribution distributions. Computational Statistics & Data analysis, 11 (2007), 4497-4509.
• [12] F. Louzada, M. Roman and V. Cancho, The complementary exponential geometric distribution: Model, properties, and a comparison with its counterpart. Manuscript submitted to Computational Statistics & Data Analysis (2011).
• [13] F. Louzada-Neto, Extended hazard regression model for reliability and survival analysis. Lifetime Data Analysis, 3 (1997), 367-381.
• [14] F. Louzada-Neto, Poly-hazard regression models for lifetime data. Biometrics, 55 (1999), 1121-1125.
• [15] G.S.C. Perdona and F. Louzada, A general hazard model for lifetime data in the presence of cure rate. Journal of Applied Statistics, 38 (2011), 1395-1405.
• [16] R Development Core Team, R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (2009).
• [17] R.A. Rigby and D.M. Stasinopoulos, Generalized additive models for location, scale and shape (with discussion). Applied Statistics, 54(3) (2005), 507-554.
• [18] J. Rodrigues, V.G. Cancho, M. de Castro and F. Louzada-Neto, On the unification of the long-term survival models. Statistics & Probability Letters, 79 (2009),753-759.
• [19] R. Tahmasbi and S. Rezaei, A two-parameter lifetime distribution with decreasing failure rate. Computational Statistics & Data Analysis, 52(8) (2008), 3889-3901.

*

Corresponding author

Publication Dates

• Publication in this collection
  28 Nov 2012
• Date of issue
  2012

History

• Received
  22 Sept 2012
• Accepted
  28 Sept 2012