Promote sign consistency in cure rate model with Weibull lifetime

1.   Introduction

Survival analysis is a method to analyze time to event of interest assuming that all subjects will eventually experience the event ^[1]. It has been widely applied in many fields such as medicine, engineering, credit scoring, etc. ^[2,3,4]. However, in practice, a proportion of subjects may not experience the event of interest. They are considered 'cured', or in other words, the long-term survivors. For instance, in medical practice, there may be a fraction of patients cured of their disease. Therefore, the cure rate model, an extension of survival analysis, is introduced to modeling data with long-term survivors ^[5].

Let ${\bf{x}}_i^{}$ be the covariate vector of subject $i$, and $S_i^{}(t|{{\bf{x}}_i})$ be the survival probability in time $t$. The cure rate model can be written as

where ${U_i}$ is the random binary variable which indicates the cure status of subject $i$. The cure status ${U_i} = 0$ denotes subject $i$ is cured and will not experience the event, and ${U_i} = 1$ otherwise. Here ${\pi _i}\left({{{\bf{x}}_i}} \right)$ is the probability function of being cured with the vector of regression coefficients ${\boldsymbol{\alpha}}$, and $S_i^{}(t|U_i^{} = 1, {{\bf{x}}_i})$ is the survival function with the vector of regression coefficients ${\boldsymbol{\beta}}$.

As shown in (1.1), the cure rate model is composed of two parts. The first part is the incident part, which predicts the probability of being cured. The second part is the latency part, which describes the conditional survival of the uncured subjects. Since the cure rate model can predict whether subjects are cured and the time to event of the uncured subjects, it is commonly adopted when a proportion of subjects have long-term survivors ^[6].

There are numerous studies in the literature regarding many extensions of the cure rate model. Cooner et al. ^[7] proposed a flexible hierarchical cure rate model to distinguish among underlying mechanisms that lead to cure. Rodrigues et al. ^[8] assumed the number of risk factors to follow the Conway-Maxwell Poisson distribution and proposed a Conway-Maxwell Poisson cure rate model to unify some cure rate models. Li et al. ^[9] considered a mixture of linear dependent tail-free processes as the prior for the distribution of the cure rate parameter to develop a latent promotion time cure rate model. Results showed that the cure rate model incorporated penalized spline has better performance in prediction ^[10]. Georgiana proposed a Bayesian spatial cure rate model with Weibull lifetime to model spatial variability in the censoring mechanism ^[11]. Pal et al. ^[12] proposed a projected non-linear conjugate gradient algorithm for the cure rate model under a competing risks scenario. Besides, many works have developed semi-parametric and nonparametric methods to investigate the effects of covariates on the outcome. For instance, Li et al. ^[13] proposed a semi-parametric additive predictor consisting of a sum of linear and nonparametric terms in the incident part. Chen et al. ^[14] modeled the covariate effects by a nonparametric form.

However, many existing works tend to pay much less attention to the relations between the vector of regression coefficients ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ in the two model parts. Generally, they assume that there are no direct constraints on the relation of the vector of regression coefficients ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ corresponding to the same covariates ^[15,16]. In other words, these works assume that the probability of being cured and the time to event are independent, and there are no direct constraints between the coefficients ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$.

Notably, in practical applications, the two model parts are quite related. The incident part describes the probability of cure (infinity survival), and the latency part describes the conditional survival of the uncured subjects (finite survival). It is desirable that there may be some connections between the vector of regression coefficients ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ corresponding to the same covariates. Theoretical derivations and case studies also suggest that relaxing the assumption of no direct constraints on regression coefficients can improve the performance of the model ^[17,18]. Liu et al. ^[19] relax the assumption of no direct constraints by establishing a joint distribution of the covariates and the logarithm of the hazard rate. Fan et al. ^[20] incorporate structural effects of ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ in the cure rate model to improve the estimation accuracy and interpretability. A joint distribution or structural effects of ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ are crude yet effective way to impose the relations between the coefficients. However, the two parts of the model still describe two different aspects. The assumption of a joint distribution or structural effects are too restrictive constraints and might not work well. So we consider a more flexible model that allows the coefficients ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ can be in different distributions and magnitudes. Sign consistency penalty is proposed to promote the similarity in sign to get more interpretable results by Zhang ^[21]. In many practical cases, it is hard to interpret the results when coefficients ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ corresponding to the same covariates have conflicting signs. In this paper, we consider a sign consistency cure rate model with a sign-based penalty. In addition, the models may suffer from bad performance due to the high dimension of the data, and grouping structure arises naturally in many practical cases. A group lasso penalty is also imposed for group variable selection.

