(NOTE: this was originally written to look at the paper it discusses. I am putting it up now after I made a comment on BlueSky about the method to quickly fit different priors using importance smapling)
Controlling populations of carnivores, such as wolves, is a contentious issue. One strategy is to cull populations, but the numbers needed to make a cull effective can be difficult to predict. This is made more difficult by illegal hunting, the effect of which is, by its own nature, difficult to quantify. One argument that has ben used in favour of legal culls is that they reduce the amount of illegal culling. In a recent paper, Chapron and Treves (2016a) used data from two US states to ask if this argument was backed up by evidence, i.e. does culling increase the base growth rate (once the effects of the legal cull itself is accounted for). Their conclusion was that culling actually increased poaching. However, I would suggest that this was an over-confident conclusion from the analysis.
The data used in the study (Chapron and Treves 2016b) are plotted in Fig. 1. We can see the increase in wolf populations, and a suggestion that the growth rate has decreased (the data are plotted on the log scale, so a constant growrh rate would equate to a straight line). The measure of culling that Chapron and Treves (2016a) used was the proportion of the year that culling was allowed. Legal culling was allowed in teh second half of the study period, so it is possible that any negative effect is related to other trends in the environmental conditions that have changed over time.
Even if we are accept that the culling efect is biologically reasonable, the size and direction is not as certain as the paper’s authors suggest. From their own results (their Fig. 2), the effect on growth was -0.03, with a confidence interval of (-0.1, 0.04). This suggests a far weaker effect: the authors’ prior assumption was that an increase or decrease was equally likely, and the posterior probability that the effect was negative was 0.87. This could mean that any effect is small, or that whilst there is a large effect, it is poorly estimated. We can get some idea of the effect: the population in Wisconsin increases by about 17% each year. So if in one yea the population is 500 wolves, in the following year it will be about 587 wolves. If culling were allowed for a whole year, the number of wolves could be between 531 and 611 animals (just using the 95% confidence interval for the culling effect, ignoring all of the sources of uncertainty). The difference is thus about 80 wolves, i.e. the growth rate could be either less than half, or a quarter as large.
Here I look at this in more detail, and ask if there is any estimable effect of culling periods on population growth rates. I slightly adapt the model of Chapron and Treves (2016a), to more explicitly incorporate the fact that the population estimates are given as minima and maxima, and then use more focussed Bayesian approaches to estimate the strengths of evidence for the different hypotheses.
First, we need to read in the data and merge a couple of states.
library(runjags)
library(plyr)
library(RColorBrewer)
library(rdryad)
File <- dryad_download("10.5061/dryad.b7d7v")
load(File$`10.5061/dryad.b7d7v`)
Data <- merge(merge(wisconsin, michigan, by="year", suffixes=c(".WI", ".MI")),
treatment, by="year")These are some function we will use later. The most interesting one is probably BF.prior(), which calculates Bayes Factors for different priors
# Function to simulate intial values for model
# Arguments
# dat - list of data send to JAGS/BUGS
# GVS - should inits for GVS model be generated? Defaults to FALSE
# JAGS - should inits for a JAGS model (rather than BUGS) be generated? Defaults to FALSE
SimInits <- function(dat, GVS=FALSE, JAGS=FALSE) {
init <- list(
Diff = runif(1,4,10),
beta0r = rnorm(dat$NStates,0,1),
gamma = rgamma(1,50,50),
muN = cbind(dat$Intervals[,1,1]+1E-3,matrix(NA, nrow = dat$NStates, ncol = dat$MaxYear-1)),
# sdN = rgamma(dat$NStates, 10,100)
sdN = rgamma(1, 10,100)
)
if(GVS) {
init$beta1rStar = rnorm(1,0,0.1)
init$Ind1r <- 2
} else {
init$beta1r = rnorm(1,0,0.1)
}
if(JAGS) {
lnNMin <- ifelse(is.na(dat$Intervals[,,1]), 6.4,dat$Intervals[,,1])
init$Intervals <- ifelse(is.na(dat$Intervals), 6.4,NA)
lnNMax <- ifelse(is.na(dat$Intervals[,,2]), lnNMin+0.9*init$Diff,
(dat$Intervals[,,1]+dat$Intervals[,,2])/2)
# lnNMax <- lnNMin+0.9*init$Diff
init$Intervals[,,2] <- ifelse(is.na(dat$Intervals[,,2]), lnNMax,NA)
init$logN <- lnNMin+0.5*(lnNMax-lnNMin)
} else {
init$lnNMin = ifelse(is.na(dat$lnNMin), 6.4,NA)
init$logN <- lnNMin+0.5*(lnNMax-lnNMin)
lnNMin <- ifelse(is.na(dat$lnNMin), 6.4,dat$lnNMin)
lnNMax <- ifelse(is.na(dat$lnNMax), lnNMin+0.9*init$Diff,dat$lnNMax)
init$lnNMax = ifelse(is.na(dat$lnNMax), lnNMax,NA)
}
init
}
# Function to calculate posterior mode, and 50% and 95% HPDIs
# Arguments
# mcmc - mcmc object, or something that can be coerced to one (e.g. a data frame)
# ... - addition argumetns to posterior.mode().
GetHPD=function(mcmc,...) {
require(MCMCglmm)
require(coda)
if(!is.mcmc(mcmc)) mcmc=as.mcmc(mcmc)
hpd=cbind(posterior.mode(mcmc,...), HPDinterval(mcmc, prob=0.5), HPDinterval(mcmc, prob=0.95))
colnames(hpd)=c("Mode", "L50", "U50", "L95", "U95")
hpd
}
# Function to re-calculate Bayes factors for different priors, using importance weights
# Arguments
# newSD - new prior standard deviation
# df - data frame of posteriors, from prior model
# priorP - prior probability (defaults to 0.5)
# newM - new prior mean (defaults to 0)
# origM - original prior mean (defaults to 0)
# origSD - original prior standard deviation (defaults to 1)
BF.prior <- function(newSD, df, priorP=0.5, newM = 0, origM = 0, origSD=1) {
Wt <- ifelse(df$Ind1r==1,1, dnorm(df$beta1r, newM,sd=newSD)/dnorm(df$beta1r, origM, sd=origSD))
Wt <- Wt/sum(Wt)
Prob <- sum(Wt[df$Ind1r!=1])
Prob*(1-priorP)/(priorP*(1-Prob))
}And the JAGS code:
cat("model {
# dinterval use from here:
# http://stats.stackexchange.com/questions/13847/how-does-dinterval-for-interval-censored-data-work-in-jags
for(s in 1:NStates) {
muN[s,1] ~ dnorm(0, 1.0E-2)
for(t in 2:MaxYear) {
muN[s,t] <- log(exp((logN[s,t-1] + r[s,t])) - gamma*H[s,t]) # dynamics
}
for(t in 1:MaxYear) {
AllIn[s,t] ~ dinterval(logN[s,t], Intervals[s,t,]) # censoring
logN[s,t] ~ dnorm(muN[s,t], tauN) # uncensored model
r[s,t] <- beta0r[s] + beta1r*(D[s,t]-mean(D[s,])) # growth rate
# for missing data
MaxN[s,t] <- Intervals[s,t,1] + Diff
Intervals[s,t,2] ~ dunif(Intervals[s,t,1], MaxN[s,t])
Intervals[s,t,1] ~ dnorm(0,1.0E-2)
}
beta0r[s] ~ dnorm(0,1.0E-6)
}
sdN ~ dunif(0,0.5);
tauN <- pow(sdN, -2)
# GVS
Pr[1] <- 0.5; Pr[2] <- 1-Pr[1]
Tau[1] <- 1.0E2; Tau[2] <- 1
Ind1r ~ dcat(Pr[])
beta1rStar ~ dnorm(0,Tau[Ind1r])
beta1r <- beta1rStar*(Ind1r-1)
gamma ~ dnorm(1.06,14) # Effect of culled animals
Diff ~ dunif(1,10)
}", file="WolvesModelGVSjags.txt")Now we can plot the data:
par(mfrow=c(1,1), mar=c(4.1,4.1,1,1), oma=rep(0,4))
plot(wisconsin$year, log(wisconsin$Nobs.Min), type="n", xlab="Year", ylab="log(population size)", main="")
apply(treatment, 1, function(yr, ylims) {
rect(yr["year"]-0.5, ylims[1], yr["year"]+0.5, ylims[2],
col=grey(1-yr["days_cull_WI"]/500), border=NA)
}, ylims=c(0,7) )
polygon(c(wisconsin$year, rev(wisconsin$year)), c(log(wisconsin$Nobs.Min), rev(log(wisconsin$Nobs.Max))),
col=3, border=NA)
lines(michigan$year, log(michigan$Nobs), col="blue", lwd=1.4)
box(bty="o")
text(1999.5,5.5, "Wisconsin", adj=1)
text(1999,5.0, "Michigan", adj=0)
The same process model is used, i.e. we assume
\[ e^{\mu_t^S} = E(N^S_{t}) = N^S_{t-1} e^{r_t^S} - \gamma H^S_t \] \[ N_t^S \sim N(\mu_t^S, \sigma^{2S}_{proc}) \] \[ r_t^S = \beta_0^{rS} + \beta_1^{r} (D_t^S - \bar{D^S}) \]
However, we use a different observation model. For Wisconsin we have maximum and minimum population sizes (\(N_t^{U1}\) and \(N_t^{L1}\)), and for Michigan we have minimum population sizes (\(N_t^{U2}\)). We therefore impute the maximum population size for Wisconsin by assuming the ceiling for the maximum is proportional to the minimum population size, i.e. \(e^C N_t^{U2}\), and the actual maximum is uniformally distbuted on the log scale between the minimum and this ceiling, i.e.
\[ log(N_t^{U2}) \sim U(log(N_t^{L2}), log(N_t^{L2})+C) \]
We then assume that \(N_t^S\) is interval censored between \(N_t^{LS}\) and \(N_t^{US}\)).
Where possible I assumed the same priors as following priors as Chapron and Treves (2016a), except for the following variables, which were either not stated in the original paper, or correspond to variables not in their model.
The effects of legal culling are the main focus here, and the question is whether there is an effect, and if there is which direction is it in. Thus, I use a variable selection technique called GVS (Dellaportas, Forster, and Ntzoufras 2000). This makes the prior a mixture of a “spike” at 0 and a “slab” with a normal prior. This is constructed by using an indicator variable, \(I^{r}\), which equals 1 if the culling effect is “in” the model, otherwise it is zero:
\[I^{r} \sim Bern(p_I) \] \[\beta_1^{r} = I^{r} \beta_1^{*}\] \[\beta_1^{*} \sim N(0, \sigma^2_{\beta I}) \]
Here I use \(p_i=0.5\), to be equivocal about whether the effect is zero or not. The prior on \(\beta_1^{*}\), is a mixture distribution. If there is a culling effect (i.e. if \(I^{r} = 1\)) I use a vague prior distribution. If \(I^{r} = 0\) the prior for \(\beta_1^{*}\) does not affect the model, so we could use any prior we want (Carlin and Chib 1995): it is called a pseudo-prior. The pseudo-prior can affect the mixing of the MCMC, so is set to something close to the posterior. Here we use a prior of \(\sigma^2_{\beta 1}=1\) and a pseuso-prior of \(\sigma^2_{\beta 1}=0.01\). Variable selection is genrally sensitive to the prior, so an uninformative prior may not be a good choice as it will penalise the parameter too much.
The posterior probability that \(I^{r}=1\) (i.e. that there is an effect of legislation) is esimated as the posterior probability from the MCMC (it is the propotion of draws for which \(I^{r}=1\)): we call this \(p_I^p\). The Bayes Factor then summarises how the evidene for a hypothesis is changed by the data (e.g. Kass and Raftery 1995):
\[B_{01} = \frac{p_I^p}{1 - p_I^p}\frac{1 - p_I}{p_I}\]
i.e. the posterior odds divided by the prior odds. This removes any effect of the prior for \(p_i\). A Bayes Factor fractor than 1 means that there is evidence (however weak) in favour of the variable having an effect. We can also calculate a probability that the parameter is positive, given it is non-zero, i.e. \(Pr(\beta_1^{r} >0 | I^{r}=1)\), and thus the Bayes Factor from this.
Because Bayes factors are known to be sensitive to prior distributions (in particular the prior for the variables that the varaible selection is applied to), I re-esimated the Bayes factor for prior distributions of \(\beta_1^{r}\) with the prior standard deviation varied between 0.05 and 5. A “legislator’s prior” was also fitted, which assumes that the legislator’s prior mean is that an annual change increases the log growth rate by 0.05, with a probability of 0.9 that there is a positive effect, which equates to a prior standard deviation of 0.039.
The model was fitted with JAGS through the runjags package Denwood (2016) in R version 4.5.1. Four chains were run in parallel, and after a burn-in of \(10^4\) iterations, a further \(10^4\) iterations were run and thinned to every fifth iteration, to give \(4 \times 10^4\) iterations in total.
Intervals <- array(c(log(Data$Nobs.Min), log(Data$Nobs),
log(Data$Nobs.Max), rep(NA, nrow(Data))), dim=c(nrow(Data),2,2),
dimnames=list(Year=Data$year,
State=c("Wisconsin" , "Michigan"),
Limit=c("Lower", "Upper")))
DataToJAGS <- list(
NStates = 2,
MaxYear = nrow(Data),
Intervals = aperm(Intervals, c(2,1,3)),
AllIn = matrix(1, ncol=nrow(Data), nrow=2),
H = rbind(Data$H.WI, Data$H.MI),
D = rbind(Data$days_cull_WI, Data$days_cull_MI)/365
)
NChains <- 4
Inits <- sapply(1:NChains, function(chain) {
init <- SimInits(DataToJAGS, GVS=TRUE, JAGS = TRUE)
init$.RNG.name <- "base::Mersenne-Twister"
init$.RNG.seed <- chain
init
}, simplify=FALSE)
Pars <- c("Diff", "beta0r", "beta1r", "gamma", "sdN", "Ind1r")
Wolves.mcmc <- run.jags("WolvesModelGVSjags.txt", data = DataToJAGS, inits=Inits, monitor=Pars,
n.chains = NChains, adapt = 1e3, burnin=9e3, sample=2e5, thin=1, silent.jags=TRUE,
method="parallel")
Wolves.df <- ldply(lapply(Wolves.mcmc$mcmc, as.matrix))
# This code looked at if there was an effect of time on population growth, but none was detected.
# The changes in the JAGS code were for muN, which was changed to:
# muN[s,t] <- log(exp((logN[s,t-1] + r[s,t]) - gamma*H[s,t])) - betat*Time[t] # dynamics
# and a prior for betat was added:
# betat ~ dnorm(0,1.0E1)
# InitsT <- lapply(Inits, function(lst) { lst$betat<-0; lst})
# Pars <- c(Pars, "betat");
# DataToJAGST <- DataToJAGS;
# DataToJAGST$Time <- 1:DataToJAGST$MaxYear
# # DataToJAGST$betat <- 0
# WolvesT.mcmc <- run.jags("WolvesModelGVSjagsQuadTime.txt", data = DataToJAGST,inits=Inits, monitor=Pars,
# n.chains = NChains, adapt = 1e3, burnin=9e3, sample=2e5, thin=1, silent.jags=TRUE,
# # n.chains = 2, adapt = 1e3, burnin=4e3, sample=1e3, thin=1, silent.jags=TRUE,
# method="parallel")
# Wolves.df <- ldply(lapply(WolvesT.mcmc$mcmc, as.matrix))
# GetHPD(Wolves.df)
#
# plot(Wolves.df$beta1r[seq(1,nrow(Wolves.df), by=100)], Wolves.df$betat[seq(1,nrow(Wolves.df), by=100)], cex=0.1)
# cor.test(Wolves.df$beta1r[Wolves.df$Ind1r==2], Wolves.df$betat[Wolves.df$Ind1r==2])Models to asses prior sensitivity were fitted using importance weights (e.g. Robert and Casella (2009), Chapter 3) to re-weight the original model (rather than running the MCMC multiple times). For a model with prior \(p(\beta_1^{ r(i) }|M_2)\), each draws from the original MCMC (with prior \(p(\beta_1^{ r(i) }|M_1\)) was given a weight
\[ w^{(i)} = \frac{p(\beta_1^{ r(i) }|M_2)}{p(\beta_1^{ r(i) }|M_1)} \]
so that the new posterior probability of an effect of legislation is
\[ p_I^p|M_2 = \frac{ \sum_{i=1}^n {w^{(i)} I^{r(i)}}}{\sum_{i=1}^n {w^{(i)}}} \]
We need a bit of code to adjust the model for the different priors:
# Bayes factors
PriorProb <- 0.5
PostProb <- mean(Wolves.df$Ind1r==2)
BayesFactor <- PostProb*(1-PriorProb)/(PriorProb*(1-PostProb))
PrPos <- mean(Wolves.df$beta1r[Wolves.df$Ind1r==2]>0)
PrPos.BF <- PrPos*(1-PriorProb)/(PriorProb*(1-PrPos))
###################
# Prior Sensitivity
# Conditional effect of legislation
beta1r.1 <- Wolves.df$beta1r[Wolves.df$Ind1r==2]
# Original prior
Wt.O <- dnorm(beta1r.1, 0,sd=1e3)/dnorm(beta1r.1, 0, sd=1)
PostProb.O <- sum(Wt.O)/(sum(Wt.O)+sum(Wolves.df$Ind1r==1))
BayesFactor.O <- PostProb.O*(1-PostProb.O)/(PriorProb*(1-PostProb))
Beta.O <- sample(x=beta1r.1, size=1e5, prob=Wt.O, replace=TRUE)
# Legislator's prior
Wt <- ifelse(Wolves.df$Ind1r==1,1,
dnorm(Wolves.df$beta1r, 0.05,sd=0.05/qnorm(0.9))/dnorm(Wolves.df$beta1r, 0, sd=1))
Wt.L <- Wt[Wolves.df$Ind1r==2]
beta1r.L <- Wolves.df$beta1r[Wolves.df$Ind1r==2]
Prob.L <- sum(Wt[Wolves.df$Ind1r!=1])/sum(Wt)
ProbPos.L <- sum(Wt.L[beta1r.L>0])/sum(Wt.L)
LegislatorsBeta <- sample(x=beta1r.L, size=1e5, prob=Wt.L, replace=TRUE)
# LegislatorsBeta.hpd <- GetHPD(LegislatorsBeta)
# Priors with different prior standard deviation for beta1r
PriorSD <- seq(0.1,5, length=30)
BFs <- sapply(PriorSD, BF.prior, df=Wolves.df, priorP=0.5)
# Now try different prior probabilities for +/-
PriorPrPlus <- 0.9
PostPrPlus <- mean(Wolves.df$beta1r[Wolves.df$Ind1r==2]>0)
BFPlus <- PostPrPlus/(1-PostPrPlus)
PostPrPlus <- PriorPrPlus*BFPlus/(1 + PriorPrPlus*(BFPlus-1))Note that in general it is ptobably better to use the widest priors for the MCMC, and importance sampling on the narrower priors, so that the importance smaples cover the range of their posterior.
Overall the evidence suggests that the culling period has no effect: the Bayes Factor in favour of no effect is 14.2, which corresponds to “Positive” evidence (but not “strong”) according to the guidelines in Kass and Raftery (1995). If we look at the estimates conditional on it being non-zero, there is a probability of 0.24, which equates to a Bayes factor of 3.2 in favour of the effect being negative: this is also “positive”, but only just (the cut-off is about 3.2). Thus the data suggest we should be doubtful about whether there is an effect, and even if we were to decide there is one, we are doubtful about the direction.
Norm <- function(x) dnorm(x, mean=0, sd=1e3)
NormG <- function(x) dnorm(x, mean=1.06, sd=sqrt(1/14))
# GamPrior <- curve(NormG, from=0, to=4, add=TRUE, col=2)
Densities <- apply(Wolves.df[,c("beta0r[1]", "beta0r[2]", "gamma", "sdN")], 2, density)
par(mfrow=c(2,2), mar=c(4,2,3,1), oma=c(0,2,0,0))
plot(Densities[["beta0r[1]"]], xlab=expression(beta[0]^{r1}), ylab="", main="Baseline Growth for Wisconsin",
ylim=c(0, 1.04*max(Densities[["beta0r[1]"]]$y)), yaxs="i")
curve(Norm, add=TRUE, col=2)
plot(Densities[["beta0r[2]"]], xlab=expression(beta[0]^{r2}), ylab="", main="Baseline Growth for Michigan",
ylim=c(0, 1.04*max(Densities[["beta0r[2]"]]$y)), yaxs="i")
curve(Norm, add=TRUE, col=2) # from = 0, to = 0.3,
plot(Densities[["gamma"]], xlab=expression(gamma), ylab="", main="Effect of Legal Culling",
ylim=c(0, 1.04*max(Densities[["gamma"]]$y)), yaxs="i")
curve(NormG, add=TRUE, col=2)
# lines(GamPrior, col=2)
plot(Densities[["sdN"]], xlab=expression(sigma[proc]), ylab="", main="Process Standard Deviation",
ylim=c(0, 1.04*max(Densities[["sdN"]]$y)), yaxs="i")
abline(h=1/0.5, col=2)
legend(0.155,21, c("Prior", "Posterior"), lty=1, col=c(1,2))
mtext("Density",2,outer=TRUE)
The estimates of the other parameters are plotted in Fig. 2. The baseline growth are very similar to those given by Chapron and Treves (2016a), but the effects of legal culling and the standard deviation of the population dynamics both have wider posterior distributions: the legal culling mostly replicates the prior distribution.
Cols <- brewer.pal(3, "Dark2")
Norm1 <- function(x) dnorm(x, mean=0, sd=1)
NormO <- function(x) dnorm(x, mean=0, sd=1e3)
NormL <- function(x) dnorm(x, mean=0.05,sd=0.05/qnorm(0.9))
par(mar=c(4.1,4.1,1,1), las=1, lwd=1.3)
plot(density(LegislatorsBeta), type="l", col=Cols[3], xlim=c(-0.2,0.15),
xlab=expression(beta[1]^{r1}), ylab="Density", main="")
curve(NormL, col=Cols[3], lty=2, add=TRUE)
lines(density(beta1r.1), col=Cols[2])
curve(Norm1, col=Cols[2], lty=2, add=TRUE)
lines(density(Beta.O), col=Cols[1])
curve(NormO, col=Cols[1], lty=2, add=TRUE)
legend(-0.19,12, c("Original", "GVS", "Legislator's"), col=Cols, lty=1, cex=0.7)
When we vary the prior standard deviation, we see that the Bayes Factor stays below 1, and becomes smaller as the prior variance increases (Fig. 3). This suggests there is still poor evidence for any effect. The Legislator’s Prior leads to a Bayes Factor of 0, and posterior probability that the effect is positive - given there is an effect - of 0.73.
If we look at the posteror densities, conditional on an effect, we see that the legislator’s prior (which assumes a positive effect) leads to a posterior distribution shrunk towards zero, whilst the prior used here and the original prior from Chapron and Treves (2016a) are almost identical.
plot(PriorSD, BFs, type="n", xlab="Prior Standard Deviation", ylab="Bayes Factor", xaxs="i", yaxs="i",
ylim=c(0,0.6), xlim=c(0,5), las=1)
rect(0,0, 5, 1/3, border=NA, col="grey70")
rect(0,0, 5, 1/20, border=NA, col="grey50")
lines(PriorSD, BFs,col=2)
text(3.9,0.2, "Positive", adj=c(0,0))
text(3.9,0.03, "Strong", adj=c(0,0))
box(bty="o")
The analyses suggest that there is little evidence for any effect of legislation on wolf dynamics, and to the extent the data point to any conclusion other than “not sure”, it suggests that there is no effect. Any support for the suggestion that allowing culling increases poaching of a large carnivore is weak, although it is stronger than the support for a positive effect. But the overall message is that we cannot be sure.
In an applied context, these result are discouraging. Capron & Treves asked an important question, but unfortunately the data were not amenable to giving a definite answer. Unfortunately this is not uncommon in ecology, where the processes that drive the dynamics of populations are noisy, and the data that is collected is often collected over too short a time scale or suffers from confounding of variables. These problems are often inevitable. Thus whilst those analysign these data can do as much as possible with it, it shoul not be surprising if we are not always able to give definitive answer to what are, after all, important questions.