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

In this vignette, we demonstrate EoAR models that contain different $g$ calculations and offsets. We estimate models with a common $g$ across all sites, and compare them to models where $g$ varies by site. We include an offset, and compare results to those from a model without an offset.

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

LBBA Analysis with Common G's and Offset

The following replicates two EoAR models reported in the EoAR manuscript. We replicate the (Intercept-only) and the top LBBA model. The Manuscript analyzes all used a common g distribution across facility-year combinations.

This section replicates results reported Table 2 of the EoAR manuscript within simulation error. Setting the random seed used for the manuscript analyses (i.e., set.seed(33134)) does not replicate manuscript results exactly because all possible models were fitted in a particular order during the manuscript analysis. We only estimate two models here.

set.seed(33134)

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

# For manuscript, we set nburns = 750000, niters = 500*500*5, nthins = 500*5.
# That is way too much sampling for LBBA, but about right for INBA. Here, 
# we cut down the iterations to speed things up.

# LBBA intercept only
lbbaFacG <- EoAR::eoar(LBBA ~ 1,
                       beta.params = gParams, 
                       data = batMortData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)

Dic <- dic.samples(lbbaFacG$jags.model, 
                   n.iter = 100*thin(lbbaFacG$out), 
                   thin = thin(lbbaFacG$out), 
                   progress.bar = "none")

lbba.Intercept <- data.frame(
  gStructure = "common",
  offset = "Turbines",
  structure = as.character(lbbaFacG$call$lambda)[3], 
  DIC = sum(Dic$deviance),
  mHat = lbbaFacG$estimates[labels(lbbaFacG,"Mtot"),"Estimate"],
  mHat.lo = lbbaFacG$intervals[labels(lbbaFacG,"Mtot"),"5%"],
  mHat.hi = lbbaFacG$intervals[labels(lbbaFacG,"Mtot"),"95%"]
)

summary(lbbaFacG)

# LBBA top model
lbbaFacG <- EoAR::eoar(LBBA ~ DistToWater + EastWestStratum + SiteAge,
                       beta.params = gParams, 
                       data = batMortData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)
Dic <- dic.samples(lbbaFacG$jags.model, 
                   n.iter = 100*thin(lbbaFacG$out), 
                   thin = thin(lbbaFacG$out), 
                   progress.bar = "none")

lbba.top <- data.frame(
  gStructure = "common",
  offset = "Turbines",
  structure = as.character(lbbaFacG$call$lambda)[3], 
  DIC = sum(Dic$deviance),
  mHat = lbbaFacG$estimates[labels(lbbaFacG,"Mtot"),"Estimate"],
  mHat.lo = lbbaFacG$intervals[labels(lbbaFacG,"Mtot"),"5%"],
  mHat.hi = lbbaFacG$intervals[labels(lbbaFacG,"Mtot"),"95%"]
)

summary(lbbaFacG)

lbba.comG.turbOffset <- rbind(lbba.top, lbba.Intercept)
print(lbba.comG.turbOffset)

LBBA Analysis with Facility-Specific G's and Offsets

The following code is exactly like that of the previous section except that a different g distribution is used for each facility-year combination.

set.seed(33134)

gParams <- batMortData %>% 
  rename("alpha" = gAlpha,
         "beta" = gBeta) %>% 
  select(alpha, beta)

# intercept only
lbbaFacG <- EoAR::eoar(LBBA ~ 1,
                       beta.params = gParams, 
                       data = batMortData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)

Dic <- dic.samples(lbbaFacG$jags.model, 
                   n.iter = 100*thin(lbbaFacG$out), 
                   thin = thin(lbbaFacG$out), 
                   progress.bar = "none")

lbba.Intercept <- data.frame(
  gStructure = "facilityYear",
  offset = "Turbines",
  structure = as.character(lbbaFacG$call$lambda)[3], 
  DIC = sum(Dic$deviance),
  mHat = lbbaFacG$estimates[labels(lbbaFacG,"Mtot"),"Estimate"],
  mHat.lo = lbbaFacG$intervals[labels(lbbaFacG,"Mtot"),"5%"],
  mHat.hi = lbbaFacG$intervals[labels(lbbaFacG,"Mtot"),"95%"]
)

summary(lbbaFacG)

# LBBA top model
lbbaFacG <- EoAR::eoar(LBBA ~ DistToWater + EastWestStratum + SiteAge,
                       beta.params = gParams, 
                       data = batMortData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)
Dic <- dic.samples(lbbaFacG$jags.model, 
                   n.iter = 100*thin(lbbaFacG$out), 
                   thin = thin(lbbaFacG$out), 
                   progress.bar = "none")

lbba.top <- data.frame(
  gStructure = "facilityYear",
  offset = "Turbines",
  structure = as.character(lbbaFacG$call$lambda)[3], 
  DIC = sum(Dic$deviance),
  mHat = lbbaFacG$estimates[labels(lbbaFacG,"Mtot"),"Estimate"],
  mHat.lo = lbbaFacG$intervals[labels(lbbaFacG,"Mtot"),"5%"],
  mHat.hi = lbbaFacG$intervals[labels(lbbaFacG,"Mtot"),"95%"]
)


