dev edits

for multiphen implementation

src/Data.cpp

@@ -3490,6 +3490,7 @@ void Data::compute_tests_mt_multiphen(int const& chrom, vector<uint64> indices,v
       MultiPhen mphen_i = mphen;
 
       // run MultiPhen test & save summary stats
+      /* cout << snpinfo[snp_index].ID << endl; */
       mphen_i.run(Gmat, pheno_data.cov_phenotypes, pheno_data.new_cov.cols() - 1, params.n_pheno); // the last 2 arg.: #cov excluding intercept; #phenotypes
 
       /* string tmp_str = mphen_i.print_sumstats(isnp, snp_index, test_string + params.condtl_suff, model_type + params.condtl_suff, block_info, snpinfo, &params); */
@@ -3543,26 +3544,31 @@ void Data::prep_multiphen()
 {
   // user parameters 
   mphen.test = params.multiphen_test;
+  mphen.optim = params.multiphen_optim;
   mphen.pval_thr = params.multiphen_thr;
-  mphen.firth_binom = (params.multiphen_firth_mult <= 0);
   mphen.firth_mult = params.multiphen_firth_mult;
+  mphen.firth_binom = (params.multiphen_firth_mult > 0);
+  mphen.firth_multinom = (params.multiphen_firth_mult > 0);
   mphen.tol = params.multiphen_tol;
   mphen.trace = params.multiphen_trace;
   mphen.verbose = params.multiphen_verbose;
   // parameters for model fitting
-  mphen.optim = "WeightHalving";
-  mphen.maxit = 250; mphen.maxit2 = 25;
-  mphen.check_step = (params.multiphen_maxstep <= 0);
+  mphen.maxit = params.multiphen_maxit; mphen.maxit2 = params.multiphen_maxit2;
+  mphen.strict = params.multiphen_strict;
+  mphen.check_step = (params.multiphen_maxstep > 0);
   mphen.max_step = params.multiphen_maxstep; 
   mphen.reuse_start = true;
-  mphen.approx_offset = params.multiphen_approx_offset;
+  mphen.mac_approx_offset = params.multiphen_approx_offset;
+  mphen.pseudo_stophalf = params.multiphen_pseudo_stophalf;
+  mphen.reset_start = params.multiphen_reset_start;
+  mphen.offset_mode = params.multiphen_offset;
 
   // prepare new matrix of covariates X + matrix of phenotypes Y
   unsigned int n_samples = pheno_data.new_cov.rows(), n_cov1 = pheno_data.new_cov.cols();
   unsigned int n_cov = n_cov1 - 1;
   unsigned int n_phen = params.n_pheno;
 
-  pheno_data.cov_phenotypes.resize(n_samples, n_cov + n_phen + 2); // +2 intercepts
+  pheno_data.cov_phenotypes.resize(n_samples, n_cov + 2*n_phen + 2); // +2 intercepts
   // column # 1 = Intercept
   pheno_data.cov_phenotypes.col(0) = ArrayXd::Constant(n_samples, 1.0);
   // next n_cov columns = covariates **without** intercept
@@ -3573,6 +3579,16 @@ void Data::prep_multiphen()
   // next n_phen columns = phenotypes (skipped here & to be filled in for each chr.)
   // next & the last column = Intercept
   pheno_data.cov_phenotypes.rightCols(1) = ArrayXd::Constant(n_samples, 1.0);
+
+  // v2
+  if(!params.strict_mode) throw std::runtime_error("--strict mode is required for MultiPhen test");
+
+  VectorXb Mask = pheno_data.masked_indivs.col(0);
+  for(unsigned int i = 1; i < pheno_data.masked_indivs.cols(); i++) {
+    Mask.col(0).array() = Mask.col(0).array() || pheno_data.masked_indivs.col(i).array();
+  }
+  pheno_data.cov_phenotypes.array().colwise() *= Mask.array().cast<double>().array();
+  mphen.setup_x(Mask, pheno_data.cov_phenotypes, n_cov, n_phen, true, false); // (ignored by MultiPhen) pos_intercept_first = true, pos_phen_first = false
 }
 
 void Data::set_multiphen()
@@ -3582,12 +3598,21 @@ void Data::set_multiphen()
   unsigned int n_phen = params.n_pheno;
 
   if(pheno_data.cov_phenotypes.rows() != n_samples) throw std::runtime_error("#rows in cov_phenotypes");
-  if(pheno_data.cov_phenotypes.cols() != n_cov + n_phen + 2) throw std::runtime_error("#rows in cov_phenotypes");
+  if(pheno_data.cov_phenotypes.cols() != n_cov + 2*n_phen + 2) throw std::runtime_error("#rows in cov_phenotypes");
 
-  pheno_data.cov_phenotypes.rightCols(n_phen + 1).leftCols(n_phen) = res;
+  // v2
+  for(unsigned i = n_cov1, k = 0; k < n_phen; i++, k++) {
+    pheno_data.cov_phenotypes.col(i) = mphen.Mask.select(res.col(k), 0.0);
+  }
+  for(unsigned i = n_cov1 + n_phen + 1, k = 0; k < n_phen; i++, k++) {
+    pheno_data.cov_phenotypes.col(i) = mphen.Mask.select(res.col(k), 0.0);
+  }
+  /* pheno_data.cov_phenotypes.rightCols(n_phen + 1).leftCols(n_phen) = res; */
+  /* pheno_data.cov_phenotypes.array().colwise() *= mphen.Mask.array().cast<double>().array(); */
 
-  if(!params.strict_mode) throw std::runtime_error("--strict mode is required for MultiPhen test");
-  mphen.setup_x(pheno_data.masked_indivs.col(0), pheno_data.cov_phenotypes, n_cov, n_phen, true, false); // (ignored by MultiPhen) pos_intercept_first = true, pos_phen_first = false
+  // v1
+  /* if(!params.strict_mode) throw std::runtime_error("--strict mode is required for MultiPhen test"); */
+  /* mphen.setup_x(pheno_data.masked_indivs.col(0), pheno_data.cov_phenotypes, n_cov, n_phen, true, false); // (ignored by MultiPhen) pos_intercept_first = true, pos_phen_first = false */
 }
 
 /////////////////////////////////////////////////