# A Fast Divide-and-Conquer Sparse Cox Regression

###### Abstract

We propose a computationally and statistically efficient divide-and-conquer (DAC) algorithm to fit sparse Cox regression to massive datasets where the sample size is exceedingly large and the covariate dimension is not small but . The proposed algorithm achieves computational efficiency through a one-step linear approximation followed by a least square approximation to the partial likelihood (PL). These sequences of linearization enable us to maximize the PL with only a small subset and perform penalized estimation via a fast approximation to the PL. The algorithm is applicable for the analysis of both time-independent and time-dependent survival data. Simulations suggest that the proposed DAC algorithm substantially outperforms the full sample-based estimators and the existing DAC algorithm with respect to the computational speed, while it achieves similar statistical efficiency as the full sample-based estimators. The proposed algorithm was applied to an extraordinarily large time-independent survival dataset and an extraordinarily large time-dependent survival dataset for the prediction of heart failure-specific readmission within 30 days among Medicare heart failure patients.Divide-and-conquer; shrinkage estimation; variable selection; Cox proportional hazards model; least square approximation.

section1[Introduction]Introduction Large datasets derived from health insurance claims and electronic health records are becoming increasingly available for health care and medical research. These datasets serve as valuable sources for the development of risk prediction models, which are the key components of precision medicine. Fitting risk prediction models to a dataset with a massive sample size (), however, is computationally challenging, especially when the number of candidate predictors () is also large and yet only a small subset of the predictors are informative. In such a setting, it is highly desirable to fit a sparse regression model to simultaneously remove non-informative predictors and estimate the effects of the informative predictors. When the outcome of interest is time-to-event and is subject to censoring, one may obtain a sparse risk prediction model by fitting a regularized Cox proportional hazards model (Cox, 1972) with penalty functions such as the adaptive least absolute shrinkage and selection operator (LASSO) penalty (Zhang and Lu, 2007).

When is extraordinarily large, directly fitting an adaptive LASSO penalized Cox model to such a dataset is not computationally feasible. To overcome the computational difficulty, one may employ the divide-and-conquer (DAC) strategy, which typically divides the full sample into subsets, solves the optimization problem using each subset, and combines the subset-specific estimates into a combined estimate. Various DAC algorithms have been proposed to fit penalized regression models. For example, Chen and Xie (2014) proposed a DAC algorithm to fit penalized generalized linear models (GLM). The algorithm obtains a sparse GLM estimate for each subset and then combines subset-specific estimates by majority voting and averaging. Tang and others (2016) proposed an alternative DAC algorithm to fit GLM with an extremely large and a large by combining de-biased LASSO estimates from each subset. While both algorithms are effective in reducing the computation burden compared to fitting a penalized regression model to the full data, they remain computationally intensive as penalized estimation procedures will be required. In addition, the DAC strategy has not been extended to the survival data analysis.

In this paper, we propose a novel DAC algorithm using sequences of linearization, denoted by , to fit adaptive LASSO penalized Cox proportional hazards models, which can further reduce the computation burden compared to the existing DAC algorithms. starts with obtaining an estimator that maximizes the partial likelihood (PL) of a subset of the full data, which is then updated using all subsets via one-step approximations. The updated estimator serves as a -consistent initial estimator for the adaptive LASSO problem and approximates the full sample-based maximum PL estimator. Subsequently, we obtain the final adaptive LASSO estimator based on an objective function applying the least square approximation (LSA) to the PL as in Wang and Leng (2007). The LSA allows us to fit the adaptive LASSO using a pseudo likelihood based on a sample of size . The penalized regression is only fit once in the proposed algorithm and the improvement in computation cost is substantial if . Our proposed algorithm can also accommodate time-dependent covariates.

The rest of the paper is organized as follows. We detail the algorithm in section A Fast Divide-and-Conquer Sparse Cox Regression. In section A Fast Divide-and-Conquer Sparse Cox Regression, we present simulation results demonstrating the superiority of compared to the existing methods when covariates are time-independent and when some covariates are time-dependent. In section A Fast Divide-and-Conquer Sparse Cox Regression, we employ the algorithm to develop risk prediction models for 30-day readmission after an index heart failure hospitalization with data from over 10 million Medicare patients by fitting regularized Cox models with (i) time-independent covariates and (ii) time-independent covariates and time-dependent environmental covariates. We conclude with some discussions in section A Fast Divide-and-Conquer Sparse Cox Regression.

section1[Methods]Methods

