ace.cpp

   1 /*
   2  *  ace.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 1/7/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "ace.h"
  11
  12 /***********************************************************************/
  13
  14 EstOutput Ace::getValues(SAbundVector* rank) {
  15         try {
  16                 cout << "abund = " << abund << "\n";
  17                 data.resize(3,0);
  18         //      vector<double> aceData(3,0);
  19                 double ace, acelci, acehci;
  20
  21                 int nrare = 0;
  22                 int srare = 0;
  23                 int sabund = 0;
  24
  25                 double Cace, term1, gamace;
  26                 int numsum = 0;
  27
  28                 double maxRank = (double)rank->getMaxRank();
  29
  30                 for(int i=1;i<=maxRank;i++){
  31                         if(i<=abund){
  32                                 srare += rank->get(i);
  33                                 nrare += i*rank->get(i);
  34                                 numsum += (i-1)*i*rank->get(i);
  35                         }
  36                         else if(i>abund)        {sabund += rank->get(i);}
  37                 }
  38                 int sobs = srare + sabund;
  39
  40                 if (nrare == 0){ Cace = 0.0000; }
  41                 else { Cace = 1.0000 -(double)rank->get(1)/(double)nrare; }
  42
  43                 double denom = Cace * (double)(nrare * (nrare-1));
  44
  45                 if(denom <= 0.0){       term1=0.0000;   } else {        term1 = (double)(srare * numsum)/(double)denom - 1.0;   }
  46                 if(term1 >= 0.0){       gamace = term1; } else {        gamace = 0.0;                                                                                   }
  47
  48                 if(gamace >= 0.64){
  49                         gamace = gamace * (1 + (nrare * (1 - Cace) * numsum) / denom);
  50                         if(gamace<0){                   gamace = 0;                     }
  51                 }
  52
  53                 if(Cace == 0.0){
  54                         ace = 0.00;}//ace
  55                 else{
  56                         ace = (double)sabund+((double)srare+(double)rank->get(1)*gamace)/Cace;//ace
  57                 }
  58
  59                 /*
  60                 The following code was obtained from Anne Chao for calculating the SE for her ACE estimator
  61                 My modification was to reset the frequencies so that a singleton is found in rank[1] insted
  62                 of in rank[0], etc.
  63
  64                 I have also added the forumlae to calculate the 95% confidence intervals.
  65                 */
  66
  67                 int j,D_s=0,nn=0,ww=0,Max_Index=rank->getMaxRank()+1;
  68                 double pp, temp1, temp2;
  69                 vector<double> Part_N_Part_F(Max_Index+1,0.0);
  70
  71                 for (j=1; j<Max_Index; j++) if(j<=abund) D_s += rank->get(j);
  72                 for (j=1; j<Max_Index; j++){
  73                         if(j<=abund){
  74                                 nn += rank->get(j) * j;
  75                                 ww += rank->get(j) * j * ( j - 1);
  76                         }
  77                 }
  78                 double C_hat = 1.-rank->get(1)/double(nn);
  79                 double Gamma = ( D_s * ww) / ( C_hat * nn * ( nn - 1.)) - 1.;
  80                 temp1 = double(nn - rank->get(1));
  81                 temp2 = double(nn - 1.);
  82
  83                 if ( Gamma > 0.){
  84                         Part_N_Part_F[1] =  ( D_s + nn) * ( 1. + rank->get(1) * ww / temp1 / temp2) / temp1 + nn * D_s * ww * ( temp1 - 1.) /
  85                         ( temp1 * temp1 * temp2 * temp2) - ( nn + rank->get(1)) / temp1;
  86                         for ( j=2; j<=Max_Index; j++){
  87                                 if(j<=abund){
  88                                         Part_N_Part_F[j] = ( nn * temp1 - j * rank->get(1) * D_s) / temp1 / temp1 * ( 1. + rank->get(1) * ww / temp1 / temp2)
  89                                         + j * rank->get(1) * D_s * nn * ( ( j - 1.) * temp1 * temp2 - ww * ( temp1 + temp2))
  90                                         / temp1 / temp1 / temp1 / temp2 / temp2 + j * rank->get(1) * rank->get(1) / temp1 / temp1;
  91                                 }
  92                         }
  93                 }
  94                 else{
  95                         Part_N_Part_F[1] = ( nn + D_s ) / temp1;
  96                         for ( j=2; j<=Max_Index; j++){
  97                                 if(j<=abund){
  98                                         Part_N_Part_F[j-1] = ( nn * temp1 - j * rank->get(1) * D_s ) / temp1 / temp1;
  99                                 }
 100                         }
 101                 }
 102                 if(Max_Index>abund){
 103                         for ( j=abund+1; j<=Max_Index; j++){
 104                                 Part_N_Part_F[j-1] = 1.;
 105                         }
 106                 }
 107                 for ( temp1=0., temp2=0., j=0; j<Max_Index; j++) {
 108                         pp = Part_N_Part_F[j];
 109                         temp1 += pp * rank->get(j);
 110                         temp2 += pp * pp * rank->get(j);
 111                 }
 112
 113                 double se = temp2 - temp1 * temp1 / ace;
 114
 115                 if(toString(se) == "nan"){
 116                         acelci = ace;
 117                         acehci = ace;
 118                 }
 119                 else if(ace==0.000){
 120                         acelci = ace;
 121                         acehci = ace;
 122                 }
 123                 else if(ace==sobs){
 124                         double ci = 1.96*pow(se,0.5);
 125                         acelci = ace-ci;                                        //ace lci
 126                         acehci = ace+ci;                                        //ace hci
 127                 }else{
 128                         double denom = pow(ace-sobs,2);
 129                         double c = exp(1.96*pow((log(1+se/denom)),0.5));
 130                         acelci = sobs+(ace-sobs)/c;                     //ace lci
 131                         acehci = sobs+(ace-sobs)*c;                     //ace hci
 132                 }
 133
 134                 data[0] = ace;
 135                 data[1] = acelci;
 136                 data[2] = acehci;
 137
 138                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
 139                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
 140                 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
 141
 142                 return data;
 143         }
 144         catch(exception& e) {
 145                 cout << "Standard Error: " << e.what() << " has occurred in the Ace class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 146                 exit(1);
 147         }
 148         catch(...) {
 149                 cout << "An unknown error has occurred in the Ace class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 150                 exit(1);
 151         }
 152 }
 153
 154 /***********************************************************************/