In this paper, we propose a cure rate model with group selection and sign consistency (CRGS), which can select important groups of covariates and promote similarity in the sign of coefficients. The proposed method can promote the similarity in the signs of coefficients ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ to improve interpretability. Compared to the individual variable selection approaches such as the sign consistency method in ^[22], the proposed method with the group selection approach takes into consideration the grouping structure and can lead to better prediction. Compared with previously employed approaches such as joint distribution or structural effects of coefficients, the CRGS method can avoid too strict constraints between the coefficients and lead to more consistency and hence more interpretable results.

The paper is organized as follows. In section 2, the sign consistency cure rate model with Weibull lifetime and the algorithm is introduced. Simulation is presented in section 3. Section 4 displays a real data application. Finally, section 5 discusses the conclusions.

2.   Methods

2.1.   Sign consistency cure rate model with Weibull lifetime

Consider data with $n$ subjects and $p$ covariates. Denote ${Y_i}$ as the time to event of subject $i$, that is, the time until the event of interest is observed to occur. Let ${C_i}$ be the time of right censoring, and ${\delta _i}{\text{ = }}{\rm I}\left({{Y_i} < {C_i}} \right)$ be the censoring indicator of subject $i$, where ${\delta _i}{\text{ = }}0$ for censored and ${\delta _i}{\text{ = 1}}$ for uncensored. Denote ${y_i} = \min \left({{T_i}, {C_i}} \right)$. Note that the censored subjects include the cured subjects and the uncured subjects for whom the event has not occurred at censoring time. So the cure status ${U_i}$ is unobservable. The observable data is $\left\{ {\left({{y_i}, {\delta _i}, {{{\bf{x}}_i}}} \right), i = 1, \cdots, n} \right\}$.

Denote ${\bf{x}}$ as the covariate vector with grouping structure. Let ${\bf{x}} = {({\bf{x}}_1^{\rm T}, \cdots, {\bf{x}}_J^{\rm T})^{\rm T}}$ be the covariate vector with $J$ nonoverlapping subgroups. ${{\bf{x}}_j} = {({x_{j1}}, \cdots, {x_{j{p_j}}})^{\rm T}}$ is the $j$ -th subgroup of covariate vector, with $\sum\nolimits_{j = 1}^J {{p_j} = p}$. Grouping structure arises naturally in many practical cases. Examples include the expression of a multi-level factor by a group of dummy covariates and the expression of an addictive model by several basis functions ^[23]. In addition, grouping structure can also be introduced into a model by taking advantage of prior knowledge ^[24]. For example, genes belong to the same biological pathway ^[25].

In the incident part of the cure rate model, we adopted logistic regression, a generalized linear model, to describe the probability of cure ${\pi _i}\left({{{\bf{x}}_i}} \right) = {1 \mathord{\left/ {\vphantom {1 {\left({1 + \exp \left({{\alpha _0} + {{\bf{x}}_i}^{\rm T}\alpha } \right)} \right)}}} \right. } {\left({1 + \exp \left({{\alpha _0} + {{\bf{x}}_i}^{\rm T}{\boldsymbol{\alpha}} } \right)} \right)}}$. Here ${\boldsymbol{\alpha}}= {({\boldsymbol{\alpha}} _1^{\rm T}, \cdots, {\boldsymbol{\alpha}}_J^{\rm T})^{\rm T}}$ is the vector of regression coefficients, ${{\boldsymbol{\alpha}} _j} = {({\alpha _{j1}}, \cdots, {\alpha _{j{p_j}}})^{\rm T}}$ is the $j$ -th subgroup of the coefficient vector, and ${\alpha _0}$ is the intercept.