subsection2[Notation and Settings]Notation and Settings Let denote the survival time and denote the vector of bounded and potentially time-dependent covariates. Due to censoring, for , we only observe , where , , and is the censoring time assumed to be independent of given . Suppose the data for analysis consist of subjects with independent realizations of , denoted by , where we assume that .

We denote the index set for the full data by . For all DAC algorithms discussed in this paper, we randomly partition into subsets with the -th subset denoted by . Without loss of generality, we assume that is an integer and that the index set for the subset is . For any index set , we denote the size of by with if and if . Throughout we assume that such that and .

To develop a risk prediction model for based on , we consider the Cox model,

(0.0) |

where is the conditional hazard function and is the baseline hazard function. Our goal is to develop a computationally and statistically efficient procedure to estimate using data in under the assumption that is sparse with the size of the active set much smaller than .

When is not extraordinarily large, we may obtain an efficient estimate, denoted by , based on the adaptive LASSO penalized PL likelihood estimator as proposed in Zhang and Lu (2007). Specifically,

(0.0) |

where for any index set ,

(0.0) |

is an initial -consistent estimator of model ( A Fast Divide-and-Conquer Sparse Cox Regression), is a tuning parameter, and . A simple choice of is , where for any set ,

Following the arguments given in Zhang and Lu (2007), when and , we can show that achieves the variable selection consistency, i.e. the estimated active set satisfies and that the oracle property holds, i.e.

where for any set , if is a vector and if is a matrix,

, , , , , , , and for any vector .

When is not too large, multiple algorithms are available to solve ( A Fast Divide-and-Conquer Sparse Cox Regression) with time-independent covariates, including a gradient descent algorithm (Simon and others, 2011), a least angle regression (LARS)-like algorithm (Park and Hastie, 2007), a combination gradient descent-Newton Raphson method (Goeman, 2010), and a modified shooting algorithm (Zhang and Lu, 2007). Unfortunately, when is extraordinarily large, none of the existing algorithms for ( A Fast Divide-and-Conquer Sparse Cox Regression) will be computationally feasible. While these algorithms may be extended to fit sparse Cox models with time-dependent covariates, the computation is even more demanding since each subject may contribute multiple observations in the fitting.

subsection2[The Algorithm]The Algorithm The goal of this paper is to develop an estimator that achieves the same asymptotic efficiency as but can be computed very efficiently.

Our proposed algorithm, , for attaining such a property is motivated by the LSA proposed in Wang and Leng (2007), with the LSA applied to the full sample-based PL. Specifically, it is not difficult to show that is asymptotically equivalent to , where

That is, will also achieve the variable selection consistency as and has the same limiting distribution as that of . This suggests that an estimator can recover the distribution of if we can construct an accurate DAC approximations to and . To this end, we propose a linearization-based DAC estimator, denoted by , which requires three main steps: (i) obtaining an estimator for the unpenalized problem based on a subset, say ; (ii) obtaining updated estimators for the unpenalized problem through one-step approximations using all subsets; and (iii) constructing an adaptive LASSO penalized estimator based on LSA. The procedure also brings a that well approximates .

Specifically, in step (i), we use subset to obtain a standard maximum PL estimator,

step (i) |

In step (ii), we obtain a DAC one-step approximation to ,

step (ii) |

where

(0.0) |

Let be our DAC approximation to . In practice, we find that it suffices to let . Finally, we apply the LSA to the PL and approximate using , where

step (iii) |

The optimization problem in step (iii) is equivalent to

(0.0) |

where is a vector and is a matrix. The linearization allows us to solve the penalized regression step using a pseudo likelihood based on a sample of size . The computation cost of this step compared to solving ( A Fast Divide-and-Conquer Sparse Cox Regression) reduces substantially when . In the Appendix, we show that . It then follows from the similar arguments given in Wang and Leng (2007) that if , , the estimated active set using achieves the variable selection consistency, i.e. and the oracle property holds, i.e. and have the same limiting distribution.

subsection2[Tuning and Standard Error Calculation]Tuning and Standard Error Calculation

The tuning parameter is chosen by minimizing the Bayesian information criteria (BIC) of the fitted model. Volinsky and Raftery (2000) showed that the exact Bayes factor can be better approximated for the Cox model if the number of uncensored cases, , is used to penalize the degrees of freedom in the Bayesian Information Criteria (BIC). Specifically, for any given tuning parameter with its corresponding estimate of , , the BIC suggested by Volinsky and Raftery (2000) is defined as

