Abstract

During the late 20th and early 21st centuries, Federal, state, and private organizations were motivated to advance the recovery of the bald eagle (Haliaeetus leucocephalus) in the Contiguous U.S. through several reintroduction methods, including hacking, an adapted falconry technique involving the transplanting, conditioning, and releasing nestlings. Beginning in New York state in 1976 and continuing in other states over subsequent decades, wildlife managers and conservationists sourced hundreds of nestlings from breeding populations in Alaska, Canada, and the portions of the Upper Midwestern USA and hacked them into local populations. Although subsequent monitoring confirmed the successful nesting of hacked birds and recent genetic analysis has suggested admixture related to hacking, the contributions of these efforts to observed eagle population growth rates have not been tested. Here, we used empirical data on bald eagle demographic rates, to simulate population growth of eagles among states to compare management scenarios with and without hacking with observed population data. We modeled state-level hacking programs that incorporated published timeframes and numbers of released birds, as well as accounting for annually-reported numbers of nesting territories. Overall, we found mixed support among states for the efficacy of bald eagle hacking. For example in Georgia, Massachusetts, and New York, hacking showed no contribution to observed population growth. Alternatively, simulation results in Indiana, New Jersey, and Vermont suggested that hacking was important for supporting observed population growth rates. Whereas, in Oklahoma and Pennsylvania, simulation results were less clear as to the contribution of hacking to population recovery with some state breeding population trends closely aligning with hacked simulations and others showing little to no difference from counterfactual simulations without hacking. In multiple instances, hacking appeared to advance achievement of state recovery goals by several years. Our results, which in some instances run counter to the generally-accepted thought that hacking has benefited eagle population recovery highlights the importance for the continual examination of efficacy of recovery methods used to support eagle populations, and species recovery programs in general.


Introduction

During the late 20th and early 21st centuries, Federal, state, and private organizations were motivated to advance the recovery of the bald eagle (Haliaeetus leucocephalus) in the Contiguous U.S. through several reintroduction methods, including hacking, an adapted falconry technique involving the transplanting, conditioning, and releasing nestlings (Sherrod 1982). Beginning in New York state in 1976 and continuing in other states over subsequent decades, wildlife managers and conservationists sourced hundreds of nestlings from breeding populations in Alaska, the Upper Midwestern USA, and Canada, as well as from captive-reared populations, and hacked them into local populations. Although subsequent monitoring confirmed the successful nesting of hacked birds (Sorenson et al. 2017, NYSDEC 2016, Hatcher 1991), and recent genetic analysis has suggested admixture related to hacking (Judkins et al. 2020), the contributions of these efforts to subsequent population recovery have not been tested. Here, we used empirical data on bald eagle demographic rates to simulate population growth of eagles among states to compare management scenarios with and without hacking with observed population data.


Methods

We constructed a simulation function in R (R Core Team 2020) to project population growth under management regimes with and without hacking. Each simulation began with a founder population based on the historic breeding population at start of a state’s hacking program. Within this founder population, we randomly assigned individual birds ages between four (age at first breeding) and 25 (maximum lifespan; USFWS 2016). We then divided the simulation population into hatch-year (HY) and after hatch-year (AHY) birds and applied annual mortality to each class according to known survival rates (USFWS 2016). In simulations with hacking, we then added the number of HY birds reported to have been hacked into the state population for that historic year and applied HY annual mortality rates to these birds. Next, for both simulations with and without hacking, we immigrated or emigrated four-year-old birds to or from the population by sampling a distribution derived from the difference between empirical population numbers in Connecticut and simulated population numbers for the state based on known mortality and productivity rates (assuming Net Immigration = Population + Births - Deaths). We then aged all birds by one year and moved to the next simulation year with this aged population, repeating this process until the end of the specified study period. To account for random variation, we then repeated this entire process 1,000 times for both simulations with and without hacking and calculated means and standard deviations for annual breeding populations under the two scenarios. Finally, we compared the simulated hacked and non-hacked population numbers against the state’s historic numbers.


Load required packages


Load state-level eagle population data

ny_pop <- read.csv("./Data/ny_pop.csv")
nj_pop <- read.csv("./Data/nj_pop.csv")
pa_pop <- read.csv("./Data/pa_pop.csv")
ok_pop <- read.csv("./Data/ok_pop.csv")
vt_pop <- read.csv("./Data/vt_pop.csv")
ma_pop <- read.csv("./Data/ma_pop.csv")
in_pop <- read.csv("./Data/in_pop.csv")
ga_pop <- read.csv("./Data/ga_pop.csv")
tn_pop <- read.csv("./Data/tn_pop.csv")
c_ca_pop <- read.csv("./Data/c_ca_pop.csv")
mo_pop <- read.csv("./Data/mo_pop.csv")
nc_pop <- read.csv("./Data/nc_pop.csv")