In the latency part, let $\lambda _i^{}$ be the link function. The survival function is the chance an individual survives to time $t$ given the individual will eventually experience the event of interest, while $\lambda _i^{}$ is the probability of an individual experiencing the event in the next instant of time $t$.

We assume the time to event $t$ follow the Weibull distribution. Weibull distribution is a considerable flexibility distribution to model lifetime data ^[26]. It had been justified to be a valid lifetime distribution within the broad family of generalized gamma models ^[27,28]. Referring to ^[1,29], the survival function for the uncured subjects following the Weibull distribution can be written as

with the probability density function $f\left({t|U_i^{} = 1} \right) = \left({{r \mathord{\left/ {\vphantom {r t}} \right. } t}} \right){\left({\lambda _i^{}t} \right)^r}\exp \left({ - {{\left({\lambda _i^{}t} \right)}^r}} \right).$ Here $r > 0$ is the shape parameter (more discussions below). The link function can be written as

where ${\beta _0}$ is the intercept, ${\boldsymbol{\beta}}= {({\boldsymbol{\beta}} _1^{\rm T}, \cdots, {\boldsymbol{\beta}} _J^{\rm T})^{\rm T}}$ is the vector of regression coefficients, and the $j$ -th subgroup of the coefficient vector is ${{\boldsymbol{\beta}} _j} = {({{\boldsymbol{\beta}} _{j1}}, \cdots, {{\boldsymbol{\beta}}_{j{p_j}}})^{\rm T}}$.

For observable data $\left\{ {\left({{y_i}, {\delta _i}, {{{\bf{x}}_i}}} \right), i = 1, \cdots, n} \right\}$, the log-likelihood function is

For promoting sign consistency and variable selection, we propose the following objective function

Here ${{\rm{P}}_1}{\boldsymbol{(\alpha}}, {\boldsymbol{\beta)}}$ is the group variable selection penalty function, and ${{\rm{P}}_2}{\boldsymbol{(\alpha}}, {\boldsymbol{\beta)}}$ is the sign consistency penalty function as follows.

where $\left\| \bullet \right\|$ is the ${l_2}$ norm, ${\text{sign}}(\bullet)$ is the sign function, ${\mu _1} > 0$ and ${\mu _2} > 0$ are tuning parameters. ${\mu _1}$ and ${\mu _2}$ control the degree of penalty.

The first penalty ${{\rm{P}}_1}{\boldsymbol{(\alpha}}, {\boldsymbol{\beta)}}$ is a group lasso penalty, which can conduct group variable selection and generate more accurate estimation. Group lasso has the advantage in group variable selection performance and is commonly used in literature. It has been justified that the group lasso is more robust to noise compared with lasso when the underlying signal is strongly group-sparse ^[30]. In addition, different from methods in some existing works ^[20,22], group lasso can consider grouping structures and it has better performance compared with the group LARS and the group non-negative garrotte ^[31]. The second penalty ${{\rm{P}}_2}{\boldsymbol{(\alpha}}, {\boldsymbol{\beta)}}$ promotes the similarity in signs of coefficients ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$, which can lead to more interpretable results. The proposed method is more flexible than the methods with structural effect in ^[20] or joint distribution in ^[19] of coefficients.

2.2.   EGCD algorithm

In this section, we propose the Expectation Group Coordinate Descent (EGCD) algorithm to optimize the objective function. In the E-step, we introduce a latent unobserved ${U_i}$ to obtain a complete log-likelihood function. In the GCD-step, group coordinate descent is adopted to iteratively update a subgroup of parameters with the remaining parameters fixed at their most recent values. Since the sign function ${\text{sign}}(\bullet)$ is not differentiable and continuous, and it is hard to optimize. So we introduce the following approximation for computational feasibility ^[21,22].