(0.0) |

where . With the LSA, we may further approximate by

(0.0) |

For the estimation of , we chose a such that is minimized. The oracle property is expected to hold in the setting where and is extraordinarily large. We may thus estimate the variance-covariance matrix for using . For , a confidence interval for can be calculated accordingly.

section1[Simulations]Simulations

subsection2[Simulation Settings]Simulation Settings We performed two sets of simulations to evaluate the performance of for the fitting of sparse Cox models, one with only time-independent covariates and the other with time-dependent covariates. For both scenarios, we let and . We consider the number of iterations and to examine the impact of on the proposed estimator.

subsubsection3[Time-independent covariates]Time-independent covariates We conducted extensive simulations to evaluate the performance of the proposed estimator relative to (a) the performance of the full sample-based adaptive LASSO estimator for the Cox model and (b) a majority voting-based DAC method for the Cox model, denoted by also with , penalized by a minimax concave penalty (MCP), which extends the majority voting-based DAC scheme for GLM proposed by Chen and Xie (2014). The reason of choosing as a comparison is that there is no other DAC method available for the Cox model and only Chen and Xie (2014) considered a similar majority voting-based DAC method for the penalized GLM with non-adaptive penalties. We set a priori that sets the estimate of a coefficient at zero, if at least 50% of the subset-specific estimates have a zero estimate for that coefficient. In addition, we compared the performance of the DAC estimator relative to the full sample maximum PL estimator .

For the penalized procedures, we selected the tuning parameter based on the BIC criterion discussed in section A Fast Divide-and-Conquer Sparse Cox Regression. The adaptive LASSO procedures were fit using the glmnet function in R with ; the MCP procedures were fit using the ncvsurv function in R.

For the covariates, we considered and . We generated from a multivariate normal distribution with mean and variance-covariance matrix , where denotes a vector with all elements being and we considered , , and to represent weak, moderate, and strong correlations among the covariates. For a given , we generated from a Weibull distribution with a shape parameter of 2 and a scale parameter of , where we considered three choices of to reflect different degrees of sparsity and signal strength:

For censoring, we generated from an exponential distribution with a rate parameter of , resulting in of censoring across different configurations.

subsubsection3[Time-dependent covariates]Time-dependent covariates We also conducted simulations for the settings where time-dependent covariates are present to evaluate the performance of . Since neither glmnet nor ncvsurv allows time-dependent survival data, we used as a benchmark to compare with. In addition, we compared the performance of relative to .

We considered consisting of time-independent covariates and time-dependent covariates. The simulation of the survival data with time-dependent covariates extended the simulation scheme of Austin (2012) from dichotomous time-dependent covariates to continuous time-dependent covariates. We considered four time intervals , , , and , where the time-dependent covariates are constant within each interval but can vary between intervals. We generated from a multivariate normal distribution with mean and variance-covariance matrix , where are the time-independent covariates and are the time-dependent for . We similarly considered to represent weak, moderate, and strong correlations.

We generated from a Weibull distribution with a shape parameter of 2 and a scale parameter of , where ,

We considered an administrative censoring with , leading to censoring under the three scenarios represented by weak, moderate, and strong correlations of the design matrix.

subsubsection3[Measures of performance]Measures of performance For any , we report (a) the average computation time for ; (b) the global mean squared error (GMSE), defined as ; (c) empirical probability of ; (d) the bias of each individual coefficient; and (e) mean squared error (MSE) of each individual coefficient. For and , we also report the empirical coverage level of the 95% normal confidence interval with standard error estimated as described in section A Fast Divide-and-Conquer Sparse Cox Regression. For any , we report (a) the average computation time for ; (b) the global mean squared error (GMSE), defined as .

The average computation time for each configuration is based on simulations using 50 simulated datasets performed on Intel^{®} Xeon^{®} E5-2620 v3 @2.40GHz. The statistical performance is evaluated based on 1000 simulated datasets for each configuration.

subsection2[Simulation Results]Simulation Results We first show in Table A Fast Divide-and-Conquer Sparse Cox Regression for the time-independent settings and in Table A Fast Divide-and-Conquer Sparse Cox Regression for the time-dependent settings the average computation time and GMSE of unpenalized estimators and . The results suggest that with two iterations () attains a GMSE comparable to the full sample-based estimator and reduced the computation time by more than 50%. The DAC estimator with two iterations () has a similar GMSE to . Across all settings, the results of are nearly identical with or and hence we summarize below the results for only for unless noted otherwise.

