Bibi

Load Packages

library(AHMbook)

## Loading required package: unmarked

## Loading required package: lattice

## Loading required package: parallel

## Loading required package: Rcpp

## Loading required package: reshape2

library(readxl)

library(rjags)

## Loading required package: coda

## Linked to JAGS 4.2.0

## Loaded modules: basemod,bugs

library(jagsUI)

## 
## Attaching package: 'jagsUI'

## The following object is masked from 'package:coda':
## 
##     traceplot

## The following object is masked from 'package:utils':
## 
##     View

library(ggeffects)
library(bayesplot)

## This is bayesplot version 1.6.0

## - Online documentation and vignettes at mc-stan.org/bayesplot

## - bayesplot theme set to bayesplot::theme_default()

##    * Does _not_ affect other ggplot2 plots

##    * See ?bayesplot_theme_set for details on theme setting

library(MCMCvis)

Load Data

habitat <- read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Landscape/habitat.xlsx")
habitat

## # A tibble: 17 x 10
##    LS_ID Tot_area2000 P_area2500 PerArea ED2500 NP2500 ENN_AVG2500
##    <dbl>        <dbl>      <dbl>   <dbl>  <dbl>  <dbl>       <dbl>
##  1   100         512.          0    26.1   29.1     10       121. 
##  2   101        1115.          0    56.8   21.9      1        60.0
##  3   102         974.          0    49.6   40.6      9        77.9
##  4   103         287.          0    14.6   17.5     10       208. 
##  5   104         234.          0    11.9   21.7     17       170. 
##  6   105         408.          0    20.8   15.7      6        88.5
##  7   106         420.          0    21.4   25.4      6        64.7
##  8   107         396.          0    20.2   35.9     15       132. 
##  9   108         609.          0    31.0   40.6     16       109. 
## 10   109        1204.          0    61.3   52.3      6        68.5
## 11   110        1254.          0    63.9   38.0      7        64.6
## 12   111         638.          0    32.5   45.7     12        95.0
## 13   112        1026.          0    52.3   27.3      3        68.3
## 14   113         780.          0    39.7   36.7      4       121. 
## 15   114         999.          0    50.9   20.9      2       111. 
## 16   115         868.          0    44.2   39.9     10        99.0
## 17   116         582.          0    29.6   25.1      7       116. 
## # ... with 3 more variables: ENN_SD2500 <dbl>, ganado <dbl>, gan1_0 <dbl>

data2 <- read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Landscape/data2.xlsx")
# str(data2)    # (Later read in as a system file of R package AHM)

Correla

cor_cov<-cor(habitat[,4:10],use="everything",method="spearman")
cor_cov

##                PerArea     ED2500       NP2500 ENN_AVG2500 ENN_SD2500
## PerArea      1.0000000  0.4583333 -0.581647883  -0.7107843 -0.6740196
## ED2500       0.4583333  1.0000000  0.268642455  -0.2598039 -0.3357843
## NP2500      -0.5816479  0.2686425  1.000000000   0.5619310  0.4633466
## ENN_AVG2500 -0.7107843 -0.2598039  0.561931006   1.0000000  0.9779412
## ENN_SD2500  -0.6740196 -0.3357843  0.463346619   0.9779412  1.0000000
## ganado       0.3441462  0.3756930 -0.007209582  -0.2724491 -0.2839207
## gan1_0       0.3015113  0.3266373  0.012632832  -0.2763854 -0.2763854
##                   ganado      gan1_0
## PerArea      0.344146243  0.30151134
## ED2500       0.375692982  0.32663729
## NP2500      -0.007209582  0.01263283
## ENN_AVG2500 -0.272449109 -0.27638540
## ENN_SD2500  -0.283920650 -0.27638540
## ganado       1.000000000  0.97019337
## gan1_0       0.970193366  1.00000000

# write.csv(cor_cov, file= "cor_cov")

library(corrplot)

## corrplot 0.84 loaded

library (ggcorrplot) # corplot in ggplot

## Loading required package: ggplot2

ggcorrplot(cor_cov, 
           hc.order = TRUE, 
           type = "lower",
           lab = TRUE)

Community model HME_Kery&Royle

data2$sp<- as.factor(data2$sp)

# Create various species lists (based on English names and systematic order)
(species.list <- levels(data2$sp))   # alphabetic list