with $\tau$ is a small positive constant.

The EGCD algorithm iteratively updates ${\alpha _0}$, ${\boldsymbol{\alpha}}$, ${\beta _0}$, and ${\boldsymbol{\beta}}$ in the $m$ -th iteration.

2.2.1.   E-step

In the E-step of the $m$ -th iteration, let the observation of the latent ${U_i}$ be $U_i^{[m]}$. Consider the complete data $\left\{ {({y_i}, {\delta _i}, U_i^{[m]}, {{\bf{x}}_i}), i = 1, \cdots, n} \right\}$ with the latent $U_i^{[m]}$. The complete log-likelihood is

where

Regarding the expectation of ${U_i}$, there are three possible situations of subjects. (a) ${\delta _i}{\text{ = }}0$ and ${U_i} = 0$ : censored and cured, indicates long-term survivors; (b) ${\delta _i}{\text{ = }}0$ and ${U_i} = 1$ : censored and uncured, indicates subjects who will eventually experience the event, and the event have not occurred in censoring time ${C_i}$; (c) ${\delta _i}{\text{ = }}1$ and ${U_i} = 1$ : uncensored and uncured, indicates subjects who have experienced the event. Therefore, the expectation of ${U_i}$ is 1 when uncensored (${\delta _i} = 1$). When censored (${\delta _i} = 0$), the expectation of ${U_i}$ is related to the probability of cure and the proportion that uncured subjects for whom the event has not occurred at time $t$. Denote the expectation of $U_i^{[m]}$ as $u_i^{[m]}$,

Given the complete data $\left\{ {({y_i}, {\delta _i}, U_i^{[m]}, {{\bf{x}}_i}), i = 1, \cdots, n} \right\}$, we take the expectation of ${L^{[m]}}$ in (2.7) with respect to $u_i^{[m]}$

where

The objective function can be written as

2.2.2.   GCD-step

In the GCD-step, group coordinate descent is adopted to iteratively update ${\alpha _0}$, ${\boldsymbol{\alpha}}$, ${\beta _0}$, and ${\boldsymbol{\beta}}$. Group coordinate descent algorithms optimize the objective function with respect to a group of parameters at a time, and iteratively cycles through the parameter groups until convergence ^[25]. The intercept $\alpha _0^{[m + 1]}$ is updated by

where $\partial l_{1}^{[m]} / \partial \alpha_{0}^{[m]}=\sum_{i=1}^{n}\left(\pi_{i}^{[m]}+u_{i}^{[m]}-1\right)$, and $\partial^{2} l_{1}^{[m]} / \partial\left(\alpha_{0}^{[m]}\right)^{2}=\sum_{i=1}^{n} \pi_{i}^{[m]}\left(\pi_{i}^{[m]}-1\right)$.

For ${\boldsymbol{\alpha}} _j^{[m + 1]} \in {R^{{p_j}}}$, we can obtain Taylor's quadratic expansion of the objective function in (2.14) respect to ${\boldsymbol{\alpha}} _j^{[m + 1]}$. Referring to the fast unified algorithm for group lasso ^[32], the upper bound of the objective function can be written as

where ${{\partial {l_1}^{[m]}} \mathord{\left/ {\vphantom {{\partial {l_1}^{[m]}} {\partial {\boldsymbol{\alpha}}_j^{[m]}}}} \right. } {\partial \alpha _j^{[m]}}}{\text{ = }}\sum\nolimits_{i = 1}^n {\left({\pi _i^{[m]} + u_i^{[m]} - 1} \right)} {x_{jk}}$. Here ${\boldsymbol{V}}_{1j}^{[m]}$ is ${p_j}$ -length vector, and $M_{1j}^{[m]}$ is a constant as follows.

