The Log-Kumaraswamy Generalized Gamma Regression Model with Application to Chemical Dependency Data

The five parameter Kumaraswamy generalized gamma model (Pascoa et al., 2011) includes some important distributions as special cases and it is very useful for modeling lifetime data. We propose an extended version of this distribution by assuming that a shape parameter can take negative values. The new distribution can accommodate increasing, decreasing, bathtub and unimodal shaped hazard functions. A second advantage is that it also includes as special models reciprocal distributions such as the reciprocal gamma and reciprocal Weibull distributions. A third advantage is that it can represent the error distribution for the log-Kumaraswamy generalized gamma regression model. We provide a mathematical treatment of the new distribution including explicit expressions for moments, generating function, mean deviations and order statistics. We obtain the moments of the log-transformed distribution. The new regression model can be used more effectively in the analysis of survival data since it includes as submodels several widely-known regression models. The method of maximum likelihood and a Bayesian procedure are used for estimating the model parameters for censored data. Overall, the new regression model is very useful to the analysis of real data.


Introduction
Standard lifetime distributions usually present very strong restrictions to produce bathtub curves, and thus appear to be unappropriate for data with this characteristic.The gamma distribution is the most popular model for analyzing skewed data.The generalized gamma (GG) (Stacy, 1962) distribution includes as special models the exponential, Weibull, gamma and Rayleigh distributions, among others.It is suitable for modeling data with hazard rate function of different forms (increasing, decreasing, bathtub and unimodal) and then it is useful for estimating individual hazard functions and both relative hazards and relative times (Cox, 2008).The GG distribution has been used in several research areas such as engineering, hydrology and survival analysis.In fact, Ortega et al. (2003) discussed influence diagnostics in GG regression models, Nadarajah and Gupta (2007) applied this distribution to drought data and Cox et al. (2007) presented a parametric survival analysis and taxonomy of GG hazard functions.For X 1 and X 2 independent GG random variables, Ali et al. (2008) derived the exact distributions of the product X 1 X 2 and quotient X 1 /X 2 and provided applications of their results to drought data from Nebraska.Further, Gomes et al. (2008) focused on parameter estimation, Ortega et al. (2008) compared three types of deviance component residuals in GG regression models under censored observations, Cox (2008) discussed and compared the F-generalized family with the GG model, Almpanidis and Kotropoulos (2008) presented a text-independent automatic phone segmentation algorithm based on this distribution and Nadarajah (2008a) analyzed some incorrect references with respect to its use in electrical and electronic engineering.More recently, Barkauskas et al. (2009) studied the noise part of a spectrum as an autoregressive moving average (ARMA) model with the innovations having the GG distribution and Malhotra et al. (2009) provided a unified analysis for wireless system over generalized fading channels that is modeled by a two parameter GG model.Also, Ortega et al. (2009) proposed a modified GG regression model to allow the possibility that long-term survivors may be presented in the data and Cordeiro et al. (2010) defined the exponentiated generalized gamma (EGG) distribution.This distribution due to its flexibility in accommodating many forms of the risk function seems to be an important model that can be used in a variety of problems in survival analysis.
The Kumaraswamy generalized gamma (KGG) distribution (Pascoa et al., 2011) can model four types of the failure rate function (i.e.increasing, decreasing, unimodal and bathtub) depending on the values of its parameters.It is also suitable for testing goodness-of-fit of some sub-models, such as the EGG, GG, EW, Weibull and GR distributions.In this paper, we define an extended form of the KGG distribution to cope with several reciprocal type distributions and study some of its structural properties.We extend the model by Pascoa et al. (2011), since a shape parameter can take negative values.The method of maximum likelihood is used for estimating the model parameters.Unless otherwise stated, all of the results presented in the paper are new and original.It is expected that they could encourage further research on the new distribution.
Different forms of regression models have been proposed in survival analysis.Among them, the location-scale regression model (Lawless, 2003) is distinguished since it is frequently used in clinical trials.Here, we propose a location-scale regression model based on the new distribution called the log-Kumaraswamy generalized gamma (LKGG) regression model, which is a feasible alternative for modeling the four types of failure rate functions.
The rest of the paper is organized as follows.In Section 2, we define an extended version of the KGG distribution.In Section 3, we derive its ordinary moments, moment generating function (mgf), order statistics and their moments.In Section 4, we define the LKGG distribution and derive an explicit expression for its moments.In Section 5, we propose the LKGG regression model for analysis of censored data.We estimate the model parameters by maximum likelihood, derive the observed information matrix and discuss a Bayesian methodology to estimate the model parameters.In Section 6, we illustrate the potentiality of the new regression model by means of a real data set in chemical dependency.Section 7 ends with some concluding remarks.et al. (2011) defined the KGG distribution with five positive parameters α, τ , k, λ and ϕ to extend the GG and EGG distributions pioneered by Stacy (1962) and Cordeiro et al. (2010), respectively.The KGG probability density function (pdf) (for t > 0) is given by

