2.1 Setup

Let’s assume that we have a case-control study of a multicategory outcome Y and a set of covariates X. The outcome Y takes integer value from 0 to K, with Y = 0 for controls, and Y > 0 for other K different outcomes. For example, in the motivation example we consider a case-control study of NHL with four subtypes. We can denote individual-level data from this study (called internal study) as {X_i, Y_i, i = 1, …, n}, with n_k subjects having outcome Y = k, and . In the following discussion, we refer Y = 1, …, K as different disease subtypes. We assume the following prospective PLR model as the underlying model for the effect of X on Y, (1) where ω_k is the intercept and M_k(X; θ_k) is a given function, such as M_k(X; θ_k) = X^⊤θ_k, with θ_k being the vector of parameters of interest. In the NHL example, we have M_k(X; θ_k) = θ_k ⋅ S(X), with S(X) being the PRS. S(X) is defined as ∑_i w_iX_i, with X_i being the genotype at the k-th SNP, and w_i being the weight defined by previous GWAS. Let be the collection of all coefficients. We are interested in estimating θ or comparing θ_k among different disease subtypes.

Beside data from the internal study, we have summary data from multiple external studies. In the NHL example, we consider the PRS defined by a set of 21 SNPs [39]. For each of those SNPs, we have SNP-level summary statistics from seven external studies. The SNP-level summary data is the estimated SNP’s marginal association (the regression coefficient) with one NHL subtype based on a standard logistic regression model (called the working model). Therefore, for an external study consisting of three different NHL subtype cases and controls, its summary data consists of 21 × 3 estimated regression coefficients, each of which is derived from a separate working model. The goal is to fit the PLR model (1) using individual-level and summary data and to utilize the fitted model for more efficient inference of θ_k and disease subtype heterogeneity. An illustration of the proposed integration framework is given in Fig 1. Before considering this complicated scenario, we first present the methodology assuming that there is one external study, and that the summary data is derived from a single working model. We will describe the method under the more general setting in Section 2.7.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Illustration of the PolyGIM framework.

We consider a setting similar to the NHL study. The outcome Y takes values from {0, 1, 2, 3, 4}. The set of covariates X comprises measures on 21 genetic markers (SNPs). We have individual-level data (Y, X) from the internal study and summary data from three external studies. Summary data consists of coefficient estimates from marginal logistic regression models established by the k-th external study (k = 1, 2, 3), with the binary outcome (D^(k)) of each external model defined by Y. The goal is to fit the polytomous regression model for the PRS effect given in the red box using individual-level and summary data.

https://doi.org/10.1371/journal.pcbi.1011236.g001

More specifically, we consider summary data consisting of regression coefficient estimates derived from a working model of a binary outcome D. The binary outcome D is derived from of the multicategory outcome Y. The working model is assumed to be a standard logistic regression model. Several versions of D can occur in practice. For example, a working model can treat several disease subtypes as one broad disease group, that is, D = 0 if Y = 0, and D = 1 if Y ∈ {c₁, …, c_L}, with 1 ≤ c_l ≤ K, 1 ≤ l ≤ L, where L is the total number of disease subtypes under study. We call this type of working model the grouped case-control (GC) model. We further assume the number of cases within each subtype group is known in the GC model. In some cases, such as in the NHL study, the GC model can consider just one subtype (i.e., L = 1). Another version of D arises from a case-only study, in which a logistic regression model is used to evaluate whether a risk factor has the same effect on two disease subtypes c₁ and c₂. We call this the case-case (CC) comparison model, with D = 0 for Y = c₁, and D = 1 for Y = c₂.

2.2 Likelihood for the internal study

Observations from the internal case-control study can be retrospectively sampled conditioning on the outcome Y from a source population. For some other case-control studies, cases are sampled from P(X|Y > 0), with their disease subtypes classified after the enrollment. As long as the sampling criterion is independent of X, we can consider both types of case-control studies as stratified samples from P(X|Y = 0) and P(X|Y = k), k = 1, …, K.

