[Genabel-commits] r1891 - pkg/OmicABEL/src

[noreply at r-forge.r-project.org]

Author: dfabregat
Date: 2014-11-18 13:13:25 +0100 (Tue, 18 Nov 2014)
New Revision: 1891

Modified:
   pkg/OmicABEL/src/REML.c
Log:
Fixed a bug in the estimation of heritability and residual
variance (h, sigma).


Modified: pkg/OmicABEL/src/REML.c
===================================================================
--- pkg/OmicABEL/src/REML.c	2014-11-10 20:53:26 UTC (rev 1890)
+++ pkg/OmicABEL/src/REML.c	2014-11-18 12:13:25 UTC (rev 1891)
@@ -189,8 +189,11 @@
             &MINUS_ONE, X, &n, beta, &iONE,
             &ONE, YmXB, &iONE );
 
+    // This was wrong
     // residual sigma and loglik
-	*res_sigma = variance( YmXB, n );
+    /**res_sigma = variance( YmXB, n );*/
+    //
+
 	// loglik = a + b
 	//  a -> log(det(M))
 	//  b -> YmXB' inv(M) YmXB
@@ -207,6 +210,10 @@
 	*loglik = 0.0;
 	for (i = 0; i < n; i++ )
 		*loglik += ZtY_upd[i] * D[i] * D[i] * ZtY_upd[i];
+
+    // res_sigma = 1/n * ( YmXB' inv(M) YmXB )
+    *res_sigma = (*loglik) / n;
+
 	*loglik = (*loglik) / (*res_sigma); // <- b
 	// a + b -> det + loglik above
 	//   det incomplete, we must take into account the sigma scaling M -> "+ n * log(*res_sigma)"