Incorporating External Risk Information with the Cox Model under Population Heterogeneity: Applications to Trans-Ancestry Polygenic Hazard Scores

Survival risk discrimination is important for identifying individuals with a high risk of disease progression. For example, polygenic hazard scores (PHS), calculated as the sum of single-nucleotide polymorphisms (SNPs) weighted by the risk allele effect sizes derived from Cox proportional hazard models, are effective tools for personalized disease risk assessment and are widely applied in cancer research (Seibert et al., 2018; Zhang et al., 2020). However, the development of PHS has often included only individuals of European ancestry (EUR). A major challenge in developing PHS models for a minority population is that non-EUR individuals are usually less represented in genetic databases (Landry et al., 2018), limiting their utility and exacerbating existing health disparities (Martin et al., 2017).

Our research is motivated by the study of prostate cancer, which is the second most prevalent cancer in males across the world (Bray et al., 2018). Recently, a PHS calculated with 46 SNPs (denoted as PHS46) has been developed to predict the age of onset of prostate cancer in EUR individuals (Huynh-Le et al., 2021). However, it is well-known that African ancestry (AFR) individuals have higher prostate cancer incidences and mortality rates than EUR individuals (DeSantis et al., 2016), which may be partially explained by disparities in gene mutations associated with prostate cancer (Huang et al., 2017). Due to these clinical and genetic heterogeneity, directly applying EUR-based PHS has limited utility for non-EUR populations (Vilhjálmsson et al., 2015). On the other hand, constructing PHS for non-EUR populations suffers from limited sample sizes, high dimensionality, and low signal-to-noise ratios. Thus, it would be desirable to develop a method that integrates risk information from published prediction models upon EUR individuals (e.g. PHS46) with individual-level data collected in non-EUR cohorts (e.g. genotype data of AFR males), automatically measure the heterogeneity across different populations, and efficiently balance the contribution of each information source.

Ideally, the access to individual-level data will enable the pooled analysis and optimize the use of each information source. However, due to restrictions on data privacy, summary-level information is commonly used for data sharing in most applications. In the context of survival analysis, several methods for integrating summary-level external information with an internal study have been proposed. For example, Huang et al. (2016) proposed a double empirical likelihood approach based on a homogeneous assumption; i.e. both the covariate distribution and the association between covariates and outcomes are the same in the internal and external cohorts. However, this strong assumption is rarely met. Ignoring heterogeneity across different studies yields substantially biased parameter estimates and a loss of efficiency (Estes et al., 2017).

To relax the homogeneity assumption, Chen et al. (2021) proposed an adaptive double empirical likelihood estimator by penalizing potential differences on the subgroup survival probabilities. Sheng et al. (2021) constructed a heterogeneity test based on the penalized empirical likelihood and proposed a semiparametric density ratio model. Although these methods are useful, they typically require external subgroup survival probabilities, which are not available through EUR-based PHS models. The published PHS contains only covariate coefficients, but with no baseline hazard information. Such risk scores cannot be given any direct probability interpretation. Furthermore, it can be computationally expensive for existing methods to handle large-scale data, especially when the number of predictors largely exceeds the sample size (e.g. non-EUR genotype data), leading to a high demand for new integration methods that are computationally efficient and enable full utilization of high-dimensional information.

Kullback-Leibler (KL) information (Kullback and Leibler, 1951) has been widely applied in the development of data integration methods for binary and ordinal outcomes under population heterogeneity (Liu and Shum, 2003; Schapire et al., 2005; Jiang et al., 2016; Hector and Martin, 2022). However, these methods apply only to generalized linear models (GLM). For the integration with survival outcomes, Wang et al. (2021) extended the GLM-based KL to discrete failure time models and developed an integration procedure to incorporate external survival probabilities across multiple time points with internal time-to-event data. However, it is not clear how to incorporate external survival risk information without probability interpretation. To fill this gap, we propose a KL-based integrated Cox model (CoxKL), which incorporates external risk scores with internal individual-level time-to-event data to improve the estimation and prediction performance in the internal cohort. A key ingredient of the proposed method is the partial-likelihood-based KL information, which measures the heterogeneity between the external and the internal populations. We first develop the CoxKL method under low-dimensional settings, and then extend it to CoxKL-LASSO with a $L_{1}$ penalization, which can be further applied to high-dimensional trans-ancestry PHS developments.