subsubsection3[Computation Time]Computation Time There are substantial differences in computation time across methods (Table A Fast Divide-and-Conquer Sparse Cox Regression- A Fast Divide-and-Conquer Sparse Cox Regression) for time-independent survival data. Across different settings, the average computation time of ranges from to seconds for and from to seconds for , with virtually all time spent on the computation of the unpenalized estimator . On the contrary, requires a substantially longer computation time with average time ranging from to seconds for and from to seconds for . This suggests that the computation time of is about of the full sample estimator when and about when . On the other hand, has a substantially longer average computation time than . This is because the MCP procedure, requiring more computational time than the adaptive LASSO, needs to be fitted times.

In the presence of time-dependent covariates, Table A Fast Divide-and-Conquer Sparse Cox Regression shows that has an average computation time of seconds for and ; has an average computation time of seconds. Virtually all computation time for and is spent on the computation of the unpenalized initial estimator , which has more observations and requires substantially more computation time compared to the setting with time-independent covariates given the same and .

subsubsection3[Statistical Performance]Statistical Performance The results for the simulation scenarios with only time-independent covariates are summarized in Table A Fast Divide-and-Conquer Sparse Cox Regression- A Fast Divide-and-Conquer Sparse Cox Regression. In general, is able to achieve a statistical performance comparable to , while generally has a worse performance, with respect to the GMSE and variable selection, bias, and MSE of individual coefficient. For example, as shown in Table A Fast Divide-and-Conquer Sparse Cox Regression, the GMSEs () for , and are respectively , , and when and ; , , when and . The relative performance of different procedures has similar patterns across different levels of correlation among the covariates. When the signals are relatively strong and sparse as for or , and have small biases and achieved perfect variable selection, while substantially excludes the signal when . For the more challenging case of where some of the signals are weak, the variable selection of and is also near perfect. Both penalized estimators for the weakest signal (0.035) exhibits a small amount of bias when and and an increased bias when . Such biases in the weak signals are expected for shrinkage estimators (Menelaos and others, 2016), especially in the presence of high correlation among covariates. However, it is important to note that and perform nearly identically, suggesting that our procedure incur negligible additional approximation errors. On the other hand, has difficulty in detecting the 0.05 and 0.035 signals and tends to produce substantially higher MSE than . The empirical coverage levels for the confidence intervals are close to the nominal level across all settings except for the very challenging set very weak signals when the correlation is . This again is due to the bias inherent in shrinkage estimators..

The results for the time-dependent survival are summarized in Table A Fast Divide-and-Conquer Sparse Cox Regression. We find that also generally has a good performance in estimating for both time-independent and time-dependent covariates. The variable selection consistency holds perfectly for all parameters of interest. The coverage of the confidence intervals also has similar patterns as the case with time-independent covariates.

section1[Application of the DAC procedure to Medicare Data]Application of the DAC procedure to Medicare Data

We applied the proposed algorithm to develop risk prediction models for heart failure-specific readmission or death within 30 days of discharge among Medicare patients who were admitted due to heart failure. The Medicare inpatient claims were assembled for all Medicare fee-for-service beneficiaries during to identify the eligible study population. The index date was defined as the discharge date of the first heart failure admission of each patient. We restricted the study population to patients who were discharged alive from the first heart failure admission. The outcome of interest is time to heart failure-specific readmission or death after the first heart failure admission. Because readmission rates within 30 days have been used to assess the quality of care at hospitals by the Centers for Medicare and Medicaid Services (CMS) (CMS, 2016), we censored the time to readmission at 30 days. For a patient who were readmitted or dead on the same day as discharge (whose claim did not indicate discharge dead), the time-to-event was set at 0.5 days. Due to the large number of ICD-9 codes, we classified each discharge ICD-9 code into disease phenotypes indexed by phenotype codes according to Denny and others (2013). A heart failure admission or readmission was identified, if the claim for that admission or readmission had a heart failure phenotype code at discharge.

