Multi-Species Occupancy model, Agroforestería para la Conservación
Diego J. Lizcano, Angélica Benítez-G, Angelica Diaz-Pulido, Maria Fernanda Ordonez-Castro, et al.
5/3/2021
Dorazio-Royle Community model with covariates
This is a Multi-Species occupancy model (Dorazio-Royle) with covariates and data augmentation.
Aunque el modelo permite el modelando de la abundancia y el número de especies, la respuesta que estamos buscando esta enfocada en el hiperparametro de la ocupación para toda la comunidad de mamíferos.
# rad packages
library (coda)
library (jagsUI)
library (tidyverse)
library (readr)
library (MCMCvis)
library (ggplot2)Read and prepare data
Estamos usando 21 especies 58 sitios.
Las covariables estan organizadas con la etiqueta de la camara.
# Modified from Chapter 11 Marc Kéry & J. Andy Royle
# Applied hierarchical modeling in ecology
#
# 11.7.2 Dorazio-Royle community model with covariates
# ------------------------------------------------------------------------
# Augment data set: choose one of two different priors on Ntotal
nz <- 5 # Use for vague prior on Ntotal: M = 395
#nz <- 215 - nspec # Use for informative prior on Ntotal: M = 215
nsp = 21
nsite= 58
#setwd(choose.dir()) # Choose folder where data was downloaded
#### SPP datA
#spp <- read.csv("D:/BoxFiles/Box Sync/CodigoR/Biodiv_Caqueta/data/S2.csv", header=TRUE,sep=",",na.strings=TRUE) # species data
spp1 <- read.csv("D:/BoxFiles/Box Sync/CodigoR/Biodiv_Caqueta/data/21spp_by12.csv", header=TRUE,sep=",") # species data
#spp.melt=reshape::melt(spp,id.var=c("species", "site", "period"), measure.var="occ")
#spp1=reshape::cast(spp.melt, site ~ period ~ species) # 3D array
invec <- gather(spp1[,2:253])
mat.names <- c("Cabassous",
"Cuniculus",
"Dasyprocta",
"Dasypus_k",
"Dasypus",
"Eira",
"Leopardus_p",
"Nasua",
"Panthera",
"Pecari",
"Puma_c",
"Puma_y",
"Leopardus_w",
"Procyon",
"Tamandua",
"Galictis",
"Didelphis",
"Tremarctos",
"Myoprocta",
"Philander_a",
"Philander_o",
"Augment_sp1",
"Augment_sp2",
"Augment_sp3",
"Augment_sp4",
"Augment_sp5")
datarray <- array(invec$value,dim = c(58,12,21),
dimnames = list(site=spp1$X,
rep=c(1:12),
sps=mat.names[1:21])) # exclude Augment_spp
yaug <- array(0, dim=c(58, 12, nsp+nz)) # array with only zeroes
yaug[,,1:nsp] <- datarray #y # copy into it the observed data
# Create same NA pattern in augmented species as in the observed species
missings <- is.na(yaug[,,1]) # e.g., third survey in high-dist_defvation quads
for(k in (nsp+1):(nsp+nz)){
yaug[,,k][missings] <- NA
}
##################
# import covs
##################
# elev = dist_def,
# forest = ndvi,
# DAT = MOON,
# DUR = RAINFALL
covs.original <- read_csv("D:/BoxFiles/Box Sync/CodigoR/Biodiv_Caqueta/data/cam_and_covs.csv")
cam_and_covs <- read_csv("D:/BoxFiles/Box Sync/CodigoR/Biodiv_Caqueta/data/cam_and_covs2.csv")
dist_def <- cam_and_covs$dist_def
roughness <- cam_and_covs$roughness#$ndvi
MOON <- as.matrix(cam_and_covs[26:37])
RAINFALL <- as.matrix(as.data.frame(scale(cam_and_covs[14:25])))
### MOON
# center: 0.5255917
# scale: 0.3251623
cam_and_covsBundle data and specify the model
##
## model {
## # Priors
## omega ~ dunif(0,1) # Data augmentation or 'occupancy' parameter
## # Priors for species-specific effects in occupancy and detection
## for(k in 1:M){ # loop over sp
## lpsi[k] ~ dnorm(mu.lpsi, tau.lpsi) # Hyperparams describe community, occupancy probability
## betalpsi1[k] ~ dnorm(mu.betalpsi1, tau.betalpsi1)
## betalpsi2[k] ~ dnorm(mu.betalpsi2, tau.betalpsi2)
## betalpsi3[k] ~ dnorm(mu.betalpsi3, tau.betalpsi3)
## lp[k] ~ dnorm(mu.lp, tau.lp) # Hyperparams describe community, detection probability
## betalp1[k] ~ dnorm(mu.betalp1, tau.betalp1)
## betalp2[k] ~ dnorm(mu.betalp2, tau.betalp2)
## betalp3[k] ~ dnorm(mu.betalp3, tau.betalp3)
## }
##
## # Hyperpriors
## # For the model of occupancy
## mu.lpsi ~ dnorm(0,0.01) # Community mean of occupancy (logit)
## tau.lpsi <- pow(sd.lpsi, -2)
## sd.lpsi ~ dunif(0,8) # heterogeneity in the logit transform of occupancy probability # as always, bounds of uniform chosen by trial and error
## mu.betalpsi1 ~ dnorm(0,0.1)
## tau.betalpsi1 <- pow(sd.betalpsi1, -2)
## sd.betalpsi1 ~ dunif(0, 4)
## mu.betalpsi2 ~ dnorm(0,0.1)
## tau.betalpsi2 <- pow(sd.betalpsi2, -2)
## sd.betalpsi2 ~ dunif(0,2)
## mu.betalpsi3 ~ dnorm(0,0.1)
## tau.betalpsi3 <- pow(sd.betalpsi3, -2)
## sd.betalpsi3 ~ dunif(0,2)
##
## # For the model of detection
## mu.lp ~ dnorm(0,0.1)
## tau.lp <- pow(sd.lp, -2)
## sd.lp ~ dunif(0, 2) # heterogeneity in the logit transform of detection probability
## mu.betalp1 ~ dnorm(0,0.1)
## tau.betalp1 <- pow(sd.betalp1, -2)
## sd.betalp1 ~ dunif(0,1)
## mu.betalp2 ~ dnorm(0,0.1)
## tau.betalp2 <- pow(sd.betalp2, -2)
## sd.betalp2 ~ dunif(0,1)
## mu.betalp3 ~ dnorm(0,0.1)
## tau.betalp3 <- pow(sd.betalp3, -2)
## sd.betalp3 ~ dunif(0,1)
##
## # Superpopulation process: Ntotal species sampled out of M available
## for(k in 1:M){
## w[k] ~ dbern(omega)
## }
##
## # Ecological model for true occurrence (process model)
## for(k in 1:M){
## for (i in 1:nsite) {
## logit(psi[i,k]) <- lpsi[k] + betalpsi1[k] * dist_def[i] +
## betalpsi2[k] * pow(dist_def[i],2) + betalpsi3[k] * roughness[i]
## mu.psi[i,k] <- w[k] * psi[i,k]
## z[i,k] ~ dbern(mu.psi[i,k])
## }
## }
##
## # Observation model for replicated detection/nondetection observations
## for(k in 1:M){
## for (i in 1:nsite){
## for(j in 1:nrep){
## logit(p[i,j,k]) <- lp[k] + betalp1[k] * MOON[i,j] +
## betalp2[k] * pow(MOON[i,j],2) + betalp3[k] * RAINFALL[i,j]
## mu.p[i,j,k] <- z[i,k] * p[i,j,k]
## y[i,j,k] ~ dbern(mu.p[i,j,k])
## }
## }
## }
##
## # Derived quantities
## for(k in 1:M){
## Nocc.fs[k] <- sum(z[,k]) # Number of occupied sites among the 21sp
## mean.psi.fs[k] <- sum(z[,k])/58
## }
## for (i in 1:nsite){
## Nsite[i] <- sum(z[i,]) # Number of occurring species at each site
## # lpsi[i]
## }
## n0 <- sum(w[(nsp+1):(nsp+nz)]) # Number of unseen species
## Ntotal <- sum(w[]) # Total metacommunity size
##
## # Vectors to save (S for <U+FFFD>save<U+FFFD>; discard posterior samples for
## # all minus 1 of the potential species to save disk space)
## # we do this for nz = 250 (i.e., M = 395)
## lpsiS[1:(nsp+1)] <- lpsi[1:(nsp+1)]
## betalpsi1S[1:(nsp+1)] <- betalpsi1[1:(nsp+1)]
## betalpsi2S[1:(nsp+1)] <- betalpsi2[1:(nsp+1)]
## betalpsi3S[1:(nsp+1)] <- betalpsi3[1:(nsp+1)]
## lpS[1:(nsp+1)] <- lp[1:(nsp+1)]
## betalp1S[1:(nsp+1)] <- betalp1[1:(nsp+1)]
## betalp2S[1:(nsp+1)] <- betalp2[1:(nsp+1)]
## betalp3S[1:(nsp+1)] <- betalp3[1:(nsp+1)]
## }Initial values MCMC settings and run JAGS model
# Initial values
wst <- rep(1, nsp+nz) # Simply set everybody at occurring
zst <- array(1, dim = c(nsite, nsp+nz)) # ditto
inits <- function() list(z = zst, w = wst, lpsi = rnorm(n = nsp+nz), betalpsi1 = rnorm(n = nsp+nz), betalpsi2 = rnorm(n = nsp+nz), betalpsi3 = rnorm(n = nsp+nz), lp = rnorm(n = nsp+nz), betalp1 = rnorm(n = nsp+nz), betalp2 = rnorm(n = nsp+nz), betalp3 = rnorm(n = nsp+nz))
# Set 1
params1 <- c( "omega", "mu.lpsi", "sd.lpsi",
"mu.betalpsi1", "sd.betalpsi1",
"mu.betalpsi2", "sd.betalpsi2",
"mu.betalpsi3", "sd.betalpsi3",
"mu.lp", "sd.lp",
"mu.betalp1", "sd.betalp1",
"mu.betalp2", "sd.betalp2",
"mu.betalp3", "sd.betalp3",
"Ntotal", "Nsite","Nocc.fs", "w",
"mean.psi.fs",
"betalpsi1",
"betalpsi2",
"betalpsi3")
# MCMC settings
ni <- 150000#00 ;
nt <- 1000#0 ;
nb <- 50000#00 ;
nc <- 3
# Run JAGS, check convergence and summarize posteriors
out101 <- jags(win.data,
inits,
params1,
model.file="D:/BoxFiles/Box Sync/CodigoR/Biodiv_Caqueta/R/model10.txt",
n.chains = nc, n.thin = nt,
n.iter = ni, n.burnin = nb,
parallel = FALSE,
verbose=TRUE)Processing function input…….
Done.
Compiling model graph Resolving undeclared variables Allocating nodes Graph information: Observed stochastic nodes: 15158 Unobserved stochastic nodes: 4697 Total graph size: 162181
Initializing model
Adaptive phase….. Adaptive phase complete
Burn-in phase, 50000 iterations x 3 chains
Sampling from joint posterior, 1e+05 iterations x 3 chains
Calculating statistics…….
Done.
# par(mfrow = c(2, 2))
# traceplot(out101, c(c("omega", "mu.lpsi", "sd.lpsi", "mu.betalpsi1", "sd.betalpsi1", "mu.betalpsi2", "sd.betalpsi2", "mu.betalpsi3", "sd.betalpsi3", "mu.lp", "sd.lp", "mu.betalp1", "sd.betalp1", "mu.betalp2", "sd.betalp2", "mu.betalp3", "sd.betalp3", "Ntotal")) )
summary(out101)Summary for model ‘D:/BoxFiles/Box Sync/CodigoR/Biodiv_Caqueta/R/model10.txt’ Saved parameters: omega mu.lpsi sd.lpsi mu.betalpsi1 sd.betalpsi1 mu.betalpsi2 sd.betalpsi2 mu.betalpsi3 sd.betalpsi3 mu.lp sd.lp mu.betalp1 sd.betalp1 mu.betalp2 sd.betalp2 mu.betalp3 sd.betalp3 Ntotal Nsite Nocc.fs w mean.psi.fs betalpsi1 betalpsi2 betalpsi3 deviance MCMC ran for 691.224 minutes at time 2021-03-09 16:05:07.
For each of 3 chains: Adaptation: 100 iterations (sufficient) Burn-in: 50000 iterations Thin rate: 1000 iterations Total chain length: 150100 iterations Posterior sample size: 100 draws
WARNING Rhat values indicate convergence failure.
DIC info: (pD = var(deviance)/2) pD = 866 and DIC = 4542.293
Trace and density plots of MCMC chains
MCMCtrace(out101,
params = c('omega',
'mu.lpsi', "sd.lpsi",
"mu.betalpsi1", "sd.betalpsi1", "mu.betalpsi2", "sd.betalpsi2", "mu.betalpsi3", "sd.betalpsi3",
"mu.lp", "sd.lp", "mu.betalp1", "sd.betalp1", "mu.betalp2", "sd.betalp2", "mu.betalp3", "sd.betalp3",
"Ntotal"),
ind = TRUE,
pdf = FALSE, Rhat = TRUE)
Occupancy per species
# Occupancy by species
MCMCplot(out101,
params = c("mean.psi.fs"),
horiz = TRUE,
rank = TRUE,
ref_ovl = FALSE,
labels = mat.names)
Model parameters
out10 <- out101
#### worm plots
# covariate effect
MCMCplot(out101,
params = c('mu.betalpsi1', "mu.betalpsi2", "mu.betalpsi3",
"mu.betalp1", "mu.betalp2","mu.betalp3"),
horiz = TRUE,
rank = TRUE,
ref_ovl = FALSE,
#xlab = 'My x-axis label',
main = 'Posterior estimate visualization',
labels = c('Deforestation distance', 'Defor distance square', 'Roughness',
'Moon ligth', 'Moon ligth square', 'Rainfall')) 
# 'Eighth param', 'Nineth param', 'Tenth param'),
#labels_sz = 1.5, med_sz = 2, thick_sz = 7, thin_sz = 3,
#ax_sz = 4, main_text_sz = 2)
par(mfrow = c(1,2))
psi.sample <- plogis(rnorm(10^6, mean = out10$mean$mu.lpsi, sd = out10$mean$sd.lpsi))
p.sample <- plogis(rnorm(10^6, mean = out10$mean$mu.lp, sd = out10$mean$sd.lp))
hist(psi.sample, freq = F, breaks = 50, col = "grey", xlab = "Species occupancy probability", ylab = "Density", main = "")
hist(p.sample, freq = F, breaks = 50, col = "grey", xlab = "Species detection probability", ylab = "Density", main = "")
summary(psi.sample) ; summary(p.sample) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0000765 0.0979857 0.2566546 0.3305747 0.5247494 0.9992624 Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0006388 0.0654152 0.1281388 0.1718877 0.2358770 0.9726464
parecieran no ser tan buenos predictores de la comunidad…
Species response to deforestation distance (square)
# Occupancy by species
MCMCplot(out101,
params = c("betalpsi2"),
horiz = TRUE,
rank = TRUE,
ref_ovl = FALSE,
labels = mat.names)
Species response to Roughness
# Occupancy by species
MCMCplot(out101,
params = c("betalpsi3"),
horiz = TRUE,
rank = TRUE,
ref_ovl = FALSE,
labels = mat.names)
Testing the linear and cuadratic effects.
# testing linear, or cuadratic effects
par(mfrow = c(2,4)) # Among-species variability in parameters (not shown)
hist(out10$sims.list$sd.lpsi, breaks = 100, col = "grey", xlim = c(0,6), main = "Occupancy: intercept")
abline(v = mean(out10$sims.list$sd.lpsi), col = "blue", lwd = 3)
hist(out10$sims.list$sd.betalpsi1, breaks = 100, col = "grey", xlim = c(0,3), main = "Occupancy: linear effect of dist_deforestation")
abline(v = mean(out10$sims.list$sd.betalpsi1), col = "blue", lwd = 3)
hist(out10$sims.list$sd.betalpsi2, breaks = 100, col = "grey", xlim = c(0,3), main = "Occupancy: quadratic effect of dist_deforestation")
abline(v = mean(out10$sims.list$sd.betalpsi2), col = "blue", lwd = 3)
hist(out10$sims.list$sd.betalpsi3, breaks = 100, col = "grey", xlim = c(0,3), main = "Occupancy: linear effect of roughness")
abline(v = mean(out10$sims.list$sd.betalpsi3), col = "blue", lwd = 3)
hist(out10$sims.list$sd.lp, breaks = 100, col = "grey", xlim = c(0,2), main = "Detection: intercept")
abline(v = mean(out10$sims.list$sd.lp), col = "blue", lwd = 3)
hist(out10$sims.list$sd.betalp1, breaks = 100, col = "grey", xlim = c(0,1), main = "Detection: linear effect of Moon ligth")
abline(v = mean(out10$sims.list$sd.betalp1), col = "blue", lwd = 3)
hist(out10$sims.list$sd.betalp2, breaks = 100, col = "grey", xlim = c(0,1), main = "Detection: quadratic linear effect of Moon ligth")
abline(v = mean(out10$sims.list$sd.betalp2), col = "blue", lwd = 3)
hist(out10$sims.list$sd.betalp3, breaks = 100, col = "grey", xlim = c(0,1), main = "Detection: linear effect of rainfall")
abline(v = mean(out10$sims.list$sd.betalp3), col = "blue", lwd = 3)
See predictions
#####
mean.dist_def <- 1076.08
sd.dist_def <- 729.7648
mean.roug <- 72.87509
sd.roug <- 26.82622
roughness <- mean(cam_and_covs$roughness)#ndvi)
sd.roughness <- sd(cam_and_covs$roughness)#ndvi)
mean.rainfall <- mean(as.matrix(data.frame(covs.original[,14:25])))
sd.rainfall <- sd(as.matrix(data.frame(covs.original[,14:25])))
mean.moon <- mean(as.matrix(data.frame(covs.original[,26:37])))
sd.moon <- sd(as.matrix(data.frame(covs.original[,26:37])))
# Visualize covariate mean relationships for the average species
o.dist_def <- seq(90, 3900,,500) # Get covariate values for prediction
o.ndvi <- seq(-1., 1.061252,,500)
o.moon <- seq(0.01, 0.99,,500)
o.rain <- seq( 1.13, 18.88,,500)
o.roug <- seq(22.91671,134.5,,500)
dist_def.pred <- (o.dist_def - mean.dist_def) / sd.dist_def
roughness.pred <- (o.roug - mean.roug) / sd.roug #o.ndvi #(o.ndvi - mean.ndvi) / sd.ndvi
moon.pred <- (o.moon - mean.moon) / sd.moon
rain.pred <- (o.rain - mean.rainfall) / sd.rainfall
# Predict occupancy for elevation and forest and detection for date and duration
# Put all fourpredictions into a single object
str( tmp <- out10$sims.list ) # grab MCMC samplesList of 26 $ omega : num [1:300] 0.907 0.723 0.74 0.718 0.932 … $ mu.lpsi : num [1:300] -0.409 -1.702 -1.394 -0.803 -1.07 … $ sd.lpsi : num [1:300] 1.02 1.69 1.79 1.67 1.89 … $ mu.betalpsi1: num [1:300] 0.157 -0.0505 0.0748 -0.1924 0.0828 … $ sd.betalpsi1: num [1:300] 0.5884 0.445 0.3321 0.1707 0.0238 … $ mu.betalpsi2: num [1:300] -0.1153 -0.381 -0.0846 0.0114 -0.0877 … $ sd.betalpsi2: num [1:300] 0.118 1.208 0.605 0.87 0.467 … $ mu.betalpsi3: num [1:300] -0.0522 -0.0983 -0.1407 0.0923 -0.2407 … $ sd.betalpsi3: num [1:300] 0.26528 0.04947 0.00598 0.11315 0.13212 … $ mu.lp : num [1:300] -2 -1.94 -1.95 -2.06 -1.9 … $ sd.lp : num [1:300] 0.903 1.453 0.941 0.989 0.921 … $ mu.betalp1 : num [1:300] 0.01589 0.11976 0.00242 -0.00144 0.08165 … $ sd.betalp1 : num [1:300] 0.0194 0.0252 0.0539 0.0508 0.056 … $ mu.betalp2 : num [1:300] -0.1103 0.0172 -0.0345 -0.0657 -0.013 … $ sd.betalp2 : num [1:300] 0.0542 0.2671 0.1266 0.2757 0.1698 … $ mu.betalp3 : num [1:300] -0.03298 -0.01635 -0.0986 -0.07851 0.00108 … $ sd.betalp3 : num [1:300] 0.0182 0.2177 0.1941 0.1757 0.1437 … $ Ntotal : num [1:300] 21 24 23 21 23 22 23 25 23 22 … $ Nsite : num [1:300, 1:58] 14 13 12 13 13 12 12 13 12 13 … $ Nocc.fs : num [1:300, 1:26] 22 21 28 22 26 19 21 31 22 33 … $ w : num [1:300, 1:26] 1 1 1 1 1 1 1 1 1 1 … $ mean.psi.fs : num [1:300, 1:26] 0.379 0.362 0.483 0.379 0.448 … $ betalpsi1 : num [1:300, 1:26] 0.4199 -0.4938 0.344 0.0194 0.0792 … $ betalpsi2 : num [1:300, 1:26] -0.4815 -0.6851 -0.5074 -0.0777 -0.2811 … $ betalpsi3 : num [1:300, 1:26] -0.581 -0.11 -0.14 0.235 -0.145 … $ deviance : num [1:300] 3654 3741 3673 3740 3712 …
nsamp <- length(tmp[[1]]) # number of mcmc samples
predC <- array(NA, dim = c(500, nsamp, 4)) # "C" for 'community mean'
for(i in 1:nsamp){
predC[,i,1] <- plogis(tmp$mu.lpsi[i] + tmp$mu.betalpsi1[i] * dist_def.pred +
tmp$mu.betalpsi2[i] * dist_def.pred^2 )
predC[,i,2] <- plogis(tmp$mu.lpsi[i] + tmp$mu.betalpsi3[i] * roughness.pred)
predC[,i,3] <- plogis(tmp$mu.lp[i] + tmp$mu.betalp1[i] * moon.pred +
tmp$mu.betalp2[i] * moon.pred^2 )
predC[,i,4] <- plogis(tmp$mu.lp[i] + tmp$mu.betalp3[i] * rain.pred)
}
# Get posterior means and 95% CRIs and plot (Fig. 11�17)
pmC <- apply(predC, c(1,3), mean)
criC <- apply(predC, c(1,3), function(x) quantile(x, prob = c(0.025, 0.975)))
par(mfrow = c(2, 2))
plot(o.dist_def, pmC[,1], col = "blue", lwd = 3, type = 'l', lty = 1, frame = F,
ylim = c(0.1, 0.69),
xlab = "Distance to deforestation (m)", ylab = "Community mean occupancy")
matlines(o.dist_def, t(criC[,,1]), col = "grey", lty = 1)
plot(o.roug, pmC[,2], col = "blue", lwd = 3, type = 'l', lty = 1, frame = F,
ylim = c(0.1, 0.69),
xlab = "Roughness", ylab = "Community mean occupancy")
matlines(o.roug, t(criC[,,2]), col = "grey", lty = 1)
plot(o.moon, pmC[,3], col = "blue", lwd = 3, type = 'l', lty = 1, frame = F,
ylim = c(0, 0.3), xlab = "Moon ligth", ylab = "Community mean detection")
matlines(o.moon, t(criC[,,3]), col = "grey", lty = 1)
plot(o.rain, pmC[,4], col = "blue", lwd = 3, type = 'l', lty = 1, frame = F,
ylim = c(0, 0.3),
#xlim =c (1,19),
xlab = "Rainfall (mm)", ylab = "Community mean detection")
matlines(o.rain, t(criC[,,4]), col = "grey", lty = 1)
Información de la sesión en R.
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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MCMCvis_0.14.0 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.4
## [5] purrr_0.3.4 readr_1.4.0 tidyr_1.1.2 tibble_3.0.6
## [9] ggplot2_3.3.3 tidyverse_1.3.0 jagsUI_1.5.1 lattice_0.20-41
## [13] coda_0.19-4
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 xfun_0.21 bslib_0.2.4 haven_2.3.1
## [5] colorspace_2.0-0 vctrs_0.3.6 generics_0.1.0 htmltools_0.5.1.1
## [9] yaml_2.2.1 rlang_0.4.10 jquerylib_0.1.3 pillar_1.4.7
## [13] withr_2.4.1 glue_1.4.2 DBI_1.1.1 dbplyr_2.1.0
## [17] modelr_0.1.8 readxl_1.3.1 lifecycle_1.0.0 munsell_0.5.0
## [21] gtable_0.3.0 cellranger_1.1.0 rvest_0.3.6 codetools_0.2-18
## [25] evaluate_0.14 knitr_1.31 parallel_4.0.3 highr_0.8
## [29] broom_0.7.5 Rcpp_1.0.6 backports_1.2.1 scales_1.1.1
## [33] jsonlite_1.7.2 fs_1.5.0 rjags_4-10 hms_1.0.0
## [37] digest_0.6.27 stringi_1.5.3 grid_4.0.3 cli_2.3.0
## [41] tools_4.0.3 magrittr_2.0.1 sass_0.3.1 crayon_1.4.1
## [45] pkgconfig_2.0.3 ellipsis_0.3.1 xml2_1.3.2 reprex_1.0.0
## [49] lubridate_1.7.9.2 rstudioapi_0.13 assertthat_0.2.1 rmarkdown_2.7
## [53] httr_1.4.2 R6_2.5.0 compiler_4.0.3