Identifying Factors Predicting Kidney Graft Survival in Chile Using Elastic-Net-Regularized Cox’s Regression

Background and Objectives: We developed a predictive statistical model to identify donor–recipient characteristics related to kidney graft survival in the Chilean population. Given the large number of potential predictors relative to the sample size, we implemented an automated variable selection mechanism that could be revised in future studies as more national data is collected. Materials and Methods: A retrospective multicenter study was conducted to analyze data from 822 adult kidney transplant recipients from adult donors between 1998 and 2018. To the best of our knowledge, this is the largest kidney transplant database to date in Chile. A procedure based on a cross-validated regularized Cox regression using the Elastic Net penalty was applied to objectively identify predictors of death-censored graft failure. Hazard ratios were estimated by adjusting a multivariate Cox regression with the selected predictors. Results: Seven variables were associated with the risk of death-censored graft failure; four from the donor: age (HR = 1.02, 95% CI: 1.00–1.03), male sex (HR = 0.64, 95% CI: 0.46–0.90), history of hypertension (HR = 1.49, 95% CI: 0.98–2.28), and history of diabetes (HR = 2.04, 95% CI: 0.97–4.29); two from the recipient: years on dialysis log-transformation (HR = 1.29, 95% CI: 0.99–1.67) and history of previous solid organ transplantation (HR = 2.02, 95% CI: 1.18–3.47); and one from the transplant: number of HLA mismatches (HR = 1.13, 95% CI: 0.99–1.28). Only the latter is considered for patient prioritization in deceased kidney allocation in Chile. Conclusions: A risk model for kidney graft failure was developed and trained for the Chilean population, providing objective criteria which can be used to improve efficiency in deceased kidney allocation.


Introduction
The shortage of deceased donor kidneys available for transplantation has increased the relevance of allocation mechanisms. To better understand the factors predicting transplant outcomes, several survival analysis methodologies have been applied in the literature [1][2][3][4].
The Cox proportional hazards model is the most widely used for assessing the risk of kidney graft failure as it provides comprehensible results that aid clinicians and decision makers. This methodology has been incorporated in the estimated post-transplant survival (EPTS) and kidney donor profile index (KDPI) models used for matching patients with good prognosis after transplant with kidneys with high predicted outcomes in the United States [5,6]. The EPTS considers four recipient factors: age, history of diabetes, history of previous solid organ transplant, and time on dialysis [5]. On the other hand, the KDPI considers ten deceased donor characteristics, including age, height, weight, race, creatinine level, and history of hypertension [6]. In Chile, the deceased kidney allocation system mainly considers HLA mismatch to prioritize patients with good transplant prognosis. Authors have suggested including factors such as donor-recipient age difference [7,8] for longevity matching and for maximizing the utilization of deceased donor kidneys. In this context, it is relevant to provide decision makers with objective information about the effect these factors have shown in historical data from the Chilean population. This may lead to better allocation mechanisms, based on statistical evidence that guides prioritization proportional to the benefit each factor has shown on transplant outcomes. Furthermore, it is of interest to identify which factors are associated with kidney graft survival and to develop a risk model for kidney graft failure capable of revealing improvement opportunities in deceased organ allocation.
Therefore, we have collected data from multiple transplant centers to construct the largest database of kidney-transplanted patients in Chile, comprising 1459 transplants in total. Using this evidence, we have undertaken an in-depth analysis to objectively identify factors influencing transplant outcomes in Chile. However, this sample size is relatively small compared to the number of factors that can be used to predict graft survival [1,6,9]. For this reason, we combined survival analysis with an objective variable selection criterion based on a cross-validated regularized Cox model using the Elastic Net penalty, which provides transparency in the model specification process and could be applied in future studies as more data are collected from transplant centers.