# The following states are 'unhacked' populations used for 
# gut checking accuracy of simulation's growth rates
ct_pop <- read.csv("./Data/ct_pop.csv")
de_pop <- read.csv("./Data/de_pop.csv")
nh_pop <- read.csv("./Data/nh_pop.csv")
va_pop <- read.csv("./Data/va_pop.csv")

#Immature empirical distribution used in model
load("./Data/imm_dist.RData")


Model description

List of model parameters:

  1. p(Adult survival) ~ rnorm(mean=0.9089, sd=0.0054)
  2. p(Hatch-year survival) ~ rnorm(mean=0.0=86, sd=0.0395)
  3. p(Hacked bird survival) ~ rnorm(mean=0.0=86, sd=0.0395)
  4. Senescence Age ~ runif(min=25,max=30)
  5. Recruitment Age ~ runif(min=4,max=5)
  6. Number of birds hacked per year ~ from states
  7. Num of simulations = 1,000
  8. Error (general) term ~ runif(min=0.9,max=1.1)
  9. Error hacked bird survival ~ runif(min=0.7,max=1)


Hacking simulation model

#load required functions to generate random model parameters
#random error term
randError<-function(){
  new.rand<-runif(n=1, min=0.9, max=1.1)
  return(new.rand)
}

#random age of adult bird senescence
randSenesce<-function(){
  new.rand<-runif(n=1, min=25, max=30)
  return(new.rand)
}

#random age of recruitment
randRecruit<-function(){
  new.rand<-runif(n=1, min=4, max=5)
  return(new.rand)
}

#modify hacked birds survival
modHack<-function(){
  new.rand<-runif(n=1, min=0.7, max=1)
  return(new.rand)
}


