coxphMIC: An R Package for Sparse Estimation of Cox Proportional Hazards Models via Approximated Information Criteria

In this paper, we describe an R package named coxphMIC , which implements the sparse estimation method for Cox proportional hazards models via approximated information criterion (Su et al., 2016). The developed methodology is named MIC which stands for “Minimizing approximated Information Criteria". A reparameterization step is introduced to enforce sparsity while at the same time keeping the objective function smooth. As a result, MIC is computationally fast with a superior performance in sparse estimation. Furthermore, the reparameterization tactic yields an additional advantage in terms of circumventing post-selection inference (Leeb and Pötscher, 2005). The MIC method and its R implementation are introduced and illustrated with the PBC data.


Introduction
Time to event (survival time) is often a primary outcome of interest in many research areas, especially in medical research such as time that takes to respond to a particular therapy, time to death, remission, or relapse.Survival times are typically right skewed and subject to censoring due to study termination, loss of follow ups, or withdrawals.Moreover, covariates may vary by time.
Cox Proportional Hazards (PH) model (Cox, 1972) is commonly used to model survival data.Given a typical survival data set that consists of {(T i , δ i , z i ) : i = 1, . . ., n}, where T i is the observed event time, δ i is the 0-1 binary censoring indicator, and z i ∈ R p is the covariate vector associated with the i-th subject, the Cox PH model formulates the hazard function h(t|z i ) for the ith subject as h(t|z i ) = h 0 (t) exp(β T z i ), where z i ∈ R p denotes the p-dimensional covariate vector associated with subject i, β = (β j ) ∈ R p is the unknown regression parameter vector, and h 0 (t) is the unspecified baseline hazard function.The vector of β can be estimated by maximizing the partial log-likelihood (Cox, 1975), which is given by Let β denote the resultant maximum partial likelihood estimator (MPLE).
Since the true β is often sparse, we need to look for methods that identify the zero components in β and at the same time estimate the nonzero ones.Best subset selection (BSS) and regularization are among two major algorithms used in survival analyses for variable selection.Both are derived from a penalized partial likelihood.Let pen(β) and λ denote the penalty function and penalty parameter, respectively.The general objective function in both of the techniques is as follows: In BSS, the penalty function is set to pen(β) = ∑ p j=1 I{β j = 0} (number of nonzero coefficients), and the penalty parameter is fixed as λ = 2 for AIC (Akaike, 1974) or λ = ln(n 0 ), where n 0 is the total number of uncensored failures, with a slight modification of BIC (Vollinsk and Raftery, 2000).In regularization, the penalty function is set to pen(β) = ∑ p j=1 |β j |, and the penalty parameter is not fixed and is appropriately chosen.The sparse estimation is reformulated into a continuous convex optimization problem.The optimization of the two techniques is a two-step process.In BSS, one needs to fit every model with the maximum partial likelihood method and then compare the fitted models according to an information criterion such as AIC (Akaike, 1974) or BIC (Schwarz, 1978).This makes the BSS infeasible for moderately large p.In regularization, one need to solve the objective function for every fixed positive value of λ to obtain a regularization path { β(λ) : λ > 0}, and then select the best λ according to an information criterion such as AIC or BIC along the path.Since such a search is only along the regularization path (a one-dimensional curve in R p ), the search space is Prepared in the R Journal style June, 2016 much reduced and hence, it may not perform as well as the estimator obtained with BSS, if AIC or BIC is used as the yardstick.Beside the computational burden, both methods face the post-selection inference challenge.A new technique is developed by Su et al. (2016) on the basis of Su (2015) for conducting sparse estimation of Cox PH models to help address the aforementioned deficiencies.

