1) Closed Population Model

Our data has three sites (North Shore, Richmond, and Vancouver) which each belong under the VCHA. There are 15 time replications (the years 2000 to 2014). Next we will need to prepare the data for model fitting in “unmarked”. We will be using 1 count = 100 reported cases. This is for computational efficiency, as these models require large computation times which increase with count size. This will reduce our model accuracy, since we will now only be able to report results down to 100’s of cases, however the trade off is that the models will be computationaly feasible on a modern home computer.

# counts
# for computational efficiency, 1 count represented 100 cases
Y <- round(dat[,2:16]/100)
# site covariates
X <- data.frame(Site=as.factor(dat[,1]))
# within site covariates (eg site/sampling occasion specific)
Z <- data.frame(Year=as.factor(rep(2000:2014,3)))

# uframe for closed population/Royle pcount model
uframeC <- unmarkedFramePCount(y=Y,siteCovs=X, obsCovs=Z)
# minimum site population size
Kmin = max(uframeC@y)

kable(as(uframeC, "data.frame"), "html", caption = "unmarkedFramePCount") %>%
  kable_styling()

Once the data has been prepared, it is easy to fit the model using the pcount function in “unmarked”. This fits an N Mixture model with covariates specified for abundance and detection rate. We will not specify any covariates for the first model, which corresponds to a constant population size and constant detection rate over sites and time.

The parameter K is an artificial upper bound on the population abundance estimate, and must be chosen such that increasing K does not increase the abundance estimate. To facilitate this choice, we will fit a sequence of models with increasing K, and choose a value of K large enough that no change to the abundance estimate is visible.

It is important to note that K must be larger than the maximum observed count, since the actual population size cannot be smaller than the observed abundance.

# checking for appropriate K
ks <- seq(1*Kmin+1, 2*Kmin,Kmin/50)
Abun <- numeric(0)
for(t in ks) {
  m <- pcount(~1 ~1, uframeC, K=t, se = FALSE, starts=c(5, -0.5))
  
  r <- ranef(m)
  a <- bup(r)
  Abun <- rbind(Abun, a)
}

# To get abundance we need to multiply by 100 to convert from 1 count to 100 cases.
# We use rowSums to get the total estimate for VCH, rather than the site estimates.
ggplot(data=tibble(ks=ks, Abun=(100*rowSums(Abun))), aes(x=ks, y=Abun)) +
  geom_vline(xintercept = 131, colour="red") +
  geom_line(colour="blue") +
  geom_point() +
  xlab("Upper Bound K") +
  ylab("Abundance Estimate") +
  theme_light() +
  ggtitle("Estimated Abundance vs K")

The plot of estimated abundance as a function of increasing K shows that for K>130 (shown by the red vertical line), the estimate is constant. This is a good indication that the maximum likelihood estimates of the parameters have converged. We can choose to use any K value larger than 130, and as such will use K=200 going forwards.

We are now ready to fit the most basic closed population model, and examine the parameter estimates. This will provide a baseline of comparison for the more complicated open population models.

# fit the model
toymodel <- pcount(~1 ~1, uframeC, K=200, se = TRUE, starts=c(5, -0.5))
s_toy <- summary(toymodel)

## 
## Call:
## pcount(formula = ~1 ~ 1, data = uframeC, K = 200, starts = c(5, 
##     -0.5), se = TRUE)
## 
## Abundance (log-scale):
##  Estimate    SE    z P(>|z|)
##      4.27 0.094 45.4       0
## 
## Detection (logit-scale):
##  Estimate   SE    z P(>|z|)
##     0.545 0.18 3.03 0.00242
## 
## AIC: 344.3374 
## Number of sites: 3
## optim convergence code: 0
## optim iterations: 28 
## Bootstrap iterations: 0

det_prob_toy <- plogis(coef(toymodel, type="det"))
det_conf_toy <- confint(toymodel, type="det", level=0.95)

# get model parameter estimates
r_toy <- ranef(toymodel)
a_toy <- 100*bup(r_toy)

# 95% confidence interval (likelihood based)
conf_toy <- 100*confint(r_toy, level=0.95) 

df_toy <- data.frame(Site=as.factor(c("North Shore","Richmond","Vancouver")), Abundance_Estimate=round(a_toy), Lower2.5=conf_toy[,1], Upper97.5=conf_toy[,2])

kable(df_toy, "html", caption = "Closed Population Model Abundance Estimates") %>%
  kable_styling()

ggplot(data=df_toy) +
  geom_errorbar(aes(x=Site, ymin=Lower2.5, ymax=Upper97.5), color="red", width=0.2) +
  geom_point(aes(y = Abundance_Estimate, x=Site ), color="black", size=1.25) +
  ggtitle("Royle Model (~1 ~1)\nAbundance Estimates By Site, with 95% CI") +
  geom_text(aes(x=Site, y=Abundance_Estimate, label=Abundance_Estimate), hjust=-1) +
  geom_text(aes(x=1, y=max(Abundance_Estimate), label=paste0("95% CI for Detection Probability:\n", signif(det_conf_toy[1],2)," < p=",signif(det_prob_toy,2)," < ",signif(det_conf_toy[2],2)))) +
  theme_light()

This model predicts a CDR of 0.63, and a total population of depressed individuals of 21400. This leads to a prediction on the unobserved depression cases of 7900. Of course, the model assumptions are not met, so these conclusions are merely for illustration of the methods, and not to be taken seriously.

