#clear all variables rm(list=ls()) library(data.table) library(ggplot2) library(cowplot) library(vegan) library(Distance) library(AICcmodavg) #set the foldername to the path where you have all scripts and data files setwd("xxxxxxxxxxxxxxxxx") # distance estimation analysis -------------------------------------------- distance0=fread("points vs sound/distance estimation data.csv") #correct distances according to Pythagorean theorem to get direct distances instead of horizontal #estimated height of recorder: 1.5m distance0[,distance_m_direct:=sqrt(distance_m^2+1.5^2)] distance0[call_type!="",Species_call_type:=paste(Species,"\n","(",call_type,")",sep="")] distance0[call_type=="",Species_call_type:=Species] distance1=distance0[exclude==""] #count total and missing recordings nrow(distance0[Estimator=="Estimator 1\nfamiliar with birds"]) nrow(distance0[exclude=="not found in recording" & Estimator=="Estimator 1\nfamiliar with birds"]) #count calls per species call type distance0[Estimator=="Estimator 1\nfamiliar with birds" & exclude!="not found in recording" ,.N,Species_call_type] #reorder factor distance1[,Estimator:=factor(Estimator,levels=c("Estimator 1\nfamiliar with birds" ,"Estimator 2\nclosest distance known" ,"Estimator 2\nclosest and farthest distance known"))] T1=distance1[,.(Slope=round(coef(lm(distance_audio~distance_m_direct))[2],2) ,Intercept=round(coef(lm(distance_audio~distance_m_direct))[1]) ,Correlation=round(cor(distance_m_direct,distance_audio),2) ,Calls=.N) ,.(`Species and call type`=Species_call_type,Estimator,Habitat=habitat)] fwrite(T1[order(`Species and call type`,Estimator)],"points vs sound/Figures field/T1.csv") ggplot(distance1,aes(distance_m_direct,distance_audio,color=Species_call_type,label=ID))+ geom_point(alpha=0.5)+ scale_color_discrete(name="Bird species and call type")+ stat_smooth(method="lm",se=F)+ geom_abline(slope=1,intercept=0,lty=2)+ labs(x="Measured distance to bird (m)",y="Estimated distance (m)\nto bird call in sound recording")+ facet_grid(.~Estimator)+ theme(legend.key.size = unit(2, 'lines')) ggsave("points vs sound/Figures field/Fig 1.tiff",width=14,height=6) # Distance to bird detections profiles ------------------------------------------------------- # import data points.sounds=fread("points vs sound/field survey data.csv") points.sounds=points.sounds[Comment==""] #check whether replicates are balanced points.sounds[,.(replicates=length(unique(method)) ,detections=.N),Observer_plot] #total detections per method points.sounds[,.(detections=.N),method] max_distance=60 bin_area=600 area_bins=seq(0,pi*(max_distance^2),bin_area) distance_bins=round(sqrt(area_bins/pi)) par(mfrow=c(1,2)) #horizontal distance histograms hist(points.sounds[Horizontal_distance<=max_distance & method=="Sound recordings" & Evidence!="Visual",Horizontal_distance],freq=T,breaks=seq(0,80,10),main="horizontal distance, sound recordings") hist(points.sounds[Horizontal_distance<=max_distance & method=="Point counts" & Evidence!="Visual",Horizontal_distance],freq=T,breaks=seq(0,80,10),main="horizontal distance, point counts") par(mfrow=c(1,1)) ggplot(points.sounds[Horizontal_distance<=max_distance & Evidence!="Visual"],aes(pi*(Horizontal_distance)^2))+ geom_histogram(binwidth=bin_area,boundary=-0.5)+ #scale_y_continuous(breaks=seq(0,1,0.05),labels=percent)+ scale_x_continuous(breaks=area_bins,labels=distance_bins)+ labs(y=paste("Frequency of detections"),x="Distance from observer, bins of equal area")+ facet_grid(method~.) # Richness and abundance comparison (preparation) ------------------------------------------------- densities=fread("points vs sound/hds birds per ha per count.csv") #quick summary densities[,.(density_km2_AR=mean(hds.ar)*100 ,density_km2_PC=mean(hds.pc)*100)] #cut-off distance to get fixed ranges between methods #choose threshold distance above which there are only 5% of observations distance.sound0=points.sounds[method=="Sound recordings"] distance.sound1=distance.sound0[order(Horizontal_distance)] distance_limit=distance.sound1[nrow(distance.sound1)-round(nrow(distance.sound1)*0.05),"direct_distance_m",with=F][[1]] dt0=rbind(points.sounds[Horizontal_distance<=distance_limit & method=="Sound recordings" ,.(Abundance=max(Individuals),range="fixed") ,.(Species,Plot,method,Evidence,Observer,Observer_plot)] ,points.sounds[method=="Sound recordings" ,.(Abundance=max(Individuals),range="unlimited") ,.(Species,Plot,method,Evidence,Observer,Observer_plot)] ,points.sounds[Horizontal_distance<=distance_limit & method=="Point counts" ,.(Abundance=max(Individuals),range="fixed") ,.(Species,Plot,method,Evidence,Observer,Observer_plot)] ,points.sounds[method=="Point counts" ,.(Abundance=max(Individuals),range="unlimited") ,.(Species,Plot,method,Evidence,Observer,Observer_plot)]) richness.abundance=dt0[ #alternative dataset with only aural detections # Evidence=="Aural" ,.(Richness=length(unique(Species)) ,Abundance=sum(Abundance)) ,.(Plot,method,Observer_plot,range)] #optionally, remove outlier # richness.abundance=richness.abundance[Plot!="F26"] #manually add plot with no detections richness.abundance=rbind(richness.abundance,data.table(Plot="F09",method="Point counts",Observer_plot="edho F09" ,range="fixed",Richness=0,Abundance=0)) #melt abundance and richness dt1.melt0=melt(richness.abundance,id=c("Plot","method","Observer_plot","range") ,measure=c("Abundance","Richness")) # extracting data for meta-analysis in JAPPL ------------------------------------------------------- # alpha=dt1.melt0[variable=="Richness" # ,.(alpha=mean(value),sd=sd(value),N=.N) # ,.(method,range)] # # # total observed richness # gamma0=dt0[,.(gamma=length(unique(Species))),.(method,range)] # # #unique in sound unlimited range # length(unique(dt0[!(Species %in% dt0[method=="Point counts" & range=="unlimited",Species]) & method=="Sound recordings" & range=="unlimited",Species])) # #unique in points unlimited range # length(unique(dt0[!(Species %in% dt0[method=="Sound recordings" & range=="unlimited",Species]) & method=="Point counts" & range=="unlimited",Species])) # #unique in sound fixed range # length(unique(dt0[!(Species %in% dt0[method=="Point counts" & range=="fixed",Species]) & method=="Sound recordings" & range=="fixed",Species])) # #unique in points fixed range # length(unique(dt0[!(Species %in% dt0[method=="Sound recordings" & range=="fixed",Species]) & method=="Point counts" & range=="fixed",Species])) # # #proportion of visually detected birds # nrow(points.sounds[Evidence=="Visual" & method=="Point counts"])/nrow(points.sounds[method=="Point counts"]) # #number of detections # sightings=points.sounds[,.(.N),method] # Richness and abundance comparison (tests and graph) ----------------------------------------------------- #adapt density data format densities.melt=melt(densities[,.(Observer_plot=site,hds.ar,hds.pc)],id="Observer_plot",variable.name="original_variable") densities.melt[grepl(".ar",original_variable),method:="Sound recordings"] densities.melt[grepl(".pc",original_variable),method:="Point counts"] densities.melt[,range:="Distance sampling\n(birds/ha)"] densities.melt[,variable:="Abundance"] #normality # par(mfrow=c(2,2)) # qqnorm(richness.abundance[method=="Sound recordings",Richness],main="Richness from AR") # qqline(richness.abundance[method=="Sound recordings",Richness]) # qqnorm(richness.abundance[method=="Point counts",Richness],main="Richness from PC") # qqline(richness.abundance[method=="Point counts",Richness]) # par(mfrow=c(1,1)) # # #are variances different? # var.test(richness.abundance[method=="Sound recordings" & range=="fixed",Richness] # ,richness.abundance[method=="Point counts" & range=="fixed",Richness]) #check number of rows # check=richness.abundance[,.N,Plot] # nrow(check[N==8]) richness.abundance=rbind(richness.abundance ,data.table(Plot="F09",method="Point counts",Observer_plot="edho F09",range="fixed","Richness"=0,"Abundance"=0)) #Wilcoxon signed-rank is appropriate for paired samples where populations do not have to be normal and do not have to have equal variances test_fixed_richness=wilcox.test(richness.abundance[method=="Sound recordings" & range=="fixed",Richness] ,richness.abundance[method=="Point counts" & range=="fixed",Richness],paired=T,conf.int=T)$p.value test_unlimited_richness=wilcox.test(richness.abundance[method=="Sound recordings" & range=="unlimited",Richness] ,richness.abundance[method=="Point counts" & range=="unlimited",Richness],paired=T,conf.int=T)$p.value test_fixed_abundance=wilcox.test(richness.abundance[method=="Sound recordings" & range=="fixed",Abundance] ,richness.abundance[method=="Point counts" & range=="fixed",Abundance],paired=T,conf.int=T)$p.value test_unlimited_abundance=wilcox.test(richness.abundance[method=="Sound recordings" & range=="unlimited",Abundance] ,richness.abundance[method=="Point counts" & range=="unlimited",Abundance],paired=T,conf.int=T)$p.value test_ds_abundance=wilcox.test(densities.melt[method=="Sound recordings",value] ,densities.melt[method=="Point counts",value],paired=T,conf.int=T)$p.value p_values0=c(test_ds_abundance,test_fixed_abundance,test_fixed_richness,test_unlimited_abundance,test_unlimited_richness) p_values1=round(p_values0,2) #keep trailing zeroes p_values2=data.table(annotation=paste("P-value =",sprintf("%.2f",p_values1)) ,range=c("Distance sampling\n(birds/ha)","Fixed radius","Fixed radius","Unlimited radius","Unlimited radius") ,variable=c("Abundance","Abundance","Richness","Abundance","Richness") ,value=c(max(densities.melt[,value]),max(richness.abundance[range=="fixed",Abundance]),max(richness.abundance[range=="fixed",Richness]),max(richness.abundance[range=="unlimited",Abundance]),max(richness.abundance[range=="unlimited",Richness])) ,method="Sound\nrecordings") #manual adjustment p_values2[annotation=="P-value = 0.00",annotation:="P-value < 0.01"] dt1.melt1=rbind(dt1.melt0,densities.melt,fill=T) dt1.melt1[,method_label0:=tstrsplit(method," ")[1]] dt1.melt1[,method_label1:=tstrsplit(method," ")[2]] dt1.melt1[,method_label:=paste(method_label0,"\n",method_label1,sep="")] dt1.melt1[,c("method_label0","method_label1"):=.(NULL,NULL)] dt1.melt1[range=="unlimited",range:="Unlimited radius"] dt1.melt1[range=="fixed",range:="Fixed radius"] #Figure 2 ggplot(dt1.melt1,aes(method_label,value,group=Observer_plot))+ geom_boxplot(width=0.5,aes(group=method),outlier.size = 0)+ geom_line(alpha=0.2)+ geom_text(data=p_values2,aes(method,value,label=annotation),group=1,nudge_x = -0.5)+ geom_point(alpha=0.2)+ facet_grid(range~variable,scales="free",switch="both")+ labs(x=NULL,y=NULL) ggsave("points vs sound/Figures field/Fig 3.tiff",width=6,height=7) # Crosscheck ------------------------------------------------------------------- #number of detections at close range for both methods points.sounds[direct_distance_m<=18,.(detections=.N),method] points.sounds[Horizontal_distance<=10,distance_class:="0-10"] points.sounds[Horizontal_distance<=20 & Horizontal_distance>10,distance_class:="10-20"] points.sounds[Horizontal_distance<=30 & Horizontal_distance>20,distance_class:="20-30"] points.sounds[Horizontal_distance<=40 & Horizontal_distance>30,distance_class:="30-40"] points.sounds[Horizontal_distance<=50 & Horizontal_distance>40,distance_class:="40-50"] points.sounds[Horizontal_distance<=60 & Horizontal_distance>50,distance_class:="50-60"] points.crosscheck.detail=points.sounds[check!="" & method=="Point counts",.(detections=.N) ,.(check,distance_class)] points.crosscheck.total=points.sounds[check!="" & method=="Point counts",.(detections_total=.N) ,.(distance_class)] #merge both tables to calculate proportions points.crosscheck=merge(points.crosscheck.detail,points.crosscheck.total ,by=c("distance_class")) points.crosscheck[,detection_proportion:=detections/detections_total] points.crosscheck[,check:=factor(check,levels=c("not heard","faint","heard"))] points.crosscheck1=points.crosscheck[order(check)] #Fig S3 ggplot(points.crosscheck1,aes(distance_class,detection_proportion,fill=check,label=detections_total))+ geom_bar(stat="identity")+ geom_text(aes(y=1.03))+ scale_fill_manual(values=c("red","orange","darkgreen"),name="Bird vocalization")+ labs(y="Proportion of detected birds",x="Horizontal distance (m)") ggsave("points vs sound/Figures field/Fig S3.png",width=6,height=4)