[GenABEL-dev] faster polygenic

[Xia Shen]

Hi Yurii,

Bad news and good news:

[BAD]
I cannot open the two links that you sent...

[GOOD]
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.

The way to fit hglm using GenABEL results is straightforward:

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
2. create L that satisfies LL' = G by
   s <- svd(G)
   L <- s$u%*%diag(sqrt(s$d))
3. call hglm like
   hglm(y = phdata(...)$..., X = design matrix of fixed effects, Z = L, family = gaussian())
   where NA has to be removed before calling, gaussian() can be changed to binomial(link = 'logit'), etc.

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.

I think the residuals from the mixed model is the most interesting stuff that will be passed into the next GWAS in GenABEL.

Hopefully this can serve some help :)

Cheers

Xia


       o   \o/  _ o        _|    \ /    |_        o_   \o/   o
      /|\   |    /\    _\o   \o   |   o/   o/_    /\    |   /|\
      / \  / \  | \  /) |   ( \  /o\ / )   |  (\  / |  / \  / \
----------------------------------------------------------------------
 Xia Shen - PhD student in Bioinformatics, Uppsala University.
          - Computational Genetics Group: www.computationalgenetics.se
          - Personal URL: www.19850911.com

 All judgements are, in their rationale, STATISTICS.        - C.R.Rao
----------------------------------------------------------------------



On Feb 10, 2011, at 3:35 PM, Yurii Aulchenko wrote:

> Hi Gulya,
> 
> As I mentioned, I am also very bothered by the fact that 'polygenic'
> works so slowly. I cc to Xia who mentioned he had some idea on how to
> make it working faster -- Xia, would be great if you could contribute
> to the discussion! I also cc to GA devel list, so people know what is
> going on.
> 
> Here I will describe how polygenic works currently (very far from optimal!)
> 
> 'polygenic'
> 
> https://r-forge.r-project.org/scm/viewvc.php/pkg/GenABEL/R/polygenic.R?view=markup&revision=623&root=genabel
> 
> is numerically maximizing likelihood function computed by procedure 'polylik'
> 
> https://r-forge.r-project.org/scm/viewvc.php/pkg/GenABEL/R/polylik.R?view=markup&revision=545&root=genabel
> 
> For maximization, standard R procedures 'nlm' or 'optim' can be used
> -- see what happens in 'polygenic' code between lines 246-300.
> 
> Any ideas on how we can speed it up are very welcome!
> 
> Yurii

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20110216/0662c9e8/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smalllogo.jpg
Type: image/jpeg
Size: 16295 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20110216/0662c9e8/attachment-0001.jpg>