[Genabel-commits] r1031 - pkg/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 28 00:59:11 CET 2012
Author: lckarssen
Date: 2012-11-28 00:59:11 +0100 (Wed, 28 Nov 2012)
New Revision: 1031
Modified:
pkg/ProbABEL/src/command_line_settings.h
pkg/ProbABEL/src/coxph_data.cpp
pkg/ProbABEL/src/data.cpp
pkg/ProbABEL/src/main.cpp
pkg/ProbABEL/src/maskedmatrix.cpp
pkg/ProbABEL/src/phedata.cpp
pkg/ProbABEL/src/phedata.h
pkg/ProbABEL/src/reg1.cpp
pkg/ProbABEL/src/regdata.cpp
Log:
ProbABEL:
- Many code layout improvements
- Fixed a couple of cppcheck warnings, mostly related to initialisation of default values in classes. For some of those cppcheck recommends using initialisation lists. Why cppcheck only recommends this for some variables is unclear to me. When adding -Weffc++ to the g++ flags you will see more warnings about initialisation lists.
Modified: pkg/ProbABEL/src/command_line_settings.h
===================================================================
--- pkg/ProbABEL/src/command_line_settings.h 2012-11-27 23:51:20 UTC (rev 1030)
+++ pkg/ProbABEL/src/command_line_settings.h 2012-11-27 23:59:11 UTC (rev 1031)
@@ -45,38 +45,38 @@
public:
cmdvars()
{
- program_name = NULL;
+ program_name = NULL;
- std::fill_n(neco, 3, 0);
- phefilename = NULL;
- mlinfofilename = NULL;
- genfilename = NULL;
- mapfilename = NULL;
- outfilename = NULL;
- inverse_filename = NULL;
+ std::fill_n(neco, 3, 0);
+ phefilename = NULL;
+ mlinfofilename = NULL;
+ genfilename = NULL;
+ mapfilename = NULL;
+ outfilename = NULL;
+ inverse_filename = NULL;
- sep = " ";
- nohead = 0;
- score = 0;
- npeople = -1;
- ngpreds = 1;
- interaction = 0;
- interaction_excluded = 0;
- is_interaction_excluded = false; //Oh Holy Matrix, forgive me for this!
- robust = 0;
- chrom = "-1";
- str_genfilename = "";
+ sep = " ";
+ nohead = 0;
+ score = 0;
+ npeople = -1;
+ ngpreds = 1;
+ interaction = 0;
+ interaction_excluded = 0;
+ is_interaction_excluded = false; //Oh Holy Matrix, forgive me for this!
+ robust = 0;
+ chrom = "-1";
+ str_genfilename = "";
- iscox = false;
- isFVF = 0;
+ isFVF = 0;
- skipd = 2;
- allcov = 0;
+ skipd = 2;
+ allcov = 0;
#if COXPH
- noutcomes = 2;
- iscox=true;
+ noutcomes = 2;
+ iscox=true;
#else
- noutcomes = 1;
+ noutcomes = 1;
+ iscox = false;
#endif
}
Modified: pkg/ProbABEL/src/coxph_data.cpp
===================================================================
--- pkg/ProbABEL/src/coxph_data.cpp 2012-11-27 23:51:20 UTC (rev 1030)
+++ pkg/ProbABEL/src/coxph_data.cpp 2012-11-27 23:59:11 UTC (rev 1031)
@@ -44,14 +44,13 @@
return -9;
}
-coxph_data::coxph_data(const coxph_data &obj)
+coxph_data::coxph_data(const coxph_data &obj) : sstat(obj.sstat)
{
nids = obj.nids;
ncov = obj.ncov;
ngpreds = obj.ngpreds;
weights = obj.weights;
stime = obj.stime;
- sstat = obj.sstat;
offset = obj.offset;
strata = obj.strata;
X = obj.X;
@@ -60,23 +59,27 @@
for (int i = 0; i < nids; i++)
masked_data[i] = 0;
}
+
coxph_data::coxph_data(phedata &phed, gendata &gend, int snpnum)
{
nids = gend.nids;
masked_data = new unsigned short int[nids];
for (int i = 0; i < nids; i++)
masked_data[i] = 0;
+
ngpreds = gend.ngpreds;
if (snpnum >= 0)
ncov = phed.ncov + ngpreds;
else
ncov = phed.ncov;
+
if (phed.noutcomes != 2)
{
std::cerr << "coxph_data: number of outcomes should be 2 (now: "
<< phed.noutcomes << ")\n";
exit(1);
}
+
// X.reinit(nids,(ncov+1));
X.reinit(nids, ncov);
stime.reinit(nids, 1);
@@ -275,21 +278,26 @@
}
void coxph_reg::estimate(coxph_data &cdatain, int verbose, int maxiter,
- double eps, double tol_chol, int model, int interaction, int ngpreds,
- bool iscox, int nullmodel)
+ double eps, double tol_chol, int model,
+ int interaction, int ngpreds, bool iscox,
+ int nullmodel)
{
coxph_data cdata = cdatain.get_unmasked_data();
mematrix<double> X = t_apply_model(cdata.X, model, interaction, ngpreds,
- iscox, nullmodel);
+ iscox, nullmodel);
// X.print();
int length_beta = X.nrow;
beta.reinit(length_beta, 1);
sebeta.reinit(length_beta, 1);
- mematrix<double> newoffset = cdata.offset;
- newoffset = cdata.offset - (cdata.offset).column_mean(0);
+ mematrix<double> newoffset = cdata.offset -
+ (cdata.offset).column_mean(0);
mematrix<double> means(X.nrow, 1);
+
for (int i = 0; i < X.nrow; i++)
+ {
beta[i] = 0.;
+ }
+
mematrix<double> u(X.nrow, 1);
mematrix<double> imat(X.nrow, X.nrow);
double work[X.ncol * 2 + 2 * (X.nrow) * (X.nrow) + 3 * (X.nrow)];
@@ -304,6 +312,7 @@
means.data.data(), beta.data.data(), u.data.data(),
imat.data.data(), loglik_int, &flag, work, &eps, &tol_chol,
&sctest);
+
for (int i = 0; i < X.nrow; i++)
{
sebeta[i] = sqrt(imat.get(i, i));
Modified: pkg/ProbABEL/src/data.cpp
===================================================================
--- pkg/ProbABEL/src/data.cpp 2012-11-27 23:51:20 UTC (rev 1030)
+++ pkg/ProbABEL/src/data.cpp 2012-11-27 23:59:11 UTC (rev 1031)
@@ -168,9 +168,8 @@
//_________________________________________Maksim_start
-InvSigma::InvSigma(const char * filename_, phedata * phe)
+InvSigma::InvSigma(const char * filename_, phedata * phe) : filename(filename_)
{
- filename = filename_;
npeople = phe->nids;
std::ifstream myfile(filename_);
char * line = new char[MAXIMUM_PEOPLE_AMOUNT];
@@ -178,7 +177,7 @@
matrix.reinit(npeople, npeople);
-//idnames[k], if (allmeasured[i]==1)
+ // idnames[k], if (allmeasured[i]==1)
if (myfile.is_open())
{
@@ -193,9 +192,9 @@
if (phe->idnames[row] != id)
{
std::cerr << "error:in row " << row << " id="
- << phe->idnames[row]
- << " in inverse variance matrix but id=" << id
- << " must be there. Wrong inverse variance matrix (only measured id must be there)\n";
+ << phe->idnames[row]
+ << " in inverse variance matrix but id=" << id
+ << " must be there. Wrong inverse variance matrix (only measured id must be there)\n";
exit(1);
}
unsigned col = 0;
Modified: pkg/ProbABEL/src/main.cpp
===================================================================
--- pkg/ProbABEL/src/main.cpp 2012-11-27 23:51:20 UTC (rev 1030)
+++ pkg/ProbABEL/src/main.cpp 2012-11-27 23:59:11 UTC (rev 1031)
@@ -74,14 +74,16 @@
}
void open_files_for_output(std::vector<std::ofstream*>& outfile,
- std::string& outfilename_str)
+ std::string& outfilename_str)
{
//create a list of filenames
const int amount_of_files = 5;
- std::string filenames[amount_of_files] = { outfilename_str + "_2df.out.txt",
- outfilename_str + "_add.out.txt", outfilename_str
- + "_domin.out.txt", outfilename_str + "_recess.out.txt",
- outfilename_str + "_over_domin.out.txt" };
+ std::string filenames[amount_of_files] =
+ { outfilename_str + "_2df.out.txt",
+ outfilename_str + "_add.out.txt",
+ outfilename_str + "_domin.out.txt",
+ outfilename_str + "_recess.out.txt",
+ outfilename_str + "_over_domin.out.txt" };
for (int i = 0; i < amount_of_files; i++)
{
@@ -90,27 +92,32 @@
if (!outfile[i]->is_open())
{
std::cerr << "Cannot open file for writing: " << filenames[i]
- << "\n";
+ << "\n";
exit(1);
}
}
}
-int create_phenoytype(phedata& phd, cmdvars& input_var)
+int create_phenotype(phedata& phd, cmdvars& input_var)
{
phd.set_is_interaction_excluded(input_var.isIsInteractionExcluded());
- phd.setphedata(input_var.getPhefilename(), input_var.getNoutcomes(),
- input_var.getNpeople(), input_var.getInteraction(),
- input_var.isIscox());
+ phd.setphedata(input_var.getPhefilename(),
+ input_var.getNoutcomes(),
+ input_var.getNpeople(),
+ input_var.getInteraction(),
+ input_var.isIscox());
+
int interaction_cox = input_var.getInteraction();
#if COXPH
interaction_cox--;
#endif
- if (input_var.getInteraction() < 0 || input_var.getInteraction() > phd.ncov
- || interaction_cox > phd.ncov)
+
+ if (input_var.getInteraction() < 0 ||
+ input_var.getInteraction() > phd.ncov ||
+ interaction_cox > phd.ncov)
{
- std::cerr << "error: Interaction parameter is out of range (ineraction="
- << input_var.getInteraction() << ") \n";
+ std::cerr << "error: Interaction parameter is out of range (interaction="
+ << input_var.getInteraction() << ") \n";
exit(1);
}
@@ -132,12 +139,15 @@
{
for (unsigned int i = 0; i < outfile.size(); i++)
{
- (*outfile[i]) << "name" << input_var.getSep() << "A1"
- << input_var.getSep() << "A2" << input_var.getSep() << "Freq1"
- << input_var.getSep() << "MAF" << input_var.getSep()
- << "Quality" << input_var.getSep() << "Rsq"
- << input_var.getSep() << "n" << input_var.getSep()
- << "Mean_predictor_allele";
+ (*outfile[i]) << "name" << input_var.getSep()
+ << "A1" << input_var.getSep()
+ << "A2" << input_var.getSep()
+ << "Freq1" << input_var.getSep()
+ << "MAF" << input_var.getSep()
+ << "Quality" << input_var.getSep()
+ << "Rsq" << input_var.getSep()
+ << "n" << input_var.getSep()
+ << "Mean_predictor_allele";
if (input_var.getChrom() != "-1")
(*outfile[i]) << input_var.getSep() << "chrom";
if (input_var.getMapfilename() != NULL)
@@ -148,9 +158,10 @@
{
for (unsigned int file = 0; file < outfile.size(); file++)
for (int i = 0; i < phd.n_model_terms - 1; i++)
- *outfile[file] << input_var.getSep() << "beta_"
- << phd.model_terms[i] << input_var.getSep() << "sebeta_"
- << phd.model_terms[i];
+ *outfile[file] << input_var.getSep()
+ << "beta_" << phd.model_terms[i]
+ << input_var.getSep() << "sebeta_"
+ << phd.model_terms[i];
}
}
@@ -160,48 +171,51 @@
create_start_of_header(outfile, input_var, phd);
*outfile[0] << input_var.getSep() << "beta_SNP_A1A2" << input_var.getSep()
- << "beta_SNP_A1A1" << input_var.getSep() << "sebeta_SNP_A1A2"
- << input_var.getSep() << "sebeta_SNP_A1A1";
+ << "beta_SNP_A1A1" << input_var.getSep() << "sebeta_SNP_A1A2"
+ << input_var.getSep() << "sebeta_SNP_A1A1";
*outfile[1] << input_var.getSep() << "beta_SNP_addA1" << input_var.getSep()
- << "sebeta_SNP_addA1";
+ << "sebeta_SNP_addA1";
*outfile[2] << input_var.getSep() << "beta_SNP_domA1" << input_var.getSep()
- << "sebeta_SNP_domA1";
+ << "sebeta_SNP_domA1";
*outfile[3] << input_var.getSep() << "beta_SNP_recA1" << input_var.getSep()
- << "sebeta_SNP_recA1";
+ << "sebeta_SNP_recA1";
*outfile[4] << input_var.getSep() << "beta_SNP_odom" << input_var.getSep()
- << "sebeta_SNP_odom";
+ << "sebeta_SNP_odom";
if (input_var.getInteraction() != 0)
{
//Han Chen
*outfile[0] << input_var.getSep() << "beta_SNP_A1A2_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "sebeta_SNP_A1A2_" << phd.model_terms[interaction_cox]
- << input_var.getSep() << "beta_SNP_A1A1_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "sebeta_SNP_A1A1_" << phd.model_terms[interaction_cox];
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep() << "sebeta_SNP_A1A2_"
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep() << "beta_SNP_A1A1_"
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep() << "sebeta_SNP_A1A1_"
+ << phd.model_terms[interaction_cox];
#if !COXPH
if (input_var.getInverseFilename() == NULL && !input_var.getAllcov())
{
*outfile[0] << input_var.getSep() << "cov_SNP_A1A2_int_SNP_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "cov_SNP_A1A1_int_SNP_"
- << phd.model_terms[interaction_cox];
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep() << "cov_SNP_A1A1_int_SNP_"
+ << phd.model_terms[interaction_cox];
}
#endif
//Oct 26, 2009
for (unsigned int file = 1; file < outfile.size(); file++)
{
*outfile[file] << input_var.getSep() << "beta_SNP_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "sebeta_SNP_" << phd.model_terms[interaction_cox];
+ << phd.model_terms[interaction_cox]
+ << input_var.getSep() << "sebeta_SNP_"
+ << phd.model_terms[interaction_cox];
//Han Chen
#if !COXPH
if (input_var.getInverseFilename() == NULL
- && !input_var.getAllcov())
+ && !input_var.getAllcov())
{
*outfile[file] << input_var.getSep() << "cov_SNP_int_SNP_"
- << phd.model_terms[interaction_cox];
+ << phd.model_terms[interaction_cox];
}
#endif
//Oct 26, 2009
@@ -215,17 +229,17 @@
}
void create_header2(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
- phedata phd, int interaction_cox)
+ phedata phd, int interaction_cox)
{
create_start_of_header(outfile, input_var, phd);
*outfile[0] << input_var.getSep() << "beta_SNP_add" << input_var.getSep()
- << "sebeta_SNP_add";
+ << "sebeta_SNP_add";
if (input_var.getInteraction() != 0)
{
*outfile[0] << input_var.getSep() << "beta_SNP_"
- << phd.model_terms[interaction_cox] << input_var.getSep()
- << "sebeta_SNP_" << phd.model_terms[interaction_cox];
+ << phd.model_terms[interaction_cox] << input_var.getSep()
+ << "sebeta_SNP_" << phd.model_terms[interaction_cox];
}
if (input_var.getInverseFilename() == NULL)
@@ -235,7 +249,7 @@
if (input_var.getInteraction() != 0 && !input_var.getAllcov())
{
*outfile[0] << input_var.getSep() << "cov_SNP_int_SNP_"
- << phd.model_terms[interaction_cox];
+ << phd.model_terms[interaction_cox];
}
#endif
*outfile[0] << input_var.getSep() << "loglik"; //"chi2_SNP";
@@ -258,7 +272,7 @@
mlinfo mli(input_var.getMlinfofilename(), input_var.getMapfilename());
int nsnps = mli.nsnps;
phedata phd;
- int interaction_cox = create_phenoytype(phd, input_var);
+ int interaction_cox = create_phenotype(phd, input_var);
//interaction--;
// if(input_var.getInverseFilename()!= NULL && phd.ncov > 1)
@@ -288,19 +302,20 @@
gendata gtd;
if (!input_var.getIsFvf())
- //use the non non filevector input format
+ // use the non-filevector input format
gtd.re_gendata(input_var.getGenfilename(), nsnps,
- input_var.getNgpreds(), phd.nids_all, phd.nids, phd.allmeasured,
- input_var.getSkipd(), phd.idnames);
+ input_var.getNgpreds(), phd.nids_all, phd.nids,
+ phd.allmeasured, input_var.getSkipd(), phd.idnames);
else
- //use the filevector input format (missing second last skipd parameter)
+ // use the filevector input format (missing second last skipd
+ // parameter)
gtd.re_gendata(input_var.getStrGenfilename(), nsnps,
- input_var.getNgpreds(), phd.nids_all, phd.nids, phd.allmeasured,
- phd.idnames);
+ input_var.getNgpreds(), phd.nids_all, phd.nids,
+ phd.allmeasured, phd.idnames);
std::cout << " loaded genotypic data ..." << std::flush;
/**
- if (input_var.getIsFvf())
+ if (input_var.getIsFvf())
gendata gtd (str_genfilename,nsnps,input_var.getNgpreds(),phd.nids_all,phd.allmeasured,phd.idnames);
else
gendata gtd (input_var.getGenfilename(),nsnps,input_var.getNgpreds(),phd.nids_all,phd.nids,phd.allmeasured,skipd,phd.idnames);
@@ -315,7 +330,9 @@
std::cout << " loaded null data ..." << std::flush;
#if LOGISTIC
logistic_reg nrd = logistic_reg(nrgd);
- nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
+ nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0,
+ input_var.getInteraction(), input_var.getNgpreds(),
+ invvarmatrix, input_var.getRobust(), 1);
#elif LINEAR
linear_reg nrd = linear_reg(nrgd);
@@ -323,11 +340,12 @@
std::cout << "[DEBUG] linear_reg nrd = linear_reg(nrgd); DONE.";
#endif
nrd.estimate(nrgd, 0, CHOLTOL, 0, input_var.getInteraction(),
- input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
+ input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
#elif COXPH
coxph_reg nrd(nrgd);
- nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0, input_var.getInteraction(), input_var.getNgpreds(), 1);
+ nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0,
+ input_var.getInteraction(), input_var.getNgpreds(), 1);
#endif
std::cout << " estimated null model ...";
@@ -358,12 +376,12 @@
} else //Only additive model. Only one output file
{
outfile.push_back(
- new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
+ new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
if (!outfile[0]->is_open())
{
std::cerr << "Cannot open file for writing: " << outfilename_str
- << "\n";
+ << "\n";
exit(1);
}
if (input_var.getNohead() != 1)
@@ -434,7 +452,8 @@
float snpdata2[gtd.nids];
if (input_var.getNgpreds() == 2)
{
- //freq = ((gtd.G).column_mean(csnp*2)*2.+(gtd.G).column_mean(csnp*2+1))/2.;
+ //freq = ((gtd.G).column_mean(csnp*2)*2. +
+ // (gtd.G).column_mean(csnp*2+1))/2.;
gtd.get_var(csnp * 2, snpdata1);
gtd.get_var(csnp * 2 + 1, snpdata2);
for (unsigned int ii = 0; ii < gtd.nids; ii++)
@@ -445,7 +464,7 @@
}
} else
{
- // freq = (gtd.G).column_mean(csnp)/2.;
+ // freq = (gtd.G).column_mean(csnp)/2.;
gtd.get_var(csnp, snpdata1);
for (unsigned int ii = 0; ii < gtd.nids; ii++)
if (!isnan(snpdata1[ii]))
@@ -467,17 +486,18 @@
//Write mlinfo to output:
for (unsigned int file = 0; file < outfile.size(); file++)
{
- *outfile[file] << mli.name[csnp] << input_var.getSep()
- << mli.A1[csnp] << input_var.getSep() << mli.A2[csnp]
- << input_var.getSep() << mli.Freq1[csnp]
- << input_var.getSep() << mli.MAF[csnp]
- << input_var.getSep() << mli.Quality[csnp]
- << input_var.getSep() << mli.Rsq[csnp]
- << input_var.getSep() << gcount << input_var.getSep()
- << freq;
+ *outfile[file] << mli.name[csnp]
+ << input_var.getSep() << mli.A1[csnp]
+ << input_var.getSep() << mli.A2[csnp]
+ << input_var.getSep() << mli.Freq1[csnp]
+ << input_var.getSep() << mli.MAF[csnp]
+ << input_var.getSep() << mli.Quality[csnp]
+ << input_var.getSep() << mli.Rsq[csnp]
+ << input_var.getSep() << gcount
+ << input_var.getSep() << freq;
if (input_var.getChrom() != "-1")
*outfile[file] << input_var.getSep()
- << input_var.getChrom();
+ << input_var.getChrom();
if (input_var.getMapfilename() != NULL)
*outfile[file] << input_var.getSep() << mli.map[csnp];
}
@@ -489,15 +509,23 @@
#if LOGISTIC
logistic_reg rd(rgd);
if (input_var.getScore())
- rd.score(nrd.residuals, rgd, 0, CHOLTOL, model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
+ rd.score(nrd.residuals, rgd, 0, CHOLTOL, model,
+ input_var.getInteraction(),
+ input_var.getNgpreds(),
+ invvarmatrix);
else
- rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model,input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust());
+ rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model,
+ input_var.getInteraction(),
+ input_var.getNgpreds(),
+ invvarmatrix,
+ input_var.getRobust());
#elif LINEAR
linear_reg rd(rgd);
if (input_var.getScore())
rd.score(nrd.residuals, rgd, 0, CHOLTOL, model,
- input_var.getInteraction(),
- input_var.getNgpreds(), invvarmatrix);
+ input_var.getInteraction(),
+ input_var.getNgpreds(),
+ invvarmatrix);
else
{
// rd.mmscore(rgd,0,CHOLTOL,model,input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
@@ -508,20 +536,22 @@
}
#elif COXPH
coxph_reg rd(rgd);
- rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model, input_var.getInteraction(), true, input_var.getNgpreds());
+ rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model,
+ input_var.getInteraction(), true,
+ input_var.getNgpreds());
#endif
if (!input_var.getAllcov() && model == 0
- && input_var.getInteraction() == 0)
+ && input_var.getInteraction() == 0)
start_pos = rd.beta.nrow - 2;
else if (!input_var.getAllcov() && model == 0
- && input_var.getInteraction() != 0)
+ && input_var.getInteraction() != 0)
start_pos = rd.beta.nrow - 4;
else if (!input_var.getAllcov() && model != 0
- && input_var.getInteraction() == 0)
+ && input_var.getInteraction() == 0)
start_pos = rd.beta.nrow - 1;
else if (!input_var.getAllcov() && model != 0
- && input_var.getInteraction() != 0)
+ && input_var.getInteraction() != 0)
start_pos = rd.beta.nrow - 2;
else
start_pos = 0;
@@ -529,13 +559,14 @@
for (int pos = start_pos; pos < rd.beta.nrow; pos++)
{
*beta_sebeta[model] << input_var.getSep()
- << rd.beta[pos] << input_var.getSep()
- << rd.sebeta[pos];
+ << rd.beta[pos]
+ << input_var.getSep()
+ << rd.sebeta[pos];
//Han Chen
#if !COXPH
if (input_var.getInverseFilename() == NULL
- && !input_var.getAllcov()
- && input_var.getInteraction() != 0)
+ && !input_var.getAllcov()
+ && input_var.getInteraction() != 0)
{
if (pos > start_pos)
{
@@ -544,9 +575,9 @@
if (pos > start_pos + 2)
{
*covvalue[model]
- << rd.covariance[pos - 3]
- << input_var.getSep()
- << rd.covariance[pos - 2];
+ << rd.covariance[pos - 3]
+ << input_var.getSep()
+ << rd.covariance[pos - 2];
}
} else
{
@@ -573,16 +604,16 @@
} else //beta, sebeta = nan
{
if (!input_var.getAllcov() && model == 0
- && input_var.getInteraction() == 0)
+ && input_var.getInteraction() == 0)
start_pos = rgd.X.ncol - 2;
else if (!input_var.getAllcov() && model == 0
- && input_var.getInteraction() != 0)
+ && input_var.getInteraction() != 0)
start_pos = rgd.X.ncol - 4;
else if (!input_var.getAllcov() && model != 0
- && input_var.getInteraction() == 0)
+ && input_var.getInteraction() == 0)
start_pos = rgd.X.ncol - 1;
else if (!input_var.getAllcov() && model != 0
- && input_var.getInteraction() != 0)
+ && input_var.getInteraction() != 0)
start_pos = rgd.X.ncol - 2;
else
start_pos = 0;
@@ -601,17 +632,17 @@
for (int pos = start_pos; pos < end_pos; pos++)
{
*beta_sebeta[model] << input_var.getSep() << "nan"
- << input_var.getSep() << "nan";
+ << input_var.getSep() << "nan";
}
//Han Chen
#if !COXPH
if (!input_var.getAllcov()
- && input_var.getInteraction() != 0)
+ && input_var.getInteraction() != 0)
{
if (model == 0)
{
- *covvalue[model] << "nan" << input_var.getSep()
- << "nan";
+ *covvalue[model] << "nan"
+ << input_var.getSep() << "nan";
} else
{
*covvalue[model] << "nan";
@@ -672,12 +703,15 @@
} else //Only additive model. Only one output file
{
//Write mlinfo to output:
- *outfile[0] << mli.name[csnp] << input_var.getSep() << mli.A1[csnp]
- << input_var.getSep() << mli.A2[csnp] << input_var.getSep();
- *outfile[0] << mli.Freq1[csnp] << input_var.getSep()
- << mli.MAF[csnp] << input_var.getSep() << mli.Quality[csnp]
- << input_var.getSep() << mli.Rsq[csnp]
- << input_var.getSep();
+ *outfile[0] << mli.name[csnp]
+ << input_var.getSep() << mli.A1[csnp]
+ << input_var.getSep() << mli.A2[csnp]
+ << input_var.getSep();
+ *outfile[0] << mli.Freq1[csnp]
+ << input_var.getSep() << mli.MAF[csnp]
+ << input_var.getSep() << mli.Quality[csnp]
+ << input_var.getSep() << mli.Rsq[csnp]
+ << input_var.getSep();
*outfile[0] << gcount << input_var.getSep() << freq;
if (input_var.getChrom() != "-1")
*outfile[0] << input_var.getSep() << input_var.getChrom();
@@ -689,9 +723,16 @@
#if LOGISTIC
logistic_reg rd(rgd);
if (input_var.getScore())
- rd.score(nrd.residuals, rgd, 0, CHOLTOL, model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
+ rd.score(nrd.residuals, rgd, 0, CHOLTOL, model,
+ input_var.getInteraction(),
+ input_var.getNgpreds(),
+ invvarmatrix);
else
- rd.estimate( rgd, 0, MAXITER, EPS, CHOLTOL, model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust());
+ rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model,
+ input_var.getInteraction(),
+ input_var.getNgpreds(),
+ invvarmatrix,
+ input_var.getRobust());
#elif LINEAR
//cout << (rgd.get_unmasked_data()).nids << " 1\n";
#if DEBUG
@@ -712,14 +753,15 @@
cout << CHOLTOL << " <-CHOLTOL\n";
cout << model << " <-model\n";
cout << input_var.getInteraction()
- << " <-input_var.getInteraction()\n";
+ << " <-input_var.getInteraction()\n";
cout << input_var.getNgpreds()
- << " <-input_var.getNgpreds()\n";
+ << " <-input_var.getNgpreds()\n";
invvarmatrix.print();
#endif
rd.score(nrd.residuals, rgd, 0, CHOLTOL, model,
- input_var.getInteraction(), input_var.getNgpreds(),
- invvarmatrix);
+ input_var.getInteraction(),
+ input_var.getNgpreds(),
+ invvarmatrix);
#if DEBUG
rd.beta.print();
cout << rd.chi2_score << " <-chi2_scoren\n";
@@ -731,28 +773,30 @@
#endif
} else
{
- // if(input_var.getInverseFilename()== NULL)
- // {
- //cout << (rgd.get_unmasked_data()).nids << " 3\n";
+ // if(input_var.getInverseFilename()== NULL)
+ // {
+ // cout << (rgd.get_unmasked_data()).nids << " 3\n";
#if DEBUG
cout << "rd.estimate\n";
cout << CHOLTOL << " <-CHOLTOL\n";
cout << model << " <-model\n";
cout << input_var.getInteraction()
- << " <-input_var.getInteraction()\n";
+ << " <-input_var.getInteraction()\n";
cout << input_var.getNgpreds()
- << " <-input_var.getNgpreds()\n";
+ << " <-input_var.getNgpreds()\n";
cout << input_var.getRobust()
- << " <-input_var.getRobust()\n";
+ << " <-input_var.getRobust()\n";
cout << "start invarmatrix\n";
invvarmatrix.print();
cout << "end invarmatrix\n";
cout << rgd.is_interaction_excluded
- << " <-rgd.is_interaction_excluded\n";
+ << " <-rgd.is_interaction_excluded\n";
#endif
rd.estimate(rgd, 0, CHOLTOL, model,
- input_var.getInteraction(), input_var.getNgpreds(),
- invvarmatrix, input_var.getRobust());
+ input_var.getInteraction(),
+ input_var.getNgpreds(),
+ invvarmatrix,
+ input_var.getRobust());
#if DEBUG
cout << "rd.beta\n";
@@ -768,22 +812,25 @@
cout << rd.sigma2 << " <-sigma2\n";
#endif
//cout << (rgd.get_unmasked_data()).nids << " 4\n";
- // }
- // else
- // {
- // rd.mmscore(rgd,0,CHOLTOL,model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/genabel -r 1031
More information about the Genabel-commits
mailing list