Group MCP for Cox Models with Time-Varying Coefficients

1 Introduction

Cox’s proportional hazard model, which was introduced by Cox^[1,2] has become the most common method to analyze time to event data. In classical setting where the number of covariates of covariates p is fixed and the sample size n is large, many variable selection techniques for linear regression models have been extended to Cox model such as best-subset selection, stepwise selection, asymptotic procedures based on score tests, Wald tests, other approximate chi-square testing procedures, bootstrap procedures^[3] and Bayesian variable selection ^[4,5]. However, the theoretical properties of these methods are generally unknown^[6]. Recently, a number of regu larization methods have been proposed for Cox model with time-independent coefficients. The lasso, proposed for simultaneous coefficient estimation and variable estimation^[7] , was extended to Cox models^[8]. Fan and Li ^[6] extended the smoothly clipped absolute deviation (SCAD) penalty^[9] to Cox models. Zhang and Lu^[10] proposed the adaptive lasso with the adaptively weighted penalty L₁ on each coefficient. For Cox models with time-varying coefficients^[11], pro viding more flexibility and challenge than the parametric linear models, the literature appears to be limited. Lin and Zhang^[12] proposed a component selection and smoothing operator in the frame of smoothing spline analysis of variance, and Leng and Zhang^[13] extended this approach to varying coefficients Cox models. Huang, Horowitz, and Wei^[14] proposed the adaptive group lasso method which approximated additive components with B-spline expansions and selected nonzero components by selecting the groups coefficients in the expansion for nonparametric additive models. This adaptive group lasso method, which not only selects important variables but also selects between time-independent and time-varying specifications, was later extended for Cox models with time-varying coefficients by Yan and Huang^[15].

Recent technological advances have made it possible to collect a huge amount of covari ate information. However, due to the high dimensional space of the predictors, the standard maximum Cox partial likelihood method cannot be applied directly to obtain the parameter estimator. In the high-dimensional and low-sample size setting for Cox models, Li and Luan^[16] investigated the L₂ penalized estimation which used the kernel tricks to reduce the computation to just involve the inversion of matrix of the sample size. Li and Gui^[17] proposed to use the LARS algorithm to obtain the solution for the Cox model with L₁ penalty. Tibshirani^[18] used the univariate shrinkage method which assumed the independence of the covariates in each risk set and the partial likelihood factors into a product, and then the variables are entered into the model based on the size of their Cox score statistics. Bradic, Fan and Jiang^[19] proposed the penalized partial likelihood approach employing a class of folded-concave penalties to the Cox parametric relative risk model for non-polynomial dimensional date. Fan, Feng and Wu^[20] extended the sure independence screening (SIS) and iterative SIS (ISIS) for Cox model and have shown its encouraging performance.

In this paper, we study the Cox models with time-varying coefficients, and then apply the group minimax concave penalty (group MCP)^[21] method that uses smooth function for coefficients as well as maximizes the log partial likelihood with nonconvex penalties on regression coefficients^[22] Although the variable selection methods involving nonconvex penalty functions, such as SCAD^[9] and MCP^[23], may introduce numerical challenges in fitting these models, they are worthwhile alternatives to the lasso in many applications and have attractive theoretical properties^[24]. Under some appropriate conditions, Yang, Huang and Zhou^[21] showed that the group MCP has the oracle selection property that it correctly selects important variables with probability converging to one. Furthermore, these concave can be applied to high-dimensional problems. In this paper, each time-varying coefficient is expanded over a B-spline basis, and then characterized by a set of bases coefficients that are further treated as a group. We select significant variables by applying the group MCP approach over these groups of bases coefficients^[25] with a fast, iterative group shooting algorithm.

The rest of the paper is organized as follows. Section 2 proposes a group MCP method with penalized partial likelihood based on B-splines. Computational details of the proposed model selection procedure are described in Section 3. The variance estimation for estimated coefficients are also presented in Section 3. Section 4 gives numerical comparisons among the group MCP and group Lasso. Section 5 gives the discussion.

2 Group MCP with B-Splines