summary(lbbaFacG)

lbba.facG.turbOffset <- rbind(lbba.top, lbba.Intercept)
print(lbba.facG.turbOffset)

NLBA Analysis with Facility-Specific G's and Offsets

set.seed(33134)

gParams <- batMortData %>% 
  rename("alpha" = gAlpha,
         "beta" = gBeta) %>% 
  select(alpha, beta)

# intercept only
fit <- EoAR::eoar(NLBA ~ 1,
                       beta.params = gParams, 
                       data = batMortData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)

Dic <- dic.samples(fit$jags.model, 
                   n.iter = 100*thin(fit$out), 
                   thin = thin(fit$out),
                   progress.bar = "none")

nlba.Intercept <- data.frame(
  gStructure = "facilityYear",
  offset = "Turbines",
  structure = as.character(fit$call$lambda)[3], 
  DIC = sum(Dic$deviance),
  mHat = fit$estimates[labels(fit,"Mtot"),"Estimate"],
  mHat.lo = fit$intervals[labels(fit,"Mtot"),"5%"],
  mHat.hi = fit$intervals[labels(fit,"Mtot"),"95%"]
)

summary(fit)

# LBBA top model
fit <- EoAR::eoar(NLBA ~ DistToWater + EastWestStratum + SiteAge,
                       beta.params = gParams, 
                       data = batMortData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)
Dic <- dic.samples(fit$jags.model, 
                   n.iter = 100*thin(fit$out), 
                   thin = thin(fit$out), 
                   progress.bar = "none")

nlba.top <- data.frame(
  gStructure = "facilityYear",
  offset = "Turbines",
  structure = as.character(fit$call$lambda)[3], 
  DIC = sum(Dic$deviance),
  mHat = fit$estimates[labels(fit,"Mtot"),"Estimate"],
  mHat.lo = fit$intervals[labels(fit,"Mtot"),"5%"],
  mHat.hi = fit$intervals[labels(fit,"Mtot"),"95%"]
)

summary(fit)

nlba.facG.turbOffset <- rbind(nlba.top, nlba.Intercept)
print(nlba.facG.turbOffset)

NLBA Analysis with Facility-Specific G's But No Offset

set.seed(33134)

gParams <- batMortData %>% 
  rename("alpha" = gAlpha,
         "beta" = gBeta) %>% 
  select(alpha, beta)

# Fit the EoAR model, intercept only
nlbaFacG <- EoAR::eoar(NLBA ~ 1,
                       beta.params = gParams, 
                       data = batMortData, 
                       offset = NULL, 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)
summary(nlbaFacG)

NLBA Analysis with Common G's and Offsets

set.seed(33134)

gParams <- EoAR::mixtureBeta(batMortData$gAlpha, 
                             batMortData$gBeta, 
                             batMortData$Turbines)

# intercept only
fit <- EoAR::eoar(NLBA ~ 1,
                       beta.params = gParams, 
                       data = batMortData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)

Dic <- dic.samples(fit$jags.model, 
                   n.iter = 100*thin(fit$out), 
                   thin = thin(fit$out),
                   progress.bar = "none")

nlba.Intercept <- data.frame(
  gStructure = "facilityYear",
  offset = "Turbines",
  structure = as.character(fit$call$lambda)[3], 
  DIC = sum(Dic$deviance),
  mHat = fit$estimates[labels(fit,"Mtot"),"Estimate"],
  mHat.lo = fit$intervals[labels(fit,"Mtot"),"5%"],
  mHat.hi = fit$intervals[labels(fit,"Mtot"),"95%"]
)

summary(fit)

# LBBA top model
fit <- EoAR::eoar(NLBA ~ DistToWater + EastWestStratum + SiteAge,
                       beta.params = gParams, 
                       data = batMortData, 
                       offset = log(Turbines), 
                       nadapt = 10000, 
                       nburns = 50000,
                       niters = 50000, 
                       nthins = 100,
                       nchains = 3,
                       quiet = TRUE)
Dic <- dic.samples(fit$jags.model, 
                   n.iter = 100*thin(fit$out), 
                   thin = thin(fit$out), 
                   progress.bar = "none")

nlba.top <- data.frame(
  gStructure = "facilityYear",
  offset = "Turbines",
  structure = as.character(fit$call$lambda)[3], 
  DIC = sum(Dic$deviance),
  mHat = fit$estimates[labels(fit,"Mtot"),"Estimate"],
  mHat.lo = fit$intervals[labels(fit,"Mtot"),"5%"],
  mHat.hi = fit$intervals[labels(fit,"Mtot"),"95%"]
)

summary(fit)

nlba.comG.turbOffset <- rbind(nlba.top, nlba.Intercept)
print(nlba.facG.turbOffset)

Final estimates from all runs

ans <- rbind(lbba.comG.turbOffset, 
             lbba.facG.turbOffset, 
             nlba.comG.turbOffset, 
             nlba.facG.turbOffset)
print(ans)

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