unmarkedHow to estimate abundance without marking animals?
In an occupancy model:
Replace detection with counts and occupancy with abundance:
This is the N-mixture model of Royle (2004)
Parameters
$N_i$: Latent abundance at site $i$
$\lambda_i$: Expected abundance at site $i$
Math
$$N_i \sim \mathrm{Poisson}(\lambda_i)$$You could also use other distributions, such as:
Parameters/data
$y_{ij}$: Observed count at site $i$ for repeated sample $j$
$p_{ij}$: Probability of detecting an individual at site $i$ in sample $j$
Math
$$y_{ij} \sim \mathrm{Binomial}(p_{ij}, N_i)$$
library(unmarked)
data(mallard)
head(mallard.y)
| y.1 | y.2 | y.3 |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 0 | 0 |
| 3 | 2 | 1 |
| 0 | 0 | 0 |
| 3 | 0 | 3 |
| 0 | 0 | 0 |
jags_data <- list(y = mallard.y, M = 239, J = 3)
Math
$$N_i \sim \mathrm{Poisson}(\lambda)$$Code
for (i in 1:M){
N[i] ~ dpois(lambda)
}
lambda ~ dunif(0, 100) # prior should be positive
Math
$$y_{ij} \sim \mathrm{Binomial}(p,N_i)$$Code
for (i in 1:M){
for (j in 1:J){
y[i,j] ~ dbinom(p, N[i])
}
}
p ~ dunif(0,1) # prior must be between 0 and 1
cat("
model{
# State model
for (i in 1:M){
N[i] ~ dpois(lambda)
}
lambda ~ dunif(0, 100)
# Detection model
for (i in 1:M){
for (j in 1:J){
y[i,j] ~ dbinom(p, N[i])
}
}
p ~ dunif(0,1)
sumN <- sum(N[])
}
", file = "nmixture_null.bugs")
pars <- c("lambda", "p", "sumN")
# initialize z at maximum count ever observed at a site
Ninit <- apply(mallard.y, 1, max, na.rm=TRUE)
Ninit[is.infinite(Ninit)] <- NA
Ninit
inits <- function() list(N = Ninit)
Warning message in FUN(newX[, i], ...): “no non-missing arguments to max; returning -Inf” Warning message in FUN(newX[, i], ...): “no non-missing arguments to max; returning -Inf” Warning message in FUN(newX[, i], ...): “no non-missing arguments to max; returning -Inf” Warning message in FUN(newX[, i], ...): “no non-missing arguments to max; returning -Inf”
library(jagsUI)
jags_null <- jags(data = jags_data,
inits = inits,
parameters.to.save = pars,
model.file = "nmixture_null.bugs",
n.chains = 3, n.iter = 10000, n.burnin = 8000,
parallel=TRUE)
Processing function input....... Done. Beginning parallel processing using 3 cores. Console output will be suppressed. Parallel processing completed. Calculating statistics....... Done.
jagsUI::traceplot(jags_null)
jags_null
JAGS output for model 'nmixture_null.bugs', generated by jagsUI.
Estimates based on 3 chains of 10000 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 8000 iterations and thin rate = 1,
yielding 6000 total samples from the joint posterior.
MCMC ran in parallel for 0.129 minutes at time 2023-07-19 12:21:33.799883.
mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
lambda 0.352 0.041 0.276 0.350 0.437 FALSE 1 1.000 6000
p 0.645 0.038 0.567 0.646 0.719 FALSE 1 1.001 2621
sumN 83.176 3.641 77.000 83.000 91.000 FALSE 1 1.002 1170
deviance 261.189 14.429 237.703 259.782 292.910 FALSE 1 1.003 1196
Successful convergence based on Rhat values (all < 1.1).
Rhat is the potential scale reduction factor (at convergence, Rhat=1).
For each parameter, n.eff is a crude measure of effective sample size.
overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.
DIC info: (pD = var(deviance)/2)
pD = 104 and DIC = 365.142
DIC is an estimate of expected predictive error (lower is better).
Including the site-level covariates
head(mallard.site)
# already standardized?
round(colMeans(mallard.site), 2)
round(apply(mallard.site, 2, sd), 2)
| elev | length | forest | |
|---|---|---|---|
| <dbl> | <dbl> | <dbl> | |
| 1 | -1.173 | 0.801 | -1.156 |
| 2 | -1.127 | 0.115 | -0.501 |
| 3 | -0.198 | -0.479 | -0.101 |
| 4 | -0.105 | 0.315 | 0.008 |
| 5 | -1.034 | -1.102 | -1.193 |
| 6 | -0.848 | 0.741 | 0.917 |
jags_data <- list(y = mallard.y, M = 239, J = 3,
elev = mallard.site$elev, forest = mallard.site$forest)
str(jags_data)
List of 5 $ y : num [1:239, 1:3] 0 0 3 0 3 0 0 0 0 0 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:3] "y.1" "y.2" "y.3" $ M : num 239 $ J : num 3 $ elev : num [1:239] -1.173 -1.127 -0.198 -0.105 -1.034 ... $ forest: num [1:239] -1.156 -0.501 -0.101 0.008 -1.193 ...
Math
$$N_i \sim \mathrm{Poisson}(\lambda_i)$$$$\mathrm{log}(\lambda_i) = \beta_0 + \beta_{ele} \cdot \mathrm{elev}_i + \beta_{for} \cdot \mathrm{forest}_i$$Code
for (i in 1:M){
N[i] ~ dpois(lambda[i])
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i]
}
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
cat("
model{
# State model
for (i in 1:M){
N[i] ~ dpois(lambda[i])
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i]
}
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
# Detection model
for (i in 1:M){
for (j in 1:J){
y[i,j] ~ dbinom(p, N[i])
}
}
p ~ dunif(0,1)
sumN <- sum(N[])
}
", file = "nmixture_covs.bugs")
pars <- c("beta_0", "beta_ele", "beta_for", "p", "sumN")
inits <- function() list(N = Ninit)
library(jagsUI)
jags_covs <- jags(data = jags_data,
inits = inits,
parameters.to.save = pars,
model.file = "nmixture_covs.bugs",
n.chains = 3, n.iter = 20000, n.burnin = 18000,
parallel=TRUE)
Processing function input....... Done. Beginning parallel processing using 3 cores. Console output will be suppressed. Parallel processing completed. Calculating statistics....... Done.
jagsUI::traceplot(jags_covs)
jags_covs
JAGS output for model 'nmixture_covs.bugs', generated by jagsUI.
Estimates based on 3 chains of 20000 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 18000 iterations and thin rate = 1,
yielding 6000 total samples from the joint posterior.
MCMC ran in parallel for 0.384 minutes at time 2023-07-19 12:21:42.594846.
mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
beta_0 -2.322 0.278 -2.943 -2.299 -1.829 FALSE 1 1.236 13
beta_ele -1.615 0.255 -2.187 -1.608 -1.146 FALSE 1 1.357 10
beta_for -0.951 0.155 -1.246 -0.954 -0.657 FALSE 1 1.008 298
p 0.613 0.045 0.521 0.616 0.696 FALSE 1 1.001 1503
sumN 88.908 5.396 80.000 88.000 102.000 FALSE 1 1.001 1482
deviance 259.127 14.175 236.036 257.557 290.180 FALSE 1 1.000 6000
**WARNING** Rhat values indicate convergence failure.
Rhat is the potential scale reduction factor (at convergence, Rhat=1).
For each parameter, n.eff is a crude measure of effective sample size.
overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.
DIC info: (pD = var(deviance)/2)
pD = 100.5 and DIC = 359.609
DIC is an estimate of expected predictive error (lower is better).
In each MCMC iteration:
fit (Pearson's $\chi^2$) on the real datafit_new for the new datasetAfter the run:
fit < fit_new = Bayesian p-valueMath
$$\hat{y}_{ij} = N_i \cdot p $$$$\mathrm{fit}_{\chi^2} = \sum{\frac{(y_{ij} - \hat{y}_{ij})^2}{\hat{y}_{ij}}}$$Code
# real data
for (i in 1:M){
for (j in 1:J){
yhat[i,j] = N[i] * p
chi2[i,j] = (y[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit = sum(chi2)
Simulate new dataset
$$y^{new}_{ij} \sim \mathrm{Binomial}(p, N_i)$$Calculate fit for new dataset
$$\hat{y}_{ij} = N_i \cdot p $$$$\mathrm{fit}^{new}_{\chi^2} = \sum{\frac{(y^{new}_{ij} - \hat{y}_{ij})^2}{\hat{y}_{ij}}}$$Code
# simulated data
for (i in 1:M){
for (j in 1:J){
y_new[i,j] ~ dbinom(p, N[i]) # simulate new datapoint
yhat[i,j] = N[i] * p
chi2_new[i,j] = (y_new[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit_new = sum(chi2_new)
cat("
model{
# State model
for (i in 1:M){
N[i] ~ dpois(lambda[i])
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i]
}
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
# Detection model
for (i in 1:M){
for (j in 1:J){
y[i,j] ~ dbinom(p, N[i])
}
}
p ~ dbeta(1,1)
# Fit statistic for real data
for (i in 1:M){
for (j in 1:J){
yhat[i,j] = N[i] * p + 0.001 # add small value to avoid divide by zero
chi2[i,j] = (y[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit = sum(chi2)
# Fit statistic for simulated data
for (i in 1:M){
for (j in 1:J){
y_new[i,j] ~ dbinom(p, N[i]) # simulate new datapoint
chi2_new[i,j] = (y_new[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit_new = sum(chi2_new)
sumN <- sum(N[])
}
", file = "nmixture_covs.bugs")
pars <- c("beta_0", "beta_ele", "beta_for", "p", "fit", "fit_new", "sumN")
inits <- function() list(N = Ninit)
library(jagsUI)
jags_covs <- jags(data = jags_data,
inits = inits,
parameters.to.save = pars,
model.file = "nmixture_covs.bugs",
n.chains = 3, n.iter = 10000, n.burnin = 8000,
parallel=TRUE)
Processing function input....... Done. Beginning parallel processing using 3 cores. Console output will be suppressed. Parallel processing completed. Calculating statistics....... Done.
Bayesian p-value
mean(jags_covs$sims.list$fit < jags_covs$sims.list$fit_new)
Posterior predictive check plot
Plot fit vs. fit_new.
A line with intercept 0 and slope 1 should approximately bisect the point cloud.
plot(jags_covs$sims.list$fit, jags_covs$sims.list$fit_new)
abline(a=0, b=1)
pp.check(jags_covs, "fit","fit_new")
Based on this test, the model appears to fit the data poorly.
Potential solutions:
More research should be done in this area to develop better GoF tests.
Use JAGS to fit the following model.
mallard.obs)Note there are missing values in the observation covariate ivel! See the dynamic occupancy model code for one way to handle this in JAGS. Hint: you can model the missing values of ivel.
Don't forget: p now becomes indexed by i and j
Check goodness-of-fit.
head(mallard.obs$ivel)
dim(mallard.obs$ivel)
| ivel.1 | ivel.2 | ivel.3 |
|---|---|---|
| -0.506 | -0.506 | -0.506 |
| -0.934 | -0.991 | -1.162 |
| -1.136 | -1.339 | -1.610 |
| -0.819 | -0.927 | -1.197 |
| 0.638 | 0.880 | 1.042 |
| -1.329 | -1.042 | -0.899 |
jags_data <- list(y = mallard.y, M = 239, J = 3,
elev = mallard.site$elev, forest = mallard.site$forest,
ivel = mallard.obs$ivel)
str(jags_data)
List of 6 $ y : num [1:239, 1:3] 0 0 3 0 3 0 0 0 0 0 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:3] "y.1" "y.2" "y.3" $ M : num 239 $ J : num 3 $ elev : num [1:239] -1.173 -1.127 -0.198 -0.105 -1.034 ... $ forest: num [1:239] -1.156 -0.501 -0.101 0.008 -1.193 ... $ ivel : num [1:239, 1:3] -0.506 -0.934 -1.136 -0.819 0.638 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr [1:3] "ivel.1" "ivel.2" "ivel.3"
cat("
model{
# State model
for (i in 1:M){
N[i] ~ dpois(lambda[i])
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i] +
beta_ele2 * elev[i] * elev[i]
}
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
beta_ele2 ~ dnorm(0, 0.01)
# Model missing values in ivel
for (i in 1:M){
for (j in 1:J){
ivel[i,j] ~ dnorm(ivel_mean, ivel_tau)
}
}
ivel_mean ~ dnorm(0, 0.01)
ivel_tau <- pow(ivel_sd, -2)
ivel_sd ~ dunif(0, 10)
# Detection model
for (i in 1:M){
for (j in 1:J){
logit(p[i,j]) <- alpha_0 + alpha_ivel * ivel[i,j]
y[i,j] ~ dbinom(p[i,j], N[i])
}
}
alpha_0 ~ dlogis(0, 1)
alpha_ivel ~ dlogis(0, 1)
# Fit statistic for real data
for (i in 1:M){
for (j in 1:J){
yhat[i,j] <- N[i] * p[i,j] + 0.001 # add small value to avoid divide by zero
chi2[i,j] <- (y[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit <- sum(chi2)
# Fit statistic for simulated data
for (i in 1:M){
for (j in 1:J){
y_new[i,j] ~ dbinom(p[i,j], N[i]) # simulate new datapoint
chi2_new[i,j] <- (y_new[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit_new <- sum(chi2_new)
sumN <- sum(N[])
}
", file = "nmixture_obscovs.bugs")
pars <- c("beta_0", "beta_ele", "beta_for", "beta_ele2", "alpha_0", "alpha_ivel", "fit", "fit_new", "sumN")
inits <- function() list(N = Ninit)
library(jagsUI)
jags_covs2 <- jags(data = jags_data,
inits = inits,
parameters.to.save = pars,
model.file = "nmixture_obscovs.bugs",
n.chains = 3, n.iter = 10000, n.burnin = 8000,
parallel=TRUE)
Processing function input....... Done. Beginning parallel processing using 3 cores. Console output will be suppressed. Parallel processing completed. Calculating statistics....... Done.
jags_covs2
jagsUI::traceplot(jags_covs2)
JAGS output for model 'nmixture_obscovs.bugs', generated by jagsUI.
Estimates based on 3 chains of 10000 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 8000 iterations and thin rate = 1,
yielding 6000 total samples from the joint posterior.
MCMC ran in parallel for 0.715 minutes at time 2023-07-19 12:22:21.342982.
mean sd 2.5% 50% 97.5% overlap0 f Rhat
beta_0 -1.916 0.292 -2.500 -1.913 -1.359 FALSE 1.000 1.017
beta_ele -1.438 0.245 -1.983 -1.420 -1.014 FALSE 1.000 1.005
beta_for -0.761 0.180 -1.131 -0.753 -0.424 FALSE 1.000 1.001
beta_ele2 -0.056 0.263 -0.595 -0.036 0.438 TRUE 0.555 1.013
alpha_0 0.402 0.219 -0.040 0.405 0.818 TRUE 0.963 1.003
alpha_ivel 0.439 0.171 0.109 0.442 0.770 FALSE 0.996 1.002
fit 78.023 12.696 59.058 76.216 109.082 FALSE 1.000 1.004
fit_new 65.582 13.616 43.622 63.966 97.240 FALSE 1.000 1.005
sumN 95.117 8.808 82.000 94.000 116.000 FALSE 1.000 1.008
deviance 2159.504 18.342 2128.796 2158.085 2200.556 FALSE 1.000 1.006
n.eff
beta_0 132
beta_ele 6000
beta_for 2215
beta_ele2 157
alpha_0 672
alpha_ivel 1117
fit 906
fit_new 800
sumN 401
deviance 491
Successful convergence based on Rhat values (all < 1.1).
Rhat is the potential scale reduction factor (at convergence, Rhat=1).
For each parameter, n.eff is a crude measure of effective sample size.
overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.
DIC info: (pD = var(deviance)/2)
pD = 167.6 and DIC = 2327.08
DIC is an estimate of expected predictive error (lower is better).
pp.check(jags_covs2, 'fit', 'fit_new')
Math
$$N_i \sim \mathrm{NegativeBinomial}(\lambda_i, \phi)$$$$\mathrm{log}(\lambda_i) = \beta_0 + \beta_{ele} \cdot \mathrm{elev}_i + \beta_{for} \cdot \mathrm{forest}_i$$$\phi$ is the dispersion parameter, as $\phi$ gets small NB converges with Poisson
Code
JAGS parameterizes the negative binomial in the success/failure way:
dnegbin(p, r)
We'll call the p parameter s (for success) since we will use p for detection.
N[i] ~ dnegbin(s[i], r)
We need to convert between our parameterization ($\lambda$) and the JAGS way ($s$)
$$s_i = \frac{r}{r+\lambda_i}$$$$\phi =\frac{1}{r}$$for (i in 1:M){
N[i] ~ dnegbin(s[i], r)
s[i] <- r / (r + lambda[i])
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i] +
beta_ele2 * elev[i] * elev[i]
}
r ~ dunif(0, 50)
phi <- 1/r
# other priors
cat("
model{
# State model
for (i in 1:M){
N[i] ~ dnegbin(s[i], r)
s[i] <- r / (r + lambda[i])
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i] +
beta_ele2 * elev[i] * elev[i]
}
r ~ dunif(0, 50)
phi <- 1/r
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
beta_ele2 ~ dnorm(0, 0.01)
# Model missing values in ivel
for (i in 1:M){
for (j in 1:J){
ivel[i,j] ~ dnorm(ivel_mean, ivel_tau)
}
}
ivel_mean ~ dnorm(0, 0.01)
ivel_tau <- pow(ivel_sd, -2)
ivel_sd ~ dunif(0, 10)
# Detection model
for (i in 1:M){
for (j in 1:J){
logit(p[i,j]) <- alpha_0 + alpha_ivel * ivel[i,j]
y[i,j] ~ dbinom(p[i,j], N[i])
}
}
alpha_0 ~ dunif(-10,10)
alpha_ivel ~ dnorm(0, 0.01)
# Fit statistic for real data
for (i in 1:M){
for (j in 1:J){
yhat[i,j] <- N[i] * p[i,j] + 0.001 # add small value to avoid divide by zero
chi2[i,j] <- (y[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit <- sum(chi2)
# Fit statistic for simulated data
for (i in 1:M){
for (j in 1:J){
y_new[i,j] ~ dbinom(p[i,j], N[i]) # simulate new datapoint
chi2_new[i,j] <- (y_new[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit_new <- sum(chi2_new)
sumN <- sum(N[])
}
", file = "nmixture_nb.bugs")
pars <- c("beta_0", "beta_ele", "beta_for", "beta_ele2", "alpha_0", "alpha_ivel", "phi", "fit", "fit_new", "sumN")
inits <- function() list(N = Ninit)
library(jagsUI)
jags_nb <- jags(data = jags_data,
inits = inits,
parameters.to.save = pars,
model.file = "nmixture_nb.bugs",
n.chains = 3, n.iter = 10000, n.burnin = 8000,
parallel=TRUE)
Processing function input....... Done. Beginning parallel processing using 3 cores. Console output will be suppressed. Parallel processing completed. Calculating statistics....... Done.
jagsUI::traceplot(jags_nb)
jags_nb
JAGS output for model 'nmixture_nb.bugs', generated by jagsUI.
Estimates based on 3 chains of 10000 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 8000 iterations and thin rate = 1,
yielding 6000 total samples from the joint posterior.
MCMC ran in parallel for 1.011 minutes at time 2023-07-19 12:23:05.539394.
mean sd 2.5% 50% 97.5% overlap0 f Rhat
beta_0 -1.803 0.396 -2.598 -1.795 -1.034 FALSE 1.000 1.007
beta_ele -1.506 0.290 -2.130 -1.488 -0.987 FALSE 1.000 1.005
beta_for -0.714 0.230 -1.171 -0.710 -0.279 FALSE 0.999 1.002
beta_ele2 0.031 0.335 -0.661 0.045 0.655 TRUE 0.553 1.008
alpha_0 -0.092 0.339 -0.824 -0.060 0.497 TRUE 0.578 1.011
alpha_ivel 0.238 0.205 -0.158 0.235 0.645 TRUE 0.881 1.002
phi 2.443 0.830 1.144 2.330 4.368 FALSE 1.000 1.000
fit 81.230 14.435 59.833 78.802 115.121 FALSE 1.000 1.014
fit_new 76.655 17.303 50.357 74.109 117.080 FALSE 1.000 1.012
sumN 119.791 23.745 89.975 114.000 179.025 FALSE 1.000 1.027
deviance 2156.426 16.525 2127.885 2155.039 2192.254 FALSE 1.000 1.009
n.eff
beta_0 295
beta_ele 534
beta_for 1503
beta_ele2 364
alpha_0 228
alpha_ivel 1724
phi 5753
fit 232
fit_new 236
sumN 163
deviance 274
Successful convergence based on Rhat values (all < 1.1).
Rhat is the potential scale reduction factor (at convergence, Rhat=1).
For each parameter, n.eff is a crude measure of effective sample size.
overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.
DIC info: (pD = var(deviance)/2)
pD = 135.6 and DIC = 2292.008
DIC is an estimate of expected predictive error (lower is better).
pp.check(jags_nb, "fit","fit_new")
Useful when there are more 0s in the dataset than you would expect from just the Poisson
Math
$$N_i \sim \mathrm{Poisson}(w_i \cdot \lambda_i)$$$$w_i \sim \mathrm{Bernoulli}(1-\psi)$$$$\mathrm{log}(\lambda_i) = \beta_0 + \beta_{ele} \cdot \mathrm{elev}_i + \beta_{for} \cdot \mathrm{forest}_i$$$\psi$ is the zero-inflation parameter (~ % of samples that are forced 0s)
In ecology we could think of this distribution like so:
Code
for (i in 1:M){
N[i] ~ dpois(w[i] * lambda[i])
w[i] ~ dbern(1-psi)
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i] +
beta_ele2 * elev[i] * elev[i]
}
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
beta_ele2 ~ dnorm(0, 0.01
psi ~ dunif(0, 1)
cat("
model{
# State model
for (i in 1:M){
N[i] ~ dpois(w[i]*lambda[i])
w[i] ~ dbern(1-psi)
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i] +
beta_ele2 * elev[i] * elev[i]
}
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
beta_ele2 ~ dnorm(0, 0.01)
psi ~ dunif(0, 1)
# Model missing values in ivel
for (i in 1:M){
for (j in 1:J){
ivel[i,j] ~ dnorm(ivel_mean, ivel_tau)
}
}
ivel_mean ~ dnorm(0, 0.01)
ivel_tau <- pow(ivel_sd, -2)
ivel_sd ~ dunif(0, 10)
# Detection model
for (i in 1:M){
for (j in 1:J){
logit(p[i,j]) <- alpha_0 + alpha_ivel * ivel[i,j]
y[i,j] ~ dbinom(p[i,j], N[i])
}
}
alpha_0 ~ dunif(-10,10)
alpha_ivel ~ dnorm(0, 0.01)
# Fit statistic for real data
for (i in 1:M){
for (j in 1:J){
yhat[i,j] <- N[i] * p[i,j] + 0.001 # add small value to avoid divide by zero
chi2[i,j] <- (y[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit <- sum(chi2)
# Fit statistic for simulated data
for (i in 1:M){
for (j in 1:J){
y_new[i,j] ~ dbinom(p[i,j], N[i]) # simulate new datapoint
chi2_new[i,j] <- (y_new[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit_new <- sum(chi2_new)
sumN <- sum(N[])
}
", file = "nmixture_zip.bugs")
pars <- c("beta_0", "beta_ele", "beta_for", "beta_ele2", "alpha_0", "alpha_ivel", "psi", "fit", "fit_new", "sumN")
inits <- function() list(N = Ninit, psi = 0)
library(jagsUI)
jags_zip <- jags(data = jags_data,
inits = inits,
parameters.to.save = pars,
model.file = "nmixture_zip.bugs",
n.chains = 3, n.iter = 10000, n.burnin = 8000,
parallel=TRUE)
Processing function input....... Done. Beginning parallel processing using 3 cores. Console output will be suppressed. Parallel processing completed. Calculating statistics....... Done.
jagsUI::traceplot(jags_zip)
jags_zip
JAGS output for model 'nmixture_zip.bugs', generated by jagsUI.
Estimates based on 3 chains of 10000 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 8000 iterations and thin rate = 1,
yielding 6000 total samples from the joint posterior.
MCMC ran in parallel for 0.792 minutes at time 2023-07-19 12:24:07.475039.
mean sd 2.5% 50% 97.5% overlap0 f Rhat
beta_0 -0.949 0.406 -1.708 -0.954 -0.138 FALSE 0.991 1.031
beta_ele -1.408 0.264 -1.982 -1.393 -0.940 FALSE 1.000 1.003
beta_for -0.578 0.201 -0.980 -0.580 -0.185 FALSE 0.999 1.022
beta_ele2 -0.093 0.273 -0.630 -0.089 0.436 TRUE 0.637 1.012
alpha_0 0.074 0.312 -0.614 0.104 0.619 TRUE 0.637 1.023
alpha_ivel 0.466 0.171 0.135 0.466 0.797 FALSE 0.997 1.002
psi 0.516 0.084 0.342 0.520 0.670 FALSE 1.000 1.037
fit 80.059 14.109 59.963 77.609 114.254 FALSE 1.000 1.011
fit_new 72.203 15.328 48.180 69.978 108.101 FALSE 1.000 1.014
sumN 112.540 20.396 87.000 108.000 166.025 FALSE 1.000 1.025
deviance 2162.634 17.492 2133.511 2160.811 2201.165 FALSE 1.000 1.014
n.eff
beta_0 75
beta_ele 1099
beta_for 102
beta_ele2 207
alpha_0 98
alpha_ivel 1036
psi 60
fit 220
fit_new 160
sumN 111
deviance 164
Successful convergence based on Rhat values (all < 1.1).
Rhat is the potential scale reduction factor (at convergence, Rhat=1).
For each parameter, n.eff is a crude measure of effective sample size.
overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.
DIC info: (pD = var(deviance)/2)
pD = 151.2 and DIC = 2313.801
DIC is an estimate of expected predictive error (lower is better).
pp.check(jags_zip, "fit","fit_new")
(also called Poisson-lognormal model)
Math
$$N_i \sim \mathrm{Poisson}(\lambda_i)$$$$\mathrm{log}(\lambda_i) = \beta_0 + \beta_{ele} \cdot \mathrm{elev}_i + \beta_{for} \cdot \mathrm{forest}_i + \epsilon_i$$$$\epsilon_i \sim \mathrm{Normal}(0, \sigma^2)$$Code
for (i in 1:M){
N[i] ~ dpois(lambda[i])
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i] +
beta_ele2 * elev[i] * elev[i] + eps[i]
eps[i] ~ dnorm(0, tau)
}
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
beta_ele2 ~ dnorm(0, 0.01)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 10)
cat("
model{
# State model
for (i in 1:M){
N[i] ~ dpois(lambda[i])
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i] +
beta_ele2 * elev[i] * elev[i] + eps[i]
eps[i] ~ dnorm(0, tau)
}
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
beta_ele2 ~ dnorm(0, 0.01)
tau <- pow(sigma, -2)
sigma ~ dunif(0, 10)
# Model missing values in ivel
for (i in 1:M){
for (j in 1:J){
ivel[i,j] ~ dnorm(ivel_mean, ivel_tau)
}
}
ivel_mean ~ dnorm(0, 0.01)
ivel_tau <- pow(ivel_sd, -2)
ivel_sd ~ dunif(0, 10)
# Detection model
for (i in 1:M){
for (j in 1:J){
logit(p[i,j]) <- alpha_0 + alpha_ivel * ivel[i,j]
y[i,j] ~ dbinom(p[i,j], N[i])
}
}
alpha_0 ~ dunif(-10,10)
alpha_ivel ~ dnorm(0, 0.01)
# Fit statistic for real data
for (i in 1:M){
for (j in 1:J){
yhat[i,j] <- N[i] * p[i,j] + 0.001 # add small value to avoid divide by zero
chi2[i,j] <- (y[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit <- sum(chi2)
# Fit statistic for simulated data
for (i in 1:M){
for (j in 1:J){
y_new[i,j] ~ dbinom(p[i,j], N[i]) # simulate new datapoint
chi2_new[i,j] <- (y_new[i,j] - yhat[i,j])^2 / yhat[i,j]
}
}
fit_new <- sum(chi2_new)
sumN <- sum(N[])
}
", file = "nmixture_pln.bugs")
pars <- c("beta_0", "beta_ele", "beta_for", "beta_ele2", "alpha_0", "alpha_ivel", "sigma", "fit", "fit_new", "sumN")
inits <- function() list(N = Ninit, psi = 0)
library(jagsUI)
jags_pln <- jags(data = jags_data,
inits = inits,
parameters.to.save = pars,
model.file = "nmixture_pln.bugs",
n.chains = 3, n.iter = 10000, n.burnin = 8000,
parallel=TRUE)
Processing function input....... Done. Beginning parallel processing using 3 cores. Console output will be suppressed. Parallel processing completed. Calculating statistics....... Done.
jagsUI::traceplot(jags_pln)
jags_pln
JAGS output for model 'nmixture_pln.bugs', generated by jagsUI.
Estimates based on 3 chains of 10000 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 8000 iterations and thin rate = 1,
yielding 6000 total samples from the joint posterior.
MCMC ran in parallel for 0.811 minutes at time 2023-07-19 12:24:56.146931.
mean sd 2.5% 50% 97.5% overlap0 f Rhat
beta_0 -3.140 0.551 -4.292 -3.122 -2.191 FALSE 1.000 1.023
beta_ele -1.623 0.309 -2.252 -1.613 -1.048 FALSE 1.000 1.028
beta_for -0.755 0.281 -1.352 -0.736 -0.247 FALSE 0.998 1.007
beta_ele2 0.197 0.363 -0.489 0.185 0.889 TRUE 0.697 1.013
alpha_0 -0.134 0.349 -0.960 -0.107 0.460 TRUE 0.630 1.108
alpha_ivel 0.173 0.230 -0.265 0.170 0.634 TRUE 0.771 1.007
sigma 1.558 0.285 1.107 1.522 2.269 FALSE 1.000 1.046
fit 81.587 15.266 60.488 78.693 121.967 FALSE 1.000 1.113
fit_new 77.232 17.521 51.540 74.137 120.784 FALSE 1.000 1.107
sumN 124.769 32.469 91.000 117.000 205.000 FALSE 1.000 1.136
deviance 2154.674 16.230 2128.274 2152.758 2191.424 FALSE 1.000 1.077
n.eff
beta_0 370
beta_ele 88
beta_for 277
beta_ele2 180
alpha_0 33
alpha_ivel 329
sigma 56
fit 32
fit_new 32
sumN 33
deviance 36
**WARNING** Rhat values indicate convergence failure.
Rhat is the potential scale reduction factor (at convergence, Rhat=1).
For each parameter, n.eff is a crude measure of effective sample size.
overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.
DIC info: (pD = var(deviance)/2)
pD = 124.3 and DIC = 2278.968
DIC is an estimate of expected predictive error (lower is better).
pp.check(jags_pln, "fit","fit_new")
stats <- t(sapply(list(jags_covs2, jags_nb, jags_zip, jags_pln),
function(x) c(mean=x$mean$beta_ele, lower=x$q2.5$beta_ele, upper=x$q97.5$beta_ele)))
plot(1:4, stats[,1], pch=19, xaxt='n', xlim=c(0.5, 4.5),
ylim=c(-2.5, -0.5), xlab="", ylab="beta_ele")
axis(1, at=1:4, labels=c("Pois","NB","ZIP","PLN"))
segments(1:4, stats[,2], 1:4, stats[,3])
stats <- t(sapply(list(jags_covs2, jags_nb, jags_zip, jags_pln),
function(x) c(mean=x$mean$sumN, lower=x$q2.5$sumN, upper=x$q97.5$sumN)))
plot(1:4, stats[,1], pch=19, xaxt='n', xlim=c(0.5, 4.5),
ylim=c(50, 210), xlab="", ylab="sumN")
axis(1, at=1:4, labels=c("Pois","NB","ZIP","PLN"))
segments(1:4, stats[,2], 1:4, stats[,3])
unmarked¶Look at the data again:
library(unmarked)
data(mallard)
head(mallard.y)
head(mallard.site)
| y.1 | y.2 | y.3 |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 0 | 0 |
| 3 | 2 | 1 |
| 0 | 0 | 0 |
| 3 | 0 | 3 |
| 0 | 0 | 0 |
| elev | length | forest | |
|---|---|---|---|
| <dbl> | <dbl> | <dbl> | |
| 1 | -1.173 | 0.801 | -1.156 |
| 2 | -1.127 | 0.115 | -0.501 |
| 3 | -0.198 | -0.479 | -0.101 |
| 4 | -0.105 | 0.315 | 0.008 |
| 5 | -1.034 | -1.102 | -1.193 |
| 6 | -0.848 | 0.741 | 0.917 |
umf <- unmarkedFramePCount(y = mallard.y, siteCovs = mallard.site)
head(umf)
Data frame representation of unmarkedFrame object. y.1 y.2 y.3 elev length forest 1 0 0 0 -1.173 0.801 -1.156 2 0 0 0 -1.127 0.115 -0.501 3 3 2 1 -0.198 -0.479 -0.101 4 0 0 0 -0.105 0.315 0.008 5 3 0 3 -1.034 -1.102 -1.193 6 0 0 0 -0.848 0.741 0.917 7 0 0 0 -0.910 0.115 -1.083 8 0 0 0 -1.003 -1.007 -0.792 9 0 0 0 -0.058 -0.913 0.553 10 0 0 0 -0.631 1.556 0.808
unmarked function pcountpcount fits the N-mixture model: does not need to be point-count data!K, an upper bound on possible site abundancenmix_null <- pcount(~1~1, umf, K = 20)
summary(nmix_null)
Warning message: “4 sites have been discarded because of missing data.”
Call:
pcount(formula = ~1 ~ 1, data = umf, K = 20)
Abundance (log-scale):
Estimate SE z P(>|z|)
-1.06 0.118 -9 2.17e-19
Detection (logit-scale):
Estimate SE z P(>|z|)
0.611 0.17 3.59 0.00033
AIC: 631.8909
Number of sites: 235
ID of sites removed due to NA: 12 69 118 146
optim convergence code: 0
optim iterations: 15
Bootstrap iterations: 0
Estimates are on the transformed scale! To get them on the original scale, use predict.
predict(nmix_null, type='state')[1,]
Warning message: “4 sites have been discarded because of missing data.”
| Predicted | SE | lower | upper | |
|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 0.3460371 | 0.04078127 | 0.2746671 | 0.4359521 |
predict(nmix_null, type='det')[1,]
Warning message: “4 sites have been discarded because of missing data.”
| Predicted | SE | lower | upper | |
|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 0.6482038 | 0.03881639 | 0.5689397 | 0.7200641 |
max(umf@y, na.rm=TRUE)
test_K <- function(K){
fit <- suppressWarnings(pcount(~1~1, umf, K=K))
exp(coef(fit)[1])
}
kvals <- c(13,14,15,20,25,30)
lams <- sapply(kvals, test_K)
plot(kvals, lams, type='o', pch=19, xlim=c(10, 30),
xlab="K", ylab="lambda estimate")
abline(v=12, col='red', lty=2)
Just need to change the formulas
nmix_covs <- pcount(~1~elev+forest, umf, K = 20)
nmix_covs
Warning message: “4 sites have been discarded because of missing data.”
Call:
pcount(formula = ~1 ~ elev + forest, data = umf, K = 20)
Abundance:
Estimate SE z P(>|z|)
(Intercept) -1.905 0.228 -8.34 7.58e-17
elev -1.317 0.222 -5.93 2.95e-09
forest -0.738 0.163 -4.53 5.79e-06
Detection:
Estimate SE z P(>|z|)
0.492 0.188 2.62 0.00878
AIC: 530.7381
The estimates are a little different from JAGS, but estimates in both cases are quite uncertain
data.frame(JAGS=unlist(jags_covs$mean[1:4]),
unmarked=coef(nmix_covs))
unm <- c(coef(nmix_covs)[1:3], plogis(coef(nmix_covs)[4]))
data.frame(JAGS=unlist(jags_covs$mean[1:4]),
unmarked=unm)
| JAGS | unmarked | |
|---|---|---|
| <dbl> | <dbl> | |
| beta_0 | -2.3535487 | -1.9051037 |
| beta_ele | -1.6609773 | -1.3167989 |
| beta_for | -0.9428308 | -0.7378125 |
| p | 0.6108338 | 0.4921321 |
| JAGS | unmarked | |
|---|---|---|
| <dbl> | <dbl> | |
| beta_0 | -2.3535487 | -1.9051037 |
| beta_ele | -1.6609773 | -1.3167989 |
| beta_for | -0.9428308 | -0.7378125 |
| p | 0.6108338 | 0.6206086 |
newdata, holding other covariates at mean or medianpredict to generate corresponding abundance values# sequence of covariate values
elev_rng <- range(mallard.site$elev)
elev_seq <- seq(elev_rng[1], elev_rng[2], length.out=100)
# insert into newdata
nd <- data.frame(elev=elev_seq, forest=median(mallard.site$forest))
# predict
pr <- predict(nmix_covs, type='state', newdata=nd)
head(pr)
| Predicted | SE | lower | upper | |
|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 1.0343542 | 0.2120818 | 0.6920552 | 1.545959 |
| 2 | 0.9824581 | 0.1955079 | 0.6651598 | 1.451116 |
| 3 | 0.9331658 | 0.1802507 | 0.6390581 | 1.362628 |
| 4 | 0.8863465 | 0.1662269 | 0.6137164 | 1.280087 |
| 5 | 0.8418763 | 0.1533586 | 0.5891027 | 1.203111 |
| 6 | 0.7996373 | 0.1415724 | 0.5651860 | 1.131344 |
plot(elev_seq, pr$Predicted, type='l')
# add CIs
plot(elev_seq, pr$Predicted, type='l', ylim=c(0, 1.5))
polygon(c(elev_seq, rev(elev_seq)), c(pr$lower, rev(pr$upper)), col='lightgrey', border=NA)
lines(elev_seq, pr$Predicted)
plotEffects(nmix_covs, type = "state", covariate = "elev")
Actually, we will use tools in AICcmodavg package
Process:
Calculate manually in R:
yhat <- fitted(nmix_covs)
sum((umf@y - yhat)^2 / yhat, na.rm=TRUE)
library(AICcmodavg)
Nmix.gof.test(nmix_covs, nsim = 30)
Chi-square goodness-of-fit for N-mixture model of 'unmarkedFitPCount' class Observed chi-square statistic = 1405.956 Number of bootstrap samples = 30 P-value = 0 Quantiles of bootstrapped statistics: 0% 25% 50% 75% 100% 459 558 625 709 987 Estimate of c-hat = 2.17
nmix_covs_nb <- pcount(~1~elev+forest, umf, K = 20, mixture = "NB")
summary(nmix_covs_nb)
Warning message: “4 sites have been discarded because of missing data.”
Call:
pcount(formula = ~1 ~ elev + forest, data = umf, K = 20, mixture = "NB")
Abundance (log-scale):
Estimate SE z P(>|z|)
(Intercept) -1.753 0.269 -6.53 6.67e-11
elev -1.378 0.266 -5.19 2.14e-07
forest -0.725 0.215 -3.37 7.40e-04
Detection (logit-scale):
Estimate SE z P(>|z|)
0.00988 0.243 0.0407 0.968
Dispersion (log-scale):
Estimate SE z P(>|z|)
-0.885 0.317 -2.79 0.00525
AIC: 481.0938
Number of sites: 235
ID of sites removed due to NA: 12 69 118 146
optim convergence code: 0
optim iterations: 41
Bootstrap iterations: 0
Here the dispersion estimate is $\mathrm{log}(r)$, not $\phi$
Nmix.gof.test(nmix_covs_nb, nsim = 30)
Chi-square goodness-of-fit for N-mixture model of 'unmarkedFitPCount' class Observed chi-square statistic = 1475.87 Number of bootstrap samples = 30 P-value = 0.1 Quantiles of bootstrapped statistics: 0% 25% 50% 75% 100% 791 898 982 1117 2413 Estimate of c-hat = 1.34
nmix_covs_zip <- pcount(~1~elev+forest, umf, K = 20, mixture = "ZIP")
summary(nmix_covs_zip)
Warning message: “4 sites have been discarded because of missing data.”
Call:
pcount(formula = ~1 ~ elev + forest, data = umf, K = 20, mixture = "ZIP")
Abundance (log-scale):
Estimate SE z P(>|z|)
(Intercept) -1.097 0.304 -3.62 3.0e-04
elev -1.274 0.238 -5.35 8.9e-08
forest -0.614 0.181 -3.38 7.2e-04
Detection (logit-scale):
Estimate SE z P(>|z|)
0.261 0.235 1.11 0.266
Zero-inflation (logit-scale):
Estimate SE z P(>|z|)
-0.000622 0.345 -0.0018 0.999
AIC: 511.0829
Number of sites: 235
ID of sites removed due to NA: 12 69 118 146
optim convergence code: 0
optim iterations: 49
Bootstrap iterations: 0
Nmix.gof.test(nmix_covs_zip, nsim = 30)
Chi-square goodness-of-fit for N-mixture model of 'unmarkedFitPCount' class Observed chi-square statistic = 1388.716 Number of bootstrap samples = 30 P-value = 0 Quantiles of bootstrapped statistics: 0% 25% 50% 75% 100% 591 717 757 831 973 Estimate of c-hat = 1.8
umf@siteCovs$id <- factor(1:numSites(umf))
nmix_covs_pln <- pcount(~1~elev+forest+(1|id), umf, K = 20)
nmix_covs_pln
Warning message: “4 sites have been discarded because of missing data.”
Call:
pcount(formula = ~1 ~ elev + forest + (1 | id), data = umf, K = 20)
Abundance:
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 2.397 1.548
Fixed effects:
Estimate SE z P(>|z|)
(Intercept) -2.814 0.398 -7.08 1.47e-12
elev -1.376 0.294 -4.69 2.76e-06
forest -0.742 0.244 -3.04 2.40e-03
Detection:
Estimate SE z P(>|z|)
0.0218 0.226 0.0966 0.923
AIC: 472.0072
fl <- fitList(null=nmix_null, covs=nmix_covs, nb=nmix_covs_nb, zip=nmix_covs_zip)
modSel(fl)
Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.”
nPars AIC delta AICwt cumltvWt nb 5 481.09 0.00 1.0e+00 1.00 zip 5 511.08 29.99 3.1e-07 1.00 covs 4 530.74 49.64 1.7e-11 1.00 null 2 631.89 150.80 1.8e-33 1.00
unmarked we have to "integrate out" the latent abundancer <- ranef(nmix_covs_nb)
r
Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.”
Mean Mode 2.5% 97.5% [1,] 0.046984840 0 0 1 [2,] 0.041395270 0 0 1 [3,] 3.420860989 3 3 5 [4,] 0.017222689 0 0 0 [5,] 4.169619341 4 3 7 [6,] 0.021924662 0 0 0 [7,] 0.042814162 0 0 1 [8,] 0.041867701 0 0 1 [9,] 0.012209217 0 0 0 [10,] 0.019017597 0 0 0 [11,] 0.004868602 0 0 0 [12,] 0.042369172 0 0 1 [13,] 0.012868910 0 0 0 [14,] 0.030768078 0 0 1 [15,] 0.013552704 0 0 0 [16,] 0.006787598 0 0 0 [17,] 1.236817812 1 1 3 [18,] 0.035243109 0 0 1 [19,] 1.570388818 1 1 3 [20,] 0.004115638 0 0 0 [21,] 0.033620157 0 0 1 [22,] 0.020767930 0 0 0 [23,] 0.005920753 0 0 0 [24,] 0.033426545 0 0 1 [25,] 17.400203488 18 14 20 [26,] 0.036955636 0 0 1 [27,] 0.018010432 0 0 0 [28,] 0.007449491 0 0 0 [29,] 0.013189178 0 0 0 [30,] 1.090162820 1 1 2 [31,] 1.289206263 1 1 3 [32,] 0.045600556 0 0 1 [33,] 0.009331223 0 0 0 [34,] 0.021110646 0 0 0 [35,] 0.029696489 0 0 1 [36,] 0.015875879 0 0 0 [37,] 0.001977142 0 0 0 [38,] 0.006563806 0 0 0 [39,] 0.016031427 0 0 0 [40,] 1.128675124 1 1 2 [41,] 0.003354656 0 0 0 [42,] 0.004217539 0 0 0 [43,] 0.032385161 0 0 1 [44,] 0.025701403 0 0 0 [45,] 0.007267892 0 0 0 [46,] 0.029316871 0 0 1 [47,] 0.024588414 0 0 0 [48,] 1.161669381 1 1 2 [49,] 0.007361064 0 0 0 [50,] 0.015598795 0 0 0 [51,] 1.023594071 1 1 1 [52,] 0.016392929 0 0 0 [53,] 0.011800553 0 0 0 [54,] 0.039073809 0 0 1 [55,] 0.016360014 0 0 0 [56,] 0.012864006 0 0 0 [57,] 0.003536612 0 0 0 [58,] 0.039583086 0 0 1 [59,] 0.017127829 0 0 0 [60,] 2.387034184 2 2 4 [61,] 0.016996642 0 0 0 [62,] 0.045253570 0 0 1 [63,] 0.003969714 0 0 0 [64,] 1.146730304 1 1 2 [65,] 0.032294674 0 0 1 [66,] 0.035034680 0 0 1 [67,] 0.025138904 0 0 0 [68,] 7.423940054 7 5 11 [69,] 0.016289323 0 0 0 [70,] 0.007188429 0 0 0 [71,] 0.006177321 0 0 0 [72,] 0.036527938 0 0 1 [73,] 0.021187202 0 0 0 [74,] 0.025523797 0 0 0 [75,] 0.001530405 0 0 0 [76,] 0.001827743 0 0 0 [77,] 0.002538851 0 0 0 [78,] 0.006844917 0 0 0 [79,] 0.009829618 0 0 0 [80,] 0.020434962 0 0 0 [81,] 0.029248528 0 0 1 [82,] 0.036507427 0 0 1 [83,] 1.506835519 1 1 3 [84,] 0.016223472 0 0 0 [85,] 0.005659877 0 0 0 [86,] 0.003538239 0 0 0 [87,] 0.029832112 0 0 1 [88,] 0.217707648 0 0 2 [89,] 1.056936591 1 1 2 [90,] 0.019822430 0 0 0 [91,] 0.035276331 0 0 1 [92,] 0.005900424 0 0 0 [93,] 0.007346588 0 0 0 [94,] 0.033663017 0 0 1 [95,] 0.006091983 0 0 0 [96,] 0.041178389 0 0 1 [97,] 0.007692091 0 0 0 [98,] 0.012568247 0 0 0 [99,] 0.035022527 0 0 1 [100,] 1.158387684 1 1 2 [101,] 0.009788777 0 0 0 [102,] 0.002055558 0 0 0 [103,] 0.002153691 0 0 0 [104,] 0.005607716 0 0 0 [105,] 0.021262460 0 0 0 [106,] 1.236582664 1 1 3 [107,] 0.023211501 0 0 0 [108,] 0.005323028 0 0 0 [109,] 0.010051920 0 0 0 [110,] 0.040827652 0 0 1 [111,] 1.445577679 1 1 3 [112,] 0.001838545 0 0 0 [113,] 0.013663635 0 0 0 [114,] 0.022028534 0 0 0 [115,] 1.136665438 1 1 2 [116,] 0.044880345 0 0 1 [117,] 0.048964632 0 0 1 [118,] 0.015838401 0 0 0 [119,] 0.022368638 0 0 0 [120,] 1.304049417 1 1 3 [121,] 5.399131779 5 4 8 [122,] 0.021463651 0 0 0 [123,] 0.003422626 0 0 0 [124,] 0.023733418 0 0 0 [125,] 0.033330007 0 0 1 [126,] 0.044823512 0 0 1 [127,] 0.046542188 0 0 1 [128,] 1.289383433 1 1 3 [129,] 0.003637699 0 0 0 [130,] 3.598665972 3 3 6 [131,] 0.040943411 0 0 1 [132,] 0.046600089 0 0 1 [133,] 0.031201219 0 0 1 [134,] 2.479394103 2 2 4 [135,] 1.254382435 1 1 3 [136,] 0.006211823 0 0 0 [137,] 0.008327797 0 0 0 [138,] 0.009882644 0 0 0 [139,] 1.189427991 1 1 2 [140,] 0.042224498 0 0 1 [141,] 0.038790389 0 0 1 [142,] 0.028642725 0 0 1 [143,] 0.009590561 0 0 0 [144,] 0.011462420 0 0 0 [145,] 2.400826404 2 2 4 [146,] 0.035874812 0 0 1 [147,] 0.003721906 0 0 0 [148,] 0.001907943 0 0 0 [149,] 0.006156528 0 0 0 [150,] 5.955966574 6 4 9 [151,] 0.021926870 0 0 0 [152,] 0.010832767 0 0 0 [153,] 0.012729012 0 0 0 [154,] 2.263263834 2 2 4 [155,] 0.021731719 0 0 0 [156,] 0.021731719 0 0 0 [157,] 0.037377593 0 0 1 [158,] 0.019970002 0 0 0 [159,] 0.004523068 0 0 0 [160,] 0.037889244 0 0 1 [161,] 0.004224435 0 0 0 [162,] 0.010942642 0 0 0 [163,] 0.012693808 0 0 0 [164,] 0.007431375 0 0 0 [165,] 0.011382292 0 0 0 [166,] 0.010954780 0 0 0 [167,] 0.007485956 0 0 0 [168,] 0.022584640 0 0 0 [169,] 0.012335548 0 0 0 [170,] 0.007220592 0 0 0 [171,] 0.003451671 0 0 0 [172,] 0.018118934 0 0 0 [173,] 4.610648264 4 3 7 [174,] 1.540941761 1 1 3 [175,] 0.010412973 0 0 0 [176,] 0.028859646 0 0 1 [177,] 0.002868679 0 0 0 [178,] 2.882512706 2 2 5 [179,] 0.015201274 0 0 0 [180,] 0.006996819 0 0 0 [181,] 0.005142620 0 0 0 [182,] 0.005795895 0 0 0 [183,] 0.009830855 0 0 0 [184,] 0.005275956 0 0 0 [185,] 0.005025610 0 0 0 [186,] 1.125139987 1 1 2 [187,] 0.018404511 0 0 0 [188,] 1.565865151 1 1 3 [189,] 0.003169450 0 0 0 [190,] 0.002021378 0 0 0 [191,] 0.014749085 0 0 0 [192,] 1.578976603 1 1 3 [193,] 0.019504834 0 0 0 [194,] 0.010452845 0 0 0 [195,] 0.016969642 0 0 0 [196,] 1.079748248 1 1 2 [197,] 0.041302157 0 0 1 [198,] 3.010439203 3 2 5 [199,] 2.071956916 2 2 3 [200,] 1.035223644 1 1 2 [201,] 0.014757212 0 0 0 [202,] 0.008609226 0 0 0 [203,] 0.014538496 0 0 0 [204,] 0.009595070 0 0 0 [205,] 0.012969280 0 0 0 [206,] 0.007582165 0 0 0 [207,] 0.008532701 0 0 0 [208,] 0.002666074 0 0 0 [209,] 1.008689515 1 1 1 [210,] 0.001487545 0 0 0 [211,] 0.001403987 0 0 0 [212,] 0.001457889 0 0 0 [213,] 0.007188429 0 0 0 [214,] 0.002117769 0 0 0 [215,] 0.002981001 0 0 0 [216,] 0.010093609 0 0 0 [217,] 0.010093609 0 0 0 [218,] 0.006760214 0 0 0 [219,] 0.008271254 0 0 0 [220,] 0.004494322 0 0 0 [221,] 0.003653496 0 0 0 [222,] 0.006760214 0 0 0 [223,] 0.017933077 0 0 0 [224,] 0.012278793 0 0 0 [225,] 0.021055879 0 0 0 [226,] 0.014881020 0 0 0 [227,] 0.021524680 0 0 0 [228,] 0.010393372 0 0 0 [229,] 0.012278793 0 0 0 [230,] 0.014881020 0 0 0 [231,] 0.017135632 0 0 0 [232,] 0.017933077 0 0 0 [233,] 0.010093609 0 0 0 [234,] 0.010093609 0 0 0 [235,] 0.008271254 0 0 0
head(bup(r))
head(bup(r, stat='mode'))
Use unmarked to fit the following model.
mallard.obs, you will need to modify the unmarkedFrame)Check goodness-of-fit.
Compare with the other models using model selection.
I(elev^2)
umf <- unmarkedFramePCount(y = mallard.y, siteCovs = mallard.site, obsCovs = mallard.obs)
nmix_obs <- pcount(~ivel~elev+forest+I(elev^2), umf, K = 20)
nmix_obs
Warning message: “4 sites have been discarded because of missing data.”
Call:
pcount(formula = ~ivel ~ elev + forest + I(elev^2), data = umf,
K = 20)
Abundance:
Estimate SE z P(>|z|)
(Intercept) -1.88792 0.283 -6.6696 2.57e-11
elev -1.35446 0.231 -5.8597 4.64e-09
forest -0.73976 0.170 -4.3539 1.34e-05
I(elev^2) 0.00943 0.245 0.0385 9.69e-01
Detection:
Estimate SE z P(>|z|)
(Intercept) 0.420 0.214 1.97 0.0493
ivel 0.423 0.171 2.47 0.0136
AIC: 528.2189
suppressWarnings(Nmix.gof.test(nmix_obs, nsim=30))
Chi-square goodness-of-fit for N-mixture model of 'unmarkedFitPCount' class Observed chi-square statistic = 1479.577 Number of bootstrap samples = 30 P-value = 0 Quantiles of bootstrapped statistics: 0% 25% 50% 75% 100% 498 587 630 718 1010 Estimate of c-hat = 2.23
# re-fit to make sure all models use the same dataset
nmix_null <- pcount(~1~1, umf, K=20)
nmix_covs <- pcount(~1~forest+elev, umf, K=20)
fl <- fitList(nmix_null=nmix_null, nmix_covs=nmix_covs, nmix_obs=nmix_obs)
modSel(fl)
Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.” Warning message: “4 sites have been discarded because of missing data.”
nPars AIC delta AICwt cumltvWt nmix_obs 6 528.22 0.00 7.8e-01 0.78 nmix_covs 4 530.74 2.52 2.2e-01 1.00 nmix_null 2 631.89 103.67 2.4e-23 1.00
A common situation:
Potential solution: "stacking"!
Potential issue:
Could ignore this if you don't expect much correlation in abundance at a site among years, or add a site-level random effect.
We don't have a convenient dataset with multiple years of data, so I'll just simulate some based on our mallard dataset and our model fit.
new_data <- unmarked::simulate(nmix_covs, 2, na.rm=FALSE)
lapply(new_data, head)
| 1 | 1 | 1 |
| 0 | 0 | 0 |
| 0 | 0 | 0 |
| 0 | 0 | 0 |
| 2 | 2 | 1 |
| 2 | 2 | 1 |
| 0 | 0 | 1 |
| 0 | 0 | 0 |
| 0 | 0 | 0 |
| 1 | 0 | 0 |
| 0 | 0 | 0 |
| 0 | 0 | 0 |
Original (real) data is year 1, and simulated datasets are years 2 and 3
y_all <- rbind(umf@y, new_data[[1]], new_data[[2]])
# "Site" (site-year) covariates
sc_all <- rbind(umf@siteCovs, umf@siteCovs, umf@siteCovs)
sc_all$site <- factor(rep(1:numSites(umf), 3))
sc_all$year <- rep(1:3, each=numSites(umf))
head(sc_all)
| elev | length | forest | site | year | |
|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <fct> | <int> | |
| 1 | -1.173 | 0.801 | -1.156 | 1 | 1 |
| 2 | -1.127 | 0.115 | -0.501 | 2 | 1 |
| 3 | -0.198 | -0.479 | -0.101 | 3 | 1 |
| 4 | -0.105 | 0.315 | 0.008 | 4 | 1 |
| 5 | -1.034 | -1.102 | -1.193 | 5 | 1 |
| 6 | -0.848 | 0.741 | 0.917 | 6 | 1 |
# Obs covariates
oc_all <- rbind(umf@obsCovs, umf@obsCovs, umf@obsCovs)
Make the combined unmarkedFrame
umf_all <- unmarkedFramePCount(y=y_all, siteCovs=sc_all, obsCovs=oc_all)
summary(umf_all)
unmarkedFrame Object
717 sites
Maximum number of observations per site: 3
Mean number of observations per site: 2.92
Sites with at least one detection: 145
Tabulation of y observations:
0 1 2 3 4 5 7 10 12 <NA>
1778 230 56 15 9 2 1 1 1 58
Site-level covariates:
elev length forest site
Min. :-1.436000 Min. :-4.945000 Min. :-1.2650000 1 : 3
1st Qu.:-0.972000 1st Qu.:-0.563000 1st Qu.:-0.9740000 2 : 3
Median :-0.198000 Median : 0.045000 Median :-0.0650000 3 : 3
Mean :-0.000046 Mean :-0.000029 Mean : 0.0000669 4 : 3
3rd Qu.: 0.994000 3rd Qu.: 0.626000 3rd Qu.: 0.8080000 5 : 3
Max. : 2.434000 Max. : 2.255000 Max. : 2.2990000 6 : 3
(Other):699
year
Min. :1
1st Qu.:1
Median :2
Mean :2
3rd Qu.:3
Max. :3
Observation-level covariates:
ivel date
Min. :-1.75300 Min. :-2.90400
1st Qu.:-0.66600 1st Qu.:-1.11900
Median :-0.13900 Median :-0.11900
Mean : 0.00002 Mean : 0.00007
3rd Qu.: 0.54900 3rd Qu.: 1.31000
Max. : 5.98000 Max. : 3.81000
NA's :156 NA's :126
Not accounting for pseudoreplicated sites:
mod_stack <- pcount(~1~elev+forest+year, umf_all, K=20)
mod_stack
Warning message: “4 sites have been discarded because of missing data.”
Call:
pcount(formula = ~1 ~ elev + forest + year, data = umf_all, K = 20)
Abundance:
Estimate SE z P(>|z|)
(Intercept) -1.7821 0.2054 -8.677 4.06e-18
elev -1.2438 0.1211 -10.275 9.17e-25
forest -0.6809 0.0914 -7.451 9.27e-14
year -0.0344 0.0819 -0.421 6.74e-01
Detection:
Estimate SE z P(>|z|)
0.488 0.107 4.57 4.83e-06
AIC: 1534.382
Accounting for pseudoreplicated sites with a random effect:
mod_stack2 <- pcount(~1~elev+forest+year+(1|site), umf_all, K=20)
mod_stack2
Warning message: “4 sites have been discarded because of missing data.”
Call:
pcount(formula = ~1 ~ elev + forest + year + (1 | site), data = umf_all,
K = 20)
Abundance:
Random effects:
Groups Name Variance Std.Dev.
site (Intercept) 0.271 0.52
Fixed effects:
Estimate SE z P(>|z|)
(Intercept) -1.8956 0.2148 -8.823 1.11e-18
elev -1.2477 0.1318 -9.470 2.80e-21
forest -0.6738 0.1056 -6.383 1.73e-10
year -0.0395 0.0819 -0.483 6.29e-01
Detection:
Estimate SE z P(>|z|)
0.449 0.111 4.06 4.98e-05
AIC: 1521.633
ubms¶library(ubms)
nmix_stan <- stan_pcount(~1~elev+forest, umf, K = 20, chains = 3, iter = 300, refresh = 0)
nmix_stan
Attaching package: ‘ubms’
The following objects are masked from ‘package:unmarked’:
fitList, projected
The following object is masked from ‘package:base’:
gamma
Warning message:
“Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess”
Warning message:
“Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess”
Warning message:
“Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
”
Call:
stan_pcount(formula = ~1 ~ elev + forest, data = umf, K = 20,
chains = 3, iter = 300, refresh = 0)
Abundance (log-scale):
Estimate SD 2.5% 97.5% n_eff Rhat
(Intercept) -1.954 0.271 -2.53 -1.490 194 1.01
elev -1.353 0.266 -1.89 -0.921 218 1.00
forest -0.756 0.176 -1.13 -0.431 383 1.00
Detection (logit-scale):
Estimate SD 2.5% 97.5% n_eff Rhat
0.48 0.191 0.123 0.828 326 1
LOOIC: 537.534
Runtime: 12.186 sec
Nest <- posterior_predict(nmix_stan, "z")
dim(Nest)
hist(apply(Nest, 1, sum, na.rm=TRUE), main="sumN", xlab="sumN")
Demonstration of how prior settings can affect coefficient estimates
Set the prior_intercept_state argument in stan_pcount
vary_prior <- function(prior){
suppressWarnings(stan_pcount(~1~elev+forest, umf, K = 20, chains = 3, iter = 300, refresh = 0,
prior_intercept_state = prior))
}
priors <- list(p5=normal(0, 5), p1=normal(0, 1), p25=normal(0, 2.5), p05=normal(0, 0.5), p25=normal(0, 0.1))
mods <- lapply(priors, vary_prior)
stats <- t(sapply(mods, function(x) as.numeric(summary(x, 'state')[1,c(1,4,8)])))
plot(1:5, stats[,1], ylim=c(-2.5, 0), pch=19, xaxt='n',
xlab="", ylab="Intercept estimate")
axis(1, at=1:5, labels=c("N(0,5)","N(0,2.5)","N(0,1)", "N(0,0.5)","N(0,0.1)"))
segments(1:5, stats[,2], 1:5, stats[,3])
NIMBLE is a general-purpose model fitting enginelibrary(nimble)
nimble version 0.13.2 is loaded.
For more information on NIMBLE and a User Manual,
please visit https://R-nimble.org.
Note for advanced users who have written their own MCMC samplers:
As of version 0.13.0, NIMBLE's protocol for handling posterior
predictive nodes has changed in a way that could affect user-defined
samplers in some situations. Please see Section 15.5.1 of the User Manual.
Attaching package: ‘nimble’
The following object is masked from ‘package:unmarked’:
simulate
The following object is masked from ‘package:stats’:
simulate
nmix_code <- nimbleCode({
# State model
for (i in 1:M){
N[i] ~ dpois(lambda[i])
log(lambda[i]) <- beta_0 + beta_ele * elev[i] + beta_for * forest[i]
}
beta_0 ~ dunif(-10, 10)
beta_ele ~ dnorm(0, 0.01)
beta_for ~ dnorm(0, 0.01)
# Detection model
for (i in 1:M){
for (j in 1:J){
y[i,j] ~ dbinom(p, N[i])
}
}
p ~ dunif(0,1)
sumN <- sum(N[])
})
pars <- c("beta_0", "beta_ele", "beta_for", "p")
inits <- function() list(N = Ninit)
nimble_covs <- nimbleMCMC(code = nmix_code,
constants = jags_data,
inits = inits,
monitors = pars,
nchains = 3,
niter=3000,
nburnin=2000,
summary = TRUE,
samplesAsCodaMCMC = TRUE)
Defining model
[Note] 'ivel' is provided in 'constants' but not used in the model code and is being ignored.
[Note] Using 'y' (given within 'constants') as data.
Building model
Setting data and initial values
Running calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.
Checking model sizes and dimensions
[Note] This model is not fully initialized. This is not an error.
To see which variables are not initialized, use model$initializeInfo().
For more information on model initialization, see help(modelInitialization).
Checking model calculations
[Note] NAs were detected in model variables: beta_0, logProb_beta_0, beta_ele, logProb_beta_ele, beta_for, logProb_beta_for, p, logProb_p, lambda, logProb_N, N, logProb_y, y, sumN.
Compiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
running chain 1...
warning: problem initializing stochastic node N[40]: logProb less than -1e12. warning: problem initializing stochastic node N[64]: logProb less than -1e12. warning: problem initializing stochastic node N[72]: logProb less than -1e12. warning: problem initializing stochastic node N[86]: logProb less than -1e12. warning: problem initializing stochastic node N[204]: logProb less than -1e12. warning: problem initializing stochastic node N[206]: logProb less than -1e12. warning: problem initializing stochastic node N[207]: logProb less than -1e12. warning: problem initializing stochastic node N[210]: logProb less than -1e12. warning: problem initializing stochastic node N[211]: logProb less than -1e12. warning: problem initializing stochastic node N[217]: logProb less than -1e12. warning: problem initializing stochastic node N[220]: logProb less than -1e12. warning: problem initializing stochastic node N[221]: logProb less than -1e12. warning: problem initializing stochastic node N[222]: logProb less than -1e12. warning: problem initializing stochastic node N[223]: logProb less than -1e12. warning: problem initializing stochastic node N[224]: logProb less than -1e12. warning: problem initializing stochastic node N[225]: logProb less than -1e12. warning: problem initializing stochastic node N[226]: logProb less than -1e12. warning: problem initializing stochastic node N[228]: logProb less than -1e12. warning: problem initializing stochastic node N[230]: logProb less than -1e12. warning: problem initializing stochastic node N[233]: logProb less than -1e12. warning: problem initializing stochastic node N[234]: logProb less than -1e12. warning: problem initializing stochastic node N[237]: logProb less than -1e12. warning: problem initializing stochastic node N[238]: logProb less than -1e12. warning: problem initializing stochastic node N[239]: logProb less than -1e12. |-------------|-------------|-------------|-------------| |Warning: slice sampler reached maximum number of contractions for 'N[224]'. Current parameter value is 16. Warning: slice sampler reached maximum number of contractions for 'N[64]'. Current parameter value is 40. Warning: slice sampler reached maximum number of contractions for 'N[64]'. Current parameter value is 40. -------------------------------------------------------|
running chain 2...
warning: problem initializing stochastic node y[69, 1]: logProb is -Inf. warning: problem initializing stochastic node y[69, 2]: logProb is -Inf. warning: problem initializing stochastic node y[69, 3]: logProb is -Inf. warning: problem initializing stochastic node y[146, 3]: logProb is -Inf. |-------------|-------------|-------------|-------------| |-------------------------------------------------------|
running chain 3...
warning: problem initializing stochastic node y[222, 3]: logProb is -Inf. |-------------|-------------|-------------|-------------| |-------------------------------------------------------|
coda::traceplot(nimble_covs$samples)
nimble_covs$summary$all.chains
| Mean | Median | St.Dev. | 95%CI_low | 95%CI_upp | |
|---|---|---|---|---|---|
| beta_0 | -1.9311467 | -1.9154466 | 0.23295557 | -2.4262519 | -1.5196816 |
| beta_ele | -1.3398020 | -1.3210108 | 0.21367156 | -1.8186875 | -0.9539358 |
| beta_for | -0.7354201 | -0.7312231 | 0.16262465 | -1.0655075 | -0.4295455 |
| p | 0.6157582 | 0.6192895 | 0.04458034 | 0.5247495 | 0.7021781 |
NIMBLE macros¶library(nimbleMacros)
nimbleOptions(enableMacroComments = FALSE)
priors <- nimbleMacros::priors
nimble_lp <- nimbleCode({
lambda[1:M] <- linPred(~elev+forest, link = log)
})
Attaching package: ‘nimbleMacros’
The following object is masked _by_ ‘.GlobalEnv’:
priors
mod <- nimbleModel(nimble_lp, constants = jags_data)
mod$getCode()
Defining model
[Note] 'y' is provided in 'constants' but not used in the model code and is being ignored.
[Note] 'ivel' is provided in 'constants' but not used in the model code and is being ignored.
[Note] 'J' is provided in 'constants' but not used in the model code and is being ignored.
Building model
Running calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.
Checking model sizes and dimensions
[Note] This model is not fully initialized. This is not an error.
To see which variables are not initialized, use model$initializeInfo().
For more information on model initialization, see help(modelInitialization).
{
for (i_1 in 1:M) {
log(lambda[i_1]) <- beta_Intercept + beta_elev * elev[i_1] +
beta_forest * forest[i_1]
}
beta_Intercept ~ dunif(-100, 100)
beta_elev ~ dnorm(0, sd = 100)
beta_forest ~ dnorm(0, sd = 100)
}
jags_data$group <- factor(sample(letters[1:5], jags_data$M, replace=TRUE))
nimble_lp2 <- nimbleCode({
lambda[1:M] <- linPred(~forest + (1|group), link = log)
})
mod <- nimbleModel(nimble_lp2, constants = jags_data)
mod$getCode()
Defining model
[Note] 'y' is provided in 'constants' but not used in the model code and is being ignored.
[Note] 'ivel' is provided in 'constants' but not used in the model code and is being ignored.
[Note] 'elev' is provided in 'constants' but not used in the model code and is being ignored.
[Note] 'J' is provided in 'constants' but not used in the model code and is being ignored.
Building model
Running calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.
Checking model sizes and dimensions
[Note] This model is not fully initialized. This is not an error.
To see which variables are not initialized, use model$initializeInfo().
For more information on model initialization, see help(modelInitialization).
{
for (i_1 in 1:M) {
log(lambda[i_1]) <- beta_Intercept + beta_forest * forest[i_1] +
beta_group[group[i_1]]
}
beta_Intercept ~ dunif(-100, 100)
beta_forest ~ dnorm(0, sd = 100)
sd_group ~ T(dt(0, 0.01, 1), 0, )
for (i_2 in 1:5) {
beta_group[i_2] ~ dnorm(0, sd = sd_group)
}
}
nimble_covs2 <- nimbleCode({
y[1:M, 1:J] ~ nmixture(~1~elev+forest+(1|group))
})
nimble_covs <- nimbleModel(nimble_covs2, constants = jags_data)
nimble_covs$getCode()
Defining model
[Note] 'ivel' is provided in 'constants' but not used in the model code and is being ignored.
[Note] Using 'y' (given within 'constants') as data.
Building model
Setting data and initial values
Running calculate on model
[Note] Any error reports that follow may simply reflect missing values in model variables.
Checking model sizes and dimensions
[Note] This model is not fully initialized. This is not an error.
To see which variables are not initialized, use model$initializeInfo().
For more information on model initialization, see help(modelInitialization).
{
for (i_1 in 1:M) {
log(lam[i_1]) <- state_Intercept + state_elev * elev[i_1] +
state_forest * forest[i_1] + state_group[group[i_1]]
}
state_Intercept ~ dunif(-10, 10)
state_elev ~ dnorm(0, 0.01)
state_forest ~ dnorm(0, 0.01)
state_sd_group ~ T(dt(0, 0.01, 1), 0, )
for (i_2 in 1:5) {
state_group[i_2] ~ dnorm(0, sd = state_sd_group)
}
for (i_3 in 1:M) {
for (i_4 in 1:J) {
logit(p[i_3, i_4]) <- det_Intercept
}
}
det_Intercept ~ dunif(-10, 10)
for (i_5 in 1:M) {
N[i_5] ~ dpois(lam[i_5])
}
for (i_6 in 1:M) {
for (i_7 in 1:J) {
y[i_6, i_7] ~ dbinom(N[i_6], p[i_6, i_7])
}
}
}