in fact you should not use the residuals but the phenofile shoud contain the original phenotype and the covariates also.
However, it is still not clear if the invmat correlation matrix should be generated from a run of polygenic_hglm where the covariates were included or when they were excluded.
So is this the proper approach?:
1) Use polygenic_hgm with trait.type = "binomial" to compute the required correlation matrix, while adjusting for the covariates:
Code: Select all
polygene <- polygenic_hglm(obesity ~ age + as.factor(sex) + I(age^2) + age:as.factor(sex), ibs2, data, trait.type = "binomial",verbose=TRUE)
invmat <- polygene$InvSigma * polygene$h2an$estimate[length(polygene$h2an$estimate)]
2) Then use palogist --mmscore to regress the binary trait against exactly the same set of covariates plus the SNP of interest.
The InvSigma matrix is computed by the polygenic_hglm program as follows:
Code: Select all
Sigma <- tVar * out$esth2 * relmat + diag(tVar * (1 - out$esth2), ncol = length(y), nrow = length(y))
out$InvSigma <- ginv(Sigma)
where
Code: Select all
tVar <- res_hglm$varRan + res_hglm$varFix
out$esth2 <- res_hglm$varRan/tVar
so it is scaled by the esth2 estimate of heritability, but if we are estimating h2 while adjusting for exactly the same set of covariates,
then I think this might be O.K.
Thank you,
Dan Weeks