##  [1] "Cabassous_chacoensis"      "Cerdocyon_thous"          
##  [3] "Chaetophractus_vellerosus" "Chaetophractus_villosus"  
##  [5] "Conepatus_chinga"          "Dasypus_novemcinctus"     
##  [7] "Didelphis_albiventris"     "Euphractus_sexcinctus"    
##  [9] "Galictis_cuja"             "Leopardus_geoffroyi"      
## [11] "Lepus_europaeus"           "Lycalopex_gymnocercus"    
## [13] "Mazama_gouazoubira"        "Myrmecophaga_tridactyla"  
## [15] "Pecari_tajacu"             "Puma_concolor"            
## [17] "Puma_yagouaroundi"         "Silvilagus_brasiliensis"  
## [19] "Tamandua_tetradactyla"     "Tolypeutes_matacus"

(spec.name.list <- tapply(data2$sp_id, data2$sp, mean)) # species ID

##      Cabassous_chacoensis           Cerdocyon_thous 
##                         1                         2 
## Chaetophractus_vellerosus   Chaetophractus_villosus 
##                         3                         4 
##          Conepatus_chinga      Dasypus_novemcinctus 
##                         5                         6 
##     Didelphis_albiventris     Euphractus_sexcinctus 
##                         7                         8 
##             Galictis_cuja       Leopardus_geoffroyi 
##                         9                        10 
##           Lepus_europaeus     Lycalopex_gymnocercus 
##                        11                        12 
##        Mazama_gouazoubira   Myrmecophaga_tridactyla 
##                        13                        14 
##             Pecari_tajacu             Puma_concolor 
##                        15                        16 
##         Puma_yagouaroundi   Silvilagus_brasiliensis 
##                        17                        18 
##     Tamandua_tetradactyla        Tolypeutes_matacus 
##                        19                        20

(spec.id.list <- unique(data2$sp_id))    # ID list

##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

(ordered.spec.name.list <- spec.name.list[order(spec.name.list)]) # ID-order list

##      Cabassous_chacoensis           Cerdocyon_thous 
##                         1                         2 
## Chaetophractus_vellerosus   Chaetophractus_villosus 
##                         3                         4 
##          Conepatus_chinga      Dasypus_novemcinctus 
##                         5                         6 
##     Didelphis_albiventris     Euphractus_sexcinctus 
##                         7                         8 
##             Galictis_cuja       Leopardus_geoffroyi 
##                         9                        10 
##           Lepus_europaeus     Lycalopex_gymnocercus 
##                        11                        12 
##        Mazama_gouazoubira   Myrmecophaga_tridactyla 
##                        13                        14 
##             Pecari_tajacu             Puma_concolor 
##                        15                        16 
##         Puma_yagouaroundi   Silvilagus_brasiliensis 
##                        17                        18 
##     Tamandua_tetradactyla        Tolypeutes_matacus 
##                        19                        20

#count to detection/no detection
COUNTS <- cbind(data2$'1', data2$'2', data2$'3', data2$'4', data2$'5', data2$'6', data2$'7')  # Counts camera traps

DET <- COUNTS
DET[DET > 1] <- 1       # now turned into detection/nondetection data

# Put detection data into 3D array: site x rep x species
nsite <- 17              # number of sites 
nrep <- 7                # number of replicate surveys per season
nspec <- length(species.list)   # 20 species
Y <- array(NA, dim = c(nsite, nrep, nspec))
for(i in 1:nspec){
  Y[,,i] <- DET[((i-1)*nsite+1):(i*nsite),]
}
dimnames(Y) <- list(NULL, NULL, names(ordered.spec.name.list))

# Get covariates and standardise them
names(habitat)

##  [1] "LS_ID"        "Tot_area2000" "P_area2500"   "PerArea"     
##  [5] "ED2500"       "NP2500"       "ENN_AVG2500"  "ENN_SD2500"  
##  [9] "ganado"       "gan1_0"

par(mfrow=c(1,3))    # set the plotting area into a 1*2 array
#Standardize covariates
area <- as.vector(habitat$PerArea)
m_area <- mean(area, na.rm=TRUE)
sd_area <- sd(area, na.rm=TRUE)
ls_area <- as.vector( (area-m_area) / sd_area )
boxplot(ls_area, main="ls_area")

ED <- as.vector(habitat$ED2500)
m_ED <- mean(ED, na.rm=TRUE)
sd_ED <- sd(ED, na.rm=TRUE)
ED <- as.vector( (ED-m_ED) / sd_ED )
boxplot(ED, main="ED")