We consider two sets of covariates: (I) time-independent covariates including baseline individual-level covariates collected at time of discharge from the index heart failure hospitalization, baseline area-level covariates at the residential ZIP code of each patient, and indicators for time trend including include dummy variables for each year and dummy variables for each months, and (II) time-dependent predictors that vary day-by-day. Baseline individual-level covariates include age, sex, race (white, black, others), calendar year and month of the discharge, Charlson co-morbidity index (CCI) (Quan and others, 2005) which describes the degree of illness of a patient, and indicators for non-rare co-morbidities (defined as prevalence among the study population). Baseline area-level covariates include socioeconomic status variables [percent black residents (ranging from 0 to 1), percent Hispanic residents (ranging from 0 to 1), median household income (per ten thousand increase), median home value (per ten thousand increase), percent below poverty (ranging from 0 to 1), percent below high school (ranging from 0 to 1), percent owned houses (ranging from 0 to 1)], population density (1000 per squared kilometer), and health status variables [percent taking hemoglobin A1C test (ranging 0-1), average BMI, percent ambulance use (ranging from 0 to 1), percent having low-density lipoprotein test (ranging from 0 to 1), and smoke rate (ranging from 0 to 1)]. The time-dependent covariates include daily fine particulate matter (PM) predicted using a neural network algorithm (Di and others, 2016), daily temperature with its quadratic form, and daily dew point temperature with its quadratic form. There were 574 time-independent covariates and 5 time-dependent covariates.

There were eligible patients with a total of heart failure readmissions or deaths, among which were readmissions and were deaths. After expanding the dataset by accounting for time-dependent variables which vary day-by-day, the time-dependent dataset contains rows of records.

We fit cause-specific Cox models for readmission due to heart failure or deaths as a composite outcome, considering two separate models: (i) a model containing only time-independent covariates and (ii) a model incorporating time-dependent covariates. In both cases, the datasets are too large for glmnet package to analyze as a whole, demonstrating the need for .

subsection2[Time-independent Covariates Only]Time-independent Covariates Only We applied with and paralleled on 25 Authentic AMD Little Endian @2.30GHz CPUs. Computing with took 1.1 hours, including the time of reading datasets from hard drives during each iteration of the update of the one-step estimator. Figure 1 shows the hazard ratio of each covariate based on with predicting heart failure-specific readmission and death within 30 days.

Multiple co-morbidities were associated with an increased risk of 30-day readmission or death with the leading factors including renal failures, cancers, malnutrition, subdural or intracerebral hemorrhage, myocardial infarction, endocarditis, respiratory failure, and cardiac arrest. CCI was also associated with an increased hazard of the outcome. These findings are generally consistent with those reported in the literature. For example, Philbin and DiSalvo (1999) reported that ischemic heart disease, diabetes, renal diseases, and idiopathic cardiomyopathy were associated with an increased risk of heart failure-specific readmission within a year. Leading factors negatively associated with readmissions included virus infections, asthma, and chronic kidney disease in earlier stages. These negative association findings are reflective of both clinical practice patterns and the biological effects, as most of the negative predictors are generally less severe than the positive predictors.

Some socioeconomic status predictors were relatively less important in predicting the outcome after accounting for the phenotypes, where percent black, median household income, and percent below poverty were dropped and dual eligibility, median home value, percent below high school had a small hazard ratio. By comparison, Philbin and others (2001) reported a decrease in readmission as neighborhood income increased. Foraker and others (2011) reported that given co-morbidity measured by CCI, the readmission or death hazard was higher for low socioeconomic status patients. The present paper considered more detailed phenotypes in addition to CCI suggested a relatively smaller impact of socioeconomic status. The difference in results is possibly because co-morbidity may be on the causal pathway between socioeconomic status and readmission or death. Adjusting for a detailed set of co-morbidities partially blocks the effect of socioeconomic status. Percent Hispanic residents was negatively associated with readmission or death. Percent occupied houses increased the risk of readmission or death, which was consistent with the strong positive prediction by population density. Most ecological health variable showed a small hazard ratio.

Black and other race groups had a lower hazard than white. Females had a lower hazard than males, which was consistent with Roger and others (2004) that females had a higher survival rate than males after heart failure. Age was associated with an increased hazard of readmission or death, as expected.

The coefficient by month suggested a higher risk of readmission or death in cold seasons than warm seasons, with a larger negative hazard ratio for summer indicators. The short-term readmission or death rate was decreasing over time, which was suggested by the negative hazard ratio of later years. The later calendar year being negatively associated with readmission risk may be an indication of improved follow-up care for patients discharged from heart failure. Consistently, (Roger and others, 2004) also suggests an improved heart failure survival rate over time.

