Post by regevsch on Jan 18, 2020 20:42:40 GMT
Hey,
I am using GCTA's Haseman-Elston regression to estimate the heritability of my trait. This is my output:
Trait #1 is included for analysis.
Non-missing phenotypes of 8974 individuals are included from [...].
8974 individuals are in common in these files.
Performing Haseman-Elston regression ...
Standardising the phenotype ...
Constructing ordinary least squares equations ...
Left-hand side of OLS equations (X'X)
4.02619e+07 -3985.72
-3985.72 1806.33
HE-CP
Coefficient Estimate SE_OLS SE_Jackknife P_OLS P_Jackknife
Intercept -4.9121e-05 0.000157576 6.60569e-06 0.755247 1.03647e-13
V(G)/Vp 0.629447 0.0235255 0.0702348 1.05823e-157 3.18756e-19
HE-SD
Coefficient Estimate SE_OLS SE_Jackknife P_OLS P_Jackknife
Intercept -0.999938 0.000232525 0.016196 0 0
V(G)/Vp 0.621837 0.034715 0.0470234 9.40868e-72 6.38067e-40
I should be happy, since I got h^2=0.62 with narrow CIs, but I was suspicious.
I tried estimating with GCTA's REML but it did not converge (V(G) exploded). I used an simpler, alternative REML program which gave me the estimate of h^2 = 0.006.
I tried recovering the HE result myself using the linear regression underlying the model (basically, scale Y to have zero mean and unit variance; linear-regress K and I on Y*Y; get the coefficients of this linear regression as sigma_g^2 and sigma_e^2). I got 0.0072.
At this point I distrust GCTA's HE output. I wonder though how to explain the discrepancy. Any clues? Is the source code available for this function?
Thanks!
I am using GCTA's Haseman-Elston regression to estimate the heritability of my trait. This is my output:
Trait #1 is included for analysis.
Non-missing phenotypes of 8974 individuals are included from [...].
8974 individuals are in common in these files.
Performing Haseman-Elston regression ...
Standardising the phenotype ...
Constructing ordinary least squares equations ...
Left-hand side of OLS equations (X'X)
4.02619e+07 -3985.72
-3985.72 1806.33
HE-CP
Coefficient Estimate SE_OLS SE_Jackknife P_OLS P_Jackknife
Intercept -4.9121e-05 0.000157576 6.60569e-06 0.755247 1.03647e-13
V(G)/Vp 0.629447 0.0235255 0.0702348 1.05823e-157 3.18756e-19
HE-SD
Coefficient Estimate SE_OLS SE_Jackknife P_OLS P_Jackknife
Intercept -0.999938 0.000232525 0.016196 0 0
V(G)/Vp 0.621837 0.034715 0.0470234 9.40868e-72 6.38067e-40
I should be happy, since I got h^2=0.62 with narrow CIs, but I was suspicious.
I tried estimating with GCTA's REML but it did not converge (V(G) exploded). I used an simpler, alternative REML program which gave me the estimate of h^2 = 0.006.
I tried recovering the HE result myself using the linear regression underlying the model (basically, scale Y to have zero mean and unit variance; linear-regress K and I on Y*Y; get the coefficients of this linear regression as sigma_g^2 and sigma_e^2). I got 0.0072.
At this point I distrust GCTA's HE output. I wonder though how to explain the discrepancy. Any clues? Is the source code available for this function?
Thanks!