Writing output local to function

2 files changed, 76 insertions, 7 deletions

diff --git a/src/gemma.cpp b/src/gemma.cpp
index c7f7ab6..f7ebdce 100644
--- a/src/gemma.cpp
+++ b/src/gemma.cpp
@@ -2846,7 +2846,7 @@ void GEMMA::BatchRun(PARAM &cPar) {
             cLmm.mdb_calc_gwa(U, eval, UtW, &UtY_col.vector, W,
                               &Y_col.vector, cPar.setGWASnps);
           }
-          else
+          else {
             if (!cPar.file_bfile.empty()) {
               // PLINK analysis
               if (cPar.file_gxe.empty()) {
@@ -2869,8 +2869,9 @@ void GEMMA::BatchRun(PARAM &cPar) {
                                       &Y_col.vector, env);
               }
             }
-          cLmm.WriteFiles(); // write output files
-          cLmm.CopyToParam(cPar);
+            cLmm.WriteFiles(); // write output files
+            cLmm.CopyToParam(cPar);
+          }
         } else {
           debug_msg("fit mvLMM (multiple phenotypes)");
           MVLMM cMvlmm;
diff --git a/src/lmm.cpp b/src/lmm.cpp
index 0a5d8e6..8a57885 100644
--- a/src/lmm.cpp
+++ b/src/lmm.cpp
@@ -1876,12 +1876,30 @@ void LMM::Analyze(std::function< SnpNameValues(size_t) >& fetch_snp,
  */
 
 
+// Results for LMM.
+class SUMSTAT2 {
+public:
+  double beta;         // REML estimator for beta.
+  double se;           // SE for beta.
+  double lambda_remle; // REML estimator for lambda.
+  double lambda_mle;   // MLE estimator for lambda.
+  double p_wald;       // p value from a Wald test.
+  double p_lrt;        // p value from a likelihood ratio test.
+  double p_score;      // p value from a score test.
+  double logl_H1;      // log likelihood under the alternative
+                       // hypothesis as a measure of goodness of fit,
+                       // see https://github.com/genetics-statistics/GEMMA/issues/81
+};
+
+
+
 void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
                       const gsl_matrix *U, const gsl_vector *eval,
                       const gsl_matrix *UtW, const gsl_vector *Uty,
                       const gsl_matrix *W, const gsl_vector *y,
                       const set<string> gwasnps,
                       size_t num_markers) {
+  vector<SUMSTAT2> sumstat2;
   clock_t time_start = clock();
 
   checkpoint_nofile("start-lmm-mdb-analyze");
@@ -1985,9 +2003,9 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
       time_opt += (clock() - time_start) / (double(CLOCKS_PER_SEC) * 60.0);
 
       // Store summary data.
-      SUMSTAT SNPs = {beta,   se,    lambda_remle, lambda_mle,
+      SUMSTAT2 marker = {beta,   se,    lambda_remle, lambda_mle,
                       p_wald, p_lrt, p_score, logl_H1};
-      sumStat.push_back(SNPs);
+      sumstat2.push_back(marker);
     }
   };
 
@@ -2005,7 +2023,7 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
   for (size_t t = 0; t < num_markers; ++t) {
   // for (size_t t = 0; t < 2; ++t) {
     if (t % progress_step == 0 || t == (num_markers - 1)) {
-      ProgressBar("Reading SNPs", t, num_markers - 1);
+      ProgressBar("Reading markers", t, num_markers - 1);
     }
     // if (indicator_snp[t] == 0)
     // continue;
@@ -2072,9 +2090,59 @@ void LMM::mdb_analyze(std::function< SnpNameValues2(size_t) >& fetch_snp,
     batch_compute(c % msize);
   ProgressBar("Reading markers", num_markers - 1, num_markers - 1);
   cout << endl;
-  cout << "Counted markers " << c << " sumStat " << sumStat.size() << endl;
+  cout << "Counted markers " << c << " sumStat " << sumstat2.size() << endl;
   checkpoint_nofile("end-lmm-mdb-analyze");
 
+  string file_str;
+  debug_msg("LMM::WriteFiles");
+  file_str = path_out + "/" + file_out;
+  file_str += ".assoc.txt";
+  checkpoint("lmm-write-files",file_str);
+
+  ofstream outfile(file_str.c_str(), ofstream::out);
+  if (!outfile) {
+    cout << "error writing file: " << file_str << endl;
+    return;
+  }
+
+  auto sumstats = [&] (SUMSTAT2 st) {
+    outfile << scientific << setprecision(6);
+
+    if (a_mode != M_LMM2) {
+      outfile << st.beta << "\t";
+      outfile << st.se << "\t";
+    }
+
+    if (a_mode != M_LMM3 && a_mode != M_LMM9)
+      outfile << st.logl_H1 << "\t";
+
+    switch(a_mode) {
+    case M_LMM1:
+      outfile << st.lambda_remle << "\t"
+              << st.p_wald << endl;
+      break;
+    case M_LMM2:
+    case M_LMM9:
+      outfile << st.lambda_mle << "\t"
+              << st.p_lrt << endl;
+      break;
+    case M_LMM3:
+      outfile << st.p_score << endl;
+      break;
+    case M_LMM4:
+      outfile << st.lambda_remle << "\t"
+              << st.lambda_mle << "\t"
+              << st.p_wald << "\t"
+              << st.p_lrt << "\t"
+              << st.p_score << endl;
+      break;
+    }
+  };
+
+  for (auto &s : sumstat2) {
+    sumstats(s);
+  }
+
   gsl_vector_safe_free(x);
   gsl_vector_safe_free(x_miss);
   gsl_vector_safe_free(Utx);