# create function to simulate management scenarios with & without hacking
hack_sim <- function(found, # number of breeding pairs in found population
                     start, # start year for simulation YYYY
                     end, # end year for simulation YYYY
                     nsims, # number of simulations
                     brd_age, # age of first breeding
                     hack_yr, # vector of years in which birds are hacked
                     hack_n, # vector of corresponding number of hacked birds
                     state_dat = FALSE, # csv containing observed state data
                     pop_goal = FALSE, # number of desired breeding pairs
                     study_area){ # name of study area
  
  

  # matrix created from combined vectors of hacking years and hacking numbers
  hack_mat <- matrix(c(hack_yr, hack_n), ncol = 2)
  
  # empty matrix to receive breeding population numbers 
  # for simulation w/ hacking
  sim_hack <- matrix(data = NA, nrow = nsims, ncol = end-start+1,
                   dimnames = list(1:nsims,start:end))
  
  # empty matrix to receive breeding population numbers 
  # for simulation w/out hacking
  sim_non <- matrix(data = NA, nrow = nsims, ncol = end-start+1,
                    dimnames = list(1:nsims,start:end)) 
  
  # simulation across study period of
  # two populations: hacked and non
  for(z in 1:nsims){
    
    #initialize random Senecence age for each iteration
    randSenesce_mod<-randSenesce()
    
    #initialize age at recruitment
    randRecruit_mod<-randRecruit()
    
    
    # create founding population for both simulations
    # age randomly selected from between age of first breeding and twenty five
    pop <- list(runif(found*2, min = randRecruit_mod, max = randSenesce_mod),runif(found*2, min = randRecruit_mod, max = randSenesce_mod))
    
    # list-vectors to receive simulated values
    pop_hy <- vector("list",2)
    pop_ahy <-  vector("list",2)
    pop_fl <-  vector("list",2)
    pop_ad <- vector("list",2)
    sim_brd <- vector("list",2)
    sim_pop <- vector("list",2)
    pop_4y <- vector("list",2)
    pop_n4y <- vector("list",2)
    pop_imm <- vector("list",2)
    
    # create sub-loops to simulate both hacked and non-hack populations 
    # through the study period
    for(y in 1:2){
      for(i in 1:(end-start+1)){
        # subset population to hatch year and after hatch birds
        pop_hy[[y]] <- pop[[y]][pop[[y]] < 2]
        pop_ahy[[y]] <- pop[[y]][pop[[y]] >= 2]
        
        # apply mortality to both age classes (options: mean=0.71, sd=0.212 or mean=0.86, sd=0.0395)
        
        #HY birds
        tryCatch({
          pop_hy[[y]] <- sample(pop_hy[[y]], length(pop_hy[[y]])*rnorm(1, mean = 0.86, sd = 0.0395))
        }, error=function(cond2=NULL){
          while(is.null(cond2)){
            pop_hy[[y]] <- sample(pop_hy[[y]], length(pop_hy[[y]])*rnorm(1, mean = 0.86, sd = 0.0395))
          }
        })
          
        #AHY birds (options mean=0.9098, sd=0.0054 or mean=0.91, sd=0.014- double check)
        tryCatch({
          pop_ahy[[y]] <- sample(pop_ahy[[y]], length(pop_ahy[[y]])*rnorm(1, mean = 0.9098, sd = 0.0054))
        }, error=function(cond2=NULL){
          while(is.null(cond2)){
            pop_ahy[[y]] <- sample(pop_ahy[[y]], length(pop_ahy[[y]])*rnorm(1, mean = 0.9098, sd = 0.0054))
          }
        })
        
          
        # eliminate birds that are senesced or older
        pop_ahy[[y]] <- pop_ahy[[y]][pop_ahy[[y]] <=randSenesce_mod] 
        
        # simulate breeding
        pop_ad[[y]] <- pop_ahy[[y]][pop_ahy[[y]]>=randRecruit_mod] # potential breeders
        
        # create 'fledglings' by taking number of breeding pairs
        # and multiplying by success rate and random
        
        # sampling of productivity distribution as taken from USFWS 2020 BAEA Population Update
        tryCatch({
          pop_fl[[y]] <- rep(0,floor(floor(length(pop_ad[[y]])/2)*rlnorm(1, meanlog = 0.12, sdlog = 0.21)))
        }, error=function(cond2=NULL){
          while(is.null(cond2)){
            pop_fl[[y]] <- rep(0,floor(floor(length(pop_ad[[y]])/2)*rlnorm(1, meanlog = 0.12, sdlog = 0.21)))
          }
        })
        
        # combine fledglings, hatch-year birds, and after hatch-year birds
        # and age all birds one year
        pop[[y]] <- c(pop_fl[[y]],pop_hy[[y]],pop_ahy[[y]])+1
        
        # simulate hacking for defined hacking period
        if(y==2 & i %in% (hack_mat[,1]-start)){
          # hack in number of birds specified in hacking matrix
          hacks <- rep(1,hack_mat[which(hack_mat[,1]-start==i),2])
          # apply mortality to hacked birds
          hacks <- sample(hacks, length(hacks)*rnorm(1, mean = 0.86, sd = 0.0395)) # mortality from NY hacking
          # add hacked birds to population
          pop[[y]] <- c(pop[[y]],hacks)
        }
        
        # simulate immigration/emigration
        # immigration/emigration is only applied to 4th-year birds in this simulation
        imm_em <- sample(imm_dist,1)
        if(imm_em<0){
          pop_4y[[y]] <- pop[[y]][pop[[y]]==4]
          pop_n4y[[y]] <- pop[[y]][pop[[y]]!=4]
          if(abs(imm_em)>length(pop_4y[[y]])){
            pop[[y]] <- pop_n4y[[y]]
          }else{
            pop_4y[[y]] <- sample(pop_4y[[y]],(length(pop_4y[[y]])-abs(imm_em)))
            pop[[y]] <- c(pop_n4y[[y]],pop_4y[[y]])
          }

        }else{
          pop_imm[[y]] <- rep(4,imm_em)
          pop[[y]] <- c(pop[[y]],pop_imm[[y]])
        }
        
        # record the number of breeding pairs at the end of each simulation
        sim_brd[[y]][i] <- floor(length(pop[[y]][pop[[y]] >= brd_age])/2)
        
    }
  }
    sim_hack[z,] <- sim_brd[[2]]*randError()
    sim_non[z,] <- sim_brd[[1]]*randError()
    
  }
  
  # calculate means and confidence intervals across simulated years
  # for hacked and non-hacked populations
  hack <- data.frame(year=as.numeric(colnames(sim_hack)), 
                     mean = colMeans(sim_hack),
                     lwr = apply(sim_hack,2,function(x)confint(lm(x~1),level = 0.95)[,1]),
                     upr = apply(sim_hack,2,function(x)confint(lm(x~1),level = 0.95)[,2]))
  
    no_hack <- data.frame(year=as.numeric(colnames(sim_non)), 
                        mean = colMeans(sim_non),
                        lwr =apply(sim_non,2,function(x)confint(lm(x~1),level = 0.95)[,1]),
                        upr = apply(sim_non,2,function(x)confint(lm(x~1),level = 0.95)[,2]))
    if(pop_goal != FALSE && state_dat != FALSE){
      obs_goal_year <- min(state_dat[which(state_dat$brd>pop_goal),]$year)
        
      no_hack_goal_year <- min(no_hack[which(no_hack$mean>pop_goal),]$year)
      
      print(paste0("With hacking, the population goal was achieved by ",obs_goal_year,". ",
            "Without hacking, the population goal would have been achieved by ",no_hack_goal_year,
            ", representing a delay of ", no_hack_goal_year - obs_goal_year, " years."))
    }
  # generate plot with hacked and non-hacked breeding populations with
  # confidence bands and actual historic breeding populations
  print(ggplot(data=hack, aes(x=year, y=mean, group = 1))+
    geom_line()+
    geom_ribbon(data = hack, aes(ymin=lwr,ymax=upr), linetype = 2, 
                color = "blue", fill = "blue", alpha = 0.4)+
    scale_x_continuous(breaks = seq(min(hack$year),max(hack$year), by = 5))+
    geom_line(data=no_hack, aes(x=year, y=mean, group = 1))+
    geom_ribbon(data = no_hack, aes(ymin=lwr,ymax=upr), linetype = 2, 
                  color = "red", fill = "red", alpha = 0.4)+
    geom_point(data = state_dat, aes(x = year, y = brd), na.rm = TRUE)+
    {if(pop_goal != FALSE)geom_hline(yintercept = pop_goal)} +  
    labs(title = study_area, x = "Year", y = "Number of Breeding Pairs")+
    scale_fill_manual(name = "Simulation", values = c(blue='blue',red='red'),
                      labels = c("Hack", "No hack"), guide = "legend"))
    
  output <- list(sim_hack,sim_non, hack, no_hack, state_dat,study_area)
}