We consider the retrospective likelihood for the case-control data as this framework is convenient for incorporating external summary data. Let π_k = P(Y = k) and , k = 1, …, K. We define the vector τ = (τ₁, …, τ_K)^⊤. By Bayes’ rule, model (1) has the following equivalent retrospective representation, (2) where Δ_k(X; ξ) = exp{τ_k + M_k(X; θ_k)}, k = 1, …, K, with ξ = (τ^⊤, θ^⊤)^⊤. We denote ξ* = (τ*^⊤, θ*^⊤)^⊤ as the true population value of ξ.

We use (2) to form the likelihood of observing X given Y. Denote as the empirical distribution of P(X|Y = 0) supported on samples from the internal study, and as the indicator function. We can estimate ξ by maximizing the following empirical log-likelihood, (3) subject to constraints and for k = 1, …, K. Using the method of Lagrange multipliers to profile out , we can infer ξ by finding the stationary point of the following profile log-likelihood, (4) where ρ_k = n_k/n₀ for k = 1, …, K. Let be the maximum likelihood estimate (MLE) based on the internal data. When M_k(X; θ_k) = X^⊤θ_k(k = 1, …, K), we can follow [40] to show that the empirical likelihood estimate of θ derived from (4) is exactly the same as the maximum likelihood estimate based on the standard prospective likelihood function specified by (1).

2.3 Properties of summary data

Here we will show the property of the summary data and its relationship with θ, the parameter of interest. We assume that the external study consists of observations {X_i, D_i, i = n + 1, …, n + N}, representing stratified samples from P(X|D = 0) and P(X|D = 1). Suppose a standard logistic regression model is used as the working model to study the effect of X on D. We can represent this model by its equivalent retrospective formation as, (5) with P(X|D = 1) and P(X|D = 0) being distributions of X in the two groups. This model can be misspecified and thus inconsistent with (2). In model (5), all unknown parameters are divided into two parts, and β. We use α to represent the set of nuisance parameters whose estimates are not given as part of the summary data. Note that the intercept term α₀ is always assumed to be part of nuisance parameters. The summary data only consists of the estimate of β.

Let N₀ and N₁ be sample sizes in groups D = 0 and D = 1, respectively. Based on the working model (5), the log-likelihood function of the external study is where ρ = N₁/N₀ and . is the solution of the following estimating equation, (6) Note that is the same as the estimate from the standard package for the logistic regression model. For simplicity, let and . According to [41], is a consistent estimate of (α*, β*), which satisfies the following stochastic constraint equation (7) where and are expectations over P(X|D = 0) and P(X|D = 1), respectively.

The estimate depends on how D is derived from Y. As mentioned in Section 2.1, two types of binary outcomes are typically encountered. One is used in the CC model, and the other is used in the GC model.

The CC model compares the distribution of X between disease subtype groups c₀ and c₁, c₀ ≠ c₁ ∈ {1, …, K}. Its binary outcome is defined as D = D_CC, with D_CC = 0 if Y = c₀, and D_CC = 1 if Y = c₁. The CC model is commonly used for the study of disease subtype heterogeneity. Let μ = (τ^⊤, θ^⊤, α^⊤, β^⊤)^⊤ be the vector of all unknown parameters, and μ* = (τ*^⊤, θ*^⊤, α*^⊤, β*^⊤)^⊤ be their true population values. Notice that and due to (2), we can express (7) as (8) with where is the expectation of over P(X|Y = 0). The asymptotic distribution of is given by (9) with A and B shown in S1 Appendix. Based on (9), we know , where Σ₀ = (A⁻¹B A⁻¹)_{β β} is the submatrix of A⁻¹B A⁻¹ corresponding to β. Notice that A and B are determined by the expectation defined by P(X|Y = 0). We will provide an empirical likelihood based estimate of P(X|Y = 0) later. According to (9), is a consistent estimate of β*, whose relationship with θ*, the parameter of interest, is governed by (8). Therefore, summary data contains information on θ*.

In the GC model, one or several disease subtypes are considered as one broad case group. The binary outcome is defined as D = D_GC, with D_GC = 0 if Y = 0, and D_GC = 1 if Y belongs to {c₁, c₂, …, c_L}. We further assume that the proportion of cases belonging to each disease subtype k is known, and denote it as q_k(k = 1, …, K). Note that some q_ks can be zero for studies that do not collect cases with certain disease subtypes. Based on (2), we have Then (8) still applies with g(X; μ) defined as Furthermore, the asymptotic distribution of takes the same form as (9), where A and B are defined in S1 Appendix. We can obtain similarly as for the CC model.

