[GenABEL-dev] Interesting idea by Gulnara Svischeva

[Yurii Aulchenko]

Ok, here are test results. I did comparison with h2-estimates obtained
by native 'polygenic', 'polygenic_hglm' and 'REML.rotate'.

Your (modified) procedure, test-script, and resulting plot attached.
In summary: very good agreement of h2-estimates between
'polygenic_hglm' and 'REML.rotate'; quite bad timing for REML.rotate.
I do not this this should discourage us though: polygenic-ML, based on
the same ideas, shows excellent timing.

-------- More details ---------------

Changes compared to your code:

1) In REML.rotate: as I ran tests on 'genomic kinship', some
eigenvalues are negative. Therefore I changed

N <- eigenvalues
to
N <- abs(eigenvalues)

2) In REML.rotate: added max.h2 to the output:

list(fitted.wls=wls, computed.h2=allh2, REMLlogL=allREML.fit)
to
list(esth2=max.h2,fitted.wls=wls, computed.h2=allh2, REMLlogL=allREML.fit)

3) During tests: used conv=1e-8 for both REML.rotate and
polygenic_hglm for better convergence.


Results:

Very good agreement between estimates obtained by polygenic_hglm and
REML.rotate (as it should be, I presume :) ). See attached pdf plot
(on X: h2 estimate from native 'polygenic'; Y: either polygenic_hglm
[black circles] or REML.rotate [red ones] ). Lack of complete
agreement native - REML is known issue we are at doubt to explain, see

http://lists.r-forge.r-project.org/pipermail/genabel-devel/2011-March/000132.html


Timing:

With conv=1e-8 mean time for estimation on 300 IDs was:

17.81 sec for REML.rotate
2.27 sec for 'polygenic'-ML
1.33 sec for polygenic_hglm

Using default conv=1e-6 will make ~5 sec for REML.rotate, so, still
rather slow.

Taking this all together I would say REML.rotate is suboptimal --
'polygenic'-ML uses the same ideas, does exhaustive numeric
maximization, and is very fast!

2011/5/18 Lars Rönnegård <lrn at du.se>:
> Hi,
>
> I attach an implementation of an idea proposed by Gulnara. It uses weighted
> least squares to compute the REML heritability estimate.
>
> The function is found in “REMLrotateFun.R” and some example code is found in
> “ExampleCode_REMLrotate.R” that uses the data in “RemlEx.RData”.
>
>
>
> If someone would be interested in trying it (and possibly also suggest
> improvements), I would appreciate it.
>
> It works for linear mixed models having a polygenic random effect and iid
> residuals. It seems to be very fast. The example I attach contains 680
> observations.
>
> Is this a new algorithm or has it been done already?
>
>
>
> Best wishes,
>
> Lars Rönnegård
>
> Dalarna University
>
> _______________________________________________
> genabel-devel mailing list
> genabel-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-devel
>
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: REMLrotateFun.R
Type: application/octet-stream
Size: 2409 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20110518/77ccdadc/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test_poly_hglm_RRot.pdf
Type: application/pdf
Size: 25633 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20110518/77ccdadc/attachment-0001.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test_REMLrotate.R
Type: application/octet-stream
Size: 1558 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20110518/77ccdadc/attachment-0003.obj>