where ${{{\partial ^2}{l_1}^{[m]}} \mathord{\left/ {\vphantom {{{\partial ^2}{l_1}^{[m]}} {\partial \alpha _{j{k_1}}^{[m]}\partial \alpha _{jk2}^{[m]}}}} \right. } {\partial \alpha _{j{k_1}}^{[m]}\partial \alpha _{jk2}^{[m]}}} = \sum\nolimits_{i = 1}^n {\left({1 - \pi _i^{[m]}} \right)} \pi _i^{[m]}x_{jk}^2$, and $\psi (·)$ is the maximum eigenvalue function. By minimizing (2.16), we obtain ${\boldsymbol{\alpha}} _j^{[m + 1]}$ :

where ${\left(a \right)_ + } = \max \left\{ {a, 0} \right\}$.

Similarly, the intercept $\beta _0^{[m + 1]}$ is updated by

where ${{\partial l_2^{[m]}} \mathord{\left/ {\vphantom {{\partial l_2^{[m]}} {\partial \beta _0^{[m]}}}} \right. } {\partial \beta _0^{[m]}}}{\text{ = }}\sum\nolimits_{i = 1}^n {\left({{\delta _i} - u_i^{[m]}r{{\left({{y_i}\lambda _i^{[m]}} \right)}^r}} \right)}$, and ${{{{\partial ^2}l_2^{[m]}} \mathord{\left/ {\vphantom {{{\partial ^2}l_2^{[m]}} {\partial \left({\beta _0^{[m]}} \right)}}} \right. } {\partial \left({\beta _0^{[m]}} \right)}}^2}{\text{ = }} - \sum\nolimits_{i = 1}^n {\left({u_i^{[m]}{r^2}{{\left({{y_i}\lambda _i^{[m]}} \right)}^r}} \right)}$.

For ${\boldsymbol{\beta}} _j^{[m + 1]} \in {R^{{p_j}}}$, consider the optimization function

where ${{\partial l_2^{[m]}} \mathord{\left/ {\vphantom {{\partial l_2^{[m]}} {\partial {\boldsymbol{\beta}} _j^{[m]}}}} \right. } {\partial \beta _j^{[m]}}} = \sum\nolimits_{i = 1}^n {\left({{\delta _i}x_{jk}^{} - u_i^{[m]}rx_{jk}^{}{{\left({{y_i}\lambda _i^{[m]}} \right)}^r}} \right)}$.

Here ${\boldsymbol{V}}_{2j}^{[m]}$ is ${p_j}$ -length vector, and $M_{2j}^{[m]}$ is a constant as follows.

where ${{{\partial ^2}l_2^{[m]}} \mathord{\left/ {\vphantom {{{\partial ^2}l_2^{[m]}} {\partial \beta _{j{k_1}}^{[m]}\partial \beta _{j{k_2}}^{[m]}}}} \right. } {\partial \beta _{j{k_1}}^{[m]}\partial \beta _{j{k_2}}^{[m]}}}{\text{ = }} - {\sum\nolimits_{i = 1}^n {u_i^{[m]}\left({rx_{jk}^{}} \right)} ^2}{\left({{y_i}\lambda _i^{[m]}} \right)^r}$.

By minimizing (2.21), we obtain ${\boldsymbol{\beta}} _j^{[m + 1]}$ :

Regarding the parameters ${\mu _1}$, ${\mu _2}$, and $r$, Wang et al ^[33] had demonstrated that the tuning parameters selected by the Bayesian information criterion (BIC) type criterion can identify the true model consistently as long as the covariate dimension is fixed. So the parameters, ${\mu _1}$, ${\mu _2}$, and $r$, are selected by Bayesian information criterion (BIC). The EGCD algorithm is shown in Table 1.

3.   Simulation

In this section, we perform a numerical study to evaluate our method, in terms of both variable selection and estimation performance. The variable selection is assessed in terms of the (1) true positive rate (TPR), and (2) false positive rate (FPR) of variable selection. Estimation is evaluated in terms of mean square error (MSE) of coefficient estimates.