Iterate model over states

#define NumSims
NumSims<-1000

ny_hack <- hack_sim(found = 1, start = 1970, end = 2009, nsims = NumSims,brd_age = 4,
                    hack_yr = c(1976:1988), hack_n = c(5,4,5,4,5,22,22,22,22,22,22,22,22),
                    state_dat = ny_pop, pop_goal = 40, study_area = "New York")
## [1] "With hacking, the population goal was achieved by 1999. Without hacking, the population goal would have been achieved by 2004, representing a delay of 5 years."

nj_hack <- hack_sim(found = 1, start = 1987, end = 2020, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1983:1990), hack_n = rep(8,8), 
                    state_dat = nj_pop, study_area = "New Jersey")

pa_hack <- hack_sim(found = 2, start = 1980, end = 2015, nsims = NumSims ,brd_age = 4,
                    hack_yr = c(1983:1988), hack_n = rep(14,6),
                    state_dat = pa_pop, pop_goal = 150, study_area = "Pennsylvania")
## [1] "With hacking, the population goal was achieved by 2008. Without hacking, the population goal would have been achieved by Inf, representing a delay of Inf years."

ok_hack <- hack_sim(found = 1, start = 1984, end = 2005, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1984:1990), hack_n = rep(12,7),
                    state_dat = ok_pop, pop_goal = 10, study_area = "Oklahoma")
## [1] "With hacking, the population goal was achieved by 1992. Without hacking, the population goal would have been achieved by 2005, representing a delay of 13 years."

vt_hack <- hack_sim(found = 1, start = 2002, end = 2035, nsims = NumSims, brd_age = 4,
                    hack_yr = c(2004:2006), hack_n = c(8,11,10),
                    state_dat = vt_pop, pop_goal = 28, study_area = "Vermont")
## [1] "With hacking, the population goal was achieved by 2019. Without hacking, the population goal would have been achieved by 2032, representing a delay of 13 years."

ma_hack <- hack_sim(found = 1, start = 1981, end = 2018, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1982:1988), hack_n = rep(6,7),
                    state_dat = ma_pop, study_area = "Massachusetts")

in_hack <- hack_sim(found = 0, start = 1985, end = 2020, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1985:1989), hack_n = c(14,14,15,15,15),
                    state_dat = in_pop, pop_goal = 50, study_area = "Indiana")
## [1] "With hacking, the population goal was achieved by 2004. Without hacking, the population goal would have been achieved by Inf, representing a delay of Inf years."

ga_hack <- hack_sim(found = 1, start = 1978, end = 2020, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1979:1995), hack_n = c(rep(5,13),6,6,6,6),
                    state_dat = ga_pop, study_area = "Georgia")

tn_hack <- hack_sim(found = 0, start = 1980, end = 2015, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1980:2012), hack_n = c(2,6,4,5,10,11,16,22,19,20,31,37,34,13,15,8,11,5,5,6,4,0,3,2,5,6,2,10,10,4,4,7,7),
                    state_dat = tn_pop, study_area = "Tennessee")

c_ca_hack <- hack_sim(found = 0, start = 1986, end = 2015, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1986:1994), hack_n = c(1,4,12,10,10,12,0,12,5),
                    state_dat = c_ca_pop, pop_goal = 4, study_area = "California")
## [1] "With hacking, the population goal was achieved by 1997. Without hacking, the population goal would have been achieved by 1999, representing a delay of 2 years."

mo_hack <- hack_sim(found = 0, start = 1980, end = 2020, nsims = NumSims, brd_age = 4,
                      hack_yr = c(1981:1990), hack_n = c(7,6,7,6,7,6,15,7,6,7),
                      state_dat = mo_pop, study_area = "Missouri")

nc_hack <- hack_sim(found = 0, start = 1980, end = 2020, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1983:1988), hack_n = c(4,4,4,4,9,4),
                    state_dat = nc_pop, study_area = "North Carolina")

#model calibration testing
ct_hack <- hack_sim(found = 1, start = 1990, end = 2019, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1994:1996), hack_n = c(2,2,2),
                    state_dat = ct_pop, study_area = "Connecticut")