In the following discussion, we will classify summary data into two distinct types, the regular and irregular summary data. Regular summary data has its corresponding g(X; μ*) to be not constantly 0 over X. Irregular summary data is the one with g(X; μ*) ≡ 0 for all X. The validity of the proposed procedure for integrating regular summary data requires some standard regularity conditions, which would not hold if g(X; μ*) ≡ 0. Therefore, a different integration procedure is needed for irregular summary data.

2.4 Procedure for integrating regular summary data

Here we extend the generalized integration model (GIM) of [38] to fit a PLR model by integrating individual-level and summary data. We call the new procedure PolyGIM.

The log-likelihood for the internal study is given by (3). The log-likelihood function of the summary data can be represented as . Since we do not know Σ₀, we can replace it with a known matrix V (e.g., the identity matrix I) as the starting point. By combining the two into a joint (pseudo) log-likelihood function, we can estimate μ via solving the following optimization problem over , (10) The last constraint equation in (10) is from (8) with a specific g.

We employ the Lagrange multiplier approach to solve (10). The Lagrange function can be written as, where κ, λ and ν are the Lagrange multipliers. It can be seen that κ = 1 and Let η = (λ^⊤, ν^⊤, μ^⊤)^⊤ be the vector of all variables. Therefore, the profiled log-likelihood function can be written as (11) Hence, solving the original problem (10) can be translated into finding the solution of a set of score equations (shown in S1 Appendix). Then we can apply the Newton-Raphson algorithm to to find the solution . In S1 Appendix, we show in Lemma 1 that under some regularity conditions is a consistent estimate of η* = (λ*^⊤, ν*^⊤, μ*^⊤)^⊤, with and . In particular, regularity condition C4 requires g(X; μ*) to be not constantly 0. More intuitively, if g(X; μ*) ≡ 0, the Lagrange multiplier ν in (11) is not identifiable. Consequently, the procedure based on (11) is only applicable to summary data that satisfy those regular conditions. We refer to such summary data as regular summary data, while summary data not meeting those conditions are called irregular summary data. Summary data derived from a GC model in general is regular. It is also regular if it is derived from a CC model that does not consider the same set of covariates as the one by the underlying PLR model, as in this setting g(X; μ) can not be constantly 0 for any μ.

Based on the estimate , we can obtain the estimated empirical distribution of P(X|Y = 0) as (12) Furthermore, we can estimate A and B in (9) by calculating the expectation over P(X|Y = 0) with . An example is shown in S1 Appendix. Therefore, can be estimated as , with given by (13) which is the submatrix of corresponding to β.

Here is the summary of theoretic properties of the estimate.

Proposition 1 Under the regularity conditions given in S1 Appendix, assuming that N/n → γ, N₁/N₀ → ρ and n_k/n → ρ_k for k = 1, …, K as n → ∞. We have is asymptotically normal, and its asymptotic variance-covariance matrix attains its minimum (in term of positive semidefinite) at V = Σ₀. Furthermore, is asymptotically more efficient than the internal data based MLE .

Proofs are given in S1 Appendix. The optimality of estimate still holds when we replace Σ₀ with its consistent estimate . We propose the following iterative Algorithm 1 to obtain this optimal estimate.

Algorithm 1 Algorithm for PolyGIM

1: Based on the internal study, we fit the PLR model to obtain a consistent estimate of ξ, and fit the working model to obtain a consistent estimate of (α, β). We denote these estimates as , where and . Then an initial estimate of is computed according to (13) at .

2: Resolve the score equation with . Let the estimates be .

3: Estimate the empirical probability via (12).

4: Update via (13).

5: Repeat Steps 2 to 4 until is converged.