3.1.   Design and assessment

In Scenario 1, the covariates ${x_{jk}}$, $j = 1, ..., J$, $k = 1, ..., {p_j}$, are generated from a multivariate normal distribution. The correlation coefficient of covariates ${x_{j{k_m}}}$ and ${x_{j{k_n}}}$ in the same group is $\rho {\text{ = }}{0.1^{|{k_m} - {k_n}|}}$, whereas $\rho {\text{ = }}0$ when in different groups. We consider the case with discrete covariates in Scenario 2. In Scenario 2, ${x_{jk}}$ is defined as follows.

The sample size is $n = 500.$ We consider low-dimension data with $p = 40,$ and high-dimension data with $p = 200$. The censoring time is generated from a Weibull distribution with the shape parameter $r\text{ = }\left\{\text{0.25, 2.5}\right\}$. We compare the proposed method (CRGS) with three alternative methods. The alternatives are the standard cure rate model without sign consistency and variable selection penalty (CR), the cure rate model with sign consistency (CRS), and the cure rate model with group lasso penalty (CRG), respectively. For comparison, we also consider the alternatives with the logistic regression in the incident part and the Weibull distribution in the latency part. The grouping structure and coefficients of the 2 scenarios are generated as listed in Table 2.

Let $\theta \in \{ {\boldsymbol{\alpha}}, {\boldsymbol{\beta}}\}$, and $\hat{\boldsymbol{\theta }}$ be the estimation of ${\boldsymbol{\theta }}$. We evaluate variable selection performance in terms of $TPR\left({\boldsymbol{\theta }} \right)$ and $FPR\left({\boldsymbol{\theta }}\right)$ :

where

Estimates are evaluated by $MSE\left({\boldsymbol{\theta }} \right)$ :

3.2.   Results

Tables 3 and 4 summarize the mean and standard deviation in parentheses of the MSEs, TPRs, and FPRs for Scenarios 1 and 2. Both the scenarios are repeated 100 times.

As shown in Tables 3 and 4, the MSEs of methods with group lasso penalty (the proposed and CRG) are smaller than the CRS and CR methods. The estimates of CRS and CR methods have increasing MSEs with higher dimensions. The results indicate that higher dimension leads to less efficient estimation, and group lasso penalty can improve the estimation performance. Comparing the TPRs and FPRs of the proposed and CRG methods, the proposed method has better performance in terms of variable selection. Compared with alternatives, the proposed method has the lowest MSEs. Results of simulation reveal that compared with alternatives, the proposed method can improve the performance in terms of variable selection as well as estimation.

4.   Analysis of credit data in China

In this section, we apply the proposed method to credit data. The data comes from a retail business of a commercial bank in China. It contains 16 covariates of 1213 customers in a personal loan from 2014 to 2019. The primary interest is to assess the credit risk of a credit loan and find the important covariates to predict the time to default of the credit loan customers. The mean observed time is 1.38 years with a standard deviation of 0.69. Customers with missing data of annual household income are removed from our analysis. After preprocessing, the covariates and their descriptions are summarized in Table 5. By transforming the multi-level covariates, we have 24 covariates in the credit model. Censoring time ${C_i}$ is the interval between the value date and either default or the end of observation (June 1, 2019). Due to the different value dates of the loans, the censoring time vary from individual to individual. Customers whose time to event ${Y_i}$ is longer than the censoring time ${C_i}$ are censored (${\delta }_{i} = 0$). In this data, 1201 out of 1213 customers are censored.

The data are randomly divided into the training set and test set by 7:3. The training data is used for fitting the model and the test data is used to verify the performance of the fitted model. The parameters ${\mu _1}$, ${\mu _2}$, and $r$ are selected by BIC. Different from the simulations, the real coefficient is unknown for the real data. Therefore, in this section, we adopt the negative log-likelihood to evaluate the performance of the methods. The mean negative log-likelihood (standard error) of the proposed method is 25.91 (11.61), compared with 56.33 (258.73), 31.04 (27.34), and 26.17 (13.02) for the CR, CRS, and CRG method respectively. For stability, all the results are based on 100 duplicates. The results indicate that the proposed method has competitive prediction performance than alternatives.