Suppose a random sample of n individuals is chosen. Let Ti∗ and C_i be the failure time and censoring time of subject i, i = 1, 2, · · · , n, respectively. Define observed time Ti=min{Ti∗,Ci} and censoring indicators δi=I(Ti∗≤Ci)i = 1,2, ··· , n. Let X_i = (X_i1, X_i2, · · · , X_ip)^T be the vector of covariate for subject i. Assume that Ti∗ and C_i are conditionally independent given X_i, and the censoring mechanism is uninformative. Our data set are independently and identically random samples that consists of the triplets {T_i, δ_i, X_i}, i = 1, 2, · · · , n.

The Cox model with time-varying coefficients is

where h₀is an unspecified baseline function, and β(t) is p × 1 vector of time-varying coefficients. Let (B₁(t), B₂(t), · · · , B_q(t))^Tbe a set of B-spline basis of q degrees of freedom without intercept on a predetermined time interval [0, τ]. This B-spline basis can be obtained from function bs in package splines from base R^[26]. In our simulation and analysis, we used function bs with quadratic splines with q degrees of freedom (df = q), with equally spaced interior knots, and intercept = FALSE. Further assume that β(t)is expanded by the B-spline basis, that is, β(t) = ΘF(t), where F(t)= (B₁(t), B₂(t), · · · , B_q(t))^Tand Θ is a p × q matrix of parameters to be estimated. Therefore, each time-varying coefficient β_j(t) = Θ_jF(t), j = 1, 2, ··· , p, is determined by the jthrow of parameter matrix Θ.

For simplicity, assume that there are no ties in the observed failure times. Present of ties we may use the technique in Breslow^[27]. Let θ = vech(Θ), the vectorization of Θ by row. Let R_i = {j : T_j < T_i}be the risk set at time T_i. The log partial likelihood function is then given

by

Following Tibshirani^[8] and Fan and Li^[6], we estimate θ by minimizing the negative penalized log partial likelihood function

where Ρ (θ; λ, γ) is a any penalty function that penalizes coefficient estimates in groups with a tuning penalty parameter λ^[25] and a regularization parameter γ.

To select significant variables, we should partition Θ first. A straightforward and useful way is to put each row in Θ into a single group, which leads to penalty function

The penalty function ρ is applied to a penalty function for the group of parameters. In generally, any penalty function that works well for individual variable selection can be considered to be extended to this group penalty function. For instance, the resulting criterion is the group lasso^[25] if we use the L₂ penalty, while it comes to be group MCP^[21] if MCP penalty function is used.

The group lasso method^[25] was defined as the value of that minimizes

The group MCP was defined as

for λ < 0 and γ ≥ 1, and the (a)₊ = aI{a ≥ 0} is the nonnegative part of a. The regularization parameter γ controls the degree of concavity, meaning that, a large value of γ makes ρ less concave^[24]. Detail discussion of MCP can be found by Zhang^[23]. Note that the group lasso can be considered a special case of group MCP by taking γ → ∞. This penalty treats Θ_j as a whole group and each covariate coefficient is penalized via a single penalty. When Θ_j are nonzero, the corresponding covariate coefficient β_j is nonzero while β_j is zero if Θ_j are zero.

3 Computation

A convex objective is desirable because convexity ensures that the algorithm converges to the unique global minimum once it converges to a critical point of the objective function^[24] The penalty function of MCP is concave, which introduces a challenge that the algorithm are not guaranteed to converge a global minimum in general. However, Breheny and Huang^[24] proposed a convexity diagnostics measure and showed that it is possible for the objection function to be convex even though it contains a nonconvex penalty component. For example, the MCP object function is convex if γ > 1/c_*, where c_* is the minimum eigenvalue of n^-1X^TX, and then the algorithm can converge to the global minimum.

4 Numerical Studies

In this section, we compare the performance of our group MCP method and group lasso method under the Cox model with time-varying coefficients.

In our simulation design, we consider three factors: Sample size (100 and 150), number of covariates (10 and 20) and the censoring rates (20% and 30%).

