In mjevans26/eaglesFWS:
library(dplyr)
library(plotly)
library(sensitivity)
library(shiny)
library(shinydashboard)
library(stringr)
library(tidyr)
library(viridis)
load("data/newdata.RData")
source("R/helper_fxns.R")
tab <- select(Bay16,SITE, TURBINES, EFFORT, FLIGHT_MIN)
prior <- curve(dgamma(x,
shape = expose$shape,
rate = expose$rate),
from = 0.01,
to = 6)
collision <- curve(dbeta(x,
shape1 = collide$shape,
#shape1 = 9.38,
shape2 = collide$rate),
#shape2 = 3224.51),
from = 0,
to = 0.05)
act <- mean(Bay_16$FLIGHT_MIN)/mean(Bay_16$EFFORT)
demo <- curve(dnorm(x, 5, 1), from = 0, to = 10)
plotList <- function(x){
a <- expose$shape + Bay_16$FLIGHT_MIN[x]
b <- expose$rate + Bay_16$EFFORT[x]
obs <- Bay_16$FLIGHT_MIN[x]/Bay_16$EFFORT[x]
post <- curve(dgamma(x, shape = a, rate = b), 0, 6)
site <- Bay_16$SITE[x]
plot_ly()%>%
add_trace(name = "Combined",
x = post$x, y = post$y,
type = "scatter", mode = "lines", fill = "tozeroy",
fillcolor = 'blue', line = list(color = 'blue'),
text = ~paste("Combined Activity<br>estimate = ",
round(post$x, 2),
sep = ""),
hoverinfo = "text")%>%
add_trace(name = "Prior",
x = prior$x, y = prior$y,
type = "scatter", mode = "lines", fill = "tozeroy",
fillcolor = 'yellow', line = list(color = 'yellow'),
text = ~paste("Prior Activity<br>estimate = ",
round(prior$x, 2),
sep = ""),
hoverinfo = "text")%>%
add_trace(x = c(obs, obs), y = c(0,
max(c(prior$y, post$y))),
type = "scatter", mode = "lines",
name = "Observed", line = list(color = 'black'),
text = ~paste("Observed Activity<br>at site = ",
round(obs,2),
sep = ""),
hoverinfo = "text")%>%
layout(#title = site,
showlegend = FALSE,
yaxis = list(showticklabels = FALSE),
xaxis = list(range = c(0,6))
)
}
s1 <- subplot(lapply(1:16, plotList),
nrows = 4, shareY = TRUE, titleY = FALSE, shareX =TRUE,titleX= FALSE)
src_discrep <- src(X = sim3[, c(2, 6)], y = sim3$UQ_F - sim3$UQ, nboot = 100)
Background
The U.S. Fish and Wildlife Service (FWS) uses a Bayesian framework to predict eagle collisions and fatalities at proposed wind energy sites. One of the appeals of Bayesian modeling is the ability to incorporate prior information and uncertainty regarding parameter values into estimates, ideally providing a better representation of reality. However, there is still debate over the appropriate choice and use of priors. Additionally, as Bayesian methodology is still gaining widespread familiarity, the ways in which prior information affects model outcomes can be unclear.
The purpose of this paper is to illustrate and quantify the effects of Bayesian priors used by FWS on the predicted outcomes of the eagle mortality model, in terms of predicted eagle fatalities. We use both empirical data from wind farms provided in a study presenting the methods used by FWS to implement the Bayesian prediction model (Bay et al. 2016), as well as a simulated dataset to demonstrate these effects both in theory and in practice.
Fully understanding the effects of priors is relevant to both energy developers and wildlife conservationists, because estimates of eagle take are used to determine permitting and mitigation costs.
Bayesian Model
The components of the Bayesian model used by FWS to predict eagle fatality are:
- Eagle Exposure Rate: minutes per hour that eagles are estimated to fly within risk areas around turbines where collisions may occur.
- Collision Rate: probability an eagle flying within a turbine risk area is struck and killed.
- ExpansionFactor: constant converting fatality rate to annual fatalities, based on the size of a wind facility.
These parameters are represented as probability distributions, and combined to produce a posterior distribution representing the relative likelihood of how many ealges will be killed annually at a site.
$Fatalities = Exposure * Collision Rate * Expansion$
fluidPage(
fluidRow(
# column(4,
# plot_ly()%>%
# add_trace(x = ~collision$x, y = ~collision$y, type = "scatter", mode = "lines",
# fill = "tozeroy", name = "Prior", line = list(color = "grey"),
# hoverinfo = "none")%>%
# layout(title = "Collision",
# xaxis = list(title = "Collision Rate"),
# yaxis = list(title = "Probability",
# showticklabels=FALSE)
# )
# ),
# column(4,
# plot_ly()%>%
# add_trace(x = prior$x, y = prior$y,
# type = "scatter", mode = "lines", fill = "tozeroy",
# name = "Prior", line = list(color = vir_col(3)[1]),
# hoverinfo = "none")%>%
# layout(title = "Exposure",
# xaxis = list(title = "Exposure Rate"),
# yaxis = list(showticklabels = FALSE,
# title = NA)
# )
# ),
column(12,
plot_ly()%>%
add_trace(x = ~demo$x, y = ~demo$y, type = "scatter",
mode = "lines", fill = "tozeroy", line = list(color = "grey"),
showlegend = FALSE,
hoverinfo = "none")%>%
add_trace(x = ~c(mean(rnorm(10000, 5, 1)), mean(rnorm(10000, 5, 1))),
y = ~ c(0, max(demo$y)),
type = "scatter", mode = "lines", name = "Mean",
hoverinfo = "none")%>%
add_trace(x = ~c(quantile(rnorm(10000, 5, 1), 0.8),
quantile(rnorm(10000, 5, 1), 0.8)),
y = ~ c(0, max(demo$y)), type = "scatter", mode = "lines",
name = "80th Percentile", hoverinfo = "none")%>%
layout(
title = "Predicted Eagle Fatalities",
xaxis = list(title = "Eagle Fatalities"),
yaxis = list(showticklabels = FALSE,
title = "Probability"),
legend = list(x = 0.7, y = 0.9)
)
)
),
fluidRow(
column(12,
p(class = "caption",
tags$b("Figure 1."), " Demonstration output of Bayesian eagle fatality model.")
)
)
)
To be conservative, FWS uses the 80th percentile of this posterior distribution in setting take limits, and we display these values unless otherwise noted.
In this analysis, we focus on eagle exposure rate, as this parameter is subject to both site-specific survey data, and a Bayesian prior distribution. The expansion factor is determined completely by the characteristics of a site. Similarly, the collision rate is estimated entirely by a prior distribution based on rates observed at existing wind facilities.
Survey Data
The data shown in Table 1 are the survey data used as example values for the collision risk model. These were derived from an empirical dataset (Appendix A in Bay et al. (2016)), which contains pre-construction eagle survey data from 26 wind energy sites. Survey effort indicates the amount of time and area covered during pre-construction surveys, and eagle observations are the duration of time over which eagles were observed flying within survey areas.
fluidPage(
fluidRow(
column(12,
p(class = "caption",
tags$b("Table 1."), " Data from Bay et al. (2016). Effort was calculated as the product of survey area (km", tags$sup("2"),") and observation time (hr),
derived from Appendix A.")
)
),
fluidRow(
div(
column(12,
DT::datatable(tab,
fillContainer = TRUE,
selection = list(mode = 'single',
#selected = which(Bay16$SITE == input$sites),
target = 'row'),
colnames = c("Site", "Turbines", "Survey Effort (hr*km2)",
"Eagle Obs (min)"),
options = list(rownames = FALSE,
pageLength = nrow(tab),
dom = 'tip'
)
)%>%
DT::formatRound("EFFORT", 2)
))
)
)
Combining Prior Information with Survey Data
Eagle flight time and survey effort are used to estimate eagle exposure at a given wind site. Figure 2 shows how the exposure measured at a site (green lines) is integrated with a consistent prior distribution of exposure probabilities (purple), to produce a probability distribution of exposure (yellow). This resulting distribution is used as the exposure rate to predict eagle fatalities at a site. The effect of the prior distribution is to modify extreme observed values that differ greatly from mean observations across a larger sample of sites.
fluidPage(
fluidRow(
column(12,
s1%>%
layout(xaxis = list(title = "Eagle Exposure (min/km<sup>3</sup>*hr)",
align = "center"),
yaxis = list(title = "Probability",
showticklabels = FALSE,
align = "center")
)
)
),
fluidRow(
column(12,
p(class = "caption",
tags$b("Figure 2."), " Eagle exposure rates used to estimate predicted fatalities
across from a sample of wind farms. A consistent prior probability distribution (purple) is combined with site-specific
survey data (green line) to estimate a distribution of exposure probabilities for a
proposed project (yellow).")
)
)
)
Using Priors Changes Estimates
At the majority of sites, the Bayesian estimate of eagle fatalities using prior exposure information was higher than a site-specific estimate (Fig. 2). These instances are indicated by points above the 1:1 line in Figure 2. The influence of eagle exposure priors on estimated eagle fatalities was less at facilities that expended greater effort on pre-construction surveys. This pattern is indicated by the sites with the most intensive surveys falling closest to a 1:1 relationship between the two estimates (Fig. 2).
fluidPage(
fluidRow(
column(12,
plot_ly(Bay16)%>%
add_trace(y = ~UQ_F, x = ~UQ, type = "scatter", mode = "markers",
size = ~((UQ_F)-(MN_F)),
marker = list(line = list(color = "black"),
sizemode = "diameter"
),
sizes = c(10,25), color = ~ EFFORT, name = "Estimated Eagle Fatalities",
text = ~paste("Site =", SITE, "<br>Effort =", round(EFFORT,2),
"hr*km<sup>3</sup>", "<br>Prediction w/Priors =",
round(UQ_F, 2),
"<br>Prediction w/o Priors =", round(UQ, 2)),
hoverinfo = "text")%>%
add_trace(x = ~c(0, max(UQ_F)), y = ~c(0, max(UQ_F)), type = "scatter",
mode = "lines", name = "1:1 Line")%>%
colorbar (y = 0.7, x = 0.80,
title = "Survey Effort<br>hr*km<sup>3</sup>")%>%
layout(hovermode = "closest", font = list(color = "black"),
yaxis = list(title = "Estimates using Priors"),
xaxis = list(title = "Estimates without Priors"),
legend = list(x = 0.4, y = 0.9, bordercolor = "black",
borderwidth = 1)
)
)
),
fluidRow(
column(12,
p(class = "caption",
tags$b("Figure 3."), "Relationships between estimated eagle fatalities produced incorporating prior information about eagle exposure at wind facilities, and those produced using only site-specific survey data. Plotted values are the 80th percentiles of the posterior distributions. Marker size is proportional to the size of the 80% CI produced from the posterior distribution using exposure priors.")
)
)
)
Source: [Bay et al. (2016)](http://www.pnas.org/content/113/13/3563.full?tab=ds)
Simulation Study
To investigate the effects of prior information and survey effort on predicted take, we created a dataset representing a range of possible pre-construction survey values for survey area, time, and eagle observations. We then produced eagle fatality estimates for these possible values using the FWS bayesian model.
Prediction Discrepancies Increase at Extreme Values
For eagle exposure rates < 1.21 min/hr*km^3^, predicted fatlity rates estimated using priors were higher than site-specific estimates (i.e. above the 1:1 line in Fig. 3). This pattern reversed when exposure was > 1.21. The threshold of 1.21 corresponds to the mean exposure rate of the prior exposure distribution. The actual value of this threshold is subject to change as prior distributions are updated with new survey data, but the pattern will be consistent. Figure 4 also illustrates that predicted fatalities estimated with and without priors are more similar when survey effort is greater (i.e., larger circles).
fluidPage(
fluidRow(
column(12,
plot_ly(sim3[!(sim3$eagle_rate %in% c(0.02, 0.03, 0.04, 0.06, 0.07, 0.08, 0.09)),])%>%
add_trace(x = ~UQ, y = ~UQ_F, type = "scatter",
mode = "markers", size = ~b,
marker = list(line = list(color = "black"),
sizemode = "diameter"
),
sizes = c(5,15), color = ~ eagle_rate, name = "Estimated Eagle Fatalities",
text = ~paste("Effort =", round(b, 3), "hr*km<sup>3</sup>",
"<br>Exposure =", a, "min"),
hoverinfo = "text")%>%
add_trace(x = ~c(0, max(UQ_F)), y = ~c(0, max(UQ_F)),
type = "scatter", mode = "lines", name = "1:1 Line")%>%
# add_trace(x = ~c(0.002, 0.006), y = ~c(0.004, 0.004),
# type = "scatter", mode = "lines", name = "Mean")%>%
colorbar(y = 1, x = 0.1,
title = "Eagle Exposure<br>Rate (min/hr)")%>%
layout(hovermode = "closest", font = list(color = "black"),
xaxis = list(title = "Estimates without Priors"),
yaxis = list(title = "Estimates using Priors"),
legend = list(x = 0.6, y = 0.3, bordercolor = "black", borderwidth = 1)
)
)
),
fluidRow(
column(12,
p(class = "caption",
tags$b("Figure 4. The mean of the exposure prior determines whether Bayesian estimates under or over predict fatalities relative to survey data"), "Plotted values are the 80th percentile of posterior fatality distributions estimated with and without exposure priors, multiplied by the mean site scaling factor from Bay et al. 2016. Circle size is proportional to survey effort")
)
)
)
The prior distribution of exposure probabilities has a greater effect on estimates when observed exposure rate at a site is farther from the mean of the prior. Figure 5 illustrates this relationship in terms of how many standard deviations an observed exposure rate was from the prior mean. At one standard deviation above or below the mean, take estimates can differ by as much as 1.5 eagles per year.
fluidPage(
fluidRow(
plot_ly(sim3[!(sim3$eagle_rate %in% c(0.02, 0.03, 0.04, 0.06, 0.07, 0.08, 0.09)),])%>%
add_trace(x = ~(eagle_rate - 1.1)/0.107,
y = ~(UQ_F-UQ), type = "scatter", mode = "markers",
color = ~ b,
text = ~paste("Effort =", round(b, 3), "<br>Z =",
round((eagle_rate - 1.1)/0.107, 2),
"<br>Discrepancy =", round((UQ_F-UQ), 2)),
hoverinfo = 'text')%>%
colorbar(x = 0.9, y = 1, title = "Survey Effort<br>(hr*km<sup>3</sup>",
lenmode = 'fraction', len = 0.4)%>%
layout(hovermode = 'closest',
font = list(color = 'black'),
xaxis = list(title = " Difference between Observed Exposure and Prior Mean (standard deviations)"
),
yaxis = list(title = "Fatality Estimate Discrepancy")
)
),
fluidRow(
column(12,
p(class = "caption",
tags$b("Figure 5. Discrepancy between Bayesian and site-specific estimates of eagle fatalities are greater at extreme values."), " Differences between the observed exposure rate and prior mean were standardized to z-scores. Discrepancies were calculated between the 80th percentile of the posterior fatality distributions. Values were generated from simulated data, and multiplied by the mean site scaling factor from Bay et al. 2016.")
)
)
)
Greater Survey Effort Reduces Effect of Priors
Increased survey effort reduced the influence of prior distributions of exposure probability on estimates of fatalities, resulting in less discrepacy between predictions obtained using priors and those obtained using site survey data only. The magnitude of the effect of effort on estimate discrepancy was contingent upon the observed eagle exposure rate. This discrepancy decreased logarithmically with greater survey effort
(Fig.6). The precision of posterior estimates, measured by the difference between the 80th percentile and mean estimate standardized by the mean, did not change as survey efforts increased.
#sim$b2 <- sim$b^2
dev <- glm(data = sim3, (UQ_F - UQ) ~ b*eagle_rate, family = gaussian(link = "inverse"))
dev_max <- glm(data = sim3, (UQ_F - UQ)~ b, subset = eagle_rate == max(eagle_rate),
family = gaussian(link = "inverse"))
dev_min <- glm(data = sim3, (UQ_F - UQ)~ sqrt(b), subset = eagle_rate == min(eagle_rate),
family = gaussian(link = "inverse"))
int <- glm(data = sim3, (UQ_F-MN_F)~ b * eagle_rate)
fluidPage(
fluidRow(
column(6,
plot_ly(sim3)%>%
add_trace(x = ~ b, y = ~(UQ_F - UQ),
type = "scatter", mode = "markers",
color = ~ eagle_rate,
text = ~paste("Effort =", round(b, 3), "hr*km<sup>2</sup>", "<br>Exposure =", eagle_rate, "min"),
hoverinfo = 'text',
showlegend = FALSE)%>%
#add_trace(x = ~c(10, 10),
# y = ~c(min(UQ_F - UQ), max(UQ_F - UQ)),
# type = "scatter", mode = "lines", name = "FWS Minimum<br>Survey Effort")%>%
colorbar(title = "Eagle Exposure<br>Rate (min/hr)", x = 0.8, y = 0.8,
lenmode = 'fraction', len = 0.35)%>%
layout(hovermode = 'closest',
font = list(color = 'black'),
xaxis = list(title = "Survey Effort (hr*km<sup>3</sup>)"
),
yaxis = list(title = "Fatality Estimate Discrepancy"),
legend = list(x = 0.5, y = 0.2)
)
),
column(6,
plot_ly(sim3[!(sim$eagle_rate %in% c(0.02, 0.03, 0.04, 0.06, 0.07, 0.08, 0.09)),])%>%
add_trace(x = ~b, y = ~(UQ_F-MN_F)/MN_F,
type = "scatter", mode = "markers",
color = ~ eagle_rate,
marker = list(colorbar = list(colorbar = list(showscale = FALSE))
), showlegend = FALSE,
text = ~paste("Effort =", round(b, 3), "hr*km<sup>3</sup>",
"Exposure =", eagle_rate , "min"
),
hoverinfo = 'text')%>%
add_trace(x = ~c(10, 10), y = ~c(min((UQ_F-MN_F)/MN_F), max((UQ_F-MN_F)/MN_F)),
type = "scatter",
mode = "lines", name = "FWS Minimum<br>Survey Effort",
showlegend = FALSE)%>%
hide_colorbar()%>%
layout(hovermode = 'closest',
font = list(color = 'black'),
xaxis = list(title = "Survey Effort (hr*km<sup>3</sup>)"
),
yaxis = list(title = "Precision ((80th% - Mean)/Mean)")
)
)
),
fluidRow(
column(12,
p(class = "caption",
tags$b("Figure 6. Estimates of eagle fatality become more similar as estimated with and without priors with greater survey effort."), " Changes in a) difference between fatalities predicted with and without prior exposure information, and b) the size of the 80% CI of the predicted fatality distribution, relative to the mean, as a funciton of survey effort.")
)
)
)
Harder to Prove Big Projects Take No Eagles
A bayesian approach to estimating fatalities guards against situations where surveys fail to detect eagles that are in fact present around the site. However, this means that sites where there are truly no birds will have some predicted take. In situations where no eagles are observed during pre-construciton surveys, greater survey effort moves predicted take closer to zero (Fig. 7). This response is highly dependent on the size of the wind farm. For the smallest sites, 48 hr.km^2^ of survey effort reduced fatality predictions < 1. Large sites with no eagle observations require more survey effort to reduce predicted take to zero (~300 hr.km^2^), however the rate of reduction is greater at large sites than at smaller sites.
#PLOT FOR ZERO OBSERVED EAGLES
fluidPage(
fluidRow(
plot_ly(data = zerosim, x = ~Effort, y = ~UQ_F,
color = ~Hectares, type = 'scatter', mode = 'markers',
marker = list(colorbar = list(x = 0.7, y = 1, title = "Project<br>Size (ha)")),
hoverinfo = 'text',
text = ~paste("Effort:", Effort, "hr*km<sup>3</sup>",
"<br>Project Size:", Hectares, "ha",
"<br>Fatalities:", round(UQ_F, 1), "eagles"))%>%
colorbar(x = 0.8, y = 1, title = 'Project Size<br>(# Turbines)')%>%
layout(legend = list(x = 0.7, y = 1),
xaxis = list(title = "Effort (hr*km2)"),
yaxis = list(title = "Predicted Annual Eagle Fatalities")
)
),
fluidRow(
column(12,
p(class = 'caption',
tags$b("Figure 7. Project size affects predicted take at sites where no eagles are observed."), " Changes in predicted fatalities, given no observed eagles, as a function of survey effort for different sized wind farms.")
)
)
)
Methods
We generated posterior distributions for predicted eagle fatalities from empirical and simulation data using the FWS Bayesian model. The prior distribution on eagle exposure used by FWS is a gamma distribution Exposure ~ Gamma (0.968, 0.552), with shape and rate parameters represented by eagle flight time in minutes, and survey effort in (hr km^3^). The prior distribution on collisions used by FWS is a beta distribution Collision ~ Beta (2.31, 396.69) with parameters representative of number of collisions and minutes of exposure.
Empirical Data
We used wind energy site survey data from Appendix A in Bay et al. (2016) to obtain example survey values. This data recorded plot area in hectares, and observation time in minutes. We converted these to survey effort consistent with the units used by FWS (hr km^3^) by rescaling plot area (ha) by 0.2 km/ha, and multiplying by the reported observation time in hours.
Simulation Data
We generated hypothetical values for survey effort based on the minimum requirements provided by U.S. FWS for pre-construction monitoring. FWS requires at least one cylindrical survey plot with radius >= 800m and height >= 200m, and that plots be surveyed for at least 12 hours per year, for two years. Thus, the minimum values for survey area and time were r min(area) km3, and r min(time) hrs, respectively. We simulated up to five plots (area = r max(area)) in increments of one, and up to r max(time) hrs in increments of 12, and calculated survey effort for all combinations.
Observed eagle flight time (min) provided in Bay et al. (2016) is a function of survey effort and eagle activity at a site. Therefore, to generate a range of potential exposure values (min/hr*km^3^), we divided flight time by effort at each survey site to obtain an effort invariant measure of eagle exposure rate, and took a random sample of 20 values. The final simulation data included all combinations of survey effort and eagle exposure rate, which we multiplied to obtain flight minutes, producing 1000 values.
flight <- sample(BAY16$MIN/Bay16$EFFORT, 20)
survey_time <- seq(1, 10, 1)*12*2
survey_area <- seq(0.402, 2.01, 0.402)
df <- expand.grid(TIME = time, AREA = area, eagle_rate = flight)
df$beta <- df$TIME*df$AREA
df$alpha <- df$eagle_rate*df$beta
Model Output
We generated eagle fatality predictions using only simulated or empirical values for observed eagle flight time (min) and survey effort (hr km3), and by updating the prior probability of exposure using these values. We obtained the posterior exposure distribution as in New et al. (2015) Exposure ~ Gamma(a + minutes, b + effort), where a and b are the shape and rate parameters of the prior exposure distribution.
# simFatal funciton from New et al. (2013) modified to perform multiple iterations.
# BMin: observed number of bird minutes
# Fatal: annual avian fatalities on an operational wind facility
# SmpHrKm: total time and area surveyed for bird minutes
# ExpFac: expansion factor
# aPriExp: alpha parameter for the prior on lambda
# bPriExp: beta parameter for the prior on lambda
# aPriCPr: alpha parameter for the prior on C
# bPriCPr: beta parameter for the prior on C
# iters: number of sampling iterations
simFatal <- function(BMin=-1, Fatal=-1, SmpHrKm, ExpFac, aPriExp=1,
bPriExp=1,aPriCPr=1, bPriCPr=1, iters){
out <- data.frame(collision = rep(NA,iters),
expose = rep(NA, iters),
fatality = rep(NA, iters)
)
# The default of a negative value for BMin or Fatal indicates no data
if(BMin>=0){
aPostExp <- aPriExp + BMin
bPostExp <- bPriExp + SmpHrKm
}else{
aPostExp <- aPriExp
bPostExp <- bPriExp}
# Update the collisions prior
if(Fatal>=0){
aPostCPr <- aPriCPr + Fatal
bPostCPr <- ((rvmean(Exp) * ExpFac) - Fatal) + bPriCPr
}else{
aPostCPr <- aPriCPr
bPostCPr <- bPriCPr}
for(i in 1:iters){
Exp <- rgamma(n=1, aPostExp, bPostExp)
CPr <- rbeta(n=1, aPostCPr, bPostCPr)
Fatalities <- ExpFac * Exp * CPr
out[i,] <- c(CPr, Exp, Fatalities)
}
return(out)
}
estimates <- function(iters, a, b){
out <- simFatal(BMin = a,
Fatal = -1,
SmpHrKm = b,
ExpFac = mean(Bay16$SCALE),
aPriExp = 11.81641,
bPriExp = 9.7656250,
aPriCPr = 1.638029,
bPriCPr = 290.0193,
iters = iters)
fatality <- mean(out$fatality)
q80 <- quantile(out$fatality, c(0.8))
out2 <- simFatal(BMin = a,
Fatal = -1,
SmpHrKm = b,
ExpFac = mean(Bay16$SCALE),
aPriExp = 0,
bPriExp = 0,
aPriCPr = 1.638029,
bPriCPr = 290.0193,
iters = iters)
fatality2 <- mean(out2$fatality)
q82 <- quantile(out2$fatality, 0.8)
return (c("MN_F" = fatality, "U_F" = q80, "MN" = fatality2, "U" = q82))
}
We used an iterative sampling procedure to produce a posterior distribution of annual eagle fatalities. During each iteration, we drew a random sample from the updated exposure distribution and the prior collision probability distribution, and multiplied these values together with a site specific expansion factor. For simulated values, we used the mean expansion factor from the empirical data. For each site and set of simulated values, we performed 100000 iterations of this sampling procedure, and calculated the mean and 80th percentile values of the resulting distribution.
#Predictions integrating Exposure prior and survey data
sim <- plyr::mdply(df[, c(5, 4)], estimates, niters = 100000)
#Predictions from survey data
emp <- vapply(1:nrow(Bay16), function(x){
estimates(100000, Bay16$FLIGHT_MIN[x], BAY16$EFFORT[x])
return (c(fatality, q80, fatality2, q82))
},
USE.NAMES = FALSE, FUN.VALUE = c(0,0,0,0))
)
fluidPage(
fluidRow(
column(2),
column(8,
div(
HTML('<a rel="license" href="http://creativecommons.org/licenses/by-nc-sa/4.0/"> <img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-nc-sa/4.0/88x31.png" /></a> <br /> This <span xmlns:dct="http://purl.org/dc/terms/" href="http://purl.org/dc/dcmitype/InteractiveResource" rel="dct:type">work</span> by <a xmlns:cc="http://creativecommons.org/ns" href="http://defenders.org" property="cc:attributionName" rel="cc:attributionURL">Defenders of Wildlife</a> is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-nc-sa/4.0/">Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License</a>. <br />'),
style = "text-align: center"
)
),
column(2)
)
)
mjevans26/eaglesFWS documentation built on Dec. 29, 2021, 1:35 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(dplyr) library(plotly) library(sensitivity) library(shiny) library(shinydashboard) library(stringr) library(tidyr) library(viridis) load("data/newdata.RData") source("R/helper_fxns.R") tab <- select(Bay16,SITE, TURBINES, EFFORT, FLIGHT_MIN) prior <- curve(dgamma(x, shape = expose$shape, rate = expose$rate), from = 0.01, to = 6) collision <- curve(dbeta(x, shape1 = collide$shape, #shape1 = 9.38, shape2 = collide$rate), #shape2 = 3224.51), from = 0, to = 0.05) act <- mean(Bay_16$FLIGHT_MIN)/mean(Bay_16$EFFORT) demo <- curve(dnorm(x, 5, 1), from = 0, to = 10) plotList <- function(x){ a <- expose$shape + Bay_16$FLIGHT_MIN[x] b <- expose$rate + Bay_16$EFFORT[x] obs <- Bay_16$FLIGHT_MIN[x]/Bay_16$EFFORT[x] post <- curve(dgamma(x, shape = a, rate = b), 0, 6) site <- Bay_16$SITE[x] plot_ly()%>% add_trace(name = "Combined", x = post$x, y = post$y, type = "scatter", mode = "lines", fill = "tozeroy", fillcolor = 'blue', line = list(color = 'blue'), text = ~paste("Combined Activity<br>estimate = ", round(post$x, 2), sep = ""), hoverinfo = "text")%>% add_trace(name = "Prior", x = prior$x, y = prior$y, type = "scatter", mode = "lines", fill = "tozeroy", fillcolor = 'yellow', line = list(color = 'yellow'), text = ~paste("Prior Activity<br>estimate = ", round(prior$x, 2), sep = ""), hoverinfo = "text")%>% add_trace(x = c(obs, obs), y = c(0, max(c(prior$y, post$y))), type = "scatter", mode = "lines", name = "Observed", line = list(color = 'black'), text = ~paste("Observed Activity<br>at site = ", round(obs,2), sep = ""), hoverinfo = "text")%>% layout(#title = site, showlegend = FALSE, yaxis = list(showticklabels = FALSE), xaxis = list(range = c(0,6)) ) } s1 <- subplot(lapply(1:16, plotList), nrows = 4, shareY = TRUE, titleY = FALSE, shareX =TRUE,titleX= FALSE) src_discrep <- src(X = sim3[, c(2, 6)], y = sim3$UQ_F - sim3$UQ, nboot = 100)
Background
The U.S. Fish and Wildlife Service (FWS) uses a Bayesian framework to predict eagle collisions and fatalities at proposed wind energy sites. One of the appeals of Bayesian modeling is the ability to incorporate prior information and uncertainty regarding parameter values into estimates, ideally providing a better representation of reality. However, there is still debate over the appropriate choice and use of priors. Additionally, as Bayesian methodology is still gaining widespread familiarity, the ways in which prior information affects model outcomes can be unclear.
The purpose of this paper is to illustrate and quantify the effects of Bayesian priors used by FWS on the predicted outcomes of the eagle mortality model, in terms of predicted eagle fatalities. We use both empirical data from wind farms provided in a study presenting the methods used by FWS to implement the Bayesian prediction model (Bay et al. 2016), as well as a simulated dataset to demonstrate these effects both in theory and in practice.
Fully understanding the effects of priors is relevant to both energy developers and wildlife conservationists, because estimates of eagle take are used to determine permitting and mitigation costs.
Bayesian Model
The components of the Bayesian model used by FWS to predict eagle fatality are:
- Eagle Exposure Rate: minutes per hour that eagles are estimated to fly within risk areas around turbines where collisions may occur.
- Collision Rate: probability an eagle flying within a turbine risk area is struck and killed.
- ExpansionFactor: constant converting fatality rate to annual fatalities, based on the size of a wind facility.
These parameters are represented as probability distributions, and combined to produce a posterior distribution representing the relative likelihood of how many ealges will be killed annually at a site.
$Fatalities = Exposure * Collision Rate * Expansion$
fluidPage( fluidRow( # column(4, # plot_ly()%>% # add_trace(x = ~collision$x, y = ~collision$y, type = "scatter", mode = "lines", # fill = "tozeroy", name = "Prior", line = list(color = "grey"), # hoverinfo = "none")%>% # layout(title = "Collision", # xaxis = list(title = "Collision Rate"), # yaxis = list(title = "Probability", # showticklabels=FALSE) # ) # ), # column(4, # plot_ly()%>% # add_trace(x = prior$x, y = prior$y, # type = "scatter", mode = "lines", fill = "tozeroy", # name = "Prior", line = list(color = vir_col(3)[1]), # hoverinfo = "none")%>% # layout(title = "Exposure", # xaxis = list(title = "Exposure Rate"), # yaxis = list(showticklabels = FALSE, # title = NA) # ) # ), column(12, plot_ly()%>% add_trace(x = ~demo$x, y = ~demo$y, type = "scatter", mode = "lines", fill = "tozeroy", line = list(color = "grey"), showlegend = FALSE, hoverinfo = "none")%>% add_trace(x = ~c(mean(rnorm(10000, 5, 1)), mean(rnorm(10000, 5, 1))), y = ~ c(0, max(demo$y)), type = "scatter", mode = "lines", name = "Mean", hoverinfo = "none")%>% add_trace(x = ~c(quantile(rnorm(10000, 5, 1), 0.8), quantile(rnorm(10000, 5, 1), 0.8)), y = ~ c(0, max(demo$y)), type = "scatter", mode = "lines", name = "80th Percentile", hoverinfo = "none")%>% layout( title = "Predicted Eagle Fatalities", xaxis = list(title = "Eagle Fatalities"), yaxis = list(showticklabels = FALSE, title = "Probability"), legend = list(x = 0.7, y = 0.9) ) ) ), fluidRow( column(12, p(class = "caption", tags$b("Figure 1."), " Demonstration output of Bayesian eagle fatality model.") ) ) )
To be conservative, FWS uses the 80th percentile of this posterior distribution in setting take limits, and we display these values unless otherwise noted.
In this analysis, we focus on eagle exposure rate, as this parameter is subject to both site-specific survey data, and a Bayesian prior distribution. The expansion factor is determined completely by the characteristics of a site. Similarly, the collision rate is estimated entirely by a prior distribution based on rates observed at existing wind facilities.
Survey Data
The data shown in Table 1 are the survey data used as example values for the collision risk model. These were derived from an empirical dataset (Appendix A in Bay et al. (2016)), which contains pre-construction eagle survey data from 26 wind energy sites. Survey effort indicates the amount of time and area covered during pre-construction surveys, and eagle observations are the duration of time over which eagles were observed flying within survey areas.
fluidPage( fluidRow( column(12, p(class = "caption", tags$b("Table 1."), " Data from Bay et al. (2016). Effort was calculated as the product of survey area (km", tags$sup("2"),") and observation time (hr), derived from Appendix A.") ) ), fluidRow( div( column(12, DT::datatable(tab, fillContainer = TRUE, selection = list(mode = 'single', #selected = which(Bay16$SITE == input$sites), target = 'row'), colnames = c("Site", "Turbines", "Survey Effort (hr*km2)", "Eagle Obs (min)"), options = list(rownames = FALSE, pageLength = nrow(tab), dom = 'tip' ) )%>% DT::formatRound("EFFORT", 2) )) ) )
Combining Prior Information with Survey Data
Eagle flight time and survey effort are used to estimate eagle exposure at a given wind site. Figure 2 shows how the exposure measured at a site (green lines) is integrated with a consistent prior distribution of exposure probabilities (purple), to produce a probability distribution of exposure (yellow). This resulting distribution is used as the exposure rate to predict eagle fatalities at a site. The effect of the prior distribution is to modify extreme observed values that differ greatly from mean observations across a larger sample of sites.
fluidPage( fluidRow( column(12, s1%>% layout(xaxis = list(title = "Eagle Exposure (min/km<sup>3</sup>*hr)", align = "center"), yaxis = list(title = "Probability", showticklabels = FALSE, align = "center") ) ) ), fluidRow( column(12, p(class = "caption", tags$b("Figure 2."), " Eagle exposure rates used to estimate predicted fatalities across from a sample of wind farms. A consistent prior probability distribution (purple) is combined with site-specific survey data (green line) to estimate a distribution of exposure probabilities for a proposed project (yellow).") ) ) )
Using Priors Changes Estimates
At the majority of sites, the Bayesian estimate of eagle fatalities using prior exposure information was higher than a site-specific estimate (Fig. 2). These instances are indicated by points above the 1:1 line in Figure 2. The influence of eagle exposure priors on estimated eagle fatalities was less at facilities that expended greater effort on pre-construction surveys. This pattern is indicated by the sites with the most intensive surveys falling closest to a 1:1 relationship between the two estimates (Fig. 2).
fluidPage( fluidRow( column(12, plot_ly(Bay16)%>% add_trace(y = ~UQ_F, x = ~UQ, type = "scatter", mode = "markers", size = ~((UQ_F)-(MN_F)), marker = list(line = list(color = "black"), sizemode = "diameter" ), sizes = c(10,25), color = ~ EFFORT, name = "Estimated Eagle Fatalities", text = ~paste("Site =", SITE, "<br>Effort =", round(EFFORT,2), "hr*km<sup>3</sup>", "<br>Prediction w/Priors =", round(UQ_F, 2), "<br>Prediction w/o Priors =", round(UQ, 2)), hoverinfo = "text")%>% add_trace(x = ~c(0, max(UQ_F)), y = ~c(0, max(UQ_F)), type = "scatter", mode = "lines", name = "1:1 Line")%>% colorbar (y = 0.7, x = 0.80, title = "Survey Effort<br>hr*km<sup>3</sup>")%>% layout(hovermode = "closest", font = list(color = "black"), yaxis = list(title = "Estimates using Priors"), xaxis = list(title = "Estimates without Priors"), legend = list(x = 0.4, y = 0.9, bordercolor = "black", borderwidth = 1) ) ) ), fluidRow( column(12, p(class = "caption", tags$b("Figure 3."), "Relationships between estimated eagle fatalities produced incorporating prior information about eagle exposure at wind facilities, and those produced using only site-specific survey data. Plotted values are the 80th percentiles of the posterior distributions. Marker size is proportional to the size of the 80% CI produced from the posterior distribution using exposure priors.") ) ) )
Source: [Bay et al. (2016)](http://www.pnas.org/content/113/13/3563.full?tab=ds)
Simulation Study
To investigate the effects of prior information and survey effort on predicted take, we created a dataset representing a range of possible pre-construction survey values for survey area, time, and eagle observations. We then produced eagle fatality estimates for these possible values using the FWS bayesian model.
Prediction Discrepancies Increase at Extreme Values
For eagle exposure rates < 1.21 min/hr*km^3^, predicted fatlity rates estimated using priors were higher than site-specific estimates (i.e. above the 1:1 line in Fig. 3). This pattern reversed when exposure was > 1.21. The threshold of 1.21 corresponds to the mean exposure rate of the prior exposure distribution. The actual value of this threshold is subject to change as prior distributions are updated with new survey data, but the pattern will be consistent. Figure 4 also illustrates that predicted fatalities estimated with and without priors are more similar when survey effort is greater (i.e., larger circles).
fluidPage( fluidRow( column(12, plot_ly(sim3[!(sim3$eagle_rate %in% c(0.02, 0.03, 0.04, 0.06, 0.07, 0.08, 0.09)),])%>% add_trace(x = ~UQ, y = ~UQ_F, type = "scatter", mode = "markers", size = ~b, marker = list(line = list(color = "black"), sizemode = "diameter" ), sizes = c(5,15), color = ~ eagle_rate, name = "Estimated Eagle Fatalities", text = ~paste("Effort =", round(b, 3), "hr*km<sup>3</sup>", "<br>Exposure =", a, "min"), hoverinfo = "text")%>% add_trace(x = ~c(0, max(UQ_F)), y = ~c(0, max(UQ_F)), type = "scatter", mode = "lines", name = "1:1 Line")%>% # add_trace(x = ~c(0.002, 0.006), y = ~c(0.004, 0.004), # type = "scatter", mode = "lines", name = "Mean")%>% colorbar(y = 1, x = 0.1, title = "Eagle Exposure<br>Rate (min/hr)")%>% layout(hovermode = "closest", font = list(color = "black"), xaxis = list(title = "Estimates without Priors"), yaxis = list(title = "Estimates using Priors"), legend = list(x = 0.6, y = 0.3, bordercolor = "black", borderwidth = 1) ) ) ), fluidRow( column(12, p(class = "caption", tags$b("Figure 4. The mean of the exposure prior determines whether Bayesian estimates under or over predict fatalities relative to survey data"), "Plotted values are the 80th percentile of posterior fatality distributions estimated with and without exposure priors, multiplied by the mean site scaling factor from Bay et al. 2016. Circle size is proportional to survey effort") ) ) )
The prior distribution of exposure probabilities has a greater effect on estimates when observed exposure rate at a site is farther from the mean of the prior. Figure 5 illustrates this relationship in terms of how many standard deviations an observed exposure rate was from the prior mean. At one standard deviation above or below the mean, take estimates can differ by as much as 1.5 eagles per year.
fluidPage( fluidRow( plot_ly(sim3[!(sim3$eagle_rate %in% c(0.02, 0.03, 0.04, 0.06, 0.07, 0.08, 0.09)),])%>% add_trace(x = ~(eagle_rate - 1.1)/0.107, y = ~(UQ_F-UQ), type = "scatter", mode = "markers", color = ~ b, text = ~paste("Effort =", round(b, 3), "<br>Z =", round((eagle_rate - 1.1)/0.107, 2), "<br>Discrepancy =", round((UQ_F-UQ), 2)), hoverinfo = 'text')%>% colorbar(x = 0.9, y = 1, title = "Survey Effort<br>(hr*km<sup>3</sup>", lenmode = 'fraction', len = 0.4)%>% layout(hovermode = 'closest', font = list(color = 'black'), xaxis = list(title = " Difference between Observed Exposure and Prior Mean (standard deviations)" ), yaxis = list(title = "Fatality Estimate Discrepancy") ) ), fluidRow( column(12, p(class = "caption", tags$b("Figure 5. Discrepancy between Bayesian and site-specific estimates of eagle fatalities are greater at extreme values."), " Differences between the observed exposure rate and prior mean were standardized to z-scores. Discrepancies were calculated between the 80th percentile of the posterior fatality distributions. Values were generated from simulated data, and multiplied by the mean site scaling factor from Bay et al. 2016.") ) ) )
Greater Survey Effort Reduces Effect of Priors
Increased survey effort reduced the influence of prior distributions of exposure probability on estimates of fatalities, resulting in less discrepacy between predictions obtained using priors and those obtained using site survey data only. The magnitude of the effect of effort on estimate discrepancy was contingent upon the observed eagle exposure rate. This discrepancy decreased logarithmically with greater survey effort (Fig.6). The precision of posterior estimates, measured by the difference between the 80th percentile and mean estimate standardized by the mean, did not change as survey efforts increased.
#sim$b2 <- sim$b^2 dev <- glm(data = sim3, (UQ_F - UQ) ~ b*eagle_rate, family = gaussian(link = "inverse")) dev_max <- glm(data = sim3, (UQ_F - UQ)~ b, subset = eagle_rate == max(eagle_rate), family = gaussian(link = "inverse")) dev_min <- glm(data = sim3, (UQ_F - UQ)~ sqrt(b), subset = eagle_rate == min(eagle_rate), family = gaussian(link = "inverse")) int <- glm(data = sim3, (UQ_F-MN_F)~ b * eagle_rate) fluidPage( fluidRow( column(6, plot_ly(sim3)%>% add_trace(x = ~ b, y = ~(UQ_F - UQ), type = "scatter", mode = "markers", color = ~ eagle_rate, text = ~paste("Effort =", round(b, 3), "hr*km<sup>2</sup>", "<br>Exposure =", eagle_rate, "min"), hoverinfo = 'text', showlegend = FALSE)%>% #add_trace(x = ~c(10, 10), # y = ~c(min(UQ_F - UQ), max(UQ_F - UQ)), # type = "scatter", mode = "lines", name = "FWS Minimum<br>Survey Effort")%>% colorbar(title = "Eagle Exposure<br>Rate (min/hr)", x = 0.8, y = 0.8, lenmode = 'fraction', len = 0.35)%>% layout(hovermode = 'closest', font = list(color = 'black'), xaxis = list(title = "Survey Effort (hr*km<sup>3</sup>)" ), yaxis = list(title = "Fatality Estimate Discrepancy"), legend = list(x = 0.5, y = 0.2) ) ), column(6, plot_ly(sim3[!(sim$eagle_rate %in% c(0.02, 0.03, 0.04, 0.06, 0.07, 0.08, 0.09)),])%>% add_trace(x = ~b, y = ~(UQ_F-MN_F)/MN_F, type = "scatter", mode = "markers", color = ~ eagle_rate, marker = list(colorbar = list(colorbar = list(showscale = FALSE)) ), showlegend = FALSE, text = ~paste("Effort =", round(b, 3), "hr*km<sup>3</sup>", "Exposure =", eagle_rate , "min" ), hoverinfo = 'text')%>% add_trace(x = ~c(10, 10), y = ~c(min((UQ_F-MN_F)/MN_F), max((UQ_F-MN_F)/MN_F)), type = "scatter", mode = "lines", name = "FWS Minimum<br>Survey Effort", showlegend = FALSE)%>% hide_colorbar()%>% layout(hovermode = 'closest', font = list(color = 'black'), xaxis = list(title = "Survey Effort (hr*km<sup>3</sup>)" ), yaxis = list(title = "Precision ((80th% - Mean)/Mean)") ) ) ), fluidRow( column(12, p(class = "caption", tags$b("Figure 6. Estimates of eagle fatality become more similar as estimated with and without priors with greater survey effort."), " Changes in a) difference between fatalities predicted with and without prior exposure information, and b) the size of the 80% CI of the predicted fatality distribution, relative to the mean, as a funciton of survey effort.") ) ) )
Harder to Prove Big Projects Take No Eagles
A bayesian approach to estimating fatalities guards against situations where surveys fail to detect eagles that are in fact present around the site. However, this means that sites where there are truly no birds will have some predicted take. In situations where no eagles are observed during pre-construciton surveys, greater survey effort moves predicted take closer to zero (Fig. 7). This response is highly dependent on the size of the wind farm. For the smallest sites, 48 hr.km^2^ of survey effort reduced fatality predictions < 1. Large sites with no eagle observations require more survey effort to reduce predicted take to zero (~300 hr.km^2^), however the rate of reduction is greater at large sites than at smaller sites.
#PLOT FOR ZERO OBSERVED EAGLES fluidPage( fluidRow( plot_ly(data = zerosim, x = ~Effort, y = ~UQ_F, color = ~Hectares, type = 'scatter', mode = 'markers', marker = list(colorbar = list(x = 0.7, y = 1, title = "Project<br>Size (ha)")), hoverinfo = 'text', text = ~paste("Effort:", Effort, "hr*km<sup>3</sup>", "<br>Project Size:", Hectares, "ha", "<br>Fatalities:", round(UQ_F, 1), "eagles"))%>% colorbar(x = 0.8, y = 1, title = 'Project Size<br>(# Turbines)')%>% layout(legend = list(x = 0.7, y = 1), xaxis = list(title = "Effort (hr*km2)"), yaxis = list(title = "Predicted Annual Eagle Fatalities") ) ), fluidRow( column(12, p(class = 'caption', tags$b("Figure 7. Project size affects predicted take at sites where no eagles are observed."), " Changes in predicted fatalities, given no observed eagles, as a function of survey effort for different sized wind farms.") ) ) )
Methods
We generated posterior distributions for predicted eagle fatalities from empirical and simulation data using the FWS Bayesian model. The prior distribution on eagle exposure used by FWS is a gamma distribution Exposure ~ Gamma (0.968, 0.552), with shape and rate parameters represented by eagle flight time in minutes, and survey effort in (hr km^3^). The prior distribution on collisions used by FWS is a beta distribution Collision ~ Beta (2.31, 396.69) with parameters representative of number of collisions and minutes of exposure.
Empirical Data
We used wind energy site survey data from Appendix A in Bay et al. (2016) to obtain example survey values. This data recorded plot area in hectares, and observation time in minutes. We converted these to survey effort consistent with the units used by FWS (hr km^3^) by rescaling plot area (ha) by 0.2 km/ha, and multiplying by the reported observation time in hours.
Simulation Data
We generated hypothetical values for survey effort based on the minimum requirements provided by U.S. FWS for pre-construction monitoring. FWS requires at least one cylindrical survey plot with radius >= 800m and height >= 200m, and that plots be surveyed for at least 12 hours per year, for two years. Thus, the minimum values for survey area and time were r min(area) km3, and r min(time) hrs, respectively. We simulated up to five plots (area = r max(area)) in increments of one, and up to r max(time) hrs in increments of 12, and calculated survey effort for all combinations.
Observed eagle flight time (min) provided in Bay et al. (2016) is a function of survey effort and eagle activity at a site. Therefore, to generate a range of potential exposure values (min/hr*km^3^), we divided flight time by effort at each survey site to obtain an effort invariant measure of eagle exposure rate, and took a random sample of 20 values. The final simulation data included all combinations of survey effort and eagle exposure rate, which we multiplied to obtain flight minutes, producing 1000 values.
flight <- sample(BAY16$MIN/Bay16$EFFORT, 20) survey_time <- seq(1, 10, 1)*12*2 survey_area <- seq(0.402, 2.01, 0.402) df <- expand.grid(TIME = time, AREA = area, eagle_rate = flight) df$beta <- df$TIME*df$AREA df$alpha <- df$eagle_rate*df$beta
Model Output
We generated eagle fatality predictions using only simulated or empirical values for observed eagle flight time (min) and survey effort (hr km3), and by updating the prior probability of exposure using these values. We obtained the posterior exposure distribution as in New et al. (2015) Exposure ~ Gamma(a + minutes, b + effort), where a and b are the shape and rate parameters of the prior exposure distribution.
# simFatal funciton from New et al. (2013) modified to perform multiple iterations. # BMin: observed number of bird minutes # Fatal: annual avian fatalities on an operational wind facility # SmpHrKm: total time and area surveyed for bird minutes # ExpFac: expansion factor # aPriExp: alpha parameter for the prior on lambda # bPriExp: beta parameter for the prior on lambda # aPriCPr: alpha parameter for the prior on C # bPriCPr: beta parameter for the prior on C # iters: number of sampling iterations simFatal <- function(BMin=-1, Fatal=-1, SmpHrKm, ExpFac, aPriExp=1, bPriExp=1,aPriCPr=1, bPriCPr=1, iters){ out <- data.frame(collision = rep(NA,iters), expose = rep(NA, iters), fatality = rep(NA, iters) ) # The default of a negative value for BMin or Fatal indicates no data if(BMin>=0){ aPostExp <- aPriExp + BMin bPostExp <- bPriExp + SmpHrKm }else{ aPostExp <- aPriExp bPostExp <- bPriExp} # Update the collisions prior if(Fatal>=0){ aPostCPr <- aPriCPr + Fatal bPostCPr <- ((rvmean(Exp) * ExpFac) - Fatal) + bPriCPr }else{ aPostCPr <- aPriCPr bPostCPr <- bPriCPr} for(i in 1:iters){ Exp <- rgamma(n=1, aPostExp, bPostExp) CPr <- rbeta(n=1, aPostCPr, bPostCPr) Fatalities <- ExpFac * Exp * CPr out[i,] <- c(CPr, Exp, Fatalities) } return(out) } estimates <- function(iters, a, b){ out <- simFatal(BMin = a, Fatal = -1, SmpHrKm = b, ExpFac = mean(Bay16$SCALE), aPriExp = 11.81641, bPriExp = 9.7656250, aPriCPr = 1.638029, bPriCPr = 290.0193, iters = iters) fatality <- mean(out$fatality) q80 <- quantile(out$fatality, c(0.8)) out2 <- simFatal(BMin = a, Fatal = -1, SmpHrKm = b, ExpFac = mean(Bay16$SCALE), aPriExp = 0, bPriExp = 0, aPriCPr = 1.638029, bPriCPr = 290.0193, iters = iters) fatality2 <- mean(out2$fatality) q82 <- quantile(out2$fatality, 0.8) return (c("MN_F" = fatality, "U_F" = q80, "MN" = fatality2, "U" = q82)) }
We used an iterative sampling procedure to produce a posterior distribution of annual eagle fatalities. During each iteration, we drew a random sample from the updated exposure distribution and the prior collision probability distribution, and multiplied these values together with a site specific expansion factor. For simulated values, we used the mean expansion factor from the empirical data. For each site and set of simulated values, we performed 100000 iterations of this sampling procedure, and calculated the mean and 80th percentile values of the resulting distribution.
#Predictions integrating Exposure prior and survey data sim <- plyr::mdply(df[, c(5, 4)], estimates, niters = 100000) #Predictions from survey data emp <- vapply(1:nrow(Bay16), function(x){ estimates(100000, Bay16$FLIGHT_MIN[x], BAY16$EFFORT[x]) return (c(fatality, q80, fatality2, q82)) }, USE.NAMES = FALSE, FUN.VALUE = c(0,0,0,0)) )
fluidPage( fluidRow( column(2), column(8, div( HTML('<a rel="license" href="http://creativecommons.org/licenses/by-nc-sa/4.0/"> <img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-nc-sa/4.0/88x31.png" /></a> <br /> This <span xmlns:dct="http://purl.org/dc/terms/" href="http://purl.org/dc/dcmitype/InteractiveResource" rel="dct:type">work</span> by <a xmlns:cc="http://creativecommons.org/ns" href="http://defenders.org" property="cc:attributionName" rel="cc:attributionURL">Defenders of Wildlife</a> is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-nc-sa/4.0/">Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License</a>. <br />'), style = "text-align: center" ) ), column(2) ) )
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.