error in rg

[marionabustamante]

I am running rg and I have the following ERROR. However I can run each of the datsets with another one...

Call:
./ldsc.py
--ref-ld-chr eur_w_ld_chr/
--out d1y_DBP
--rg /home/mbustamante/ldsc/diarrhea/d1y.sumstats.gz,DBP.sumstats.gz
--w-ld-chr eur_w_ld_chr/

Beginning analysis at Tue Dec 15 16:09:48 2015
Reading summary statistics from /home/mbustamante/ldsc/diarrhea/d1y.sumstats.gz ...
Read summary statistics for 1148956 SNPs.
Reading reference panel LD Score from eur_w_ld_chr/[1-22] ...
Read reference panel LD Scores for 1293150 SNPs.
Removing partitioned LD Scores with zero variance.
Reading regression weight LD Score from eur_w_ld_chr/[1-22] ...
Read regression weight LD Scores for 1293150 SNPs.
After merging with reference panel LD, 1147988 SNPs remain.
After merging with regression SNP LD, 1147988 SNPs remain.
Computing rg for phenotype 2/2
Reading summary statistics from DBP.sumstats.gz ...
Read summary statistics for 1217311 SNPs.
After merging with summary statistics, 1147988 SNPs remain.
973560 SNPs with valid alleles.
ERROR computing rg for phenotype 2/2, from file DBP.sumstats.gz.
Traceback (most recent call last):
File "/home/mbustamante/ldsc/ldscore/sumstats.py", line 343, in estimate_rg
rghat = _rg(loop, args, log, M_annot, ref_ld_cnames, w_ld_cname, i)
File "/home/mbustamante/ldsc/ldscore/sumstats.py", line 465, in _rg
intercept_gencov=intercepts[2], n_blocks=n_blocks, twostep=args.two_step)
File "/home/mbustamante/ldsc/ldscore/regressions.py", line 699, in init
np.multiply(hsq1.tot_delete_values, hsq2.tot_delete_values))
FloatingPointError: invalid value encountered in sqrt

Summary of Genetic Correlation Results
p1 p2 rg se z p h2_obs h2_obs_se h2_int h2_int_se gcov_int gcov_int_se
/home/mbustamante/ldsc/diarrhea/d1y.sumstats.gz DBP.sumstats.gz NA NA NA N A NA NA NA NA NA NA

[rkwalters]

Hello,
Thanks for including the full error log, it's very helpful.

This error appears to be caused by negative h2 estimates for one of the datasets when running LD score regression with the set of SNPs overlapping between the two datasets. It's likely that the problematic dataset has a very small (but positive) h2 estimate in your other analysis where there's a somewhat different SNP set, and shifts to negative h2 here in the available overlapping SNPs.

@tpoterba: Can we look at improving error messaging for this case? There's a check on the sign of h2 upstream, but it only checks the univariate h2 estimate and not h2 in overlapping SNPs.

Thanks,
Raymond

[bulik]

I think the issue is in the jackknife delete-a-block-values for sqrt(h21*h22). Sometimes one of the delete-a-block jackknife values for one of the traits will be negative even when both point estimates of h2 are positive (typically this means that one of the h2 estimates is positive but close to zero).

Here are the lines throwing the error; the software solution is to check that all entries of hsq1.tot_delete_values and hsq2.tot_delete_values are positive. I haven't thought about the best statistical solution. In principle, the distribution of a ratio will be very skew when the denominator has a lot of mass close to zero, so the mean and SE are poor summary statistics. On the other hand, the SE will be really big, which suffices to convey the key point, namely that the results are too noisy to be useful.

@hilaryfinucane what do you think?

[hilaryfinucane]

Moving ldsc to bcc.

I think the best solution is to return NA if any of the delete-a-block
values of h21*h22 are negative. Like you say, the results are too noisy to
be useful. We should still be able to report estimates of genetic
covariance though, and in particular test for non-zero genetic covariance,
instead of having the whole analysis crash, though.

On Tue, Dec 15, 2015 at 10:56 AM, Brendan Bulik-Sullivan <
notifications@github.com> wrote:

I think the issue is in the jackknife delete-a-block-values for
sqrt(h21*h22). Sometimes one of the delete-a-block jackknife values for one
of the traits will be negative even when both point estimates of h2 are
positive (typically this means that one of the h2 estimates is positive but
close to zero).

Here are the lines throwing the error

ldsc/ldscore/regressions.py

Lines 698 to 699 in 78874c7

	denom_delete_values = np.sqrt(
	np.multiply(hsq1.tot_delete_values, hsq2.tot_delete_values))

;
the software solution is to check that all entries of
hsq1.tot_delete_values and hsq2.tot_delete_values are positive. I haven't
thought about the best statistical solution. In principle, the posterior
distribution of a ratio will be very skew when the posterior of the
denominator has a lot of mass close to zero, so the mean and SE are poor
summary statistics. On the other hand, the SE will be really big, which
suffices to convey the key point, namely that the results are too noisy to
be useful.

@hilaryfinucane https://github.com/hilaryfinucane what do you think?

—
Reply to this email directly or view it on GitHub
#40 (comment).

[marionabustamante]

Sorry for my ignorance, but I cannot follow your comments. Do you mean that as the data is too noisy it is not possible to estimate the genetic correlation? No solution...?

In fact, the sample size of one of the datasets is quite small (N total=5758, N cases=2693) and when I estimate heritability I get a message saying this:
Total Observed scale h2: 0.0139 (0.0735)
Lambda GC: 1.0046
Mean Chi^2: 0.9982
Intercept: 0.9966 (0.0064)
Ratio: NA (mean chi^2 < 1).

And when I "munge" the data:
Metadata:
Mean chi^2 = 0.998
WARNING: mean chi^2 may be too small. !!!!
Lambda GC = 1.005
Max chi^2 = 36.187
12 Genome-wide significant SNPs (some may have been removed by filtering).

Mariona

[hilaryfinucane]

Hi Mariona,

Right, there is not enough signal in the data to estimate the genetic
correlation between these two traits. It will still be possible to test
whether there is non-zero genetic covariance, but that is not in the
current version of the code.

Best,

Hilary

On Tue, Dec 15, 2015 at 4:17 PM, marionabustamante <notifications@github.com

wrote:

Sorry for my ignorance, but I cannot follow your comments. Do you mean
that as the data is too noisy it is not possible to estimate the genetic
correlation? No solution...?

In fact, the sample size of one of the datasets is quite small (N
total=5758, N cases=2693) and when I estimate heritability I get a message
saying this:
Total Observed scale h2: 0.0139 (0.0735)
Lambda GC: 1.0046
Mean Chi^2: 0.9982
Intercept: 0.9966 (0.0064)
Ratio: NA (mean chi^2 < 1).

And when I "munge" the data:
Metadata:
Mean chi^2 = 0.998
WARNING: mean chi^2 may be too small. !!!!
Lambda GC = 1.005
Max chi^2 = 36.187
12 Genome-wide significant SNPs (some may have been removed by filtering).

Mariona

—
Reply to this email directly or view it on GitHub
#40 (comment).