np<- as.vector(habitat$NP2500)
m_np <- mean(np, na.rm=TRUE)
sd_np <- sd(np, na.rm=TRUE)
np <- as.vector( (np-m_np) / sd_np)
boxplot(np, main="np")

Dorazio-Royle community model with covariates and without Data augmentation

# Bundle and summarize data
win.data <- list(Y = Y, nsite = dim(Y)[1], nrep = dim(Y)[2], nspec = dim(Y)[3], ls_area = ls_area, ED = ED, np = np)

# # Specify model in BUGS language
# sink("C:/Users/diego.lizcano/Box Sync/CodigoR/Landscape/model10.txt")
# cat("
#     model {
#     
#     # Priors
#     # Priors for species-specific effects in occupancy and detection
#     for(k in 1:nspec){
#     lpsi[k] ~ dnorm(mu.lpsi, tau.lpsi)    # Hyperparams describe community
#     lp[k] ~ dnorm(mu.lp, tau.lp)
# 
#     betalpsi1[k] ~ dnorm(mu.betalpsi1, tau.betalpsi1)
#     betalpsi2[k] ~ dnorm(mu.betalpsi2, tau.betalpsi2)
#     betalpsi3[k] ~ dnorm(mu.betalpsi3, tau.betalpsi3)
# 
#     #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) #en el otro modelo dunif(0,1)
#     
#     mu.betalpsi1 ~ dnorm(0,0.1)
#     mu.betalpsi2 ~ dnorm(0,0.1)
#     mu.betalpsi3 ~ dnorm(0,0.1)
# 
#     tau.lpsi <- pow(sd.lpsi, -2)
#     tau.betalpsi1 <- pow(sd.betalpsi1, -2)
#     tau.betalpsi2 <- pow(sd.betalpsi2, -2)
#     tau.betalpsi3 <- pow(sd.betalpsi3, -2)
#     
#     sd.lpsi ~ dunif(0,8)   # as always, bounds of uniform chosen by trial and error
#     sd.betalpsi1 ~ dunif(0,2)
#     sd.betalpsi2 ~ dunif(0,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, 3)    # Was U(0,2) for first run of model
# 
#     #mu.betalp1 ~ dnorm(0,0.1)
#     #mu.betalp2 ~ dnorm(0,0.1)
#     #mu.betalp3 ~ dnorm(0,0.1)
# 
#     #tau.betalp1 <- pow(sd.betalp1, -2)
#     #tau.betalp2 <- pow(sd.betalp2, -2)
#     #tau.betalp3 <- pow(sd.betalp3, -2)
#     
#     #sd.betalp1 ~ dunif(0,1)
#     #sd.betalp2 ~ dunif(0,1)
#     #sd.betalp3 ~ dunif(0,1)
#     
#     # Ecological model for true occurrence (process model)
#     for(k in 1:nspec){
#     for (i in 1:nsite) {
#     logit(psi[i,k]) <- lpsi[k] + betalpsi1[k] *  ls_area[i] + betalpsi2[k] * ED[i] + betalpsi3[k] * np[i]
#     z[i,k] ~ dbern(psi[i,k])
#     }
#     }
#     
#     # Observation model for replicated detection/nondetection observations
#     for(k in 1:nspec){
#     for (i in 1:nsite){
#     for(j in 1:nrep){
#     logit(p[i,j,k]) <- lp[k] 
#     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:nspec){
#     Nocc.fs[k] <- sum(z[,k])       # Number of occupied sites among the 17
#     }
#     for (i in 1:nsite){
#     Nsite[i] <- sum(z[i,])          # Number of occurring species at each site
#     }
#     }
#     ",fill = TRUE)
# sink()


# Initial values
zst <- array(1, dim = c(nsite, nspec))
inits <- function() list(z = zst, lpsi = rnorm(n = nspec),  betalpsi1 = rnorm(n = nspec), betalpsi2 = rnorm(n = nspec), betalpsi3 = rnorm(n = nspec), lp = rnorm(n = nspec))


# Set 1
params1 <- c("mu.lpsi", "sd.lpsi", "mu.betalpsi1", "sd.betalpsi1", "mu.betalpsi2", "sd.betalpsi2", "mu.betalpsi3", "sd.betalpsi3","mu.lp", "sd.lp", "Nocc.fs", "Nsite")