Materials and Methods
A major data collection effort was conducted in collaboration with five transplant centers in Chile: Hospital Barros Luco Trudeau, Hospital del Salvador, Hospital Las Higueras, Hospital Sótero del Río, and Clínica Santa María, which jointly comprise approximately 14% of the total number of transplants conducted in Chile. To the best of our knowledge, this is the largest kidney transplant dataset that has been collected to date in Chile. This process required the digitalization and standardization of thousands of physical files stored in each medical center. The resulting database was implemented in REDCap and registered 1459 kidney transplant cases. Specifically, it contains pre-transplant information of donorrecipient characteristics and follow-up diagnosis post-transplant, including information on graft failure, with a total of 25 pre-transplant characteristics, which we sought to evaluate as potential predictors of graft survival. Anonymity was preserved for clinical records used in the analysis.
To select the final sample for the study, we applied the criteria summarized in Figure 1. From the 1459 original transplants, 637 entries were discarded sequentially for the following reasons: 244 had transplant dates outside the study period 1998-2018, 12 due to inconsistent recorded survival time or dialysis time, 333 with missing values on variables for which the imputation would be complicated-survival time, HLA mismatch, end-stage renal disease cause, and recipient comorbidities-(later we show how we were able to use these observations to conduct an out-of-sample test), and 48 transplants with pediatric donor or recipient (age < 18). As a result, the final sample comprised 822 transplants, of which: (i) 140 graft failures were observed; (ii) 58 recipients died before graft failure was reported; and (iii) 624 had a functioning graft in the last follow-up. Table 1 reports summary statistics and the percentage of missing values of the main characteristics collected in the study. Five variables presented missing information that needs to be addressed. To avoid a relevant reduction in sample size and the deletion of relevant information provided by other variables, missing values were replaced by reasonable estimates using MissForest [10]. This nonparametric machine learning algorithm based on random forest has shown good performance imputing clinical data [11]. reasonable estimates using MissForest [10]. This nonparametric machine learning algorithm based on random forest has shown good performance imputing clinical data [11].    Two definitions of graft failure are typically used in survival analysis after a kidney transplant [12][13][14]. Overall, graft failure considers the time of transplantation to the date of irreversible graft failure (return to dialysis or a new transplant) or the time of death, i.e., death is considered graft failure. This time may be censored by the date of the last follow-up with a functioning graft. Graft survival censored for death with functioning graft (death-censored graft failure) treats patient death as a censoring event equivalent to the last follow-up, and hence considers a functioning graft at the time of death. The present study used death-censored survival analysis, which more adequately accounts for varying death rates from other causes in the patient population [14].
Preliminary evidence was evaluated through different analyses to better understand the influence between potential predictors and transplant outcomes. First, an exploratory data analysis and a statistical description were computed. Distribution plots for survival and censorship times were examined, and variables were described by their mean and standard deviation or by their frequency and proportion in categories. Second, survival curves were estimated using the nonparametric Kaplan-Meier method for different subsets of patients [15], and the null hypothesis of equality between curves was contrasted using the Mantel-Cox test, considering a 10% significance level [16]. This approximates the survival function, allowing us to assess transplant outcomes in different groups of patients. Finally, penalized splines with 5 degrees of freedom were adjusted in univariate Cox regression for non-binary variables to determine functional transformations that may better fit Cox's log-linearity assumption [17]. We considered a 10% significance level for detecting statistical evidence of nonlinear effects [18], and we examined complementary graphical evidence for determining adequate functional transformations.
Despite our significant data collecting effort, the sample size was still too limited to evaluate the predictive power of the 15+ factors that needed to be tested, including alternative functional forms. Consequently, we applied an automatized procedure to objectively identify predictors for death-censored graft failure based on a cross-validated regularized Cox model using the Elastic Net penalty [19]. Regularization is a technique for creating generalized models which consists of applying a penalty over the estimated coefficients in order to reduce overfitting and assist variable selection. The Elastic Net penalty is a convex combination of the Lasso and Ridge penalties. Therefore, it combines the strength of the two approaches. On one hand, the Lasso penalty chooses only a few nonzero coefficients, performing variable selection but causing undesirable behavior in the presence of correlated predictors [20]. On the other hand, the Ridge penalty reduces the value of the coefficients towards zero proportionally but sets none to exactly zero; therefore, it better handles correlated predictors but fails to exclude irrelevant variables from the model. The value of α ∈ [0, 1] determines the trade-off between Lasso and Ridge penalties in the Elastic Net. With α = 0.95, the Elastic Net behaves similarly to the Lasso, only removing degenerate behavior due to extreme correlations [21]. We considered this approach (α = 0.95) to prioritize variable selection while maintaining stability in the presence of highly correlated predictors. The optimal penalty in the regularized model was determined by applying 10-fold cross-validation using the Harrell's Concordance Index ("C-statistic" or "C-index") as the ensemble metric [22]. To obtain a simplified model, we considered the highest penalty with a C-index no more than one standard deviation from the optimum. This procedure was repeated 1000 times to reduce group selection noise in the 10-fold cross-validation. The selected variables (i.e., those associated with nonzero coefficients) in more than half of the scenarios were defined as the identified predictors by the procedure (see Algorithm S1). These were used as covariables in a multivariate Cox model to estimate hazard ratios with their respective confidence intervals. Finally, to communicate model predictions beyond hazard ratios, survival curves were estimated for different transplant cases using the Breslow estimator to approximate the cumulative baseline hazard function [23].
Predictive accuracy was assessed by computing the C-statistic, which indicates the proportion of comparable pairs ranked correctly by the model. In addition to in-sample predictions, predictive accuracy was also assessed out-of-sample using 76 kidney transplants not included in the training data due to missing information. These transplants had missing data in recipient comorbidities assessed in the variable selection process. Since some of these were not selected in the final model, this sub-sample was used for evaluating the out-of-sample prediction of the final model specification.
The proportionality assumption in the model was tested using the scaled Schoenfeld residuals [24] and considering the chi-squared distributed statistic proposed by Grambsch and Therneau (1994) [25]. We contrasted the null hypothesis of proportionality considering a 5% significance level, and we examined the graphical evidence of the evolution in time of the scaled residuals for each predictor. On the other hand, the log-linearity assumption was examined in non-binary predictors adjusting penalized splines with five degrees of freedom independently on multivariate Cox regressions [18]. We considered a 5% significance level for null hypothesis rejection and used graphical evidence as a complement.
The statistical analyses were conducted using R (v4.1.2; R Core Team 2021; Vienna, Austria), while data processing was carried out using Python (v3.8.8; Python Software Foundation 2021; Fredericksburg, VA, USA). Detailed information about the libraries and functions used can be found in the Supplementary Materials.