subsection2[Incorporating Time-Dependent Covariates]Incorporating Time-Dependent Covariates The analysis has two goals. First, the covariates serve as the risk predictors of the hazard of heart failure-specific readmission. Second, all covariates other than PM serve as the potential confounders of the association between PM and readmission, particularly time trend and area-level covariates. The procedure is a variable selection technique to drop non-informative confounders given the high dimensionality of confounders. This goal aligns with Belloni and others (2014), which constructs separate penalized regressions for the propensity score model and outcome regression model to identify confounders. We herein focused on building a penalized regression for the outcome regression model.

We applied algorithm with to this time-dependent survival dataset. The procedure was paralleled on 10 Authentic AMD Little Endian @2.30GHz CPUs. The estimation of with took 36.5 hours, including the time of loading the datasets into memory. The result suggests each 10 increase in daily PM was associated with increase of risk (95% confidence interval ) adjusting for individual-level, area-level covariates, and temperature. Because there is rare evidence on whether air pollution is associated with heart failure-specific readmission or death among heart failure patients and it is rare to estimate the health effect of daily air pollution using a time-dependent Cox model, this model provides a novel approach to address a new research question. While evidence is rare on the association between daily PM and heart failure-specific readmission, some studies used case-crossover design to estimate the effect of short-term PM on the incidence of heart failure admissions. Pope and others (2008) found that a 10 increase in 14-day moving average PM was associated with a 13.1% () increase in the incidence of heart failure admissions among elderly patients; Zanobetti and others (2009) reported that each 10 increase in 2-day averaged PM was associated with a 1.85% () increase in the incidence of congestive heart failure admission. There is also a large body of literature suggesting that short-term exposure to PM is associated with an increased risk of death. For example, (Di and others, 2017) shows among the Medicare population during that each 10 increase in PM was associated with an 1.05% (0.95%, 1.15%) increase in mortality risk. In addition, Figure 2 shows the covariate-specific estimates of the hazard ratio for all the covariates, with the estimates consistent with the analysis of time-independent dataset.

section1[Discussions]Discussions The proposed procedure for fitting adaptive LASSO penalized Cox model reduces the computation cost, while it maintains the precision of estimation and accuracy in variable selection with an extraordinarily large and a numerically large . The use of makes it feasible to obtain the -consistent estimator required by penalized step (e.g. when there is a constraint in RAM) and shortens the computation time of the initial estimator by . The improvement in the computation time was substantial in the regularized regression step. The LSA converted the fitting of regularized regression from using a dataset of size to a dataset of size .

The majority voting-based method with MCP (Chen and Xie, 2014) had a substantially longer computation time than . The difference primarily comes from (a) the fact that the Cox model with MCP is fitted times and (b) the computation efficiency between glmnet algorithm which is more efficient than the MCP algorithm in ncvsurv (Breheny and Huang, 2011).

The difference in variable selection between and (Chen and Xie, 2014) is primarily due to the majority voting. The simulations in Chen and Xie (2014) have shown that an increase in the percentage for the majority voting decreased the sensitivity and increased the specificity of variable selection. Similarly in the simulations of the present study with a of the majority vote, Chen and Xie’s procedure showed a high specificity but the sensitivity was low for weaker signals as demonstrated in the simulation studies.

For non-weak signals, the oracle property appears to hold well as evidenced by the simulation results for and shown in Tables A Fast Divide-and-Conquer Sparse Cox Regression- A Fast Divide-and-Conquer Sparse Cox Regression. For weak signals such as 0.035 in , the oracle property does not appear to hold even with and the bias in the shrinkage estimators results in confidence intervals with low coverage. This is consistent with the the impossibility result shown in Potscher and Schneider (2009), which suggests difficulty in deriving precise interval estimators when adaptive LASSO penalty is employed.

Acknowledgments
The study was partially supported by an USEPA grant
RD-83587201, an NIH grant U54 HG007963, and NIEHS grants R01 ES024332 and P30 ES000002. Its contents are solely the responsibility of the grantee and do not necessarily represent the official views of the USEPA. Further, USEPA does not endorse the purchase of any commercial products or services mentioned in the publication. The analysis of the Medicare data was performed on the secured clusters of Research Computing, Faculty of Arts and Sciences at Harvard University.

Conflict of Interest: None declared.

section1[Appendix]Appendix Throughout we assume all regularity conditions required in Zhang and Lu (2007), is fixed and belongs to a compact support. We next establish the asymptotic equivalence of and in that . To this end, we note that , , and , where . It follows that . From a taylor series expansion, it’s straightforward to see that for ,

On the other hand, . Therefore,

Furthermore, from the convergence rate of and the fact that , we have , where

It follows that

and therefore

Similarly, we may show that

This implies that