Cats Ucumari
Juan Camilo Cepada-Duque, Diego J. Lizcano, Bibiana Gomez-Valencia, John Harold Castaño
1 de marzo de 2019
R Load Packages and Data
Load Packages
library (readxl) # read excel data
library (tidyverse) # Data handling
library (tmap) # nice maps
library (unmarked) # occu analitics
library (ggcorrplot) # corplot in ggplot
library (raster)
library (rasterVis)
library (knitr) #
library (kableExtra) # nice tables
library (sf) #sf
library (mapview)Load data
#Covariables de ocupacion
siteCovs<- read.table('C:/Users/diego.lizcano/Box Sync/CodigoR/Camilo/data/covariates.txt', sep='\t', header=T)
#########Subir los datos e la matriz de presencia-ausencia
y<- read.table('C:/Users/diego.lizcano/Box Sync/CodigoR/Camilo/data/PUMA_MATRIX.txt', sep = '\t', header=T )
######Indexar matriz de manera que solo me considere los 1's y 0's
y<-y[ ,2:19]
# see covs
kable(siteCovs) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
scroll_box(width = "100%", height = "500px")| CAMERA_NAME | Y | X | POPDEN_MEAN | CANOPY_MEAN | HFP | EVI_MEAN | NDWI_.MEAN | PRECIPITATION | NDVI_MEAN | NPP_MEAN | GPP_MEAN |
|---|---|---|---|---|---|---|---|---|---|---|---|
| CT-UC-01-02 | 4.730583 | -75.53442 | 0.00000 | 33 | 21 | 0.46940 | 0.02013 | 2430 | 0.70748 | 12057 | 21984 |
| CT-UC-01-03 | 4.713333 | -75.54614 | 5.42337 | 34 | 28 | 0.52759 | 0.01823 | 2670 | 0.75493 | 12552 | 23565 |
| CT-UC-01-07 | 4.732139 | -75.53922 | 0.00000 | 33 | 21 | 0.46940 | 0.02013 | 2430 | 0.70748 | 12057 | 21984 |
| CT-UC-01-09 | 4.774028 | -75.57292 | 4.66756 | 17 | 21 | 0.44093 | 0.01825 | 2823 | 0.66879 | 10216 | 18684 |
| CT-UC-01-11 | 4.797861 | -75.58356 | 26.81536 | 6 | 26 | 0.47913 | 0.01013 | 2835 | 0.68954 | 9508 | 16716 |
| CT-UC-01-13 | 4.768861 | -75.58131 | 49.43168 | 22 | 21 | 0.52437 | 0.00925 | 2754 | 0.72856 | 11698 | 21486 |
| CT-UC-01-16 | 4.721056 | -75.54000 | 0.12275 | 34 | 28 | 0.49353 | 0.02253 | 2034 | 0.71808 | 11641 | 21053 |
| CT-UC-01-17 | 4.728083 | -75.56258 | 1.61679 | 42 | 28 | 0.53002 | 0.01678 | 2745 | 0.73281 | 13351 | 24862 |
| CT-UC-01-21 | 4.735333 | -75.53575 | 0.20614 | 28 | 21 | 0.53724 | 0.01972 | 2471 | 0.74700 | 12254 | 23090 |
| CT-UC-01-24 | 4.724750 | -75.55342 | 1.99266 | 33 | 28 | 0.53012 | 0.02375 | 2692 | 0.74128 | 12618 | 23530 |
| CT-UC-01-25 | 4.723778 | -75.54900 | 1.99266 | 33 | 21 | 0.53012 | 0.02375 | 2379 | 0.74128 | 11823 | 21950 |
| CT-UC-01-27 | 4.720111 | -75.55761 | 8.55517 | 34 | 28 | 0.52556 | 0.01994 | 2692 | 0.71446 | 12618 | 23530 |
| CT-UC-01-31 | 4.776083 | -75.55733 | 2.97775 | 7 | 21 | 0.49805 | 0.01140 | 2077 | 0.68476 | 3182 | 6306 |
| CT-UC-02-01 | 4.733611 | -75.56972 | 229.13182 | 39 | 28 | 0.53166 | 0.01719 | 2697 | 0.72762 | 13074 | 25265 |
| CT-UC-02-02 | 4.738389 | -75.56989 | 100.70332 | 35 | 28 | 0.56691 | 0.01508 | 2697 | 0.71679 | 13074 | 25265 |
| CT-UC-02-04 | 4.738500 | -75.56094 | 9.44532 | 41 | 28 | 0.53153 | 0.02332 | 2670 | 0.72343 | 12954 | 24696 |
| CT-UC-02-06 | 4.751250 | -75.58108 | 28.46955 | 33 | 28 | 0.52114 | 0.01739 | 2818 | 0.69203 | 12238 | 22649 |
| CT-UC-02-08 | 4.751056 | -75.57167 | 1.71999 | 28 | 28 | 0.52473 | 0.02604 | 2196 | 0.70765 | 10666 | 18988 |
| CT-UC-02-11 | 4.759333 | -75.56731 | 6.78518 | 29 | 21 | 0.45973 | 0.01302 | 2755 | 0.70926 | 10318 | 18094 |
| CT-UC-02-13 | 4.762722 | -75.55858 | 1.44236 | 15 | 21 | 0.49707 | 0.02643 | 1988 | 0.69914 | 6035 | 12023 |
| CT-UC-02-14 | 4.700472 | -75.53417 | 13.56040 | 38 | 28 | 0.50501 | 0.02164 | 2550 | 0.75940 | 13205 | 24863 |
| CT-UC-02-18 | 4.705389 | -75.51744 | 3.70714 | 39 | 28 | 0.48151 | 0.01089 | 2179 | 0.71408 | 12259 | 22858 |
| CT-UC-02-21 | 4.708250 | -75.51331 | 1.17315 | 34 | 28 | 0.45813 | 0.01449 | 2071 | 0.75737 | 12434 | 22340 |
| CT-UC-02-23 | 4.705111 | -75.50372 | 1.48404 | 35 | 28 | 0.48113 | 0.00474 | 2020 | 0.73521 | 10880 | 19067 |
| CT-UC-02-25 | 4.723056 | -75.56886 | 1.69512 | 36 | 28 | 0.51264 | 0.01749 | 2734 | 0.72597 | 12773 | 24305 |
| CT-UC-02-26 | 4.709361 | -75.48961 | 57.04740 | 33 | 28 | 0.50781 | 0.00129 | 1945 | 0.76352 | 11636 | 21599 |
| CT-UC-02-28 | 4.700222 | -75.51733 | 3.70714 | 39 | 28 | 0.48151 | 0.01089 | 2179 | 0.71408 | 12259 | 22858 |
| CT-UC-02-31 | 4.698225 | -75.51256 | 3.70714 | 39 | 28 | 0.48151 | 0.01089 | 2059 | 0.71408 | 12529 | 22925 |
| CT-UC-02-32 | 4.709861 | -75.55786 | 41.67831 | 27 | 28 | 0.49391 | 0.02116 | 2711 | 0.69654 | 12816 | 24068 |
| CT-UC-02-33 | 4.708861 | -75.54864 | 3.29384 | 28 | 28 | 0.52814 | 0.01874 | 2670 | 0.72200 | 12552 | 23565 |
See Cameras
# put coords
covs_sf <- st_as_sf(siteCovs, coords = c("X", "Y"), crs = "+proj=longlat +datum=WGS84")
cams_plot <- mapview(covs_sf, zcol = "CAMERA_NAME",
color = "black",
legend = FALSE,
map.types = "Esri.WorldImagery")
cams_plotcams_sp <- as(covs_sf, "Spatial")
centroid <- c(mean(siteCovs$X), mean(siteCovs$Y))
##################################
# set extent get SRTM
##################################
clip_window <- extent(-75.61 , -75.46, 4.65, 4.85)
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)
siteCovs$Elev <- extract(srtm_crop, cams_sp)
plot(srtm_crop)
plot(cams_sp, add=TRUE)
See Correlations in Site Covs
#### See correlation in Covs
# correlation mat
corr <- round(cor(siteCovs[,4:13]), 1)
ggcorrplot(corr,
#method = "circle",
hc.order = TRUE,
type = "lower",
lab = TRUE)
#### Estandarizar Covs
proms<-apply(siteCovs[,4:13], MARGIN = 2,FUN = mean)
des<-apply(siteCovs[,4:13], MARGIN = 2, FUN = sd)
cov_stand<-t(apply(siteCovs[,4:13], MARGIN = 1, FUN = function(x) (x-proms)/des))
Covariables<-as.data.frame(cov_stand)Unmarked frame object
#Correr el modelo de ocupacion
umf<- unmarkedFrameOccu(y= y, siteCovs = Covariables)
#######Graficar umf
plot(umf)
#Construir modelos
#Correr el modelo nulo donde la detecci?n (primer~1) y la ocupaci?n son constantes (segundo~1)
#Build Models and Model Selection
#Construir modelos
#Correr el modelo nulo donde la detecci?n (primer~1) y la ocupaci?n son constantes (segundo~1)
#
# PRIMERO HAGO MODELO NULO CAMBIANDO LA DETECCION CON LAS COVARIABLES DE SITIO QUE LA ALTEREN
# HACER CRITERIO DE AKAIKE PARA TODOS LOS MODELOS DE DETECCION
# ES MEJOR QUEDARME INICIALMENTE CON EL MEJOR MODELO DE DETECCI?N Y LUEGO PASAR A MODELAR OCUPACI?N
# PRIMERO SALDRIAN 4 MODELOS TODOS CON OCUPACION CONSTANTE.
m0<-occu(~1 ~1, umf)
m1<-occu(~PRECIPITATION ~1, umf)
m2<-occu(~1 ~POPDEN_MEAN + HFP, umf)
m3<-occu(~PRECIPITATION ~Elev+NDWI_.MEAN, umf)
#m4<-occu (~Elev ~NDWI_.MEAN+Elev, umf)
m5<-occu(~PRECIPITATION ~POPDEN_MEAN, umf)
m6<-occu (~CANOPY_MEAN ~NDWI_.MEAN, umf)
m7<-occu (~PRECIPITATION ~CANOPY_MEAN, umf)
m8<-occu (~1 ~POPDEN_MEAN + NDWI_.MEAN, umf)
m9<-occu (~1 ~Elev, umf)
m10<-occu (~CANOPY_MEAN ~NDWI_.MEAN , umf)
#m11<-occu (~PRECIPITATION ~NDWI_.MEAN , umf)
m12<-occu (~Elev ~NDWI_.MEAN , umf)
# Model names
fms<-fitList("p(.) Ocu(.)"=m0,
"p(precipitation) Ocu(.)"=m1,
"p(.) Ocu(POPDEN+HFP)"=m2,
"p(precipitation) Ocu(NDVI+Elev)"=m3,
#"p(Elev) Ocu(NDVI+Elev)"=m4,
"p(precipitation) Ocu(POPDEN)"=m5,
"p(Canopy) Ocu(NDWI)"=m6,
"p(precipitation) Ocu(Canopy)"=m7,
"p(.) Ocu(POPDEN+NDWI)"=m8,
"p(.) Ocu(Elev)"=m9,
"p(Canopy) Ocu(NDVI)"=m10,
#"p(precipitation) Ocu(NDVI)"=m11,
"p(Elev) Ocu(NDVI)"=m12)
seleccion_modelos<-modSel(fms)
seleccion_modelos## nPars AIC delta AICwt cumltvWt
## p(precipitation) Ocu(NDVI+Elev) 5 147.24 0.00 0.4385 0.44
## p(Elev) Ocu(NDVI) 4 149.71 2.47 0.1277 0.57
## p(Canopy) Ocu(NDWI) 4 150.05 2.80 0.1080 0.67
## p(Canopy) Ocu(NDVI) 4 150.05 2.80 0.1080 0.78
## p(.) Ocu(POPDEN+NDWI) 4 150.40 3.15 0.0905 0.87
## p(.) Ocu(Elev) 3 151.56 4.32 0.0507 0.92
## p(.) Ocu(.) 2 152.92 5.68 0.0256 0.95
## p(precipitation) Ocu(.) 3 153.38 6.13 0.0204 0.97
## p(precipitation) Ocu(POPDEN) 4 154.36 7.11 0.0125 0.98
## p(precipitation) Ocu(Canopy) 4 154.50 7.25 0.0117 0.99
## p(.) Ocu(POPDEN+HFP) 4 155.69 8.44 0.0064 1.00Best model is m3
### m3 is best
summary (m3)##
## Call:
## occu(formula = ~PRECIPITATION ~ Elev + NDWI_.MEAN, data = umf)
##
## Occupancy (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -0.164 0.540 -0.303 0.7618
## Elev 1.179 0.717 1.643 0.1003
## NDWI_.MEAN 1.463 0.665 2.199 0.0279
##
## Detection (logit-scale):
## Estimate SE z P(>|z|)
## (Intercept) -1.139 0.268 -4.24 0.000022
## PRECIPITATION -0.314 0.270 -1.16 0.244859
##
## AIC: 147.2429
## Number of sites: 30
## optim convergence code: 0
## optim iterations: 22
## Bootstrap iterations: 0pb <- parboot(m3, nsim=500, report=10) # goodness of fit## t0 = 20.43962
## Running in parallel. Bootstrapped statistics not reported.plot (pb)
Seems to be a good model to predict.
Predictions
# mat as data frame
Cov_stand <- as.data.frame(cov_stand)
# new data
elevrange<-data.frame(Elev=seq(min(Cov_stand$Elev),
max(Cov_stand$Elev),length=100))
NDVIrange <- data.frame(NDWI_.MEAN=seq(min(Cov_stand$NDWI_.MEAN),
max(Cov_stand$NDWI_.MEAN),length=100))
Precrange <- data.frame(PRECIPITATION=seq(min(Cov_stand$PRECIPITATION),
max(Cov_stand$PRECIPITATION),length=100))
# new dat in a single dataframe
newdat <- cbind(elevrange, NDVIrange, Precrange)
# plot and predict Occu
pred_psi <-predict(m3,type="state",newdata=newdat,appendData=TRUE)
plot(Predicted~Elev, pred_psi,type="l",col="blue",
xlab="Elev",
ylab="psi")
lines(lower~Elev, pred_psi,type="l",col=gray(0.5))
lines(upper~Elev, pred_psi,type="l",col=gray(0.5))
# plot and predict det
pred_p <-predict(m3,type="det",newdata=newdat,appendData=TRUE)
plot(Predicted~PRECIPITATION, pred_p,type="l",col="blue",
xlab="Elev",
ylab="p")
lines(lower~PRECIPITATION, pred_p,type="l",col=gray(0.5))
lines(upper~PRECIPITATION, pred_p,type="l",col=gray(0.5))
En conclusión
Un modelo con menos interacciones entre las covariables es mas facil de interpretar.
La precipitación y elevación estan bastante correlacionadas, vale la pena dejar una de las dos… yo pensaria que es mejor quedarse con la elevación.
En resumen
- Estoy trabajando ya en el análisis de patrones de actividad.
Bien!… porfa mandeme el codigo.
- Ya transformé los modelos de escala logit a escala de probabilidad.
- No en nesesario, por ahora. Mas adelante se hace, por ahora falta incluir una mejor covariable de probabilidad de deteccion.
- Necesito incluir la covariable de días de muestreo como función de la detección en los modelos.
- No se requiere. Ese tipo de covariables no son utiles y no dicen mucho.
- Necesito realizar la predicción de mis modelos y elaborar las gráficas y el mapa correspondiente para poder trabajar de lleno en el manuscrito.
- ya lo hice.
Session Info
Details and pakages used
print(sessionInfo(), locale = FALSE)## R version 3.5.2 (2018-12-20)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17134)
##
## Matrix products: default
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] mapview_2.6.3 sf_0.7-2 kableExtra_1.0.1.9000
## [4] knitr_1.21 rasterVis_0.45 latticeExtra_0.6-28
## [7] RColorBrewer_1.1-2 raster_2.8-19 sp_1.3-1
## [10] ggcorrplot_0.1.2 unmarked_0.12-3 reshape2_1.4.3
## [13] Rcpp_1.0.0 lattice_0.20-38 tmap_2.2
## [16] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.0.1
## [19] purrr_0.3.0 readr_1.3.1 tidyr_0.8.2
## [22] tibble_2.0.1 ggplot2_3.1.0 tidyverse_1.2.1
## [25] readxl_1.2.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-137 satellite_1.0.1 lubridate_1.7.4
## [4] webshot_0.5.1 httr_1.4.0 tools_3.5.2
## [7] backports_1.1.3 rgdal_1.3-6 R6_2.4.0
## [10] KernSmooth_2.23-15 rgeos_0.4-2 DBI_1.0.0
## [13] lazyeval_0.2.1 colorspace_1.4-0 withr_2.1.2
## [16] tidyselect_0.2.5 leaflet_2.0.2 compiler_3.5.2
## [19] cli_1.0.1 rvest_0.3.2 xml2_1.2.0
## [22] labeling_0.3 scales_1.0.0 hexbin_1.27.2
## [25] classInt_0.3-1 digest_0.6.18 rmarkdown_1.11
## [28] base64enc_0.1-3 dichromat_2.0-0 pkgconfig_2.0.2
## [31] htmltools_0.3.6 highr_0.7 htmlwidgets_1.3
## [34] rlang_0.3.1 rstudioapi_0.9.0 shiny_1.2.0
## [37] generics_0.0.2 zoo_1.8-4 jsonlite_1.6
## [40] crosstalk_1.0.0 magrittr_1.5 Matrix_1.2-15
## [43] munsell_0.5.0 stringi_1.3.1 yaml_2.2.0
## [46] tmaptools_2.0-1 plyr_1.8.4 grid_3.5.2
## [49] promises_1.0.1 crayon_1.3.4 haven_2.1.0
## [52] hms_0.4.2 pillar_1.3.1 codetools_0.2-16
## [55] stats4_3.5.2 XML_3.98-1.16 glue_1.3.0
## [58] evaluate_0.13 modelr_0.1.4 png_0.1-7
## [61] httpuv_1.4.5.1 cellranger_1.1.0 gtable_0.2.0
## [64] assertthat_0.2.0 xfun_0.5 mime_0.6
## [67] lwgeom_0.1-5 xtable_1.8-3 broom_0.5.1
## [70] e1071_1.7-0.1 later_0.7.5 class_7.3-15
## [73] viridisLite_0.3.0 units_0.6-2