# MCMC settings
# Same MCMC settings
ni <-10000   ;   nt <- 5   ;   nb <- 1000   ;   nc <- 3

# Run JAGS, check convergence and summarize posteriors: ART = 15 min

out101 <- jags(win.data, inits, params1, 
               "C:/Users/diego.lizcano/Box Sync/CodigoR/Landscape/model1.txt", 
               n.chains = nc, 
               n.thin = nt, 
               n.iter = ni, 
               n.burnin = nb)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2380
##    Unobserved stochastic nodes: 450
##    Total graph size: 4949
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 9000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

summary(out101)

## Summary for model 'C:/Users/diego.lizcano/Box Sync/CodigoR/Landscape/model1.txt'
## Saved parameters: mu.lpsi sd.lpsi mu.betalpsi1 sd.betalpsi1 mu.betalpsi2 sd.betalpsi2 mu.betalpsi3 sd.betalpsi3 mu.lp sd.lp Nocc.fs Nsite deviance 
## MCMC ran for 2.49 minutes at time 2019-07-19 20:08:13.
## 
## For each of 3 chains:
## Adaptation:            100 iterations (sufficient)
## Burn-in:               1000 iterations
## Thin rate:             5 iterations
## Total chain length:    10100 iterations
## Posterior sample size: 1800 draws
## 
## Successful convergence based on Rhat values (all < 1.1). 
## 
## DIC info: (pD = var(deviance)/2) 
## pD = 251.2 and DIC = 1647.233

par(mfrow = c(2, 2))
traceplot(out101, c(c("mu.lpsi", "sd.lpsi", "mu.betalpsi1", "sd.betalpsi1", "mu.betalpsi2", "sd.betalpsi2",  "mu.betalpsi3", "sd.betalpsi3","mu.lp", "sd.lp")) )

# Set 2
params2 <- c("mu.lpsi", "sd.lpsi", "mu.betalpsi1", "sd.betalpsi1", "mu.betalpsi2", "sd.betalpsi2","mu.betalpsi3", "sd.betalpsi3", "lpsi", "betalpsi1", "betalpsi2","betalpsi3", "lp", "z")

# Same MCMC settings
ni <- 10000   ;   nt <- 5   ;   nb <- 1000   ;   nc <- 3

# ART 15 min
out102 <- jags(win.data, inits, params2, 
               "C:/Users/diego.lizcano/Box Sync/CodigoR/Landscape/model1.txt", 
               n.chains = nc, 
               n.thin = nt, 
               n.iter = ni, 
               n.burnin = nb)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2380
##    Unobserved stochastic nodes: 450
##    Total graph size: 4949
## 
## Initializing model
## 
## Adaptive phase..... 
## Adaptive phase complete 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 9000 iterations x 3 chains 
##  
## 
## Calculating statistics.......

## Warning in doTryCatch(return(expr), name, parentenv, handler): At least one
## Rhat value could not be calculated.

## 
## Done.

summary(out102)

## Summary for model 'C:/Users/diego.lizcano/Box Sync/CodigoR/Landscape/model1.txt'
## Saved parameters: mu.lpsi sd.lpsi mu.betalpsi1 sd.betalpsi1 mu.betalpsi2 sd.betalpsi2 mu.betalpsi3 sd.betalpsi3 lpsi betalpsi1 betalpsi2 betalpsi3 lp z deviance 
## MCMC ran for 3.483 minutes at time 2019-07-19 20:10:44.
## 
## For each of 3 chains:
## Adaptation:            100 iterations (sufficient)
## Burn-in:               1000 iterations
## Thin rate:             5 iterations
## Total chain length:    10100 iterations
## Posterior sample size: 1800 draws
## 
## Successful convergence based on Rhat values (all < 1.1). 
## 
## DIC info: (pD = var(deviance)/2) 
## pD = 244.3 and DIC = 1642.289

library(coda)
all10 <- as.matrix(out102) # Put output from 3 chains into a matrix
summary(out102)            # May take a loooong time

