Commit 2864a7fe authored by zhanxw's avatar zhanxw
Browse files

fix merge conflicts, remove compiler flag -Wno-error=int-in-bool-context

parents b629ab4d cdd07dc6
......@@ -3,6 +3,10 @@
* Improve speed of covariance calculations
* Clean codes
2017-07-13 Xiaowei Zhan <zhanxw@fantasia.csgstat.sph.umich.edu>
* Fix beta estimation in binary score tests for unrelated individuals
2017-06-13 Xiaowei Zhan <zhanxw@fantasia.csgstat.sph.umich.edu>
* Fix a bug related Wald test outputs (fix #30)
......
......@@ -55,7 +55,7 @@ GONCALO_LIB_DBG = $(ROOT)/libsrc/lib-dbg-goncalo.a
##################################################
## Common compiler options
CXX ?= g++
DEFAULT_CXXFLAGS = -D__STDC_LIMIT_MACROS -std=c++0x -Wall -Wno-unused-function $(OPENMP_FLAG) -DGIT_VERSION="\"$(GIT_VERSION)\"" -Wno-error=int-in-bool-context
DEFAULT_CXXFLAGS = -D__STDC_LIMIT_MACROS -std=c++0x -Wall -Wno-unused-function $(OPENMP_FLAG) -DGIT_VERSION="\"$(GIT_VERSION)\""
INCLUDE = $(THIRD_INC) $(REGRESSION_INC) $(VCF_INC) $(BASE_INC) $(GONCALO_INC)
LIB = $(REGRESSION_LIB) $(VCF_LIB) $(BASE_LIB) $(GONCALO_LIB) $(THIRD_LIB)
......
......@@ -195,6 +195,14 @@ class VCFRecord {
return;
}
for (unsigned int i = 9; i < sa.size(); i++) {
// do not allow empty name
if (sa[i].empty()) {
fprintf(stderr,
"One inddividual (column %d, or extra tab at the line end) has "
"an empty column header, please check file format\n",
(int)i);
exit(1);
}
int idx = i - 9;
VCFIndividual* p = new VCFIndividual;
this->allIndv[idx] = p;
......
......@@ -1345,6 +1345,20 @@ class BoltLMM::BoltLMMImpl {
}
fprintf(stderr, "ratio = %f\n", x_V_inv_x.sum() / x_x.sum());
}
fprintf(stderr, "infStatCalibration_ = %f\n", infStatCalibration_);
// calculate empirically X'HX / X'X
xVx_xx_ratio_ = x_V_inv_x.sum() / x_x.sum();
if (!std::isfinite(xVx_xx_ratio_)) {
xVx_xx_ratio_ = 1.0;
}
fprintf(stderr, "\ni\tx_v_inv_x\tx_x\tratio\n");
for (int i = 0; i < nSnp; ++i) {
fprintf(stderr, "%d\t%f\t%f\t%f\n", i, x_V_inv_x(i), x_x(i),
x_x(i) == 0 ? 0.0 : x_V_inv_x(i) / x_x(i));
}
fprintf(stderr, "ratio = %f\n", x_V_inv_x.sum() / x_x.sum());
return 0;
}
......
......@@ -30,7 +30,7 @@
Logger* logger = NULL;
const char* VERSION = "20170731";
const char* VERSION = "20170818";
void banner(FILE* fp) {
const char* string =
......@@ -41,7 +41,7 @@ void banner(FILE* fp) {
"| Bingshan Li, Dajiang Liu | \n"
"| Goncalo Abecasis | \n"
"| zhanxw@umich.edu | \n"
"| July 2017 | \n"
"| August 2017 | \n"
"| zhanxw.github.io/rvtests | \n"
"|----------------------------------------+ \n"
" \n";
......
......@@ -36,7 +36,7 @@ BASE += SingleDummy
OBJ = $(BASE:%=%.o)
OBJ_DBG = $(BASE:%=%_dbg.o)
DEFAULT_CXXFLAGS ?= -std=c++0x -I. $(addprefix -I, $(THIRD_INC)) -DGIT_VERSION="\"$(GIT_VERSION)\"" -fopenmp -Wno-error=int-in-bool-context
DEFAULT_CXXFLAGS ?= -std=c++0x -I. $(addprefix -I, $(THIRD_INC)) -DGIT_VERSION="\"$(GIT_VERSION)\"" -fopenmp
# enable read over HTTP and FTP
DEFAULT_CXXFLAGS += -D_USE_KNETFILE
......
......@@ -376,6 +376,7 @@ void makeVariableThreshodlGenotype(
makeVariableThreshodlGenotype(dc, in, freqIn, out, freqOut, collapseFunc);
}
#if 0
void SingleVariantScoreTest::calculateConstant(Matrix& phenotype) {
int nCase = 0;
int nCtrl = 0;
......@@ -393,14 +394,17 @@ void SingleVariantScoreTest::calculateConstant(Matrix& phenotype) {
alpha = 500.;
}
obtainB(alpha, &this->b);
fprintf(stderr, "alpha = %g, b = %g\n", alpha, b);
}
#endif
void MetaScoreTest::MetaFamBinary::calculateB() {
obtainB(this->alpha, &this->b);
// fprintf(stderr, "alpha = %g, b = %g\n", alpha, b);
return;
}
void MetaScoreTest::MetaUnrelatedBinary::calculateB() {
void MetaCovTest::MetaCovFamBinary::calculateB() {
obtainB(this->alpha, &this->b);
return;
}
......
......@@ -296,7 +296,7 @@ class SingleVariantScoreTest : public ModelFitter {
warnOnce("Single variant score test failed in fitting null model.");
return -1;
}
calculateConstant(phenotype);
// calculateConstant(phenotype);
needToFitNullModel = false;
}
fitOK = logistic.TestCovariate(cov, pheno, genotype);
......@@ -347,19 +347,20 @@ class SingleVariantScoreTest : public ModelFitter {
result.updateValue("DIRECTION",
logistic.GetU()[0][0] > 0 ? "+" : "-");
}
if (v > 0 && b > 0) {
result.updateValue("EFFECT", u / v / b);
result.updateValue("SE", 1.0 / sqrt(v) / b);
if (v > 0) {
result.updateValue("EFFECT", u / v);
result.updateValue("SE", 1.0 / sqrt(v));
}
result.updateValue("PVALUE", logistic.GetPvalue());
}
}
result.writeValueLine(fp);
}
void calculateConstant(Matrix& phenotype);
// don't need this
// void calculateConstant(Matrix& phenotype);
private:
double b; // a constant
// double b; // a constant
double af;
int nSample;
Vector pheno;
......@@ -3677,12 +3678,15 @@ class MetaScoreTest : public ModelFitter {
++nCtrl;
}
}
/*
no need to adjust for logistic regression
if (nCtrl > 0) {
alpha = log(1.0 * nCase / nCtrl);
} else {
alpha = 500.;
}
calculateB();
*/
}
// fit null model
bool fitOK = logistic.FitNullModel(cov, pheno, 100);
......@@ -3742,8 +3746,8 @@ class MetaScoreTest : public ModelFitter {
double GetV() { return logistic.GetV()[0][0]; }
double GetEffect() {
if (!useMLE) {
if (logistic.GetU()[0][0] != 0.0 && this->b != 0.0) {
return logistic.GetU()[0][0] / logistic.GetV()[0][0] / b;
if (logistic.GetU()[0][0] != 0.0) {
return logistic.GetU()[0][0] / logistic.GetV()[0][0];
}
} else {
return logisticAlt.GetCovEst()[1];
......@@ -3752,19 +3756,19 @@ class MetaScoreTest : public ModelFitter {
}
double GetEffectSE() {
const double v = GetV();
return v != 0.0 ? 1.0 / sqrt(v) / b : 0.0;
return v != 0.0 ? 1.0 / sqrt(v) : 0.0;
}
double GetPvalue() { return logistic.GetPvalue(); }
private:
void calculateB();
// private:
// void calculateB();
private:
LogisticRegressionScoreTest logistic;
Vector pheno;
Matrix X; // intercept, cov(optional) and genotype
double alpha;
double b;
// double alpha;
// double b;
bool useMLE;
LogisticRegression logisticAlt;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment