################################################################################################### ## CUSTOM CODE TO PERFORM A SENSITIVITY ANALYSIS OF THE BIRD-FRIENDLINESS INDEX ## ## This code performs sensitivity analyses of BFI components. Specifically, it evaluates how ## species density, conservation weight and functional evenness affect the overall output of the BFI ## ## AUTHOR: Nicole Michel, National Audubon Society #################################################################################################### ##read in required libraries library(raster) library(rgdal) library(MASS) library(fitdistrplus) library(overlapping) library(psych) # set working directory setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # define the path to where your input and output data files are stored pathtooutput <- "/Outputs/DensityPrediction/" pathtofiles <- "/Outputs/" pathtofd <- paste0(pathtofiles, "/FunctionalDiversity/") pathtodensity <- "/Outputs/DensityPrediction/Predictions/" pathtobfi <- paste0(pathtofiles, "BFI_estimated/") # create species list speclist <- read.csv(file="/CountData/BFI_NumberOfDetectionsBySpeciesYear.csv") # extract species list with sufficient detections speclist <- speclist$BirdCode[which(speclist$DetectTot >= 60)] # create yearlist yearlist <- c(2009:2014) ############################################################################################ ####THIS FIRST SECTION WILL CALCULATE THE EFFECT OF CHANGING THE STANDARD DEVIATION ABOUT THE MEAN ####DENSITY OF EACH SPECIES BY 10%; FINAL RESULT IS CALLED 'ES_DENSITY' (EFFECT SIZE DENSITY) ############################################################################################ #read in predicted density for each year, subset to include just the species density estimates #and then change the column names to include just the species alpha codes for (year in 2009:2014){ AbundFD <- read.delim(file=paste(pathtofd,"FinalDensities_FD_BFI_",year,".txt", sep=""),header=T) PDyr <- AbundFD[,grep("*DensityKm",names(AbundFD))] colnames(PDyr) <- substr(colnames(PDyr),1,4) # read in functional trait data functrt <- read.delim(file=paste(pathtofd,"EltonianTraits_BFISpecies.txt", sep = ""), header=T) # get list of species names spp_names <- functrt$BirdCode #reorder the rows of SNAB so that they match the order of the column #names of 'PD09' functrt <- functrt[(match(colnames(PDyr), functrt$BirdCode)),] #create new vector called 'CW'(conservation weight) using the breeding CCS score for each species CW <- functrt$CCSb #determine the column means for each species and log transform MDenyr <- colMeans(PDyr, na.rm=T) SDDenyr <- apply(PDyr,2,sd, na.rm=T) SDDenyr10 <- SDDenyr*1.10 meanvals <- cbind(MDenyr, SDDenyr, SDDenyr10) #determine sample size needed for resampling distribution N <- nrow(PDyr) #write resampFunc; use log normal distribution because density has to be greater than/equal to 0 # use code here to correctly determine parameters for lognormal distribution: https://msalganik.wordpress.com/2017/01/21/making-sense-of-the-rlnorm-function-in-r/ resampFunc <- function(x) { location <- log((x[1])^2 / sqrt((x[3])^2 + (x[1])^2)) shape <- sqrt(log(1 + ((x[3])^2 / (x[1])^2))) rlnorm(n=N, location, shape) } #draw all new density values for each species from a log normal distribution and put in dataframe; for species that had 0 value for mean and SD replace NaN with 0 SensitivityDens <- as.data.frame(apply(meanvals,1,resampFunc), by.row = TRUE) SensitivityDens[is.na(SensitivityDens)]<-NA SensitivityDens[is.na(SensitivityDens)] <- 0 summary(SensitivityDens) #create new dataframe called 'weightDen' by multiplying CW by each row #of SensitivityDen and original densities weightDen_samp <- data.frame(mapply(`*`,SensitivityDens,CW,SIMPLIFY=FALSE)) weightDen_orig <- data.frame(mapply(`*`,PDyr,CW,SIMPLIFY=FALSE)) #get the summed weighted density for each site(row)and put in 'SumWD' AbundFD$SumWD_samp <- rowSums(weightDen_samp) #calculate the variants of the final BFI by multiplying 'SumWD' #by the appropriate column in the Functional Diversity metric AbundFD$BFI_samp <- AbundFD$FunctShDiv * AbundFD$SumWD_samp #BFI calculated using Functional Shannon's Diversity for the sensitivity analysis AbundFD$BFI_logistic <- logistic(scale(AbundFD$BFI_samp)) # run 1000 bootstrap samples and calculate mean, sd, and overlap for each WDsampdf <- data.frame() for (i in 1:1000){ SensitivityDens <- as.data.frame(apply(meanvals,1,resampFunc), by.row = TRUE) SensitivityDens[is.na(SensitivityDens)]<-NA SensitivityDens[is.na(SensitivityDens)] <- 0 summary(SensitivityDens) weightDen_samp <- data.frame(mapply(`*`,SensitivityDens,CW,SIMPLIFY=FALSE)) AbundFD$SumWD_samp <- rowSums(weightDen_samp) AbundFD$BFI_samp <- AbundFD$FunctShDiv * AbundFD$SumWD_samp #BFI calculated using Functional Shannon's Diversity for the sensitivity analysis AbundFD$BFI_samp_logistic <- logistic(scale(AbundFD$BFI_samp)) WDov <- overlap(list(BFI_WD=AbundFD$BFI_samp_logistic[!(is.na(AbundFD$BFI_samp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), nbins = 1000, plot = FALSE, partial.plot = FALSE)$OV WDsampdf <- rbind(WDsampdf, data.frame(WDmean=mean(AbundFD$BFI_samp_logistic[!(is.na(AbundFD$BFI_samp_logistic))]), WDsd=sd(AbundFD$BFI_samp_logistic[!(is.na(AbundFD$BFI_samp_logistic))]), WDoverlap=WDov)) } write.csv(WDsampdf, paste(pathtofiles, "Sensitivity_Analysis_Results/WD_SensitityAnalysis_",year,".csv", sep="") ) WDsampdf <- read.csv(paste(pathtofiles, "Sensitivity_Analysis_Results/WD_SensitityAnalysis_",year,".csv", sep="") ) summary(WDsampdf) png(file=paste(pathtofiles, "Sensitivity_Analysis_Results/SensitivityAnalysis_Dens_",year,".png", sep=""), width=1500, height=1000, units="px" , res=300) final.plot2(list(BFI_WD=AbundFD$BFI_samp_logistic[!(is.na(AbundFD$BFI_samp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), facet.label="Avian density") dev.off() ############################################################################################ ####THIS SECTION WILL CALCULATE THE EFFECT OF CHANGING THE CONSERVATION WEIGHT #### OF EACH SPECIES; FINAL RESULT IS CALLED 'ES_CW' (EFFECT SIZE CONSERVATION WEIGHT) ############################################################################################ #determine sample size needed for resampling distribution and boundaries of uniform distribution N <- length(CW) maxCW <- max(CW) minCW <- min(CW) #draw new conservation weights from uniform distribution bounded by the min and max possible values for conservation weights CW_SA<-round(runif(N,minCW,maxCW),digits=0) #create new dataframe called 'weightDen' by multiplying CW and CW_SA by each row #of PD09 weightDen_orig <- data.frame(mapply(`*`,PDyr,CW,SIMPLIFY=FALSE)) weightDen_CWSA <- data.frame(mapply(`*`,PDyr,CW_SA,SIMPLIFY=FALSE)) #get the summed weighted density for each site(row)and put in 'SumWD' AbundFD$SumWD_CWsamp=rowSums(weightDen_CWSA) #calculate the variants of the final BFI by multiplying 'SumWD' #by the appropriate column in the Functional Diversity metric AbundFD$BFI_CWsamp <- AbundFD$FunctShDiv * AbundFD$SumWD_CWsamp #BFI calculated using Functional Shannon's Diversity for the sensitivity analysis # run 1000 bootstrap samples and calculate mean, sd, and overlap for each CWsampdf <- data.frame() for (i in 1:1000){ CW_SA<-round(runif(N,minCW,maxCW),digits=0) weightDen_CWSA=data.frame(mapply(`*`,PDyr,CW_SA,SIMPLIFY=FALSE)) AbundFD$SumWD_CWsamp=rowSums(weightDen_CWSA) AbundFD$BFI_CWsamp <- AbundFD$FunctShDiv * AbundFD$SumWD_CWsamp #BFI calculated using Functional Shannon's Diversity for the sensitivity analysis AbundFD$BFI_CWsamp_logistic <- logistic(scale(AbundFD$BFI_CWsamp)) CWov <- overlap(list(BFI_CW=AbundFD$BFI_CWsamp_logistic[!(is.na(AbundFD$BFI_CWsamp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), nbins = 1000, plot = FALSE, partial.plot = FALSE)$OV CWsampdf <- rbind(CWsampdf, data.frame(CWmean=mean(AbundFD$BFI_CWsamp_logistic[!(is.na(AbundFD$BFI_CWsamp_logistic))]), CWsd=sd(AbundFD$BFI_CWsamp_logistic[!(is.na(AbundFD$BFI_CWsamp_logistic))]), CWoverlap=CWov)) } write.csv(CWsampdf, paste(pathtofiles, "Sensitivity_Analysis_Results/CW_SensitityAnalysis_",year,".csv", sep="") ) summary(CWsampdf) png(file=paste(pathtofiles, "Sensitivity_Analysis_Results/SensitivityAnalysis_CW_",year,".png", sep=""), width=1500, height=1000, units="px" , res=300) final.plot2(list(BFI_CW=AbundFD$BFI_CWsamp_logistic[!(is.na(AbundFD$BFI_CWsamp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), facet.label="Conservation score") dev.off() ############################################################################################ ####THIS SECTION WILL CALCULATE THE EFFECT OF CHANGING THE FUNCTIONAL DIVERSITY SCORE #### OF EACH SITE; FINAL RESULT IS CALLED 'ES_FD' (EFFECT SIZE FUNCTIONAL DIVERSITY) ############################################################################################ #look at the scale, shape parameters of beta distribution from calculated functional diversity scores fitdist(AbundFD$FunctShDiv, "norm") N <- nrow(AbundFD) AbundFD$FD_SA<-rnorm(N,mean=mean(AbundFD$FunctShDiv),sd=(sd(AbundFD$FunctShDiv)*1.1)) #calculate the variants of the final BFI by multiplying 'SumWD' #by the appropriate column in the Functional Diversity metric AbundFD$BFI_FDsamp <- AbundFD$FD_SA * AbundFD$SumWD #BFI calculated using Functional Evenness for the sensitivity analysis # run 1000 bootstrap samples and calculate mean, sd, and overlap for each FDsampdf <- data.frame() for (i in 1:1000){ AbundFD$FD_SA<-rnorm(N,mean=mean(AbundFD$FunctShDiv),sd=(sd(AbundFD$FunctShDiv)*1.1)) AbundFD$BFI_FDsamp <- AbundFD$FD_SA * AbundFD$SumWD #BFI calculated using Functional Evenness for the sensitivity analysis AbundFD$BFI_FDsamp_logistic <- logistic(scale(AbundFD$BFI_FDsamp)) FDov <- overlap(list(BFI_FD=AbundFD$BFI_FDsamp_logistic[!(is.na(AbundFD$BFI_FDsamp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), nbins = 1000, plot = FALSE, partial.plot = FALSE)$OV FDsampdf <- rbind(FDsampdf, data.frame(FDmean=mean(AbundFD$BFI_FDsamp_logistic[!(is.na(AbundFD$BFI_FDsamp_logistic))]), FDsd=sd(AbundFD$BFI_FDsamp_logistic[!(is.na(AbundFD$BFI_FDsamp_logistic))]), FDoverlap=FDov)) } write.csv(FDsampdf, paste(pathtofiles, "Sensitivity_Analysis_Results/FD_SensitityAnalysis_",year,".csv", sep="") ) summary(FDsampdf) png(file=paste(pathtofiles, "Sensitivity_Analysis_Results/SensitivityAnalysis_FD_",year,".png", sep=""), width=1500, height=1000, units="px" , res=300) final.plot2(list(BFI_FD=AbundFD$BFI_FDsamp_logistic[!(is.na(AbundFD$BFI_FDsamp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), facet.label="Functional diversity") dev.off() } ############################################################################## # compile all years and produce plots # compile mean overlap values for each sample (6000 samples = 1000 x 6 years) from all years WDsampdf <- data.frame() CWsampdf <- data.frame() FDsampdf <- data.frame() for (year in 2009:2014){ temp <- read.csv(paste(pathtofiles, "Sensitivity_Analysis_Results/WD_SensitityAnalysis_",year,".csv", sep="") ) temp$X <- NULL WDsampdf <- rbind(WDsampdf, temp) temp <- read.csv(paste(pathtofiles, "Sensitivity_Analysis_Results/CW_SensitityAnalysis_",year,".csv", sep="") ) temp$X <- NULL CWsampdf <- rbind(CWsampdf, temp) temp <- read.csv(paste(pathtofiles, "Sensitivity_Analysis_Results/FD_SensitityAnalysis_",year,".csv", sep="") ) temp$X <- NULL FDsampdf <- rbind(FDsampdf, temp) } # read in and compile estimated BFIs AbundFD <- data.frame() for (year in 2009:2014){ temp <- read.delim(file=paste(pathtofd,"FinalDensities_FD_BFI_",year,".txt", sep=""),header=T) AbundFD <- rbind(AbundFD, temp) } AbundFD.all <- AbundFD #################### # species densities # get estimated mean overlap across all years summary(WDsampdf) # estimate mean densities for each species in 2014 (for plotting purposes) AbundFD <- read.delim(file=paste(pathtofd,"FinalDensities_FD_BFI_2014.txt", sep=""),header=T) PDyr <- AbundFD[,grep("*DensityKm",names(AbundFD))] colnames(PDyr) <- substr(colnames(PDyr),1,4) # read in functional trait data functrt <- read.delim(file=paste(pathtofd,"Priority Species Info_Eltonian1.0.txt", sep = ""),header=T) functrt <- functrt[which(functrt$BirdCode %in% speclist),] # get list of species names spp_names=functrt$BirdCode #reorder the rows of SNAB so that they match the order of the column functrt <- functrt[(match(colnames(PDyr),functrt$BirdCode)),] #create new vector called 'CW'(conservation weight) using the breeding CCS score for each species CW <- functrt$CCSb #determine the column means for each species and log transform MDenyr <- colMeans(PDyr, na.rm=T) SDDenyr <- apply(PDyr,2,sd, na.rm=T) SDDenyr10 <- SDDenyr*1.10 meanvals <- cbind(MDenyr, SDDenyr, SDDenyr10) # resample 2014 mean densities to get simulated densities N <- nrow(PDyr) SensitivityDens <- as.data.frame(apply(meanvals,1,resampFunc), by.row = TRUE) SensitivityDens[is.na(SensitivityDens)]<-NA SensitivityDens[is.na(SensitivityDens)] <- 0 # multiply simulated densities by Conservation Weights weightDen_samp <- data.frame(mapply(`*`,SensitivityDens,CW,SIMPLIFY=FALSE)) AbundFD$SumWD_samp <- rowSums(weightDen_samp) # calculate BFI AbundFD$BFI_samp <- AbundFD$FunctShDiv * AbundFD$SumWD_samp AbundFD$BFI_samp_logistic <- logistic(scale(AbundFD$BFI_samp)) WDov <- overlap(list(BFI_WD=AbundFD$BFI_samp_logistic[!(is.na(AbundFD$BFI_samp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), nbins = 1000, plot = FALSE, partial.plot = FALSE)$OV # check that overlap is similar to Mean WDoverlap from all years WDov # compare to estimated mean overlap across all years summary(WDsampdf) # repeat above code as needed until WDov (2014 only) is similar to estimated mean overlap across all years #Median :0.8722 #Mean :0.8388 write.csv(AbundFD, paste0(pathtobfi, "SensitivityAnalysis_SummarizedDensities_ForFigure.csv")) png(file=paste(pathtobfi, "SensitivityAnalysis_Dens_All.png", sep=""), width=1500, height=1000, units="px" , res=300) final.plot2(list(BFI_WD=AbundFD$BFI_samp_logistic[!(is.na(AbundFD$BFI_samp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), facet.label="Avian density") dev.off() ################################################## ## conservation scores # get estimated mean overlap across all years summary(CWsampdf) # estimated mean CWoverlap: N <- length(CW) maxCW <- max(CW) minCW <- min(CW) #draw new conservation weights from uniform distribution bounded by the min and max possible values for conservation weights CW_SA<-round(runif(N,minCW,maxCW),digits=0) weightDen_CWSA <- data.frame(mapply(`*`,PDyr,CW_SA,SIMPLIFY=FALSE)) AbundFD$SumWD_CWsamp <- rowSums(weightDen_CWSA) AbundFD$BFI_CWsamp <- AbundFD$FunctShDiv * AbundFD$SumWD_CWsamp #BFI calculated using Functional Shannon's Diversity for the sensitivity analysis AbundFD$BFI_CWsamp_logistic <- logistic(scale(AbundFD$BFI_CWsamp)) CWov <- overlap(list(BFI_CW=AbundFD$BFI_CWsamp_logistic[!(is.na(AbundFD$BFI_CWsamp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), nbins = 1000, plot = FALSE, partial.plot = FALSE)$OV CWov # compare with estimated mean overlap across all years summary(CWsampdf) write.csv(AbundFD, paste0(pathtobfi, "SensitivityAnalysis_SummarizedConservationWeights_ForFigure.csv")) png(file=paste(pathtobfi, "SensitivityAnalysis_CW_All.png", sep=""), width=1500, height=1000, units="px" , res=300) final.plot2(list(BFI_CW=AbundFD$BFI_CWsamp_logistic[!(is.na(AbundFD$BFI_CWsamp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), facet.label="Conservation score") dev.off() ################################### ## functional diversity # need to use 2014 data for FD year <- 2012 AbundFD <- read.delim(file=paste(pathtobfi,"FinalDensities_FD_BFI_",year,".txt", sep=""),header=T) PDyr <- AbundFD[,grep("*DensityKm",names(AbundFD))] colnames(PDyr) <- substr(colnames(PDyr),1,4) #create new vector called 'CW'(conservation weight) using the breeding CCS score for each species CW <- functrt$CCSb #determine the column means for each species and log transform MDenyr <- colMeans(PDyr, na.rm=T) SDDenyr <- apply(PDyr,2,sd, na.rm=T) SDDenyr10 <- SDDenyr*1.10 meanvals <- cbind(MDenyr, SDDenyr, SDDenyr10) summary(FDsampdf) AbundFD$FD_SA<-rnorm(nrow(AbundFD),mean=mean(AbundFD$FunctShDiv),sd=(sd(AbundFD$FunctShDiv)*1.1)) AbundFD$BFI_FDsamp <- AbundFD$FD_SA * AbundFD$SumWD #BFI calculated using Functional Evenness for the sensitivity analysis AbundFD$BFI_FDsamp_logistic <- logistic(scale(AbundFD$BFI_FDsamp)) FDov <- overlap(list(BFI_FD=AbundFD$BFI_FDsamp_logistic[!(is.na(AbundFD$BFI_FDsamp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), nbins = 1000, plot = FALSE, partial.plot = FALSE)$OV FDov # compare with estimated mean overlap summary(FDsampdf) write.csv(AbundFD, paste0(pathtobfi, "SensitivityAnalysis_SummarizedFunctionalDiversity_ForFigure.csv")) png(file=paste(pathtobfi, "SensitivityAnalysis_FD_All.png", sep=""), width=1500, height=1000, units="px" , res=300) final.plot2(list(BFI_FD=AbundFD$BFI_FDsamp_logistic[!(is.na(AbundFD$BFI_FDsamp_logistic))], BFI=AbundFD$BFIlogistic[!(is.na(AbundFD$BFIlogistic))]), facet.label="Functional diversity") dev.off()