omerwe / polyfun

PolyFun (POLYgenic FUNctionally-informed fine-mapping)
MIT License
85 stars 21 forks source link

Some doubts about the standard error of calculating R-square using the jackknife method #197

Closed Y-Isaac closed 2 months ago

Y-Isaac commented 2 months ago

Hi, I apologize for bothering you again. I would like to follow your manuscript and apply the jackknife method to calculate whether the R² improvement of polypred compared to bolt-lmm is statistically significant. However, I have encountered some difficulties in the process. I have two questions I hope you can help me with:

1.Based on the .jk file and using Equation 1 (Figure 1), I have calculated the R² s.e. for both polypred and bolt-lmm (I believe this process is correct, maybe?).

image

I then attempted to use Equation 2 (Figure 2) to calculate the R² diff s.e. However, I found that the results do not match those in your Supplementary Tables 4-5. 1717776083751

After carefully reviewing your manuscript, I suspect that the R² diff s.e. might also be directly calculated using the jackknife method, rather than using Equation 2. Is my understanding correct? If this is the case, could you please confirm if my understanding is correct: For each PRS subset obtained using the jackknife method, we calculate the R² difference between polypred and bolt-lmm. After obtaining 200 pairs of differences, we use Equation 1 again to calculate the standard error of these differences. Additionally, the PRS subset indices must be consistent for each calculation of the difference (i.e., the column names in the .jk file)?

2.You mentioned that you used inverse-variance weighting for meta-analysis of multiple phenotypic results. I would like to confirm if the correct approach is to first use the jackknife method to obtain the R² differences and their standard errors for each phenotype under both methods (polypred and bolt-lmm), and then apply inverse-variance weighting to these differences? However, when I attempted to apply this method to the results in Supplementary Table 5, my weighted calculation did not match your reported results. For example, in Supplementary Table 5, the R² for PolyPred in non-British Europeans for the seven phenotypes are: 0.1307, 0.3823, 0.1972, 0.1679, 0.1129, 0.1127, and 0.1123. The corresponding standard errors are: 0.0104, 0.0107, 0.0197, 0.0161, 0.0047, 0.0122, and 0.0119. My inverse-variance weighted R² for these seven phenotypes is 0.148, whereas your table reports a Meta-analysis R² of 0.174. I do not understand where my calculation went wrong and would appreciate your guidance.

My statistical foundation is not very well, and some of my questions might seem a bit silly. I apologize for any inconvenience this might cause, and thank you in advance for your help!

Y-Isaac commented 2 months ago

Some new ideas about the second question: In the method, you said 'We estimated statistics (for example, relative R2) for meta-analyzed traits via an inverse-variance-weighted average, using weights inversely proportional to the s.e. of the R2 of BOLT-LMM in the target population (as estimated via genomic block-jackknife).'

But, I don't quite understand, why did you say that the weight is inversely proportional to s.e? In the inverse variance weighting method, shouldn't the weight be inversely proportional to the variance?

Moreover, the meaning here is that when I want to perform a meta-analysis on R2 diff with multiple phenotypes, the weight is not determined by R2 diff s.e., but by R2 s.e. obtained from the corresponding phenotype-population pair from bolt-lmm?

I'm really sorry for asking you so many questions, and I hope you can forgive me...

omerwe commented 2 months ago

@Y-Isaac these are good questions, and no worries.

I'll try to answer as best as I can, but given that I've worked on this more than five years ago I can't guarantee to get all the details right.

  1. Yes, you basically need to compute the R2 difference for each jackknife estimate (each one corresponds to 199 genomic regions). Then you need to compute the standard deviation of these 200 different estimates (with zero degrees of freedom) and multiply by sqrt(199).

In Python, this is written as follows: np.std(arr_estimates, ddof=0) * np.sqrt(199) where arr_estimates is an array of the 200 different estimates (each one from a different jackknife estimate).

  1. Here is Python code that shows how we obtain a value of 0.174: (if you're more comfortable using R, ChatGPT can translate this to R quite well I believe).
import numpy as np
arr_r2 = np.array([0.1307, 0.3823, 0.1972, 0.1679, 0.1129, 0.1127, 0.1123])
arr_se = np.array([0.0104, 0.0107, 0.0197, 0.0161, 0.0047, 0.0122, 0.0119])
arr_w = arr_r2 / arr_r2.sum()
se = np.sqrt(np.mean(arr_r2**2 / arr_w)) / np.sqrt(len(arr_r2))

Hope this is clear, please let me know if not!

Y-Isaac commented 2 months ago

@omerwe Thanks!! I fully understand the first question. I also understand how 0.174 was obtained regarding the issue of meta. But the only question I don't quite understand, it doesn't seem like it was obtained using the inverse-variance-weighted average?

omerwe commented 2 months ago

@Y-Isaac you're correct, this is a mistake in the text of the Methods section, my apologies. the words "the standard error of" are mistaken... The correct text should read "using weights inversely proportional to the R2 of BOLT-LMM in the target population". The idea was to upweight traits for which we can get strong prediction performance (otherwise if we had a trait with zero predictive capability, the s.e. for that trait would be zero and thus we would place all weight on that trait, although this is the least interesting trait).

Y-Isaac commented 2 months ago

@omerwe I agree with you! The phenotypes with higher heritability have more causal loci, so we should give higher weights, which is very coherent.