#R code accompanying Tyack and Thomas dose:response impact paper #Code calculates expected impact using Miller et al. (2014) as and example #Also reproduces the figures in the paper #Last updated: 17th August 2018 #Set to true if you want figures produced and saved to png files make.figures<-TRUE #variables used to scale the figure sizes figure.x.multiplier<-7 figure.y.multiplier<-5.5 figure.pointsize<-23 #Input data from Table 4 of Miller et al. #Posterior mean p(response) p<-c(0,0.01,0.03,0.06,0.10,0.17,0.26,0.37,0.49,0.62,0.74,0.83,0.91,0.96,1) #Quantiles of p(response) - 25, 75, 0.025 and 0.975 respectively #Used for Figure 1 q<-matrix(0,4,length(p)) q[1,]<-c(0,0,0.01,0.02,0.05,0.09,0.17,0.26,0.38,0.52,0.66,0.78,0.88,0.95,1) q[2,]<-c(0,0.015,0.04,0.08,0.14,0.23,0.34,0.47,0.60,0.73,0.84,0.91,0.96,0.99,1) q[3,]<-c(0,0,0,0,0.01,0.03,0.06,0.12,0.2,0.3,0.44,0.58,0.73,0.87,1) q[4,]<-c(0,0.05,0.11,0.20,0.30,0.43,0.56,0.68,0.79,0.87,0.94,0.97,0.99,1,1) levels<-seq(60,200,by=10) quants<-seq(0,1,length=length(levels)) pred.levels<-seq(60,200,by=0.01) #Interpolate input data with splines so we can estimate p(response) and # quantiles for any dose #Choose which quantile to use as input - # e.g., input.p<-p for posterior mean, or input.p<-q[3,] for 0.025 quantile input.p<-p #q[4,] mod.spline.p<-splinefun(levels,input.p) pred.p<-mod.spline.p(pred.levels) pred.q<-matrix(0,4,length(pred.p)) for(i in 1:4){ mod.spline.q<-splinefun(levels,q[i,]) pred.q[i,]<-mod.spline.q(pred.levels) } #Work out RL.p50 and RL.p10 p50.index<-which.max(pred.p>0.5) RL.p50<-pred.levels[p50.index] p10.index<-which.max(pred.p>0.1) RL.p10<-pred.levels[p10.index] if(make.figures){ png("fig1.png", width=480*figure.x.multiplier, height=480*figure.y.multiplier, pointsize=figure.pointsize,res=300) #Plot predicted levels vs p response plot(pred.levels,pred.p,type="l",xlab="Received Level (dB re 1 microPascal)",ylab="Probability of Response",lwd=3,col="black") red.dot<-which.max(pred.p>0.5) lines(pred.levels,pred.q[3,],lty=3,lwd=2) lines(pred.levels,pred.q[4,],lty=3,lwd=2) points(RL.p50,0.5,pch=16,cex=2,col="red") arrows(60,0.5,RL.p50-4,0.5,col="red",lwd=2) arrows(RL.p50,0.5-0.05,RL.p50,0,col="red",lwd=2) points(RL.p10,0.1,pch=16,cex=2,col="blue") arrows(60,0.1,RL.p10-4,0.1,col="blue",lwd=2) arrows(RL.p10,0.1-0.05,RL.p10,0,col="blue",lwd=2) dev.off() } #Get a transmission loss curve #Note: assumes input is in kilometers #a is absorption coefficient (/km, default assumes 3khZ signal) TL<-function(r,a=0.185){ loss<-20*log10(r*1000) loss[loss<0]<-0 loss<-loss+a*r return(loss) } #Code used to find range as a function of TL, but now I first find the max range # and then just work across equally spaced ranges #Find max range to use obj.fun<-function(r,SL,min.RL=60){ #Return squared difference between target.TL (SL-60) and actual TL target.TL<-SL-min.RL return((target.TL - TL(r))^2) } SL<-210 #Assumed source level ret<-optimize(obj.fun,interval=c(0,30000),SL=SL) #Found up max range to nearest higher 10km max.r<-ceiling(ret$minimum/10)*10 #Also use the same technique to get the p50 and p10 ranges ret<-optimize(obj.fun,interval=c(0,30000),SL=SL,min.RL=RL.p50) r.p50<-ret$minimum ret<-optimize(obj.fun,interval=c(0,30000),SL=SL,min.RL=RL.p10) r.p10<-ret$minimum #Work out RL in a set of bins n.bins<-10000 range.interval<-max.r/n.bins cutpoints<-seq(0,max.r,length=(n.bins+1)) r<-cutpoints[-1]-range.interval/2 TLs<-TL(r) RL<-SL-TLs #Get the corresponding p new.p<-mod.spline.p(RL) #constrain in the 0-1 range, in case it goes outside new.p[new.p<0]<-0; new.p[new.p>1]<-1 if(make.figures){ png("fig2.png", width=480*figure.x.multiplier, height=480*figure.y.multiplier, pointsize=figure.pointsize,res=300) #Plot recevied level as a function of range plot(r,RL,type="l",xlab="Range (km)",ylab="Received Level (db re 1 microPascal)",lwd=2,yaxt="n") axis(2,at=seq(60,200,by=20),labels=c("60","","100","","140","","180","")) points(r.p50,RL.p50,pch=16,cex=1.5,col="red") arrows(-7,RL.p50,r.p50-5,RL.p50,col="red",lwd=2) arrows(r.p50,RL.p50-5,r.p50,60,col="red",lwd=2) text(r.p50+5,RL.p50,"50% respond",pos=4,col="red") points(r.p10,RL.p10,pch=16,cex=1.5,col="blue") arrows(-7,RL.p10,r.p10-5,RL.p10,col="blue",lwd=2) arrows(r.p10,RL.p10-5,r.p10,60,col="blue",lwd=2) text(r.p10+5,RL.p10,"10% respond",pos=4,col="blue") dev.off() } #Work out animal density as a function of range D<-1 #Animal density - per square kilometer #Area of each bin bin.area<-pi*(cutpoints[-1]^2-cutpoints[1:n.bins]^2) bin.N<-D*bin.area bin.n.respond<-bin.N*new.p n.respond<-sum(bin.n.respond) #Print number responding print(n.respond) #Work out effective response range n.cum<-cumsum(bin.N) n.respond.cum<-cumsum(bin.n.respond) n.norespond.cum<-n.cum-n.respond.cum #get the difference between the number that don't respond at smaller ranges # and the number that do at larger ranges diff<-(n.norespond.cum-(sum(n.respond)-n.respond.cum)) #look for the place where this difference is 0 - that's the ERR mod.spline.ERR<-splinefun(diff,r) ERR<-mod.spline.ERR(0) #Note - another way to work this out (easier to compute) is sqrt(D*n.respond/pi). #Print Effective Response Radius print(ERR) #Get the effective received level associated with this ERL<-SL-TL(ERR) #Print Effective Received Level print(ERL) #Get the detection probability associated with this ERp<-mod.spline.p(ERL) if(make.figures){ #work out which of the r's is more than r.p50 and ERR r.p50.index<-which.max(r>=r.p50) ERR.index<-which.max(r>=ERR) png("fig3.png", width=480*figure.x.multiplier, height=480*figure.y.multiplier, pointsize=figure.pointsize,res=300) #Plot number of animals and number responding against range plot(r,bin.n.respond,type="l",xlab="Range bin midpoint (km)",ylab="Number of animals",lwd=2,ylim=c(0,bin.N[ERR.index]*1.05)) polygon(x=c(r[1:ERR.index],r[seq(ERR.index,1,by=-1)]), y=c(bin.n.respond[1:ERR.index],bin.N[seq(ERR.index,1,by=-1)]), col=gray(0.1),border=NA,density=20,angle=-45) l.r<-length(r) polygon(x=c(r[ERR.index:l.r],r[seq(l.r,ERR.index,by=-1)]), y=c(bin.n.respond[ERR.index:l.r],rep(0,l.r-ERR.index+1)), col=gray(0.1),border=NA,density=20,angle=45) lines(r,bin.N,lty=2,lwd=2) points(r[r.p50.index],bin.n.respond[r.p50.index],pch=16,cex=1.5,col="red") points(r[ERR.index],bin.n.respond[ERR.index],pch=16,cex=1.5,col="green") dev.off() } if(make.figures){ #This code produces Figure 4 - it is not needed to calculate ERR, ERL, etc #Work out the number of animals (rounded) responding at each range #Note the actual number of animals is scaled down (N<- line) so the plot is clear set.seed(14318) N<-round(sum(bin.N))/30 animal.pos<-data.frame(id=1:N,r=0,phi=0,x=0,y=0,p=0,respond=FALSE) #Assign distances deterministicly according to a triangular distribution cum.area.per.animal<-cumsum(rep(pi*max.r^2/N,N)) animal.pos$r<-sqrt(cum.area.per.animal/pi) #Random angles animal.pos$phi<-runif(N,0,2*pi) #Turn from polar coords into cartesian animal.pos$x<-animal.pos$r*sin(animal.pos$phi) animal.pos$y<-animal.pos$r*cos(animal.pos$phi) #Simulate whether it responded animal.pos$p<-mod.spline.p(SL-TL(animal.pos$r)) animal.pos$p[animal.pos$p<0]<-0;animal.pos$p[animal.pos$p>1]<-1 animal.pos$respond<-(runif(N)1]<-1 color.ramp.scale<-0.5 for(i in 1:(n.ranges.for.p-1)){ my.color<-hsv(2/3,ps[i]^color.ramp.scale,1) filledcircle(r2=ranges[i],r1=ranges[i+1],mid=c(0,0),col=my.color) } #Draw marker lines at 10km intervals for(ri in seq(0,max.r,by=10)){ plotcircle(r=ri,mid=c(0,0),lcol="black",lwd=1) } #animals points(animal.pos$x,animal.pos$y,pch=20,cex=1) #animals that respond are in red points(animal.pos$x[animal.pos$respond],animal.pos$y[animal.pos$respond],col="red",pch=20,lwd=2,cex=1) legend("topright",title="p(response)",legend=seq(1,0,by=-0.2),fill=hsv(2/3,(seq(1,0,by=-0.2))^color.ramp.scale,1)) par(mar=c(4,4,0,0)+0.1) plot(0,0,xlim=c(-max.r2,max.r2),ylim=c(-max.r2,max.r2),type="n",xlab="x (km)",ylab="y (km)") #Draw marker lines at 10km intervals for(ri in seq(0,max.r,by=10)){ plotcircle(r=ri,mid=c(0,0),lcol="black",lwd=1) } #Draw ERR plotcircle(r=ERR,mid=c(0,0),lcol="green",lwd=5) #animals points(animal.pos$x,animal.pos$y,pch=20,cex=1) #animals within ERR are in red points(animal.pos$x[animal.pos$r<=ERR],animal.pos$y[animal.pos$r<=ERR],col="red",pch=20,lwd=2,cex=1) dev.off() par(pars) } #Calculate uncertainty on number of takes, by reading in 1000 replicate dose-response functions reps<-read.csv(file="dose_response_replicates.csv", row.names=1, stringsAsFactors = FALSE) n.reps<-dim(reps)[1] rep.ERR<-rep.n.respond<-numeric(n.reps) for(i in 1:n.reps){ mod.spline.p<-splinefun(levels,reps[i,]) rep.p<-mod.spline.p(RL) rep.n.respond[i]<-sum(bin.N*rep.p) rep.ERR[i]<-sqrt(D*rep.n.respond[i]/pi) } #95% interval on number affected print(quantile(rep.n.respond,c(0.025,0.975))) #95% interval on ERR print(quantile(rep.ERR,c(0.025,0.975))) #95% interval on ERL print(SL-TL(quantile(rep.ERR,c(0.975,0.025))))