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)
library(raster)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", "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 %>% dplyr::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")
Aves_sf %>% dplyr::select (occurrenceID, recordedBy, Sistema, locality) %>%
mapview(zcol="Sistema")load distance to forest
#### load mapa bosque
dist_forest <- raster("d://BoxFiles/Box Sync/CodigoR/Aves_Andrea/defore_raster/dist_forest_2019.tif")
plot(dist_forest)
sites <- Aves %>% distinct (Sistema, decimalLatitude , decimalLongitude) %>%
arrange(Sistema) %>% # ordena por sistema luego pone site number
mutate(site=c(1:164)) # pone site number
sites_sf <- st_as_sf(sites, coords = c("decimalLongitude", "decimalLatitude"), # convert to sf
crs = "+proj=longlat +datum=WGS84 +no_defs")
# Covs non scaled
forest.ovr <- raster::extract(dist_forest, sites_sf, method='bilinear')
forest.ovr[1:18] <- 0 # en bosque es distancia cero
fore_dist <- wiqid::standardize(forest.ovr)# standarizePrepare data for 1 species (Setophaga striata)
##### 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 cov1
Setophaga_striata <- Setophaga_striata %>% left_join(hab, by = c("site"))
Vanellus_chilensis <- Vanellus_chilensis %>% left_join(hab, by = c("site"))
Crotophaga_ani <- Crotophaga_ani %>% left_join(hab, by = c("site"))
# add cov2
Setophaga_striata$dist_forest <- fore_dist
Vanellus_chilensis$dist_forest <- fore_dist
Crotophaga_ani$dist_forest <- fore_dist
#################################
# The species
#################################
C <- Setophaga_striata[,1:4]
# Load unmarked, format data in unmarked data frame and summarize
library(unmarked)
umf_ss <- unmarkedFramePCount(
y=Setophaga_striata[,1:4], # Counts matrix
siteCovs= data.frame(hab = Setophaga_striata$Sistema,
dist_forest = Setophaga_striata$dist_forest) # , # 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: 20
##
## Tabulation of y observations:
## 0 1 2
## 618 32 6
##
## Site-level covariates:
## hab dist_forest
## Bosque : 18 Min. :-1.0283
## Pastizal: 33 1st Qu.:-0.6331
## SSP :113 Median :-0.2257
## Mean : 0.0000
## 3rd Qu.: 0.3454
## Max. : 3.7614Models
# 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), K=364)## initial value 223.264280
## iter 2 value 211.697867
## iter 3 value 207.484542
## iter 4 value 180.323448
## iter 5 value 141.932611
## iter 6 value 139.894789
## iter 7 value 138.841334
## iter 8 value 138.760121
## iter 9 value 138.666581
## iter 10 value 138.666320
## iter 10 value 138.666319
## iter 10 value 138.666319
## final value 138.666319
## convergedfm.Nmix1 <- pcount(~1 ~hab, data=umf_ss, control=list(trace=T, REPORT=1), K=364)## initial value 223.264280
## iter 2 value 141.740720
## iter 3 value 141.522996
## iter 4 value 139.265436
## iter 5 value 136.998879
## iter 6 value 134.864852
## iter 7 value 134.726776
## iter 8 value 134.669747
## iter 9 value 134.648313
## iter 10 value 134.631578
## iter 11 value 134.631499
## iter 12 value 134.631271
## iter 13 value 134.631264
## iter 14 value 134.631257
## iter 14 value 134.631257
## iter 14 value 134.631257
## final value 134.631257
## convergedfm.Nmix2 <- pcount(~1 ~hab, data=umf_ss, mixture="NB",control=list(trace=T, REPORT=1), K=364)## initial value 190.834330
## iter 2 value 172.800338
## iter 3 value 157.947224
## iter 4 value 150.801568
## iter 5 value 136.441964
## iter 6 value 132.316969
## iter 7 value 131.491846
## iter 8 value 130.912461
## iter 9 value 130.642648
## iter 10 value 130.535512
## iter 11 value 130.501067
## iter 12 value 130.500107
## iter 13 value 130.499967
## iter 14 value 130.499809
## iter 15 value 130.499772
## iter 16 value 130.499766
## iter 17 value 130.499721
## iter 17 value 130.499720
## iter 17 value 130.499720
## final value 130.499720
## convergedfm.Nmix3 <- pcount(~1 ~hab, data=umf_ss, mixture="ZIP",control=list(trace=T, REPORT=1), K=364)## initial value 154.354436
## iter 2 value 149.448852
## iter 3 value 140.065021
## iter 4 value 138.351844
## iter 5 value 135.406701
## iter 6 value 132.580396
## iter 7 value 129.935724
## iter 8 value 129.313392
## iter 9 value 129.142333
## iter 10 value 129.092294
## iter 11 value 129.085361
## iter 12 value 129.049562
## iter 13 value 129.044402
## iter 14 value 129.040816
## iter 15 value 129.037363
## iter 16 value 129.029256
## iter 17 value 129.012945
## iter 18 value 128.986971
## iter 19 value 128.963566
## iter 20 value 128.936066
## iter 21 value 128.924768
## iter 22 value 128.918582
## iter 23 value 128.911626
## iter 24 value 128.911475
## iter 25 value 128.911374
## iter 26 value 128.911355
## iter 27 value 128.911251
## iter 28 value 128.910939
## iter 29 value 128.910921
## iter 30 value 128.910887
## iter 31 value 128.910745
## iter 32 value 128.910459
## iter 33 value 128.909953
## iter 34 value 128.909836
## iter 35 value 128.909767
## iter 36 value 128.909614
## iter 37 value 128.909560
## iter 38 value 128.909408
## iter 39 value 128.909389
## iter 40 value 128.909386
## iter 41 value 128.909375
## iter 42 value 128.909373
## iter 42 value 128.909373
## iter 42 value 128.909373
## final value 128.909373
## convergedfm.Nmix4 <- pcount(~hab ~hab, data=umf_ss, control=list(trace=T, REPORT=1), K=364)## initial value 223.264280
## iter 2 value 141.415679
## iter 3 value 138.757737
## iter 4 value 137.873439
## iter 5 value 133.911423
## iter 6 value 133.291405
## iter 7 value 132.951969
## iter 8 value 132.833847
## iter 9 value 132.809195
## iter 10 value 132.801637
## iter 11 value 132.791081
## iter 12 value 132.783134
## iter 13 value 132.775396
## iter 14 value 132.769457
## iter 15 value 132.769297
## iter 16 value 132.768690
## iter 17 value 132.757244
## iter 18 value 132.752180
## iter 19 value 132.750781
## iter 20 value 132.743025
## iter 21 value 132.735024
## iter 22 value 132.717403
## iter 23 value 132.710291
## iter 24 value 132.705968
## iter 25 value 132.703992
## iter 26 value 132.702947
## iter 27 value 132.702437
## iter 27 value 132.702436
## final value 132.702436
## convergedfm.Nmix5 <- pcount(~hab ~hab, data=umf_ss, mixture="NB", control=list(trace=T, REPORT=1), K=364)## initial value 190.834330
## iter 2 value 181.889997
## iter 3 value 173.553889
## iter 4 value 151.088537
## iter 5 value 139.882837
## iter 6 value 134.835260
## iter 7 value 132.349402
## iter 8 value 130.469524
## iter 9 value 129.250612
## iter 10 value 129.041222
## iter 11 value 128.907613
## iter 12 value 128.848887
## iter 13 value 128.753139
## iter 14 value 128.701810
## iter 15 value 128.673714
## iter 16 value 128.659993
## iter 17 value 128.659639
## iter 18 value 128.659577
## iter 19 value 128.659505
## iter 20 value 128.659340
## iter 21 value 128.659165
## iter 22 value 128.658621
## iter 23 value 128.655822
## iter 24 value 128.652401
## iter 25 value 128.651471
## iter 26 value 128.651388
## iter 27 value 128.649670
## iter 28 value 128.649510
## iter 29 value 128.649466
## iter 30 value 128.649463
## iter 31 value 128.649460
## iter 31 value 128.649459
## final value 128.649459
## convergedfm.Nmix6 <- pcount(~hab ~hab, data=umf_ss, mixture="ZIP",control=list(trace=T, REPORT=1), K=364)## initial value 154.354436
## iter 2 value 153.802810
## iter 3 value 148.136897
## iter 4 value 139.597598
## iter 5 value 134.343128
## iter 6 value 133.928142
## iter 7 value 129.936346
## iter 8 value 128.511897
## iter 9 value 128.258050
## iter 10 value 128.150702
## iter 11 value 127.943866
## iter 12 value 127.844325
## iter 13 value 127.775812
## iter 14 value 127.728892
## iter 15 value 127.677332
## iter 16 value 127.628020
## iter 17 value 127.627782
## iter 18 value 127.624615
## iter 19 value 127.617987
## iter 20 value 127.617279
## iter 21 value 127.615828
## iter 22 value 127.614275
## iter 23 value 127.613746
## iter 24 value 127.612475
## iter 25 value 127.609456
## iter 26 value 127.603848
## iter 27 value 127.595021
## iter 28 value 127.585791
## iter 29 value 127.584612
## iter 30 value 127.583448
## iter 31 value 127.582558
## iter 32 value 127.582509
## iter 33 value 127.582454
## iter 34 value 127.582372
## iter 35 value 127.582343
## iter 36 value 127.582319
## iter 37 value 127.582302
## iter 38 value 127.582297
## iter 39 value 127.582287
## iter 40 value 127.582256
## iter 41 value 127.582237
## iter 42 value 127.582221
## iter 43 value 127.582215
## iter 43 value 127.582215
## iter 43 value 127.582215
## final value 127.582215
## convergedfm.Nmix7 <- pcount(~1 ~1, data=umf_ss, mixture="ZIP",control=list(trace=T, REPORT=1), K=364)## initial value 154.354436
## iter 2 value 139.799120
## iter 3 value 138.086648
## iter 4 value 134.085123
## iter 5 value 132.822679
## iter 6 value 131.998663
## iter 7 value 131.917394
## iter 8 value 131.872583
## iter 9 value 131.867167
## iter 10 value 131.854241
## iter 11 value 131.829917
## iter 12 value 131.825274
## iter 13 value 131.825213
## iter 13 value 131.825211
## iter 13 value 131.825211
## final value 131.825211
## convergedfm.Nmix8 <- pcount(~hab ~dist_forest, data=umf_ss, mixture="ZIP",control=list(trace=T, REPORT=1), K=364)## initial value 154.354436
## iter 2 value 141.536862
## iter 3 value 137.419236
## iter 4 value 135.710527
## iter 5 value 131.474175
## iter 6 value 127.250629
## iter 7 value 126.893378
## iter 8 value 126.637816
## iter 9 value 126.609562
## iter 10 value 126.603826
## iter 11 value 126.602977
## iter 12 value 126.602328
## iter 13 value 126.602301
## iter 13 value 126.602300
## iter 13 value 126.602300
## final value 126.602300
## converged#fm.Nmix9 <- pcount(~hab ~dist_forest +I(dist_forest^2), data=umf_ss, mixture="ZIP",control=list(trace=T, REPORT=1))
fm.Nmix10 <- pcount(~hab ~hab+dist_forest, data=umf_ss, mixture="ZIP",control=list(trace=T, REPORT=1), K=364)## initial value 154.354436
## iter 2 value 152.203757
## iter 3 value 141.285025
## iter 4 value 137.780366
## iter 5 value 133.203685
## iter 6 value 128.167407
## iter 7 value 126.757170
## iter 8 value 126.311499
## iter 9 value 126.128457
## iter 10 value 125.942518
## iter 11 value 125.852870
## iter 12 value 125.762629
## iter 13 value 125.729647
## iter 14 value 125.709536
## iter 15 value 125.697965
## iter 16 value 125.691941
## iter 17 value 125.688604
## iter 18 value 125.687024
## iter 19 value 125.686998
## iter 20 value 125.686988
## iter 21 value 125.686975
## iter 22 value 125.686973
## iter 23 value 125.686967
## iter 23 value 125.686965
## iter 23 value 125.686964
## final value 125.686964
## convergedrbind(AIC.null=fm.Nmix0@AIC,
AIC.P_hab=fm.Nmix1@AIC,
AIC.NB_hab=fm.Nmix2@AIC,
AIC.ZIP_hab=fm.Nmix3@AIC, # best
AIC.P_hab_hab=fm.Nmix4@AIC,
AIC.NB_hab_hab=fm.Nmix5@AIC,
AIC.ZIP_hab_hab=fm.Nmix6@AIC,
AIC.ZIP_hab_dist_forest=fm.Nmix7@AIC
) # fm.Nmix2 best## [,1]
## AIC.null 281.3326
## AIC.P_hab 277.2625
## AIC.NB_hab 270.9994
## AIC.ZIP_hab 267.8187
## AIC.P_hab_hab 277.4049
## AIC.NB_hab_hab 271.2989
## AIC.ZIP_hab_hab 269.1644
## AIC.ZIP_hab_dist_forest 269.6504#
models <- fitList(
'p(.)lambda(.)ZIP' = fm.Nmix7,
'p(hab)lambda(dist_forest)ZIP' = fm.Nmix8,
#'p(hab)lambda(dist_forest^2)ZIP' = fm.Nmix9,
'p(hab)lambda(hab+dist_forest)ZIP' = fm.Nmix10
)
modSel(models)## nPars AIC delta AICwt cumltvWt
## p(hab)lambda(dist_forest)ZIP 6 265.20 0.00 0.691 0.69
## p(hab)lambda(hab+dist_forest)ZIP 8 267.37 2.17 0.234 0.93
## p(.)lambda(.)ZIP 3 269.65 4.45 0.075 1.00LRT(fm.Nmix7, fm.Nmix8)# Predictions of lambda for specified values of vegHt, say 0.2 and 2.1
newdat <- data.frame(dist_forest=c(seq(min(fore_dist), max(fore_dist), length= 100),
seq(min(fore_dist), max(fore_dist), length= 100),
seq(min(fore_dist), max(fore_dist), length= 100)
),
hab= c(rep("SSP", 100),
rep("Bosque", 100),
rep("Pastizal", 100))
)
# data.frame(hab=hab$Sistema) #uniques also
# result <- predict(fm.Nmix8, 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.Nmix8, nsim=500, report=10) # goodness of fit ## t0 = 51.63815plot (pb) # plot goodness of fit
Ok!… model pass!
Predicting abundance in habitat types using best model and random effect
##### ##### ##### ##### ##### #####
############ Detection probability
plogis(coef(fm.Nmix8, type="det")) ## p(Int) p(habPastizal) p(habSSP)
## 0.16860618 0.08225214 0.61589764##### ##### ##### ##### ##### #####
################## Per site
# Estimates of conditional abundance distribution at each site
re <- ranef(fm.Nmix8)
# Best Unbiased Predictors
bup(re, stat="mean") # Posterior mean## [1] 1.901375782 0.062888079 0.150323950 0.150323950 0.150323950 1.508554300
## [7] 0.150323950 1.286099261 1.194171606 0.150323950 2.282665872 0.150323950
## [13] 1.557184095 0.150323950 0.150323950 0.150323950 1.356771023 1.356771023
## [19] 2.218685890 0.051849842 1.048682994 0.015548343 0.103432944 2.161825056
## [25] 1.145276261 0.036411420 1.196257040 0.025699798 0.055680040 0.072355925
## [31] 0.059955166 0.073451966 1.145548121 2.160556045 2.222255622 1.449865638
## [37] 2.201200005 0.038260346 0.096512141 0.112850208 0.101912690 0.095139852
## [43] 0.098704009 0.105281823 0.110512176 1.385808320 0.081063219 0.116091469
## [49] 0.033107353 0.038794360 0.035557799 0.081145675 0.010416228 0.015420324
## [55] 0.012717079 0.029351555 0.015177425 0.053942085 0.032331674 0.058277247
## [61] 0.049392990 0.021693249 0.019556981 0.022194915 0.043251146 0.033344070
## [67] 0.029111801 0.045632705 0.023557485 0.030354449 0.046413759 0.028947611
## [73] 0.026760102 0.033630387 0.039511665 0.037955597 0.053374586 0.029650361
## [79] 0.027407061 0.024843554 0.033442547 0.036636042 0.037811018 0.037103158
## [85] 0.036753363 0.039139025 0.033966399 0.029780721 0.089848830 0.023007070
## [91] 0.069801942 0.019681958 0.087316153 0.032779097 0.028703044 0.033649589
## [97] 0.029395331 0.071139644 0.019404395 0.018833904 0.019401563 0.023996772
## [103] 1.091718025 0.019335383 0.019021206 0.022211193 0.032114824 0.038263711
## [109] 0.030693250 0.024040537 0.031234509 0.033583185 0.034672902 0.057839685
## [115] 0.037104597 0.025093211 0.033668941 0.033008050 0.036360744 0.033902690
## [121] 0.027857750 0.031769161 0.036552631 0.035145163 0.031016262 0.031561269
## [127] 0.031627008 0.026165966 0.029706077 0.033852247 0.037226218 0.025695514
## [133] 0.029605461 0.029446739 0.028508672 0.029451913 0.025704982 0.021371603
## [139] 0.018811371 0.019692092 0.022369902 0.023432848 0.002832504 0.002559653
## [145] 0.003415543 0.004195030 0.004290567 0.005198051 0.005195632 0.005199800
## [151] 0.004588663 0.004086839 0.003040955 0.018509904 0.021484142 0.023500629
## [157] 0.020037632 0.017184847 0.014610246 0.013173960 0.014252056 0.015610709
## [163] 0.016293879 0.018019500bup(re, stat="mode") # Posterior mode## [1] 2 0 0 0 0 1 0 1 1 0 2 0 1 0 0 0 1 1 2 0 1 0 0 2 1 0 1 0 0 0 0 0 1 2 2 1 2
## [38] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [149] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0# confint(re, level=0.9) # 90% CI
# TOTAL # Estimated population size
sum(bup(re))## [1] 38.0782sistema.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=164)
#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))Predict with distance to forest
library(AICcmodavg)
pred_lambda<-predict(fm.Nmix8, type="state",newdata=newdat, appendData=TRUE) # or state
pred_lambda <- pred_lambda[1:100,]
pred1 <- predictSE(fm.Nmix8, #parm.type="lambda",
newdata=newdat[1:100,], c.hat = 2.47)
pred2 <- modavgPred(cand.set = list(fm.Nmix8), newdata=newdat[1:100,], parm.type =
"lambda", type = "response", c.hat = 2.47)
# 1-SE bounds
bounds <- as.data.frame(cbind(pred1[[1]]-pred1[[2]], pred1[[1]]+pred1[[2]]))
names(bounds) <- c("inf", "sup")
pred_lambda <- cbind(pred_lambda, bounds)
##########################################
########### plot Lambda
##########################################
# png(file = "C:/Users/diego.lizcano/Box Sync/CodigoR/Toshiba/Pacoche_unmakerd/fig/p_L_wiedii.png", width = 1200, height = 1200, res = 300)
par(mar = c(5, 4, 2, 2) + 0.1)
plot(Predicted ~ dist_forest, pred_lambda, type="l", ylim=c(-0.1,1), col="blue",
xlab="distance to forest",
ylab="Expected abundance",
xaxt="n")#
lines(inf ~ dist_forest, pred_lambda, type="l", col=gray(0.5))
lines(sup ~ dist_forest, pred_lambda, type="l", col=gray(0.5))
#matlines(pred_lambda$Predicted, cbind(pred1[[1]]-pred1[[2]], pred1[[1]]+pred1[[2]]), type = "l",
#lty = 1, lwd = 1, col = "gray")
xticks <- -2:3.9 # see min(cam.and.covs.scaled$elev) # see wich number equals to var name using: colnames(cam.and.covs)
xlabs <- xticks*sd(forest.ovr) + mean(forest.ovr) #Use the mean and sd of original value to change label name
axis(1, at=xticks, labels=round(xlabs, 1))
##turn off graphics device
#dev.off( )