## Summary for model 'C:/Users/diego.lizcano/Box Sync/CodigoR/Landscape/model1.txt'
## Saved parameters: mu.lpsi sd.lpsi mu.betalpsi1 sd.betalpsi1 mu.betalpsi2 sd.betalpsi2 mu.betalpsi3 sd.betalpsi3 lpsi betalpsi1 betalpsi2 betalpsi3 lp z deviance 
## MCMC ran for 3.483 minutes at time 2019-07-19 20:10:44.
## 
## For each of 3 chains:
## Adaptation:            100 iterations (sufficient)
## Burn-in:               1000 iterations
## Thin rate:             5 iterations
## Total chain length:    10100 iterations
## Posterior sample size: 1800 draws
## 
## Successful convergence based on Rhat values (all < 1.1). 
## 
## DIC info: (pD = var(deviance)/2) 
## pD = 244.3 and DIC = 1642.289

# gelman.diag(out102)        # ditto

out10 <- out101
out10

## JAGS output for model 'C:/Users/diego.lizcano/Box Sync/CodigoR/Landscape/model1.txt', generated by jagsUI.
## Estimates based on 3 chains of 10000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 1000 iterations and thin rate = 5,
## yielding 5400 total samples from the joint posterior. 
## MCMC ran for 2.49 minutes at time 2019-07-19 20:08:13.
## 
##                  mean     sd     2.5%      50%    97.5% overlap0     f
## mu.lpsi         1.076  0.455    0.252    1.043    2.080    FALSE 0.994
## sd.lpsi         0.956  0.465    0.198    0.897    2.010    FALSE 1.000
## mu.betalpsi1    0.429  0.504   -0.482    0.398    1.538     TRUE 0.811
## sd.betalpsi1    0.909  0.448    0.110    0.876    1.833    FALSE 1.000
## mu.betalpsi2    0.115  0.437   -0.697    0.100    1.028     TRUE 0.598
## sd.betalpsi2    0.680  0.456    0.033    0.615    1.730    FALSE 1.000
## mu.betalpsi3   -0.277  0.401   -1.105   -0.268    0.516     TRUE 0.769
## sd.betalpsi3    0.399  0.319    0.017    0.318    1.183    FALSE 1.000
## mu.lp          -1.723  0.391   -2.512   -1.715   -0.976    FALSE 1.000
## sd.lp           1.615  0.340    1.062    1.571    2.398    FALSE 1.000
## Nocc.fs[1]     11.131  3.201    4.000   11.000   17.000    FALSE 1.000
## Nocc.fs[2]     10.623  0.890   10.000   10.000   13.000    FALSE 1.000
## Nocc.fs[3]     10.324  3.174    4.000   10.000   16.000    FALSE 1.000
## Nocc.fs[4]      9.838  1.931    7.000   10.000   14.000    FALSE 1.000
## Nocc.fs[5]     11.881  0.966   11.000   12.000   14.000    FALSE 1.000
## Nocc.fs[6]      8.587  0.906    8.000    8.000   11.000    FALSE 1.000
## Nocc.fs[7]     14.054  1.966   10.000   14.000   17.000    FALSE 1.000
## Nocc.fs[8]     12.725  1.870    9.000   13.000   16.000    FALSE 1.000
## Nocc.fs[9]     10.541  3.481    3.000   11.000   16.000    FALSE 1.000
## Nocc.fs[10]    15.107  0.333   15.000   15.000   16.000    FALSE 1.000
## Nocc.fs[11]    10.729  3.197    4.000   11.000   17.000    FALSE 1.000
## Nocc.fs[12]     9.180  0.473    9.000    9.000   10.000    FALSE 1.000
## Nocc.fs[13]    14.321  0.577   14.000   14.000   16.000    FALSE 1.000
## Nocc.fs[14]    15.098  0.314   15.000   15.000   16.000    FALSE 1.000
## Nocc.fs[15]     9.936  3.583    2.000   10.000   16.000    FALSE 1.000
## Nocc.fs[16]     9.813  3.156    4.000   10.000   16.000    FALSE 1.000
## Nocc.fs[17]    12.493  2.415    7.000   13.000   17.000    FALSE 1.000
## Nocc.fs[18]    11.951  1.536   10.000   12.000   15.000    FALSE 1.000
## Nocc.fs[19]    10.570  3.570    2.000   11.000   16.000    FALSE 1.000
## Nocc.fs[20]    10.391  0.660   10.000   10.000   12.000    FALSE 1.000
## Nsite[1]       15.216  1.572   12.000   15.000   18.000    FALSE 1.000
## Nsite[2]       14.457  1.641   11.000   15.000   17.000    FALSE 1.000
## Nsite[3]       14.355  1.672   11.000   14.000   17.000    FALSE 1.000
## Nsite[4]       11.005  1.769    8.000   11.000   15.000    FALSE 1.000
## Nsite[5]        7.956  1.979    5.000    8.000   12.000    FALSE 1.000
## Nsite[6]       13.828  1.825   10.000   14.000   17.000    FALSE 1.000
## Nsite[7]       13.793  1.624   11.000   14.000   17.000    FALSE 1.000
## Nsite[8]        8.668  1.816    5.000    9.000   12.000    FALSE 1.000
## Nsite[9]       14.893  1.351   12.000   15.000   18.000    FALSE 1.000
## Nsite[10]      15.278  1.824   11.000   15.000   18.000    FALSE 1.000
## Nsite[11]      15.914  1.296   13.000   16.000   18.000    FALSE 1.000
## Nsite[12]      13.412  1.628   10.000   13.000   17.000    FALSE 1.000
## Nsite[13]      16.288  1.560   13.000   16.000   19.000    FALSE 1.000
## Nsite[14]      14.212  2.147   10.000   14.000   18.000    FALSE 1.000
## Nsite[15]      14.919  1.741   11.000   15.000   18.000    FALSE 1.000
## Nsite[16]      14.431  1.591   11.000   15.000   17.000    FALSE 1.000
## Nsite[17]      10.669  1.851    7.000   11.000   14.000    FALSE 1.000
## deviance     1396.015 22.425 1354.558 1395.166 1442.009    FALSE 1.000
##               Rhat n.eff
## mu.lpsi      1.004  1072
## sd.lpsi      1.020   114
## mu.betalpsi1 1.004  1894
## sd.betalpsi1 1.001  1373
## mu.betalpsi2 1.001  3294
## sd.betalpsi2 1.010   220
## mu.betalpsi3 1.000  5400
## sd.betalpsi3 1.004   470
## mu.lp        1.001  2346
## sd.lp        1.002  1068
## Nocc.fs[1]   1.004   751
## Nocc.fs[2]   1.001  5083
## Nocc.fs[3]   1.001  1882
## Nocc.fs[4]   1.000  5400
## Nocc.fs[5]   1.000  5400
## Nocc.fs[6]   1.001  4277
## Nocc.fs[7]   1.002   806
## Nocc.fs[8]   1.001  3253
## Nocc.fs[9]   1.002  4338
## Nocc.fs[10]  1.005  1231
## Nocc.fs[11]  1.001  5400
## Nocc.fs[12]  1.002  5400
## Nocc.fs[13]  1.003  1132
## Nocc.fs[14]  1.000  5400
## Nocc.fs[15]  1.004  1824
## Nocc.fs[16]  1.005   431
## Nocc.fs[17]  1.001  2535
## Nocc.fs[18]  1.002  2068
## Nocc.fs[19]  1.004   858
## Nocc.fs[20]  1.000  5400
## Nsite[1]     1.001  1614
## Nsite[2]     1.000  5400
## Nsite[3]     1.001  5400
## Nsite[4]     1.001  1425
## Nsite[5]     1.001  1874
## Nsite[6]     1.001  2240
## Nsite[7]     1.000  5400
## Nsite[8]     1.001  1816
## Nsite[9]     1.000  3285
## Nsite[10]    1.002  5400
## Nsite[11]    1.001  3314
## Nsite[12]    1.002  1298
## Nsite[13]    1.001  2831
## Nsite[14]    1.000  5400
## Nsite[15]    1.000  5400
## Nsite[16]    1.000  5400
## Nsite[17]    1.000  5400
## deviance     1.002  1552
## 
## 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 = 251.2 and DIC = 1647.233 
## DIC is an estimate of expected predictive error (lower is better).

