[Genabel-commits] r1297 - in pkg/ProbABEL: examples src

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

Author: lckarssen
Date: 2013-08-19 10:07:06 +0200 (Mon, 19 Aug 2013)
New Revision: 1297

Modified:
   pkg/ProbABEL/examples/example_mms.sh
   pkg/ProbABEL/src/main.cpp
Log:
More ProbABEL chi^2. In case of mmscore we now output chi^2 as well. However, since mmscore is a REML method we can't use the LRT. Consequently, the chi^2 in this case is calculated from the Wald test.

src/main.cpp:
 - Implementation of chi^2 in output in case of mmscore.  
examples/example_mms.sh: 
 - Uncomment checks of prob vs. dose (with and without filevector input)


Copyright for this change lies with the Erasmus MC.


Modified: pkg/ProbABEL/examples/example_mms.sh
===================================================================
--- pkg/ProbABEL/examples/example_mms.sh	2013-08-19 07:08:54 UTC (rev 1296)
+++ pkg/ProbABEL/examples/example_mms.sh	2013-08-19 08:07:06 UTC (rev 1297)
@@ -64,14 +64,12 @@
         "mmscore check ($model model): prob vs. prob_fv"
 done
 
-# The following checks are disabled because of the missing LogLik
-# column in the prob data
-# run_diff mmscore_prob_add.out.txt \
-#     mmscore_add.out.txt \
-#     "mmscore check: prob vs. dose" \
-#     -I beta_SNP
+run_diff mmscore_prob_add.out.txt \
+    mmscore_dose_add.out.txt \
+    "mmscore check: prob vs. dose" \
+    -I SNP
 
-# run_diff mmscore_prob_fv_add.out.txt \
-#     mmscore_fv_add.out.txt \
-#     "mmscore check: prob_fv vs. dose_fv" \
-#     -I beta_SNP
+run_diff mmscore_prob_fv_add.out.txt \
+    mmscore_dose_fv_add.out.txt \
+    "mmscore check: prob_fv vs. dose_fv" \
+    -I SNP

Modified: pkg/ProbABEL/src/main.cpp
===================================================================
--- pkg/ProbABEL/src/main.cpp	2013-08-19 07:08:54 UTC (rev 1296)
+++ pkg/ProbABEL/src/main.cpp	2013-08-19 08:07:06 UTC (rev 1297)
@@ -216,9 +216,8 @@
                             << phd.model_terms[interaction_cox];
             }
 #endif
-            *outfile[0] << input_var.getSep() << "chi2_SNP"; //"loglik";
         }
-        //Oct 26, 2009
+        *outfile[0] << input_var.getSep() << "chi2_SNP";
         *outfile[0] << "\n";
     } // ngpreds == 1
     else if (input_var.getNgpreds() == 2) // prob data: all models
@@ -692,11 +691,8 @@
                 //________________________________
                 //cout <<  rd.loglik<<" "<<input_var.getNgpreds() << "\n";
 
-                if ((input_var.getInverseFilename() == NULL
-                        && input_var.getNgpreds() != 2)
-                        || input_var.getNgpreds() == 2)
-                {
-
+                if (input_var.getInverseFilename() == NULL)
+                { // Only if we don't have an inv.sigma file can we use LRT
                     if (input_var.getScore() == 0)
                     {
                         double loglik = rd.loglik;
@@ -754,6 +750,13 @@
                         // We want score test output
                         *chi2[model] << rd.chi2_score;
                     }
+                } // END if( inv.sigma == NULL )
+                else if (input_var.getInverseFilename() != NULL)
+                {
+                    // We can't use the LRT here, because mmscore is a
+                    // REML method. Therefore go for the Wald test
+                    double Z = rd.beta[start_pos] / rd.sebeta[start_pos];
+                    *chi2[model] << Z * Z;
                 }
             } // END first part of if(poly); allele not too rare
             else
@@ -882,22 +885,14 @@
         }
         else // Dose data: only additive model. Only one output file
         {
-            if (input_var.getInverseFilename() == NULL)
-            {
-                //Han Chen
-                *outfile[0] << beta_sebeta[0]->str() << input_var.getSep();
+            *outfile[0] << beta_sebeta[0]->str() << input_var.getSep();
 #if !COXPH
-                if (!input_var.getAllcov() && input_var.getInteraction() != 0)
-                {
-                    *outfile[0] << covvalue[0]->str() << input_var.getSep();
-                }
-#endif
-                *outfile[0] << chi2[0]->str() << "\n";
-                //Oct 26, 2009
-            } else
+            if (!input_var.getAllcov() && input_var.getInteraction() != 0)
             {
-                *outfile[0] << beta_sebeta[0]->str() << "\n";
+                *outfile[0] << covvalue[0]->str() << input_var.getSep();
             }
+#endif
+            *outfile[0] << chi2[0]->str() << "\n";
         }  // End ngpreds == 1 when writing output files