de_hack <- hack_sim(found = 6, start = 1990, end = 2020, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1994:1996), hack_n = c(2,2,2),
                    state_dat = de_pop, study_area = "Delaware")

nh_hack <- hack_sim(found = 1, start = 1990, end = 2020, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1994:1996), hack_n = c(2,2,2),
                    state_dat = nh_pop, study_area = "New Hampshire")

va_hack <- hack_sim(found = 106, start = 1990, 2020, nsims = NumSims, brd_age = 4,
                    hack_yr = c(1994:1996), hack_n = c(2,2,2),
                    state_dat = va_pop, study_area = "Virginia")


Function to compile simulation results and add cumulative distribution values

#function to clean up results
build.df<-function(simuResults){
  
  new.simu.results<-simuResults
  
  #make hack.df
  hack.df<-as.data.frame(new.simu.results[3])
  row.names(hack.df)<-NULL
  
  #add Simulation_Type column
  hack.df$Simulation_Type<-"Hacking"
  
  #add State column
  hack.df$State<-unique(unlist(new.simu.results[6]))
  
  #modify column names
  names(hack.df)<-c("Year","Pop_size","lwr","upr","Simulation_Type","State")
  
  
  #make non-hack.df
  no_hack.df<-as.data.frame(new.simu.results[4])
  row.names(no_hack.df)<-NULL
  
  #add Simulation_Type column
  no_hack.df$Simulation_Type<-"No_Hacking"
  
  #add State column
  no_hack.df$State<-unique(unlist(new.simu.results[6]))
  
  #modify column names
  names(no_hack.df)<-c("Year","Pop_size","lwr","upr","Simulation_Type","State")
  
  
  #pop.df
  pop.df<-as.data.frame(new.simu.results[5])
  
  #add lwr and upr
  pop.df$lwr<-NA
  pop.df$upr<-NA
  
  #add Simulation_Type
  pop.df$Simulation_Type<-"Observed"
  
  #add State
  pop.df$State<-unique(unlist(new.simu.results[6]))
  
  #modify names
  names(pop.df)<-c("Year","Pop_size","lwr","upr","Simulation_Type","State")
  
  ############################################
  #stack hacking, no-hacking, and observed to get greatest pop size
  all.df<-plyr::rbind.fill(pop.df, hack.df,no_hack.df)
  
  #get vector of max values per year
  years.max<-aggregate(Pop_size~Year, FUN=max, data=all.df)
  names(years.max)[names(years.max)=="Pop_size"]<-"Max_pop_size"

  
  #get cumulative distributions for each State and treatment
  
  #now merge each data.frame with years.max
  
  #observed population estimates
  pop.df.merge<-merge(pop.df, years.max, by="Year",all.x=TRUE)
  
  pop.df.merge$CDF<- mltools::empirical_cdf(pop.df.merge$Pop_size,ubounds=pop.df.merge$Max_pop_size)
   #clean up
  pop.df.merge$CDF$UpperBound<-NULL
  pop.df.merge$CDF$N.cum<-NULL
  pop.df.merge$CDF<-pop.df.merge$CDF$CDF

  #hacking simulation
   hack.df.merge<-merge(hack.df, years.max, by="Year",all.x=TRUE)
   
   hack.df.merge$CDF<-mltools::empirical_cdf(hack.df.merge$Pop_size, ubounds=hack.df.merge$Max_pop_size)
   #clean up
   hack.df.merge$CDF$UpperBound<-NULL
   hack.df.merge$CDF$N.cum<-NULL
   hack.df.merge$CDF<-hack.df.merge$CDF$CDF
   
   #no_hacking simulation
   no_hack.df.merge<-merge(no_hack.df, years.max, by="Year",all.x=TRUE)
   no_hack.df.merge$CDF<-mltools::empirical_cdf(no_hack.df.merge$Pop_size,ubounds=no_hack.df.merge$Max_pop_size)
   #clean up
   no_hack.df.merge$CDF$UpperBound<-NULL
   no_hack.df.merge$CDF$N.cum<-NULL
   no_hack.df.merge$CDF<-no_hack.df.merge$CDF$CDF
   

  ###########################################
  #now rbind all data.frames
  data.all<-rbind(pop.df.merge, hack.df.merge, no_hack.df.merge)
  

  return(data.all)
  
}


Compile all simulation results
#create list of state-level simulation results
simuResultsList<-list(ny_hack, nj_hack, pa_hack, ok_hack, vt_hack, ma_hack, in_hack, ga_hack, tn_hack, c_ca_hack,mo_hack,nc_hack)


#create empty list to gather results
save.results<-list()

for(i in 1:length(simuResultsList)){
  
  #print(i)
  
  new.sim<-simuResultsList[[i]]
  
  save.df<-build.df(new.sim)
  
  save.results<-rbind(save.results, save.df)
  
}

save.results.df<-as.data.frame(save.results)

#bring in other covariates 
myCovs<-readxl::read_excel("./Data/Data Sources.xlsx",sheet="Sheet1")