Event times are generated from varying-coefficient Cox model with covariate vector X and coefficient β(t), and the nonzero components are β₂(t) = {l + cos(πt)}I(0 < t < 1), β₃(t)= {1 + 0.5sin(πt)} and β₈(t) = l, for t ∊ (0,2), the baseline hazard function is λ(t) = exp{sin(πt/2)}. That is, the second and the third ones have time-varying coefficients, the eighth one is a constant coefficient out of 10 or 20 covariates, and all the rest have coefficient zero. Note that β₂(t) diminished to zero when t ≥ 1, which makes model selection and estimation harder. Covariate vector X is generated from a multivariate normal distribution whose marginals are all N(0,0.5), and whose pairwise correlation coefficients are ρ^|j–k| for pair (j, k), we take ρ = 0.5, which mean that covariates have moderate correlation. Censoring times are generated from a mixture of uniform distribution over (0, 2) and a point mass at 2, with the mixing probability calibrated to yield desire censoring percentage c_p, this way was used by Yan and Huang^[15].

For each scenario, 100 datasets are generated. Given a simulated dataset, we use quadratic B-splines with 5 degrees of freedom, with equally spaced knots in time window (0, 2), for each covariate coefficient. This gives two equal distant interior knots in (0,2). Generalized cross validation is used to estimate the tuning parameter γ in the group lasso and group MCP method as described above.

To measure prediction accuracy, we report the average number of groups selected (NG), and the average mean squared error (MSE) over 100 runs in each scenario. The “correct” NGs are all 3 for group MCP and group lasso, respectively. Following Tibshirani^[8], the MSE is calculated as {β̂(t) − β(t)}^TV{β̂(t)— β(t)}, where V is the population covariance matrix of the covariates, and we report the average of pointwise MSE.

In our simulation study, both group MCP method and group lasso method work analogously in important covariates X₂, X₃and X₈. Due to a diminishing effect in covariate X₂, it is difficult to be selected, especially in low-sample size. It is easily found that group lasso select more frequently than group MCP method when choosing unimportant covariates. However, the group lasso method has smaller MSE relative to group MCP method.

Several observations such as the variable selection results, NG and MSE from 100 runs with 10 covariates can be obtained from Table 1. Covariates X₃and X₈are selected most of the times, and covariate X₂is selected less. When the sample size increases, both methods seem to perform better as expected. For example, in the scenario of sample size n = 100 and censoring rate c_p = 30%, the covariate X₂is selected almost half of the times for both methods. As sample size increase to 150, the covariate X₂is selected most of the times with smaller MSE and less overselection. The MSE of group MCP method in the scenario of the n = 100, c_p = 30% with 10 covariates is 0.526, and it is reduced to 0.349 in the scenario of n = 150, c_p = 30%, which means that the group MCP method selects important variables more accurately as the sample size increasing. On the other hand, the unimportant covariates are also selected less as a whole. As censoring rate increase, both methods appear to perform worse, especially in the small size setting. In the scenario of n = 100 and c_p = 20%, both methods select covariate X₂about three/four of the times, and they select it only half of the times for c_p = 30%. The covariates X₃and X₈ are also select less at the same time.

It is worth noting that we observed such a phenomenon: On the one hand the MSE in small sample size with a high censoring rate is significantly higher than that in other settings, on the other hand the group MCP has a smaller MSE than the group lasso method. For example, with regard to mean squared error, although the group lasso perform better than the group MCP in most situations, in the scenario of n = 100, c_p = 30%, the MSE of the group lasso method is 0.554, While the MSE of the group MCP method is 0.526. That is, the group MCP method performs better than group lasso method in such low sample size with high censoring rate setting. We conduct more simulation in higher censoring rate setting to study the performance of these two methods, Table 2 summarizes the results. Both scenarios indicate that the group MCP method and group lasso method done poorly, but the group MCP method performs better than the competing method. When censoring rate increases, it is difficult to reflect the complete information in such setting. Moreover, the number of parameters to be estimated in Θ is large, so a large sample size is necessary.