Results
A visual description of the study cohort is shown in Figure S1. MissForest was used to impute missing values in donor history of diabetes (9.8% missing), donor creatinine > 1.5 mg/dL (9.6%), recipient weight (6.8%), recipient years on dialysis (2.5%), and cold ischemia time (0.9%). The statistical comparison between imputed and pre-existing values showed a reasonable estimation by the machine learning algorithm (see Table S1). Table 1 shows a statistical summary of donor-recipient characteristics studied as potential predictors for death-censored graft failure. Most kidneys come from deceased donors (80%), male sex (58%), and without history of hypertension (81%). On the other hand, most recipients present history of hypertension (83%), a maximum panel-reactive antibodies (PRA) less than or equal to 10% (64%), and at least two HLA mismatches (88%). In addition, the average recipient is 45 years old, weighs 67 kg, and has spent 3 years on dialysis therapy, while the average donor is approximately 44 years old.
The preliminary log-linearity analysis in non-binary variables using penalized splines is shown in Figure S5. A nonlinear effect is observed for recipient age; the risk of graft failure increases at a higher pace beyond 50 years old. The nonlinear component of the penalized spline is significant (p = 0.044) at a 10% significance level, while the linear component is not (p = 0.732). Considering this evidence, a functional transformation was included to consider an effect for recipients beyond 50 years old. On the other hand, the rest of examined variables showed a reasonable log-linearity adjustment. Nonetheless, three extra functional forms were incorporated considering graphical evidence and medical criteria: a binary variable that becomes active when maximum recipient PRA is greater than 50%, and log-transformations for recipient weight and recipient years on dialysis. These four functional transformations included for variable selection are shown in Table S5.
The variable selection mechanism is executed, including 25 potential predictors of death-censored graft failure (see Table 2, which includes nonlinear transformations of the variables described in Table 1). Figure 2 shows the number of selected predictors for different penalty values in one iteration of the algorithm. The vertical right line indicates the selected penalty, and the number of selected variables (i.e., with nonzero coefficients in the regularized model) is shown in the top row. The number of times each variable is selected in 1000 iterations is shown in Table 2. Seven variables were selected as predictors by the algorithm: donor age, donor male sex, donor history of hypertension, donor history of diabetes, previous solid organ transplantation in the recipient, the recipient's years on dialysis log-transformation, and the number of HLA mismatches. The Cox multivariate model adjusted to these factors is presented in Table 3 and presents a C-index of 0.659 in the training data (n = 822) and 0.733 in testing out-of-sample data (n = 76).  The proportionality test shows insufficient evidence to reject the null hypothesis and, therefore, supports the chosen model specification (see Table S6). Specifically, the p-value for the overall test is 0.4, and when evaluating each variable independently, the lowest p-value is 0.092 (corresponding to the number of HLA mismatches). In addition, graphical evidence shows a reasonable adjustment of the proportionality assumption (see Figure  S7). Correspondingly, the test for log-linearity in non-binary predictors shows a reasonable adjustment (see Figure S6). The nonlinear components were insignificant for donor age, mismatch HLA, and recipient years on dialysis log-transformation.
ferent penalty values in one iteration of the algorithm. The vertical right line indicates the selected penalty, and the number of selected variables (i.e., with nonzero coefficients in the regularized model) is shown in the top row. The number of times each variable is selected in 1000 iterations is shown in Table 2. Seven variables were selected as predictors by the algorithm: donor age, donor male sex, donor history of hypertension, donor history of diabetes, previous solid organ transplantation in the recipient, the recipient's years on dialysis log-transformation, and the number of HLA mismatches. The Cox multivariate model adjusted to these factors is presented in Table 3 and presents a C-index of 0.659 in the training data (n = 822) and 0.733 in testing out-of-sample data (n = 76). The proportionality test shows insufficient evidence to reject the null hypothesis and, therefore, supports the chosen model specification (see Table S6). Specifically, the p-value To further explain estimated effects, survival curves were derived from model estimates using Breslow's estimator to approximate the cumulative baseline hazard. Figure 3 shows the estimated survival curves for different values in each predictor while keeping the rest constant in the reference values. This graphical evidence provides comprehensible information of model specification under proportionality and log-linearity assumptions.  Table  3)