#fix columm types
myCovs$State<-as.factor(myCovs$State)
myCovs$Type<-as.factor(myCovs$Type)
myCovs$Source<-as.factor(myCovs$Source)

#merge with save.results.df
save.results.merge<-merge(save.results.df, myCovs, by="State",all.x=TRUE)

#save compiled simulation results
#write.csv(save.results.merge, file="./Data/all_simulation_results_compiled.csv",row.names=FALSE)

all.simu.data<-save.results.merge

#fix columns again
all.simu.data$State<-as.factor(all.simu.data$State)
all.simu.data$Type<-as.factor(all.simu.data$Type)
all.simu.data$Source<-as.factor(all.simu.data$Source)
all.simu.data$Simulation_Type<-factor(all.simu.data$Simulation_Type,levels=c("Observed","Hacking","No_Hacking"))

#subset data to include only 1980-2020
all.simu.data.sub<-subset(all.simu.data, c(Year > 1979 & Year < 2021))


Results

None of the states we examined showed population growth rates lower than what was predicted in simulations without hacking, suggesting our demographic rates and modeling assumptions accurately the range of population growth. In two states (MA, Central CA), the observed population growth matched what was predicted under simulations without hacking, and in two others (NY, TN), observed population growth fell between what was predicted under simulations with and without hacking. The observed population growth of three states (PA, OK, GA) closely tracked what was predicted in simulations with hacking. Finally, in five states (NJ, VT, IN, MS, NC), observed population growth exceeded the predicted growth of simulated hacked populations. Overall, historic hacking efforts appeared to improve the rate of recovery among state-level breeding bald eagle populations.


Fit linear mixed-effects model with Year as random effect

#fit some models (with Year as random effect)
mod.1<-lmer(Pop_size~Simulation_Type*State+Source+ (1|Year), data=all.simu.data.sub)

#look at model summary
model.summary<-summary(mod.1)

#make data.frame
summary.mod.df<-data.frame(model.summary$coefficients)

#add column of model parameter
summary.mod.df<-data.frame("Parameter"=row.names(summary.mod.df), summary.mod.df)

#remove now.names
row.names(summary.mod.df)<-NULL


panderOptions('table.alignment.default', function(df) ifelse(sapply(df, is.numeric), 'right', 'left'))
panderOptions("table.split.table", Inf)

pander(summary.mod.df)
Parameter Estimate Std..Error t.value
(Intercept) 14.17 16.88 0.8398
Simulation_TypeHacking 38.24 16.98 2.251
Simulation_TypeNo_Hacking -0.862 16.98 -0.05075
StateGeorgia 64.58 19.24 3.357
StateIndiana 49.4 17.83 2.77
StateMassachusetts 14.79 19.21 0.77
StateMissouri 77.56 20.67 3.752
StateNew Jersey 34.83 16.6 2.098
StateNew York 64.03 19.83 3.229
StateNorth Carolina 41.58 19.88 2.092
StateOklahoma 42.66 21.13 2.019
StatePennsylvania 61.61 19.17 3.213
StateTennessee 60.26 17.78 3.389
StateVermont -61.72 19.98 -3.089
Simulation_TypeHacking:StateGeorgia 26.6 23.98 1.11
Simulation_TypeNo_Hacking:StateGeorgia -51.96 23.98 -2.167
Simulation_TypeHacking:StateIndiana -22.22 23.13 -0.9608
Simulation_TypeNo_Hacking:StateIndiana -61.08 23.13 -2.641
Simulation_TypeHacking:StateMassachusetts 7.075 24.1 0.2935
Simulation_TypeNo_Hacking:StateMassachusetts -7.063 24.1 -0.293
Simulation_TypeHacking:StateMissouri -5.709 25.15 -0.227
Simulation_TypeNo_Hacking:StateMissouri -71.6 25.15 -2.847
Simulation_TypeHacking:StateNew Jersey -66.84 22.29 -2.999
Simulation_TypeNo_Hacking:StateNew Jersey -51.49 22.29 -2.31
Simulation_TypeHacking:StateNew York 90.73 25.12 3.612
Simulation_TypeNo_Hacking:StateNew York -19.8 25.12 -0.7884
Simulation_TypeHacking:StateNorth Carolina -40.24 24.5 -1.642
Simulation_TypeNo_Hacking:StateNorth Carolina -35.97 24.5 -1.468
Simulation_TypeHacking:StateOklahoma -27.65 26.9 -1.028
Simulation_TypeNo_Hacking:StateOklahoma -13.53 26.9 -0.5029
Simulation_TypeHacking:StatePennsylvania 2.955 24.19 0.1222
Simulation_TypeNo_Hacking:StatePennsylvania -41.03 24.19 -1.696
Simulation_TypeHacking:StateTennessee 92.47 23.03 4.016
Simulation_TypeNo_Hacking:StateTennessee -44.88 23.03 -1.949
Simulation_TypeHacking:StateVermont -39.05 26.42 -1.478
Simulation_TypeNo_Hacking:StateVermont -7.676 26.42 -0.2906


