Hi Xia,<br><br>just gave your procedure a try -- and time-wise it works _incredibly_ fast (code attached FYI)! So, this sounds like a good way to go to speed up 'polygenic' -- I thinking of writing a wrapper (say, 'polygenic_hglm') which will require(hglm), prepare the data, call hglm and the format the output to make a valid 'polygenic' object.<br>
<br>Few things I could not see answer right away -- can you please help? -- it will save me quite some time!<br><br>1) do I understand correctly that fixed effects estimates are in hglm(...)$fixef<br>2) do I understand correctly that 'breeding values' are in hglm(...)$ranef<br>
3) do I understand correctly that environmental residuals are in hglm(...)$resid<br>4) how can I get hold on the estimates of polygenic and environmental (residual) variances, and find out what was heritability estimate<br>
5) is there an easy way to get hold on inverse of variance-covariance matrix. I mean given parameters estimates, it is no problem to compute that, but just in case it is already hanging around would be stupid to do so<br>
<br>best, and many thanks for you efforts,<br>Yurii<br><br>[code]<br>library(GenABEL)<br>library(hglm)<br>data(ge03d2.clean)<br>df <- ge03d2.clean[1:200,autosomal(ge03d2.clean)]<br>gkin <- ibs(df,w="freq")<br>
<br>formula <- as.formula(height ~ sex + age)<br>t0 <- system.time(res0 <- polygenic(formula,kin=gkin,df))<br><br>grel <- gkin<br>grel[upper.tri(grel)] <- t(grel)[upper.tri(grel)]<br>data <- phdata(df)<br>
mf <- model.frame(formula,data,na.action=na.omit,drop.unused.levels=TRUE)<br>y <- model.response(mf)<br>desmat <- model.matrix(formula,mf)<br>phids <- rownames(data)[rownames(data) %in% rownames(mf)]<br>relmat <- grel[phids,phids]*2.0<br>
s <- svd(grel)<br>L <- s$u %*% diag(sqrt(s$d))<br>t1 <- system.time(res1 <- hglm(y = y, X = desmat, Z = L, family = gaussian()))<br><br>[/code]<br><br><div class="gmail_quote">On Wed, Feb 16, 2011 at 2:11 PM, Xia Shen <span dir="ltr"><<a href="mailto:xia.shen@lcb.uu.se">xia.shen@lcb.uu.se</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><div style="word-wrap: break-word;">Hi Yurii,<div><br></div><div>Bad news and good news:</div>
<div><br></div><div>[BAD]</div><div>I cannot open the two links that you sent...</div><div><br></div><div>[GOOD]</div><div>One of my colleague happens to be playing with your 'polygenic' function, and has got stuck by some binary trait because the numerical algorithm didn't converge. However, treating the trait as normal works. To fit the same mixed model, I tried my 'hglm' package (published on the recent issue of The R Journal), it worked pretty well, and fitting binary, poisson, gamma, ... traits will all work.</div>
<div><br></div><div>The way to fit hglm using GenABEL results is straightforward:</div><div><br></div><div>1. create the IBS matrix G, has to be symmetric, just need to replicate the lower triangle of the IBS from GenABEL and overwrite the upper triangle</div>
<div>2. create L that satisfies LL' = G by</div><div> s <- svd(G)</div><div> L <- s$u%*%diag(sqrt(s$d))</div><div>3. call hglm like</div><div> hglm(y = phdata(...)$..., X = design matrix of fixed effects, Z = L, family = gaussian())</div>
<div> where NA has to be removed before calling, gaussian() can be changed to binomial(link = 'logit'), etc.</div><div><br>
The hglm function will return all the parameter estimates (fixed effects, variance components) and random effects predicts and residuals, etc. The algorithm is based on inter-connected GLMs (h-likelihood theory, the work by Lee & Nelder in the last decade), efficient and stable, without establishing the likelihood and calling numeric algorithm like FGLS.</div>
<div><br></div><div>I think the residuals from the mixed model is the most interesting stuff that will be passed into the next GWAS in GenABEL.</div><div><br></div><div>Hopefully this can serve some help :)</div><div><br>
</div><div>Cheers</div><div><br></div><div>Xia</div></div><br></blockquote></div>