################################################################################ ################################################################################ # R-Code for the simulations in # 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 #------------------------------------------------------------------------------- # set working directory setwd(".....") #------------------------------------------------------------------------------- # creat functions lower<-function(x) quantile(x, prob=0.025, na.rm=TRUE) # lower limit of 95% interval upper<-function(x) quantile(x, prob=0.975, na.rm=TRUE) # upper limit of 95% interval ## function for the new formula presented in the article det.bat<-function(s, f, d, n){ # s = probability that a dead bat is still there after 24 hours. # f = searchers efficiency, probability that a dead bat is detected by a searcher # d = (average) number of days between two searches # n = number of searches (n * d = study period) x<-(1-f)*s^d A<-s*(1-s^d)/(1-s) summep<-numeric(n+1) for(k in 0:n) summep[k+1]<-(n-k)*x^k p<-A*f*sum(summep)/(d*n) p } ## function for the decreasing-searcher efficiency-version of the new formula presented in the article productsef<-function(f, k, n, j){ res<-1 for(i in 0:(n-j-1)) res<-res*(1-f*k^i) res } det.batdecrease<-function(s, f, d, n, k=0.25){ # s = probability that a dead bat is still there after 24 hours. # f = searchers efficiency, probability that a dead bat is detected by a searcher # d = (average) number of days between two searches # n = number of searches (n * d = study period) # k = factor by which serachers efficiency decreases per search, default = 0.25 as assumed by Manuela Huso A<-s*(1-s^d)/(1-s) summepfound<-A*f for(x in 2:n){ summep<-1+k*s^d*(1-f) for(j in 1:(x-1)) summep<-summep+k^(x-j) * s^((x-j)*d) * productsef(f, k, n=x, j) summepfound<-summepfound+A*f*summep } p<-summepfound/(d*n) p } # fnction presented in Huso (2010) phuso<-function(s, f, d){ # s = probability that a dead bat is still there after 24 hours. # f = searchers efficiency, probability that a dead bat that is present is detected by a searcher # d = (average) number of days between two searches t.bar<-1/(-log(s)) # expected number of days until dissapearance (for exponential distribution of laying time) d.tilde<- -log(0.01)*t.bar d.star <- min(d,d.tilde) r<-t.bar * (1-exp(-d.star/t.bar))/d.star k<-min(1, d.tilde/d) phuso<-r*f*k phuso } #------------------------------------------------------------------------------- ################################################################################ # 1. simulations: E(m) constant f.vector <-c(0.8) # searcher efficiency m.vector <- c(0.1) # average number of fatalities per night d.vector<-c(1, 7, 14) # search interval I.vector <- c(100) # study period (in days) T.vector <- c(3, 30) # expected persistence time of the carcasses results <- expand.grid(f=f.vector, m=m.vector, T=T.vector, I=I.vector, d=d.vector) results$s <- exp(1/-results$T) nsim <- 10000 for(i in 1:nrow(results)){ d <- results$d[i] f <- results$f[i] start.t <- d # day of first search I<-results$I[i] n<-floor(I/d) R<-start.t+I-d s<-results$s[i] # persistence probability # mean persistence time t_quer<-1/(-log(s)) countsvec <- numeric(nsim) diff0<-numeric(nsim) diff1<-numeric(nsim) diff2<-numeric(nsim) diff3<-numeric(nsim) diff4<-numeric(nsim) muN <- results$m[i] for(b in 1:nsim){ n.neu_wahr<-rpois(R, muN) # new fatalities per night n_wahr<-numeric(R) n_wahr[1]<-rbinom(1, n.neu_wahr[1],s) # number of carcasses present at day 1 count<-rep(NA, R) sampling<-rep(FALSE, R) # identify days with searches index<-seq(start.t, start.t+I-d, by=d) sampling[index]<-TRUE sum(sampling) # control: should equal n if(start.t>1) for(k in 2:start.t) n_wahr[k]<-rbinom(1, n_wahr[k-1], s)+rbinom(1, n.neu_wahr[k],s) # number of carcasses present on the following days (without searches) count[start.t]<-rbinom(1, n_wahr[start.t], prob=f) # first search n_wahr[start.t]<-n_wahr[start.t]-count[start.t] for(k in (start.t+1):R){ # true number of carcasses and counted number of carcasses for days with searches n_wahr[k]<-rbinom(1, n_wahr[k-1], s)+rbinom(1, n.neu_wahr[k],s) if(sampling[k]){ count[k]<-rbinom(1, n_wahr[k], prob=f) n_wahr[k]<-n_wahr[k]-count[k] } } N.wahr<-sum(n.neu_wahr[(R-I+1):R]) # true number of fatalities during the study period count<-count[!is.na(count)] N.est0 <- sum(count) sint <- 0 for(o in 1:d){ sint <- sint+s^(d-o+1) } sint <- sint/d N.est1<-sum(count)/(sint*f) # naive estimates # new estimate p2<-det.bat(s, f, d, n) N.est2<- sum(count)/p2 # Erickson 2004 N.est3<-sum(count)/(t_quer*f/d*((exp(d/t_quer)-1)/(exp(d/t_quer)-1+f))) # Erickson 2004 # Huso 2010 p4<-phuso(s,f,d) N.est4<-sum(count)/p4 countsvec[b] <- sum(count, na.rm=TRUE) diff0[b]<-ifelse(N.wahr>0, (N.est0-N.wahr)/N.wahr, NA) diff1[b]<-ifelse(N.wahr>0, (N.est1-N.wahr)/N.wahr, NA) diff2[b]<-ifelse(N.wahr>0, (N.est2-N.wahr)/N.wahr, NA) diff3[b]<-ifelse(N.wahr>0, (N.est3-N.wahr)/N.wahr, NA) diff4[b]<-ifelse(N.wahr>0, (N.est4-N.wahr)/N.wahr, NA) } # close b results$antnullcount[i] <- sum(countsvec==0)/nsim results$mcount[i] <- mean(countsvec) results$biascount[i] <- mean(diff0, na.rm=TRUE) results$ugrcount[i] <- lower(diff0) results$ogrcount[i] <- upper(diff0) results$sdcount[i] <- sd(diff0, na.rm=TRUE) results$biaseinf[i] <- mean(diff1, na.rm=TRUE) results$ugreinf[i] <- lower(diff1) results$ogreinf[i] <- upper(diff1) results$sdeinf[i] <- sd(diff1, na.rm=TRUE) results$biasoikostat[i] <- mean(diff2, na.rm=TRUE) results$ugroikostat[i] <- lower(diff2) results$ogroikostat[i] <- upper(diff2) results$sdoikostat[i] <- sd(diff2, na.rm=TRUE) results$biaserick[i] <- mean(diff3, na.rm=TRUE) results$ugrerick[i] <- lower(diff3) results$ogrerick[i] <- upper(diff3) results$sderick[i] <- sd(diff3, na.rm=TRUE) results$biashuso[i] <- mean(diff4, na.rm=TRUE) results$ugrhuso[i] <- lower(diff4) results$ogrhuso[i] <- upper(diff4) results$sdhuso[i] <- sd(diff4, na.rm=TRUE) write.table(results, file="simresults110404mkonst.txt", row.names=FALSE) # save results } # close i #------------------------------------------------------------------------------- # 2. simulation: E(m) proportional to emprical acoustic activity measurements f.vector <-c(0.8) # searcher efficiency m.vector <- c(0.1) # average number of fatalities per night T.vector <- c(3, 30) d.vector<-c(1, 7, 14) # search interval I.vector <- c(100) # study period in days results <- expand.grid(f=f.vector, m=m.vector, T=T.vector, I=I.vector, d=d.vector) results$s <- exp(1/-results$T) nsim <- 10000 # original data on bat activity (see Behr et al. 2011) activity <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 5, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 36, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 4, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0, 35, 13, 1, 0, 6, 0, 25, 16, 5, 0, 0, 3, 2, 3, 16, 8, 379, 10, 16, 68, 7, 1, 161, 23, 39, 3, 16, 20, 0, 15, 51, 1, 1, 0, 29, 0, 0, 23, 0, 0, 0, 4, 1, 0, 1, 8, 64, 1, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 4, 8, 0, 0, 0, 21, 1, 19, 0, 0, 1, 1, 16, 0, 0, 118, 0, 0, 224, 0, 8, 167, 2, 11, 0, 27, 23, 0, 26, 3, 3, 3, 1, 7, 10, 0, 6, 1, 0, 0, 0, 0, 3, 1, 0, 44, 0, 0, 0, 6, 0, 0, 0, 0, 5, 0, 0, 2, 0, 0, 0, 0, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 4, 2, 1, 0, 12, 0, 4, 6, 2, 0, 0, 0, 24, 59, 58, 128, 28, 34, 41, 36, 146, 69, 8, 0, 1, 9, 0, 0, 1, 363, 55, 37, 0, 31, 0, 0, 3, 4, 41, 4, 136, 23, 4, 2, 16, 60, 66, 4, 24, 5, 10, 7, 27, 15, 11, 12, 15, 110, 1, 14, 6, 0, 0, 0, 0, 11, 0, 1, 0, 1, 0, 7, 0, 0, 2, 0, 0, 8, 5, 0, 0, 0, 2, 1, 0, 0, 1, 0, 0, 0, 16, 18, 46, 10, 1, 6, 8, 4, 76, 69, 3, 1, 0, 7, 0, 1, 6, 87, 12, 46, 3, 14, 0, 0, 3, 0, 30, 4, 59, 22, 0, 0, 4, 35, 7, 1, 30, 1, 2, 13, 2, 5, 4, 5, 3, 30, 0, 5, 1, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 6, 0, 0, 5, 4, 0, 0, 0, 5, 0, 0, 0, 0, 0, 168, 182, 59, 15, 64, 1, 50, 244, 8, 0, 0, 0, 11, 0, 0, 1, 484, 37, 10, 0, 44, 0, 0, 17, 8, 113, 5, 28, 22, 0, 0, 72, 54, 118, 8, 13, 5, 84, 0, 35, 8, 48, 38, 29, 36, 22, 24, 14, 15, 0, 0, 0, 0, 0, 6, 15, 0, 0, 0, 0, 5, 3, 2, 0, 5, 2, 9, 0, 111, 185, 99, 123, 30, 20, 209, 64, 5, 1, 0, 4, 50, 0, 0, 0, 278, 24, 13, 0, 14, 0, 0, 0, 26, 87, 19, 164, 26, 12, 1, 27, 37, 118, 10, 6, 17, 28, 6, 17, 13, 8, 1, 53, 35, 21, 39, 37, 1, 0, 1, 9, 0, 0, 0, 5, 0, 0, 1, 0, 12, 0, 5, 2, 6, 1, 0, 0, 1, 0, 0, 0, 5, 0, 0, 70, 216, 6, 0, 1, 229, 0, 6, 0, 0, 0, 0, 0, 3, 8, 0, 43, 0, 56, 77, 0, 3, 0, 0, 0, 0, 30, 7, 1, 0, 0, 0, 1, 0, 28, 259, 13, 2, 0, 0, 0, 2, 0, 0, 7, 0, 79, 0, 7, 19, 0, 8, 0, 0, 0, 7, 0, 0, 4, 7, 16, 1, 33, 9, 0, 0, 0, 3, 0, 5, 12, 0, 6, 8, 1, 2, 0, 3, 13, 13, 0, 0, 2, 0, 74, 55, 22, 0, 0, 10, 4, 18, 2, 0, 0, 0, 1, 25, 0, 5, 19, 21, 0, 0, 37, 11, 69, 11, 82, 36, 38, 19, 44, 164, 173, 44, 50, 35, 13, 0, 44, 63, 0, 0, 0, 0, 0, 4, 5, 12, 0, 0, 0, 0, 0, 0, 1, 0, 2, 14, 0, 0, 0, 0, 0, 0, 0, 26, 50, 18, 62, 1, 5, 0, 67, 134, 13, 0, 0, 41, 41, 0, 5, 14, 3, 0, 2, 40, 0, 38, 11, 12, 0, 0, 51, 0, 16, 0, 139, 26, 45, 15, 91, 19, 28, 111, 44, 97, 132, 46, 295, 12, 58, 374, 348, 157, 130, 226, 29, 9, 87, 41, 0, 8, 22, 1, 3, 2, 4, 60, 5, 3, 0, 6, 0, 0, 0, 0, 0, 9, 7, 0, 0, 2, 0, 9, 5, 0, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 1, 10, 0, 0, 9, 0, 0, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 11, 3, 0, 0, 0, 0, 2, 0, 0, 0, 10, 14, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 5, 3, 5, 0, 22, 0, 43, 0, 0, 11, 0, 0, 0, 0, 0, 0, 0, 0, 14, 0, 0, 0, 5, 3, 0, 0, 0, 0, 0, 0, 8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 39, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 21, 0, 0, 1, 0, 5, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 32, 0, 0, 9, 0, 9, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 22, 0, 0, 1, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 23, 3, 3, 0, 0, 2, 0, 3, 1, 10, 1, 3, 3, 2, 1, 0, 10, 36, 26, 3, 90, 19, 13, 66, 5, 3, 34, 57, 303, 5, 23, 46, 55, 46, 5, 42, 36, 1, 0, 2, 1, 0, 0, 0, 10, 28, 7, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 1, 2, 4, 2, 0, 16, 0, 0, 0, 0, 0, 0, 5, 0, 1, 0, 0, 2, 13, 5, 13, 9, 1, 9, 0, 1, 0, 3, 0, 3, 3, 0, 0, 0, 1, 0, 0, 0, 2, 2, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) N <- length(activity) N*m.vector activity.scaled <- activity/sum(activity) *N*m.vector # scale activity so that the mean is m for(i in 1:nrow(results)){ d <- results$d[i] f <- results$f[i] start.t <- d # day if first search I<-results$I[i] n<-floor(I/d) R<-start.t+I-d s<-results$s[i] # persistence probability # average persistence time t_quer<-1/(-log(s)) countsvec <- numeric(nsim) diff0<-numeric(nsim) diff1<-numeric(nsim) diff2<-numeric(nsim) diff3<-numeric(nsim) diff4<-numeric(nsim) for(b in 1:nsim){ starttime <- sample(1:(N-R), 1) # choose at random a start point of the interval used n.neu_wahr<-rpois(R, activity.scaled[starttime:(starttime+R)]) # new fatalities per day n_wahr<-numeric(R) n_wahr[1]<-rbinom(1, n.neu_wahr[1],s) # number of carcasses at day 1 count<-rep(NA, R) sampling<-rep(FALSE, R) # identify days with searches index<-seq(start.t, start.t+I-d, by=d) sampling[index]<-TRUE sum(sampling) # control: should equal n if(start.t>1) for(k in 2:start.t) n_wahr[k]<-rbinom(1, n_wahr[k-1], s)+rbinom(1, n.neu_wahr[k],s) # number of carcasses for the following days without searches count[start.t]<-rbinom(1, n_wahr[start.t], prob=f) # first search n_wahr[start.t]<-n_wahr[start.t]-count[start.t] for(k in (start.t+1):R){ # true number of carcasses and counts for days with searches n_wahr[k]<-rbinom(1, n_wahr[k-1], s)+rbinom(1, n.neu_wahr[k],s) if(sampling[k]){ count[k]<-rbinom(1, n_wahr[k], prob=f) n_wahr[k]<-n_wahr[k]-count[k] } } N.wahr<-sum(n.neu_wahr[(R-I+1):R]) # true number of fatalities during the study period count<-count[!is.na(count)] N.est0 <- sum(count) sint <- 0 for(o in 1:d){ sint <- sint+s^(d-o+1) } sint <- sint/d N.est1<-sum(count)/(sint*f) # naive estimate # new method presented in the article p2<-det.bat(s, f, d, n) N.est2<- sum(count)/p2 # Erickson 2004 N.est3<-sum(count)/(t_quer*f/d*((exp(d/t_quer)-1)/(exp(d/t_quer)-1+f))) # Erickson 2004 # Huso 2010 p4<-phuso(s,f,d) N.est4<-sum(count)/p4 countsvec[b] <- sum(count, na.rm=TRUE) diff0[b]<-ifelse(N.wahr>0, (N.est0-N.wahr)/N.wahr, NA) diff1[b]<-ifelse(N.wahr>0, (N.est1-N.wahr)/N.wahr, NA) diff2[b]<-ifelse(N.wahr>0, (N.est2-N.wahr)/N.wahr, NA) diff3[b]<-ifelse(N.wahr>0, (N.est3-N.wahr)/N.wahr, NA) diff4[b]<-ifelse(N.wahr>0, (N.est4-N.wahr)/N.wahr, NA) } # close b results$antnullcount[i] <- sum(countsvec==0)/nsim results$mcount[i] <- mean(countsvec) results$biascount[i] <- mean(diff0, na.rm=TRUE) results$ugrcount[i] <- lower(diff0) results$ogrcount[i] <- upper(diff0) results$sdcount[i] <- sd(diff0, na.rm=TRUE) results$biaseinf[i] <- mean(diff1, na.rm=TRUE) results$ugreinf[i] <- lower(diff1) results$ogreinf[i] <- upper(diff1) results$sdeinf[i] <- sd(diff1, na.rm=TRUE) results$biasoikostat[i] <- mean(diff2, na.rm=TRUE) results$ugroikostat[i] <- lower(diff2) results$ogroikostat[i] <- upper(diff2) results$sdoikostat[i] <- sd(diff2, na.rm=TRUE) results$biaserick[i] <- mean(diff3, na.rm=TRUE) results$ugrerick[i] <- lower(diff3) results$ogrerick[i] <- upper(diff3) results$sderick[i] <- sd(diff3, na.rm=TRUE) results$biashuso[i] <- mean(diff4, na.rm=TRUE) results$ugrhuso[i] <- lower(diff4) results$ogrhuso[i] <- upper(diff4) results$sdhuso[i] <- sd(diff4, na.rm=TRUE) write.table(results, file="simresults110404mpropA.txt", row.names=FALSE) } # close i ################################################################################ # 3. simulation: E(m) constant, removal rate increasing with time library(survival) f.vector <-c(0.8) # searcher efficiency m.vector <- c(0.1) # average number of fatalities per night T.vector <- c(3, 30) # mean persistence time d.vector<-c(1, 7, 14) # search interval I.vector <- c(100) # study period in days alpha.vector <- c(0.7, 1.3) #alpha >1 -> increasing removal, alpha <1-> decreasing removal results <- expand.grid(f=f.vector, m=m.vector, T=T.vector, I=I.vector, d=d.vector, alpha=alpha.vector) results$lambda <- gamma(1+1/results$alpha)/results$T # show the simulation scenarios Fig S1 par(mfrow=c(2,2), mar=c(4,6,3,1)) x <- seq(0, 14, by=0.1) ntest <- 1000 index <- results$d==1 for(i in 1:4){ t_quer <- results$T[index][i] alpha <- results$alpha[index][i] lambda <- gamma(1+1/alpha)/t_quer sx <- exp(-(lambda*(x+1))^alpha)/exp(-(lambda*(x))^alpha) Sx <- exp(-(lambda*x)^alpha) testt <- rweibull(ntest, shape=alpha, scale=1/lambda) status <- rep(1, ntest) status[testt>14] <- 0 testt[testt>14] <- 14 testt <- ceiling(testt) survmod <- survreg(Surv(testt, status)~1, dist="exponential") results$s[results$T==t_quer & results$alpha==alpha] <- exp(1/(-predict(survmod)[1])) Sxass <- results$s[index][i]^x plot(x, 1-sx, type="l", lty=3, lwd=2, ylim=c(0,1), las=1, ylab="PROPORTION REMAINING\nREMOVAL PROBABILITY", main=paste("Mean persistence time = ", t_quer, ",\n", ifelse(alpha<1, "decreasing", "increasing"), " removal rate", sep=""), xlab="", cex.lab=1.3, font.main=1, bty="l") abline(h=1-results$s[index][i], lwd=2, lty=3, col=grey(0.5)) lines(x,Sxass, lty=1, lwd=2, col=grey(0.5)) lines(x, Sx, lty=1, lwd=2) mtext("DAY", side=1, line=2) } legend(0, 0.47, lwd=2, lty=c(1,1, 3, 3),col=c(1,grey(0.5), 1,grey(0.5)), cex=0.8, legend=c("True proportion remaining", "Assumed proportion remaining", "True removal probability", "Assumed removal probability")) # end of Fig S1-------------------------------------------------------------------- nsim <- 1000 for(i in 1:nrow(results)){ d <- results$d[i] f <- results$f[i] start.t <- d # day of first search I<-results$I[i] n<-floor(I/d) R<-start.t+I-d t_quer<-results$T[i] s <- results$s[i] countsvec <- numeric(nsim) diff0<-numeric(nsim) diff1<-numeric(nsim) diff2<-numeric(nsim) diff3<-numeric(nsim) diff4<-numeric(nsim) muN <- results$m[i] t.time <- 1:(R+1) s.time <- exp(-(results$lambda[i]*(t.time+1))^results$alpha[i])/exp(-(results$lambda[i]*(t.time))^results$alpha[i]) for(b in 1:nsim){ n.neu_wahr<-rpois(R, muN) # freshly dead bats per night count<-rep(NA, R) sampling<-rep(FALSE, R+1) # identify days with searches index<-seq(start.t, start.t+I-d, by=d) sampling[index]<-TRUE sum(sampling) # control: must equal the number of searches n_true.mat <- matrix(ncol=R+1,nrow=R) count.mat <- matrix(ncol=R+1,nrow=R) diag(n_true.mat) <- rbinom(R, n.neu_wahr, prob=s.time[1]) # number of bats still there after the first night for(z in 1:R){ if(sampling[z]){ # first day count.mat[z, z] <- rbinom(1, n_true.mat[z,z], prob=f) n_true.mat[z,z] <- n_true.mat[z,z] - count.mat[z, z] } for(zz in (z+1):(R+1)){ # following days n_true.mat[z,zz] <- rbinom(1, n_true.mat[z, zz-1], prob=s.time[zz-z+1]) if(sampling[zz]){ count.mat[z, zz] <- rbinom(1, n_true.mat[z,zz], prob=f) n_true.mat[z,zz] <- n_true.mat[z,zz] - count.mat[z, zz] } } } count <- apply(count.mat, 2, sum, na.rm=TRUE) N.wahr<-sum(n.neu_wahr[(R-I+1):R]) # true number of fatalities during the study period count<-count[!is.na(count)] N.est0 <- sum(count) sint <- 0 for(o in 1:d){ sint <- sint+s^(d-o+1) } sint <- sint/d N.est1<-sum(count)/(sint*f) # naive estimate # new method presented in the article p2<-det.bat(s, f, d, n) N.est2<- sum(count)/p2 # Erickson 2004 N.est3<-sum(count)/(t_quer*f/d*((exp(d/t_quer)-1)/(exp(d/t_quer)-1+f))) # Erickson 2004 # Huso 2010 p4<-phuso(s,f,d) N.est4<-sum(count)/p4 countsvec[b] <- sum(count, na.rm=TRUE) diff0[b]<-ifelse(N.wahr>0, (N.est0-N.wahr)/N.wahr, NA) diff1[b]<-ifelse(N.wahr>0, (N.est1-N.wahr)/N.wahr, NA) diff2[b]<-ifelse(N.wahr>0, (N.est2-N.wahr)/N.wahr, NA) diff3[b]<-ifelse(N.wahr>0, (N.est3-N.wahr)/N.wahr, NA) diff4[b]<-ifelse(N.wahr>0, (N.est4-N.wahr)/N.wahr, NA) } # close b results$antnullcount[i] <- sum(countsvec==0)/nsim results$mcount[i] <- mean(countsvec) results$biascount[i] <- mean(diff0, na.rm=TRUE) results$ugrcount[i] <- lower(diff0) results$ogrcount[i] <- upper(diff0) results$sdcount[i] <- sd(diff0, na.rm=TRUE) results$biaseinf[i] <- mean(diff1, na.rm=TRUE) results$ugreinf[i] <- lower(diff1) results$ogreinf[i] <- upper(diff1) results$sdeinf[i] <- sd(diff1, na.rm=TRUE) results$biasoikostat[i] <- mean(diff2, na.rm=TRUE) results$ugroikostat[i] <- lower(diff2) results$ogroikostat[i] <- upper(diff2) results$sdoikostat[i] <- sd(diff2, na.rm=TRUE) results$biaserick[i] <- mean(diff3, na.rm=TRUE) results$ugrerick[i] <- lower(diff3) results$ogrerick[i] <- upper(diff3) results$sderick[i] <- sd(diff3, na.rm=TRUE) results$biashuso[i] <- mean(diff4, na.rm=TRUE) results$ugrhuso[i] <- lower(diff4) results$ogrhuso[i] <- upper(diff4) results$sdhuso[i] <- sd(diff4, na.rm=TRUE) write.table(results, file="simresults110416svariable.txt", row.names=FALSE) } # close i #------------------------------------------------------------------------------- # Draw figures ################################################################################ results.konstm <- read.table("simresults110404mkonst.txt", header=TRUE) results.varm <- read.table("simresults110404mpropA.txt", header=TRUE) results.vars <- read.table("simresults110416svariable.txt", header=TRUE) # Grafikparameter: versch <- c(-0.4, -0.2, 0, 0.2, 0.4) symbole <- c(20, 21, 22, 23, 24) oblim <- 3 cex.main <- 1.2 col.oikostat <- grey(0.6) cex.axis <- 1.4 cex.lab <- 1.7 windows(width=14, height=8) par(mar=c(0.1, 0.5, 0.5, 0.1), mgp=c(1, 0.3, 0), oma=c(6, 8, 4, 6), mfrow=c(2,4)) index <- results.konstm$T==3 plot(1:3, results.konstm$antnullcount[index], type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8)) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.konstm$antnullcount[index][i], border=NA, col=grey(0.8)) axis(4, las=1, labels=NA, tck=0.015) par(new=TRUE) plot(c(1:3)+versch[1], results.konstm$biascount[index], ylim=c(-1, oblim), las=1, tck=0.015, cex.lab=1.2, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.konstm$ugrcount[index], c(1:3)+versch[1], results.konstm$ogrcount[index], lwd=3, lend="butt") points(c(1:3)+versch[1], results.konstm$biascount[index], pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.konstm$ugreinf[index], c(1:3)+versch[2], results.konstm$ogreinf[index], lwd=3, lend="butt") points(c(1:3)+versch[2], results.konstm$biasein[index], pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.konstm$ugroikostat[index], c(1:3)+versch[3], results.konstm$ogroikostat[index], lwd=3, lend="butt") points(c(1:3)+versch[3], results.konstm$biasoikostat[index], pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.konstm$ugrerick[index], c(1:3)+versch[4], results.konstm$ogrerick[index], lwd=3, lend="butt") points(c(1:3)+versch[4], results.konstm$biaserick[index], pch=symbole[4], bg="white", cex=2) segments(c(1:3)+versch[5], results.konstm$ugrhuso[index], c(1:3)+versch[5], results.konstm$ogrhuso[index], lwd=3, lend="butt") points(c(1:3)+versch[5], results.konstm$biashuso[index], pch=symbole[5], bg="white", cex=2) axis(2, at=seq(-1, oblim, by=0.5), las=1, cex.axis=cex.axis*0.8, tck=0.015) box() mtext(side=2, expression(bar(t) == 3), line=2, cex=cex.lab) mtext("Mortality and\nremoval constant", side=3, line=1, cex=cex.main) par(mfg=c(2,1)) index <- results.konstm$T==30 plot(1:3, results.konstm$antnullcount[index], type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8)) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.konstm$antnullcount[index][i], border=NA, col=grey(0.8)) axis(4, las=1, labels=NA, tck=0.015) par(new=TRUE) plot(c(1:3)+versch[1], results.konstm$biascount[index], ylim=c(-1, oblim), las=1, tck=0.015, cex.lab=1.2, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.konstm$ugrcount[index], c(1:3)+versch[1], results.konstm$ogrcount[index], lwd=3, lend="butt") points(c(1:3)+versch[1], results.konstm$biascount[index], pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.konstm$ugreinf[index], c(1:3)+versch[2], results.konstm$ogreinf[index], lwd=3, lend="butt") points(c(1:3)+versch[2], results.konstm$biasein[index], pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.konstm$ugroikostat[index], c(1:3)+versch[3], results.konstm$ogroikostat[index], lwd=3, lend="butt") points(c(1:3)+versch[3], results.konstm$biasoikostat[index], pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.konstm$ugrerick[index], c(1:3)+versch[4], results.konstm$ogrerick[index], lwd=3, lend="butt") points(c(1:3)+versch[4], results.konstm$biaserick[index], pch=symbole[4], bg="white", cex=2) segments(c(1:3)+versch[5], results.konstm$ugrhuso[index], c(1:3)+versch[5], results.konstm$ogrhuso[index], lwd=3, lend="butt") points(c(1:3)+versch[5], results.konstm$biashuso[index], pch=symbole[5], bg="white", cex=2) axis(2, at=seq(-1, oblim, by=0.5), las=1, cex.axis=cex.axis*0.8, tck=0.015, labels=c(seq(-1, oblim, by=0.5)[-length(seq(-1, oblim, by=0.5))], NA)) axis(1, at=1:3, labels=results.konstm$d[index], tck=0, cex.axis=cex.axis) box() legend(0.5, 2.9, pch=c(20:24), pt.bg=c("white", "white", col.oikostat, "white", "white"), legend=c("Count", "Simple formula", "This article", "Erickson et al. (2004)", "Huso (2010)"), cex=2) mtext(side=2, expression(bar(t) == 30), line=2.5, cex=cex.lab) # m prop Activity par(mfg=c(1,2)) index <- results.varm$T==3 plot(1:3, results.varm$antnullcount[index], type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8)) mtext("Mortality variable,\nremoval const.", side=3, line=1, cex=cex.main) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.varm$antnullcount[index][i], border=NA, col=grey(0.8)) axis(4, las=1, labels=NA, tck=0.015) par(new=TRUE) plot(c(1:3)+versch[1], results.varm$biascount[index], ylim=c(-1, oblim), las=1, tck=0.015, cex.lab=1.2, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.varm$ugrcount[index], c(1:3)+versch[1], results.varm$ogrcount[index], lwd=3, lend="butt") points(c(1:3)+versch[1], results.varm$biascount[index], pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.varm$ugreinf[index], c(1:3)+versch[2], results.varm$ogreinf[index], lwd=3, lend="butt") points(c(1:3)+versch[2], results.varm$biasein[index], pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.varm$ugroikostat[index], c(1:3)+versch[3], results.varm$ogroikostat[index], lwd=3, lend="butt") points(c(1:3)+versch[3], results.varm$biasoikostat[index], pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.varm$ugrerick[index], c(1:3)+versch[4], results.varm$ogrerick[index], lwd=3, lend="butt") points(c(1:3)+versch[4], results.varm$biaserick[index], pch=symbole[4], bg="white", cex=2) segments(c(1:3)+versch[5], results.varm$ugrhuso[index], c(1:3)+versch[5], results.varm$ogrhuso[index], lwd=3, lend="butt") points(c(1:3)+versch[5], results.varm$biashuso[index], pch=symbole[5], bg="white", cex=2) axis(2, at=seq(-1, oblim, by=0.5), las=1, labels=NA, cex.axis=cex.axis, tck=0.015) box() par(mfg=c(2,2)) index <- results.varm$T==30 plot(1:3, results.varm$antnullcount[index], type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8)) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.varm$antnullcount[index][i], border=NA, col=grey(0.8)) axis(4, las=1, labels=NA, tck=0.015) par(new=TRUE) plot(c(1:3)+versch[1], results.varm$biascount[index], ylim=c(-1, oblim), las=1, tck=0.015, cex.axis=cex.axis, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.varm$ugrcount[index], c(1:3)+versch[1], results.varm$ogrcount[index], lwd=3, lend="butt") points(c(1:3)+versch[1], results.varm$biascount[index], pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.varm$ugreinf[index], c(1:3)+versch[2], results.varm$ogreinf[index], lwd=3, lend="butt") points(c(1:3)+versch[2], results.varm$biasein[index], pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.varm$ugroikostat[index], c(1:3)+versch[3], results.varm$ogroikostat[index], lwd=3, lend="butt") points(c(1:3)+versch[3], results.varm$biasoikostat[index], pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.varm$ugrerick[index], c(1:3)+versch[4], results.varm$ogrerick[index], lwd=3, lend="butt") points(c(1:3)+versch[4], results.varm$biaserick[index], pch=symbole[4], bg="white", cex=2) segments(c(1:3)+versch[5], results.varm$ugrhuso[index], c(1:3)+versch[5], results.varm$ogrhuso[index], lwd=3, lend="butt") points(c(1:3)+versch[5], results.varm$biashuso[index], pch=symbole[5], bg="white", cex=2) axis(2, at=seq(-1, oblim, by=0.5),labels=NA, las=1, cex.axis=cex.axis, tck=0.015) axis(1, at=1:3, labels=results.varm$d[index], tck=0, cex.axis=cex.axis*1.2) box() # decreasing removal par(mfg=c(1,3)) index <- results.vars$T==3 & results.vars$alpha<1 # decreasing removal plot(1:3, results.vars$antnullcount[index], type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8)) mtext("Mortality constant,\nremoval decreasing", side=3, line=1, cex=cex.main) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.vars$antnullcount[index][i], border=NA, col=grey(0.8)) axis(4, las=1, labels=NA, tck=0.015) par(new=TRUE) plot(c(1:3)+versch[1], results.vars$biascount[index], ylim=c(-1, oblim), las=1, tck=0.015, cex.lab=1.2, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.vars$ugrcount[index], c(1:3)+versch[1], results.vars$ogrcount[index], lwd=3, lend="butt") points(c(1:3)+versch[1], results.vars$biascount[index], pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.vars$ugreinf[index], c(1:3)+versch[2], results.vars$ogreinf[index], lwd=3, lend="butt") points(c(1:3)+versch[2], results.vars$biasein[index], pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.vars$ugroikostat[index], c(1:3)+versch[3], results.vars$ogroikostat[index], lwd=3, lend="butt") points(c(1:3)+versch[3], results.vars$biasoikostat[index], pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.vars$ugrerick[index], c(1:3)+versch[4], results.vars$ogrerick[index], lwd=3, lend="butt") points(c(1:3)+versch[4], results.vars$biaserick[index], pch=symbole[4], bg="white", cex=2) segments(c(1:3)+versch[5], results.vars$ugrhuso[index], c(1:3)+versch[5], results.vars$ogrhuso[index], lwd=3, lend="butt") points(c(1:3)+versch[5], results.vars$biashuso[index], pch=symbole[5], bg="white", cex=2) axis(2, at=seq(-1, oblim, by=0.5), las=1, labels=NA, cex.axis=cex.axis, tck=0.015) box() par(mfg=c(2,3)) index <- results.vars$T==30 & results.vars$alpha<1 # decreasing removal plot(1:3, results.vars$antnullcount[index], type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8)) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.vars$antnullcount[index][i], border=NA, col=grey(0.8)) axis(4, las=1, labels=NA, tck=0.015) par(new=TRUE) plot(c(1:3)+versch[1], results.vars$biascount[index], ylim=c(-1, oblim), las=1, tck=0.015, cex.axis=cex.axis, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.vars$ugrcount[index], c(1:3)+versch[1], results.vars$ogrcount[index], lwd=3, lend="butt") points(c(1:3)+versch[1], results.vars$biascount[index], pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.vars$ugreinf[index], c(1:3)+versch[2], results.vars$ogreinf[index], lwd=3, lend="butt") points(c(1:3)+versch[2], results.vars$biasein[index], pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.vars$ugroikostat[index], c(1:3)+versch[3], results.vars$ogroikostat[index], lwd=3, lend="butt") points(c(1:3)+versch[3], results.vars$biasoikostat[index], pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.vars$ugrerick[index], c(1:3)+versch[4], results.vars$ogrerick[index], lwd=3, lend="butt") points(c(1:3)+versch[4], results.vars$biaserick[index], pch=symbole[4], bg="white", cex=2) segments(c(1:3)+versch[5], results.vars$ugrhuso[index], c(1:3)+versch[5], results.vars$ogrhuso[index], lwd=3, lend="butt") points(c(1:3)+versch[5], results.vars$biashuso[index], pch=symbole[5], bg="white", cex=2) axis(2, at=seq(-1, oblim, by=0.5),labels=NA, las=1, cex.axis=cex.axis, tck=0.015) axis(1, at=1:3, labels=results.vars$d[index], tck=0, cex.axis=cex.axis*1.2) box() # increasing removal par(mfg=c(1,4)) index <- results.vars$T==3 & results.vars$alpha>1 # decreasing removal plot(1:3, results.vars$antnullcount[index], type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8)) mtext("Mortality constant,\nremoval increasing", side=3, line=1, cex=cex.main) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.vars$antnullcount[index][i], border=NA, col=grey(0.8)) axis(4, las=1, tck=0.015, cex.axis=cex.axis) par(new=TRUE) plot(c(1:3)+versch[1], results.vars$biascount[index], ylim=c(-1, oblim), las=1, tck=0.015, cex.lab=1.2, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.vars$ugrcount[index], c(1:3)+versch[1], results.vars$ogrcount[index], lwd=3, lend="butt") points(c(1:3)+versch[1], results.vars$biascount[index], pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.vars$ugreinf[index], c(1:3)+versch[2], results.vars$ogreinf[index], lwd=3, lend="butt") points(c(1:3)+versch[2], results.vars$biasein[index], pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.vars$ugroikostat[index], c(1:3)+versch[3], results.vars$ogroikostat[index], lwd=3, lend="butt") points(c(1:3)+versch[3], results.vars$biasoikostat[index], pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.vars$ugrerick[index], c(1:3)+versch[4], results.vars$ogrerick[index], lwd=3, lend="butt") points(c(1:3)+versch[4], results.vars$biaserick[index], pch=symbole[4], bg="white", cex=2) segments(c(1:3)+versch[5], results.vars$ugrhuso[index], c(1:3)+versch[5], results.vars$ogrhuso[index], lwd=3, lend="butt") points(c(1:3)+versch[5], results.vars$biashuso[index], pch=symbole[5], bg="white", cex=2) axis(2, at=seq(-1, oblim, by=0.5), las=1, labels=NA, cex.axis=cex.axis, tck=0.015) box() par(mfg=c(2,4)) index <- results.vars$T==30 & results.vars$alpha>1 # decreasing removal plot(1:3, results.vars$antnullcount[index], type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8)) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.vars$antnullcount[index][i], border=NA, col=grey(0.8)) axis(4, at=seq(0, 1, by=.2), las=1, tck=0.015, cex.axis=cex.axis, labels=c(seq(0, 0.8, by=.2), NA)) par(new=TRUE) plot(c(1:3)+versch[1], results.vars$biascount[index], ylim=c(-1, oblim), las=1, tck=0.015, cex.axis=cex.axis, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.vars$ugrcount[index], c(1:3)+versch[1], results.vars$ogrcount[index], lwd=3, lend="butt") points(c(1:3)+versch[1], results.vars$biascount[index], pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.vars$ugreinf[index], c(1:3)+versch[2], results.vars$ogreinf[index], lwd=3, lend="butt") points(c(1:3)+versch[2], results.vars$biasein[index], pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.vars$ugroikostat[index], c(1:3)+versch[3], results.vars$ogroikostat[index], lwd=3, lend="butt") points(c(1:3)+versch[3], results.vars$biasoikostat[index], pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.vars$ugrerick[index], c(1:3)+versch[4], results.vars$ogrerick[index], lwd=3, lend="butt") points(c(1:3)+versch[4], results.vars$biaserick[index], pch=symbole[4], bg="white", cex=2) segments(c(1:3)+versch[5], results.vars$ugrhuso[index], c(1:3)+versch[5], results.vars$ogrhuso[index], lwd=3, lend="butt") points(c(1:3)+versch[5], results.vars$biashuso[index], pch=symbole[5], bg="white", cex=2) axis(2, at=seq(-1, oblim, by=0.5),labels=NA, las=1, cex.axis=cex.axis, tck=0.015) axis(1, at=1:3, labels=results.vars$d[index], tck=0, cex.axis=cex.axis*1.2) box() mtext(side=1, line=3.1,outer=TRUE, text="SEARCH INTERVAL (days)", cex=cex.lab) mtext(side=2, line=5,outer=TRUE, text="MEAN RELATIVE ERROR (symbols)", cex=cex.lab) mtext(side=4, line=3,outer=TRUE, text="PROPORTION OF COUNT = 0 (grey bars)", cex=cex.lab) savePlot("relErr_constm.varm.vars.bmp", type="bmp") #------------------------------------------------------------------------------- # influence of f ################################################################################ # 4. simulation: varying f f.vector <-c(0.5, 0.8) # searcher efficiency m.vector <- c(0.1) # average number of fatalities per night s.vector <- c(0.8) # persistence probability d.vector<-c(1, 7, 14) # search interval I.vector <- c(100) # study period in days results <- expand.grid(f=f.vector, m=m.vector, s=s.vector, I=I.vector, d=d.vector) nsim <- 1000 for(i in 1:nrow(results)){ d <- results$d[i] f <- results$f[i] start.t <- d # day of first search I<-results$I[i] n<-floor(I/d) R<-start.t+I-d s<-results$s[i] # persistence probability # mean persistence time t_quer<-1/(-log(s)) countsvec <- numeric(nsim) diff0<-numeric(nsim) diff1<-numeric(nsim) diff2<-numeric(nsim) diff2a <- numeric(nsim) diff3<-numeric(nsim) diff4<-numeric(nsim) muN <- results$m[i] for(b in 1:nsim){ n.neu_wahr<-rpois(R, muN) # new fatalities per day n_wahr<-numeric(R) n_wahr[1]<-rbinom(1, n.neu_wahr[1],s) # number of carcasses present at day 1 count<-rep(NA, R) sampling<-rep(FALSE, R) # identify days with searches index<-seq(start.t, start.t+I-d, by=d) sampling[index]<-TRUE sum(sampling) # control: should equal n if(start.t>1) for(k in 2:start.t) n_wahr[k]<-rbinom(1, n_wahr[k-1], s)+rbinom(1, n.neu_wahr[k],s) # number of carcasses for the following days without searches count[start.t]<-rbinom(1, n_wahr[start.t], prob=f) # first search n_wahr[start.t]<-n_wahr[start.t]-count[start.t] for(k in (start.t+1):R){ # true number of carcasses and counts for days with searches n_wahr[k]<-rbinom(1, n_wahr[k-1], s)+rbinom(1, n.neu_wahr[k],s) if(sampling[k]){ count[k]<-rbinom(1, n_wahr[k], prob=f) n_wahr[k]<-n_wahr[k]-count[k] } } N.wahr<-sum(n.neu_wahr[(R-I+1):R]) # true number of fatalities per study period count<-count[!is.na(count)] N.est0 <- sum(count) sint <- 0 for(o in 1:d){ sint <- sint+s^(d-o+1) } sint <- sint/d N.est1<-sum(count)/(sint*f) # naive estimate # new method presented in the article p2<-det.bat(s, f, d, n) N.est2<- sum(count)/p2 p2a<-det.batdecrease(s, f, d, n, k=0.25) N.est2a<- sum(count)/p2a # Erickson 2004 N.est3<-sum(count)/(t_quer*f/d*((exp(d/t_quer)-1)/(exp(d/t_quer)-1+f))) # Erickson 2004 # Huso 2010 p4<-phuso(s,f,d) N.est4<-sum(count)/p4 countsvec[b] <- sum(count, na.rm=TRUE) diff0[b]<-ifelse(N.wahr>0, (N.est0-N.wahr)/N.wahr, NA) diff1[b]<-ifelse(N.wahr>0, (N.est1-N.wahr)/N.wahr, NA) diff2[b]<-ifelse(N.wahr>0, (N.est2-N.wahr)/N.wahr, NA) diff2a[b]<-ifelse(N.wahr>0, (N.est2a-N.wahr)/N.wahr, NA) diff3[b]<-ifelse(N.wahr>0, (N.est3-N.wahr)/N.wahr, NA) diff4[b]<-ifelse(N.wahr>0, (N.est4-N.wahr)/N.wahr, NA) } # close b results$antnullcount[i] <- sum(countsvec==0)/nsim results$mcount[i] <- mean(countsvec) results$biascount[i] <- mean(diff0, na.rm=TRUE) results$ugrcount[i] <- lower(diff0) results$ogrcount[i] <- upper(diff0) results$sdcount[i] <- sd(diff0, na.rm=TRUE) results$biaseinf[i] <- mean(diff1, na.rm=TRUE) results$ugreinf[i] <- lower(diff1) results$ogreinf[i] <- upper(diff1) results$sdeinf[i] <- sd(diff1, na.rm=TRUE) results$biasoikostat[i] <- mean(diff2, na.rm=TRUE) results$ugroikostat[i] <- lower(diff2) results$ogroikostat[i] <- upper(diff2) results$sdoikostat[i] <- sd(diff2, na.rm=TRUE) results$biasoikostat.a[i] <- mean(diff2a, na.rm=TRUE) results$ugroikostat.a[i] <- lower(diff2a) results$ogroikostat.a[i] <- upper(diff2a) results$sdoikostat.a[i] <- sd(diff2a, na.rm=TRUE) results$biaserick[i] <- mean(diff3, na.rm=TRUE) results$ugrerick[i] <- lower(diff3) results$ogrerick[i] <- upper(diff3) results$sderick[i] <- sd(diff3, na.rm=TRUE) results$biashuso[i] <- mean(diff4, na.rm=TRUE) results$ugrhuso[i] <- lower(diff4) results$ogrhuso[i] <- upper(diff4) results$sdhuso[i] <- sd(diff4, na.rm=TRUE) write.table(results, file="simresults110414varf.txt", row.names=FALSE) } # close i ################################################################################ # 5. simulation: decreasubg f: f(t) = f*0.25^t f.vector <-c(0.8) # searcher efficiency m.vector <- c(0.1) # average number of fatalities per night s.vector <- c(0.8) # persistence probability d.vector<-c(1, 7, 14) # search interval I.vector <- c(100) # study period in days results <- expand.grid(f=f.vector, m=m.vector, s=s.vector, I=I.vector, d=d.vector) nsim <- 1000 for(i in 1:nrow(results)){ d <- results$d[i] f <- results$f[i] start.t <- d # day of first search I<-results$I[i] n<-floor(I/d) R<-start.t+I-d s<-results$s[i] # persistence probability # average persistence time t_quer<-1/(-log(s)) countsvec <- numeric(nsim) diff0<-numeric(nsim) diff1<-numeric(nsim) diff2<-numeric(nsim) diff2a<-numeric(nsim) diff3<-numeric(nsim) diff4<-numeric(nsim) muN <- results$m[i] for(b in 1:nsim){ n.neu_wahr<-rpois(R, muN) # number of new fatalities per night count<-rep(NA, R) sampling<-rep(FALSE, R+1) # identify days with searches index<-seq(start.t, start.t+I-d, by=d) sampling[index]<-TRUE sum(sampling) # control: should equal n n_true.mat <- matrix(ncol=R+1,nrow=R) count.mat <- matrix(ncol=R+1,nrow=R) diag(n_true.mat) <- rbinom(R, n.neu_wahr, prob=s) # number of bats still there after the first night for(z in 1:R){ if(sampling[z]){ # first day count.mat[z, z] <- rbinom(1, n_true.mat[z,z], prob=f) n_true.mat[z,z] <- n_true.mat[z,z] - count.mat[z, z] } for(zz in (z+1):(R+1)){ # following days n_true.mat[z,zz] <- rbinom(1, n_true.mat[z, zz-1], prob=s) if(sampling[zz]){ count.mat[z, zz] <- rbinom(1, n_true.mat[z,zz], prob=f*0.25^(sum(sampling[z:zz])-1)) n_true.mat[z,zz] <- n_true.mat[z,zz] - count.mat[z, zz] } } } count <- apply(count.mat, 2, sum, na.rm=TRUE) N.wahr<-sum(n.neu_wahr[(R-I+1):R]) # true number of fatalities count<-count[!is.na(count)] N.est0 <- sum(count) sint <- 0 for(o in 1:d){ sint <- sint+s^(d-o+1) } sint <- sint/d N.est1<-sum(count)/(sint*f) # simple estimator p2<-det.bat(s, f, d, n) # new method presented in the article N.est2<- sum(count)/p2 p2a<-det.batdecrease(s, f, d, n, k=0.25) N.est2a<- sum(count)/p2a # Erickson 2004 N.est3<-sum(count)/(t_quer*f/d*((exp(d/t_quer)-1)/(exp(d/t_quer)-1+f))) # Erickson 2004 # Huso 2010 p4<-phuso(s,f,d) N.est4<-sum(count)/p4 countsvec[b] <- sum(count, na.rm=TRUE) diff0[b]<-ifelse(N.wahr>0, (N.est0-N.wahr)/N.wahr, NA) diff1[b]<-ifelse(N.wahr>0, (N.est1-N.wahr)/N.wahr, NA) diff2[b]<-ifelse(N.wahr>0, (N.est2-N.wahr)/N.wahr, NA) diff2a[b]<-ifelse(N.wahr>0, (N.est2a-N.wahr)/N.wahr, NA) diff3[b]<-ifelse(N.wahr>0, (N.est3-N.wahr)/N.wahr, NA) diff4[b]<-ifelse(N.wahr>0, (N.est4-N.wahr)/N.wahr, NA) } # close b results$antnullcount[i] <- sum(countsvec==0)/nsim results$mcount[i] <- mean(countsvec) results$biascount[i] <- mean(diff0, na.rm=TRUE) results$ugrcount[i] <- lower(diff0) results$ogrcount[i] <- upper(diff0) results$sdcount[i] <- sd(diff0, na.rm=TRUE) results$biaseinf[i] <- mean(diff1, na.rm=TRUE) results$ugreinf[i] <- lower(diff1) results$ogreinf[i] <- upper(diff1) results$sdeinf[i] <- sd(diff1, na.rm=TRUE) results$biasoikostat[i] <- mean(diff2, na.rm=TRUE) results$ugroikostat[i] <- lower(diff2) results$ogroikostat[i] <- upper(diff2) results$sdoikostat[i] <- sd(diff2, na.rm=TRUE) results$biasoikostat.a[i] <- mean(diff2a, na.rm=TRUE) results$ugroikostat.a[i] <- lower(diff2a) results$ogroikostat.a[i] <- upper(diff2a) results$sdoikostat.a[i] <- sd(diff2a, na.rm=TRUE) results$biaserick[i] <- mean(diff3, na.rm=TRUE) results$ugrerick[i] <- lower(diff3) results$ogrerick[i] <- upper(diff3) results$sderick[i] <- sd(diff3, na.rm=TRUE) results$biashuso[i] <- mean(diff4, na.rm=TRUE) results$ugrhuso[i] <- lower(diff4) results$ogrhuso[i] <- upper(diff4) results$sdhuso[i] <- sd(diff4, na.rm=TRUE) write.table(results, file="simresults110414fdecreasing.txt", row.names=FALSE) } # close i #------------------------------------------------------------------------------- # Draw figures results.varf <- read.table("simresults110414varf.txt", header=TRUE) results.decf <- read.table("simresults110414fdecreasing.txt", header=TRUE) versch <- seq(-0.42, 0.38, length=6) symbole <- c(20, 21, 22, 22, 23, 24) oblim <- 3 cex.main <- 1.5 col.oikostat <- grey(0.6) col.oikostat.a <- grey(0.9) cex.axis <- 1.4 cex.lab <- 1.7 windows(width=10, height=5) par(mar=c(0.1, 0.5, 3, 0.1), mgp=c(1, 0.3, 0), oma=c(6, 5, 0, 6), mfrow=c(1,3)) for(x in 1:2){ index <- results.varf$f==c(0.5, 0.8)[x] plot(1:3, results.varf$antnullcount[index], type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8), main=paste("f =", c(0.5, 0.8)[x], "constant"), cex.main=cex.main) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.varf$antnullcount[index][i], border=NA, col=grey(0.8)) axis(4, las=1, labels=NA, tck=0.015) par(new=TRUE) plot(c(1:3)+versch[1], results.varf$biascount[index], ylim=c(-1, oblim), las=1, tck=0.015, cex.lab=1.2, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.varf$ugrcount[index], c(1:3)+versch[1], results.varf$ogrcount[index], lwd=3, lend="butt") points(c(1:3)+versch[1], results.varf$biascount[index], pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.varf$ugreinf[index], c(1:3)+versch[2], results.varf$ogreinf[index], lwd=3, lend="butt") points(c(1:3)+versch[2], results.varf$biasein[index], pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.varf$ugroikostat[index], c(1:3)+versch[3], results.varf$ogroikostat[index], lwd=3, lend="butt") points(c(1:3)+versch[3], results.varf$biasoikostat[index], pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.varf$ugroikostat.a[index], c(1:3)+versch[4], results.varf$ogroikostat.a[index], lwd=3, lend="butt") points(c(1:3)+versch[4], results.varf$biasoikostat.a[index], pch=symbole[4], bg=col.oikostat.a, cex=2) segments(c(1:3)+versch[5], results.varf$ugrerick[index], c(1:3)+versch[5], results.varf$ogrerick[index], lwd=3, lend="butt") points(c(1:3)+versch[5], results.varf$biaserick[index], pch=symbole[5], bg="white", cex=2) segments(c(1:3)+versch[6], results.varf$ugrhuso[index], c(1:3)+versch[6], results.varf$ogrhuso[index], lwd=3, lend="butt") points(c(1:3)+versch[6], results.varf$biashuso[index], pch=symbole[6], bg="white", cex=2) if(x==1) axis(2, at=seq(-1, oblim, by=0.5), las=1, cex.axis=cex.axis*0.8, tck=0.015) if(x==1) mtext(side=2, text="MEAN RELATIVE ERROR (symbols)", line=2.5, cex=cex.axis, outer=TRUE) axis(1, at=1:3, labels=results.varf$d[index], cex.axis=cex.axis, tck=0) box() } plot(1:3, results.decf$antnullcount, type="n", ylim=c(0, 1), xlim=c(0.4,3.6), xlab="", ylab="", axes=FALSE, las=1, cex.axis=cex.axis, xaxs="i", yaxs="i", col=grey(0.8), main="f decreasing", cex.main=cex.main) for(i in 1:3) rect(c(0.51, 1.51, 2.51)[i], 0, c(1.49, 2.49, 3.49)[i], results.decf$antnullcount[i], border=NA, col=grey(0.8)) axis(4, las=1, tck=0.015, cex.axis=cex.axis) par(new=TRUE) plot(c(1:3)+versch[1], results.decf$biascount, ylim=c(-1, oblim), las=1, tck=0.015, cex.lab=1.2, xaxs="i", yaxs="i", axes=FALSE, ylab="", xlab="", xlim=c(0.4,3.6)) abline(h=0, lty=3, lwd=2) segments(c(1:3)+versch[1], results.decf$ugrcount, c(1:3)+versch[1], results.decf$ogrcount, lwd=3, lend="butt") points(c(1:3)+versch[1], results.decf$biascount, pch=symbole[1], bg="white", cex=2) segments(c(1:3)+versch[2], results.decf$ugreinf, c(1:3)+versch[2], results.decf$ogreinf, lwd=3, lend="butt") points(c(1:3)+versch[2], results.decf$biasein, pch=symbole[2], bg="white", cex=2) segments(c(1:3)+versch[3], results.decf$ugroikostat, c(1:3)+versch[3], results.decf$ogroikostat, lwd=3, lend="butt") points(c(1:3)+versch[3], results.decf$biasoikostat, pch=symbole[3], bg=col.oikostat, cex=2) segments(c(1:3)+versch[4], results.decf$ugroikostat.a, c(1:3)+versch[4], results.decf$ogroikostat.a, lwd=3, lend="butt") points(c(1:3)+versch[4], results.decf$biasoikostat.a, pch=symbole[4], bg=col.oikostat.a, cex=2) segments(c(1:3)+versch[5], results.decf$ugrerick, c(1:3)+versch[5], results.decf$ogrerick, lwd=3, lend="butt") points(c(1:3)+versch[5], results.decf$biaserick, pch=symbole[5], bg="white", cex=2) segments(c(1:3)+versch[6], results.decf$ugrhuso, c(1:3)+versch[6], results.decf$ogrhuso, lwd=3, lend="butt") points(c(1:3)+versch[6], results.decf$biashuso, pch=symbole[6], bg="white", cex=2) if(x==1) axis(2, at=seq(-1, oblim, by=0.5), las=1, cex.axis=cex.axis*0.8, tck=0.015) if(x==1) mtext(side=2, text="MEAN RELATIVE ERROR (symbols)", line=2.5, cex=cex.axis, outer=TRUE) axis(1, at=1:3, labels=results.decf$d, cex.axis=cex.axis, tck=0) box() legend(0.5, 2.9, pch=c(20, 21, 22, 22, 23, 24), pt.bg=c("white", "white", col.oikostat, col.oikostat.a, "white", "white"), legend=c("Count", "Simple formula", "This article", "This article, adapted for decr. f", "Erickson et al. (2004)", "Huso (2010)"), cex=1.3) mtext(side=4, line=2.5, text="PROPORTION OF COUNTS = 0 (grey bars)", cex=cex.lab*0.8) mtext(side=1, line=2.5, text="DAY", cex=cex.lab, outer=TRUE) #------------------------------------------------------------------------------- # simulation 6: show how precision changes with the number of counted carcasses f.vector <-seq(0.05, 0.95, by=0.1) # searcher efficiency m.vector <- seq(0.01, 1, length=10) # average number of fatalities per night d.vector<-c(7) # search interval I.vector <- c(100) # study period in days T.vector <- c(30) # average persistence time results <- expand.grid(f=f.vector, m=m.vector, T=T.vector, I=I.vector, d=d.vector) results$s <- exp(1/-results$T) results$n<-floor(results$I/results$d) for(i in 1:nrow(results)) results$p[i] <- det.bat(results$s[i], results$f[i], results$d[i], results$n[i]) nsim <- 1000 for(i in 1:nrow(results)){ d <- results$d[i] f <- results$f[i] start.t <- d # day of first search I<-results$I[i] n <- results$n[i] R<-start.t+I-d s<-results$s[i] # persistence probability countsvec <- numeric(nsim) diff2<-numeric(nsim) muN <- results$m[i] for(b in 1:nsim){ n.neu_wahr<-rpois(R, muN) # new fatalities per day n_wahr<-numeric(R) n_wahr[1]<-rbinom(1, n.neu_wahr[1],s) # number of carcasses at day 1 count<-rep(NA, R) sampling<-rep(FALSE, R) # identify days with searches index<-seq(start.t, start.t+I-d, by=d) sampling[index]<-TRUE if(start.t>1) for(k in 2:start.t) n_wahr[k]<-rbinom(1, n_wahr[k-1], s)+rbinom(1, n.neu_wahr[k],s) count[start.t]<-rbinom(1, n_wahr[start.t], prob=f) # first search n_wahr[start.t]<-n_wahr[start.t]-count[start.t] for(k in (start.t+1):R){ n_wahr[k]<-rbinom(1, n_wahr[k-1], s)+rbinom(1, n.neu_wahr[k],s) if(sampling[k]){ count[k]<-rbinom(1, n_wahr[k], prob=f) n_wahr[k]<-n_wahr[k]-count[k] } } N.wahr<-sum(n.neu_wahr[(R-I+1):R]) # true number of fatalites per study period count<-count[!is.na(count)] N.est2<- sum(count)/results$p[i] countsvec[b] <- sum(count, na.rm=TRUE) diff2[b]<-ifelse(N.wahr>0, (N.est2-N.wahr)/N.wahr, NA) } # close b results$antnullcount[i] <- sum(countsvec==0)/nsim results$mcount[i] <- mean(countsvec) results$biasoikostat[i] <- mean(diff2, na.rm=TRUE) results$ugroikostat[i] <- lower(diff2) results$ogroikostat[i] <- upper(diff2) results$sdoikostat[i] <- sd(diff2, na.rm=TRUE) write.table(results, file="simresults110417_ncounts.txt", row.names=FALSE) } # close i results <- read.table("simresults110417_ncounts.txt", header=TRUE) plot( results$mcount, results$p) par(mar=c(5,7,1,1)) plot(results$mcount, results$sdoikostat, type="n", las=1, xlab="NUMBER OF CARCASSES FOUND", ylab="STANDARD DEVIATION OF\nTHE RELATIVE ERROR", cex.axis=1.2, cex.lab=1.4, ylim=c(0,max(results$sdoikostat)), xlim=c(-5, max(results$mcount)), bty="l") for(i in unique(results$p)){ lines(results$mcount[results$p==i], results$sdoikostat[results$p==i], lwd=1.5) if(!is.element(round(i,2),c(0.72, 0.81, 0.84))) text(-4, results$sdoikostat[results$p==i][1], round(i, 2)) }