|
|
Post by tfwillems on Jul 21, 2014 20:23:26 GMT
My understanding is that GREML is supposed to produce unbiased estimates of the variance explained by tagging SNPs.
Is this true for both the --reml and --rem-no-constrain options, or does the constraint introduce a potential bias?
On a similar note, do you recommend using --reml-no-constrain under certain scenarios, or is the default option (--reml) generally superior?
|
|
|
|
Post by Zhihong Zhu on Sept 6, 2014 2:55:16 GMT
Hi,
The "--reml" let the program know REML algorithm will be used in the analysis. And I would like to use "--rem-no-constrain" in the REML analysis. Usually, the results are the same, whether "--rem-no-constrain" is used or not. In some particular scenarios, the vg (genetic variance) may be negative if trying "--rem-no-constrain". You might constrain the vg to zero, because it should be >= 0 as the definition. And I would like to report the later one(ie. constrain to zero) , since it might be more easily accepted by people. In the scenarios, before reporting the results, I will double check if there is any difference except for vg.
Cheers,
Zhihong
|
|
ruomu
New Member
Posts: 3 
|
Post by ruomu on Apr 9, 2015 11:47:07 GMT
Hi,
I came across a problem. If I didn't use "rem-no-constrain", the V(e) was infinitely near zero, and V(G)/Vp_L (The value I want) was 0.53 maybe. If I used, the V(e) was negative, and V(G)/Vp_L value was 0.68. And the two values varied widely of course. I was wondering which one I should use, the former one?
If the former, then what is the purpose of "rem-no-contrain" command? Does it only for those situations when we even cannot get a value?
|
|
|
|
Post by Zhihong Zhu on Apr 11, 2015 2:18:42 GMT
Hi,
It is weird that the V(e) is nearly zero or negative. I guess you fit multiple grms or covariates. Please check the hg^2 ( V(G)/V(P) ) and V(e) if fitting only one grms, genotypes and phenotypes in the data.
Cheers,
Zhihong
|
|
ruomu
New Member
Posts: 3
|
Post by ruomu on Apr 13, 2015 2:52:12 GMT
Thank you for your reply.
But I had only one grm and not many covariates. Thank you for taking some time to see this: "./gcta64 --reml --grm commonsnps --pheno plink.fam --mpheno 4 --covar discrete.covar --qcovar quantity.covar --prevalence 0.007 --thread-num 20 --reml-maxit 1000 --reml-no-constrain --out commonsnps_s". This is the code I used. The file of "--grm commonsnps" came from "--make-grm" of all the common snps. "discrete.covar" only includes gender, and "quantity.covar" includes 20 pcas. This is a case-control data, and there is only one trait: disease or not.
Thank you.
|
|
|
|
Post by Zhihong Zhu on Apr 13, 2015 23:00:56 GMT
Hi,
The command looks fine. But I'm not sure whether there is enough power to detect the hg, ie, whether the sample size is large enough or not. Please check the SE and power by GCTA power calculator (http://spark.rstudio.com/ctgg/gctaPower/). If the SE is very large, or even almost equal to the vg, it is hard to say vg!= 0.
Cheers,
Zhihong
|
|
ruomu
New Member
Posts: 3
|
Post by ruomu on Apr 17, 2015 2:54:54 GMT
Thank you very much for your advice.
I tried the GCTA power calculator and got a SE of 0.0086 and power of 1 (in case the "Variance of the SNP-derived genetic relationships" means V11 from the matrix of --reml output if Vg was the first component, I am not sure of this). My sample size was 944: 1195 and Vg was 0.23. Does this mean a proper SE ?
Thank you.
|
|
|
|
Post by Zhihong Zhu on Apr 25, 2015 5:27:49 GMT
Hi,
I was wondering the estimate of VG and SE(VG), because the negative VE suggests the VG could be >VP, suppose only one GRM fit in the model. Please compare the SE(VG) from GCTA power calculator with that from your estimates. And please double check the SE(VG) got from GCTA power calculator, since 0.0086 seems to be a little small, given the sample size of 944 : 1195 and the risk prevalence 0.007.
PS: There must be V(1)/VP and V(1)_L/VP_L in the .hsq file if you specify the risk prevalence. V(1)/VP - proportion of variance explained by the first GRM in observed scale. Because in your case the GRM must be estimated from SNPs, V(1)/VP is the SNP heritability in the observed scale. V(1)_L / VP_L - the SNP heritability in liability scale.
Cheers,
Zhihong
|
|