Pascoa
where is the incomplete gamma function.The incomplete gamma function defined here for positive real k and x can be developed into the context of holomorphic functions for any complex k and x.Complex analysis shows how properties of the real incomplete gamma functions extend to their holomorphic counterparts.In the density function (1), α is a scale parameter and τ , k, λ and ϕ are shape parameters.The Weibull, GG and EGG distributions are special models of (1) when λ = k = 1, λ = ϕ = 1 and ϕ = 1, respectively.The KGG distribution approaches the log-normal distribution when λ = ϕ = 1 and k → ∞.Now, we define an extended form of the density function (1) (for t > 0) given by where τ is not zero and the other parameters are positive.The terminology Kumaraswamy generalized gamma (KGG) distribution is maintained for (2).For a random variable T having this density, we write T ∼KGG(α, τ, k, λ, ϕ).For τ > 0, (2) reduces to (1) and for τ > 0 and λ = ϕ = 1, it is identical to the GG distribution.(2011).Plots of the KGG density function for selected values of τ > 0 and τ < 0 are displayed in Figure 1.
This equation holds for any ϕ > 0. From expansion (27) (given in Appendix A), we obtain For any real λj + m > 0, we can write .
Applying the binomial expansion, (6) can be rewritten as .
From (31) (given in Appendix A) and replacing ∞ l=0 l q=0 by ∞ q=0 ∞ l=q , the density function f (t) can be expressed as a double linear combination of GG density functions where g α,τ,k(1+q)+d (t) is the GG(α, τ, k(1 + q) + d) density function given by ( 4) and the coefficients w d,q are where the quantity s m (λ − 1) is given by (28) and the quantities c q,d can be determined from the recurrence (30) (see Appendix A).Clearly, ∞ d,q=0 w d,q = 1.( 7) is the main result of this section.Some KGG mathematical properties (for example, the ordinary, inverse, factorial and incomplete moments and generating function can be obtained from (7) and those GG properties.

General Properties
Henceforth, let T and Y be random variables having the densities (2) and (4), respectively.

Moments
Here, we obtain an infinite representation for the rth ordinary moment of T , say µ r = E(T r ).The rth moment of the GG(α, τ, k) distribution is Based on (7), µ r can be expressed as where the quantity s m (λ) is given by ( 28) and the coefficients c q,d can be determined from (30).

Generating Function
Here, we provide an explicit expression for the mgf of Y , say M α,τ,k (s) = E(e sY ), for τ > 1, based on the Wright function (Wright, 1935).We are unable to obtain a closed-form expression for M α,τ,k (s) when τ < 0. We can write Setting u = t/α, we have Now, we assume τ > 0. Expanding the first exponential in Taylor series and using ∞ 0 u kτ +m−1 e −u τ du = τ −1 Γ(m/τ + k), we obtain (9) holds for any real τ > 0. However, for τ > 1, it can be simplified by considering the Wright generalized hypergeometric function (Wright, 1935) defined by x m m! .
This function exists if 1 + q j=1 B j − p j=1 A j > 0. By combining the last two equations, we have Finally, the mgf of T (for τ > 1) can be written from ( 7) and (10) as (9)-( 11) are the main result of this section.

Order Statistics
The density function f i:n (t) of the ith order statistic, say where f (t) and F (t) are the pdf and cdf of the KGG distribution, respectively and B(p, q) = Γ(p)Γ(q)/Γ(p + q) denotes the beta function.We readily obtain using the binomial expansion Now, we derive an expression for the density function of the KGG order statistics in terms of the baseline density function multiplied by a power series of G(t) = γ 1 (k, (t/α) τ ).Setting ϕ = λ = 1 in (3), it is clear that the cdf of the GG(α, τ, k) distribution is G(t) = γ 1 (k, (t/α) τ ) for τ > 0 and 1 − G(t) = 1 − γ 1 (k, (t/α) τ ) for τ < 0. This result enables us to obtain the ordinary moments of the KGG order statistics as infinite weighted sums of convenient quantities defined by These quantities δ s,r can be easily integrated numerically.They represent the probability weighted moments (PWMs) of Y when τ > 0. We shall consider two distinct cases: τ > 0 and τ < 0.
• For τ > 0 An expansion for F (t) i+j 1 −1 follows from (3) as Using the series expansion for 1 − G(t) λ l 1 ϕ , we obtain , where the coefficients s r (m 1 λ) can be obtained from (28) as where Interchanging the sums, we obtain and then The last equation can be rewritten as where Then, for τ > 0, E (T s i:n ) can be expressed as where δ s+τ (kq+d),r was defined before.
• For τ < 0 From (3), we can write The density function of the KGG order statistics becomes and then The last equation can be rewritten as where m(d, q, i, j 1 , n) is given by ( 13).So, the moments of the KGG order statistics for τ < 0 can be expressed as m(d, q, i, j 1 , n)q r (i, j 1 , λ, ϕ) δ s+τ (kq+d),r .
(12) and ( 14) are the main results of this section.

The LKGG Distribution
Henceforth, let T be a random variable following the KGG density function ( 2) where −∞ < y, µ < ∞, σ > 0, λ > 0, ϕ > 0 and q is different from zero.We consider an extended form including the case q = 0 (Lawless, 2003).Thus, the density of Y can be expressed as where Φ(•) is the standard normal cumulative distribution.We refer to (15) as the log-Kumaraswamy generalized gamma (LKGG) distribution, say Y ∼ LKGG(µ, σ, q, λ, ϕ), where µ ∈ R is the location parameter, σ > 0 is the scale parameter and q, λ and ϕ are shape parameters.For q = 0 and λϕ = 1 and q = 1, we have from ( 15) the skew normal and extreme value distributions, respectively.For ϕ = 1 and λ = ϕ = 1, we obtain the log-exponentiated generalized gamma (Ortega et al., 2012) and log-gamma generalized distributions, respectively.The case λ = ϕ = 1 and q = −1 corresponds to the log-inverse Weibull distribution.Thus, Setting µ = 0 and σ = 1, the plots of the density function ( 15) for selected values of λ and ϕ, when q < 0, q > 0 and q = 0, are displayed in Figure 2.These plots clearly indicate that the LKGG distribution could be very flexible for modeling its kurtosis.• For q > 0, the survival function of Y is given by ; • For q < 0, the survival function of Y is given by .
• For q = 0, we can write Using (5), we have The PWM (for j and p non-negative integers) of the standardized normal distribution is where φ(z) is the standard normal density.Setting y = µ + σz and using (27) in Appendix A for Φ (λ(m+1)−1) (z), we obtain The standard cumulative normal can be expressed as Using the binomial expansion and interchanging terms, we obtain From the power series for the error function erf(•) the last integral follows from ( 9)-(11) of Nadarajah (2008b).When j + p − l is even, we have where is given in terms of the Lauricella function of type A (Exton, 1978;Aarts, 2000) defined by and (a) i = a(a + 1) • • • (a + i − 1) is the ascending factorial with the convention that (a) 0 = 1.Numerical routines for the direct computation of the Lauricella function of type A are available, see Exton (1978) and Mathematica (Trott, 2006).Plots of the skewness and kurtosis for selected values of µ, σ, q, λ and ϕ are displayed in Figures 3 and 4, respectively.for some values of λ with µ = 1.5, σ = 2 and q = 0.5

The LKGG Regression Model
In many practical applications, the lifetimes are affected by explanatory variables such as the cholesterol level, blood pressure, weight and many others.Parametric models to estimate univariate survival functions for censored data regression problems are widely used.A parametric model that provides a good fit to lifetime data tends to yield more precise estimates of the quantities of interest.If Y ∼ LKGG(µ, σ, q, λ, ϕ), we define the standardized random variable Z = (Y − µ)/σ having density function given by We write Z ∼ LKGG(0, 1, q, λ, ϕ).Further, we propose a linear location-scale regression model linking the response variable y i and the explanatory variable vector where the random error z i has density function ( 19), T is a known model matrix.The LKGG model (20) opens new possibilities for fitted many different types of data.It contains as special models some well-known regression models.For ϕ = 1, we obtain the log-exponentiated generalized gamma (LEGG) regression model (Ortega et al., 2012).The case λ = 1 leads to the log-exponentiated complement generalized gamma (LECGG) regression model, whereas ϕ = λ = 1 yields the log-gamma generalized (LGG) regression model.For ϕ = λ = q = 1, we obtain the classical log-Weibull (or extreme value) regression model (see, Lawless, 2003).If σ = 1 and σ = 0.5, in addition to ϕ = λ = q = 1, the new regression model reduces to the exponential and Rayleigh regression models, respectively.The case λ = ϕ = q = 1 refers to the log-exponentiated Weibull (LEW) regression model (Mudholkar et al., 1995).See, also, Cancho et al. (1999Cancho et al. ( , 2009)), Ortega et al. (2006) and Hashimoto et al. (2010).If σ = 1, in addition to q = 1, it coincides with the log-exponentiated exponential regression model.If σ = 0.5, in addition to q = 1, it gives the log-generalized Rayleigh regression model.For λ = 1, we have the log-gamma generalized (LGG) regression model (Lawless, 2003).Recently, the LGG distribution has been used in several research areas.

Maximum Likelihood Estimation
Consider a sample (y 1 , x 1 ), • • • , (y n , x n ) of n independent observations, where each random response is defined by y i = min{log(t i ), log(c i )}.We assume noninformative censoring such that the observed lifetimes and censoring times are independent.Let F and C be the sets of individuals for which y i is the log-lifetime or log-censoring, respectively.Conventional likelihood estimation techniques can be applied here.The log-likelihood function for the vector of parameters θ = (q, λ, ϕ, σ, β T ) T from model ( 20) has the form l(θ) is the density function (15) and S(y i ) is the survival function ( 16) of Y i .The total log-likelihood function for θ can be partitioned as where v is the number of uncensored observations (failures) and z i = (y i −x T i β)/σ.The maximum likelihood estimate (MLE) θ of the vector of the model parameters can be calculated by maximizing the log-likelihood (21).Initial values for β and σ are obtained by fitting the Weibull regression model with λ = 1, ϕ = 1 and q = 1.If ẑi = (y i − x T i β)/σ, the fit of the LKGG model yields the estimated survival function for y i
Further, we can use the likelihood ratio (LR) statistic for comparing some sub-models with the LKGG model.We consider the partition θ = (θ T 1 , θ T 2 ) T , where θ 1 is a subset of parameters of interest and θ 2 is a subset of remaining parameters.The LR statistic for testing the null hypothesis H 0 : 1 is given by w = 2{ ( θ) − ( θ)}, where θ and θ are the estimates under the null and alternative hypotheses, respectively.The statistic w is asymptotically (as n → ∞) distributed as χ 2 k , where k is the dimension of the subset of parameters θ 1 of interest.

A Bayesian Analysis
As an alternative analysis, we use the Bayesian method which allows for the incorporation of previous knowledge of the parameters through informative priori density functions.When this information is not available, we can consider a noninformative prior.In the Bayesian approach, the information referring to the model parameters is obtained through a posterior marginal distribution.In this way, two difficulties usually arise.The first refers to attaining marginal posterior distribution, and the second to the calculation of the moments of interest.Both cases require numerical integration that, many times, do not present an analytical solution.Here, we use the simulation method of Markov Chain Monte Carlo (MCMC), such as the Gibbs sampler and Metropolis-Hastings algorithm.
Since we have no prior information from historical data or from previous experiment, we assign conjugate but weakly informative prior distributions to the parameters.Since we assume informative (but weakly) prior distribution, the posterior distribution is a well-defined proper distribution.Here, we consider that the elements of the parameter vector are independent and that the joint prior distribution of all unknown parameters has density function given by π(λ, ϕ, σ, q, β) ∝ π(λ) × π(ϕ) × π(σ) × π(q) × π(β). (23) , where Γ(a i , b i ) denotes a gamma distribution with mean a i /b i , variance a i /b 2 i and density function where υ > 0, a i > 0 and b i > 0, N (µ s , σ 2 s ) denotes the normal distribution with mean µ s , variance σ 2 s and density function given by where x, µ s ∈ R and σ 2 s > 0. All hyper-parameters are specified.Combining the likelihood function ( 21) and the prior distribution ( 23), the joint posterior distribution for λ, ϕ, σ, q and β T s reduces to π(λ, ϕ, σ, q, β|y) ∝ λϕ(−q)(q −2 ) q −2 σΓ(q −2 ) The joint posterior density above is analytically intractable because the integration of the joint posterior density is not easy to perform.So, the inference can be based on MCMC simulation methods such as the Gibbs sampler and Metropolis-Hastings algorithm, which can be used to draw samples, from which features of the marginal distributions of interest can be inferred.In this direction, we first obtain the full conditional distributions of the unknown parameters given by π(λ|y, ϕ, σ, q, β) Since the full conditional distributions for λ, ϕ, σ, q and β do not have closedform, we require the use of the Metropolis-Hastings algorithm.The joint posterior density for q > 0 is defined analogously, considering the likelihood function ( 21) for q > 0 and the prior distribution (23).Combining the likelihood function ( 21) and the prior distribution (23), the joint posterior distribution for λ, ϕ, σ and β T s for q = 0 reduces to The conditional distributions for these parameters taking q = 0 are given by The MCMC computations were implemented in the statistical software package R.

Application: Chemical Dependency Data
We demonstrate the potentiality of the LKGG regression model applied to chemical dependency patients as investigated by Pascoa (2008).The data were provided by the Association Mother Brave, localized in the city of Caratinga, MG, Brazil, where n = 141 residents were considered drug addicts from 2000 to 2005.The response variable was the time spent in the community until the withdraw of the treatment, considering that each resident remains in the community for a maximum period of 270 days without any contact with drugs or free drugs.The resident who achieves this goal was regarded as a censored observation.Chemical dependency is a chronic disease of the brain (Kalivas and Volkow, 2005).One of its main characteristics is the relapse phenomenon, whereby sufferers return to their consumption of the problem substance, followed by a renewed attempt to stop or reduce this consumption (Brandon et al., 2007;Koob and Le Moal, 1997).Relapse does not mean failure of treatment.Rather it is a part of the rehabilitation process (Marlatt, 2001).The patient is in relapse when he or she presents all or part of the dysfunctional behavior patterns that were present before treatment.Thus, relapse starts before actual renewed consumption.It is marked by the occurrence of facilitating stimuli with cognitive, emotional, physical and social aspects (Maisto et al., 2003).Behavioral changes and exposure to high-risk situations generally precede resumed consumption of the substance (Miller et al., 1996).Even the most motivated patients can present relapse episodes (Baker et al., 2004).Nevertheless, relapse is predictable and avoidable (Brandon et al., 2007).As the time of abstinence increases, the chance of relapse decreases, although this risk is never totally eliminated during the lifetime of a person with an addiction.Studies of the risk factors associated with relapse in chemical dependency are very important, because their findings can increase the probability of suitable early clinical interventions.Furthermore, risk studies allow identification of protective factors to reduce vulnerability and favor resistance to the temptation to relapse.The variables involved in the study are: • t i : the time spent in the community until the withdrawal of treatment (days); • δ i : the censoring indicator (0 = censoring, 1 = lifetime observed); • x i1 : the marital status (0 = single, 1 = married); • x i2 : schooling (0 = no education or incomplete bases, 1 = elementary education and higher education); • x i3 : age (0 =< 30 years, 1 =≥ 30 years).
We fit the LKGG regression model where the errors z 1 , • • • , z 141 are independent random variables having density function ( 19).Table 1 lists the MLEs of the parameters for the LKGG regression model fitted to the current data using the NLMixed procedure in SAS.This model involves an extra parameter which gives it more flexibility to fit these data.Initial values for β, q and σ are taken from the fit of the LGG regression model with λ = 1 and ϕ = 1.The explanatory variables x 1 , x 2 and x 3 are marginally significant for the LKGG model at the significance level of 5%.Cox (1972) proposed a very useful regression model for analyzing censoring failure times, where the random variable of interest represents failure time and the failures times are assumed identically distributed in some specified form.He noted that if the proportional hazards assumption holds (or, is assumed to hold) then it is possible to estimate the effect parameter(s) without any consideration of the hazard function (non-parametric approach).This approach to survival data is called proportional hazards model.The Cox model may be specialized if a reason exists to assume that the baseline hazard follows a parametric form.In this case, the baseline hazard can be replaced by a parametric density.Typically, we can then maximize the full likelihood which greatly simplifies model-fitting and provides interpretability at the cost of flexibility.
Let R(t i ) be the set of individuals at risk at time t i .Conditionally on the risk sets, the required likelihood L(β) can be expressed as where δ i is the censoring indicator.
The MLE β of β can be calculated by maximizing the likelihood function (25) using the R software.In Table 2, we list the estimates, corresponding standard errors and p-values for the fitted Cox regression model.Explanatory variables x 2 and x 3 are marginally significant at the 5% significance level.Further, Table 3 gives the Akaike Information Criterion (AIC), Consistent Akaike Information Criterion (CAIC) and Bayesian Information Criterion (BIC) to compare the LKGG model and Cox proportional hazard regression models.The LKGG regression model outperforms the other models irrespective of the criteria and it can be used effectively in the analysis of these data.So, the proposed model is a great alternative to model survival data.In order to assess if the model is appropriate, we fit the LKGG regression model for each explanatory variable.In Figures 5a, b, c, d, we plot the empirical survival function and the estimated survival function ( 22) for each explanatory variable.We can conclude that the LKGG regression model provides a good fit to these data.

Bayesian analysis
The following independent priors were considered to perform the Metropolis-Hastings algorithm: λ ∼ Γ(0.01, 0.01), ϕ ∼ Γ(0.01, 0.01), σ ∼ Γ(0.01, 0.01), q ∼ N (0, 10) and β s ∼ N (0, 10), so that we have a vague prior distribution.Considering these prior density functions, we generate two parallel independent runs of the Metropolis-Hastings with size 100, 000 for each parameter, disregarding the first 10, 000 iterations to eliminate the effect of the initial values and, to avoid correlation problems, we consider a spacing of size 10, obtaining a sample of size 9, 000 from each chain.To monitor the convergence of the Metropolis-Hastings, we perform the methods suggested by Cowles and Carlin (1996).We use the between and within sequence information, following the approach developed in Gelman and Rubin (1992), to obtain the potential scale reduction, R. In all cases, these values were close to one, indicating the convergence of the chain.The approximate posterior marginal density functions for the parameters are presented in Figure 6.In Table 4, we report posterior summaries for the parameters of the LKGG model.We note that the values for the means a posteriori (Table 4) are quite close (as expected) to the MLEs given in Table 1.SD denotes the standard deviation from the posterior distributions of the parameters and HPD represents the 95% highest posterior density (HPD) intervals.

Concluding Remarks
We introduce an extended form of the Kumaraswamy generalized gamma (KGG) distribution (Pascoa et al., 2011) for which the hazard rate function accommodates the four types of shape forms, i.e. increasing, decreasing, bathtub and unimodal.The KGG model includes as special cases several useful lifetime models.We derive explicit expressions for its moments, generating function, order statistics and their moments.Further, we also introduce the called log-Kumaraswamy generalized gamma (LKGG) distribution and obtain explicit expressions for its moments.Based on this new distribution, we define the LKGG regression model which is very suitable for modeling censored and uncensored lifetime data.The new regression model allow us to test as special models the goodness of fit of some widely known regression models.Hence, it represents a good alternative for lifetime data analysis.We estimate the model parameters using maximum likelihood and a Bayesian approach.We demonstrate by means of an application to real data that the LKGG model can produce better fits than some well-known models.
which always converges since 0 < γ 1 (k, x) < 1.Hence, We can substitute where We use the power series for the incomplete gamma ratio function given by By application a result by Gradshteyn and Ryzhik (2000, Section 0.314) for a power series raised to a positive integer q, we obtain Here, the coefficients c q,d (for d = 1, 2, • • • ) are determined from the recurrence equation where c q,0 = a q 0 and a p = (−1) p /[(k + p)p!].So, the coefficient c q,d can be calculated from c q,0 , • • • , c q,d−1 , and then it can be expressed as a function of the quantities a 0 , • • • , a d , although it is not necessary for programming numerically the expansions.Further, using (29), we obtain Appendix B: Matrix of second derivatives − L(θ) Here, we provide formulas for the second-order partial derivatives of the loglikelihood function.After some algebraic manipulations, we obtain × exp(−w i ) + q −2 (q −2 − 1)w q −2 −1 i exp(qz i ) + q −1 Γ(q −2 ) w q −2 −1 i z i × exp(qz i − 2w i )[γ 1 (q −2 , w i )] −1 − q −1 Γ(q −2 ) w q −2 −1 i z i exp(qz i − w i ) where z i = (y i − x T i β)/σ and w i = q −2 exp(qz i ).

Figure 2 :
Figure 2: The LKGG density curves: (a) For some values of q > 0. (b) For some values of q < 0. (c) For some values of q = 0

Figure 3 :Figure 4 :
Figure 3: Skewness and kurtosis of the LKGG distribution as a function of λ for some values of ϕ with µ = 1.5, σ = 2 and q = 0.5

Figure 5 :
Figure 5: Kaplan-Meier curves stratified by explanatory variable and estimated survival functions to the chemical dependency data: (a) The marital status explanatory variable.(b) Schooling explanatory variable.(c) Age explanatory variable

Figure 6 :
Figure 6: Approximate posterior marginal densities for the parameters from the LKGG model for the chemical dependency data

Table 1 :
MLEs of the parameters, standard errors, p-values and 95% confidence intervals for the LKGG model fitted to the chemical dependency data

Table 2 :
Estimates of the Cox model fitted to the data of chemical dependency and corresponding hazard ratios (HR j )

Table 3 :
The AIC, CAIC and BIC statistics for comparing the LKGG model and Cox proportional hazard regression models

Table 4 :
Posterior summaries for the parameters from the LKGG model for the chemical dependency data