We use the following strategy to choose the initial point in Step 1. We obtain by fitting the PLR model with the internal data. Since we formulate the PLR model with the empirical likelihood representation (4), we adjust estimates of intercept terms from the standard package for the PLR model by subtracting log(ρ_i), i = 1, …, K. To obtain the initial estimate of (α, β), we first adjust for the sample size difference between the internal and external studies by assigning each subject from the internal study an appropriate weight. Then we fit a weighted logistic regression model to obtain the initial estimate of (α, β). More specifically, suppose that the external study fits a GC model based on N₀ controls and N₁ grouped cases, with r_i proportion of them having disease subtype i(i = 1, …, K). We assign N₀/n₀ as the weight for each control and (N₁r_i)/n_i as the weight for each subtype i case in the internal study. We can obtain by fitting a weighted logistic regression model with those weights. Again, we need to adjust the estimate of the intercept term due to the use of the empirical likelihood representation.

2.5 Integrating irregular summary data

When summary data is derived from a working CC model that is consistent with the underlying risk model, it may become irregular with g(X; μ) ≡ 0. More specifically, if M_k(⋅) = X^⊤θ_k in the underlying risk model (2) and m(⋅) = X^⊤β in the working CC model (5) with two subtypes c₀ and c₁, based on the definitions of {φ₀, φ₁, Δ_k, δ} in Section 2.2 and 2.3, it can be shown that So we have g(X; μ ≡ 0 if we let (14) Notice that the true value of μ satisfies (14) due to the consistency between the working model and the underlying risk model. Under these constraints (14), we can eliminate from (10) and use the Lagrange multiplier approach to solve a modified version of (10). The resultant profile log-likelihood function can be written as (15) We can obtain the estimate based on (15). Under regularity conditions given in S1 Appendix, we can show that is consistent for any given positive definite V.

The asymptotic distribution of and the optimal choice of V are summarized by the following result, with proofs given in S1 Appendix.

Proposition 2 Under model (15) and regularity conditions given in S1 Appendix, we have is asymptotically normal, and its asymptotic variance-covariance matrix attains its minimum at V = Σ₀. In particular, the asymptotic variance-covariance matrix of has the following form, where is defined in S1 Appendix, and with and 1 being a vector of 1’s.

Similar to the regular PolyGIM procedure, we can obtain by an iterative algorithm. Even though we use different procedures for integrating regular and irregular summary data, in the following discussion we still call them the PolyGIM procedure when there is no confusion.

We can also use a restricted MLE (RMLE) approach to incorporate this irregular summary data. Let be the solution of the following constraint optimization problem, (16)

This setup is different from the standard RMLE, as in the constraint equation has variability. We need to account for this uncertainty when estimating the variance-covariance matrix of . In S1 Appendix, we prove the following result that shows the PolyGIM estimate is more efficient than other considered estimates.

Proposition 3 The estimate based on (15) is asymptotically more efficient than both the internal data based MLE and the RMLE .

2.6 Test for disease subtype heterogeneity

As mentioned in the Introduction, researchers are often interested in testing whether a risk factor X_j has the same effect on different disease subtypes, with the null hypothesis being H₀ : θ_j1 = θ_j2 = … = θ_jK. We can use PolyGIM to combine data from multiple sources for a more efficient test.

By Proposition 1, , the estimated coefficients corresponding to X_j, has a multivariate normal distribution , where is extracted from the estimated variance-covariance matrix of . So its log-likelihood function can be written as Under the null, follows . The MLE of a under the null can be derived as . Thus, we can construct the following likelihood ratio test, Under the null, this test follows a chi-square distribution with K − 1 degrees of freedom.

2.7 Summary data from multiple models

So far, we have considered summary data from a single working model based on one external study. Here we show how to incorporate summary data from multiple working models, some of which can be fitted with overlapping samples. For example, in the NHL example, we have summary data from seven independent external studies. From one study consisting of cases of three NHL subtypes, we have three summary statistics on each SNP, which are estimated from three GC models sharing a common set of controls. For notation simplicity, we focus on regular summary data in the following discussion.

First, we provide some theoretical insights on the PolyGIM procedure with multiple summary data. Suppose that we are given the summary data , from two external models. Let α_i (i = 1, 2) be the corresponding nuisance parameters. Similar to Section 2.3, we can establish the following asymptotics where are the true values. The specific form of Σ depends on the two external models and whether there are overlapped samples used for fitting them. For example, if the two external models are fitted with data from two different studies, we have Σ₁₂ = 0. If overlapped samples are used for fitting the two models, Σ₁₂ ≠ 0. Later we will provide more details on how to estimate Σ.

Let be the optimal PolyGIM estimate of ξ using summary data , be the optimal estimate using summary data . We can show theoretically that it is always beneficial by using more summary data (Proposition 4), see S1 Appendix.

More technical details on integrating summary data from multiple studies are presented in S1 Appendix.

3.1 Simulation studies with summary data from one external study

To verify the theoretic properties of the proposed PolyGIM method, we first considered a simple scenario when summary data comes from one external study. Assume that X = (X₁, X₂) was a vector of two binary biomarkers and that the outcome Y had three classes, with 0 for the control group and 1/2 for the two disease subtypes. In the study population, the two biomarkers were correlated, with the joint probability of (X₁, X₂) = (0, 0), (0, 1), (1, 0), and (1, 1) specified as 0.28, 0.12, 0.18 and 0.42, respectively. The true disease risk model was chosen as We let (ω₁, ω₂) = (−0.2, 0.1) and (θ₁₁, θ₁₂, θ₂₁, θ₂₂) = (1.5, −1.0, −0.5, 1.2). From this population, case-control studies were retrospectively generated conditioning on the outcome Y, with the internal study consisting of 2000 controls, 500 subtype 1 cases and 500 subtype 2 cases, and the external study consisting of 1000 controls, 600 subtype 1 cases and 400 subtype 2 cases.

We considered two classes of summary data, one from the CC model and the other from the GC model. The binary outcome D used in a CC model was defined as D = 0 for subtype 1 and D = 1 for subtype 2. Summary data from two types of CC working models (CC1 and CC2) were studied. The CC1 working model only included covariate X₁ and thus was inconsistent with the underlying risk model. It generated regular summary data , which was the estimated coefficient of X₁. The CC2 working model was consistent with the underlying risk model and included covariates (X₁, X₂). It generated irregular summary data , which were estimated coefficients of X₁ and X₂.

The GC model merged the two disease subtypes into one group and compared it with the control group. Two GC working models (GC1 and GC2) were used to generate regular summary data. GC1 model included only X₁ as a covariate and generated summary data . GC2 model had X₁ and X₂ as covariates and generated summary data .

For each simulated dataset (including an internal study and summary data), we applied three versions of PolyGIM, one based on the most optimal estimate (called GIM_opt), the others based on the estimate with two options of V, one called GIM_I with V = I, and the other call with V = V_σ, a diagonal matrix that has each diagonal element being , with being the variance of the i-th summary statistic from the summary data . We applied different types of PolyGIM procedures depending on whether the summary data was regular or irregular. As a comparison, we also analyzed the internal study using the standard PLR model (called MLE_int). We simulated 2000 datasets under each scenario to evaluate the performance of all considered methods. Tables 1 and 2 summarized simulation results in situations when summary data is generated from CC and GC models, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Simulation results in situations when summary data is derived from one external study based on case-case (CC) comparison models.

https://doi.org/10.1371/journal.pcbi.1011236.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Simulation results in situations when summary data is derived from one external study based on grouped case-control (GC) models.

https://doi.org/10.1371/journal.pcbi.1011236.t002

Table 1 for the CC models shows that all considered methods have expected performances, with unbiased estimates and their estimated standard errors matching well with corresponding empirical standard errors. When using irregular summary data , we can see that GIM_opt provides more efficient estimates of {θ_k1, θ_k2, k = 1, 2} compared to GIM_I, and MLE_int. These findings align with the conclusions of Propositions 2 and 3. If only is used as summary data, GIM_opt and are more effective for estimating coefficients of X₁, whereas all methods have similar efficiency levels for estimating coefficients of X₂. Notably, results from Table 1 reveal that has a similar level of efficiency as GIM_opt when integrating only one summary statistic, but it becomes less efficient when integrating two correlated summary statistics. This efficiency loss arises because selects the matrix V = V_σ in (10), assuming all summary statistics to be independent. This chosen V_σ is suboptimal when summary statistics are correlated, since the optimal matrix should be a consistent estimate of the variance-covariance matrix for these statistics. However, when the summary statistics are indeed independent, V_σ can serve as a satisfactory practical approximation for the optimal matrix. Although they differ theoretically in their diagonal terms when external models are mis-specified (e.g., the CC1 model), this discrepancy is more theoretical than practical in our context, as evident from S1 Table. From Table 2 and S2 Table. for the GC models, we can reach similar conclusions as those from Table 1 and S1 Table.

3.2 Simulation studies with summary data from multiple external studies

We conducted additional simulation studies under a more complex setting to mimic the NHL study. We considered 21 SNPs X = (X₁, …, X₂₁), and the outcome Y had five classes, with Y = 0 representing controls and Y = 1 to 4 representing four different disease subtypes. We selected the 21 SNPs identical to those used in the NHL example, with the exception that we assumed they were independent to simplify the simulation procedure. The characteristics of these SNPs are summarized in S3 Table. The true underlying risk model was defined as with the PRS . We considered the following two set of parameters for the true model:

1. Null PRS Model: (ω₁, ω₂, ω₃, ω₄) = (−3.9, −4.1, −3.6, −3.8) and (θ₁, θ₂, θ₃, θ₄) = (0, 0, 0, 0).
2. Alternative PRS Model: ω_i = −3.8, i = 1, …, 4, and (θ₁, θ₂, θ₃, θ₄) = (0.019, 0.092, −0.12, 0.047).

Both models were designed with intercept terms to ensure the rarity of each disease subtype in the study population (each with a prevalence of less than 2%). For Alternative PRS Model, the effects of the PRS were chosen to match those observed in the NCI study of the NHL example.

In practice the weights (w_i) used in PRS are estimated with uncertainty from other studies, leading to a PRS with measurement error. We aim to evaluate the performance of procedures under consideration while accounting for the measurement error. To accomplish this, we generated estimates () of the true weights (w_i) for each simulated dataset, assuming that the true weights w_i were identical to those used in the real NHL example (S3 Table). We randomly generated from a normal distribution , where se_i was the standard error reported in the published GWAS from which the SNP was genome-wide significantly detected (see S3 Table). We used to calculate the PRS with measurement error, denoted as , and varied the level of measurement error by choosing the scaling factor c from {0, 1, 16}. In the analysis of each simulated dataset, we used instead of S(X) to reflect the measurement error.

First, we considered performance of MLE_int and GIM_opt under the null model. We used summary data from External Studies 1–5 listed in Table 3. All these five external studies, as well as the internal study, were generated from the same source population (the internal study population). S3 Table provides the allele frequency of the effect allele “1” for each of the 21 SNPs in the study population. For each external study, we used a logistic regression model to estimate the marginal effect of each SNP (i.e., the regression coefficient) on the risk of a specific disease subtype. The summary data consisted of these coefficient estimates for all considered 21 SNPs. We simulated 2000 datasets under Null PRS model and analyzed each dataset with MLE_int and GIM_opt. To assess the impact of measurement error, we employed three sets of PRS in the analysis of each simulated dataset, including the PRS without measurement error (i.e., ), the PRS with measurement error at the same level as the real NHL study (i.e., ), and the PRS with elevated measurement error (i.e., ). The simulation results, presented in Table 4, indicate that the level of uncertainty in PRS does not have an impact on the statistical properties of MLE_int and GIM_opt. This is expected since likelihood models for both procedures remain valid under the null model (i.e., θ = 0), even when using PRS with high levels of measurement error. Consequently, the estimate of the effect of the PRS remains consistent. Additionally, we evaluated the type I errors of GIM_opt at various significance levels, as shown in S4 Table, and generated the Q-Q plots of estimated Z-scores for the PRS effect on different subtypes in S1 Fig. These results demonstrate that GIM_opt has well calibrated type I errors and P-values.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Simulation results in situations when summary data is derived from one external study based on grouped case-control (GC) models.

https://doi.org/10.1371/journal.pcbi.1011236.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Simulation results on the impact of measurement error under the null PRS model.

https://doi.org/10.1371/journal.pcbi.1011236.t004

Next, we assessed the performance of the considered procedures under Alternative PRS Model with non-zero PRS effects, and we summarized the simulation results in Table 5. Here are some notable observations. Firstly, when the PRS has either no or relatively low measurement error, both MLE_int and GIM_opt exhibit desirable statistical properties in terms of consistency, the accuracy of standard error estimation, and 95% confidence interval coverage probability. In these cases, GIM_opt is more efficient than MLE_int. Secondly, when the PRS has increased measurement error with , estimates obtained using MLE_int and GIM_opt become inconsistent, with a noticeable increase in bias.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Simulation results on the impact of measurement error under the alternative PRS model.

https://doi.org/10.1371/journal.pcbi.1011236.t005

We further evaluated how the sample sizes of external studies affect the efficiency of GIM_opt by increasing the sample size of the five external studies by 5 and 10 times, with a focus on PRS that have no or relatively low measurement error (i.e., ). The results, which are summarized in Table 6, demonstrate that the efficiency of GIM_opt improves as the sample sizes of external studies increase. This same pattern is observed when using an internal study with only 10% of the original sample size (see S5 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Simulation results on the impact of sample sizes of external studies under the alternative PRS model.

https://doi.org/10.1371/journal.pcbi.1011236.t006

Next, we conducted simulations to compare strategies for integrating summary data when there were differences in the joint distribution of X between the internal and external study populations. Seven external studies were considered, as listed in Table 3. The first five were generated from the same source population as the internal study, while External Studies 6 and 7 were generated from two different external study populations. The allele frequency of allele “1” in External Study 6 was 0.2 lower than that in the internal study population for each of the first 10 SNPs listed in S3 Table, with a lower bound of 0.05. Similarly, in External Study 7, the allele frequency of allele “1” was 0.2 lower than that in the internal study population for each of the last 11 SNPs listed in S3 Table, with a lower bound of 0.05. All other SNPs in these two external studies had the same distribution as in the internal study population. It is worth noting that the discrepancy in the SNP distribution between the two studies is substantially higher than the discrepancy observed in the NHL study. We examined three strategies for integrating summary data from the seven external studies. The first strategy utilized summary data from only the first five external studies. The second strategy employed complete summary data, including summary statistics on the 21 SNPs across all seven studies. The third strategy used summary data from the first five external studies and partial summary data from the remaining two external studies, which included summary data only on SNPs having the same allele frequency as in the internal study population. This excluded summary data on the first 10 SNPs from External Study 6 and the last 11 SNPs from External Study 7. To investigate the impact of varying allele frequencies between the internal and external study populations, we conducted additional simulations. Specifically, we increased the sample sizes of External Studies 6 and 7 by a factor of five to further explore this effect.

Table 7 summarizes the results of our simulation study, from which several observations can be made. The GIM_opt method, when using additional partial summary data from External Studies 6 and 7, shows desirable statistical properties and is more efficient than using only summary data from the first five external studies. This advantage becomes more pronounced as we increase the sample size of the two additional external studies. In contrast, GIM_opt using complete summary data could lead to erroneous standard error estimates, especially when the sample size is increased in External Studies 6 and 7. This is an expected outcome, as the PolyGIM procedure is tailored to the scenario where the internal and external studies are conducted in the same source population. Regarding the use of partial summary data, we can show, using arguments similar to those presented in [42], that it is valid to use partial summary data on a subset of SNPs if the following conditions are met. First, the joint distribution of this subset of SNPs in the external study population must be the same as in the internal study population. Second, these SNPs must be independent of the remaining SNPs in both study populations. Third, the disease prevalence must be relatively rare in both populations so that the SNP distribution is similar between the control group and the general population. Our simulation study demonstrated that the use of partial summary data is valid and beneficial when these underlying assumptions are satisfied.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 7. Simulation results comparing different summary data integration strategies under the alternative PRS model.

https://doi.org/10.1371/journal.pcbi.1011236.t007

Finally, we carried out supplementary experiments to assess the computational efficiency and memory requirements of the PolyGIM package. In these experiments, we focused on the internal study and the external studies listed in Table 3, and analyzed a PRS with an increasing number of SNPs, ranging from 21 to 105. We considered summary data derived from the first 1, 3, or 5 external studies provided in Table 3. The computational time and memory requirements were recorded over 100 replications and are presented in S6 and S7 Tables. These tables offer a glimpse into the performance of the PolyGIM package under various conditions.

3.3 Real data application: The NHL study

NHL is the most common hematological malignancy with many subtypes with distinct molecular and clinical features. It has been hypothesized that various NHL subtypes and other lymphoid malignancies might share some degree of genetic susceptibility. Recently, [39] used eight GWAS within the InterLymph Consortium to study the genetic heritability in four major NHL subtypes, including CLL, DLBCL, FL, and MZL. To explore pleiotropy between NHL subtypes and other lymphoid malignancies (e.g., acute lymphoblastic leukemia and Hodgkin lymphoma), they generated PRS using SNPs that had been established as being associated with each lymphoid malignancy and tested their associations with risk for the four NHL subtypes.

We utilized PolyGIM to analyze the project described by [39], using the PRS for Hodgkin lymphoma (HL) as our primary example. This PRS was derived from 21 HL-associated SNPs that were selected by [39] due to their identification as genome-wide significant SNPs, each with a P-value of less than 5 × 10⁻⁸ from seven previously published GWAS [43–49]. Each SNP was selected as an index SNP to represent a nearby gene, and they were mostly independent of each other, with only two pairs having R² > 0.01 (one pair with R² = 0.024 and the other with R² = 0.038. The PRS was calculated as , where X_j represents the genotype, and w_j (weight) is the estimated effect by the j-th SNP on HL. Further information on the 21 SNPs and their estimated weights used in the PRS calculation can be found in S3 Table.

[39] collected data from eight GWAS with European ancestry, including six US-based studies and two European studies. Table 3 shows the sample sizes of cases with specific NHL subtypes and controls from each study. S2 Fig displays the minor allele frequency (MAF) of each of the 21 SNPs in the control groups of the eight studies. The plot exhibits a consistent distribution of MAF across the studies, with a maximum range of approximately 0.1. This suggests that it is reasonable to assume that all eight studies were conducted on the same population.

For our analysis, we obtained individual-level data from the US-based NCI study, which collected cases from all four NHL subtypes and controls, and treated it as the internal study. We obtained SNP-level summary statistics from the other seven studies (external studies), consisting of each SNP’s marginal effect (i.e., the estimated regression coefficient) on one NHL subtype. We assumed the following PLR model to assess the effect of PRS on the four NHL subtype, We utilized PolyGIM to fit the model in two different ways. The first approach involved integrating summary data solely from the five US-based studies, while the second approach integrated summary data from all seven studies, which included the five US-based studies and two European studies.

Table 8 summarizes the estimates of the effects of PRS on the four subtypes of NHL. The two PolyGIM approaches provide consistent results, indicating that the PRS has a significant effect on DLBCL and FL subtypes. Using summary data from all seven studies yields slightly more significant results compared to using only summary data from the five US-based studies, as the two European studies contribute additional cases on DLBCL and FL subtypes. It is interesting to notice that the PRS is positively associated with the risk of DLBCL but inversely associated with the risk of FL. This observation indicates that SNPs that increase the risk of HL tend to be associated with an increased risk of DLBCL but a reduced risk of FL. On the other hand, the PRS has no significant effect on the risk of CLL or MZL, suggesting that it is less likely to have SNPs with pleiotropic effects on both HL and CLL/MZL. In Table 9, we compared PRS effects on each pair of NHL subtypes. A global test for disease subtype heterogeneity based on the one given in Section 2.6 is significant, with P-value = 6.97 × 10⁻¹¹ based on PolyGIM using summary data from all seven external studies. We also present results based on the PLR model fitted with the internal study (the NCI study). As expected, PolyGIM estimates are more efficient than internal study-based MLE since summary data from external studies provide additional helpful information. However, the overall improvement is somewhat limited as the internal study has a much larger sample size in controls and each NHL subtype group.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 8. Estimated PRS effects on four NHL subtypes based on the NHL study.

https://doi.org/10.1371/journal.pcbi.1011236.t008

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 9. Disease subtype heterogeneity testing P-values based on the NHL study.

https://doi.org/10.1371/journal.pcbi.1011236.t009

As an experiment, we reduced the internal study sample size by randomly removing 2/3 of the control and each subtype group. Results based on this downsized internal study and the original summary data are given in S8 and S9 Tables. Compared with the results in Tables 8 and 9, the advantage of PolyGIM over MLE_int becomes more evident. Finally, we considered the following model by allowing the PRS had a nonlinear effect on each subtype, By applying PolyGIM, we found that the PRS had no significant nonlinear effect on any subtype, with P-values associated θ_k2 all being larger than 0.05.