The proposed CoxKL procedure attains the following advantages: (1) it only requires summary-level risk scores to be obtained from external studies, but with no use of any individual-level external data or survival probabilities, and therefore achieves greater flexibility in data sharing and integrative analysis; (2) it accounts for population heterogeneity, automatically detects the quality of different information sources, and provides an optimal integrated model based on the aggregated information; and (3) its objective function retains a similar form of the log-partial likelihood of the classical Cox model; thus, it is computationally efficient for high-dimensional problems, which is crucial for the motivating trans-ancestry PHS analysis; and (4) it provides a unified integration framework to incorporate risk information from various types of external prediction models or multiple external models.

The rest of the paper is organized as follows: In Section 2, we introduce the notations and model setup. In Section 3, we propose the CoxKL method. Asymptotic properties are discussed in Section 4. Section 5 conducts simulations to evaluate the proposed method. In Section 6, we applied CoxKL method to develop a trans-ancestry PHS model for prostate cancer in AFR males. Section 7 concludes with a brief discussion.

In Section 2.1, we introduce the notations of the Cox proportional hazards model for the target internal cohort, which collects individual-level data with time-to-event outcomes and a set of covariates from a limited-size sample (denoted as internal data). In Section 2.2, we describe various types of external risk scores that can be integrated by the CoxKL method.

We propose using KL information to evaluate the disparity between the internal data and external risk scores, and consequently integrate the optimally weighted external information with the internal data to obtain the best integrated model for the target internal population. We present CoxKL under low-dimensional settings in Section 3.1, and then present CoxKL-LASSO, which extends CoxKL to high-dimensional trans-ancestry PHS problems, in Section 3.3.

In this section, we study the asymptotic distribution of the CoxKL estimator under the low-dimensional setting (i.e. $p$ is fixed and $p\ll n$) and show that the proposed method improves the asymptotic mean squared error (aMSE) given a certain range of $\eta$.

We first provide some prerequisite notations. For a given $\eta$, let $\hat{\bm{\beta}}(\eta)$ be the solution obtained by maximizing the penalized objective function in (7), and let $\bm{\beta}^{*}(\eta)$ be the maximum of the expectation of the penalized objective function $\mathbb{E}\{\ell_{\eta}(\bm{\beta}))\}$. Denote $\hat{\bm{\beta}}$ as the maximum log-partial likelihood estimator of the internal model in (2). The true value of the regression parameter for the Cox model using the internal cohort only is denoted by $\bm{\beta}_{0}$. Define $\bm{\Delta}=\bm{\beta}_{0}-\widetilde{\bm{\beta}}$, where $\widetilde{{\bm{\beta}}}$ are the estimated coefficients from the external model. For the ease of notations, we assume that the external risk scores can be constructed as $\bm{Z}_{i}^{\top}\widetilde{\boldsymbol{\beta}}$. Denote

where $\bm{b}^{\otimes 2}=\bm{b}\bm{b}^{\top}$. Then the score function of $\ell_{\eta}(\bm{\beta})$ in (7) is given by

The corresponding information matrix is

To derive the theoretical properties of the CoxKL method under the low-dimensional setting, we impose the following regularity conditions for the internal cohort:

• (C1).

  $(X_{i},\delta_{i},\bm{Z}_{i})$, $i=1,\ldots,n$, are independent and identically distributed random vectors.

• (C2).

  $P(X_{i}\geq\tau)=0$, where $\tau$ is a pre-specified time point.

• (C3).

  $Z_{ij}<\kappa$ for $i=1,\ldots,n$, and $j=1,\dots,p$, where $\kappa$ is a constant and $Z_{ij}$ is the $j$-th component of $\bm{Z}_{i}$.

• (C4).

  $\int_{0}^{\tau}\lambda_{0}(t)<\infty$.

• (C5).

  Continuity of $s^{(d)}(t;\bm{\beta})$, $d=0,1,2$, for any $\bm{\beta}$ between $\bm{\beta}_{0}$ and $\widetilde{\bm{\beta}}$, where $s^{(d)}(t;\bm{\beta})$ is the limiting value of $S^{(d)}(t;\bm{\beta})$, with $s^{(1)}(t;\bm{\beta})$ and $s^{(2)}(t;\bm{\beta})$ bounded, and $s^{(0)}(t;\bm{\beta})$ bounded away from $0$ for $t\in[0,\tau]$.

• (C6).

  Positive-definiteness of the matrix $\bm{h}(\bm{\beta})=\int_{0}^{\tau}v(t;\bm{\beta})s^{(0)}(t;\bm{\beta})\lambda_{0}(t)dt$, for any $\bm{\beta}$ between $\bm{\beta}_{0}$ and $\widetilde{\bm{\beta}}$, where

  $\displaystyle v(t;\bm{\beta})=\left[\frac{s_{k}^{(2)}(\bm{\beta})}{s_{k}^{(0)}(\bm{\beta})}-\left\{\frac{s_{k}^{(1)}(\bm{\beta})}{s_{k}^{(0)}(\bm{\beta})}\right\}^{\otimes 2}\right].$