#First, the community, which is described by the following quantities in the model: mu.lpsi, sd.lpsi, mu.betalpsi2, sd.betalpsi2, mu.betalpsi3, sd.betalpsi3, mu.lp, sd.lp, sd.betalp1, mu.betalp2, sd.betalp2, mu.betalp3, sd.betalp3

par(mfrow = c(1,2))       # Fig. 11-16
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) 

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.02448 0.60610 0.74559 0.71234 0.84807 0.99634

summary(p.sample)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000733 0.0567946 0.1521114 0.2325167 0.3472456 0.9987299

We can also look at the interspecific variability in every parameter, which appears in our model as the standard deviations of the normal priors for the species-specific effects. We simply plot histograms of the posteriors for each. Executing the next body of code you will see that the difference among the species in the Swiss breeding bird community varies a great deal among the different parameters in the occupancy and the detection model (note different scales of abscissa).

par(mfrow = c(2,2))  # 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 percentage of area in landscape")
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: linear effect edge density")
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 edge density")
abline(v = mean(out10$sims.list$sd.betalpsi3), col = "blue", lwd = 3)

Next, we look at the community average response to the modeled covariates, where we predict the mean community response of occupancy probability to elevation and forest cover and of species detection probability to survey date and survey duration.

#values for prediction
o.ls_area <- seq(0, 100,,500)
o.ed <- seq(0.1, 100,,500)
o.np <- seq(1, 100,,500)
ls_area.pred <- (o.ls_area - m_area) / sd_area
o.ed.pred <- (o.ed - m_ED) / sd_ED
o.np.pred <- (o.np - m_np) / sd_np


