[Genabel-commits] r1741 - in branches/ProbABEL-v0.4.3-coxfix/ProbABEL: . src

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

Author: lckarssen
Date: 2014-05-26 19:28:26 +0200 (Mon, 26 May 2014)
New Revision: 1741

Modified:
   branches/ProbABEL-v0.4.3-coxfix/ProbABEL/configure.ac
   branches/ProbABEL-v0.4.3-coxfix/ProbABEL/src/coxph_data.cpp
Log:
First stab at solving bug #5726.
- Checks for infinite beta are now per beta component, instead of the norm of the beta vector.
- Explicit setting of the loglik, beta and se_beta to NAN is commented, so the actual values are printed.

This version is sent to a tester for feedback.


Modified: branches/ProbABEL-v0.4.3-coxfix/ProbABEL/configure.ac
===================================================================
--- branches/ProbABEL-v0.4.3-coxfix/ProbABEL/configure.ac	2014-05-26 14:55:52 UTC (rev 1740)
+++ branches/ProbABEL-v0.4.3-coxfix/ProbABEL/configure.ac	2014-05-26 17:28:26 UTC (rev 1741)
@@ -2,7 +2,7 @@
 # Process this file with autoconf to produce a configure script.
 
 AC_PREREQ([2.67])
-AC_INIT(ProbABEL, 0.4.3, genabel-devel at r-forge.wu-wien.ac.at)
+AC_INIT(ProbABEL, 0.4.3-coxfix1, genabel-devel at r-forge.wu-wien.ac.at)
 AM_INIT_AUTOMAKE([silent-rules subdir-objects])
 AM_SILENT_RULES([yes])
 AC_CONFIG_SRCDIR([src/data.h])

Modified: branches/ProbABEL-v0.4.3-coxfix/ProbABEL/src/coxph_data.cpp
===================================================================
--- branches/ProbABEL-v0.4.3-coxfix/ProbABEL/src/coxph_data.cpp	2014-05-26 14:55:52 UTC (rev 1740)
+++ branches/ProbABEL-v0.4.3-coxfix/ProbABEL/src/coxph_data.cpp	2014-05-26 17:28:26 UTC (rev 1741)
@@ -439,14 +439,19 @@
     niter = maxiter;
 
     // Check the results of the Cox fit; mirrored from the same checks
-    // in coxph_fit.S from the R survival package
+    // in coxph.fit.S and coxph.R from the R survival package.
 
     bool setToZero = false;
 
+    // Based on coxph.fit.S lines with 'which.sing' and coxph.R line
+    // with if(any(is.NA(coefficients))). These lines set coefficients
+    // to NA if flag < nvar (with nvar = ncol(x)) and maxiter >
+    // 0. coxph.R then checks for any NAs in the coefficients and
+    // outputs the warning message if NAs were found.
     if (flag < X.nrow && maxiter > 0) {
         cerr << "Warning for " << snpinfo.name[cursnp]
              << ": X matrix deemed to be singular,"
-             << " setting beta and se to 'nan'\n";
+             << " setting beta and se to 'NaN'\n";
         setToZero = true;
     }
 
@@ -461,21 +466,36 @@
     {
         cerr << "Warning for " << snpinfo.name[cursnp]
              << ": Cox regression ran out of iterations and did not converge,"
-             << " setting beta and se to 'nan'\n";
+             << " setting beta and se to 'NaN'\n";
         setToZero = true;
     } else {
 #if EIGEN
         VectorXd ueigen = u.data;
         MatrixXd imateigen = imat.data;
         VectorXd infs = ueigen.transpose() * imateigen;
+        infs = infs.cwiseAbs();
         VectorXd betaeigen = beta.data;
-        if ( infs.norm() > eps ||
-            infs.norm() > sqrt(eps) * betaeigen.norm() )
-        {
+        bool problems = false;
+
+        assert(betaeigen.size() == infs.size());
+
+        // We check the beta's for all coefficients
+        // (incl. covariates), maybe stick to only checking the SNP
+        // coefficient?
+        for (int i = 0; i < infs.size(); i++) {
+            if (infs[i] > eps && infs[i] > sqrt(eps) * abs(betaeigen[i])) {
+                problems = true;
+            }
+        }
+
+        if (problems) {
             cerr << "Warning for " << snpinfo.name[cursnp]
                  << ": beta may be infinite,"
-                 << " setting beta and se to 'nan'\n";
+                 << " setting beta and se to 'NaN'\n";
 
+            // cout << "beta values for SNP: " << snpinfo.name[cursnp]
+            //      << " are: " << betaeigen << std::endl;
+
             setToZero = true;
         }
 #else
@@ -491,9 +511,9 @@
         if (setToZero)
         {
             // Cox regression failed
-            sebeta[i] = NAN;
-            beta[i]   = NAN;
-            loglik    = NAN;
+            // sebeta[i] = NAN;
+            // beta[i]   = NAN;
+            // loglik    = NAN;
         } else {
             sebeta[i] = sqrt(imat.get(i, i));
             loglik = loglik_int[1];