Examine post-hoc pair-wise comparisons among Observed, Hacking, and No_Hacking scenearios
#use lsmeans to get pair-wise comparisons
Simulation_Type_compare<-lsmeans(mod.1, pairwise~Simulation_Type, adjust="tukey")

#make data.frame
Simulation_Type_Compare.df<-data.frame(Simulation_Type_compare$contrasts)

#plot
pander(Simulation_Type_Compare.df)
contrast estimate SE df t.ratio p.value
Observed - Hacking -40.97 5.596 967.7 -7.321 1.618e-12
Observed - No_Hacking 30.31 5.596 967.7 5.417 2.294e-07
Hacking - No_Hacking 71.28 4.65 965 15.33 0


Summary table of scenario mean, var, sd, se, CV, and lower and upper 95% CIs

#get means
simu.type.summary<-summaryFunction(dataIn=all.simu.data.sub, factor="Simulation_Type",response = "Pop_size")
simu.type.summary$Simulation_Type<-factor(simu.type.summary$Simulation_Type, levels=c("Observed","Hacking","No_Hacking"))

pander(simu.type.summary)
Simulation_Type n mean var SD SE CV lwr upr
Observed 320 42.95 3672 60.59 3.387 1.411 36.31 49.59
Hacking 404 94.12 17359 131.8 6.555 1.4 81.28 107
No_Hacking 404 15.17 354 18.81 0.9361 1.24 13.34 17.01


Figure 1. Boxplot showing comparisons of number of breeding eagle pairs among scenarios
#define myColors
myColors<-c("gray50","royalblue3","tomato")

#boxplot of Simulation_Type
simulationTypeBoxPlot<-ggplot(data=all.simu.data.sub, aes(x=Simulation_Type, y=Pop_size))+
  geom_boxplot(aes(fill=Simulation_Type, color=Simulation_Type),alpha=0.7)+
  scale_fill_manual(values=myColors)+
  scale_color_manual(values=myColors)+
  theme(panel.background = element_rect(fill="white",color="black"),
        panel.border = element_rect(fill="transparent",color="black"),
        legend.position="right")+
  guides(fill=guide_legend(title="Scenario"),color=FALSE)+
  labs(x="Scenario",y="Number of Breeding Pairs")+
  ylim(c(0,750))+
  geom_text(data=simu.type.summary, aes(x=Simulation_Type,y=c(720, 720, 720)), label=c("b","a","c"))


#save plot
# ggsave(simulationTypeBoxPlot, file="./Figures/Eagle_Simulation_Type_Boxplot.png",width=5, height=4,
#        dpi=1200)

print(simulationTypeBoxPlot)


Figure 2. State-level comparison of number of breeding eagle pairs among Observed, Hacking, and No_Hacking model scenarios.
#order States by maximum pop size in 2021
states.df<-subset(all.simu.data, c(Simulation_Type=="Observed" & Year<2022))

#get max
states.max<-aggregate(Pop_size~State, FUN=max, data=states.df)

#now order
states.order<-states.max[order(states.max$Pop_size,decreasing=TRUE),]

#set factor level order
all.simu.data.sub$State<-factor(all.simu.data.sub$State, levels=unique(states.order$State))


#remove NAs from data
all.simu.data.sub.no.NAs<-all.simu.data.sub[!is.na(all.simu.data.sub$Pop_size),]

#make Year an integer
all.simu.data.sub.no.NAs$Year<-as.integer(as.character(all.simu.data.sub.no.NAs$Year))

eagleAreaPlots<-ggplot(data=all.simu.data.sub.no.NAs, aes(x=Year, y=Pop_size,group=Simulation_Type))+
  geom_area(aes(fill=Simulation_Type,color=Simulation_Type),position=position_dodge(width=0.5),alpha=0.5)+
  #geom_path(aes(color=Simulation_Type))+
  scale_fill_manual(values=myColors)+
  scale_color_manual(values=myColors)+
  theme(panel.background = element_rect(fill="white",color="black"),
        panel.border = element_rect(fill="transparent",color="black"),
        #legend.key.size= unit(1, 'cm'),
        legend.title=element_text(size=10),
        legend.text=element_text(size=8),
        axis.text=element_text(size=8),
        axis.title=element_text(size=10))+
  labs(x="Year",y="Number of Breeding Pairs")
#eagleAreaPlots

#now facet
eagleAreaPlot_facet<-eagleAreaPlots+facet_wrap(~State, scales="free_y")+
  theme(legend.position = "top",
        strip.background = element_rect(fill="gray50"),
        strip.text=element_text(size=14, color="white"))

# ggsave(eagleAreaPlot_facet, file="./Figures/Eagle_Simulation_Type_Area_Plots_States.png",width=12, height=7,dpi=1200)