Conditions (C1)-(C6) are analogous to the regularity conditions employed in Andersen and Gill (1982) and Andersen et al. (1993) to prove the consistency and asymptotic normality of estimators from classical Cox models.

Lemma 1 is a standard results showing that the CoxKL estimator $\hat{\bm{\beta}}(\eta)$ is consistent for the unique maximum for the expectation of the CoxKL objective function $\bm{\beta}^{*}(\eta)$. By Proposition 1, the CoxKL objective function (7) retains a similar form as that of a classical Cox model, thus, Lemma 1 can be easily verified following standard arguments on the asymptotic results of maximum partial-likelihood estimator for classical Cox model under standard regularity conditions (Andersen et al., 1993).

Lemma 2 shows the asymptotic distribution of $\hat{\bm{\beta}}(\eta)$. When $\bm{\Delta}=\bm{0}$, $\hat{\bm{\beta}}(\eta)$ is an unbiased estimator for $\bm{\beta}_{0}$. In contrast, when $\bm{\Delta}\neq\bm{0}$ under population heterogeneity, as $n\rightarrow\infty$ and with the order of $\eta=o(n^{-1/2})$, $\hat{\bm{\beta}}(\eta)\overset{p}{\rightarrow}\bm{\beta}_{0}$. Moreover, the CoxKL estimator $\hat{\bm{\beta}}(\eta)$ yields a considerable efficiency gain over $\hat{\bm{\beta}}$ by incorporating external risk information, which can be easily seen from the form of asymptotic variance of $\hat{\bm{\beta}}(\eta)$. We then have the following Theorem 1, which states the advantage on the estimation of the CoxKL method.

Theorem 1 implies that the CoxKL method improves the aMSE given an appropriate range of $\eta$ regardless of the potential heterogeneity. When $\bm{\Delta}=0$, $\hat{\bm{\beta}}(\eta)$ is unbiased, and it’s intuitive that $\text{aMSE}(\hat{\bm{\beta}}(\eta))<\text{aMSE}(\hat{\bm{\beta}})$ for all $\eta>0$ since the CoxKL method improves the estimation efficiency as proved in Lemma 2. Moreover, even if the external and internal information are heterogeneous and deviate from each other, there exists a range of $\eta>0$ such that the CoxKL estimator $\hat{\bm{\beta}}(\eta)$ achieves a smaller aMSE for the model regression parameters compared to the maximum partial-likelihood estimator $\hat{\bm{\beta}}$ using the internal data only. Furthermore, as mentioned in Section 2.2, we do not require that the internal data and the external risk score share the same set of covariates. Since $\widetilde{\bm{\beta}}$ can be a sparse vector (i.e. some entries of $\widetilde{\bm{\beta}}$ are zero), we allow internal data to include additional predictors.

To evaluate the performance of the proposed CoxKL method, we conducted two sets of simulation studies. In Setting I, we focused on the evaluation of the estimation performance of CoxKL, while in Setting II, we evaluated the prediction performance of CoxKL under different scenarios, including the integration with various types of external risk score models (CoxPH and XGBoost (Chen and Guestrin, 2016)) and multiple external risk scores. The tuning parameter $\eta$ was selected by the 5-fold cross validation with the V&VH cross-validated partial likelihood as the measure of model performance. In addition to the internal and the external models, we compared the performance of the CoxKL method with the stacked regression method (Debray et al., 2014), which includes the prediction of the external model as a covariate in the integration model. We evaluated the model performance on an independent testing dataset, which was simulated from the same distribution as the internal data. The simulation studies were replicated 500 times.

Prostate cancer is the most prevalent cancer and second leading cause of cancer death in American males (Siegel et al., 2022). Identifying high-risk population and offering personalized suggestions on when to initialize prostate cancer screening is important to reduce the cancer-related morbidity and mortality (Schröder et al., 2014). PHS46 has shown accurate discrimination of onset of prostate cancer in EUR males but reduced performance in non-EUR males, such as AFR males (Huynh-Le et al., 2021). To improve the risk discrimination of prostate cancer in AFR males, we applied CoxKL to develop a trans-ancestry PHS by integrating PHS46 with a small-sized internal dataset of AFR individuals.

