[Genabel-commits] r1298 - in pkg/ProbABEL: doc src

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

Author: lckarssen
Date: 2013-08-19 13:04:24 +0200 (Mon, 19 Aug 2013)
New Revision: 1298

Modified:
   pkg/ProbABEL/doc/ChangeLog
   pkg/ProbABEL/doc/ProbABEL_manual.tex
   pkg/ProbABEL/src/main.cpp
Log:
ProbABEL: don't calculate chi^2 for the 2df model (when using --mmscore). We need a more complicated solution here (see equation just below Eq. (4) in the ProbABEL paper). 

Also updated the documentation to reflect this and the previous change.

The ErasmusMC holds copyright for this change.


Modified: pkg/ProbABEL/doc/ChangeLog
===================================================================
--- pkg/ProbABEL/doc/ChangeLog	2013-08-19 08:07:06 UTC (rev 1297)
+++ pkg/ProbABEL/doc/ChangeLog	2013-08-19 11:04:24 UTC (rev 1298)
@@ -1,6 +1,8 @@
 ***** v.0.4.0 (2013.08)
 * The output files now contain a chi^2 column with the chi^2 value based
-  on the LRT.
+  on the LRT when not using --mmscore. When using -mmscore, the chi^2
+  values are calculated from the Wald statistic (not implemented for the
+  2df model yet).
 * Fixed long-standing bug #1130: The Cox regression module of ProbABEL
   crashes (when using filevector files). Actually, the Cox PH module was
   broken ever since ProbABEL v0.1-9e. Thanks to the Erasmus Medical

Modified: pkg/ProbABEL/doc/ProbABEL_manual.tex
===================================================================
--- pkg/ProbABEL/doc/ProbABEL_manual.tex	2013-08-19 08:07:06 UTC (rev 1297)
+++ pkg/ProbABEL/doc/ProbABEL_manual.tex	2013-08-19 11:04:24 UTC (rev 1298)
@@ -33,6 +33,9 @@
 % get the links to the figures and tables right:
 \usepackage[all]{hypcap} % to be loaded after hyperref package
 
+% lowercase letters as footnote numerals (to avoid confusion with powers).
+\renewcommand{\thefootnote}{\alph{footnote}}
+
 \DeclareMathOperator{\var}{\mathbf{var}}
 
 \makeindex
@@ -488,10 +491,14 @@
 standard errors ($\text{SE}_\beta$), and the $\chi^2$ test value based
 on the likelihood ratio test. Note that for the additive, recessive,
 dominant and overdominant genetic models this is a $\chi^2$ of 1
-degree of freedom, whereas for the genotypic model this is a $\chi^2$
-of 2df.
+degree of freedom (df), whereas for the genotypic model this is a
+$\chi^2$ of 2df. If the \texttt{--mmscore} option is used, the
+$\chi^2$ values are calculated from the Wald statistic
+(1df)\footnote{For the genotypic (2df) model the $\chi^2$ values can't
+  be simply calculated using the Wald statistic, consequently the
+  $\chi^2$ values are set to \texttt{nan}. A fix for this needs to be
+  implemented still.}.
 
-
 \section{Preparing input files}
 After installing \PA{} you can find the \texttt{prepare\_data.R} file
 in the \texttt{scripts} directory. It is an R script that arranges

Modified: pkg/ProbABEL/src/main.cpp
===================================================================
--- pkg/ProbABEL/src/main.cpp	2013-08-19 08:07:06 UTC (rev 1297)
+++ pkg/ProbABEL/src/main.cpp	2013-08-19 11:04:24 UTC (rev 1298)
@@ -755,8 +755,20 @@
                 {
                     // 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;
+                    if (input_var.getNgpreds() == 2 && model == 0)
+                    {
+                        /* For the 2df model we can't simply use the
+                         * Wald statistic. This can be fixed using the
+                         * equation just below Eq.(4) in the ProbABEL
+                         * paper. TODO LCK
+                         */
+                        *chi2[model] << "nan";
+                    }
+                    else
+                    {
+                        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