Discussion
Regardless of the limitations given by the sample size of the Chilean transplanted population, the results of this study are reasonable considering previous findings in the literature. Donor age and history of hypertension are significant predictors for kidney graft survival identified in populations from the United States, the United Kingdom, and Thailand [6,9,26]. In addition, recipient history of previous solid organ transplant, recipi-  Table 3

Discussion
Regardless of the limitations given by the sample size of the Chilean transplanted population, the results of this study are reasonable considering previous findings in the literature. Donor age and history of hypertension are significant predictors for kidney graft survival identified in populations from the United States, the United Kingdom, and Thailand [6,9,26]. In addition, recipient history of previous solid organ transplant, recipient years on dialysis log-transformation and mismatched HLA have been presented as relevant predictors in several cohort studies [2][3][4] and are included in the EPTS score used for assessing the patient risk of graft failure in the United States [5]. However, regarding the estimated effect for donor sex, which suggests a higher survival in male donor kidneys, this effect could be confounded through other relevant predictors not included in our study due to lack of data-such as donor weight, height, and donor body mass index-which have shown significant associations with transplant outcome in the literature [6,27]. These factors were not included in this study because of missing values in more than 50% of the entries. To analyze whether these omitted variables could be generating a bias in the gender coefficient, we compared means and proportions between donor genders using the Wilcoxon's test (see Table S7). Donor weight, height, and cause of death by CVA are statistically different across genders; in particular, male donors have a higher value for weight and height and a higher proportion of deaths caused by CVA. On the other hand, the donor body mass index did not show significant differences between genders. Overall, this analysis suggests that the effect of donor gender could be attributed to omitted variables characterizing the donor; hence, the result should be interpreted with caution. More information about the donor is needed to assess whether donor gender is associated with differences in death-censored graft survival.
Of the six predictors identified, only donor-recipient HLA mismatch is currently considered in the score used to prioritize patients with good transplant prognosis in Chile. In addition, the estimated effects suggest that the current score distinctions between patients may need adjustment to be proportional to the benefits of transplant survival, which may be an opportunity to improve deceased kidney allocation. Furthermore, the five prognostic factors not considered for organ allocation open room for discussion about the mechanisms, based on national evidence, that increase patients' overall years of therapy.
We compared the predictive power of the model against the C-index that could be attained using the Kidney Donor Profile Index (KDPI) [6]. The KDPI is a numerical measure that combines 10 donor factors to predict graft survival, which is used to match with adequate recipients for kidney allocation. The KDPI is derived from the Kidney Donor Risk Index (KDRI), which is computed using ten factors' coefficients estimated through a multivariable Cox proportional hazards regression using a sample of deceased donors in the United States. These factors include donor's age, height, weight, ethnicity, history of hypertension, history of diabetes, cause of death, serum creatinine, hepatitis C virus (HCV) status, and donation after circulatory death (DCD) status.
We applied the KDRI score with the patient sample analyzed in this study, calculating the C-index to assess its predictive power. Since height, weight, HCV status, and DCD exhibited significant missing values in the data, these factors were set at reference levels to compute the KDRI. The predictive power of this adjusted KDRI score yields a concordance index (C-index) of 0.606, which is lower than the C-index calculated in this study (C-index = 0.659). This comparison provided further support that the methodology developed in this work is useful to improve predictions of graft survival, which can potentially be used to construct a "local" KDPI score in Chile's kidney allocation system.

