lukejoconnor / LCV

Software implementing the Latent Causal Variable Model
60 stars 17 forks source link

negative heritability estimates #13

Open humanpaingeneticslab opened 2 years ago

humanpaingeneticslab commented 2 years ago

Dear LCV users,

when using LCV as shown in the tutorial: LCV = RunLCV(data[,"L2"],data[,"Z.x"],data[,"Z.y"] ) I would often get the error "NaNs produced, probably due to negative heritability estimates."

However, when providing sample sizes, I got LCV to run without errors: LCV = RunLCV(data[,"L2"],data[,"Z.x"],data[,"Z.y"], n.1=n.1, n.2=n.2, ldsc.intercept=0 )

Notice that ldsc.intercept is set to 0, not 1 as stated in the header of RunLCV.R

humanpaingeneticslab commented 2 years ago

Setting ldsc.intercept=0 has other side effect... So I patched the code - not sure how Kosher this is though;

I now run it like this: RunLCV(data[,"L2"],data[,"Z.x"],data[,"Z.y"], n.1=n.1, n.2=n.2, ldsc.intercept=1 )

In the EstimateK4 function: In the code block for when 'ldsc.intercept=1', if the estimated heritability is negative then I resort to the way it is estimated when 'ldsc.intercept=0'; e.g.


  if(ldsc.intercept==0){ ... }
  else{
    keep.snps.1<-z.1^2 <= sig.threshold*mean(z.1^2)
    temp<WeightedRegression( ... );
    intercept.1<-temp[2];

    temp<-WeightedRegression( ... );
    h2g.1<-temp[1];

    # deal with negative heritability
    if( ( h2g.1 < 0 ) & ( n.1 > 1 ) ) {
      # like when ldsc.intercept==0
      intercept.1<-1/n.1
      temp<-WeightedRegression(ell,z.1^2-intercept.1,weights)
      h2g.1<-temp[1]
    }