2) Open Population Model

The open population models are more complicated, allowing for population dynamics. The model will be fit using the function pcountOpen in the “unmarked” package.

The model itself is one where population size for site \(i\) at time \(t+1\) is given by the survival \((S_{it})\) and the population recruitment \((G_{it})\) from time \(t\).

\[ N_{i(t+1)} = S_{it}+G_{it} \] \[ S_{it} \sim Binomial(N_{it}, \omega_{it}) \] \[ G_{it} \sim Poisson(\gamma_{it}) \] These combine with the thinned detection process to give observed counts \(y_{it}\): \[ y_{it} \sim Binomial(N_{it}, p_{it})\]

We will allow initial abundance to vary by site \(\lambda_i=\lambda_{it}\), which seems reasonable as the data suggest significantly different populations at each site. The other parameters will be held constant across time and site so that \(\gamma=\gamma_{it}\), \(\omega=\omega_{it}\), and \(p=p_{it}\). This corresponds to the assumption that the recruitment rate, recovery rate and the detection probability are independent of site and constant over time. The link functions for this model are as follows:

Initial Site Abundance: \[ log(\lambda)=\beta_{0\lambda}+\beta_{1\lambda} * Site \]

Recruitment Rate: \[ log(\gamma)=\beta_{0\gamma} \]

Apparent Survival Rate: \[ logit(\omega)=\beta_{0\omega} \]

Detection Probability: \[ logit(p)=\beta_{0p} \]

# uframe for open population pcountOpen model
uframeO <- unmarkedFramePCO(y=Y, siteCovs=X, yearlySiteCovs=Z, numPrimary = 15)

openmodel <- pcountOpen(lambdaformula = ~Site, gammaformula = ~1, omegaformula = ~1, pformula = ~1, uframeO, K=200, se = TRUE)

# these models take forever to run, so it is good practice to save the model to disk for
# future use, that way you don't need to re-compute from scratch every time
saveRDS(openmodel, here::here("/data/bccdc/openmodel_2018_01_21.RDS"))
openmodel <- readRDS(here::here("/data/bccdc/openmodel_2018_01_21.RDS"))

summary(openmodel)

# get model parameter estimates
r_open <- ranef(openmodel)
a_open <- 100*bup(r_open)

rownames(a_open) <- c("North Shore", "Richmond", "Vancouver")
colnames(a_open) <- 2000:2014

# 95% confidence interval (likelihood based)
conf_open <- 100*confint(r_open, level=0.95) 

df_open <- data.frame(Abundance_Estimate=round(as.numeric(a_open),-2), Year=rep(2000:2014, times=rep(3, times=15)), Lower2.5=round(as.numeric(conf_open[,1,]),-2), Upper97.5=round(as.numeric(conf_open[,2,]),-2), Site=as.factor(rep(rownames(a_open), times=15)))

kable(round(a_open, -2), "html", caption = "Open Population Model Abundance Estimates") %>% kable_styling()

This model predicts a CDR of 0.986, which intuitively seems too large. Reasons for this could be our lack of covariates (both Site and Time based) for the detection probability. Essentially we have assumed that detection probability is constant over time and sites, which is probably not the case. If this article was more than an exploration of the data and the “unmarked” package, we would fit many more models and select those with best performance.

Of course, data are easier to understand when they are presented graphically rather than as tables. So let’s display the results using ggplot2.

# total abundance and 95% confidence intervals by year 
df_open_tot <- df_open %>% group_by(Year) %>% summarise(
  Nhat=sum(Abundance_Estimate), 
  Lower2.5=sum(Abundance_Estimate)-sum((Abundance_Estimate-Lower2.5)^2)^(0.5),
  Upper97.5=sum(Abundance_Estimate)+sum((Abundance_Estimate-Upper97.5)^2)^(0.5))

# ggplot to display total abundance estimate over time
ggp_open_tot <- ggplot(data=df_open_tot) +
  geom_errorbar(aes(x=Year, ymin=Lower2.5, ymax=Upper97.5), color="red", width=0.2) +
  geom_point(aes(y = Nhat, x=Year), color="black", size=1.25) +
  ggtitle("Open Pop Model (~Site ~1 ~1 ~1)\nAbundance Estimates with 95% CI") +
  geom_text(aes(x=Year, y=Nhat, label=Nhat), vjust=-1) +
  xlab(label = "Year") +
  ylab(label = "Total Estimated Depressed Individuals") +
  scale_x_continuous(breaks=seq(2000,2014,2)) +
  scale_y_continuous(breaks=seq(12000,18000,1000)) +
  theme_light()
ggp_open_tot

# ggplot to display site abundance estimates
ggp_open_sites <- ggplot(data=df_open) + 
  geom_errorbar(aes(x=Year, ymin=Lower2.5, ymax=Upper97.5), color="red", width=0.2) +
  geom_point(aes(y = Abundance_Estimate, x=Year), color="black", size=1.25) +
  ggtitle("Open Pop Model (~Site ~1 ~1 ~1)\nAbundance Estimates with 95% CI") +
  geom_text(aes(x=Year, y=Abundance_Estimate, label=Abundance_Estimate), vjust=-0.6) +
  xlab(label = "Year") +
  ylab(label = "Estimated Depressed Individuals") +
  scale_x_continuous(breaks=seq(2000,2014,2)) +
  facet_grid(Site~., scales = "free_y") +
  theme_light()
ggp_open_sites