[Genabel-commits] r1274 - pkg/ProbABEL/src

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

Author: lckarssen
Date: 2013-07-23 00:17:21 +0200 (Tue, 23 Jul 2013)
New Revision: 1274

Modified:
   pkg/ProbABEL/src/main.cpp
Log:
Replaced loglik with chi^2 based on LRT in the ProbABEL output. This implementation is not complete yet, in case of missing genotype data for a given SNP I simply print NaN because the null model likelihood should in that case be estimated using only the individuals that have genotype data for that SNP. 



Modified: pkg/ProbABEL/src/main.cpp
===================================================================
--- pkg/ProbABEL/src/main.cpp	2013-07-22 21:50:17 UTC (rev 1273)
+++ pkg/ProbABEL/src/main.cpp	2013-07-22 22:17:21 UTC (rev 1274)
@@ -259,11 +259,11 @@
             //Oct 26, 2009
         }
     }
-    *outfile[0] << input_var.getSep() << "loglik\n"; //"chi2_SNP_2df\n";
-    *outfile[1] << input_var.getSep() << "loglik\n"; //"chi2_SNP_A1\n";
-    *outfile[2] << input_var.getSep() << "loglik\n"; //"chi2_SNP_domA1\n";
-    *outfile[3] << input_var.getSep() << "loglik\n"; //"chi2_SNP_recA1\n";
-    *outfile[4] << input_var.getSep() << "loglik\n"; //"chi2_SNP_odom\n";
+    *outfile[0] << input_var.getSep() << "chi2_SNP_2df\n";  // "loglik\n";
+    *outfile[1] << input_var.getSep() << "chi2_SNP_A1\n";   // "loglik\n";
+    *outfile[2] << input_var.getSep() << "chi2_SNP_domA1\n";// "loglik\n";
+    *outfile[3] << input_var.getSep() << "chi2_SNP_recA1\n";// "loglik\n";
+    *outfile[4] << input_var.getSep() << "chi2_SNP_odom\n"; // "loglik\n";
 }
 
 void create_header2(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
@@ -296,7 +296,7 @@
                         << phd.model_terms[interaction_cox];
         }
 #endif
-        *outfile[0] << input_var.getSep() << "loglik"; //"chi2_SNP";
+        *outfile[0] << input_var.getSep() << "chi2_SNP"; //"loglik";
     }
     //Oct 26, 2009
     *outfile[0] << "\n";
@@ -443,6 +443,7 @@
     nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0,
                  input_var.getInteraction(), input_var.getNgpreds(), true, 1);
 #endif
+    double null_loglik = nrd.loglik;
 
     std::cout << " estimated null model ...";
     // end null
@@ -503,6 +504,7 @@
         covvalue[i]->precision(9);
         //Oct 26, 2009
         chi2.push_back(new std::ostringstream());
+        chi2[i]->precision(9);
     }
 
 
@@ -674,12 +676,23 @@
 
                     if (input_var.getScore() == 0)
                     {
-                        //*chi2[model] << 2.*(rd.loglik-null_loglik);
-                        *chi2[model] << rd.loglik;
+                        if (gcount != gtd.nids)
+                        {
+                            // If SNP data is missing we didn't
+                            // correctly compute the null likelihood
+                            *chi2[model] << "NaN";
+                        }
+                        else
+                        {
+                            // No missing SNP data, we can compute the LRT
+                            *chi2[model] << 2. * (rd.loglik - null_loglik);
+                        }
+                        //*chi2[model] << rd.loglik;
                     } else
                     {
-                        //*chi2[model] << rd.chi2_score;
-                        *chi2[model] << "nan";
+                        // We want score test output
+                        *chi2[model] << rd.chi2_score;
+                        //*chi2[model] << "nan";
                     }
                 }
             } // END first part of if(poly); allele not too rare