Pcount N-Mixture model. Abundance, One species
Diego J. Lizcano, Andrea Morales (Unillanos), GAICA
23/3/2021
#Load packages
library(readxl)
library (coda)
library (jagsUI)
library (tidyverse)
library (MCMCvis)
library(reshape2)
library(lubridate)
library(sf)
library(mapview)Load Data
Joining Unillanos and GAICA data sets
##########################################################################
# Full data set
# Unilllanos and GAICA
# 164 sitios
# 224 observed species
##########################################################################
Aves <- read_excel("D:/BoxFiles/Box Sync/CodigoR/Aves_Andrea/data/UNILLANOS_GAICA_Aves.xlsx",
sheet = "UNILLANOS_GAICA_Aves", col_types = c("text",
"text", "text", "numeric", "text",
"text", "text", "text", "text", "text",
"numeric", "text", "text", "text",
"text", "text", "text", "text", "text",
"text", "text", "text", "text", "text",
"text", "text", "text", "text", "text",
"text", "text", "text", "text", "text"))
sites <- Aves %>% distinct (decimalLatitude , decimalLongitude) %>% mutate(site=c(1:164))
Aves_site <- inner_join(Aves,sites,by="decimalLatitude")
Aves_site$date <- yday(ymd(Aves_site$eventDate))See the species
spp.melt <- reshape::melt(Aves_site, id.var=c("date" , "eventTime", "Ocasion_de_muestreo",
"Sistema", "decimalLatitude" , "decimalLongitude.x",
"scientificName", "site"), measure.var="individualCount")
spp.melt2 <- reshape::melt(Aves_site, id.var=c("date", "eventTime", "Ocasion_de_muestreo",
"Sistema", "decimalLatitude" , "decimalLongitude.x",
"scientificName", "site"), measure.var="date")
spp1 <- reshape::cast(spp.melt, site ~ Ocasion_de_muestreo ~ scientificName) # 3D arr)ay
spp2 <- reshape::cast(spp.melt, scientificName ~ site ) #2D especie por sitio
spp3 <- reshape::cast(spp.melt2, site ~ Ocasion_de_muestreo + date ~ scientificName) #2D especie por sitio
# Ojo Falta!
#obs.cov <- filter (scientificName =="Crotophaga_ani") %>% dcast(spp.melt, date ~ Ocasion_de_muestreo) # 3D array
hab <- Aves_site %>% select (site, Sistema) %>% unique()
#### sp names
spnames <- as.vector(dimnames(spp1)$scientificName)
spnames## [1] "Actitis macularius" "Amazona amazonica"
## [3] "Amazona farinosa" "Amazona mercenarius"
## [5] "Amazona ochrocephala" "Ammodramus aurifrons"
## [7] "Ammodramus humeralis" "Anthracothorax nigricollis"
## [9] "Ara severus" "Aramides cajaneus"
## [11] "Aramus guarauna" "Ardea alba"
## [13] "Ardea cocoi" "Arundinicola leucocephala"
## [15] "Atalotriccus pilaris" "Athene cunicularia"
## [17] "Brachygalba lugubris" "Brotogeris cyanoptera"
## [19] "Bubulcus ibis" "Burhinus bistriatus"
## [21] "Buteo albonotatus" "Buteogallus meridionalis"
## [23] "Buteogallus urubitinga" "Butorides striata"
## [25] "Cacicus cela" "Campephilus melanoleucos"
## [27] "Camptostoma obsoletum" "Cantorchilus leucotis"
## [29] "Capito auratus" "Caracara cheriway"
## [31] "Cardellina canadensis" "Cathartes aura"
## [33] "Cathartes burrovianus" "Catharus minimus"
## [35] "Catharus ustulatus" "Ceratopipra erythrocephala"
## [37] "Cercomacroides parkeri" "Cercomacroides tyrannina"
## [39] "Chaetura brachyura" "Chionomesa fimbriata"
## [41] "Chloroceryle amazona" "Chloroceryle americana"
## [43] "Chlorophonia cyanea" "Chlorostilbon mellisugus"
## [45] "Circus buffoni" "Coccyzus americanus"
## [47] "Coereba flaveola" "Colinus cristatus"
## [49] "Columbina minuta" "Columbina squammata"
## [51] "Columbina talpacoti" "Contopus cinereus"
## [53] "Contopus virens" "Coragyps atratus"
## [55] "Crotophaga ani" "Crotophaga major"
## [57] "Crypturellus cinereus" "Crypturellus soui"
## [59] "Crypturellus undulatus" "Cyanerpes caeruleus"
## [61] "Cyanocorax violaceus" "Cyclarhis gujanensis"
## [63] "Dacnis albiventris" "Dacnis cayana"
## [65] "Daptrius ater" "Dendrocincla fuliginosa"
## [67] "Dendrocygna autumnalis" "Dendrocygna viduata"
## [69] "Dryocopus lineatus" "Egretta caerulea"
## [71] "Elaenia chiriquensis" "Elaenia flavogaster"
## [73] "Elaenia frantzii" "Empidonomus varius"
## [75] "Eucometis penicillata" "Eudocimus ruber"
## [77] "Euphonia chlorotica" "Euphonia chrysopasta"
## [79] "Euphonia laniirostris" "Euphonia xanthogaster"
## [81] "Eupsittula pertinax" "Eurypyga helias"
## [83] "Falco femoralis" "Falco sparverius"
## [85] "Forpus conspicillatus" "Galbula tombacea"
## [87] "Gampsonyx swainsonii" "Geothlypis aequinoctialis"
## [89] "Geotrygon montana" "Geranoaetus albicaudatus"
## [91] "Geranospiza caerulescens" "Glaucis hirsutus"
## [93] "Glyphorynchus spirurus" "Gymnomystax mexicanus"
## [95] "Heliomaster longirostris" "Herpetotheres cachinnans"
## [97] "Hirundo rustica" "Hylophilus flavipes"
## [99] "Hypnelus ruficollis" "Ictinia plumbea"
## [101] "Ixothraupis guttata" "Legatus leucophaius"
## [103] "Leistes militaris" "Leptopogon amaurocephalus"
## [105] "Leptotila rufaxilla" "Leptotila verreauxi"
## [107] "Loriotus luctuosus" "Manacus manacus"
## [109] "Mecocerculus minor" "Megaceryle torquata"
## [111] "Megarynchus pitangua" "Megascops choliba"
## [113] "Melanerpes cruentatus" "Mesembrinibis cayennensis"
## [115] "Milvago chimachima" "Mimus gilvus"
## [117] "Mionectes oleagineus" "Molothrus bonariensis"
## [119] "Molothrus oryzivorus" "Momotus momota"
## [121] "Mycteria americana" "Myiarchus apicalis"
## [123] "Myiarchus cephalotes" "Myiarchus ferox"
## [125] "Myiarchus swainsoni" "Myiarchus tuberculifer"
## [127] "Myiarchus tyrannulus" "Myiozetetes cayanensis"
## [129] "Myiozetetes similis" "Myrmophylax atrothorax"
## [131] "Nyctibius griseus" "Nyctidromus albicollis"
## [133] "Ochthornis littoralis" "Opisthocomus hoazin"
## [135] "Ortalis guttata" "Orthopsittaca manilatus"
## [137] "Pachyramphus polychopterus" "Pandion haliaetus"
## [139] "Parkesia noveboracensis" "Patagioenas cayennensis"
## [141] "Patagioenas plumbea" "Patagioenas subvinacea"
## [143] "Penelope jacquacu" "Phaethornis augusti"
## [145] "Phaethornis griseogularis" "Phaethornis hispidus"
## [147] "Phimosus infuscatus" "Phyllomyias nigrocapillus"
## [149] "Piaya cayana" "Picumnus squamulatus"
## [151] "Pilherodius pileatus" "Pionus menstruus"
## [153] "Piranga rubra" "Pitangus sulphuratus"
## [155] "Poecilotriccus sylvia" "Progne chalybea"
## [157] "Progne tapera" "Psarocolius angustifrons"
## [159] "Psarocolius bifasciatus" "Psarocolius decumanus"
## [161] "Pteroglossus castanotis" "Pteroglossus inscriptus"
## [163] "Pteroglossus pluricinctus" "Ramphastos tucanus"
## [165] "Ramphastos vitellinus" "Ramphocelus carbo"
## [167] "Riparia riparia" "Rostrhamus sociabilis"
## [169] "Rupornis magnirostris" "Saltator coerulescens"
## [171] "Saltator maximus" "Schistochlamys melanopis"
## [173] "Setophaga caerulescens" "Setophaga castanea"
## [175] "Setophaga fusca" "Setophaga petechia"
## [177] "Setophaga ruticilla" "Setophaga striata"
## [179] "Sicalis flaveola" "Sicalis luteola"
## [181] "Sporophila angolensis" "Sporophila intermedia"
## [183] "Sporophila lineola" "Sporophila nigricollis"
## [185] "Stelgidopteryx ruficollis" "Stilpnia cayana"
## [187] "Streptoprocne rutila" "Streptoprocne zonaris"
## [189] "Sturnella magna" "Synallaxis albescens"
## [191] "Syrigma sibilatrix" "Tachornis squamata"
## [193] "Tangara mexicana" "Tapera naevia"
## [195] "Tersina viridis" "Thamnophilus amazonicus"
## [197] "Thamnophilus punctatus" "Thraupis episcopus"
## [199] "Thraupis palmarum" "Thryophilus rufalbus"
## [201] "Tigrisoma fasciatum" "Tityra cayana"
## [203] "Tityra inquisitor" "Todirostrum cinereum"
## [205] "Tolmomyias flaviventris" "Tolmomyias sulphurescens"
## [207] "Touit huetii" "Touit stictopterus"
## [209] "Tringa solitaria" "Troglodytes aedon"
## [211] "Trogon viridis" "Turdus ignobilis"
## [213] "Turdus leucomelas" "Turdus nudigenis"
## [215] "Tyrannulus elatus" "Tyrannus melancholicus"
## [217] "Tyrannus savana" "Tyrannus tyrannus"
## [219] "Vanellus chilensis" "Veniliornis passerinus"
## [221] "Vireo olivaceus" "Volatinia jacarina"
## [223] "Xiphorhynchus guttatus" "Zenaida auriculata"See a map
Aves_sf <- st_as_sf(Aves_site, coords = c("decimalLongitude.x", "decimalLatitude"), crs = "+proj=longlat +datum=WGS84 +no_defs")
mapview(Aves_sf, zcol="Sistema")Prepare data for 1 species (Crotophaga ani)
##### select 1 species
Vanellus_chilensis <- as.data.frame(spp1[, , 219] )
Setophaga_striata <- as.data.frame(spp1[, , 178] )
Crotophaga_ani <- as.data.frame(spp1[, , 55] )
# replace NA by 0
Vanellus_chilensis[is.na(Vanellus_chilensis)] <- 0
Setophaga_striata[is.na(Setophaga_striata)] <- 0
Crotophaga_ani[is.na(Crotophaga_ani)] <- 0
Vanellus_chilensis$site <- as.numeric(rownames(Vanellus_chilensis))
Setophaga_striata$site <- as.numeric(rownames(Setophaga_striata))
Crotophaga_ani$site <- as.numeric(rownames(Crotophaga_ani))
# add cov
Setophaga_striata %>% left_join(hab, by = c("site"))Vanellus_chilensis %>% left_join(hab, by = c("site"))Crotophaga_ani %>% left_join(hab, by = c("site"))#################################
# The species
#################################
C <- Crotophaga_ani[,1:4]
# Load unmarked, format data in unmarked data frame and summarize
library(unmarked)
umf_ss <- unmarkedFramePCount(
y=Crotophaga_ani[,1:4], # Counts matrix
siteCovs= data.frame(hab = hab$Sistema) # , # Site covariates
#obsCovs = list(day = as.data.frame(obs.cov[,2:5]))
) # Observation covs
summary(umf_ss)## unmarkedFrame Object
##
## 164 sites
## Maximum number of observations per site: 4
## Mean number of observations per site: 4
## Sites with at least one detection: 59
##
## Tabulation of y observations:
## 0 1 2 3
## 553 90 12 1
##
## Site-level covariates:
## hab
## Bosque : 18
## Pastizal: 33
## SSP :113Models
# Fit model and extract estimates
# detection and abundance, in that order
# linear model for p follows first tilde, then comes linear model for lambda
fm.Nmix0 <- pcount(~1 ~1, data=umf_ss, control=list(trace=T, REPORT=1))## initial value 337.223488
## iter 2 value 300.122602
## iter 3 value 296.063229
## iter 4 value 296.026954
## iter 5 value 296.021365
## iter 6 value 296.021287
## iter 6 value 296.021285
## final value 296.021285
## convergedfm.Nmix1 <- pcount(~1 ~hab, data=umf_ss, control=list(trace=T, REPORT=1))## initial value 337.223488
## iter 2 value 305.579882
## iter 3 value 296.796972
## iter 4 value 294.789309
## iter 5 value 290.393314
## iter 6 value 290.135479
## iter 7 value 290.052510
## iter 8 value 290.037586
## iter 9 value 290.034272
## iter 10 value 290.034173
## iter 10 value 290.034172
## final value 290.034172
## convergedfm.Nmix2 <- pcount(~1 ~hab, data=umf_ss, mixture="NB",control=list(trace=T, REPORT=1))## initial value 329.394118
## iter 2 value 306.167477
## iter 3 value 294.312293
## iter 4 value 291.597062
## iter 5 value 291.075648
## iter 6 value 290.301278
## iter 7 value 289.729213
## iter 8 value 289.686524
## iter 9 value 289.683058
## iter 10 value 289.681851
## iter 11 value 289.681467
## iter 12 value 289.681358
## iter 12 value 289.681356
## final value 289.681356
## convergedfm.Nmix3 <- pcount(~1 ~hab, data=umf_ss, mixture="ZIP",control=list(trace=T, REPORT=1))## initial value 317.764014
## iter 2 value 307.365421
## iter 3 value 303.476253
## iter 4 value 298.290989
## iter 5 value 294.049919
## iter 6 value 291.313812
## iter 7 value 290.172622
## iter 8 value 290.022193
## iter 9 value 289.997727
## iter 10 value 289.983236
## iter 11 value 289.967688
## iter 12 value 289.959130
## iter 13 value 289.958935
## iter 14 value 289.958819
## iter 15 value 289.957845
## iter 16 value 289.957522
## iter 17 value 289.957063
## iter 18 value 289.956825
## iter 19 value 289.956730
## iter 20 value 289.956719
## iter 20 value 289.956718
## iter 20 value 289.956718
## final value 289.956718
## convergedfm.Nmix4 <- pcount(~hab ~hab, data=umf_ss, control=list(trace=T, REPORT=1))## initial value 337.223488
## iter 2 value 305.544367
## iter 3 value 298.619991
## iter 4 value 295.890067
## iter 5 value 292.981290
## iter 6 value 290.648784
## iter 7 value 289.567436
## iter 8 value 289.462966
## iter 9 value 289.343541
## iter 10 value 289.211341
## iter 11 value 289.080515
## iter 12 value 288.989981
## iter 13 value 288.950374
## iter 14 value 288.930184
## iter 15 value 288.929874
## iter 16 value 288.929571
## iter 17 value 288.929545
## iter 18 value 288.929342
## iter 19 value 288.929206
## iter 20 value 288.928190
## iter 21 value 288.923230
## iter 22 value 288.917120
## iter 23 value 288.910740
## iter 24 value 288.907685
## iter 25 value 288.906649
## iter 26 value 288.906629
## iter 27 value 288.906607
## iter 28 value 288.906597
## iter 29 value 288.906587
## iter 30 value 288.906582
## iter 30 value 288.906581
## iter 30 value 288.906581
## final value 288.906581
## convergedfm.Nmix5 <- pcount(~hab ~hab, data=umf_ss, mixture="NB", control=list(trace=T, REPORT=1))## initial value 329.394118
## iter 2 value 302.615955
## iter 3 value 295.359652
## iter 4 value 295.142227
## iter 5 value 291.385183
## iter 6 value 290.695013
## iter 7 value 289.286825
## iter 8 value 289.239493
## iter 9 value 289.124922
## iter 10 value 289.022728
## iter 11 value 288.918128
## iter 12 value 288.812611
## iter 13 value 288.754587
## iter 14 value 288.727569
## iter 15 value 288.706998
## iter 16 value 288.685871
## iter 17 value 288.685137
## iter 18 value 288.684814
## iter 19 value 288.675381
## iter 20 value 288.673366
## iter 21 value 288.670211
## iter 22 value 288.669675
## iter 23 value 288.669274
## iter 24 value 288.669217
## iter 25 value 288.669126
## iter 26 value 288.668746
## iter 27 value 288.667971
## iter 28 value 288.666150
## iter 29 value 288.662877
## iter 30 value 288.662607
## iter 31 value 288.660485
## iter 32 value 288.660324
## iter 33 value 288.660320
## iter 34 value 288.660302
## iter 35 value 288.660240
## iter 36 value 288.660217
## iter 37 value 288.660169
## iter 38 value 288.660162
## iter 39 value 288.660101
## iter 40 value 288.659989
## iter 41 value 288.659687
## iter 42 value 288.659196
## iter 43 value 288.658891
## iter 44 value 288.658679
## iter 45 value 288.658577
## iter 45 value 288.658577
## iter 45 value 288.658577
## final value 288.658577
## convergedfm.Nmix6 <- pcount(~hab ~hab, data=umf_ss, mixture="ZIP",control=list(trace=T, REPORT=1))## initial value 317.764014
## iter 2 value 304.251754
## iter 3 value 301.077475
## iter 4 value 298.462984
## iter 5 value 294.541673
## iter 6 value 291.081920
## iter 7 value 290.103110
## iter 8 value 289.584594
## iter 9 value 289.422891
## iter 10 value 289.347624
## iter 11 value 289.165179
## iter 12 value 289.068411
## iter 13 value 288.968790
## iter 14 value 288.933861
## iter 15 value 288.909547
## iter 16 value 288.900457
## iter 17 value 288.899199
## iter 18 value 288.898903
## iter 19 value 288.898751
## iter 20 value 288.898120
## iter 21 value 288.896834
## iter 22 value 288.894770
## iter 23 value 288.892013
## iter 24 value 288.888964
## iter 25 value 288.885141
## iter 26 value 288.884601
## iter 27 value 288.884399
## iter 28 value 288.884287
## iter 29 value 288.883282
## iter 30 value 288.881788
## iter 31 value 288.878025
## iter 32 value 288.877033
## iter 33 value 288.876874
## iter 34 value 288.876493
## iter 35 value 288.876393
## iter 36 value 288.876232
## iter 37 value 288.876113
## iter 38 value 288.875766
## iter 39 value 288.875608
## iter 40 value 288.875521
## iter 41 value 288.875477
## iter 42 value 288.875373
## iter 43 value 288.875131
## iter 44 value 288.874569
## iter 45 value 288.873465
## iter 46 value 288.871680
## iter 47 value 288.871458
## iter 48 value 288.871210
## iter 49 value 288.871192
## iter 50 value 288.871089
## iter 51 value 288.870829
## iter 52 value 288.870473
## iter 53 value 288.869759
## iter 54 value 288.869605
## iter 55 value 288.869570
## iter 55 value 288.869566
## iter 55 value 288.869566
## final value 288.869566
## convergedrbind(AIC.null=fm.Nmix0@AIC,
AIC.P_hab=fm.Nmix1@AIC, #best
AIC.NB_hab=fm.Nmix2@AIC,
AIC.ZIP_hab=fm.Nmix3@AIC,
AIC.P_hab_hab=fm.Nmix4@AIC,
AIC.NB_hab_hab=fm.Nmix5@AIC,
AIC.ZIP_hab_hab=fm.Nmix6@AIC
) # fm.Nmix2 best## [,1]
## AIC.null 596.0426
## AIC.P_hab 588.0683
## AIC.NB_hab 589.3627
## AIC.ZIP_hab 589.9134
## AIC.P_hab_hab 589.8132
## AIC.NB_hab_hab 591.3172
## AIC.ZIP_hab_hab 591.7391# Predictions of lambda for specified values of vegHt, say 0.2 and 2.1
newdat <- data.frame(hab=hab$Sistema) #uniques also
predict(fm.Nmix1, type="state", newdata=newdat, append = T)Empirical Bayes methods can underestimate the variance of the posterior distribution because they do not account for uncertainty in the hyperparameters (lambda or psi). Simulation studies indicate that the posterior mode can exhibit (3-5 percent) negatively bias as a point estimator of site-specific abundance.
Assess goodness-of-fit
pb <- parboot(fm.Nmix1, nsim=250, report=10) # goodness of fit## t0 = 123.1331plot (pb) # plot goodness of fit
Ok!… model pass!
predicting using best model
################## Per site
# Estimates of conditional abundance distribution at each site
re <- ranef(fm.Nmix1)
# Best Unbiased Predictors
bup(re, stat="mean") # Posterior mean## [1] 0.22110866 0.04245567 0.11110099 0.11110099 0.11110099 1.04245567
## [7] 1.21109279 0.04245567 0.22110866 0.11110099 0.22110866 1.11110099
## [13] 0.11110099 1.11110099 1.11110099 0.11110099 1.40218072 0.22110866
## [19] 0.04245567 0.11110099 0.22110866 0.22110866 0.11110099 0.22110866
## [25] 0.22110866 0.04245567 0.04245567 1.40218072 1.37713657 0.11110099
## [31] 0.11110099 0.11110099 0.04245567 1.40218072 1.40218072 0.22110866
## [37] 0.22110866 0.04245567 0.11110099 0.11110099 0.11110099 0.11110099
## [43] 0.11110099 0.11110099 0.11110099 0.22110866 0.11110099 0.11110099
## [49] 0.22110866 0.04245567 0.04245567 0.11110099 0.22110866 1.22110866
## [55] 0.22110866 0.22110866 0.22110866 0.11110099 0.22110866 1.11110099
## [61] 2.33971563 0.22110866 1.22110866 0.22110866 0.04245567 0.22110866
## [67] 0.22110866 0.04245567 1.22110866 2.22110866 0.04245567 1.22110866
## [73] 1.22110866 0.22110866 0.04245567 1.22110866 1.04245567 0.22110866
## [79] 1.66562297 0.22110866 1.22110866 1.22110866 0.22110866 1.40218072
## [85] 1.66562297 1.22110866 0.22110866 1.22110866 0.11110099 0.22110866
## [91] 0.11110099 0.22110866 0.11110099 1.40218072 1.22110866 0.22110866
## [97] 1.40218072 1.11110099 0.22110866 1.22110866 1.40218072 1.22110866
## [103] 0.22110866 1.22110866 0.22110866 0.22110866 0.04245567 0.04245567
## [109] 1.22110866 1.66562297 0.22110866 0.22110866 0.22110866 0.04245567
## [115] 0.22110866 0.22110866 0.22110866 0.22110866 0.22110866 0.22110866
## [121] 1.22110866 1.66562297 0.22110866 1.22110866 0.22110866 0.22110866
## [127] 0.22110866 0.22110866 0.22110866 0.22110866 0.22110866 0.22110866
## [133] 1.40218072 3.71808239 1.22110866 1.40218072 0.22110866 2.45456238
## [139] 1.66562297 1.40218072 1.22110866 2.62527772 0.22110866 0.22110866
## [145] 1.22110866 0.22110866 3.30182330 1.40218072 1.22110866 2.32065744
## [151] 1.40218072 0.22110866 0.22110866 0.22110866 0.22110866 0.22110866
## [157] 1.22110866 2.62527772 0.22110866 0.22110866 1.22110866 0.22110866
## [163] 2.32065744 0.22110866bup(re, stat="mode") # Posterior mode## [1] 0 0 0 0 0 1 1 0 0 0 0 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0 0
## [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 2 0 1 0 0 0 0 0 1 2 0 1 1 0
## [75] 0 1 1 0 1 0 1 1 0 1 1 1 0 1 0 0 0 0 0 1 1 0 1 1 0 1 1 1 0 1 0 0 0 0 1 1 0
## [112] 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 1 3 1 1 0 2 1 1 1 2 0 0 1 0 3 1
## [149] 1 2 1 0 0 0 0 0 1 2 0 0 1 0 2 0confint(re, level=0.9) # 90% CI## 5% 95%
## [1,] 0 1
## [2,] 0 0
## [3,] 0 1
## [4,] 0 1
## [5,] 0 1
## [6,] 1 1
## [7,] 1 2
## [8,] 0 0
## [9,] 0 1
## [10,] 0 1
## [11,] 0 1
## [12,] 1 2
## [13,] 0 1
## [14,] 1 2
## [15,] 1 2
## [16,] 0 1
## [17,] 1 3
## [18,] 0 1
## [19,] 0 0
## [20,] 0 1
## [21,] 0 1
## [22,] 0 1
## [23,] 0 1
## [24,] 0 1
## [25,] 0 1
## [26,] 0 0
## [27,] 0 0
## [28,] 1 3
## [29,] 1 2
## [30,] 0 1
## [31,] 0 1
## [32,] 0 1
## [33,] 0 0
## [34,] 1 3
## [35,] 1 3
## [36,] 0 1
## [37,] 0 1
## [38,] 0 0
## [39,] 0 1
## [40,] 0 1
## [41,] 0 1
## [42,] 0 1
## [43,] 0 1
## [44,] 0 1
## [45,] 0 1
## [46,] 0 1
## [47,] 0 1
## [48,] 0 1
## [49,] 0 1
## [50,] 0 0
## [51,] 0 0
## [52,] 0 1
## [53,] 0 1
## [54,] 1 2
## [55,] 0 1
## [56,] 0 1
## [57,] 0 1
## [58,] 0 1
## [59,] 0 1
## [60,] 1 2
## [61,] 2 3
## [62,] 0 1
## [63,] 1 2
## [64,] 0 1
## [65,] 0 0
## [66,] 0 1
## [67,] 0 1
## [68,] 0 0
## [69,] 1 2
## [70,] 2 3
## [71,] 0 0
## [72,] 1 2
## [73,] 1 2
## [74,] 0 1
## [75,] 0 0
## [76,] 1 2
## [77,] 1 1
## [78,] 0 1
## [79,] 1 3
## [80,] 0 1
## [81,] 1 2
## [82,] 1 2
## [83,] 0 1
## [84,] 1 3
## [85,] 1 3
## [86,] 1 2
## [87,] 0 1
## [88,] 1 2
## [89,] 0 1
## [90,] 0 1
## [91,] 0 1
## [92,] 0 1
## [93,] 0 1
## [94,] 1 3
## [95,] 1 2
## [96,] 0 1
## [97,] 1 3
## [98,] 1 2
## [99,] 0 1
## [100,] 1 2
## [101,] 1 3
## [102,] 1 2
## [103,] 0 1
## [104,] 1 2
## [105,] 0 1
## [106,] 0 1
## [107,] 0 0
## [108,] 0 0
## [109,] 1 2
## [110,] 1 3
## [111,] 0 1
## [112,] 0 1
## [113,] 0 1
## [114,] 0 0
## [115,] 0 1
## [116,] 0 1
## [117,] 0 1
## [118,] 0 1
## [119,] 0 1
## [120,] 0 1
## [121,] 1 2
## [122,] 1 3
## [123,] 0 1
## [124,] 1 2
## [125,] 0 1
## [126,] 0 1
## [127,] 0 1
## [128,] 0 1
## [129,] 0 1
## [130,] 0 1
## [131,] 0 1
## [132,] 0 1
## [133,] 1 3
## [134,] 3 5
## [135,] 1 2
## [136,] 1 3
## [137,] 0 1
## [138,] 2 4
## [139,] 1 3
## [140,] 1 3
## [141,] 1 2
## [142,] 2 4
## [143,] 0 1
## [144,] 0 1
## [145,] 1 2
## [146,] 0 1
## [147,] 2 5
## [148,] 1 3
## [149,] 1 2
## [150,] 2 3
## [151,] 1 3
## [152,] 0 1
## [153,] 0 1
## [154,] 0 1
## [155,] 0 1
## [156,] 0 1
## [157,] 1 2
## [158,] 2 4
## [159,] 0 1
## [160,] 0 1
## [161,] 1 2
## [162,] 0 1
## [163,] 2 3
## [164,] 0 1# TOTAL Pop
sum(bup(re))## [1] 106.4055sistema.func <- function(x){ #x will be length 20 since M=20
#Mean of first 10 sites
SSP <- mean(x[which(hab$Sistema=="SSP")])
#Mean of sites 11-20
Bosque <- mean(x[which(hab$Sistema=="Bosque")])
Pastizal <- mean(x[which(hab$Sistema=="Pastizal")])
#Naming elements of the output is optional but helpful
return(c(SSP=SSP, Bosque=Bosque, Pastizal=Pastizal))
}
#Get 100 samples of the values calculated in your function
pr <- predict(re, func=sistema.func, nsims=500)
#Summarize posteriors abundances in sistemas
data.frame(Abundance=rowSums(pr), #rowMeans()
se=apply(pr, 1, stats::sd),
lower=apply(pr, 1, stats::quantile, 0.025),
upper=apply(pr, 1, stats::quantile, 0.975))