improve SKAT and documentation

ChangeLog

@@ -1,6 +1,8 @@
 2015-05-09  Xiaowei Zhan  <zhanxw@fantasia.csgstat.sph.umich.edu>
 
 	* Fix and verified '--meta cov' outputs
+	* Improve SKAT permutation
+	* Improve documentation
 
 2015-04-30  Xiaowei Zhan  <zhanxw@fantasia.csgstat.sph.umich.edu>
 

regression/Makefile

@@ -23,6 +23,7 @@ BASE = LogisticRegression \
        MultivariateNormalDistribution \
        MultivariateVT \
        FamSkat \
+       cdflib \
 
 OBJ = $(BASE:%=%.o)
 

regression/MixtureChiSquare.cpp

@@ -1,6 +1,7 @@
 #include "MixtureChiSquare.h"
 
 #include <stdlib.h>
+#include "cdflib.h"
 #include "qfc.c"
 
 double MixtureChiSquare::getPvalue(double Q) {
@@ -21,6 +22,60 @@ double MixtureChiSquare::getPvalue(double Q) {
   return pValue;
 }
 
+double sum(double* d, int n, int power) {
+  double r = 0.0;
+  double tmp;
+  for (int i = 0; i < n; ++i) {
+    tmp = d[i];
+    for (int j = 1; j < power; ++j) {
+      tmp *= d[i];
+    }
+    r += tmp;
+  }
+  return r;
+}
+
+double MixtureChiSquare::getLiuPvalue(double Q) {
+  const double c1 = sum(lambda, lambda_size, 1);
+  const double c2 = sum(lambda, lambda_size, 2);
+  const double c3 = sum(lambda, lambda_size, 3);
+  const double c4 = sum(lambda, lambda_size, 4);
+  double s1 = c3 / c2 / sqrt(c2);
+  double s2 = c4 / c2 / c2;
+  const double muQ = c1;
+  const double sigmaQ = sqrt(2.0 * c2);
+  const double tstar = (Q - muQ) / sigmaQ;
+
+  double a;
+  double delta;
+  double l;
+  if (s1 * s1 > s2) {
+    a = 1 / (s1 - sqrt(s1 * s1 - s2));
+    delta = (s1 * a - 1) * a * a;
+    l = a * a - 2.0 * delta;
+  } else {
+    a = 1.0 / s1;
+    delta = 0.0;
+    l = c2 * c2 * c2 / c3 / c3;
+  }
+  const double muX = l + delta;
+  const double sigmaX = sqrt(2) * a;
+
+  int which = 1;
+  double p;
+  double q;
+  double x = tstar * sigmaX + muX;
+  double ncp = delta;
+  int status;
+  double bound;
+
+  cdfchn(&which, &p, &q, &x, &l, &ncp, &status, &bound);
+  if (status != 0) {
+    return 1;
+  }
+  return q;
+}
+
 void MixtureChiSquare::dumpLambda() const {
   for (int i = 0; i < lambda_size; ++i) {
     fprintf(stderr, "lambda[%d] = %g\n", i, lambda[i]);

regression/MixtureChiSquare.h

@@ -52,12 +52,15 @@ class MixtureChiSquare {
     df = newDf;
     lambda_cap = newCap;
   };
+  // use Davies's method to calculate P-value
   double getPvalue(double Q);
+  // use Liu's method to calculate P-value
+  double getLiuPvalue(double Q);  
   void dumpLambda() const;
 
  private:
   const double sigma;  // coefficient of standard normal variable
-  const int lim;       // maximum number of terms in tegration
+  const int lim;       // maximum number of terms integration
   const double acc;    // accuracy
   int lambda_cap;      // capacity of *lambda
 

regression/Skat.cpp

@@ -3,7 +3,6 @@
 #include "EigenMatrixInterface.h"
 #include <Eigen/Core>
 #include <Eigen/Eigenvalues>
-//#include <Eigen/Dense>
 #include "MixtureChiSquare.h"
 
 // #define DEBUG
@@ -88,6 +87,9 @@ class Skat::SkatImpl {
     }
     // calculate p-value
     this->pValue = this->mixChiSq.getPvalue(this->Q);
+    if (this->pValue == 0.0 || this->pValue == 1.0) {
+      this->pValue = this->mixChiSq.getLiuPvalue(this->Q);      
+    }
     return 0;
   };
 

regression/test/Makefile

@@ -9,7 +9,8 @@ EXE = testGSLIntegration \
       testLogisticVT \
       testLinearVT \
       testMultivariateNormalDistribution \
-      testLMM
+      testLMM \
+      testMixtureChiSquare
 
 all: $(EXE)
 debug: $(EXE)

regression/test/testLMM.R

@@ -45,8 +45,9 @@ getEst <- function(param) {
 
   f <- function(param, mask) {
     print(param)
-    print(mask)
+   print(mask)
     Sigma <- param[1] * k1 * mask[1] + param[2] * k2 * mask[2] + param[3] * I * mask[3]
+    # Sigma <- param[1] * k1  + param[2] * k2  + param[3] * I * 0
     SigmaInv <- solve(Sigma)
     beta <- solve(t(Cov) %*% SigmaInv %*% Cov, t(Cov) %*% SigmaInv %*% y)
     resid <- y - Cov %*% beta
@@ -61,6 +62,7 @@ getEst <- function(param) {
 
   #optim(param, f, method = "Nelder-Mead", lower = rep(1e-8, 3))
   optim(param, f, mask = c(1, 1, 0), method = "Nelder-Mead")
+  optim(param, f, mask = c(1, 0, 1), method = "Nelder-Mead")
 
   getUV <- function(param) {
     print(param)

src/Main.cpp

@@ -28,7 +28,7 @@
 
 Logger* logger = NULL;
 
-const char* VERSION = "20150430";
+const char* VERSION = "20150509";
 
 void banner(FILE* fp) {
   const char* string =
@@ -593,9 +593,9 @@ int main(int argc, char** argv) {
   ADD_STRING_PARAMETER(pl, cov, "--covar", "Specify covariate file")
   ADD_STRING_PARAMETER(
       pl, covName, "--covar-name",
-      "Specify the column name in coavriate file to be included in analysis")
+      "Specify the column name in covariate file to be included in analysis")
   ADD_BOOL_PARAMETER(pl, sex, "--sex",
-                     "Include sex (5th) as covaraite from PED file")
+                     "Include sex (5th) as covariates from PED file")
 
   ADD_PARAMETER_GROUP(pl, "Specify Phenotype")
   ADD_STRING_PARAMETER(pl, pheno, "--pheno", "Specify phenotype file")
@@ -616,7 +616,7 @@ int main(int argc, char** argv) {
   // ADD_STRING_PARAMETER(pl, glTag, "--gl", "Specify which genotype likelihood
   // tag to use. (e.g. GL)")
 
-  ADD_PARAMETER_GROUP(pl, "Chromsome X Options")
+  ADD_PARAMETER_GROUP(pl, "Chromosome X Options")
   ADD_STRING_PARAMETER(pl, xLabel, "--xLabel",
                        "Specify X chromosome label (default: 23|X")
   ADD_STRING_PARAMETER(pl, xParRegion, "--xParRegion",
@@ -658,7 +658,7 @@ int main(int argc, char** argv) {
                     "Specify minimum Minor Allele Count(inclusive) to be "
                     "included in analysis")
   ADD_STRING_PARAMETER(pl, annoType, "--annoType",
-                       "Specify annotation type that is follwed by ANNO= in "
+                       "Specify annotation type that is followed by ANNO= in "
                        "the VCF INFO field, regular expression is allowed ")
 
   ADD_PARAMETER_GROUP(pl, "Genotype Filter")
@@ -735,7 +735,7 @@ int main(int argc, char** argv) {
   ADD_STRING_PARAMETER(pl, condition, "--condition",
                        "Specify markers to be conditions (specify range)")
 
-  ADD_PARAMETER_GROUP(pl, "Auxilliary Functions")
+  ADD_PARAMETER_GROUP(pl, "Auxiliary Functions")
   ADD_BOOL_PARAMETER(pl, help, "--help", "Print detailed help message")
   END_PARAMETER_LIST(pl);
 

src/Permutation.h

@@ -78,7 +78,7 @@ class Permutation {
    */
   bool next() {
     if (this->actualPerm >= this->numPerm) return false;
-    if (numX + numEqual > threshold) {
+    if (numX + numEqual >= threshold) {
       return false;
     }
     return true;