Figure 1 shows the 100 estimated coefficient curves overlaid with the true curves for the scenario of n = 100 with 10 covariates, c_p = 20% and c_p = 30%, respectively. To compare the group MCP method and group lasso method in the performance of recovering of the coefficient, we also plot the average of the 100 estimated coefficients curves and their pointwise 95% confidence intervals constructed using the method in Section 3.2. In these scenarios, both method perform well and the estimated curves are closer to the true curves for all three coefficients, especially for coefficients β₃ and β₈. For coefficient β₂, the estimator curves are also close to the true curves in their nonzero part, and when they diminishes to zero, the estimator curves are slightly diverge. It is not a surprise because it remains zero afterward which makes model estimation harder. The group MCP method perform well that the estimated curves for coefficient β₂ are tighter around the true curves compared with the group lasso method. As we observe in Figure 1, the standard errors seem to underestimate the true variation. Overall, the group MCP can fit better than the group Lasso.

Table 3 shows the model selection results from 100 runs with 20 covariates and the conclusion are similar to those in Table 1. However, in low sample size with a high censoring rate setting, the group MCP method seems to select more covariate X₂and less covariates X₃and X₈ compared to the same setting with 10 covariates.

Figure 2 shows the 100 estimated coefficient curves overlaid with the true curves for the scenario of n = 100 with 20 covariates, c_p = 20% and c_p = 30%, respectively. The conclusion are also similar to those in Figure 1 that the group MCP can fit better than the group Lasso.

Figure 1

Estimated curves (gray) of the three nonzero coefficients from 100 runs in the scenario of sample size n = 100 with 10 covariates, c_p = 20% and c_p = 30%, respectively. The dark lines are the true curves. The dashed lines are the average of 100 estimated coefficients curves. The dotted lines are the pointwise 95% confidence intervals

Figure 2

Estimated curves (gray) of the three nonzero coefficients from 100 runs in the scenario of sample size n = 100 with 20 covariates, c_p = 20% and c_p = 30%, respectively. The dark lines are the true curves. The dashed lines are the average of 100 estimated coefficients curves. The dotted lines are the pointwise 95% confidence intervals

These simulation results indicate that the group MCP and group lasso method have good selection and estimation performance in the simulated model. Overall, the performances of these two methods are similar. However, the group MCP method selects fewer unimportant covariates than the group lasso method, yet the MSE of the group Lasso is smaller. This was also observed in the study of Yang, Huang and Zhou^[21]. In the low sample size with high censoring rate setting, the group MCP have better selection and estimation performance in the simulation study. On the whole, the group MCP can improve the selection and estimation results over the group Lasso.

5 Discussion

In this article, we studied the model selection and estimation of the group Lasso and group MCP in Cox models with time-varying coefficients on right-censored failure times. In our study, our nonparametric coefficients are fitted with smooth functions which are expanded using B-spline basis. Also, a set of bases coefficients further are treated as a group and applied the group MCP and group lasso approach. Our simulation studies indicate that the group MCP outperforms the group Lasso in terms of selection and estimation.

Although the selection consistency and asymptotic oracle properties of group MCP have been proofed by Yang, Huang and Zhou^[21] in high-dimensional varying coefficient models, a rigorous proof of our group MCP estimator’s properties in Cox model is not given in this article. This is a challenging problem when extend this approach to Cox models with time-varying coefficients. Research to the theoretical properties of group MCP to such model is necessary.

The methods in this article raise several questions. As shown in Figure 1 and Figure 2, although the coefficient β₈ is a time-independent coefficient, it is selected to have a time-varying effect. This may be explained by the fact that we put each row in Θ into a single group so that both covariates coefficients are treated as time varying coefficients. Yan and Huang^[15] penalized a time-independent part and a time-varying part separately for each coefficient, and then it can select significant variables and the temporal dynamics of their effects. This is an efficient and desirable way that can be considered to solve our problem.

In principle, the group MCP method can be applied in high-dimensional Cox models with time-varying or time-independent coefficients. However, the iterative group shooting algorithm involves matrix calculation and the computational cost of this algorithm in high-dimensional setting is relatively high. Therefore, more work is needed to study such as the theoretical properties and computational procedures in these more complicated models.