Bat insect abundance using exclosures. Nmixture model using JAGS
Francisco Sanchez, Diego J. Lizcano, et al.
12/05/2022
We used hierarchical Binomial Mixture Model to analyzed a two categorical covariates that act on the ecological process and one continuous covariate on the observation process.
We model the effect of the covariates on the scale of the log and the logit, as is customary for generalized linear models.
Code based on the BPA book chapter 12.2.2. Model Introducing covariates. Binomial-mix model data.
Our model is a variation of the original code by Kery´s book, but incorporating variation in the ecological process in the form of intercept and slope in two levels varying per site. Our response variable was the insect abundance in a treatment which consisted of an enclosure experiment to keep out bats and controls.
we sampled five consecutive nights as repeated counts.
library(jagsUI)
library(parallel)
library(MCMCvis)
library(AHMbook)
library (readxl) # read excel data
library (tidyverse) # Data handling
library (bayesplot)
# source("D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/R/DBDA2E-utilities.R")Algebra of the model
Our hierarchical model included in the ecological process the counts of expected abundance (\(\lambda\)) as a realization of a Poisson process affected by the Treatment (Excluded, no excluded) and varying in the intercept and slope with the two levels of the type of system (silvopastoral vs conventional). The observation process model included the the detection probability (\(p\)) affected by the percentage of moon light as a binomial distribution.
\[ N_{i} \sim {\sf Pois}(\lambda_{i}) \tag{1: ecological process} \] \[ log(\lambda_{i}) = \beta_0[sistema_{i}] + \beta_1[sistema_{i}] * Tratamiento_{i} \tag{1: 1} \]
\[ y_{ij} | N_{i} \sim {\sf Binom} (N_{i} *p) \tag{2: observational process} \]
\[ logit(p_{ij}) = \alpha_0 * Moon_{ij} + \alpha_1 * (Moon_{ij})^{2} \tag{2: 2} \]
Cow parasites
Load Data
bat_exclusure <- read_excel("D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/data/datos_analisis_servicios2022.xlsx",
sheet = "vaca")
luna <- read_excel("D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/data/Datos_Bat_exclusion_12dic2019_tydy_fix_unmarked_no_na.xlsx",
sheet = "luna")JAGS model
# Specify model in BUGS language
sink("model_7_factors.jags")
cat("
model {
# Priors
for(k in 1:2){ # Loop over 2 levels of sistema (silvopastoril and convencional)
# alpha0[k] ~ dunif(-10, 10) # Detection intercepts
# alpha1[k] ~ dunif(-10, 10) # Detection slopes
# alpha2[k] ~ dunif(-10, 10) # Detection slopes
beta0[k] ~ dunif(-10, 10) # Abundance intercepts
beta1[k] ~ dunif(-10, 10) # Abundance slopes
# beta2[k] ~ dunif(-10, 10) # Abundance slopes
}
# # Priors
alpha0 ~ dunif(-10, 10)
alpha1 ~ dunif(-10, 10)
# alpha2 ~ dunif(-10, 10)
# beta0 ~ dunif(-10, 10)
# beta1 ~ dunif(-10, 10)
# beta2 ~ dunif(-10, 10)
# beta3 ~ dunif(-10, 10)
# Likelihood
# Ecological model for true abundance
for (i in 1:R){
N[i] ~ dpois(lambda[i]) # 1. Distribution for abundance part
log(lambda[i]) <- beta0[Sistema[i]] + beta1[Sistema[i]] * Tratamiento[i] # + beta2 * Tratamiento[i] * Sistema[i]
site.abundance[i] <- N[i] # Abundance per site
# Observation model for replicated counts
for (j in 1:T){
y[i,j] ~ dbin(p[i,j], N[i])
p[i,j] <- exp(lp[i,j])/(1+exp(lp[i,j]))
lp[i,j] <- alpha0 * luna[i,j] + alpha1 * pow(luna[i,j],2)
} #j
} #i
# Posterior predictive distributions of chi2 discrepancy
for (i in 1:R) {
for (j in 1:T) {
C.sim[i,j] ~ dbin(p[i,j], N[i]) # Create new data set under model
e.count[i,j] <- N[i] * p[i,j] # Expected datum
# Chi-square discrepancy for the actual data
chi2.actual[i,j] <- pow((y[i,j] - e.count[i,j]),2) / (e.count[i,j]+0.0001)
# Chi-square discrepancy for the simulated ('perfect') data
chi2.sim[i,j] <- pow((C.sim[i,j] - e.count[i,j]),2) / (e.count[i,j]+0.0001)
# Add small value e to denominator to avoid division by zero
}
}
# Add up individual chi2 values for overall fit statistic
fit.actual <- sum(chi2.actual[,]) # Fit statistic for actual data set
fit.sim <- sum(chi2.sim[,]) # Fit statistic for a fitting model
c.hat <- fit.actual / fit.sim # c-hat estimate
bpv <- step(fit.sim-fit.actual) # Bayesian p-value
# Derived quantities: functions of latent variables and predictions
totalN <- sum(N[]) # Total pop. size across all sites
mean.abundance[1] <- sum(site.abundance[1:44])/44
mean.abundance[2] <- sum(site.abundance[45:84])/44
Dif.mean.abundance <- (sum(site.abundance[1:44]))/44 - (sum(site.abundance[45:84]))/44
Ntratamiento[1] <- sum(N[1:44]) # Total pop for sites in encierro
Ntratamiento[2] <- sum(N[45:88]) # Total pop for sites in sin encierro
N_Excl_Conv <- (sum(N[1:4]) + sum(N[9:11]) + sum(N[15:19]) + sum(N[25:29]) + sum(N[35:39]))/24
N_Excl_Silv <- (sum(N[5:8]) + sum(N[12:14]) + sum(N[20:24]) + sum (N[30:34]) + sum(N[40:44]))/20
N_noExcl_Conv <- (sum(N[45:48]) + sum(N[53:55]) + sum(N[59:63]) + sum(N[69:73]) + sum(N[79:83]))/23
N_noExcl_Silv <- (sum(N[49:52]) + sum(N[56:58]) + sum(N[64:68]) + sum(N[74:78]) + sum(N[84:88]))/21
D_Silv <- N_Excl_Silv - N_noExcl_Silv
D_Conv <- N_Excl_Conv - N_noExcl_Conv
Difftratamiento <- sum(N[1:44]) - sum(N[45:88]) # Diferencia encierro-sin
for(s in 1:2){ # Predictions of lambda Sistema
for(t in 1:2){ # abundance for each Tratamiento of sistema factor
lam.pred.exclosure[s, t] <- exp(beta0[s] + beta1[s] * Tratamiento[t])
# logit(p.pred[s, t]) <- alpha0[s] + alpha1[s] * luna[t]
} #t
} #s
}
# ",fill = TRUE)
sink()Run model
# Bundle data
y <- as.matrix(bat_exclusure[,9:13])
win.data <- list(y = y, R = nrow(y), T = ncol(y),
luna = as.matrix(luna[,8:12]),
Sistema = as.numeric(factor(bat_exclusure$Sistema)),
Tratamiento=as.numeric(factor(bat_exclusure$Tratamiento)))
# Initial values
Nst <- apply(y, 1, max) + 1 # Important to give good inits for latent N
inits <- function() list(N = Nst,
alpha0 = rnorm(1),
alpha1 = rnorm(1),
# alpha2 = rnorm(1),
beta0 = rnorm(2),
beta1 = rnorm(2)#,
# beta2 = runif(1, -1, 1),
# beta3 = runif(1, -1, 1)
)
# Parameters monitored
params <- c("totalN",
"alpha0", "alpha1", #"alpha2",
"beta0", "beta1", # "beta2", #"beta3",
"lam.pred",
"lam.pred.exclosure",
"N_Excl_Conv",
"N_Excl_Silv",
"N_noExcl_Conv",
"N_noExcl_Silv",
"D_Silv",
"D_Conv",
#"site.abundance",
"mean.abundance",
"Dif.mean.abundance",
"Ntratamiento",
"Difftratamiento",
"fit.actual", "fit.sim", "bpv", "c.hat")
# MCMC settings
ni <- 260000
nt <- 50
nb <- 10000
nc <- 3
# Call JAGS from R (BRT 4 min)
out3 <- jags(win.data,
inits,
params,
"D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/R/model_7_factors.jags",
n.chains = nc,
n.thin = nt,
n.iter = ni,
n.burnin = nb,
parallel = TRUE,
n.cores=3)##
## Processing function input.......
##
## Done.
##
## Beginning parallel processing using 3 cores. Console output will be suppressed.
##
## Parallel processing completed.
##
## Calculating statistics.......
##
## Done.# Save
saveRDS(out3, "D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/result_models/model_7a_factors.RData")
# Summarize posteriors
print(out3, dig = 3)## JAGS output for model 'D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/R/model_7_factors.jags', generated by jagsUI.
## Estimates based on 3 chains of 260000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 10000 iterations and thin rate = 50,
## yielding 15000 total samples from the joint posterior.
## MCMC ran in parallel for 9.72 minutes at time 2022-05-20 23:53:57.
##
## mean sd 2.5% 50% 97.5% overlap0
## totalN 727.907 36.744 660.000 726.000 803.000 FALSE
## alpha0 -4.912 0.419 -5.720 -4.917 -4.074 FALSE
## alpha1 3.708 0.438 2.837 3.713 4.555 FALSE
## beta0[1] 1.946 0.227 1.504 1.946 2.388 FALSE
## beta0[2] 2.711 0.190 2.335 2.712 3.076 FALSE
## beta1[1] 0.019 0.141 -0.258 0.019 0.292 TRUE
## beta1[2] -0.333 0.123 -0.578 -0.331 -0.093 FALSE
## lam.pred.exclosure[1,1] 7.172 0.771 5.781 7.138 8.779 FALSE
## lam.pred.exclosure[2,1] 10.828 0.981 9.002 10.792 12.851 FALSE
## lam.pred.exclosure[1,2] 7.172 0.771 5.781 7.138 8.779 FALSE
## lam.pred.exclosure[2,2] 10.828 0.981 9.002 10.792 12.851 FALSE
## N_Excl_Conv 6.574 0.473 5.708 6.542 7.583 FALSE
## N_Excl_Silv 11.915 0.759 10.550 11.850 13.500 FALSE
## N_noExcl_Conv 6.991 0.514 6.043 6.957 8.043 FALSE
## N_noExcl_Silv 8.145 0.576 7.095 8.143 9.333 FALSE
## D_Silv 3.770 0.749 2.324 3.752 5.271 FALSE
## D_Conv -0.417 0.572 -1.558 -0.415 0.699 TRUE
## mean.abundance[1] 9.002 0.503 8.091 8.977 10.045 FALSE
## mean.abundance[2] 6.679 0.404 5.932 6.659 7.500 FALSE
## Dif.mean.abundance 2.322 0.451 1.455 2.318 3.227 FALSE
## Ntratamiento[1] 396.067 22.148 356.000 395.000 442.000 FALSE
## Ntratamiento[2] 331.840 19.752 295.000 331.000 372.000 FALSE
## Difftratamiento 64.227 20.278 25.000 64.000 104.000 FALSE
## fit.actual 804.771 30.195 750.143 803.361 868.001 FALSE
## fit.sim 325.473 22.403 283.831 324.728 371.110 FALSE
## bpv 0.000 0.000 0.000 0.000 0.000 FALSE
## c.hat 2.484 0.184 2.144 2.475 2.862 FALSE
## deviance 1840.218 16.357 1809.387 1839.729 1873.460 FALSE
## f Rhat n.eff
## totalN 1.000 1.000 15000
## alpha0 1.000 1.000 15000
## alpha1 1.000 1.000 15000
## beta0[1] 1.000 1.000 15000
## beta0[2] 1.000 1.000 15000
## beta1[1] 0.555 1.000 12469
## beta1[2] 0.998 1.000 15000
## lam.pred.exclosure[1,1] 1.000 1.000 15000
## lam.pred.exclosure[2,1] 1.000 1.000 15000
## lam.pred.exclosure[1,2] 1.000 1.000 15000
## lam.pred.exclosure[2,2] 1.000 1.000 15000
## N_Excl_Conv 1.000 1.000 15000
## N_Excl_Silv 1.000 1.000 15000
## N_noExcl_Conv 1.000 1.000 15000
## N_noExcl_Silv 1.000 1.000 15000
## D_Silv 1.000 1.000 15000
## D_Conv 0.767 1.000 15000
## mean.abundance[1] 1.000 1.000 15000
## mean.abundance[2] 1.000 1.000 15000
## Dif.mean.abundance 1.000 1.000 15000
## Ntratamiento[1] 1.000 1.000 15000
## Ntratamiento[2] 1.000 1.000 15000
## Difftratamiento 1.000 1.000 15000
## fit.actual 1.000 1.000 6967
## fit.sim 1.000 1.000 15000
## bpv 1.000 NaN 1
## c.hat 1.000 1.000 15000
## deviance 1.000 1.001 7171
##
## 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 = 133.8 and DIC = 1973.972
## DIC is an estimate of expected predictive error (lower is better). MCMCtrace(out3,
params = c("alpha0", "alpha1", "alpha2"),#, "beta0", "beta1", "beta2"),
#ISB = FALSE,
#exact = TRUE,
pdf = FALSE)#,
# main_den= c("Inclusion probability",
# "Mean capture",
# "Variation of phi"))
MCMCtrace(out3,
params = c("beta0", "beta1", "beta2"),
#ISB = FALSE,
#exact = TRUE,
pdf = FALSE)#,
# main_den= c("Inclusion probability",
# "Mean capture",
# "Variation of phi"))
################ GOF test
# A fitting model has a Bayesian p-value around 0.5 and a model that does not fit has a p-value close to or equal to 0 or 1. However, there are no rules to decide on a threshold value that distinguishes a fitting from a nonfitting model.
#c.hat is estimated at only a little more than 1, and the Bayesian p-value of 0.070 almost suggests a fitting model.
ppc.plot(out3) # Results from a posterior predictive check of GoF (based on Chi-square discrepancy) for simulated data under an N-mixture model variant
################
MCMCplot(out3,
params = c("alpha0", "alpha1", "alpha2", "beta0", "beta1", "beta2", "beta3"),
ci = c(50, 90),
ref_ovl = TRUE,
guide_lines = TRUE,
main= "Parameters")
MCMCplot(out3,
params = c("Ntratamiento", "Difftratamiento" ),
ci = c(50, 90),
ref_ovl = TRUE,
guide_lines = TRUE,
main= "Parameters Total Population")
MCMCplot(out3,
params = c("mean.abundance", "Dif.mean.abundance" ),
ref_ovl = TRUE,
ci = c(50, 90), #credible intercval
guide_lines = TRUE,
main= "Abundance per trap")
MCMCplot(out3,
params = c( "N_noExcl_Conv","N_Excl_Conv", "N_noExcl_Silv", "N_Excl_Silv"),
labels=c("No Exclosure Conventional",
"Exclosure Conventional",
"No Exclosure Silvopastoral",
"Exclosure Silvopastoral"),
ci = c(50, 90),
ref_ovl = TRUE,
guide_lines = TRUE,
main= "Abundance per trap")
# # Plot posteriors
# par(mfrow = c(2, 1))
# hist(log(out$sims.list$Ntratamiento[,1]), col = "gray", main = "", xlab = "Exclosure", las = 1)
# abline(v = median(log(out$sims.list$Ntratamiento[,1])), lwd = 3, col = "red")
# hist(log(out$sims.list$Ntratamiento[,2]), col = "gray", main = "", xlab = "No Eclosure", las = 1)
# abline(v = median(log(out$sims.list$Ntratamiento[,2])), lwd = 3, col = "red")
#
#
df <- data.frame(
Tratamiento=factor(rep(c("Exclosure", "No Exclosure"), each=length(out3$sims.list$N_Excl_Conv) )),
Abundance=c(((out3$sims.list$N_Excl_Conv)),
((out3$sims.list$N_noExcl_Conv)))
)
df2 <- data.frame(
Tratamiento=factor(rep(c("Diferencia"), each=length(out3$sims.list$Difftratamiento))),
Abundance=c((out3$sims.list$D_Conv))# ,
#((out$sims.list$Ntratamiento[,2])))
)
df3 <- data.frame(
Tratamiento=factor(rep(c("Exclosure", "No Exclosure"), each=length(out3$sims.list$N_Excl_Silv) )),
Abundance=c(((out3$sims.list$N_Excl_Silv)),
((out3$sims.list$N_noExcl_Silv)))
)
df4 <- data.frame(
Tratamiento=factor(rep(c("Diferencia"), each=length(out3$sims.list$Difftratamiento))),
Abundance=c((out3$sims.list$D_Silv))# ,
#((out$sims.list$Ntratamiento[,2])))
)
df5 <- data.frame(
Tratamiento=factor(rep(c("Exclosure", "No Exclosure"), each=length(out3$sims.list$mean.abundance)/2 )),
Abundance=c(((out3$sims.list$mean.abundance[,1])),
((out3$sims.list$mean.abundance[,2])))
)
df6 <- data.frame(
Tratamiento=factor(rep(c("Diferencia"), each=length(out3$sims.list$Dif.mean.abundance))),
Abundance=c((out3$sims.list$Dif.mean.abundance))# ,
#((out$sims.list$Ntratamiento[,2])))
)
library(plyr)
mu <- ddply(df, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu2 <- ddply(df2, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu3 <- ddply(df3, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu4 <- ddply(df4, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu5 <- ddply(df5, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu6 <- ddply(df6, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
library(rphylopic)
bat_pic <- image_data("18bfd2fc-f184-4c3a-b511-796aafcc70f6", size=128)[[1]] # get actual
cow_pic <- image_data("415714b4-859c-4d1c-9ce0-9e1081613df7", size=128)[[1]]
grasshoper_pic <- image_data("7b3f79f4-9d30-42c2-bb04-6b7d2a83da04", size=128)[[1]]
fly_pic <- image_data("67bb75ec-8906-4910-871c-2758415b139a", size=128)[[1]]
library(ggbreak)
# Add mean lines
p<-ggplot(df, aes(x=Abundance, color=Tratamiento, fill=Tratamiento)) +
geom_density(alpha=0.3) +
geom_histogram(alpha=0.3, bins = 50, position="identity") +
# geom_density(alpha=0.6) +
geom_vline(data=mu, aes(xintercept=grp.median, color=Tratamiento),
linetype="dashed") +
# scale_x_break(c(40000,50000), scales = 0.25) +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Posteriors of insect abundance Conventional",x="Insect abundance per trap", y = "Frecuency") +
theme(legend.position="top")
p + # add_phylopic(cow_pic, 1, 9, 200, ysize = 300) +
add_phylopic(fly_pic, 1, 9, 250, ysize = 250) 
# Add mean lines
q<-ggplot(df2, aes(x=Abundance, fill=Tratamiento)) +
geom_density( alpha=0.3) + # geom_density(alpha=0.6) +
geom_vline(data=mu2, aes(xintercept=grp.median),
linetype="dashed") +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Difference of insect abundance in Conventional - removal by bats",x="Total Insect removal", y = "Density") +
theme(legend.position="top")
q + add_phylopic(fly_pic, 1, 1, 0.1, ysize = 0.4) 
p2<-ggplot(df3, aes(x=Abundance, color=Tratamiento, fill=Tratamiento)) +
geom_density(alpha=0.3) +
geom_histogram(alpha=0.3, bins = 50, position="identity") +
# geom_density(alpha=0.6) +
geom_vline(data=mu3, aes(xintercept=grp.median, color=Tratamiento),
linetype="dashed") +
# scale_x_break(c(40000,50000), scales = 0.25) +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Posteriors of insect abundance Silvopastoral",x="Insect abundance per trap", y = "Frecuency") +
theme(legend.position="top")
p2 + add_phylopic(fly_pic, 1, 14.5, 250, ysize = 300) 
# Add mean lines
q2<-ggplot(df4, aes(x=Abundance, fill=Tratamiento)) +
geom_density( alpha=0.3) + # geom_density(alpha=0.6) +
geom_vline(data=mu4, aes(xintercept=grp.median),
linetype="dashed") +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Difference of insect abundance in Silvopastoral - removal by bats",x="Total Insect removal", y = "Density") +
theme(legend.position="top")
q2 + add_phylopic(fly_pic, 1, 6, 0.1, ysize = 0.4) 
l<-ggplot(df5, aes(x=Abundance, color=Tratamiento, fill=Tratamiento)) +
geom_density(alpha=0.3) +
geom_histogram(alpha=0.3, bins = 50, position="identity") +
# geom_density(alpha=0.6) +
geom_vline(data=mu5, aes(xintercept=grp.median, color=Tratamiento),
linetype="dashed") +
# scale_x_break(c(40000,50000), scales = 0.25) +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Posteriors of insect overall abundance",x="Mean insect per Trap", y = "Frecuency") +
theme(legend.position="top")
l
# Add mean lines
m<-ggplot(df6, aes(x=Abundance, fill=Tratamiento)) +
geom_density( alpha=0.3) + # geom_density(alpha=0.6) +
geom_vline(data=mu6, aes(xintercept=grp.median),
linetype="dashed") +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Difference of insect overall abundance - removal by bats per trap",x="Insect removal", y = "Density") +
theme(legend.position="top")
m
#
#
# hist(out$sims.list$beta0, col = "gray", main = "", xlab = "beta0", las=1)
# abline(v = data$beta0, lwd = 3, col = "red")
# hist(out$sims.list$beta1, col = "gray", main = "", xlab = "beta1", las = 1)
# abline(v = data$beta1, lwd = 3, col = "red")
#
# # Plot predicted covariate relationship with abundance
# plot(data$luna, data$N, main = "", xlab = "Covariate", ylab = "Abundance", las = 1, ylim = c(0, max(data$N)), frame.plot = FALSE)
# lines(data$luna, data$lam, type = "l", col = "red", lwd = 3)
# GLM.pred <- exp(predict(glm(apply(data$y, 1, max) ~ luna + I(luna^2), family = poisson, data = data)))
# lines(data$luna, GLM.pred, type = "l", lty = 2, col = "blue", lwd = 3)
# Nmix.pred <- exp(out$mean$alpha0 + out$mean$alpha1 * data$luna)
# points(data$luna, Nmix.pred, type = "l", col = "blue", lwd = 3)
# Grass parasites
Load Data
bat_exclusure_grass <- read_excel("D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/data/datos_analisis_servicios2022.xlsx",
sheet = "pasto")Run model
# Bundle data
y2 <- as.matrix(bat_exclusure_grass[,9:13])
win.data2 <- list(y = y2, R = nrow(y2), T = ncol(y2),
luna = as.matrix(luna[,8:12]),
Sistema = as.numeric(factor(bat_exclusure_grass$Sistema)),
Tratamiento=as.numeric(factor(bat_exclusure_grass$Tratamiento)))
# Initial values
Nst2 <- apply(y2, 1, max) + 1 # Important to give good inits for latent N
inits2 <- function() list(N = Nst2,
alpha0 = rnorm(1),
alpha1 = rnorm(1),
# alpha2 = rnorm(1),
beta0 = rnorm(2),
beta1 = rnorm(2)#,
# beta2 = runif(1, -1, 1),
# beta3 = runif(1, -1, 1)
)
# Parameters monitored
params <- c("totalN",
"alpha0", "alpha1", #"alpha2",
"beta0", "beta1", # "beta2", #"beta3",
"lam.pred",
"lam.pred.exclosure",
"N_Excl_Conv",
"N_Excl_Silv",
"N_noExcl_Conv",
"N_noExcl_Silv",
"D_Silv",
"D_Conv",
#"site.abundance",
"mean.abundance",
"Dif.mean.abundance",
"Ntratamiento",
"Difftratamiento",
"fit.actual", "fit.sim", "bpv", "c.hat")
# MCMC settings
ni <- 260000
nt <- 50
nb <- 10000
nc <- 3
# Call JAGS from R (BRT 4 min)
out3b <- jags(win.data2,
inits2,
params,
"D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/R/model_7_factors.jags",
n.chains = nc,
n.thin = nt,
n.iter = ni,
n.burnin = nb,
parallel = TRUE,
n.cores=3)##
## Processing function input.......
##
## Done.
##
## Beginning parallel processing using 3 cores. Console output will be suppressed.
##
## Parallel processing completed.
##
## Calculating statistics.......
##
## Done.# Save
saveRDS(out3b, "D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/result_models/model_7b_factors.RData")
# Summarize posteriors
print(out3b, dig = 3)## JAGS output for model 'D:/BoxFiles/Box Sync/CodigoR/Unillanos-Bat2019/R/model_7_factors.jags', generated by jagsUI.
## Estimates based on 3 chains of 260000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 10000 iterations and thin rate = 50,
## yielding 15000 total samples from the joint posterior.
## MCMC ran in parallel for 9.622 minutes at time 2022-05-21 00:03:47.
##
## mean sd 2.5% 50% 97.5% overlap0
## totalN 580.879 30.244 525.000 580.000 645.000 FALSE
## alpha0 -2.834 0.440 -3.691 -2.838 -1.971 FALSE
## alpha1 1.456 0.474 0.537 1.453 2.382 FALSE
## beta0[1] 1.758 0.232 1.289 1.762 2.199 FALSE
## beta0[2] 2.012 0.222 1.567 2.014 2.442 FALSE
## beta1[1] 0.051 0.144 -0.233 0.049 0.338 TRUE
## beta1[2] -0.057 0.140 -0.332 -0.057 0.220 TRUE
## lam.pred.exclosure[1,1] 6.140 0.682 4.876 6.114 7.533 FALSE
## lam.pred.exclosure[2,1] 7.104 0.746 5.734 7.073 8.675 FALSE
## lam.pred.exclosure[1,2] 6.140 0.682 4.876 6.114 7.533 FALSE
## lam.pred.exclosure[2,2] 7.104 0.746 5.734 7.073 8.675 FALSE
## N_Excl_Conv 5.631 0.396 4.917 5.625 6.458 FALSE
## N_Excl_Silv 7.799 0.543 6.800 7.750 8.950 FALSE
## N_noExcl_Conv 6.172 0.443 5.391 6.130 7.130 FALSE
## N_noExcl_Silv 7.037 0.501 6.143 7.000 8.095 FALSE
## D_Silv 0.762 0.572 -0.350 0.757 1.905 TRUE
## D_Conv -0.541 0.470 -1.498 -0.534 0.362 TRUE
## mean.abundance[1] 6.617 0.386 5.909 6.591 7.410 FALSE
## mean.abundance[2] 5.998 0.358 5.341 5.977 6.750 FALSE
## Dif.mean.abundance 0.619 0.354 -0.068 0.614 1.318 TRUE
## Ntratamiento[1] 291.132 16.967 260.000 290.000 326.025 FALSE
## Ntratamiento[2] 289.747 17.240 258.000 289.000 326.000 FALSE
## Difftratamiento 1.385 15.983 -30.000 1.000 33.000 TRUE
## fit.actual 680.909 27.840 630.682 679.538 739.190 FALSE
## fit.sim 311.919 21.779 271.293 311.507 356.070 FALSE
## bpv 0.000 0.000 0.000 0.000 0.000 FALSE
## c.hat 2.192 0.161 1.900 2.186 2.528 FALSE
## deviance 1676.041 17.731 1643.123 1675.521 1713.079 FALSE
## f Rhat n.eff
## totalN 1.000 1.001 2382
## alpha0 1.000 1.002 869
## alpha1 0.999 1.002 884
## beta0[1] 1.000 1.000 15000
## beta0[2] 1.000 1.000 15000
## beta1[1] 0.641 1.000 15000
## beta1[2] 0.661 1.000 15000
## lam.pred.exclosure[1,1] 1.000 1.000 12208
## lam.pred.exclosure[2,1] 1.000 1.000 8489
## lam.pred.exclosure[1,2] 1.000 1.000 12208
## lam.pred.exclosure[2,2] 1.000 1.000 8489
## N_Excl_Conv 1.000 1.000 8277
## N_Excl_Silv 1.000 1.001 2198
## N_noExcl_Conv 1.000 1.000 4383
## N_noExcl_Silv 1.000 1.000 6155
## D_Silv 0.913 1.000 11590
## D_Conv 0.876 1.000 15000
## mean.abundance[1] 1.000 1.001 2572
## mean.abundance[2] 1.000 1.000 4931
## Dif.mean.abundance 0.964 1.000 15000
## Ntratamiento[1] 1.000 1.001 2572
## Ntratamiento[2] 1.000 1.000 3645
## Difftratamiento 0.549 1.000 15000
## fit.actual 1.000 1.000 15000
## fit.sim 1.000 1.000 15000
## bpv 1.000 NaN 1
## c.hat 1.000 1.000 15000
## deviance 1.000 1.000 15000
##
## 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 = 157.2 and DIC = 1833.255
## DIC is an estimate of expected predictive error (lower is better).################ GOF test
#c.hat is estimated at only a little more than 1, and the Bayesian p-value of 0.070 almost suggests a fitting model.
ppc.plot(out3b) # Results from a posterior predictive check of GoF (based on Chi-square discrepancy) for simulated data under an N-mixture model variant
MCMCtrace(out3b,
params = c("alpha0", "alpha1", "alpha2"),#, "beta0", "beta1", "beta2"),
#ISB = FALSE,
#exact = TRUE,
pdf = FALSE)#,
# main_den= c("Inclusion probability",
# "Mean capture",
# "Variation of phi"))
MCMCtrace(out3b,
params = c("beta0", "beta1", "beta2"),
#ISB = FALSE,
#exact = TRUE,
pdf = FALSE)#,
# main_den= c("Inclusion probability",
# "Mean capture",
# "Variation of phi"))
MCMCplot(out3b,
params = c("alpha0", "alpha1", "alpha2", "beta0", "beta1", "beta2", "beta3"),
ci = c(50, 90),
ref_ovl = TRUE,
guide_lines = TRUE,
main= "Parameters")
MCMCplot(out3b,
params = c("Ntratamiento", "Difftratamiento" ),
ci = c(50, 90),
ref_ovl = TRUE,
guide_lines = TRUE,
main= "Parameters Total Population")
MCMCplot(out3b,
params = c("mean.abundance", "Dif.mean.abundance" ),
ref_ovl = TRUE,
ci = c(50, 90), #credible intercval
guide_lines = TRUE,
main= "Abundance per trap")
MCMCplot(out3b,
params = c( "N_noExcl_Conv","N_Excl_Conv", "N_noExcl_Silv", "N_Excl_Silv"),
labels=c("No Exclosure Conventional",
"Exclosure Conventional",
"No Exclosure Silvopastoral",
"Exclosure Silvopastoral"),
ci = c(50, 90),
ref_ovl = TRUE,
guide_lines = TRUE,
main= "Abundance per trap")
# # Plot posteriors
# par(mfrow = c(2, 1))
# hist(log(out$sims.list$Ntratamiento[,1]), col = "gray", main = "", xlab = "Exclosure", las = 1)
# abline(v = median(log(out$sims.list$Ntratamiento[,1])), lwd = 3, col = "red")
# hist(log(out$sims.list$Ntratamiento[,2]), col = "gray", main = "", xlab = "No Eclosure", las = 1)
# abline(v = median(log(out$sims.list$Ntratamiento[,2])), lwd = 3, col = "red")
#
#
df <- data.frame(
Tratamiento=factor(rep(c("Exclosure", "No Exclosure"), each=length(out3b$sims.list$N_Excl_Conv) )),
Abundance=c(((out3b$sims.list$N_Excl_Conv)),
((out3b$sims.list$N_noExcl_Conv)))
)
df2 <- data.frame(
Tratamiento=factor(rep(c("Diferencia"), each=length(out3b$sims.list$Difftratamiento))),
Abundance=c((out3b$sims.list$D_Conv))# ,
#((out$sims.list$Ntratamiento[,2])))
)
df3 <- data.frame(
Tratamiento=factor(rep(c("Exclosure", "No Exclosure"), each=length(out3b$sims.list$N_Excl_Silv) )),
Abundance=c(((out3b$sims.list$N_Excl_Silv)),
((out3b$sims.list$N_noExcl_Silv)))
)
df4 <- data.frame(
Tratamiento=factor(rep(c("Diferencia"), each=length(out3b$sims.list$Difftratamiento))),
Abundance=c((out3b$sims.list$D_Silv))# ,
#((out$sims.list$Ntratamiento[,2])))
)
df5 <- data.frame(
Tratamiento=factor(rep(c("Exclosure", "No Exclosure"), each=length(out3b$sims.list$mean.abundance)/2 )),
Abundance=c(((out3b$sims.list$mean.abundance[,1])),
((out3b$sims.list$mean.abundance[,2])))
)
df6 <- data.frame(
Tratamiento=factor(rep(c("Diferencia"), each=length(out3b$sims.list$Dif.mean.abundance))),
Abundance=c((out3b$sims.list$Dif.mean.abundance))# ,
#((out$sims.list$Ntratamiento[,2])))
)
library(plyr)
mu <- ddply(df, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu2 <- ddply(df2, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu3 <- ddply(df3, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu4 <- ddply(df4, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu5 <- ddply(df5, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
mu6 <- ddply(df6, "Tratamiento", summarise, grp.median=median(Abundance),
grp.mean=mean(Abundance))
library(rphylopic)
bat_pic <- image_data("18bfd2fc-f184-4c3a-b511-796aafcc70f6", size=128)[[1]] # get actual
cow_pic <- image_data("415714b4-859c-4d1c-9ce0-9e1081613df7", size=128)[[1]]
grasshoper_pic <- image_data("7b3f79f4-9d30-42c2-bb04-6b7d2a83da04", size=128)[[1]]
fly_pic <- image_data("67bb75ec-8906-4910-871c-2758415b139a", size=128)[[1]]
library(ggbreak)
# Add mean lines
p<-ggplot(df, aes(x=Abundance, color=Tratamiento, fill=Tratamiento)) +
geom_density(alpha=0.3) +
geom_histogram(alpha=0.3, bins = 50, position="identity") +
# geom_density(alpha=0.6) +
geom_vline(data=mu, aes(xintercept=grp.median, color=Tratamiento),
linetype="dashed") +
# scale_x_break(c(40000,50000), scales = 0.25) +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Posteriors of insect abundance Conventional",x="Insect abundance per trap", y = "Frecuency") +
theme(legend.position="top")
p + # add_phylopic(cow_pic, 1, 9, 200, ysize = 300) +
add_phylopic(grasshoper_pic, 1, 9.5, 50, ysize = 35) 
# Add mean lines
q<-ggplot(df2, aes(x=Abundance, fill=Tratamiento)) +
geom_density( alpha=0.3) + # geom_density(alpha=0.6) +
geom_vline(data=mu2, aes(xintercept=grp.median),
linetype="dashed") +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Difference of insect abundance in Conventional - removal by bats",x="Total Insect removal", y = "Density") +
theme(legend.position="top")
q + add_phylopic(grasshoper_pic, 1, -3, 0.2, ysize = 0.2) 
p2<-ggplot(df3, aes(x=Abundance, color=Tratamiento, fill=Tratamiento)) +
geom_density(alpha=0.3) +
geom_histogram(alpha=0.3, bins = 50, position="identity") +
# geom_density(alpha=0.6) +
geom_vline(data=mu3, aes(xintercept=grp.median, color=Tratamiento),
linetype="dashed") +
# scale_x_break(c(40000,50000), scales = 0.25) +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Posteriors of insect abundance Silvopastoral",x="Insect abundance per trap", y = "Frecuency") +
theme(legend.position="top")
p2 + add_phylopic(grasshoper_pic, 1, 8.5, 25, ysize = 20) 
# Add mean lines
q2<-ggplot(df4, aes(x=Abundance, fill=Tratamiento)) +
geom_density( alpha=0.3) + # geom_density(alpha=0.6) +
geom_vline(data=mu4, aes(xintercept=grp.median),
linetype="dashed") +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Difference of insect abundance in Silvopastoral - removal by bats",x="Total Insect removal", y = "Density") +
theme(legend.position="top")
q2 + add_phylopic(grasshoper_pic, 1, 1, 0.2, ysize = 0.2) 
l<-ggplot(df5, aes(x=Abundance, color=Tratamiento, fill=Tratamiento)) +
geom_density(alpha=0.3) +
geom_histogram(alpha=0.3, bins = 50, position="identity") +
# geom_density(alpha=0.6) +
geom_vline(data=mu5, aes(xintercept=grp.median, color=Tratamiento),
linetype="dashed") +
# scale_x_break(c(40000,50000), scales = 0.25) +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Posteriors of insect overall abundance",x="Mean insect per Trap", y = "Frecuency") +
theme(legend.position="top")
l
# Add mean lines
m<-ggplot(df6, aes(x=Abundance, fill=Tratamiento)) +
geom_density( alpha=0.3) + # geom_density(alpha=0.6) +
geom_vline(data=mu6, aes(xintercept=grp.median),
linetype="dashed") +
#scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
#scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9"))+
labs(title="Difference of insect overall abundance - removal by bats per trap",x="Insect removal", y = "Density") +
theme(legend.position="top")
m
Perhaps we should use the code in 7.6 BAYESIAN ANALYSIS IN BUGS USING THE CONDITIONAL MULTINOMIAL (for removal sampling) pg335 in AHM book?
Session Info
Details and packages used
print(sessionInfo(), locale = FALSE)## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
##
## Matrix products: default
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggbreak_0.0.9 rphylopic_0.3.0 plyr_1.8.7 bayesplot_1.9.0
## [5] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4
## [9] readr_1.4.0 tidyr_1.1.4 tibble_3.1.7 ggplot2_3.3.6
## [13] tidyverse_1.3.1 readxl_1.3.1 AHMbook_0.2.3 MCMCvis_0.15.3
## [17] jagsUI_1.5.2
##
## loaded via a namespace (and not attached):
## [1] RandomFieldsUtils_1.2.5 fs_1.5.0 lubridate_1.7.10
## [4] httr_1.4.2 tools_4.0.3 backports_1.4.1
## [7] bslib_0.2.5.1 utf8_1.2.2 R6_2.5.1
## [10] DBI_1.1.2 colorspace_2.0-3 raster_3.5-15
## [13] withr_2.5.0 sp_1.4-5 tidyselect_1.1.1
## [16] curl_4.3.2 compiler_4.0.3 RandomFields_3.3.8
## [19] cli_3.2.0 rvest_1.0.2 xml2_1.3.3
## [22] labeling_0.4.2 triebeard_0.3.0 unmarked_1.1.1
## [25] sass_0.4.0 scales_1.2.0 mvtnorm_1.1-3
## [28] ggridges_0.5.3 yulab.utils_0.0.4 digest_0.6.29
## [31] rmarkdown_2.14 pkgconfig_2.0.3 htmltools_0.5.1.1
## [34] plotrix_3.8-1 highr_0.9 dbplyr_2.1.1
## [37] rlang_1.0.2 rstudioapi_0.13 httpcode_0.3.0
## [40] farver_2.1.0 gridGraphics_0.5-1 jquerylib_0.1.4
## [43] generics_0.1.2 jsonlite_1.8.0 magrittr_2.0.3
## [46] ggplotify_0.1.0 patchwork_1.1.1 Rcpp_1.0.8.3
## [49] munsell_0.5.0 fansi_1.0.3 lifecycle_1.0.1
## [52] terra_1.5-21 stringi_1.7.6 yaml_2.3.5
## [55] MASS_7.3-54 grid_4.0.3 crayon_1.5.1
## [58] lattice_0.20-44 haven_2.4.1 hms_1.1.0
## [61] knitr_1.39 pillar_1.7.0 codetools_0.2-18
## [64] crul_1.2.0 reprex_2.0.0 glue_1.6.2
## [67] evaluate_0.15 ggfun_0.0.6 modelr_0.1.8
## [70] vctrs_0.4.1 png_0.1-7 urltools_1.7.3
## [73] cellranger_1.1.0 gtable_0.3.0 assertthat_0.2.1
## [76] xfun_0.31 gridBase_0.4-7 broom_0.7.8
## [79] coda_0.19-4 rjags_4-13 aplot_0.1.4
## [82] ellipsis_0.3.2