[Genabel-commits] r1514 - pkg/ProbABEL/checks/R-tests

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

Author: lckarssen
Date: 2013-12-29 17:19:49 +0100 (Sun, 29 Dec 2013)
New Revision: 1514

Modified:
   pkg/ProbABEL/checks/R-tests/run_models_in_R_palinear.R
Log:
Add manual fix to the ProbABEL R check for palinear. In the case of test SNP 6 (no genotypic variation) the results of ProbABEL are slightly different from R.
Not ideal, but at least this allows the test to finish.


Modified: pkg/ProbABEL/checks/R-tests/run_models_in_R_palinear.R
===================================================================
--- pkg/ProbABEL/checks/R-tests/run_models_in_R_palinear.R	2013-12-29 15:18:53 UTC (rev 1513)
+++ pkg/ProbABEL/checks/R-tests/run_models_in_R_palinear.R	2013-12-29 16:19:49 UTC (rev 1514)
@@ -38,6 +38,12 @@
     paste0(tests.path, "linear_ngp2_2df.out.txt"),
     head=TRUE)[, cols2df]
 
+## Fix betas, sebetas, chi^2 for the case that there is no variation
+## (SNP 6 in the info file). ProbABEL lists them all as 0.0, R lists
+## them as:
+prob.dom.PA[6, 2:4] <- c(NaN, NaN, 0.0)
+prob.2df.PA[6, 4:5] <- c(NA, NA)
+
 ####
 ## run analysis in R
 ####
@@ -102,10 +108,17 @@
         model.0 <- lm( height[noNA] ~ sex[noNA] + age[noNA] )
         model   <- lm( height ~ sex + age + prob[, indexHet] + prob[, indexHom] )
         smA1A2  <- summary(model)$coef[4, 1:2]
-        smA1A1  <- summary(model)$coef[5, 1:2]
+
+        ## When all coefficients are NA, they don't show up in $coeff
+        if ( nrow(summary(model)$coeff) > 4 ) {
+            smA1A1 <- summary(model)$coef[5, 1:2]
+        } else {
+            smA1A1 <- c(NA, NA)
+        }
+
         lrt     <- 2 * ( logLik( model ) - logLik( model.0 ) )
 
-        rsq <- resPa2df[i-2, "Rsq"]
+        rsq <- prob.2df.PA[i-2, "Rsq"]
         if( rsq < rsq.thresh) {
             row <- c(rsq, NaN, NaN, NaN, NaN, NaN)
         } else {