# Predict occupancy 
# Put all fourpredictions into a single

str( tmp <- out10$sims.list )              # grab MCMC samples

## List of 13
##  $ mu.lpsi     : num [1:5400] 1.138 0.532 1.172 0.488 1.246 ...
##  $ sd.lpsi     : num [1:5400] 0.874 1.142 0.726 0.794 0.777 ...
##  $ mu.betalpsi1: num [1:5400] 0.3861 0.4232 0.4021 0.2795 -0.0145 ...
##  $ sd.betalpsi1: num [1:5400] 1.04 0.871 0.905 0.721 0.819 ...
##  $ mu.betalpsi2: num [1:5400] -0.153 0.16 0.244 0.326 0.749 ...
##  $ sd.betalpsi2: num [1:5400] 0.456 0.713 0.589 0.342 0.355 ...
##  $ mu.betalpsi3: num [1:5400] 0.0822 -0.2324 -0.3485 -0.6031 -0.2928 ...
##  $ sd.betalpsi3: num [1:5400] 0.582 0.587 0.454 0.605 0.977 ...
##  $ mu.lp       : num [1:5400] -1.76 -2.34 -1.74 -1.75 -1.7 ...
##  $ sd.lp       : num [1:5400] 1.62 1.48 1.7 1.61 1.19 ...
##  $ Nocc.fs     : num [1:5400, 1:20] 9 14 15 13 13 13 15 13 12 12 ...
##  $ Nsite       : num [1:5400, 1:17] 18 15 14 14 13 17 16 15 16 17 ...
##  $ deviance    : num [1:5400] 1393 1373 1398 1385 1390 ...

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,2] <- plogis(tmp$mu.lpsi[i] + tmp$mu.betalpsi1[i] * ls_area.pred)
  predC[,i,3] <- plogis(tmp$mu.lpsi[i] + tmp$mu.betalpsi2[i] * o.ed.pred )
  predC[,i,4] <- plogis(tmp$mu.lpsi[i] + tmp$mu.betalpsi3[i] * o.np.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), na.rm = TRUE))

par(mfrow = c(2, 2))
plot(o.ls_area, pmC[,2], col = "blue", lwd = 3, type = 'l', lty = 1, frame = F, ylim = c(0, 1), xlab = "Percentage of forest cover", ylab = "Community mean occupancy")
matlines(o.ls_area, t(criC[,,2]), col = "grey", lty = 1)


plot(o.ed, pmC[,3], col = "blue", lwd = 3, type = 'l', lty = 1, frame = F, ylim = c(0.0, 1), xlab = "Edge density", ylab = "Community mean occupancy")
matlines(o.ed, t(criC[,,3]), col = "grey", lty = 1)

plot(o.np, pmC[,4], col = "blue", lwd = 3, type = 'l', lty = 1, frame = F, ylim = c(0.0, 1), xlab = "Number of patches", ylab = "Community mean occupancy")
matlines(o.np, t(criC[,,4]), col = "grey", lty = 1)

Second, moving to inferences at the community (i.e., site-specific) level.

#Second, moving to inferences at the community (i.e., site-specific) level, we can again obtain the estimates of species richness at each site (Nsite) and then do so for the same sample of nine sites that we already plotted in Fig. 11-14. Of course, you may plot estimates of Nsite against any site covariate in your model or even outside, to generate hypotheses about drivers of variation in community richness. When plotting the Nsite estimates against site elevation (Fig. 11???19), we see that incorporating covariates increases the Nsite estimates at lower elevations but decreases them at higher elevations, compared to the covariate-free estimates under model 9. This sensible pattern is due to the elevation relationships for most species that are estimated in the model (we will see these species-specific estimates below).

