+vector<double> MothurMetastats::permuted_pvalues(vector< vector<double> >& Imatrix, vector<double>& tstats, vector< vector<double> >& Fmatrix) {
+ try {
+ //# matrix stores tstats for each taxa(row) for each permuted trial(column)
+ vector<double> ps; ps.resize(row, 0.0); //# to store the pvalues
+ vector< vector<double> > permuted_ttests; permuted_ttests.resize(numPermutations);
+ for (int i = 0; i < numPermutations; i++) { permuted_ttests[i].resize(row, 0.0); }
+
+ //# calculate null version of tstats using B permutations.
+ for (int i = 0; i < numPermutations; i++) {
+ permuted_ttests[i] = permute_and_calc_ts(Imatrix);
+ }
+
+ //# calculate each pvalue using the null ts
+ if ((secondGroupingStart) < 8 || (column-secondGroupingStart) < 8){
+ vector< vector<double> > cleanedpermuted_ttests; cleanedpermuted_ttests.resize(numPermutations); //# the array pooling just the frequently observed ts
+ //# then pool the t's together!
+ //# count how many high freq taxa there are
+ int hfc = 1;
+ for (int i = 0; i < row; i++) { // # for each taxa
+ double group1Total = 0.0; double group2Total = 0.0;
+ for(int j = 0; j < secondGroupingStart; j++) { group1Total += Fmatrix[i][j]; }
+ for(int j = secondGroupingStart; j < column; j++) { group2Total += Fmatrix[i][j]; }
+
+ if (group1Total >= secondGroupingStart || group2Total >= (column-secondGroupingStart)){
+ hfc++;
+ for (int j = 0; j < numPermutations; j++) { cleanedpermuted_ttests[j].push_back(permuted_ttests[j][i]); }
+ }
+ }
+
+ //#now for each taxa
+ for (int i = 0; i < row; i++) {
+ //number of cleanedpermuted_ttests greater than tstat[i]
+ int numGreater = 0;
+ for (int j = 0; j < numPermutations; j++) {
+ for (int k = 0; k < hfc; k++) {
+ if (cleanedpermuted_ttests[j][k] > abs(tstats[i])) { numGreater++; }
+ }
+ }
+
+ ps[i] = (1/(double)(numPermutations*hfc))*numGreater;
+ }
+ }else{
+ for (int i = 0; i < row; i++) {
+ //number of permuted_ttests[i] greater than tstat[i] //(sum(permuted_ttests[i,] > abs(tstats[i]))+1)
+ int numGreater = 1;
+ for (int j = 0; j < numPermutations; j++) { if (permuted_ttests[j][i] > abs(tstats[i])) { numGreater++; } }
+ ps[i] = (1/(double)(numPermutations+1))*numGreater;
+ }
+ }
+
+ return ps;
+
+ }catch(exception& e) {
+ m->errorOut(e, "MothurMetastats", "permuted_pvalues");
+ exit(1);
+ }
+}
+/***********************************************************/
+vector<double> MothurMetastats::permute_and_calc_ts(vector< vector<double> >& Imatrix) {
+ try {
+ vector< vector<double> > permutedMatrix = Imatrix;
+
+ //randomize columns, ie group abundances.
+ for (int i = 0; i < permutedMatrix.size(); i++) { random_shuffle(permutedMatrix[i].begin(), permutedMatrix[i].end()); }
+
+ //calc ts
+ vector< vector<double> > C1; C1.resize(row);
+ for (int i = 0; i < row; i++) { C1[i].resize(3, 0.0); } // statistic profiles for class1 and class 2
+ vector< vector<double> > C2; C2.resize(row); // mean[1], variance[2], standard error[3]
+ for (int i = 0; i < row; i++) { C2[i].resize(3, 0.0); }
+ vector<double> Ts; Ts.resize(row, 0.0); // a place to store the true t-statistics
+
+ //#*************************************
+ //# generate statistics mean, var, stderr
+ //#*************************************
+ for(int i = 0; i < row; i++){ // for each taxa
+ //# find the mean of each group
+ double g1Total = 0.0; double g2Total = 0.0;
+ for (int j = 0; j < secondGroupingStart; j++) { g1Total += permutedMatrix[i][j]; }
+ C1[i][0] = g1Total/(double)(secondGroupingStart);
+ for (int j = secondGroupingStart; j < column; j++) { g2Total += permutedMatrix[i][j]; }
+ C2[i][0] = g2Total/(double)(column-secondGroupingStart);
+
+ //# find the variance of each group
+ double g1Var = 0.0; double g2Var = 0.0;
+ for (int j = 0; j < secondGroupingStart; j++) { g1Var += pow((permutedMatrix[i][j]-C1[i][0]), 2); }
+ C1[i][1] = g1Var/(double)(secondGroupingStart-1);
+ for (int j = secondGroupingStart; j < column; j++) { g2Var += pow((permutedMatrix[i][j]-C2[i][0]), 2); }
+ C2[i][1] = g2Var/(double)(column-secondGroupingStart-1);
+
+ //# find the std error of each group -std err^2 (will change to std err at end)
+ C1[i][2] = C1[i][1]/(double)(secondGroupingStart);
+ C2[i][2] = C2[i][1]/(double)(column-secondGroupingStart);
+ }
+
+ //#*************************************
+ //# two sample t-statistics
+ //#*************************************
+ for(int i = 0; i < row; i++){ // # for each taxa
+ double xbar_diff = C1[i][0] - C2[i][0];
+ double denom = sqrt(C1[i][2] + C2[i][2]);
+ Ts[i] = abs(xbar_diff/denom); // calculate two sample t-statistic
+ }
+
+ return Ts;
+
+
+ }catch(exception& e) {
+ m->errorOut(e, "MothurMetastats", "permuted_ttests");
+ exit(1);
+ }
+}
+/***********************************************************/
+vector<double> MothurMetastats::calc_qvalues(vector<double>& pValues) {