################################################################################ # R-Code given in Appendix II of # Korner-Nievergelt F, Korner-Nievergelt P, Behr O, Niermann I, Brinkmann R & Hellriegel B (2011) # A new method to determine bird and bat fatality at wind energy turbines from carcass searches. Wildl. Biol. 17: 1-14 #------------------------------------------------------------------------------- # R-Code written and tested for R. 2.12.0 #------------------------------------------------------------------------------- #Estimating the number of fatalities with a credible interval #Step 1: # Describe the uncertainty in the estimates for searcher efficiency and carcass persistence probability by a beta-distribution, i.e. transform the 95 % CI into the shape parameters of the beta-distributions. f <- 0.72; f.lower <- 0.62; f.upper <- 0.81 s <- 0.84; s.lower <- 0.64; s.upper <- 0.94 #-------------------------------------------------------------------- # function to transform the 95% CI into shape parameters of a beta # distribution shapeparameter<-function(m, lwr, upr){ # m = estimate # lwr, upr = lower and upper limit of the 95 % credible or confidence # interval #-------------------------------------------------------------------- ci <- upr - lwr sigma2 <-(ci/4)^2 a <- m*(m*(1-m)/sigma2-1) b <- (1-m)*(m*(1-m)/sigma2-1) list(a=a,b=b) } # end of function shapeparameter-------------------------------------- f.a <- shapeparameter(f, f.lower, f.upper)$a f.b <- shapeparameter(f, f.lower, f.upper)$b s.a <- shapeparameter(s, s.lower, s.upper)$a s.b <- shapeparameter(s, s.lower, s.upper)$b # Step 2 # Define the parameters of the simulations and prepare the vector for the resulting posterior distribution of the number of fatalities. Define a function to obtain the detection probability from searcher efficiency, persistence probability and search interval. Here, we use the new formula presented in this article. Define a function to obtain the posterior distribution of the number of fatalities based on the number of observed carcasses and the detection probability. This formula is based on the theorem of Bayes. Start the loop over step 3 and step 4. maxn <- 500 # define a maximum for the number of fatalities nsim <- 1000 # number of Monte Carlo simulations # prepare a vector for the posterior density distribution # of the estimated number of fatalities: Npostdist <- numeric(maxn+1) #-------------------------------------------------------------- # function to obtain the probability of detecting a carcass # given the searcher efficiency (f), persistence probability # (s), search interval (d) and the total number of searches # (n). pcarcass <- function(s, f, d, n){ # s = probability that a carcass remains 24 hours # f = probability that a carcass is detected by a # searcher during a search given it persisted to the search # d = (average) number of days between two searches # n = number of searches (n * d = length of study period) #-------------------------------------------------------------- x <- (1-f)*s^d A <- s*(1-s^d)/(1-s) summep <- numeric(n) for(k in 0:(n-1)) summep[k+1] <- (n-k)*x^k p <- A*f*sum(summep)/(d*n) p } # end of function pcarcass----------------------------------- #------------------------------------------------------------ # function to obtain the posterior density distribution of the # number of fatalities based on a (known) probability that a # carcass is detected (p) and the number of observed carcasses (nf) posterior.N <- function(p, nf=0, maxN=50, ci.int=0.95, plot=TRUE, dist=FALSE){ # p = probability that a killed animal is detected by a seacher # nf = number of carcasses found # maxN = maximal possible number of fatalities # ci.int = size of the credible interval that should be calculated # plot: posterior distribution is plotted if TRUE # dist: posterior distribution is given if TRUE #--------------------------------------------------------------- N <- nf:maxN if(nf==0) pN <- p*(1-p)^(N-nf) if(nf>0) { denom <- sum(choose(N, nf) * (1-p)^(N-nf)) pN <- choose(N, nf)*(1-p)^(N-nf)/denom pN <- c(rep(0, nf), pN) N <- c(rep(0, nf), N) } if(plot) plot(N, pN, type="h", lwd=5, lend="butt", xlab="Number of fatalities", ylab="Posterior density") index <- cumsum(pN)0) interval <- c(min(N[!indexLower]), min(N[!indexUpper])) if(interval[2]==Inf) cat("Upper limit of CI larger than maxN! -> increase maxN\n") expected <- min(N[!cumsum(pN)<0.5]) results <- list(interval=interval, expected=expected) if(dist==TRUE) results <- list(interval=interval, expected=expected, pN=pN) results } # end of function posterior.N------------------------------------------ for(i in 1:nsim){ #Step 3 #Draw a searcher efficiency f at random from the beta distribution defined by f.a and f.b. fr <- rbeta(1, f.a, f.b) #Draw a persistence probability s at random from the beta distribution defined by s.a and s.b. sr <- rbeta(1, s.a, s.b) #Calculate the detection probability given sr, fr, search interval d=2 and number of searches n=100. pr <- pcarcass(sr, fr, d=2, n=100) #Step 4 #Compute the posterior density distribution of the number of fatalities based on pr and the observed number of carcasses (number found =12) using the function posterior.N. postNtemp <- posterior.N(nf=12, p=pr, maxN=maxn, plot=FALSE, dist=TRUE) #Sum the posterior densities over all nsim simulations. Npostdist <- Npostdist + postNtemp$pN } # close loop i #Step 5 #Scale the summed posterior distribution and extract median and 95 % credible interval. Npostdist.sc <- Npostdist/nsim indexLower <- cumsum(Npostdist.sc) < 0.025 indexMedian <- cumsum(Npostdist.sc) < 0.5 indexUpper <- cumsum(Npostdist.sc) < 0.975 lower <- min(c(0:maxn)[!indexLower]) estimate <- min(c(0:maxn)[!indexMedian]) upper <- min(c(0:maxn)[!indexUpper]) lower; estimate; upper #As a result we receive an estimate of 17 fatalities with a 95 % credible interval of 12-31 fatalities. #The posterior distribution of the number of fatalities is plotted: plot(0:maxn, Npostdist.sc , type="h", lwd=5, lend="butt", xlab="Number of fatalities", ylab="Posterior density", xlim=c(0,50))