The MIC method
A new method, named MIC for "Minimizing approximated Information Criteria", is developed to conduct sparse estimation of Cox PH models.MIC borrows strength from both BSS and regularization.The main issue with BSS is the indicator function, I(β = 0), involved in the ℓ 0 penalty function, leading to a discrete optimization problem.To overcome this difficulty, MIC proposes to approximate the indicator function by a continuous or smooth unit dent function.One reasonable approximation is the hyperbolic tangent function given by where a is a nonnegative scale parameter that controls the sharpness of the approximation.As shown in Figure 1(a), w(β) provides a smooth approximation to the discrete function I{β = 0}.However, the curve does not have zero as a singular point.If we estimate β by minimizing −2 l(β) + ln(n 0 ) ∑ p j=1 w(β j ), we will not obtain sparse estimates.To enforce sparsity, MIC devises a reparameterization step.The reparameterization is based on the decomposition β = β I{β = 0}.Set γ = β and approximate I{β = 0} with w(γ) = tanh(aγ 2 ).This leads to a reparameterization of Prepared in the R Journal style June, 2016 β = γw(γ).The objective function in MIC is given by where the penalty parameter λ 0 is fixed as ln(n 0 ) for BIC (Vollinsk and Raftery, 2000) and matrix W is p × p diagonal with diagonal elements w j = w(γ j ) and hence trace tr(W) = ∑ p j=1 tanh(aγ 2 j ).With this notation, it follows that β = Wγ.
The above reparameterization offers several important conveniences: 1. Sparsity now becomes achievable in estimating β.The penalty w(γ) as a function of β = γw(γ) = γ tanh(γ) is a unit dent function that is smooth everywhere except at β = 0, as shown in Figure 1(c).This is a necessary condition to ensure sparsity as indicated by Fan and Li (2002).On this basis, the oracle properties of the MIC estimator β obtained by minimizing Q n (β) in ( 1) can be established under regularity conditions.The asymptotic results entails a n = O(n).For this reason, we fix a n = n 0 , the number of non-censored failures.In practice, the empirical performance of MIC is large stable with respect to the choice of a, as demonstrated in Su (2015).
Thus simply fixing a at a reasonably large value (say, a ≥ 10) could do as well practically.
2. In terms of practical optimization, it is preferable to consider γ as the decision vector.Namely, we minimize Q n (γ) with respect to γ by treating it as a function of γ.Let γ be the resultant MIC estimator of One immediate advantage of doing so is that Q n (γ) is smooth in γ and hence many optimization routines can be applied directly.Since no selection of tuning parameters is involved, MIC is computationally efficient.
3. One consequence of post-selection inference is that no standard error formula is available for zero estimates of β j .As depicted in Figure 1(b), β j and γ j have a one-to-one correspondence with β = 0 iff γ = 0.This motivates us to test H 0 : β j = 0 by equivalently testing H 0 : γ j = 0.The MIC estimator γ can be viewed as an M-estimator with smooth objective function Q n (γ) and hence standard arguments can be used to make inference.

Implementation in R
The R package coxphMIC implements MIC on the basis of R package survival (Therneau and Grambsch, 2000) and is hosted at CRAN.Type the following command in R console in order to install the package:

> install.packages("coxphMIC")
To summarize, MIC can be simply formulated as the following optimization problem Owing to the non-convex nature, a global optimization method is helpful in solving (3).While other R routines (Mullen, 2014) are available, we have found using the SANN method combined with the BFGS method in R function optim() is fast and quite effective.The simulated annealing (SA) implemented by SANN helps locate a nearly minimum point globally.Then the quasi-Newton BFGS method makes sure that the algorithm stops at a critical point.
There are two functions included in the coxphMIC package: an internal function LoglikPen() that computes the partial log-likelihood and a wrapper function coxphMIC() that does the MIC sparse estimation.The function coxphMIC() has the following usage:  We briefly explain some of the important options.The formula argument is a formula object similar to that in survival, with the response on the left of the ~operator being a survival object as returned by the Surv function, and the terms on the right being predictors.The arguments method.beta0,beta0, and theta0 pertains to the initial starting values.By default, the maximum partial likelihood estimator with the option MPLE is used.Otherwise, one can use the ridge estimator with option ridge.The theta0 corresponds to the tuning parameter in ridge estimation.User defined starting values can also be used such as β = γ = 0 by specifying beta0.By default, the approximated BIC (Vollinsk and Raftery, 2000) is recommended.However, one can use AIC (Akaike, 1974).Alternatively, user-specified penalty is allowed by specifying lambda0.The default value for a is n 0 .The option maxit.globalallows for specification of the maximal iteration steps in SANN while maxit.localspecifies the maximal iteration steps for BFGS.MIC computes the standard errors (SE) for both β and γ.For β, the SE computation is only applicable for its nonzero components.The option maxit.global asks whether the user wants to output the confidence intervals for γ j at the confidence level specified by conf.level(with 95% as default).
The output of Function coxphMIC() is an object of S3 class coxphMIC, which is essentially a list of detailed objects that can be used for other purposes.In particular, the item result presents the most important results, where one can see the selected model and inference based on testing γ.Two generic functions, print and plot, are made available for exploring a coxphMIC object.

Other R packages for variable selection in Cox PH models
Several other R packages are available for variable selection of Cox PH models.The best subset selection (BSS) is available in the R Package glmulti (Calcagno and de Mazancourt, 2010) with AIC only, but it is very slow owing to the intensive computation involved.For large p, a stepwise selection procedure could be used as a surrogate.LASSO (Tibshirani, 1997) can be computed via R Package glmnet (Friedman et al., 2010).Zhang and Lu (2007) have made their R codes for implementing ALASSO for Cox models available at http://www4.stat.ncsu.edu/~hzhang/paper/cox_new.tar.But the program was written without resorting well to available R routines and it takes an unnecessarily long running time.One alternative way to compute ALASSO is first transform the design matrix Z := Z diag(| β|) so that LASSO could be applied and then transform the resultant estimates back.SCAD for Cox PH models (Fan and Li, 2002;Fan et al., 2010) can be computed with an earlier version of the R package SIS (Saldana and Feng, 2016), but it is no longer available in its current version.One is referred to Table 1, which is presented as Table B1 in Su et al. (2016), for a comparison study of these above-mentioned methods.MIC clearly stands out as the top or among-the-top performer in both sparse estimation and computing time.