par(mfrow = c(2, 2))

ls_area <- habitat$PerArea[1:17]
plot(ls_area, out101$mean$Nsite, xlab = "Landscape area (percentage)", ylab = "Community size estimate (Nsite)", frame = F, xlim= c(10,70), ylim = c(0, 25), pch = 16) # black: model 101
lines(smooth.spline(out10$mean$Nsite ~ ls_area), lwd = 3, col = "blue")

ed <- habitat$ED2500[1:17]
plot(ed, out101$mean$Nsite, xlab = "Edge density", ylab = "Community size estimate (Nsite)", frame = F, xlim= c(10,60), ylim = c(0, 20), pch = 16) # black: model 101
lines(smooth.spline(out10$mean$Nsite ~ ed), lwd = 3, col = "blue")

np <- habitat$NP2500[1:17]
plot(np, out101$mean$Nsite, xlab = "Number of patches", ylab = "Community size estimate (Nsite)", frame = F, xlim= c(0,20),ylim = c(0, 20), pch = 16) # black: model 101
lines(smooth.spline(out10$mean$Nsite ~ np), lwd = 3, col = "blue")

Finally, and third, we present inferences about the individual species under our community occupancy model.

# all10 <- as.matrix(out102)

# str(all10)# look at the MCMC output

# names(out102)
#[1] "sims.list"   "mean"        "sd"          "q2.5"       
#[5] "q25"         "q50"         "q75"         "q97.5"      
#[9] "overlap0"    "f"           "Rhat"        "n.eff"      
#[13] "pD"          "DIC"         "summary"     "samples"    
#[17] "modfile"     "model"       "parameters"  "mcmc.info"  
#[21] "run.date"    "parallel"    "bugs.format" "calc.DIC" 

# str(all10[1])
# str(all10[2])
# print(all10[2]) #aqui estoy copiando los parametros

################################
#### aqui el error
################################

# pm <- apply(all10, 2, mean)    # Get posterior means and 95% CRIs
# str(pm)
# cri <- apply(all10[2], 2, function(x) quantile(x, prob = c(0.05, 0.95))) # CRIs
# print(pm)
# 
 

ex <- MCMCchains(out101, params = 'Nocc.fs')

cri <- apply(ex, 2, function(x) quantile(x, prob = c(0.025, 0.975))) # CRIs


MCMCplot(out102, 
         params = c("mu.betalpsi1", "mu.betalpsi2","mu.betalpsi3"), 
         horiz = TRUE,
         rank = TRUE,
         ref_ovl = FALSE)

MCMCplot(out102, 
         params = c('lp'), 
         horiz = TRUE,
         rank = TRUE,
         ref_ovl = FALSE)

MCMCplot(out102, 
         params = c('lpsi'), 
         horiz = TRUE,
         rank = TRUE,
         ref_ovl = FALSE)

MCMCplot(out101, 
         params = c('Nocc.fs'), 
         horiz = TRUE,
         rank = TRUE,
         ref_ovl = FALSE)

#Now we produce plots of the species-specific parameter estimates in the occupancy and the detection model for all 20 observed species and compare them with the community average (the hyperparameters). Beyond the interest of inspecting parameter estimates for each species, these plots should emphasize again that the community average may be quite different from what an individual species does.spec.name.list.


# Effect of patch (linear) on occupancy probability 
plot(ex[1:20], 1:20, xlim = c(-20, 20), xlab = "Parameter estimate", ylab = "Species number ID", main = "Effect of forest cover in landscape (linear) on occupancy", pch = 16)
abline(v = 0, lwd = 2, col = "black")
segments(cri[1, 1:20], 1:20, cri[2, 1:20], 1:20, col = "grey", lwd = 1)
sig4 <- (cri[1, 1:20] * cri[2, 1:20]) > 0
segments(cri[1, 1:20][sig4 == 1], (1:20)[sig4 == 1], cri[2, 1:20][sig4 == 1], (1:20)[sig4 == 1], col = "blue", lwd = 2)
abline(v = out101$summary[3,1], lwd = 3, col = "red")
abline(v = out101$summary[3,c(3,7)], lwd = 3, col = "red", lty = 2)