The coefficients are estimated based on 100 duplicates. With the median estimates of coefficients, we can compute the probability of cure (non-default) ${\pi _i}\left({{{\bf{x}}_i}} \right)$ for all the customers. We dichotomize ${\pi _i}\left({{{\bf{x}}_i}} \right)$ at the median and get two different groups of customers. One group with lower ${\pi _i}\left({{{\bf{x}}_i}} \right)$ is denoted as "high risk", whereas another with higher ${\pi _i}\left({{{\bf{x}}_i}} \right)$ is denoted as "low risk". Figure 1 presents the Kaplan-Meier curves of the survival of the customers belonging to different groups. Kaplan-Meier curve describes the change of the survival probability over time. It is commonly used in survival analysis, see Rodrigues et al ^[8], and Pal ^[34]. As indicated in Figure 1, the "low risk" group has higher survival probability than the "high risk" group.

The median estimates of coefficients based on 100 duplicates are listed in Table 6. A positive coefficient ${\boldsymbol{\alpha}}$ indicates that the covariate is positively related to the probability of default, and a positive coefficient ${\boldsymbol{\beta}}$ indicates that the covariate is negatively related to default time. Both the probability of default and default time are two quite relevant credit aspects. Customers with a higher probability of default are likely to default earlier. Compared with the alternative method, the signs of the ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ of the proposed method are promoted to be more consistent, whereas many covariates such as the housing status in the CR method have an opposite effect on the probability and time to default.

The coefficient results of the proposed method reveal that loan line, early repayment, gender, housing status, annual household income, employment status, type of workplace, and occupation are important covariates for credit risk assessment. The impact of business type and education on credit is not clear.

The loan line has a positive effect. One possible explanation is that customers with better credit status are more likely to obtain a higher loan line. Customers with housing loans are more likely to default. The increasing annual household income leads to better credit status. Employed customers are less likely to default. Compared with other employment groups such as self-employed, freelance, and unemployed, the employed group has a more stable income and is less likely to default. Customers who work in a government organization and institution have better credit status. Customers who are managers or commercial and service workers are less likely to default. Customers with early payment records tend to maintain good credit records and are less likely to default. Compared with women, men are more likely to default. This is consistent with the results of ^[35] and the personality characteristics of men's risk preference ^[36].

5.   Conclusions

The cure rate model is commonly used when the data has long-term survivors. The model is composed of two parts. The incident part describes the probability of cure and the latency part describes the survival function of the uncured group. The drawback of the standard cure rate model is that it assumes that there are no direct constraints between the coefficients corresponding to the same covariates in the two parts of the model. This may lead to conflicting results of covariate effects on the probability of cure and conditional survival of the uncured group. In fact, the two parts of the model describe quite related aspects. It is desirable that there may be some connections between coefficients corresponding to the same covariates. Existing works have considered joint distribution or structural effect of the two sets of covariates, which is too strict.

In this paper, we consider a more flexible cure rate model that allows the two sets of covariates can be in different distributions and magnitudes. In the proposed method, a sign consistency cure rate model is proposed to promote the similarity in the sign of coefficients in the two model parts to improve interpretability. In addition, we also impose a group lasso penalty for variable selection. Simulation results show that compared with alternatives, the proposed method has better performance in terms of variable selection and estimation. An analysis of credit data in China illustrates that the proposed method can improve prediction performance as well as interpretability.

Acknowledgments

We are grateful to the reviewers and the editor for their helpful comments and suggestions. This work was supported by the National Office for Philosophy and Social Sciences of China under Grant 20 & ZD137.

Conflict of interest

The authors declare there is no conflict of interest.