Examples
We illustrate the usage of coxphMIC via an analysis of the PBC (primary biliary cirrhosis) data, available from the survival package (Therneau and Grambsch, 2000).

Data preparation
To proceed, some minor data preparation is needed.First of all, we want to make sure that the censoring indicator is 0-1 binary.
> library(survival); data(pbc); > dat <-pbc; dim(dat); [1] 418 20 > dat$status <-ifelse(pbc$status == 2, 1, 0) Next, we explicitly created dummy variable for categorical variables.The factor() function could be used instead.Also, grouped sparsity could be used to handle these dummy variables so that they are either all selected or all excluded.We plan to explore this possibility in future research.
> dat$sex <-ifelse(pbc$sex == "f", 1, 0) Another necessary step is to handle missing values.This current version does not automatically treat missings.Here, the listwise deletion is used so that only the 276 subjects with complete records are used for further analysis.The data set now contains 20 variables.Except id, time, and status, there are a total of 17 predictors.

MIC starting with MPLE
To apply coxphMIC, one simply proceeds in the usual way of using coxph formula.By default, all predictors are standardized; the approximated BIC (λ 0 = ln(n 0 ) is used with a = n 0 ; and the MPLE is used as the starting point.
> fit.mic <-coxphMIC(formula = Surv(time, status)~.-id,data = dat, CI.gamma = FALSE) > names(fit.mic) [1] "opt.global" "opt.local""min.Q" "gamma" "beta" "VCOV.gamma"[7] "se.gamma" "se.beta" "BIC" "result" "call" The output of coxphMIC contains the minimized Q n value, the final estimates of γ and β, the variancecovariance matrix and SE for γ, SE for nonzero β, BIC value for the final model, and a summary table result.In order for the user to be able to inspect the convergence and other detailed info of the optimization algorithms, we also output two objects The above results are presented as Table 4 in Su et al. (2016).In this example, MIC started with MPLE given by the first column named beta0.Columns 2-5 present estimation of γ and the hypothesis testing results on H 0 : γ j = 0.The estimates of β are given in the last two columns.It can be seen that eight variables are selected in the final model, which are age, edema, bili, albumin, copper, ast, protime, and stage.
The plot function provides error bar plots based on the MIC estimator of both β and the reparameterized γ : > plot(fit.mic,conf.level= 0.95) as shown in Figure 2. Essentially, the 95% confidence intervals (CI) are plotted.One can modify the confidence level with the conf.leveloption.To compare two plots conveniently, they are made with the same range on the vertical y-axis.Note that CI is not available for any zero β j estimate in Panel (b), which corresponds to an unselected variable.Those selected variables are highlighted in green color in Panel.Prepared in the R Journal style June, 2016

Multiple starting points
Trying out multiple starting point is a common strategy in facing global optimization problems.We may consider starting with the 0 vector, which corresponds to the null model.Having beta0 = 0 is actually the default option if method.beta0 is neither 'MPLE' nor 'ridge' and a specific value for beta0 is not given, i.e., setting beta0 = NULL.
Concerning sparse estimation, the vectors with 0/+1/-1 values obtained by applying a threshold to the MPLE | β| could be reasonable choices for the starting point too, i.e., where c 0 > 0 is a threshold close to 0. For example, setting c = 0.06 yields 0.3881 1 In the above example, we applied a threshold of 0.06 to the MPLE to obtain a 0/+1/-1 valued vector.
To start with this user-supplied starting point, one proceeds as follows.

Different a values
We may consider obtaining the regularization path with respect to a.According to asymptotic results, a = O(n) is desirable and the recommended value is a = n 0 the number of uncensored deaths, which is n 0 = 111 in the PBC data under study.We try out a spread of a values that range from 10 to 200, as prescribed by the R object A0.

Figure 2 :
Figure 2: Error bar plots for MIC estimates of γ in (a) and β in (b).The 95% confidence intervals (CI) are plotted.The selected variables are highlighted in green in Panel (b).

>Figure 3 :
Figure 3: The regularization path with respect to a with the PBC data.The a value varies from 10 to 200.The recommended value is a = n 0 = 111.

Table 1 :
Comparison of computation time: CPU time (in seconds) averaged over three runs.
opt.global and opt.local, which result from the global (SANN by default) and local optimization (BFGS by default) algorithms.The output fit.mic is a S3 object of coxphMIC class.Two generic functions, print and plot, are available.The print function provides a summary table as below: