# CONDOR-20-044 - A large-scale experiment demonstrates line marking reduces power line # collision mortality for large terrestrial birds, but not bustards, in the Karoo, South Africa # Jessica Shaw shawmjessica@gmail.com # October 2020 # R code for the analysis: # 1. Code for testing effectiveness of markers # Dataset 'MarkingOct2020' is the number of collisions per section check in the # reduced experiment for the 4 bird groups, with the number of spans in each section check # 2. Code for generating collision rates of the top 10 species pre-marking # Dataset 'BootstrapOct2020' is the number of collisions per span check by top 10 most commonly recovered species #1. Testing marking effectiveness library(lme4) library(MASS) hydraC<-MarkingOct2020 hydraC <- hydraC[1:767,] hydraC$Block<-as.factor(hydraC$Block) hydraC$ID<-as.factor(hydraC$ID) hydraC$Experiment <- factor(hydraC$Experiment) hydraC$TreatSimple <- factor(hydraC$TreatSimple) #DOES LINE MARKING WORK OVERALL? #Blue Cranes bcC1.nb <- glmer.nb(BC~TreatSimple+Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000))) bcC2.nb <- glmer.nb(BC~TreatSimple*Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 100000))) qqnorm(resid(bcC2.nb)) anova(bcC1.nb,bcC2.nb) summary(bcC2.nb) #Ludwig's Bustards lbC1.nb <- glmer.nb(LB~TreatSimple+Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000))) lbC2.nb <- glmer.nb(LB~TreatSimple*Experiment+offset(log(count_of_spans))+(1|Block),data=hydraC, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1000000000))) qqnorm(resid(lbC2.nb)) anova(lbC1.nb,lbC2.nb) summary(lbC2.nb) #Total bustards tbC1.nb <- glmer.nb(TotalBustards~TreatSimple+Experiment+offset(log(count_of_spans))+(1|Block),data=hydraC, glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000))) tbC2.nb <- glmer.nb(TotalBustards~TreatSimple*Experiment+offset(log(count_of_spans))+(1|Block),data=hydraC, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1000000000))) qqnorm(resid(tbC2.nb)) anova(tbC1.nb,tbC2.nb) summary(tbC2.nb) #Total large birds tlC1.nb <- glmer.nb(TotalLarge~TreatSimple+Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000))) tlC2.nb <- glmer.nb(TotalLarge~TreatSimple*Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000))) qqnorm(resid(tlC2.nb)) anova(tlC1.nb,tlC2.nb) summary(tlC2.nb) #For significant models, then calculate the percentage decrease in collisions on marked spans and #estimate confidence intervals with simulations from those models #BC (exp(fixef(bcC2.nb)[1])-exp(fixef(bcC2.nb)[1]+fixef(bcC2.nb)[2]))/exp(fixef(bcC2.nb)[1]) tbcC20 <- rnorm(1000,-3.8848,0.6660) tbcC21 <- rnorm(1000,-2.5567,0.5461) tbcC2_res <- (exp(tbcC20)-exp(tbcC20+tbcC21))/exp(tbcC20) summary(tbcC2_res) quantile(tbcC2_res,c(.025,.5,.975)) #TL (exp(fixef(tlC2.nb)[1])-exp(fixef(tlC2.nb)[1]+fixef(tlC2.nb)[2]))/exp(fixef(tlC2.nb)[1]) ttlC20 <- rnorm(1000,-2.8188,0.2622) ttlC21 <- rnorm(1000,-0.7173,0.2185) ttlC2_res <- (exp(ttlC20)-exp(ttlC20+ttlC21))/exp(ttlC20) summary(ttlC2_res) quantile(ttlC2_res,c(.025,.5,.975)) #WHICH MARKER IS BETTER? #BC bcC3.nb <- glmer.nb(BC~Treatment+Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 10000000))) bcC4.nb <- glmer.nb(BC~Treatment*Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000))) qqnorm(resid(bcC4.nb)) anova(bcC3.nb,bcC4.nb) summary(bcC4.nb) #LB lbC3.nb <- glmer.nb(LB~Treatment+Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1000000000))) lbC4.nb <- glmer.nb(LB~Treatment*Experiment+offset(log(count_of_spans))+(1|Block),data=hydraC, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1000000000))) qqnorm(resid(lbC4.nb)) anova(lbC3.nb,lbC4.nb) summary(lbC4.nb) #TB tbC3.nb <- glmer.nb(TotalBustards~Treatment+Experiment+offset(log(count_of_spans))+(1|Block),data=hydraC, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1000000000))) tbC4.nb <- glmer.nb(TotalBustards~Treatment*Experiment+offset(log(count_of_spans))+(1|Block),data=hydraC, glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1000000000))) qqnorm(resid(tbC4.nb)) anova(tbC3.nb,tbC4.nb) summary(tbC4.nb) #TL tlC3.nb <- glmer.nb(TotalLarge~Treatment+Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000))) tlC4.nb <- glmer.nb(TotalLarge~Treatment*Experiment+offset(log(count_of_spans))+(1|Block), data=hydraC, glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 10000000))) qqnorm(resid(tlC4.nb)) anova(tlC3.nb,tlC4.nb) summary(tlC4.nb) #For significant models, then calculate the percentage decrease in collisions on marked spans and #estimate confidence intervals with simulations from those models #BCs #Percentage decrease with BFD (exp(fixef(bcC4.nb)[1])-exp(fixef(bcC4.nb)[1]+fixef(bcC4.nb)[2]))/exp(fixef(bcC4.nb)[1]) tbcC40 <- rnorm(1000,-3.8766,0.6613) tbcC41 <- rnorm(1000,-2.2937,0.6256) tbcC4_res <- (exp(tbcC40)-exp(tbcC40+tbcC41))/exp(tbcC40) summary(tbcC4_res) quantile(tbcC4_res,c(.025,.5,.975)) #Percentage decrease with flapper (exp(fixef(bcC4.nb)[1])-exp(fixef(bcC4.nb)[1]+fixef(bcC4.nb)[3]))/exp(fixef(bcC4.nb)[1]) tbcC40 <- rnorm(1000,-3.8766,0.6613) tbcC42 <- rnorm(1000,-3.0802,1.0360) tbcC4b_res <- (exp(tbcC40)-exp(tbcC40+tbcC42))/exp(tbcC40) summary(tbcC4b_res) quantile(tbcC4b_res,c(.025,.5,.975)) #2. COLLISION RATES OF TOP 10 SPECIES PRE-MARKING #Bootstrapping pre-marking data for the top 10 most commonly recovered species at De Aar #Bootstrapping hierarchical and correlated data library(boot) BootstrapOct2020$Year<-as.factor(BootstrapOct2020$Year) BootstrapOct2020$ID<-as.factor(BootstrapOct2020$ID) resample<-function(dat,cluster,replace){ if(nrow(dat)==1||all(replace==FALSE)) return(dat) cls<-sample(unique(dat[[cluster[1]]]),replace=replace[1]) sub<-lapply(cls,function(b) subset(dat,dat[[cluster[1]]]==b)) if(length(cluster)>1) sub<-lapply(sub,resample,cluster=cluster[-1],replace=replace[-1]) do.call(rbind,sub) } allfactors<-c("Year","Season","Powerline","ID") #Ludwig's Bustard system.time(lb10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$LB))) quantile(lb10,c(0.025,0.5,0.975)) mean(lb10) #Blue Crane system.time(bc10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$BC))) quantile(bc10,c(0.025,0.5,0.975)) mean(bc10) #Northern Black Korhaan system.time(nbk10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$NBK))) quantile(nbk10,c(0.025,0.5,0.975)) mean(nbk10) #White Stork system.time(ws10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$WS))) quantile(ws10,c(0.025,0.5,0.975)) mean(ws10) #Secretarybird system.time(sec10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$Sec))) quantile(sec10,c(0.025,0.5,0.975)) mean(sec10) #Karoo Korhaan system.time(kk10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$KK))) quantile(kk10,c(0.025,0.5,0.975)) mean(kk10) #Pied Crow system.time(pc10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$PC))) quantile(pc10,c(0.025,0.5,0.975)) mean(pc10) #Blue Korhaan system.time(bk10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$BK))) quantile(bk10,c(0.025,0.5,0.975)) mean(bk10) #Spur-winged Goose system.time(swg10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$SWG))) quantile(swg10,c(0.025,0.5,0.975)) mean(swg10) #Kori Bustard system.time(kb10<-replicate(500,mean(resample(BootstrapOct2020,allfactors,c(T,T,T,T))$KB))) quantile(kb10,c(0.025,0.5,0.975)) mean(kb10)