print(eagleAreaPlot_facet)


Discussion

Our simulation results demonstrate that in several states, hacking programs effectively augmented breeding populations. However, in an even greater number of instances, the observed rate of population growth during and after state hacking efforts differed from what was predicted under simulated conditions. These differences may be due to invalid assumptions or general limitations of our simulation model. For example, our use of state boundaries in defining populations was ecologically arbitrary and we did not account for density dependence even though breeding recruitment rate has decreased since 2005 (USFWS 2020).

Nonetheless, we maintain that a significant proportion of the variation between observed and simulated populations in states where hacking was applied is due to biological factors. This conclusion is supported by the differing degree of genetic admixture among U.S. regions with historic hacking programs (Judkins et al. 2020). States with higher than anticipated breeding populations may be due to additional management efforts such as artificial incubation and fostering programs (NJDEP 2020). Conversely, states with lower-than-expected breeding populations may be the result of varied hacking techniques (Sorenson et al. 2017), source and quantity of hacked birds (Judkins et al. 2020), or dispersal to neighboring states (Hatcher 1991). Given hacking’s apparent contribution to bald eagle recovery, and the potential need to apply this technique to other species and populations, we recommend further research into the differences between outcomes of historic hacking programs.

Where state observed breeding population numbers differ from simulated values There are several potential explanations for the difference between observed and simulated populations with hacking. including disparity between modeled and state-level demographic rates, emigration of hacked birds to adjacent states, hacking technique, source population and fitness of hacked birds, or undercounting of actual historic breeding population. Additionally, because our simulations did not incorporate the influence of density dependence, they may be unreflective of population dynamics after 2005, when national modeling indicates an increasing proportion of adult eagles were no longer recruited into the breeding population (USFWS 2020). Our simulation exercise provides general evidence that hacking often augmented bald eagle breeding populations and advanced the species’ rate of recovery at the state level. These findings argue in favor of hacking as an effective technique for recovering imperiled raptor species. Further, the variance in performance of hacking programs suggests underlying factors which, if once investigated and understood, could lead to improved implementation in the future. • Model does not incorporate density dependence. Rather, assumes full recruitment, which national models indicate has dropped from 99% in 2005 to 66% in 2018. • Simulations relied on arbitrary state boundaries for population definition; ample evidence bald eagles generally and hacked birds specifically dispersing to non-natal states • Results suggest varied success in hacking, supporting use as recovery technique for raptors, but also arguing for further investigation/optimization of techniques • Likely positive ‘spillover’ effect to adjacent states from hacking programs • Historic state nesting population numbers likely undercount of actual populations • Model assumes demographic rates do not vary by state –doubtful, particularly for HY survival, which is key demographic rate in population growth population growth


References

Gilbertson, M., 1992. Proceedings of the Third Expert Consultation Meeting on Bald Eagles in the Great Lakes Basin, February 25 to 26, 1992, International Joint Commission, Windsor, Ontario.

Hatcher, R.M., 1991. Computer model projections of Bald Eagle nesting in Tennessee. Journal of the Tennessee Academy of Science, 66(4), pp.225-228.

Jenkins, M.A. and Sherrod, S.K., 2005. Growth and recovery of the bald eagle population in Oklahoma. Wildlife Society Bulletin, 33(3), pp.810-813.

Judkins, M.E., Couger, B.M., Warren, W.C. and Van Den Bussche, R.A., 2020. A 50K SNP array reveals genetic structure for bald eagles (Haliaeetus leucocephalus). Conservation Genetics, 21(1), pp.65-76.

Murphy, T.M., 1989. Southeastern states bald eagle recovery plan. US Fish and Wildlife Service, Southeast Region.

New Jersey Department of Environmental Protection. 2020. New Jersey Bald Eagle Project, 2020. New Jersey Department of Environmental Protection, Division of Fish and Wildlife, Trenton, NJ, USA.

New York State Department of Environmental Conservation. 2016. Conservation Plan for Bald Eagles in New York State. New York State Department of Environmental Conservation, Division of Fish, Wildlife & Marine Resources, Bureau of Wildlife, Bald Eagle Team, Albany, NY, U.S.A.

R Core Team. 2020. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Sherrod, S.K., 1982. Hacking: a method for releasing peregrine falcons and other birds of prey. Peregrine Fund.

Sorenson, K.J., Burnett, L.J. and Stake, M.M., 2017. Restoring a bald eagle breeding population in central California and monitoring 25 years of regional population growth. Journal of Raptor Research, 51(2), pp.145-152.

U.S. Fish and Wildlife Service. 2016. Bald and Golden Eagles: Population demographics and estimation of sustainable take in the United States, 2016 update. Division of Migratory Bird Management, Washington D.C., USA.

U.S. Fish and Wildlife Service. 2020. Final Report: Bald Eagle Population Size: 2020 Update. U.S. Fish and Wildlife Service, Division of Migratory Bird Management, Washington, D.C., U.S.A.