Conclusions
This study analyzed the largest sample of kidney transplants collected up to date in Chile, comprising 1459 transplants from five transplant centers accounting for 14% of the transplants performed from 1998 to 2018. An automated procedure based on the Elastic-Net-regularized Cox's regression was applied to objectively select predictors for kidney death-censored graft failure in the Chilean population, and a risk model that provides comprehensible results was developed to aid clinicians and decision markers. It is relevant that the data collection process continues to increase the accuracy of the selected variables and the precision of the estimated effects through the proposed methodology. In this regard, the procedure can be used as more information is collected if proportionality and log-linearity conditions for the Cox model are adequately validated.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/medicina58101348/s1, Algorithm S1: Training algorithm for variable selection and model fitting, Table S1: Statistical comparison between imputed and preexisting values, Table S2: Categories used for numeric variables in Kaplan-Meier analysis, Table S3: Descriptive summary of recipient characteristics, Table S4: Descriptive summary of donor and transplant characteristics, Table S5: Transformations included as study variables, Table S6: Proportionality assumption tests in each predictor and in the overall model, Table S7: Comparison of means and proportions between donor genders, Figure S1: Exploratory data analysis, Figure S2: Graphic representation of Spearman's correlation matrix, Figure S3: Kaplan-Meier analysis in donor characteristics, Figure S4: Kaplan-Meier analysis in recipient characteristics, Figure S5: Preliminary log-linearity assumption analysis, Figure S6: Log-linearity assumption analysis via penalized smoothing splines, Figure S7: Proportionality assumption analysis via scaled Schoenfeld residuals.