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                 data.resize(3,0);
  17                 double ace, acelci, acehci;
  18
  19                 int nrare = 0;
  20                 int srare = 0;
  21                 int sabund = 0;
  22
  23                 double Cace, term1, gamace;
  24                 int numsum = 0;
  25
  26                 double maxRank = (double)rank->getMaxRank();
  27
  28                 for(int i=1;i<=maxRank;i++){
  29                         if(i<=abund){
  30                                 srare += rank->get(i);
  31                                 nrare += i*rank->get(i);
  32                                 numsum += (i-1)*i*rank->get(i);
  33                         }
  34                         else if(i>abund)        {sabund += rank->get(i);}
  35                 }
  36                 int sobs = srare + sabund;
  37
  38                 if (nrare == 0){ Cace = 0.0000; }
  39                 else { Cace = 1.0000 -(double)rank->get(1)/(double)nrare; }
  40
  41                 double denom = Cace * (double)(nrare * (nrare-1));
  42
  43                 if(denom <= 0.0){       term1=0.0000;   } else {        term1 = (double)(srare * numsum)/(double)denom - 1.0;   }
  44                 if(term1 >= 0.0){       gamace = term1; } else {        gamace = 0.0;                                                                                   }
  45
  46                 if(gamace >= 0.64){
  47                         gamace = gamace * (1 + (nrare * (1 - Cace) * numsum) / denom);
  48                         if(gamace<0){                   gamace = 0;                     }
  49                 }
  50
  51                 if(Cace == 0.0){
  52                         ace = 0.00;}//ace
  53                 else{
  54                         ace = (double)sabund+((double)srare+(double)rank->get(1)*gamace)/Cace;//ace
  55                 }
  56
  57                 /*
  58                 The following code was obtained from Anne Chao for calculating the SE for her ACE estimator
  59                 My modification was to reset the frequencies so that a singleton is found in rank[1] insted
  60                 of in rank[0], etc.
  61
  62                 I have also added the forumlae to calculate the 95% confidence intervals.
  63                 */
  64
  65                 int j,D_s=0,nn=0,ww=0,Max_Index=rank->getMaxRank()+1;
  66                 double pp, temp1, temp2;
  67                 vector<double> Part_N_Part_F(Max_Index+1,0.0);
  68
  69                 for (j=1; j<Max_Index; j++) if(j<=abund) D_s += rank->get(j);
  70                 for (j=1; j<Max_Index; j++){
  71                         if(j<=abund){
  72                                 nn += rank->get(j) * j;
  73                                 ww += rank->get(j) * j * ( j - 1);
  74                         }
  75                 }
  76                 double C_hat = 1.-rank->get(1)/double(nn);
  77                 double Gamma = ( D_s * ww) / ( C_hat * nn * ( nn - 1.)) - 1.;
  78                 temp1 = double(nn - rank->get(1));
  79                 temp2 = double(nn - 1.);
  80
  81                 if ( Gamma > 0.){
  82                         Part_N_Part_F[1] =  ( D_s + nn) * ( 1. + rank->get(1) * ww / temp1 / temp2) / temp1 + nn * D_s * ww * ( temp1 - 1.) /
  83                         ( temp1 * temp1 * temp2 * temp2) - ( nn + rank->get(1)) / temp1;
  84                         for ( j=2; j<=Max_Index; j++){
  85                                 if(j<=abund){
  86                                         Part_N_Part_F[j] = ( nn * temp1 - j * rank->get(1) * D_s) / temp1 / temp1 * ( 1. + rank->get(1) * ww / temp1 / temp2)
  87                                         + j * rank->get(1) * D_s * nn * ( ( j - 1.) * temp1 * temp2 - ww * ( temp1 + temp2))
  88                                         / temp1 / temp1 / temp1 / temp2 / temp2 + j * rank->get(1) * rank->get(1) / temp1 / temp1;
  89                                 }
  90                         }
  91                 }
  92                 else{
  93                         Part_N_Part_F[1] = ( nn + D_s ) / temp1;
  94                         for ( j=2; j<=Max_Index; j++){
  95                                 if(j<=abund){
  96                                         Part_N_Part_F[j-1] = ( nn * temp1 - j * rank->get(1) * D_s ) / temp1 / temp1;
  97                                 }
  98                         }
  99                 }
 100                 if(Max_Index>abund){
 101                         for ( j=abund+1; j<=Max_Index; j++){
 102                                 Part_N_Part_F[j-1] = 1.;
 103                         }
 104                 }
 105                 for ( temp1=0., temp2=0., j=0; j<Max_Index; j++) {
 106                         pp = Part_N_Part_F[j];
 107                         temp1 += pp * rank->get(j);
 108                         temp2 += pp * pp * rank->get(j);
 109                 }
 110
 111                 double se = temp2 - temp1 * temp1 / ace;
 112
 113                 if(toString(se) == "nan"){
 114                         acelci = ace;
 115                         acehci = ace;
 116                 }
 117                 else if(ace==0.000){
 118                         acelci = ace;
 119                         acehci = ace;
 120                 }
 121                 else if(ace==sobs){
 122                         double ci = 1.96*pow(se,0.5);
 123                         acelci = ace-ci;                                        //ace lci
 124                         acehci = ace+ci;                                        //ace hci
 125                 }else{
 126                         double denom = pow(ace-sobs,2);
 127                         double c = exp(1.96*pow((log(1+se/denom)),0.5));
 128                         acelci = sobs+(ace-sobs)/c;                     //ace lci
 129                         acehci = sobs+(ace-sobs)*c;                     //ace hci
 130                 }
 131
 132                 data[0] = ace;
 133                 data[1] = acelci;
 134                 data[2] = acehci;
 135
 136                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
 137                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
 138                 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
 139
 140                 return data;
 141         }
 142         catch(exception& e) {
 143                 m->errorOut(e, "Ace", "getValues");
 144                 exit(1);
 145         }
 146 }
 147
 148 /***********************************************************************/