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                 double nrare = 0;
  20                 double srare = 0;
  21                 double sabund = 0;
  22
  23                 double Cace, term1, gamace;
  24                 double 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                 double 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                 double j,D_s=0,nn=0,ww=0;
  66                 int Max_Index=rank->getMaxRank()+1;
  67                 double pp, temp1, temp2;
  68                 vector<double> Part_N_Part_F(Max_Index+1,0.0);
  69
  70                 for (j=1; j<Max_Index; j++) if(j<=abund) D_s += rank->get(j);
  71                 for (j=1; j<Max_Index; j++){
  72                         if(j<=abund){
  73                                 nn += rank->get(j) * j;
  74                                 ww += rank->get(j) * j * ( j - 1);
  75                         }
  76                 }
  77                 double C_hat = 1.-rank->get(1)/double(nn);
  78                 double Gamma = ( D_s * ww) / ( C_hat * nn * ( nn - 1.)) - 1.;
  79                 temp1 = double(nn - rank->get(1));
  80                 temp2 = double(nn - 1.);
  81
  82                 if ( Gamma > 0.){
  83                         Part_N_Part_F[1] =  ( D_s + nn) * ( 1. + rank->get(1) * ww / temp1 / temp2) / temp1 + nn * D_s * ww * ( temp1 - 1.) /
  84                         ( temp1 * temp1 * temp2 * temp2) - ( nn + rank->get(1)) / temp1;
  85                         for ( j=2; j<=Max_Index; j++){
  86                                 if(j<=abund){
  87                                         Part_N_Part_F[j] = ( nn * temp1 - j * rank->get(1) * D_s) / temp1 / temp1 * ( 1. + rank->get(1) * ww / temp1 / temp2)
  88                                         + j * rank->get(1) * D_s * nn * ( ( j - 1.) * temp1 * temp2 - ww * ( temp1 + temp2))
  89                                         / temp1 / temp1 / temp1 / temp2 / temp2 + j * rank->get(1) * rank->get(1) / temp1 / temp1;
  90                                 }
  91                         }
  92                 }
  93                 else{
  94                         Part_N_Part_F[1] = ( nn + D_s ) / temp1;
  95                         for ( j=2; j<=Max_Index; j++){
  96                                 if(j<=abund){
  97                                         Part_N_Part_F[j-1] = ( nn * temp1 - j * rank->get(1) * D_s ) / temp1 / temp1;
  98                                 }
  99                         }
 100                 }
 101                 if(Max_Index>abund){
 102                         for ( j=abund+1; j<=Max_Index; j++){
 103                                 Part_N_Part_F[j-1] = 1.;
 104                         }
 105                 }
 106                 for ( temp1=0., temp2=0., j=0; j<Max_Index; j++) {
 107                         pp = Part_N_Part_F[j];
 108                         temp1 += pp * rank->get(j);
 109                         temp2 += pp * pp * rank->get(j);
 110                 }
 111
 112                 double se = temp2 - temp1 * temp1 / ace;
 113
 114                 if(toString(se) == "nan"){
 115                         acelci = ace;
 116                         acehci = ace;
 117                 }
 118                 else if(ace==0.000){
 119                         acelci = ace;
 120                         acehci = ace;
 121                 }
 122                 else if(ace==sobs){
 123                         double ci = 1.96*pow(se,0.5);
 124                         acelci = ace-ci;                                        //ace lci
 125                         acehci = ace+ci;                                        //ace hci
 126                 }else{
 127                         double denom = pow(ace-sobs,2);
 128                         double c = exp(1.96*pow((log(1+se/denom)),0.5));
 129                         acelci = sobs+(ace-sobs)/c;                     //ace lci
 130                         acehci = sobs+(ace-sobs)*c;                     //ace hci
 131                 }
 132
 133                 data[0] = ace;
 134                 data[1] = acelci;
 135                 data[2] = acehci;
 136
 137                 if (isnan(data[0]) || isinf(data[0])) { data[0] = 0; }
 138                 if (isnan(data[1]) || isinf(data[1])) { data[1] = 0; }
 139                 if (isnan(data[2]) || isinf(data[2])) { data[2] = 0; }
 140
 141                 return data;
 142         }
 143         catch(exception& e) {
 144                 m->errorOut(e, "Ace", "getValues");
 145                 exit(1);
 146         }
 147 }
 148
 149 /***********************************************************************/