The PHS46 model is a Cox proportional hazards model developed using a genotype dataset consisting of EUR males. By analyzing associations between genetic factors and the age at diagnosis of prostate cancer, the model includes 46 SNPs and predicts prostate cancer risk in EUR males. Based on the external validation results, the PHS46 model performs poorly with a C-index of 0.55 in AFR males (Karunamuni et al., 2021). Thus, to improve the risk discrimination for prostate cancer in AFR males, an ancestry specific PHS model is desired.

The internal dataset obtained from the Michigan Genomics Initiative (MGI) included 1,036 AFR males with both genotype and phenotype data (Zawistowski et al., 2021). The outcome was defined as the age at diagnosis of prostate cancer. For subjects who were not diagnosed with prostate cancer, the study follow-up was censored at the age of the last visit, and the overall censoring rate was $82.5\%$. Standard quality control was performed on the genotype data. Eventually, a total of 161,861 imputed SNPs passing a quality threshold of $R^{2}>0.3$, with a minor allele frequency $>0.05$, and LD pruned with $R^{2}>0.1$, $+/-500kb$ and $\text{step size}=1$ were considered in the analysis. In addition, we conducted a GWAS analysis to further identify the SNPs that may be associated with the outcome using univariate Cox regression models. Top 10,000 SNPs identified based on p-values from the GWAS results were included in the PHS model development.

With the moderate sample size, high censoring rate and high-dimensional covariate space, we considered incorporating PHS46 with the internal data by utilizing the proposed CoxKL-LASSO procedure (introduced in Section 3.2). The internal data was randomly split into training and testing sets with a ratio of 3:1. We developed PHS models on the training set and evaluated the model performance on the testing set. A total of 25 independent random splits were implemented. As comparisons, results of the PHS46 model and the internal model by applying a penalized Cox regression model with LASSO on the internal data are also presented.

As shown in Figure 4, the proposed method had the best prediction performance compared to the internal model and PHS46 model. Specifically, when the dimension of the covariate space was much larger than the sample size, it was difficult to develop a proper PHS model based on the internal small-sized data only. Therefore, the CoxKL method incorporated more information from PHS46 with an optimal $\eta$ at around $3$. Furthermore, as shown in Figure 5, the CoxKL method also achieved better risk discrimination of prostate cancer in AFR males than both the internal and PHS46 models. High-risk populations (subjects with PHS score ranging from $81-100\%$) were clearly identified by the proposed CoxKL method, which indicates that our method can be used for better personalized suggestions on when to initialize prostate cancer screening for AFR males.

In this paper, we proposed the CoxKL method to incorporate external risk scores derived from published prediction models with internal time-to-event data. Partial-likelihood-based KL information is utilized to provide an overall assessment of heterogeneity across populations, and a penalized objective function is developed to balance the contribution of different information sources. Extensive simulation studies show the advantages of the CoxKL method, including improvements on estimation efficiency and prediction accuracy.

The CoxKL method is computationally efficient and, thus, can be easily extended to integration problems under high-dimensional settings. We presented CoxKL-LASSO for developing trans-ancestry PHS prediction models in non-EUR populations. We applied our method to develop a trans-ancestry PHS model for prostate cancer by integrating PHS46, a previously published EUR-based PHS, with a small-sized AFR genotype data. The prostate cancer risk discrimination in AFR males has been substantially improved (C-Index from 0.49 [Internal] and 0.54 [PHS46] to 0.68 [CoxKL]). The proposed AFR-based trans-ancestry PHS can clearly distinguish subjects with high prostate cancer risk.

As the CoxKL method summarizes external risk information into personalized risk scores, it provides a unified framework for integrating risk information from various types of external models, including prediction models by machine learning algorithms, as long as the external models can provide personalized risk scores for each subject in the internal data. Moreover, as shown by both theoretical properties and simulation studies, the CoxKL method allows internal data include additional covariates than that used in external risk score models. The CoxKL method can also be easily extended to the integration with risk information from multiple external models, details are given in Supporting Information.

In summary, the proposed CoxKL method serves as a general and flexible integration tool for time-to-event data, with an important application in developing trans-ancestry PHS models and enhancing risk discrimination in non-EUR populations.

The authors acknowledge the Michigan Genomics Initiative participants, Precision Health at the University of Michigan, the University of Michigan Medical School Central Biorepository, the University of Michigan Medical School Data Office for Clinical and Translational Research, the University of Michigan Advanced Genomics Core for providing data storage and specimen storage, management, processing, and distribution services, and the Center for Statistical Genetics in the Department of Biostatistics at the School of Public Health for genotype data curation, imputation, and management in support of the research reported in this publication. This work was partially supported by the National Institute of Diabetes and Digestive and Kidney Diseases under Grant R01DK129539.