In tmcd82070/EoAR: Evidence of Absence Regression (EoAR)

In this vignette, I demonstrate EoAR models fitted to no targets, but I include an offset.

My example data is the bat mortality dataset from the EoAR manuscript.

library(tidyverse)
library(EoAR)
data("batMortData", package="EoAR")
batMortData <- batMortData %>%
  mutate(NLBA = 0) # no carcasses found
set.seed(82001)
#knitr::kable(batMortData)

G distribution

# Compute one average g 
gParams <- EoAR::mixtureBeta(
  batMortData$gAlpha, 
  batMortData$gBeta, 
  batMortData$Turbines)

# increase width of g distribution to illustrate the point
# If we don't do this, narrowness of g makes changes 
# in uncertainty difficult to see
gParams$alpha <- gParams$alpha / 1000
gParams$beta <- gParams$beta / 1000
gParams

One Site, One Turbine

# Input dataset
zeroData <- data.frame(
  NLBA = rep(0,1),
  Turbines = 2000
)

fit <- EoAR::eoar(NLBA ~ 1,
                       beta.params = gParams, 
                       data = zeroData, 
                       offset = NULL, 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)
r1 <- data.frame(
  nSites = nrow(zeroData),
  offset = 1,
  model = paste("Y ~",as.character(fit$call$lambda)[3]), 
  mMedian = fit$estimates[labels(fit,"Mtot"),"Estimate"],
  mMean = mean(unlist(fit$out[,"Mtot"])),
  lambdaMean = mean(unlist(fit$out[,"lambda"]))
)

summary(fit)

One Site, 2000 Turbines

fit <- EoAR::eoar(NLBA ~ 1,
                       beta.params = gParams, 
                       data = zeroData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)
r2 <- data.frame(
  nSites = nrow(zeroData),
  offset = zeroData$Turbines,
  model = paste("Y ~",as.character(fit$call$lambda)[3]), 
  mMedian = fit$estimates[labels(fit,"Mtot"),"Estimate"],
  mMean = mean(unlist(fit$out[,"Mtot"])),
  lambdaMean = mean(unlist(fit$out[,"lambda"]))
)

summary(fit)

Five Sites Averaging 400 Turbines

# Input dataset
zeroData <- data.frame(
  NLBA = rep(0,4),
  Turbines = rpois(4, 400)
)
zeroData <- zeroData %>% 
  bind_rows(data.frame(NLBA = 0, Turbines = 2000 - sum(.$Turbines)))
knitr::kable(zeroData)

fit <- EoAR::eoar(NLBA ~ 1,
                       beta.params = gParams, 
                       data = zeroData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)
r3 <- data.frame(
  nSites = nrow(zeroData),
  offset = mean(zeroData$Turbines),
  model = paste("Y ~",as.character(fit$call$lambda)[3]), 
  mMedian = fit$estimates[labels(fit,"Mtot"),"Estimate"],
  mMean = mean(unlist(fit$out[,"Mtot"])),
  lambdaMean = mean(unlist(fit$out[,"lambda[1]"]))
)

summary(fit)

Twenty Sites Averaging 100 Turbines

# Input dataset
zeroData <- data.frame(
  NLBA = rep(0,19),
  Turbines = rpois(19, 2000/20)
)
zeroData <- zeroData %>% 
  bind_rows(data.frame(NLBA = 0, Turbines = 2000 - sum(.$Turbines)))
knitr::kable(zeroData)

fit <- EoAR::eoar(NLBA ~ 1,
                       beta.params = gParams, 
                       data = zeroData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)
r4 <- data.frame(
  nSites = nrow(zeroData),
  offset = mean(zeroData$Turbines),
  model = paste("Y ~",as.character(fit$call$lambda)[3]), 
  mMedian = fit$estimates[labels(fit,"Mtot"),"Estimate"],
  mMean = mean(unlist(fit$out[,"Mtot"])),
  lambdaMean = mean(unlist(fit$out[,"lambda[1]"]))
)

summary(fit)

Final Table

ans <- rbind(r1, r2, r3, r4)
ans <- ans %>% 
  mutate(lambdaScaled = lambdaMean * nSites * offset)
knitr::kable(ans)

Note that lambdaMean is mortalities per turbine per year, and is constant across rows with the same number of turbines. Multiplying mortatlities per turbine by number of turbines is approximately the total number killed (i.e., lambdaScaled = $Mtot = T\lambda$ where $T$ = number of turbines).

tmcd82070/EoAR documentation built on July 13, 2021, 5:52 p.m.