Venado Quindio, Risaralda
Vanessa Diaz, Diego J. Lizcano, Hugo Mantilla, et al.
22/01/2020
Cargar Paquetes y Datos
Paquetes
Cargar Datos
# Quindio_data <- read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Vanessa/data/matriz_v3.xlsx")
Quindio_data <- read_excel("D:/BoxFiles/Box Sync/CodigoR/Vanessa/data/matriz_v4.xlsx")
Quindio_data$Longitud <- Quindio_data$Longitud * (-1)
# Risaralda_data <- as.data.frame(read_excel("C:/Users/diego.lizcano/Box Sync/CodigoR/Vanessa/data/matriz_v2.xlsx", sheet = "Ucumari"))
# covaribles
Covs <- read.delim("D:/BoxFiles/Box Sync/CodigoR/Vanessa/data/Covariables.csv")
Covs_q <- data.frame(scale(Covs[,5:12])) # estandarizadas
Risaralda_data <- as.data.frame(read_excel("D:/BoxFiles/Box Sync/CodigoR/Vanessa/data/matriz_v2.xlsx", sheet = "Ucumari"))
Risaralda_data [,4] <- as.vector(as.numeric(Risaralda_data [,4]))
Risaralda_data [,5] <- as.vector(as.numeric(Risaralda_data [,5]))
Risaralda_data [,17] <- as.vector(as.numeric(Risaralda_data [,17]))
Risaralda_data$binomial <- paste(Risaralda_data$Genus, Risaralda_data$Species , sep = "_")
Risaralda_data$camera_trap_start_date <- ymd(Risaralda_data$camera_trap_start_date)
Risaralda_data$camera_trap_end_date <- ymd(Risaralda_data$camera_trap_end_date)Ver camaras
Risaralda y Quindio
Risaralda_sf <- st_as_sf(Risaralda_data, coords = c("Longitude", "Latitude"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
cams_plot <- mapview(Risaralda_sf, zcol = "Photo Type",
color = c("black","red"),
legend = FALSE,
map.types = "Esri.WorldImagery")
Quindio_sf <- st_as_sf(Quindio_data, coords = c("Longitud", "Latitud"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
# cams_plot
### Risaralda
cam_location_R <- unique(Risaralda_data[,c('Latitude','Longitude','Camera_Trap')])
#as.data.frame(unique(cbind(Risaralda_data$Longitude,
#Risaralda_data$Latitude)))
#### Quindio #### las Camaras se sobrelapan
cam_location_Q <- unique(Quindio_data[,c('Latitud','Longitud','camara')])
cams_sp_R <- st_as_sf(cam_location_R, coords = c("Longitude", "Latitude"), crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
centroid <- c(mean(Risaralda_data$Longitude), mean(Risaralda_data$Latitude))
cams_sp_Q <- as(Quindio_sf, "Spatial")
##################################
# set extent get SRTM
##################################
clip_window <- extent(-75.60 , -75.39, 4.59, 4.81)
srtm <- raster::getData('SRTM', lon=centroid[1], lat=centroid[2])
# plot(clip_window, border = "blue")
# plot(cams_sp, add=TRUE)
# crop the raster using the vector extent
srtm_crop <- crop(srtm, clip_window)
# extract altitudes
cam_location_R$Elev_R <- raster::extract(srtm_crop, cams_sp_R)
cam_location_Q$Elev_Q <- raster::extract(srtm_crop, cams_sp_Q)
cam_location_Q <- as.data.frame(cam_location_Q) # convert to Data frame
# unify names
names(cam_location_Q) <- c("Latitud", "Longitud", "camera", "Elev")
names(cam_location_R) <- c("Latitud", "Longitud", "camera", "Elev")
full_covs <- rbind(cam_location_R, cam_location_Q)
# plot(srtm_crop)
# plot(cams_sp, add=TRUE)
# view
pal = mapviewPalette("mapviewTopoColors")
Risaralda_sf %>% group_by(Camera_Trap) %>%
summarise(mean = mean("Number of Animals", na.omit=T), fotos = n()) %>%
filter(fotos > 1) %>%
mapview (cex = "fotos", zcol = "fotos", legend = TRUE) +
mapview (srtm_crop, col.regions = pal(400), at = seq(1000, 5000, 10), alpha=0.7) +
mapview(Risaralda_sf, cex = 1) + mapview(st_jitter(Quindio_sf, factor = 0.03), alpha = 0, cex = 4, zcol = "camara")Extrae matriz venado
Datos de Risaralda
# source("C:/Users/diego.lizcano/Box Sync/CodigoR/Vanessa/R/team.R")
source("D:/BoxFiles/Box Sync/CodigoR/Vanessa/R/team.R")
mat.per.sp<-f.matrix.creator2(data = Risaralda_data, year = c(2017, 2016))
sp.names<-names(mat.per.sp) # species names
# counting how many (total) records per species by all days
cont.per.sp<-data.frame(row.names = sp.names)
for (i in 1:length(mat.per.sp)){
cont.per.sp[i,1]<-sum(apply(as.data.frame(mat.per.sp [[i]]),FUN=sum,na.rm=T, MARGIN = 1))
}
cont.per.sp$especie<-rownames(cont.per.sp)
colnames(cont.per.sp)<-c("Numero_de_fotos","especie")
# nombre de la matriz de cada especie
# names(mat.per.sp)
print(as.data.frame(arrange(df = cont.per.sp, desc(Numero_de_fotos))))## Numero_de_fotos especie
## 1 53 NA_NA
## 2 16 Mazama_rufinaAnalisis de Ocupación
Crea objetos UMF
Quindio_data$Elev <- cam_location_Q$Elev_Q # pega elevacion
y_Q <- as.data.frame(Quindio_data[,4:111])
y_R <- as.data.frame(mat.per.sp[[1]])
# Covariables
#Correr el modelo de ocupacion
umf_Q<- unmarkedFrameOccu(y= y_Q, siteCovs = Covs_q) # Quindio
umf_R<- unmarkedFrameOccu(y= y_R, siteCovs = NULL) # Risaralda
#######Graficar umf
plot(umf_Q, main="Quindio")
## unmarkedFrame Object
##
## 28 sites
## Maximum number of observations per site: 108
## Mean number of observations per site: 34.25
## Sites with at least one detection: 11
##
## Tabulation of y observations:
## 0 1 <NA>
## 938 21 2065
##
## Site-level covariates:
## Helechos Frailejones Cob_dosel Cob_arb
## Min. :-0.6320 Min. :-0.6687 Min. :-0.9933 Min. :-1.9423
## 1st Qu.:-0.6320 1st Qu.:-0.6687 1st Qu.:-0.9933 1st Qu.:-0.4645
## Median :-0.5303 Median :-0.6687 Median : 0.1256 Median : 0.2252
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.2833 3rd Qu.: 0.4377 3rd Qu.: 0.8049 3rd Qu.: 0.8163
## Max. : 3.0291 Max. : 2.0548 Max. : 1.5642 Max. : 1.6045
## Cobert_her DAP Alti Temp
## Min. :-1.7738 Min. :-0.9793 Min. :-2.1097 Min. :-1.2904
## 1st Qu.:-0.8955 1st Qu.:-0.7335 1st Qu.:-0.6012 1st Qu.:-0.7277
## Median : 0.5418 Median :-0.2162 Median : 0.2183 Median :-0.1926
## Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.7814 3rd Qu.: 0.3658 3rd Qu.: 0.7265 3rd Qu.: 0.2417
## Max. : 0.7814 Max. : 3.8062 Max. : 1.4488 Max. : 2.3453
## unmarkedFrame Object
##
## 59 sites
## Maximum number of observations per site: 96
## Mean number of observations per site: 45.32
## Sites with at least one detection: 7
##
## Tabulation of y observations:
## 0 1 <NA>
## 2658 16 2990Selección de modelos
# build models first detection
m1<-occu(~1 ~ 1, umf_Q)
m2<-occu(~1 ~ Helechos, umf_Q)
m3<-occu(~1 ~ Frailejones, umf_Q)
m4<-occu(~1 ~ Cob_arb, umf_Q )
m5<-occu(~1 ~ Cobert_her, umf_Q )
m6<-occu(~1 ~ DAP, umf_Q )
m7<-occu(~1 ~ Alti, umf_Q )
m8<-occu(~1 ~ Temp, umf_Q )
m9<-occu(~1 ~ DAP + Alti, umf_Q )
m10<-occu(~1 ~ DAP + I(Alti^2), umf_Q )
# fit list
fms1<-fitList("p(.) Ocu(.)"=m1,
"p(.) Ocu(Helechos)"=m2,
"p(.) Ocu(Frailejones)"=m3,
"p(.) Ocu(Cob_arb)"=m4,
"p(.) Ocu(Cobert_her)"=m5,
"p(.) Ocu(DAP)"=m6,
"p(.) Ocu(Alti)"=m7,
"p(.) Ocu(Temp)"=m8,
"p(.) Ocu(DAP + Alt)"=m9,
"p(.) Ocu(Alti2)"=m10
)
modSel(fms1)## nPars AIC delta AICwt cumltvWt
## p(.) Ocu(DAP + Alt) 4 189.85 0.00 0.94361 0.94
## p(.) Ocu(DAP) 3 197.56 7.71 0.01999 0.96
## p(.) Ocu(Alti2) 4 198.97 9.12 0.00988 0.97
## p(.) Ocu(Alti) 3 199.03 9.18 0.00959 0.98
## p(.) Ocu(.) 2 200.08 10.22 0.00569 0.99
## p(.) Ocu(Temp) 3 200.71 10.86 0.00414 0.99
## p(.) Ocu(Cob_arb) 3 201.06 11.21 0.00348 1.00
## p(.) Ocu(Helechos) 3 201.76 11.91 0.00245 1.00
## p(.) Ocu(Cobert_her) 3 203.59 13.73 0.00098 1.00
## p(.) Ocu(Frailejones) 3 206.86 17.00 0.00019 1.00Que tan bueno es el ajuste del modelo?
## t0 = 20.26135
##
## Call:
## occu(formula = ~1 ~ DAP + Alti, data = umf_Q)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) 7.56 20.5 0.370 0.712
## DAP -140.75 320.8 -0.439 0.661
## Alti -59.39 142.2 -0.418 0.676
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## -3.31 0.222 -14.9 3.04e-50
##
## AIC: 189.8538
## Number of sites: 28
## optim convergence code: 0
## optim iterations: 121
## Bootstrap iterations: 0newdat_range<-data.frame(DAP=seq(min(Covs_q$DAP),max(Covs_q$DAP),length=100),
Alti=seq(min(Covs_q$Alt),max(Covs_q$Alt),length=100))
### plot Detection en escala estandarizada
# pred_det <-predict(m9, type="det", newdata=newdat_range, appendData=TRUE)
# plot(Predicted~Alt, pred_det,type="l",col="blue",
# xlab="Altitud",
# ylab="Predicted detection prob.")
# lines(lower~Alt, pred_det,type="l",col=gray(0.5))
# lines(upper~Alt, pred_det,type="l",col=gray(0.5))Usamos el modelo m9 para predecir.
### plot occupancy en escala estandarizada
pred_psi <-predict(m9, type="state", newdata=newdat_range, appendData=TRUE)
plot(Predicted~DAP, pred_psi,type="l",col="blue",
xlab="DAP",
ylab="Predicted Occupancy prob.")
lines(lower~DAP, pred_psi,type="l",col=gray(0.5))
lines(upper~DAP, pred_psi,type="l",col=gray(0.5))
### plot occupancy en escala estandarizada
pred_psi <-predict(m9, type="state", newdata=newdat_range, appendData=TRUE)
plot(Predicted~Alti, pred_psi,type="l",col="blue",
xlab="Altitud",
ylab="Predicted Occupancy prob.")
lines(lower~Alti, pred_psi,type="l",col=gray(0.5))
lines(upper~Alti, pred_psi,type="l",col=gray(0.5))
Problema! si bien DAP y Altitud son las covariables que mejor predicen la ocupacion, hay un error muy grande en la predicción. La posible causa son los pocos datos (28 sitios)
Juntemos los datos de Risaralda y Quindio
#Construir modelos
# nombre camara en quindio
row.names(y_Q) <- Quindio_data$camara
# nuevo Risaralda con mas NAs en columnas
y_Rn <- cbind(y_R, as.data.frame(matrix(NA,59,12))) # Nuevo Risaralda con 12 NAs
#names(y_Rn) <- NULL
#names(y_Q) <- NULL
# Loop to convert logical to numeric
for (i in 1:length(y_Rn)){
y_Rn[,i] <- as.numeric(y_Rn[,i])
}
names(y_Rn) <- c(1:108)
# Loop to convert logical to numeric
for (i in 1:length(y_Q)){
y_Q[,i] <- as.numeric(y_Q[,i])
}
names(y_Q) <- c(1:108)
# unir Risaralda y Quindio
y_full <- rbind (y_Rn, y_Q)
full_covs_s <- data.frame(scale(data.frame(Elev=full_covs$Elev)))
# Make unmarked frame
umf_y_full<- unmarkedFrameOccu(y= y_full)
siteCovs(umf_y_full) <- full_covs_s # data.frame(Elev=full_covs$Elev) # Full
#######Graficar umf
plot(umf_y_full)
Un modelo simple con la unica covariable comun a Quindio y Risaralda
Altitud en graficas
# build models
mf0<-occu(~1 ~ 1, umf_y_full)
mf1<-occu(~1 ~ Elev, umf_y_full)
mf2<-occu(~1 ~ Elev +I(Elev^2), umf_y_full)
mf3<-occu(~1 ~ Elev +I(Elev^3), umf_y_full)
# fit list
fms1<-fitList("p(.) Ocu(.)"=mf0,
"p(.) Ocu(Elev)"=mf1,
"p(.) Ocu(Elev^2)"=mf2,
"p(.) Ocu(Elev^3)"=mf3
)
modSel(fms1)## nPars AIC delta AICwt cumltvWt
## p(.) Ocu(Elev^2) 4 368.33 0.00 6.6e-01 0.66
## p(.) Ocu(Elev) 3 369.73 1.40 3.3e-01 0.99
## p(.) Ocu(Elev^3) 4 376.33 8.00 1.2e-02 1.00
## p(.) Ocu(.) 2 417.05 48.72 1.7e-11 1.00## t0 = 36.43499
newdat_range<-data.frame(Elev=seq(min(full_covs_s$Elev),max(full_covs_s$Elev),length=100))
### plot Detection en escala estandarizada
# pred_det <-predict(m4, type="det", newdata=newdat_range, appendData=TRUE)
# plot(Predicted~Elev, pred_det,type="l",col="blue",
# xlab="Altitud",
# ylab="Predicted detection prob.")
# lines(lower~Elev, pred_det,type="l",col=gray(0.5))
# lines(upper~Elev, pred_det,type="l",col=gray(0.5))
### plot occupancy en escala estandarizada
pred_psi <-predict(mf2, type="state", newdata=newdat_range, appendData=TRUE)
plot(Predicted~Elev, pred_psi,type="l",col="blue",
xlab="Altitud",
ylab="Predicted Occupancy prob.")
lines(lower~Elev, pred_psi,type="l",col=gray(0.5))
lines(upper~Elev, pred_psi,type="l",col=gray(0.5))
### Plot detection en escala original
plot(Predicted ~ Elev, pred_psi, type="l", ylim=c(0,1), col="blue",
xlab="Altitud",
ylab="Predicted Occupancy prob.",
xaxt="n")
xticks <- -1:2
xlabs <- xticks*sd(full_covs$Elev) + mean(full_covs$Elev) #Use the mean and sd of original value to change label name
axis(1, at=xticks, labels=round(xlabs, 1))
lines(lower ~ Elev, pred_psi, type="l", col=gray(0.5))
lines(upper ~ Elev, pred_psi, type="l", col=gray(0.5))
La altitud es buena predictora de la ocupación. la ocupación aumente con la altitud hasta los 3100 metros, luego de esta altitud la ocupación disminuye lentamente.
Modelo de Ocupación explicado por altitud espacialmente explicito.
library(RColorBrewer)
srtm_crop_s <- stack(scale(srtm_crop)) # scale altitud
names(srtm_crop_s) <- "Elev"
crs(srtm_crop_s) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
pred_psi_s <-predict(mf2, type="state", newdata=srtm_crop_s) ## doing row 1000 of 66528
## doing row 2000 of 66528
## doing row 3000 of 66528
## doing row 4000 of 66528
## doing row 5000 of 66528
## doing row 6000 of 66528
## doing row 7000 of 66528
## doing row 8000 of 66528
## doing row 9000 of 66528
## doing row 10000 of 66528
## doing row 11000 of 66528
## doing row 12000 of 66528
## doing row 13000 of 66528
## doing row 14000 of 66528
## doing row 15000 of 66528
## doing row 16000 of 66528
## doing row 17000 of 66528
## doing row 18000 of 66528
## doing row 19000 of 66528
## doing row 20000 of 66528
## doing row 21000 of 66528
## doing row 22000 of 66528
## doing row 23000 of 66528
## doing row 24000 of 66528
## doing row 25000 of 66528
## doing row 26000 of 66528
## doing row 27000 of 66528
## doing row 28000 of 66528
## doing row 29000 of 66528
## doing row 30000 of 66528
## doing row 31000 of 66528
## doing row 32000 of 66528
## doing row 33000 of 66528
## doing row 34000 of 66528
## doing row 35000 of 66528
## doing row 36000 of 66528
## doing row 37000 of 66528
## doing row 38000 of 66528
## doing row 39000 of 66528
## doing row 40000 of 66528
## doing row 41000 of 66528
## doing row 42000 of 66528
## doing row 43000 of 66528
## doing row 44000 of 66528
## doing row 45000 of 66528
## doing row 46000 of 66528
## doing row 47000 of 66528
## doing row 48000 of 66528
## doing row 49000 of 66528
## doing row 50000 of 66528
## doing row 51000 of 66528
## doing row 52000 of 66528
## doing row 53000 of 66528
## doing row 54000 of 66528
## doing row 55000 of 66528
## doing row 56000 of 66528
## doing row 57000 of 66528
## doing row 58000 of 66528
## doing row 59000 of 66528
## doing row 60000 of 66528
## doing row 61000 of 66528
## doing row 62000 of 66528
## doing row 63000 of 66528
## doing row 64000 of 66528
## doing row 65000 of 66528
## doing row 66000 of 66528Session Info
Details and pakages used
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
##
## Matrix products: default
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] RColorBrewer_1.1-2 plyr_1.8.5 lubridate_1.7.4 tmaptools_2.0-2
## [5] tmap_2.3-2 mapview_2.7.0 readxl_1.3.1 forcats_0.4.0
## [9] stringr_1.4.0 dplyr_0.8.3 purrr_0.3.3 readr_1.3.1
## [13] tidyr_1.0.0 tibble_2.1.3 ggplot2_3.2.1 tidyverse_1.3.0
## [17] unmarked_0.13-1 reshape2_1.4.3 Rcpp_1.0.3 lattice_0.20-38
## [21] sf_0.8-0 raster_3.0-7 sp_1.3-2
##
## loaded via a namespace (and not attached):
## [1] leafem_0.0.1 colorspace_1.4-1 class_7.3-15
## [4] leaflet_2.0.3 rgdal_1.4-8 satellite_1.0.2
## [7] base64enc_0.1-3 fs_1.3.1 dichromat_2.0-0
## [10] rstudioapi_0.10 farver_2.0.3 fansi_0.4.1
## [13] xml2_1.2.2 codetools_0.2-16 knitr_1.27
## [16] zeallot_0.1.0 jsonlite_1.6 broom_0.5.3
## [19] dbplyr_1.4.2 png_0.1-7 rgeos_0.5-2
## [22] shiny_1.4.0 compiler_3.6.1 httr_1.4.1
## [25] backports_1.1.5 assertthat_0.2.1 Matrix_1.2-17
## [28] fastmap_1.0.1 lazyeval_0.2.2 cli_2.0.1
## [31] later_1.0.0 leaflet.providers_1.9.0 htmltools_0.4.0
## [34] tools_3.6.1 gtable_0.3.0 glue_1.3.1
## [37] cellranger_1.1.0 vctrs_0.2.1 svglite_1.2.2
## [40] nlme_3.1-140 leafsync_0.1.0 crosstalk_1.0.0
## [43] lwgeom_0.1-7 xfun_0.12 rvest_0.3.5
## [46] mime_0.8 lifecycle_0.1.0 XML_3.98-1.20
## [49] MASS_7.3-51.4 scales_1.1.0 hms_0.5.3
## [52] promises_1.1.0 yaml_2.2.0 gdtools_0.2.1
## [55] stringi_1.4.4 leafpop_0.0.5 e1071_1.7-3
## [58] rlang_0.4.2 pkgconfig_2.0.3 systemfonts_0.1.1
## [61] evaluate_0.14 htmlwidgets_1.5.1 tidyselect_0.2.5
## [64] magrittr_1.5 R6_2.4.1 generics_0.0.2
## [67] DBI_1.1.0 pillar_1.4.3 haven_2.2.0
## [70] withr_2.1.2 units_0.6-5 modelr_0.1.5
## [73] crayon_1.3.4 uuid_0.1-2 KernSmooth_2.23-15
## [76] rmarkdown_2.1 grid_3.6.1 reprex_0.3.0
## [79] digest_0.6.23 classInt_0.4-2 webshot_0.5.2
## [82] xtable_1.8-4 httpuv_1.5.2 brew_1.0-6
## [85] stats4_3.6.1 munsell_0.5.0 viridisLite_0.3.0