From 5a7b1d45e9604becce2410333fcc71de818e3782 Mon Sep 17 00:00:00 2001 From: pschloss Date: Tue, 31 May 2011 12:23:23 +0000 Subject: [PATCH] changes to beta diversity calculators --- memchi2.cpp | 12 +++++++++--- mempearson.cpp | 4 ++-- sharedanderbergs.cpp | 2 +- sharedbraycurtis.cpp | 2 +- sharedjabund.cpp | 4 ++-- sharedjclass.cpp | 2 +- sharedjest.cpp | 2 +- sharedkulczynski.cpp | 2 +- sharedkulczynskicody.cpp | 2 +- sharedlennon.cpp | 2 +- sharedmorisitahorn.cpp | 4 ++-- sharedsorabund.cpp | 1 + sharedsorclass.cpp | 2 +- sharedsorest.cpp | 2 +- sharedthetan.cpp | 2 +- sharedthetayc.cpp | 1 + spearman.cpp | 38 +++++++++++++++++++++++++++++++------- structpearson.cpp | 2 +- whittaker.cpp | 2 +- 19 files changed, 60 insertions(+), 28 deletions(-) diff --git a/memchi2.cpp b/memchi2.cpp index 07edd66..595b54f 100644 --- a/memchi2.cpp +++ b/memchi2.cpp @@ -25,7 +25,6 @@ EstOutput MemChi2::getValues(vector shared) { if (shared[1]->getAbundance(i) != 0) { nonZeroB++; } } - double totalTerm = 1 / (float) totalGroups; double sum = 0.0; for (int i = 0; i < shared[0]->getNumBins(); i++) { int A = shared[0]->getAbundance(i); @@ -36,11 +35,18 @@ EstOutput MemChi2::getValues(vector shared) { double Aterm = A / (float) nonZeroA; double Bterm = B / (float) nonZeroB; + + int incidence = 0; + for(int j=0;jgetAbundance(i) != 0){ incidence++; } + } - sum += (totalTerm * ((Aterm-Bterm)*(Aterm-Bterm))); + if(incidence != 0){ + sum += (((Aterm-Bterm)*(Aterm-Bterm))/incidence); + } } - data[0] = sqrt((totalOtus * sum)); + data[0] = sqrt(totalOtus * sum); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/mempearson.cpp b/mempearson.cpp index 934bd01..2503b5d 100644 --- a/mempearson.cpp +++ b/mempearson.cpp @@ -46,9 +46,9 @@ EstOutput MemPearson::getValues(vector shared) { denomTerm2 = sqrt(denomTerm2); double denom = denomTerm1 * denomTerm2; - + if (denom != 0) { - data[0] = 1.0 - (numTerm / denom); + data[0] = (numTerm / denom); }else { data[0] = 1.0; } diff --git a/sharedanderbergs.cpp b/sharedanderbergs.cpp index 09ec692..10dfbd5 100644 --- a/sharedanderbergs.cpp +++ b/sharedanderbergs.cpp @@ -33,7 +33,7 @@ EstOutput Anderberg::getValues(vector shared) { if ((tempA != 0) && (tempB != 0)) { S12++; } } - data[0] = S12 / ((float)((2 * S1) + (2 * S2) - (3 * S12))); + data[0] = 1.0 - S12 / ((float)((2 * S1) + (2 * S2) - (3 * S12))); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/sharedbraycurtis.cpp b/sharedbraycurtis.cpp index 566df3b..0182e13 100644 --- a/sharedbraycurtis.cpp +++ b/sharedbraycurtis.cpp @@ -37,7 +37,7 @@ EstOutput BrayCurtis::getValues(vector shared) { else { sumSharedAB += tempB; } } - data[0] = (2 * sumSharedAB) / (float)( sumSharedA + sumSharedB); + data[0] = 1.0 - (2 * sumSharedAB) / (float)( sumSharedA + sumSharedB); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/sharedjabund.cpp b/sharedjabund.cpp index 79c9f52..d86fcb8 100644 --- a/sharedjabund.cpp +++ b/sharedjabund.cpp @@ -20,8 +20,8 @@ EstOutput JAbund::getValues(vector shared) { UVest = uv->getUVest(shared); //UVest[0] is Uest UVest[1] is Vest - data[0] = (UVest[0] * UVest[1]) / ((float)(UVest[0] + UVest[1] - (UVest[0] * UVest[1]))); - if(data[0] > 1){data[0] = 1; } + data[0] = 1.0-(UVest[0] * UVest[1]) / ((float)(UVest[0] + UVest[1] - (UVest[0] * UVest[1]))); + if(data[0] > 1){data[0] = 0; } if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } return data; diff --git a/sharedjclass.cpp b/sharedjclass.cpp index ede4e31..ac3e94e 100644 --- a/sharedjclass.cpp +++ b/sharedjclass.cpp @@ -34,7 +34,7 @@ EstOutput Jclass::getValues(vector shared) { if ((tempA != 0) && (tempB != 0)) { S12++; } } - data[0] = S12 / (float)(S1 + S2 - S12); + data[0] = 1.0 - S12 / (float)(S1 + S2 - S12); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/sharedjest.cpp b/sharedjest.cpp index 504dd7d..9a1cfbb 100644 --- a/sharedjest.cpp +++ b/sharedjest.cpp @@ -39,7 +39,7 @@ EstOutput Jest::getValues(vector shared) { S1 = chaoS1->getValues(chaoS1Sabund); S2 = chaoS2->getValues(chaoS2Sabund); - data[0] = S12[0] / (float)(S1[0] + S2[0] - S12[0]); + data[0] = 1.0 - S12[0] / (float)(S1[0] + S2[0] - S12[0]); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/sharedkulczynski.cpp b/sharedkulczynski.cpp index c19df53..6aad5df 100644 --- a/sharedkulczynski.cpp +++ b/sharedkulczynski.cpp @@ -33,7 +33,7 @@ EstOutput Kulczynski::getValues(vector shared) { if ((tempA != 0) && (tempB != 0)) { S12++; } } - data[0] = S12 / (float)(S1 + S2 - (2 * S12)); + data[0] = 1.0 - S12 / (float)(S1 + S2 - (2 * S12)); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/sharedkulczynskicody.cpp b/sharedkulczynskicody.cpp index 57fcb24..de90252 100644 --- a/sharedkulczynskicody.cpp +++ b/sharedkulczynskicody.cpp @@ -33,7 +33,7 @@ EstOutput KulczynskiCody::getValues(vector shared) { if ((tempA != 0) && (tempB != 0)) { S12++; } } - data[0] = 0.5 * ((S12 / (float)S1) + (S12 / (float)S2)); + data[0] = 1.0 - 0.5 * ((S12 / (float)S1) + (S12 / (float)S2)); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/sharedlennon.cpp b/sharedlennon.cpp index 0c0f833..f3ed5ec 100644 --- a/sharedlennon.cpp +++ b/sharedlennon.cpp @@ -39,7 +39,7 @@ EstOutput Lennon::getValues(vector shared) { if (tempA < tempB) { min = tempA; } else { min = tempB; } - data[0] = S12 / (float)(S12 + min); + data[0] = 1.0 - S12 / (float)(S12 + min); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/sharedmorisitahorn.cpp b/sharedmorisitahorn.cpp index d52ead5..bd4923f 100644 --- a/sharedmorisitahorn.cpp +++ b/sharedmorisitahorn.cpp @@ -39,9 +39,9 @@ EstOutput MorHorn::getValues(vector shared) { d += relA * relB; } - morhorn = (2 * d) / (a + b); + morhorn = 1- (2 * d) / (a + b); - if (isnan(morhorn) || isinf(morhorn)) { morhorn = 0; } + if (isnan(morhorn) || isinf(morhorn)) { morhorn = 1; } data[0] = morhorn; diff --git a/sharedsorabund.cpp b/sharedsorabund.cpp index 90b825a..9a1bd63 100644 --- a/sharedsorabund.cpp +++ b/sharedsorabund.cpp @@ -23,6 +23,7 @@ EstOutput SorAbund::getValues(vector shared) { data[0] = (2 * UVest[0] * UVest[1]) / ((float)(UVest[0] + UVest[1])); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } + data[0] = 1-data[0]; return data; } diff --git a/sharedsorclass.cpp b/sharedsorclass.cpp index 2cebc7d..32728f5 100644 --- a/sharedsorclass.cpp +++ b/sharedsorclass.cpp @@ -34,7 +34,7 @@ EstOutput SorClass::getValues(vector shared) { if ((tempA != 0) && (tempB != 0)) { S12++; } } - data[0] = (2 * S12) / (float)(S1 + S2); + data[0] = 1.0-(2 * S12) / (float)(S1 + S2); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/sharedsorest.cpp b/sharedsorest.cpp index d3ae59e..b9ab050 100644 --- a/sharedsorest.cpp +++ b/sharedsorest.cpp @@ -38,7 +38,7 @@ EstOutput SorEst::getValues(vector shared) { S1 = chaoS1->getValues(chaoS1Sabund); S2 = chaoS2->getValues(chaoS2Sabund); - data[0] = (2 * S12[0]) / (float)(S1[0] + S2[0]); + data[0] = 1.0-(2 * S12[0]) / (float)(S1[0] + S2[0]); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/sharedthetan.cpp b/sharedthetan.cpp index 361c7b2..7c42a18 100644 --- a/sharedthetan.cpp +++ b/sharedthetan.cpp @@ -46,7 +46,7 @@ EstOutput ThetaN::getValues(vector shared) { if (isnan(thetaN) || isinf(thetaN)) { thetaN = 0; } - data[0] = thetaN; + data[0] = 1.0 - thetaN; return data; } diff --git a/sharedthetayc.cpp b/sharedthetayc.cpp index ed18d92..315a61f 100644 --- a/sharedthetayc.cpp +++ b/sharedthetayc.cpp @@ -74,6 +74,7 @@ EstOutput ThetaYC::getValues(vector shared) { if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; } if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; } + data[0] = 1.0 - data[0]; return data; } catch(exception& e) { diff --git a/spearman.cpp b/spearman.cpp index 72d5667..49d5265 100644 --- a/spearman.cpp +++ b/spearman.cpp @@ -10,6 +10,7 @@ #include "spearman.h" /***********************************************************************/ + EstOutput Spearman::getValues(vector shared) { try { data.resize(1,0); @@ -20,28 +21,51 @@ EstOutput Spearman::getValues(vector shared) { double sumRanks = 0.0; int numOTUS = shared[0]->getNumBins(); - //calc the 2 denominators - for (int i = 0; i < shared[0]->getNumBins(); i++) { + vector rankVectorA(savA.getMaxRank()+1, 0); + int currentRankA = 0; + for(int i=savA.getMaxRank();i>0;i--){ + int numWithAbundanceI = savA.get(i); + + if(numWithAbundanceI > 1){ rankVectorA[i] = (currentRankA + 1 + currentRankA + numWithAbundanceI) / 2.0; } + else { rankVectorA[i] = currentRankA+numWithAbundanceI; } + currentRankA += numWithAbundanceI; + } + rankVectorA[0] = (numOTUS + currentRankA + 1) / 2.0; + + + vector rankVectorB(savB.getMaxRank()+1, 0); + int currentRankB = 0; + for(int i=savB.getMaxRank();i>0;i--){ + int numWithAbundanceI = savB.get(i); + if(numWithAbundanceI > 1){ rankVectorB[i] = (currentRankB + 1 + currentRankB + numWithAbundanceI) / 2.0; } + else { rankVectorB[i] = currentRankB+numWithAbundanceI; } + currentRankB += numWithAbundanceI; + } + rankVectorB[0] = (numOTUS + currentRankB + 1) / 2.0; + + + + for (int i = 0; i < shared[0]->getNumBins(); i++) { int Aij = shared[0]->getAbundance(i); int Bij = shared[1]->getAbundance(i); - int rankA = savA.get(Aij); - int rankB = savB.get(Bij); + float rankA = rankVectorA[Aij]; + float rankB = rankVectorB[Bij]; sumRanks += ((rankA - rankB) * (rankA - rankB)); } - - data[0] = 1.0 - ((6 * sumRanks) / (float) (numOTUS * (numOTUS-1))); + data[0] = 1.0 - ((6 * sumRanks) / (float) (numOTUS * ((numOTUS*numOTUS)-1))); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } return data; } catch(exception& e) { - m->errorOut(e, "Soergel", "getValues"); + m->errorOut(e, "Spearman", "getValues"); exit(1); } } + /***********************************************************************/ diff --git a/structpearson.cpp b/structpearson.cpp index bc2ac76..a21eb99 100644 --- a/structpearson.cpp +++ b/structpearson.cpp @@ -37,7 +37,7 @@ EstOutput StructPearson::getValues(vector shared) { double denom = denomTerm1 * denomTerm2; - data[0] = 1.0 - (numTerm / denom); + data[0] = (numTerm / denom); if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; } diff --git a/whittaker.cpp b/whittaker.cpp index 8688edd..5ab1e5d 100644 --- a/whittaker.cpp +++ b/whittaker.cpp @@ -23,7 +23,7 @@ EstOutput Whittaker::getValues(vector shared){ if(shared[1]->getAbundance(i) != 0){ countB++; } } - data[0] = 2*sTotal/(float)(countA+countB)-1; + data[0] = 2-2*sTotal/(float)(countA+countB); return data; } catch(exception& e) { -- 2.39.2