Bat insect exclosures Nmix
Francisco Sanchez, Diego J. Lizcano, et al.
24/12/2019
#Libraries used (Packages)
library(ggplot2) # Used for data visualisation (i.e barcharts)
library(reshape2) # Manipulation of data frame
library(RColorBrewer) # Colours in bar charts
library(knitr) # Complies this dynamic report
library (readxl) # read excel data
library (tidyverse) # Data handling
library (plyr)
library (hrbrthemes)
# library(dplyr)
# library(tidyr)
library (viridis)
library (unmarked) # occu analitics
# library (ggcorrplot) # corplot in ggplot
library (raster)
library (rasterVis)
library (kableExtra) # nice tables
library (stringr)
library (lubridate)
library (ggeffects) # test and plot glm
library (bayesplot)Introduction
This document gives an introduction to the statistical analysis of insect data gathered by exclosure experiments.
include_url("https://www.gaiagps.com/public/4UcmM8lScsTnQxxG8OyZ8MPV/?slideshow=true&layer=GaiaTopoRasterFeet", height = "600px")The data was analized by an N Mixture model.
Fitting the multinomial-Poisson mixture model to data collected using survey methods such as removal sampling.
Load data
bat_exclusure <- read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Unillanos-Bat2019/data/Datos_Bat_exclusion_12dic2019_tydy_fix_unmarked_no_na.xlsx",
sheet = "Full_plagas")
luna <- read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Unillanos-Bat2019/data/Datos_Bat_exclusion_12dic2019_tydy_fix_unmarked_no_na.xlsx",
sheet = "luna")Fit the N-mixture model of Royle (2004) in unmarked.
Multinomial-Poisson mixture model of Royle (2004) fit by multinomPois.
# Organize data
CountMat <- as.matrix(bat_exclusure[,7:11])
umf <- unmarkedFrameMPois(y=CountMat, siteCovs=data.frame(bat_exclusure[,1:5])
, obsCovs=list(luna=as.matrix(luna[,8:12])),
type = "removal"
)
summary(umf)## unmarkedFrame Object
##
## 88 sites
## Maximum number of observations per site: 5
## Mean number of observations per site: 5
## Sites with at least one detection: 88
##
## Tabulation of y observations:
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 30
## 101 59 69 54 36 23 23 23 12 13 2 6 7 1 4 3 1 2 1
##
## Site-level covariates:
## dia Finca Sistema
## Min. :2019-11-10 00:00:00 Length:88 Length:88
## 1st Qu.:2019-11-10 00:00:00 Class :character Class :character
## Median :2019-11-17 00:00:00 Mode :character Mode :character
## Mean :2019-11-16 05:27:16
## 3rd Qu.:2019-11-18 00:00:00
## Max. :2019-11-23 00:00:00
## Tratamiento Trampa_n
## Length:88 Min. :1.000
## Class :character 1st Qu.:1.000
## Mode :character Median :2.000
## Mean :2.455
## 3rd Qu.:3.000
## Max. :5.000
##
## Observation-level covariates:
## luna
## Min. :0.0000
## 1st Qu.:0.3180
## Median :0.6075
## Mean :0.5715
## 3rd Qu.:0.9790
## Max. :1.0000# Fit a models, multinomPois order of formulas: detection, abundance
fm0 <- multinomPois(~1 ~ 1, umf, se = TRUE)
fm1 <- multinomPois(~luna ~ 1, umf, se = TRUE)
fm2 <- multinomPois(~luna + I(luna^2) ~ 1, umf, se = TRUE)
fm3 <- multinomPois(~luna + I(luna^2) ~ Sistema, umf, se = TRUE)
fm4 <- multinomPois(~luna + I(luna^2) ~ Tratamiento, umf, se = TRUE)
fm5 <- multinomPois(~luna + I(luna^2) ~ Sistema+Tratamiento, umf, se = TRUE)
fm5 <- multinomPois(~luna + I(luna^2) ~ Sistema+Tratamiento, umf, se = TRUE)
fmList <- fitList("p(.) lamdba(.)"=fm0,
"p(luna) lamdba(.)"=fm1,
"p(luna + I(luna^2)) lamdba(.)"=fm2,
"p(luna + I(luna^2)) lamdba(Sistema)"=fm3,
"p(luna + I(luna^2)) lamdba(Tratamiento)"=fm4,
"p(luna + I(luna^2)) lamdba(Sistema+Tratamiento)"=fm5
)
# Model selection
modSel(fmList)## nPars AIC delta AICwt
## p(luna + I(luna^2)) lamdba(Sistema+Tratamiento) 6 2598.78 0.00 1.0e+00
## p(luna + I(luna^2)) lamdba(Tratamiento) 5 2619.81 21.03 2.7e-05
## p(luna + I(luna^2)) lamdba(Sistema) 5 2628.66 29.88 3.3e-07
## p(luna + I(luna^2)) lamdba(.) 4 2649.75 50.97 8.5e-12
## p(luna) lamdba(.) 3 2668.60 69.82 6.9e-16
## p(.) lamdba(.) 2 2668.71 69.93 6.5e-16
## cumltvWt
## p(luna + I(luna^2)) lamdba(Sistema+Tratamiento) 1.00
## p(luna + I(luna^2)) lamdba(Tratamiento) 1.00
## p(luna + I(luna^2)) lamdba(Sistema) 1.00
## p(luna + I(luna^2)) lamdba(.) 1.00
## p(luna) lamdba(.) 1.00
## p(.) lamdba(.) 1.00## [1] 0.0003200154 0.0003200154 0.0003200154 0.0003200154 0.0003200154
## [6] 0.0003200154 0.0003200154 0.0003200154 0.0003200154 0.0003200154
## [11] 0.0003200154 0.0003200154 0.0003200154 0.0003200154 0.0003877568
## [16] 0.0003877568 0.0003877568 0.0003877568 0.0003877568 0.0003877568
## [21] 0.0003877568 0.0003877568 0.0003877568 0.0003877568 0.0002802511
## [26] 0.0002802511 0.0002802511 0.0002802511 0.0002802511 0.0002802511
## [31] 0.0002802511 0.0002802511 0.0002802511 0.0002802511 0.0003877568
## [36] 0.0003877568 0.0003877568 0.0003877568 0.0003877568 0.0003877568
## [41] 0.0003877568 0.0003877568 0.0003877568 0.0003877568 0.0003200154
## [46] 0.0003200154 0.0003200154 0.0003200154 0.0003200154 0.0003200154
## [51] 0.0003200154 0.0003200154 0.0003200154 0.0003200154 0.0003200154
## [56] 0.0003200154 0.0003200154 0.0003200154 0.0003877568 0.0003877568
## [61] 0.0003877568 0.0003877568 0.0003877568 0.0003877568 0.0003877568
## [66] 0.0003877568 0.0003877568 0.0003877568 0.0002802511 0.0002802511
## [71] 0.0002802511 0.0002802511 0.0002802511 0.0002802511 0.0002802511
## [76] 0.0002802511 0.0002802511 0.0002802511 0.0003877568 0.0003877568
## [81] 0.0003877568 0.0003877568 0.0003877568 0.0003877568 0.0003877568
## [86] 0.0003877568 0.0003877568 0.0003877568# Empirical Bayes estimation of random effects
fm5re <- ranef(fm5, K=80)
# plot(fm5re, subset=site %in% 1:25)
sum(bup(fm5re)) # Estimated population size## [1] NaN
## lambda(Int)
## 50493.32## lambda(SistemaSilvopastoril) lambda(TratamientoSin_Encierro)
## 0.5621364 0.4271014# Estimates of random effects
re <- ranef(fm5, K=80)
#ltheme <- canonical.theme(color = FALSE)
#lattice.options(default.theme = ltheme)
plot(re, subset=site %in% 1:25)
newdat1 <- cbind(luna=as.numeric(seq(0,1, length=100)), Sistema="Silvopastoril", Tratamiento="Encierro")
newdat2 <- cbind(luna=as.numeric(seq(0,1, length=100)), Sistema="Silvopastoril", Tratamiento="Sin_Encierro")
newdat3 <- cbind(luna=as.numeric(seq(0,1, length=100)), Sistema="Convencional", Tratamiento="Encierro")
newdat4 <- cbind(luna=as.numeric(seq(0,1, length=100)), Sistema="Convencional", Tratamiento="Sin_Encierro")
newdat <- as.data.frame(rbind(newdat1, newdat2, newdat3, newdat4))
newdat$luna <- rep(as.numeric(seq(0,1, length=100)),4)
pred_lambda <-predict(fm5,type="state",newdata=newdat,appendData=TRUE)
pred_lambda$Predicted <- log(pred_lambda$Predicted)
pred_p <-predict(fm5,type="det",newdata=newdat,appendData=TRUE)
pred_p$Predicted <- plogis(pred_p$Predicted)
ggplot(pred_lambda, aes(x=Sistema, y=Predicted)) +
geom_boxplot(fill='#A4A4A4', color="black") + geom_jitter()
ggplot(pred_lambda, aes(x=Tratamiento, y=Predicted)) +
geom_boxplot(fill='#A4A4A4', color="black") + geom_jitter()
umf2 <- unmarkedFrameGMM(y=CountMat, siteCovs=data.frame(bat_exclusure[,1:5])
, obsCovs=list(luna=as.matrix(luna[,8:12])),
numPrimary=1,
type = "removal")
fm0 <- gmultmix(lambdaformula = ~1, phiformula = ~1, pformula = ~1,
data=umf2)
# Fit Poisson models
fm1 <- gmultmix(~ Sistema, ~ 1, ~ 1, data = umf2)
fm2 <- gmultmix(~ Tratamiento, ~ 1, ~ 1, data = umf2)
fm3 <- gmultmix(~ Sistema + Tratamiento, ~ 1, ~ 1, data = umf2)
# Maybe p also depends on luna?
fm5 <- gmultmix(~ Sistema + Tratamiento, ~ 1, ~ luna, data = umf2)
# Fit analogous NegBin models
fm0nb <- gmultmix(~ 1, ~ 1, ~ 1, mixture = "NB", data = umf2)
fm1nb <- gmultmix(~ Sistema, ~ 1, ~ 1, mixture = "NB", data = umf2)
fm2nb <- gmultmix(~ Tratamiento, ~ 1, ~ 1, mixture = "NB", data = umf2)
fm3nb <- gmultmix(~ Sistema + Tratamiento , ~ 1, ~ 1, mixture = "NB", data = umf2)
# maybe p also depends on luna?
fm5nb <- gmultmix(~ Sistema + Tratamiento, ~ 1, ~ luna, mixture = "NB",
data = umf2)
# Rank models by AIC
gms <- fitList(
"lam(.)p(.)" = fm0,
"lam(Sistema)p(.)" = fm1,
"lam(Tratamiento)p(.)" = fm2,
"lam(Sistema+Tratamiento)p(.)" = fm3,
"lam(Sistema+Tratamiento)p(luna)" = fm5,
"NB,lam(.)p(.)" = fm0nb,
"NB,lam(Sistema)p(.)" = fm1nb,
"NB,lam(Tratamiento)p(.)" = fm2nb,
"NB,lam(Sistema+Tratamiento)p(.)" = fm3nb,
"NB,lam(Sistema+Tratamiento)p(luna)" = fm5nb
)
gms1 <- modSel(gms)
gms1## nPars AIC delta AICwt cumltvWt
## NB,lam(Sistema+Tratamiento)p(luna) 6 -765.33 0.00 6.8e-01 0.68
## NB,lam(Sistema+Tratamiento)p(.) 5 -762.88 2.45 2.0e-01 0.87
## NB,lam(Tratamiento)p(.) 4 -761.42 3.91 9.6e-02 0.97
## NB,lam(Sistema)p(.) 4 -758.42 6.91 2.1e-02 0.99
## NB,lam(.)p(.) 3 -756.93 8.40 1.0e-02 1.00
## lam(Sistema+Tratamiento)p(.) 4 -665.33 100.00 1.3e-22 1.00
## lam(Sistema+Tratamiento)p(luna) 5 -663.41 101.92 5.0e-23 1.00
## lam(Tratamiento)p(.) 3 -652.36 112.97 2.0e-25 1.00
## lam(Sistema)p(.) 3 -640.24 125.09 4.6e-28 1.00
## lam(.)p(.) 2 -626.85 138.48 5.7e-31 1.00## [1] 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724
## [8] 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724
## [15] 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550
## [22] 0.2849550 0.2849550 0.2849550 0.3199744 0.3199744 0.3199744 0.3199744
## [29] 0.3199744 0.3199744 0.3199744 0.3199744 0.3199744 0.3199744 0.2849550
## [36] 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550
## [43] 0.2849550 0.2849550 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724
## [50] 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724 0.2565724
## [57] 0.2565724 0.2565724 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550
## [64] 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550 0.3199744 0.3199744
## [71] 0.3199744 0.3199744 0.3199744 0.3199744 0.3199744 0.3199744 0.3199744
## [78] 0.3199744 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550 0.2849550
## [85] 0.2849550 0.2849550 0.2849550 0.2849550# Empirical Bayes estimation of random effects
fm5nbre <- ranef(fm5nb, K=80)
# plot(fm5re, subset=site %in% 1:25)
sum(bup(fm5nbre)) # Estimated population size## [1] 5206.467
## lambda(Int)
## 61.04965## lambda(SistemaSilvopastoril) lambda(TratamientoSin_Encierro)
## 0.5459331 0.4364225# Estimates of random effects
re <- ranef(fm5nb, K=80)
#ltheme <- canonical.theme(color = FALSE)
#lattice.options(default.theme = ltheme)
plot(re, subset=site %in% 1:25)
pred_lambda_nb <-predict(fm5nb,type="lambda",newdata=newdat,appendData=TRUE)
pred_lambda_nb$Predicted <- log(pred_lambda$Predicted)
pred_p_nb <-predict(fm5nb,type="det",newdata=newdat,appendData=TRUE)
#pred_p_nb$Predicted <- plogis(pred_p$Predicted)
#pred_p_nb$lower <- plogis(pred_p$lower)
#pred_p_nb$upper <- plogis(pred_p$upper)
ggplot(pred_lambda_nb, aes(x=Sistema, y=Predicted)) +
geom_boxplot(fill='#A4A4A4', color="black") + geom_jitter()
ggplot(pred_lambda_nb, aes(x=Tratamiento, y=Predicted)) +
geom_boxplot(fill='#A4A4A4', color="black") + geom_jitter()
ggplot(pred_p_nb) +
geom_line( aes(x=luna, y=Predicted), color="black") +
geom_line(aes(x=luna, y=Predicted + SE), colour="blue", linetype="dotted") +
geom_line(aes(x=luna, y=Predicted - SE), colour="blue", linetype="dotted") # + ylim(0.6, 0.65)
Session Info
Details and pakages used
## R version 3.5.3 (2019-03-11)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
##
## Matrix products: default
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] bayesplot_1.7.1 ggeffects_0.14.0 lubridate_1.7.4
## [4] kableExtra_1.1.0 rasterVis_0.46 latticeExtra_0.6-28
## [7] raster_3.0-7 sp_1.3-2 unmarked_0.13-1
## [10] Rcpp_1.0.3 lattice_0.20-38 viridis_0.5.1
## [13] viridisLite_0.3.0 hrbrthemes_0.6.0 plyr_1.8.5
## [16] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3
## [19] purrr_0.3.3 readr_1.3.1 tidyr_1.0.0.9000
## [22] tibble_2.1.3 tidyverse_1.3.0 readxl_1.3.1
## [25] knitr_1.26 RColorBrewer_1.1-2 reshape2_1.4.3
## [28] ggplot2_3.2.1
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.1 jsonlite_1.6 modelr_0.1.5 assertthat_0.2.1
## [5] highr_0.8 cellranger_1.1.0 yaml_2.2.0 gdtools_0.1.8
## [9] Rttf2pt1_1.3.7 pillar_1.4.2 backports_1.1.5 glue_1.3.1
## [13] extrafontdb_1.0 digest_0.6.21 rvest_0.3.5 colorspace_1.4-1
## [17] htmltools_0.4.0 Matrix_1.2-17 pkgconfig_2.0.3 broom_0.5.3
## [21] haven_2.2.0 webshot_0.5.1 scales_1.0.0 generics_0.0.2
## [25] sjlabelled_1.1.1 withr_2.1.2 hexbin_1.27.3 lazyeval_0.2.2
## [29] cli_1.1.0 magrittr_1.5 crayon_1.3.4 evaluate_0.14
## [33] fs_1.3.1 nlme_3.1-140 MASS_7.3-51.1 xml2_1.2.2
## [37] tools_3.5.3 hms_0.5.2 lifecycle_0.1.0 munsell_0.5.0
## [41] reprex_0.3.0 compiler_3.5.3 rlang_0.4.2 ggridges_0.5.1
## [45] grid_3.5.3 rstudioapi_0.10 labeling_0.3 rmarkdown_1.17
## [49] gtable_0.3.0 codetools_0.2-16 DBI_1.0.0 sjmisc_2.8.2
## [53] R6_2.4.0 gridExtra_2.3 zoo_1.8-5 extrafont_0.17
## [57] zeallot_0.1.0 insight_0.7.1 stringi_1.4.3 vctrs_0.2.0
## [61] dbplyr_1.4.2 tidyselect_0.2.5 xfun_0.11