Multispecies occupancy model from Rota et al. (2016), for two (or more) interacting species. The model generalizes the standard single-species occupancy model from MacKenzie et al. (2002). We used Leopardus tigrinus and domestic dog data from 6 sites sampled with arrays of camera traps.

cargar paquetes

codigo R

library(knitr)
library(tidyverse) # maneja datos
library(mapview) # mapas facil
library(readxl) #leer datos
library(sf) # vector map
library(geodata) # replace getData de raster para Terra
library(raster) # mapas raster
library(spatstat) # interpola mapa
library(maptools) # to coerce to ppp. note that 'maptools' will be retired by the end of 2023
library(rgdal) # rgdal will be retired during 2023 #some tricks to change projection
# library(stars)
# library(unmarked) # occu models
library(DT) # html table


source("C:/CodigoR/tigrinus2/R/organizadato.R")

UCUMARI

cargar datos

codigo R

Full_data_ucu <- read_excel("D:/BoxFiles/Box Sync/CodigoR/tigrinus/data/Full_data_Ucumari_Huila_Cocha1_Cocha2.xlsx", 
    sheet = "ucumari", col_types = c("numeric", 
        "text", "text", "text", "text", "text", 
        "text", "text", "text", "text", "text", 
        "text", "text", "text", "text", "numeric", 
        "text", "text", "text", "text", "text", 
        "text", "text", "text", "text", "numeric", 
        "numeric", "numeric", "numeric", 
        "text"))

crerar matrices para unmarked

codigo R

casas <- read_csv("C:/CodigoR/tigrinus2/data/casas.csv")
casas_sf <- st_as_sf(casas, coords = c("lon", "lat"), crs = "EPSG:4326")




############# start spatial part
#### make sf object
ucumari_sf <- st_as_sf(Full_data_ucu, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")


# get centroid. 1st bbox, make pol
centroid_ucu <- st_centroid(st_as_sfc(st_bbox(ucumari_sf)))
# get altitude
elev_ucu_full <- elevation_3s(centroid_ucu[[1]][1], centroid_ucu[[1]][2], 
                         path="D:/BoxFiles/Box Sync/CodigoR/tigrinus/raster")

# elev_ucu_full_ras <- elev_ucu_full %>% raster()

# st_bbox(ucumari_sf) # notice order xmin, xmax, ymin, ymax
ext_ucu <- ext(-75.59, -75.47,  4.68,    4.81 )
elev_ucu <- crop(elev_ucu_full, ext_ucu) 

# convert from terra to raster
elev_ucu_ras <-  elev_ucu %>% raster()

# get uniques
cams_ucu <-Full_data_ucu %>% dplyr::select(c("Longitude", "Latitude", "camera_trap")) %>% distinct()
#### make sf object
cams_ucu_sf <- st_as_sf(cams_ucu, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")

# extract values from raster using altitude
cams_ucu_sf$elev <- extract(elev_ucu_ras, cams_ucu_sf)

############## make distance map
# Convert points to sp spatialpointdatafram
casas_points <- as(casas_sf, "Spatial")
# Projection
# Be sure to have rgdal first installed.
casas_points_utm <- spTransform(casas_points, CRS('+init=epsg:21818'))
# convert to ppp
casas_points_ppp <- as(as(casas_points_utm, "SpatialPoints"), "ppp")
# distance
casas_distance <- distmap(casas_points_ppp)

####### extract distance 

casas_distance_ras<- raster(casas_distance) # convert raster
crs(casas_distance_ras) <- '+init=epsg:21818' # add crs

# project raster
casa_distance <- projectRaster(casas_distance_ras, crs = crs(casas_points))

# cams_ucu_sf$dist_casa <- raster::extract(casas_distance_ras, cams_ucu_sf)
cams_ucu_sf$dist_casa <- raster::extract(casa_distance, cams_ucu_sf) # also works


# plot map
pal = mapviewPalette("mapviewTopoColors") #color palete

mapview(elev_ucu_ras) + 
  mapview(casa_distance, col.regions = pal(100), at = seq(600, 240000, 100), legend = TRUE) +
  mapview(cams_ucu_sf["camera_trap"])

codigo R

############### end spatial part

Full_data_ucu$binomial <- str_c (Full_data_ucu$Genus, "_", Full_data_ucu$Species)

#funcion para crear todas las tablas de datos
all_data_ucu <-  f.matrix.creator2 (Full_data_ucu)

# names(all_data) # ver lass especies y en que lista esta cada una
# kable(names(all_data)) # html table
# Tigrinus es lista 8

datatable(
data = as.data.frame(names(all_data_ucu)),
caption = "Especies Ucumari",
filter = "top"
)

unmarked

codigo R

library(unmarked)
# tabla con solo tiginus
tigrinus_ucu <- all_data_ucu[[8]]
# cargar paquete
library(unmarked)
# crear objeto umf
umf_tigrinus_ucu <- unmarkedFrameOccu(y=tigrinus_ucu)
# verificar datos en grafica
# plot(umf_tigrinus_ucu)

# tabla con solo perros
perros_ucu <- all_data_ucu[[75]]
# crear objeto umf
umf_perros_ucu <- unmarkedFrameOccu(y=perros_ucu)
# verificar datos en grafica
# plot(umf_perros_ucu)

# tabla con solo ocelote
ocelote_ucu <- all_data_ucu[[16]]
# crear objeto umf
umf_ocelote_ucu <- unmarkedFrameOccu(y=ocelote_ucu)
# verificar datos en grafica
# plot(umf_ocelote_ucu)

Modelo nulo

codigo R

# modelo nulo tigrinus
fm_tig_ucu <- occu(~1 ~1, umf_tigrinus_ucu)  # fit a model

backTransform(fm_tig_ucu, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.0184 0.00977 -3.98 1

Transformation: logistic

codigo R

backTransform(fm_tig_ucu, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.198 0.0968 -1.4 1

Transformation: logistic

codigo R

# modelo nulo perro
fm_perros_ucu <- occu(~1 ~1, umf_perros_ucu)  # fit a model

backTransform(fm_perros_ucu, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.000356 0.000614 -7.94 1

Transformation: logistic

codigo R

backTransform(fm_perros_ucu, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.985 1.38 4.17 1

Transformation: logistic

codigo R

# modelo nulo ocelote
fm_ocelote_ucu <- occu(~1 ~1, umf_ocelote_ucu)  # fit a model

backTransform(fm_ocelote_ucu, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.0109 0.00736 -4.51 1

Transformation: logistic

codigo R

backTransform(fm_ocelote_ucu, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.288 0.181 -0.904 1

Transformation: logistic

MODELOS DE CO-OCURRENCIA

codigo R

detformulas <- c( "~1", "~1")#, "~1")
#stateformulas <- c('~elev','~elev', '~elev', "~1", "~1", "~1", "~0")# 3 sp
stateformulas <- c('~elev','~elev', "~0")
y <- list(tigrinus_ucu, perros_ucu)# , ocelote_ucu)
names(y) <- c("tigrinus", "perros")#, "ocelote")


obs_covs <-as.data.frame(scale(cams_ucu_sf$dist_casa))
names(obs_covs) <- "dist_casa"

site_covs_ucu <- data.frame(cams_ucu_sf[,c('elev','dist_casa')])[,1:2]
site_covs_ucu <-as.data.frame(apply(site_covs_ucu,2,scale))
names(site_covs_ucu) <- c("elev", "dist_casa")


umf <-  unmarkedFrameOccuMulti(y=y, 
                              siteCovs=site_covs_ucu,
                              obsCovs = NULL)
plot(umf)

codigo R

#umf

# occFormulas Length should match number/order of columns in fDesign
umf@fDesign

    f1[tigrinus] f2[perros] f3[tigrinus:perros]

psi[11] 1 1 1 psi[10] 1 0 0 psi[01] 0 1 0 psi[00] 0 0 0

codigo R

fit1 <- occuMulti(detformulas, stateformulas, umf,    
        method="BFGS", se=TRUE, engine=c("C"), silent=FALSE)

fit1

Call: occuMulti(detformulas = detformulas, stateformulas = stateformulas, data = umf, method = “BFGS”, se = TRUE, engine = c(“C”), silent = FALSE, maxOrder = 2L)

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) -1.3964 0.611 -2.286 0.0223 [tigrinus] elev 0.0983 0.456 0.215 0.8295 [perros] (Intercept) -48.0649 100.952 -0.476 0.6340 [perros] elev 38.5090 77.990 0.494 0.6215

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -3.98 0.543 -7.33 2.22e-13 [perros] (Intercept) -5.77 1.105 -5.22 1.79e-07

AIC: 152.5147

codigo R

# update model
# occFormulas2 <- c('~dist_casa', '~dist_casa', '~dist_casa', "~1", "~1", "~1", "~0")
occFormulas2 <- c('~dist_casa', '~dist_casa', "~0")
fit2 <- update(fit1, stateformulas=occFormulas2)
fit2

Call: occuMulti(detformulas = c(“~1”, “~1”), stateformulas = c(“~dist_casa”, “~dist_casa”, “~0”), data = umf, maxOrder = 2L, method = “BFGS”, se = TRUE, engine = c(“C”), silent = FALSE)

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) -1.473 0.648 -2.274 0.023 [tigrinus] dist_casa 0.606 0.596 1.016 0.309 [perros] (Intercept) -22.080 35.939 -0.614 0.539 [perros] dist_casa -35.460 55.936 -0.634 0.526

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -4.00 0.552 -7.24 4.38e-13 [perros] (Intercept) -6.53 1.011 -6.46 1.04e-10

AIC: 152.6737

Model Selection

codigo R

#List of fitted models
fl <- fitList(elev=fit1, dist=fit2)
# coef(fl)


###################
# Model selection #
###################

modSel(fl)

 nPars    AIC delta AICwt cumltvWt

elev 6 152.51 0.00 0.52 0.52 dist 6 152.67 0.16 0.48 1.00

codigo R

#############
# Model fit #
#############

# bt <- parboot(fit1, nsim=50) # takes time best model
# plot(bt)

plot predicted marginal occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_ucu_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_ucu_sf$elev))/sd(cams_ucu_sf$elev)

r2 <- range(cams_ucu_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_ucu_sf$dist_casa))/sd(cams_ucu_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)



tigrinus_pred <- predict(fit1, "state", species="tigrinus", newdata=nd)
tigrinus_pred$Species <- "tigrinus"

perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Elevation", y="Marginal occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted co-occurrence occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_ucu_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_ucu_sf$elev))/sd(cams_ucu_sf$elev)

r2 <- range(cams_ucu_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_ucu_sf$dist_casa))/sd(cams_ucu_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)


tigrinus_perros_pred <- predict(fit1, "state", 
                         species=c("tigrinus", "perros"), newdata=nd)

# tigrinus_pred$Species <- c("tigrinus", "perros")

# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- tigrinus_perros_pred #rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = "grey50"), alpha=0.3) +
  geom_line(aes(y=Predicted), col="blue") +
  labs(x="Elevation", y="co-occurence") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted conditional occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_ucu_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_ucu_sf$elev))/sd(cams_ucu_sf$elev)

r2 <- range(cams_ucu_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_ucu_sf$dist_casa))/sd(cams_ucu_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)

##### conditional

tigrinus_perro_no <- predict(fit1, "state", 
                         species="tigrinus", 
                         cond='-perros',
                         newdata=nd)

tigrinus_perro_no$Species <- "perro ausente"

tigrinus_perro_si <- predict(fit1, "state", 
                         species="tigrinus", 
                         cond='perros',
                         newdata=nd)

tigrinus_perro_si$Species <- "perro presente"





# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_perro_si, tigrinus_perro_no)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Elevation", y="tigrinus conditional occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

HUILA

codigo R

huila_data <- read_excel("C:/CodigoR/tigrinus2/data/Full_data_Ucumari_Huila_Cocha1_Cocha2.xlsx", 
    sheet = "pitalito")



#huila_data$binomial <- str_c (huila_data$Genus, "_", huila_data$Species)

crerar matrices para unmarked

codigo R

############# start spatial part
#### make sf object
huila_sf <- st_as_sf(huila_data, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")


# get centroid. 1st bbox, make pol
centroid_huila <- st_centroid(st_as_sfc(st_bbox(huila_sf)))
# get altitude
elev_huila_full <- elevation_3s(centroid_huila[[1]][1], centroid_huila[[1]][2], 
                         path="D:/BoxFiles/Box Sync/CodigoR/tigrinus/raster")

# elev_ucu_full_ras <- elev_ucu_full %>% raster()

# st_bbox(ucumari_sf) # notice order xmin, xmax, ymin, ymax
ext_huila <- ext(-76.39, -76.282,  1.71,    1.84 )
elev_huila <- crop(elev_huila_full, ext_huila) 

# convert from terra to raster
elev_huila_ras <-  elev_huila %>% raster()

# get uniques
cams_huila <-huila_data %>% dplyr::select(c("Longitude", "Latitude", "camera_trap")) %>% distinct()
#### make sf object
cams_huila_sf <- st_as_sf(cams_huila, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")

# extract values from raster using altitude
cams_huila_sf$elev <- extract(elev_huila_ras, cams_huila_sf)

############## make distance map
# # Convert points to sp spatialpointdatafram
# casas_points <- as(casas_sf, "Spatial")
# # Projection
# # Be sure to have rgdal first installed.
# casas_points_utm <- spTransform(casas_points, CRS('+init=epsg:32718'))
# # convert to ppp
# casas_points_ppp <- as(as(casas_points_utm, "SpatialPoints"), "ppp")
# # distance
# casas_distance <- distmap(casas_points_ppp)

####### extract distance 
# 
# casas_distance_ras<- raster(casas_distance) # convert raster
# crs(casas_distance_ras) <- '+init=epsg:32718' # add crs
# 
# # project raster
# casa_distance <- projectRaster(casas_distance_ras, crs = crs(casas_points))

# cams_ucu_sf$dist_casa <- raster::extract(casas_distance_ras, cams_ucu_sf)
cams_huila_sf$dist_casa <- raster::extract(casa_distance, cams_huila_sf) # also works


# plot map
mapview(elev_huila_ras) + mapview(cams_huila_sf["camera_trap"])

codigo R

############### end spatial part

#huila_data$binomial <- str_c (huila_data$Genus, "_", huila_data$Species)

#funcion para crear todas las tablas de datos
all_data_huila <-  f.matrix.creator2 (huila_data)

# names(all_data) # ver lass especies y en que lista esta cada una
# kable(names(all_data)) # html table
# Tigrinus es lista 6
# perro es 43

datatable(
data = as.data.frame(names(all_data_huila)),
caption = "Especies Huila",
filter = "top"
)

unmarked

codigo R

# tabla con solo tigrinus
tigrinus_huila <- all_data_huila[[6]]

# crear objeto umf
umf_tigrinus_huila <- unmarkedFrameOccu(y=tigrinus_huila)
# verificar datos en grafica
# plot(umf_tigrinus_ucu)

# tabla con solo perros
perros_huila <- all_data_huila[[43]]
# crear objeto umf
umf_perros_huila <- unmarkedFrameOccu(y=perros_huila)
# verificar datos en grafica
# plot(umf_perros_ucu)

# # tabla con solo ocelote
# ocelote_huila <- all_data[[16]]
# # crear objeto umf
# umf_ocelote_huila <- unmarkedFrameOccu(y=ocelote_ucu)
# # verificar datos en grafica
# # plot(umf_ocelote_ucu)

Modelo nulo

codigo R

# modelo nulo tigrinus
fm_tig_huila <- occu(~1 ~1, umf_tigrinus_huila)  # fit a model

backTransform(fm_tig_huila, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.00746 0.00512 -4.89 1

Transformation: logistic

codigo R

backTransform(fm_tig_huila, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.66 0.424 0.663 1

Transformation: logistic

codigo R

# modelo nulo perro
fm_perros_huila <- occu(~1 ~1, umf_perros_huila)  # fit a model

backTransform(fm_perros_huila, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.00152 0.00076 -6.49 1

Transformation: logistic

codigo R

backTransform(fm_perros_huila, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 1 NaN 23.8 1

Transformation: logistic

codigo R

# # modelo nulo ocelote
# fm_ocelote_huila <- occu(~1 ~1, umf_ocelote_huila)  # fit a model
# 
# backTransform(fm_ocelote_huila, type="det") #estimado lineal de deteccion
# backTransform(fm_ocelote_huila, type="state") # estimado linel de ocupacion
#

MODELOS DE CO-OCURRENCIA

codigo R

detformulas <- c( "~1", "~1")#, "~1")
#stateformulas <- c('~elev','~elev', '~elev', "~1", "~1", "~1", "~0")# 3 sp
stateformulas <- c('~elev','~elev', "~0")
y <- list(tigrinus_huila, perros_huila)# , ocelote_huila)
names(y) <- c("tigrinus", "perros")#, "ocelote")


obs_covs <-as.data.frame(scale(cams_huila_sf$dist_casa))
names(obs_covs) <- "dist_casa"

site_covs_huila <- data.frame(cams_huila_sf[,c('elev','dist_casa')])[,1:2]
site_covs_huila <-as.data.frame(apply(site_covs_huila,2,scale))
names(site_covs_huila) <- c("elev", "dist_casa")


umf <-  unmarkedFrameOccuMulti(y=y, 
                              siteCovs=site_covs_huila,
                              obsCovs = NULL)
plot(umf)

codigo R

#umf

# occFormulas Length should match number/order of columns in fDesign
umf@fDesign

    f1[tigrinus] f2[perros] f3[tigrinus:perros]

psi[11] 1 1 1 psi[10] 1 0 0 psi[01] 0 1 0 psi[00] 0 0 0

codigo R

fit1 <- occuMulti(detformulas, stateformulas, umf,    
        method="BFGS", se=TRUE, engine=c("C"), silent=FALSE)

fit1

Call: occuMulti(detformulas = detformulas, stateformulas = stateformulas, data = umf, method = “BFGS”, se = TRUE, engine = c(“C”), silent = FALSE, maxOrder = 2L)

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) 96.274 149.319 0.645 0.5191 [tigrinus] elev 80.638 126.032 0.640 0.5223 [perros] (Intercept) -2.324 1.159 -2.005 0.0449 [perros] elev -0.979 0.774 -1.265 0.2060

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -5.15 0.278 -18.51 1.75e-76 [perros] (Intercept) -4.39 0.961 -4.57 4.84e-06

AIC: 227.2404

codigo R

# update model
# occFormulas2 <- c('~dist_casa', '~dist_casa', '~dist_casa', "~1", "~1", "~1", "~0")
occFormulas2 <- c('~dist_casa', '~dist_casa', "~0")
fit2 <- update(fit1, stateformulas=occFormulas2)
fit2

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) 15.69 124.97 0.126 0.900 [tigrinus] dist_casa 15.99 129.90 0.123 0.902 [perros] (Intercept) -4.09 2.98 -1.372 0.170 [perros] dist_casa -3.38 3.61 -0.934 0.350

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -5.18 0.305 -17.02 6.22e-65 [perros] (Intercept) -4.39 0.971 -4.52 6.19e-06

AIC: 228.9213

Model Selection

codigo R

#List of fitted models
fl <- fitList(elev=fit1, dist=fit2)
# coef(fl)

###################
# Model selection #
###################

modSel(fl)

 nPars    AIC delta AICwt cumltvWt

elev 6 227.24 0.00 0.70 0.70 dist 6 228.92 1.68 0.30 1.00

codigo R

#############
# Model fit #
#############

# bt <- parboot(fit1, nsim=50) # takes time best model
# plot(bt)

plot predicted marginal occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_huila_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_huila_sf$elev))/sd(cams_huila_sf$elev)

r2 <- range(cams_huila_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_huila_sf$dist_casa))/sd(cams_huila_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)



tigrinus_pred <- predict(fit1, "state", species="tigrinus", newdata=nd)
tigrinus_pred$Species <- "tigrinus"

perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Elevation", y="Marginal occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted co-occurrence occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_huila_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_huila_sf$elev))/sd(cams_huila_sf$elev)

r2 <- range(cams_huila_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_huila_sf$dist_casa))/sd(cams_huila_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)


tigrinus_perros_pred <- predict(fit1, "state", 
                         species=c("tigrinus", "perros"), newdata=nd)

# tigrinus_pred$Species <- c("tigrinus", "perros")

# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- tigrinus_perros_pred #rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = "grey50"), alpha=0.3) +
  geom_line(aes(y=Predicted), col="blue") +
  labs(x="Elevation", y="co-occurence") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted conditional occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_huila_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_huila_sf$elev))/sd(cams_huila_sf$elev)

r2 <- range(cams_huila_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_huila_sf$dist_casa))/sd(cams_huila_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)

##### conditional

tigrinus_perro_no <- predict(fit1, "state", 
                         species="tigrinus", 
                         cond='-perros',
                         newdata=nd)

tigrinus_perro_no$Species <- "perro ausente"

tigrinus_perro_si <- predict(fit1, "state", 
                         species="tigrinus", 
                         cond='perros',
                         newdata=nd)

tigrinus_perro_si$Species <- "perro presente"





# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_perro_si, tigrinus_perro_no)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Elevation", y="tigrinus conditional occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

LAGUNA DE LA COCHA

codigo R

library(readr)
cocha1_data <- read_excel("C:/CodigoR/tigrinus2/data/Full_data_Ucumari_Huila_Cocha1_Cocha2.xlsx", 
    sheet = "cocha1")

cocha2_data <- read_excel("C:/CodigoR/tigrinus2/data/Full_data_Ucumari_Huila_Cocha1_Cocha2.xlsx", 
    sheet = "cocha2")

cocha2_data$binomial <- str_c (cocha2_data$Genus, "_", cocha2_data$Species)

crerar matrices para unmarked

codigo R

############# start spatial part
#### make sf object
cocha1_sf <- st_as_sf(cocha1_data, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")
#### make sf object
cocha2_sf <- st_as_sf(cocha2_data, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")


# get centroid. 1st bbox, make pol
centroid_cocha <- st_centroid(st_as_sfc(st_bbox(cocha1_sf)))
# get altitude
elev_cocha_full <- elevation_3s(centroid_cocha[[1]][1], centroid_cocha[[1]][2], 
                         path="C:/CodigoR/tigrinus2/raster")

# elev_ucu_full_ras <- elev_ucu_full %>% raster()

# st_bbox(ucumari_sf) # notice order xmin, xmax, ymin, ymax
ext_cocha <- ext(-77.17, -77.05,  0.886,    1.085 )
elev_cocha <- crop(elev_cocha_full, ext_cocha) 

# convert from terra to raster
elev_cocha_ras <-  elev_cocha %>% raster()

##### get uniques
cams_cocha1 <-cocha1_data %>% dplyr::select(c("Longitude", "Latitude", "camera_trap")) %>% distinct()
#### make sf object
cams_cocha1_sf <- st_as_sf(cams_cocha1, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")

cams_cocha2 <-cocha2_data %>% dplyr::select(c("Longitude", "Latitude", "camera_trap")) %>% distinct()
#### make sf object
cams_cocha2_sf <- st_as_sf(cams_cocha2, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")



# extract values from raster using altitude
cams_cocha1_sf$elev <- extract(elev_cocha_ras, cams_cocha1_sf)
cams_cocha2_sf$elev <- extract(elev_cocha_ras, cams_cocha2_sf)


############## make distance map
# # Convert points to sp spatialpointdatafram
# casas_points <- as(casas_sf, "Spatial")
# # Projection
# # Be sure to have rgdal first installed.
# casas_points_utm <- spTransform(casas_points, CRS('+init=epsg:32718'))
# # convert to ppp
# casas_points_ppp <- as(as(casas_points_utm, "SpatialPoints"), "ppp")
# # distance
# casas_distance <- distmap(casas_points_ppp)

####### extract distance 
# 
# casas_distance_ras<- raster(casas_distance) # convert raster
# crs(casas_distance_ras) <- '+init=epsg:32718' # add crs
# 
# # project raster
# casa_distance <- projectRaster(casas_distance_ras, crs = crs(casas_points))


# extracr raster values
# cams_ucu_sf$dist_casa <- raster::extract(casas_distance_ras, cams_ucu_sf)
cams_cocha1_sf$dist_casa <- raster::extract(casa_distance, cams_cocha1_sf) # also works
cams_cocha2_sf$dist_casa <- raster::extract(casa_distance, cams_cocha2_sf) # also works


# plot map
mapview(elev_cocha_ras) + mapview(cams_cocha1_sf["camera_trap"]) + mapview(cams_cocha2_sf["camera_trap"])

codigo R

############### end spatial part

# fix cocha2 data
# cocha2_data$binomial <- str_c (cocha2_data$Genus, "_", cocha2_data$Species)

#funcion para crear todas las tablas de datos
all_data_cocha1 <-  f.matrix.creator2 (cocha1_data)
all_data_cocha2 <-  f.matrix.creator2 (cocha2_data)
# names(all_data) # ver lass especies y en que lista esta cada una


# kable(names(all_data_cocha1)) # html table
# Tigrinus es lista 11 y perro 2

# kable(names(all_data_cocha2)) # html table
# Tigrinus es 30
# Perro no hay

datatable(
data = as.data.frame(names(all_data_cocha1)),
caption = "Especies La Cocha",
filter = "top"
)

unmarked

codigo R

# tabla con solo tiginus
tigrinus_cocha1 <- all_data_cocha1[[11]]
tigrinus_cocha2 <- all_data_cocha2[[30]]
# cargar paquete
library(unmarked)
# crear objeto umf
umf_tigrinus_cocha1 <- unmarkedFrameOccu(y=tigrinus_cocha1)
umf_tigrinus_cocha2 <- unmarkedFrameOccu(y=tigrinus_cocha2)
# verificar datos en grafica
plot(umf_tigrinus_cocha1)

codigo R

plot(umf_tigrinus_cocha2)

codigo R

# tabla con solo perro
perros_cocha1 <- all_data_cocha1[[2]]
# cargar paquete
library(unmarked)
# crear objeto umf
umf_perros_cocha1 <- unmarkedFrameOccu(y=perros_cocha1)
# verificar datos en grafica
plot(umf_perros_cocha1)

Modelo nulo

codigo R

# modelo nulo
fm_tig_cocha1 <- occu(~1 ~1, umf_tigrinus_cocha1)  # fit a model
fm_tig_cocha2 <- occu(~1 ~1, umf_tigrinus_cocha2)  # fit a model

backTransform(fm_tig_cocha1, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.0171 0.0092 -4.05 1

Transformation: logistic

codigo R

backTransform(fm_tig_cocha1, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.263 0.131 -1.03 1

Transformation: logistic

codigo R

backTransform(fm_tig_cocha2, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.0112 0.0052 -4.48 1

Transformation: logistic

codigo R

backTransform(fm_tig_cocha2, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.236 0.105 -1.17 1

Transformation: logistic

codigo R

# modelo nulo
fm_perros_cocha1 <- occu(~1 ~1, umf_perros_cocha1)  # fit a model

backTransform(fm_perros_cocha1, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.00179 0.000896 -6.32 1

Transformation: logistic

codigo R

backTransform(fm_perros_cocha1, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 1 NaN 7.83 1

Transformation: logistic

MODELOS DE CO-OCURRENCIA

codigo R

detformulas <- c( "~1", "~1")#, "~1")
#stateformulas <- c('~elev','~elev', '~elev', "~1", "~1", "~1", "~0")# 3 sp
stateformulas <- c('~elev','~elev', "~0")
y <- list(tigrinus_cocha1, perros_cocha1)# , ocelote_ucu)
names(y) <- c("tigrinus", "perros")#, "ocelote")


obs_covs <-as.data.frame(scale(cams_cocha1_sf$dist_casa))
names(obs_covs) <- "dist_casa"

site_covs_cocha1 <- data.frame(cams_cocha1_sf[,c('elev','dist_casa')])[,1:2]
site_covs_cocha1 <-as.data.frame(apply(site_covs_cocha1,2,scale))
names(site_covs_cocha1) <- c("elev", "dist_casa")


umf <-  unmarkedFrameOccuMulti(y=y, 
                              siteCovs=site_covs_cocha1,
                              obsCovs = NULL)
plot(umf)

codigo R

#umf

# occFormulas Length should match number/order of columns in fDesign
umf@fDesign

    f1[tigrinus] f2[perros] f3[tigrinus:perros]

psi[11] 1 1 1 psi[10] 1 0 0 psi[01] 0 1 0 psi[00] 0 0 0

codigo R

fit1 <- occuMulti(detformulas, stateformulas, umf,    
        method="BFGS", se=TRUE, engine=c("C"), silent=FALSE)

fit1

Call: occuMulti(detformulas = detformulas, stateformulas = stateformulas, data = umf, method = “BFGS”, se = TRUE, engine = c(“C”), silent = FALSE, maxOrder = 2L)

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) -1.019 0.695 -1.466 0.1426 [tigrinus] elev -0.254 0.523 -0.486 0.6271 [perros] (Intercept) -2.349 1.070 -2.194 0.0282 [perros] elev -0.961 0.917 -1.048 0.2946

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -4.06 0.553 -7.35 2.04e-13 [perros] (Intercept) -4.05 0.948 -4.28 1.90e-05

AIC: 189.2644

codigo R

# update model
# occFormulas2 <- c('~dist_casa', '~dist_casa', '~dist_casa', "~1", "~1", "~1", "~0")
occFormulas2 <- c('~dist_casa', '~dist_casa', "~0")
fit2 <- update(fit1, stateformulas=occFormulas2)
fit2

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) -1.0281 0.678 -1.516 0.1295 [tigrinus] dist_casa -0.0526 0.472 -0.111 0.9113 [perros] (Intercept) -2.1771 0.969 -2.247 0.0247 [perros] dist_casa -0.5529 0.865 -0.639 0.5226

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -4.05 0.546 -7.41 1.23e-13 [perros] (Intercept) -4.03 0.930 -4.33 1.47e-05

AIC: 190.4627

Model Selection

codigo R

#List of fitted models
fl <- fitList(elev=fit1, dist=fit2)
# coef(fl)


###################
# Model selection #
###################
modSel(fl)

 nPars    AIC delta AICwt cumltvWt

elev 6 189.26 0.00 0.65 0.65 dist 6 190.46 1.20 0.35 1.00

codigo R

#############
# Model fit #
#############
# bt <- parboot(fit2, nsim=50) # takes time
# plot(bt)

plot predicted marginal occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_cocha1_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_cocha1_sf$elev))/sd(cams_cocha1_sf$elev)

r2 <- range(cams_cocha1_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_cocha1_sf$dist_casa))/sd(cams_cocha1_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)



tigrinus_pred <- predict(fit2, "state", species="tigrinus", newdata=nd)
tigrinus_pred$Species <- "tigrinus"

perros_pred <- predict(fit2, "state", species="perros", newdata=nd)
perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x2,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Distancia casa", y="Marginal occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted co-occurrence occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_cocha1_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_cocha1_sf$elev))/sd(cams_cocha1_sf$elev)

r2 <- range(cams_ucu_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_cocha1_sf$dist_casa))/sd(cams_cocha1_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)


tigrinus_perros_pred <- predict(fit2, "state", 
                         species=c("tigrinus", "perros"), newdata=nd)

# tigrinus_pred$Species <- c("tigrinus", "perros")

# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- tigrinus_perros_pred #rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = "grey50"), alpha=0.3) +
  geom_line(aes(y=Predicted), col="blue") +
  labs(x="Distancia casa", y="co-occurence") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted conditional occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_cocha1_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_cocha1_sf$elev))/sd(cams_cocha1_sf$elev)

r2 <- range(cams_cocha1_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_cocha1_sf$dist_casa))/sd(cams_cocha1_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)

##### conditional

tigrinus_perro_no <- predict(fit2, "state", 
                         species="tigrinus", 
                         cond='-perros',
                         newdata=nd)

tigrinus_perro_no$Species <- "perro ausente"

tigrinus_perro_si <- predict(fit2, "state", 
                         species="tigrinus", 
                         cond='perros',
                         newdata=nd)

tigrinus_perro_si$Species <- "perro presente"





# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_perro_si, tigrinus_perro_no)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x2,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Distancia casa", y="tigrinus conditional occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

LA FE

codigo R

#seleccionar solo 2021

lafe_data_raw <- read_excel("C:/CodigoR/tigrinus2/data/cuencaverde.xlsx", 
    sheet = "LaFe_completo")

library(lubridate)
lafe_data_raw$year <- year(lafe_data_raw$Photo_Date)

# flter by 2021
lafe_data <- lafe_data_raw %>% filter(year >= "2021")

#huila_data$binomial <- str_c (huila_data$Genus, "_", huila_data$Species)

crerar matrices para unmarked

codigo R

############# start spatial part
#### make sf object
lafe_sf <- st_as_sf(lafe_data, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")


# get centroid. 1st bbox, make pol
centroid_lafe <- st_centroid(st_as_sfc(st_bbox(lafe_sf)))
# get altitude
elev_lafe_full <- elevation_3s(centroid_lafe[[1]][1], centroid_lafe[[1]][2], 
                         path="C:/CodigoR/tigrinus2/raster")

# elev_ucu_full_ras <- elev_ucu_full %>% raster()

# st_bbox(ucumari_sf) # notice order xmin, xmax, ymin, ymax
ext_lafe <- ext(-75.49, -75.179,  5.81,    6.03 )
elev_lafe <- crop(elev_lafe_full, ext_lafe) 

# convert from terra to raster
elev_lafe_ras <-  elev_lafe %>% raster()

# get uniques
cams_lafe <-lafe_data %>% dplyr::select(c("Longitude", "Latitude", "camera_trap")) %>% distinct()
#### make sf object
cams_lafe_sf <- st_as_sf(cams_lafe, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")

# extract values from raster using altitude
cams_lafe_sf$elev <- extract(elev_lafe_ras, cams_lafe_sf)

############## make distance map
# # Convert points to sp spatialpointdatafram
# casas_points <- as(casas_sf, "Spatial")
# # Projection
# # Be sure to have rgdal first installed.
# casas_points_utm <- spTransform(casas_points, CRS('+init=epsg:32718'))
# # convert to ppp
# casas_points_ppp <- as(as(casas_points_utm, "SpatialPoints"), "ppp")
# # distance
# casas_distance <- distmap(casas_points_ppp)

####### extract distance 
# 
# casas_distance_ras<- raster(casas_distance) # convert raster
# crs(casas_distance_ras) <- '+init=epsg:32718' # add crs
# 
# # project raster
# casa_distance <- projectRaster(casas_distance_ras, crs = crs(casas_points))

# cams_ucu_sf$dist_casa <- raster::extract(casas_distance_ras, cams_ucu_sf)
cams_lafe_sf$dist_casa <- raster::extract(casa_distance, cams_lafe_sf) # also works


# plot map
mapview(elev_lafe_ras) + mapview(cams_lafe_sf["camera_trap"])

codigo R

############### end spatial part

#lafe_data$binomial <- str_c (lafe_data$Genus, "_", lafe_data$Species)

#funcion para crear todas las tablas de datos
all_data_lafe <-  f.matrix.creator2 (lafe_data)

# names(all_data) # ver lass especies y en que lista esta cada una
# kable(names(all_data)) # html table
# Tigrinus es lista 3
# perro es 2

datatable(
data = as.data.frame(names(all_data_lafe)),
caption = "Especies La Fe",
filter = "top"
)

unmarked

codigo R

# tabla con solo tiginus
tigrinus_lafe <- all_data_lafe[[3]]

# colapsa a una semana take time
tigrinus_lafe_48<-f.shrink.matrix.to26(matrix = all_data_lafe[[3]])

# crear objeto umf
umf_tigrinus_lafe <- unmarkedFrameOccu(y=tigrinus_lafe_48)
# verificar datos en grafica
# plot(umf_tigrinus_ucu)

# tabla con solo perros
perros_lafe <- all_data_lafe[[2]]

# colapsa a una semana take time
perros_lafe_48<-f.shrink.matrix.to26(matrix = all_data_lafe[[2]])
# crear objeto umf
umf_perros_lafe <- unmarkedFrameOccu(y=perros_lafe_48)
# verificar datos en grafica
# plot(umf_perros_ucu)

# # tabla con solo ocelote
# ocelote_lafe <- all_data[[16]]
# # crear objeto umf
# umf_ocelote_lafe <- unmarkedFrameOccu(y=ocelote_ucu)
# # verificar datos en grafica
# # plot(umf_ocelote_ucu)

Modelo nulo

codigo R

# modelo nulo tigrinus
fm_tig_lafe <- occu(~1 ~1, umf_tigrinus_lafe)  # fit a model

backTransform(fm_tig_lafe, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.129 0.0797 -1.91 1

Transformation: logistic

codigo R

backTransform(fm_tig_lafe, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.513 0.293 0.0508 1

Transformation: logistic

codigo R

# modelo nulo perro
fm_perros_lafe <- occu(~1 ~1, umf_perros_lafe)  # fit a model

backTransform(fm_perros_lafe, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.247 0.06 -1.11 1

Transformation: logistic

codigo R

backTransform(fm_perros_lafe, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.748 0.171 1.09 1

Transformation: logistic

codigo R

# # modelo nulo ocelote
# fm_ocelote_lafe <- occu(~1 ~1, umf_ocelote_lafe)  # fit a model
# 
# backTransform(fm_ocelote_lafe, type="det") #estimado lineal de deteccion
# backTransform(fm_ocelote_lafe, type="state") # estimado linel de ocupacion
#

MODELOS DE CO-OCURRENCIA

codigo R

detformulas <- c( "~1", "~1")#, "~1")
#stateformulas <- c('~elev','~elev', '~elev', "~1", "~1", "~1", "~0")# 3 sp
stateformulas <- c('~elev','~elev', "~0")
y <- list(tigrinus_lafe, perros_lafe)# , ocelote_lafe)
names(y) <- c("tigrinus", "perros")#, "ocelote")


obs_covs <-as.data.frame(scale(cams_lafe_sf$dist_casa))
names(obs_covs) <- "dist_casa"

site_covs_lafe <- data.frame(cams_lafe_sf[,c('elev','dist_casa')])[,1:2]
site_covs_lafe <-as.data.frame(apply(site_covs_lafe,2,scale))
names(site_covs_lafe) <- c("elev", "dist_casa")


umf <-  unmarkedFrameOccuMulti(y=y, 
                              siteCovs=site_covs_lafe,
                              obsCovs = NULL)
plot(umf)

codigo R

#umf

# occFormulas Length should match number/order of columns in fDesign
umf@fDesign

    f1[tigrinus] f2[perros] f3[tigrinus:perros]

psi[11] 1 1 1 psi[10] 1 0 0 psi[01] 0 1 0 psi[00] 0 0 0

codigo R

fit1 <- occuMulti(detformulas, stateformulas, umf,     
        method="BFGS", se=TRUE, engine=c("C"), silent=FALSE)

fit1

Call: occuMulti(detformulas = detformulas, stateformulas = stateformulas, data = umf, method = “BFGS”, se = TRUE, engine = c(“C”), silent = FALSE, maxOrder = 2L)

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) -0.362 0.738 -0.491 0.624 [tigrinus] elev 0.179 0.589 0.304 0.761 [perros] (Intercept) 12.205 NaN NaN NaN [perros] elev -2.326 NaN NaN NaN

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -3.65 0.472 -7.73 1.06e-14 [perros] (Intercept) -3.30 0.167 -19.69 2.46e-86

AIC: 449.8858

codigo R

# update model
# occFormulas2 <- c('~dist_casa', '~dist_casa', '~dist_casa', "~1", "~1", "~1", "~0")
occFormulas2 <- c('~dist_casa', '~dist_casa', "~1")
fit2 <- update(fit1, stateformulas=occFormulas2)
fit2

Call: occuMulti(detformulas = c(“~1”, “~1”), stateformulas = c(“~dist_casa”, “~dist_casa”, “~1”), data = umf, maxOrder = 2L, method = “BFGS”, se = TRUE, engine = c(“C”), silent = FALSE)

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) -2.991 NaN NaN NaN [tigrinus] dist_casa -0.545 0.693 -0.786 0.432 [perros] (Intercept) 12.387 NaN NaN NaN [perros] dist_casa 2.229 NaN NaN NaN [tigrinus:perros] (Intercept) 2.674 NaN NaN NaN

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -3.69 0.488 -7.56 3.99e-14 [perros] (Intercept) -3.30 0.167 -19.69 2.48e-86

AIC: 451.2344

Model Selection

codigo R

#List of fitted models
fl <- fitList(elev=fit1, dist=fit2)
# coef(fl)

###################
# Model selection #
###################

modSel(fl)

 nPars    AIC delta AICwt cumltvWt

elev 6 449.89 0.00 0.66 0.66 dist 7 451.23 1.35 0.34 1.00

codigo R

#############
# Model fit #
#############

# bt <- parboot(fit1, nsim=50) # takes time best model
# plot(bt)

plot predicted marginal occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_lafe_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_lafe_sf$elev))/sd(cams_lafe_sf$elev)

r2 <- range(cams_lafe_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_lafe_sf$dist_casa))/sd(cams_lafe_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)



tigrinus_pred <- predict(fit1, "state", species="tigrinus", newdata=nd)
tigrinus_pred$Species <- "tigrinus"

perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Elev", y="Marginal occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted co-occurrence occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_huila_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_huila_sf$elev))/sd(cams_huila_sf$elev)

r2 <- range(cams_huila_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_huila_sf$dist_casa))/sd(cams_huila_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)


tigrinus_perros_pred <- predict(fit1, "state", 
                         species=c("tigrinus", "perros"), newdata=nd)

# tigrinus_pred$Species <- c("tigrinus", "perros")

# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- tigrinus_perros_pred #rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = "grey50"), alpha=0.3) +
  geom_line(aes(y=Predicted), col="blue") +
  labs(x="Elev", y="co-occurence") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted conditional occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_lafe_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_lafe_sf$elev))/sd(cams_lafe_sf$elev)

r2 <- range(cams_lafe_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_lafe_sf$dist_casa))/sd(cams_lafe_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)

##### conditional

tigrinus_perro_no <- predict(fit1, "state", 
                         species="tigrinus", 
                         cond='-perros',
                         newdata=nd)

tigrinus_perro_no$Species <- "perro ausente"

tigrinus_perro_si <- predict(fit1, "state", 
                         species="tigrinus", 
                         cond='perros',
                         newdata=nd)

tigrinus_perro_si$Species <- "perro presente"





# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_perro_si, tigrinus_perro_no)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Distancia Casa", y="tigrinus conditional occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

RIO GRANDE

codigo R

riogrande_data <- read_excel("C:/CodigoR/tigrinus2/data/cuencaverde.xlsx", 
    sheet = "Riogrande_completo")



#huila_data$binomial <- str_c (huila_data$Genus, "_", huila_data$Species)

crerar matrices para unmarked

codigo R

############# start spatial part
#### make sf object
riogrande_sf <- st_as_sf(riogrande_data, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")


# get centroid. 1st bbox, make pol
centroid_riogrande <- st_centroid(st_as_sfc(st_bbox(riogrande_sf)))
# get altitude
elev_riogrande_full <- elevation_3s(centroid_riogrande[[1]][1], centroid_riogrande[[1]][2], 
                         path="C:/CodigoR/tigrinus2/raster")

# elev_ucu_full_ras <- elev_ucu_full %>% raster()

# st_bbox(ucumari_sf) # notice order xmin, xmax, ymin, ymax
ext_riogrande <- ext(-75.94, -75.275,  6.27,    6.896 )
elev_riogrande <- crop(elev_riogrande_full, ext_riogrande) 

# convert from terra to raster
elev_riogrande_ras <-  elev_riogrande %>% raster()

# get uniques
cams_riogrande <-riogrande_data %>% dplyr::select(c("Longitude", "Latitude", "camera_trap")) %>% distinct()
#### make sf object
cams_riogrande_sf <- st_as_sf(cams_riogrande, coords = c("Longitude", "Latitude"), crs = "EPSG:4326")

# extract values from raster using altitude
cams_riogrande_sf$elev <- extract(elev_riogrande_ras, cams_riogrande_sf)

############## make distance map
# # Convert points to sp spatialpointdatafram
# casas_points <- as(casas_sf, "Spatial")
# # Projection
# # Be sure to have rgdal first installed.
# casas_points_utm <- spTransform(casas_points, CRS('+init=epsg:32718'))
# # convert to ppp
# casas_points_ppp <- as(as(casas_points_utm, "SpatialPoints"), "ppp")
# # distance
# casas_distance <- distmap(casas_points_ppp)

####### extract distance 
# 
# casas_distance_ras<- raster(casas_distance) # convert raster
# crs(casas_distance_ras) <- '+init=epsg:32718' # add crs
# 
# # project raster
# casa_distance <- projectRaster(casas_distance_ras, crs = crs(casas_points))

# cams_ucu_sf$dist_casa <- raster::extract(casas_distance_ras, cams_ucu_sf)
cams_riogrande_sf$dist_casa <- raster::extract(casa_distance, cams_riogrande_sf) # also works


# plot map
mapview(elev_riogrande_ras) + mapview(cams_riogrande_sf["camera_trap"])

codigo R

############### end spatial part

#riogrande_data$binomial <- str_c (riogrande_data$Genus, "_", riogrande_data$Species)

#funcion para crear todas las tablas de datos
all_data_riogrande <-  f.matrix.creator2 (riogrande_data)

# names(all_data) # ver lass especies y en que lista esta cada una
# kable(names(all_data)) # html table
# Tigrinus es lista 11
# perro es 1

datatable(
data = as.data.frame(names(all_data_riogrande)),
caption = "Especies Rio Grande",
filter = "top"
)

unmarked

codigo R

# colapsa a una semana take time
tigrinus_riogrande_48<-f.shrink.matrix.to26(matrix = all_data_riogrande[[11]])

# tabla con solo tiginus
tigrinus_riogrande <- tigrinus_riogrande_48 #all_data_riogrande[[11]]



# crear objeto umf
umf_tigrinus_riogrande <- unmarkedFrameOccu(y=tigrinus_riogrande)
# verificar datos en grafica
# plot(umf_tigrinus_ucu)

# colapsa a una semana take time
perros_riogrande_48<-f.shrink.matrix.to26(matrix = all_data_riogrande[[1]])

# tabla con solo perros
perros_riogrande <- perros_riogrande_48 #all_data_riogrande[[1]]
# crear objeto umf
umf_perros_riogrande <- unmarkedFrameOccu(y=perros_riogrande)
# verificar datos en grafica
# plot(umf_perros_ucu)

# # tabla con solo ocelote
# ocelote_riogrande <- all_data[[16]]
# # crear objeto umf
# umf_ocelote_riogrande <- unmarkedFrameOccu(y=ocelote_ucu)
# # verificar datos en grafica
# # plot(umf_ocelote_ucu)

Modelo nulo

codigo R

# modelo nulo tigrinus
fm_tig_riogrande <- occu(~1 ~1, umf_tigrinus_riogrande)  # fit a model

backTransform(fm_tig_riogrande, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.225 0.102 -1.24 1

Transformation: logistic

codigo R

backTransform(fm_tig_riogrande, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.236 0.111 -1.18 1

Transformation: logistic

codigo R

# modelo nulo perro
fm_perros_riogrande <- occu(~1 ~1, umf_perros_riogrande)  # fit a model

backTransform(fm_perros_riogrande, type="det") #estimado lineal de deteccion

Backtransformed linear combination(s) of Detection estimate(s)

Estimate SE LinComb (Intercept) 0.118 0.0936 -2.01 1

Transformation: logistic

codigo R

backTransform(fm_perros_riogrande, type="state") # estimado linel de ocupacion

Backtransformed linear combination(s) of Occupancy estimate(s)

Estimate SE LinComb (Intercept) 0.482 0.33 -0.0723 1

Transformation: logistic

codigo R

# # modelo nulo ocelote
# fm_ocelote_riogrande <- occu(~1 ~1, umf_ocelote_riogrande)  # fit a model
# 
# backTransform(fm_ocelote_riogrande, type="det") #estimado lineal de deteccion
# backTransform(fm_ocelote_riogrande, type="state") # estimado linel de ocupacion
#

MODELOS DE CO-OCURRENCIA

codigo R

detformulas <- c( "~1", "~1")#, "~1")
#stateformulas <- c('~elev','~elev', '~elev', "~1", "~1", "~1", "~0")# 3 sp
stateformulas <- c('~elev','~elev', "~0")
y <- list(tigrinus_riogrande, perros_riogrande)# , ocelote_riogrande)
names(y) <- c("tigrinus", "perros")#, "ocelote")


obs_covs <-as.data.frame(scale(cams_riogrande_sf$dist_casa))
names(obs_covs) <- "dist_casa"

site_covs_riogrande <- data.frame(cams_riogrande_sf[,c('elev','dist_casa')])[,1:2]
site_covs_riogrande <-as.data.frame(apply(site_covs_riogrande,2,scale))
names(site_covs_riogrande) <- c("elev", "dist_casa")


umf <-  unmarkedFrameOccuMulti(y=y, 
                              siteCovs=site_covs_riogrande,
                              obsCovs = NULL)
plot(umf)

codigo R

#umf

# occFormulas Length should match number/order of columns in fDesign
umf@fDesign

    f1[tigrinus] f2[perros] f3[tigrinus:perros]

psi[11] 1 1 1 psi[10] 1 0 0 psi[01] 0 1 0 psi[00] 0 0 0

codigo R

fit1 <- occuMulti(detformulas, stateformulas, umf,    
        method="BFGS", se=TRUE, engine=c("C"), silent=FALSE)

fit1

Call: occuMulti(detformulas = detformulas, stateformulas = stateformulas, data = umf, method = “BFGS”, se = TRUE, engine = c(“C”), silent = FALSE, maxOrder = 2L)

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) -0.888 0.763 -1.164 0.245 [tigrinus] elev -2.289 1.549 -1.478 0.140 [perros] (Intercept) -0.498 0.861 -0.578 0.564 [perros] elev -0.579 0.655 -0.884 0.377

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -1.48 0.592 -2.49 0.0126 [perros] (Intercept) -1.70 0.713 -2.39 0.0170

AIC: 143.037

codigo R

# update model
# occFormulas2 <- c('~dist_casa', '~dist_casa', '~dist_casa', "~1", "~1", "~1", "~0")
occFormulas2 <- c('~dist_casa', '~dist_casa', "~0")
fit2 <- update(fit1, stateformulas=occFormulas2)
fit2

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) -1.211 0.632 -1.916 0.0554 [tigrinus] dist_casa -0.364 1.029 -0.354 0.7237 [perros] (Intercept) -0.315 1.114 -0.283 0.7772 [perros] dist_casa 1.859 2.907 0.639 0.5225

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -1.23 0.584 -2.11 0.0349 [perros] (Intercept) -1.69 0.721 -2.34 0.0190

AIC: 143.4188

Model Selection

codigo R

#List of fitted models
fl <- fitList(elev=fit1, dist=fit2)
# coef(fl)

###################
# Model selection #
###################

modSel(fl)

 nPars    AIC delta AICwt cumltvWt

elev 6 143.04 0.00 0.55 0.55 dist 6 143.42 0.38 0.45 1.00

codigo R

#############
# Model fit #
#############

# bt <- parboot(fit1, nsim=50) # takes time best model
# plot(bt)

plot predicted marginal occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_riogrande_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_riogrande_sf$elev))/sd(cams_riogrande_sf$elev)

r2 <- range(cams_riogrande_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_riogrande_sf$dist_casa))/sd(cams_riogrande_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)



tigrinus_pred <- predict(fit2, "state", species="tigrinus", newdata=nd)
tigrinus_pred$Species <- "tigrinus"

perros_pred <- predict(fit2, "state", species="perros", newdata=nd)
perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x2,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Distancia Casa", y="Marginal occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted co-occurrence occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_lafe_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_lafe_sf$elev))/sd(cams_lafe_sf$elev)

r2 <- range(cams_lafe_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_lafe_sf$dist_casa))/sd(cams_lafe_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)


tigrinus_perros_pred <- predict(fit2, "state", 
                         species=c("tigrinus", "perros"), newdata=nd)

# tigrinus_pred$Species <- c("tigrinus", "perros")

# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- tigrinus_perros_pred #rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = "grey50"), alpha=0.3) +
  geom_line(aes(y=Predicted), col="blue") +
  labs(x="Distancia Casa", y="co-occurence") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted conditional occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_riogrande_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_riogrande_sf$elev))/sd(cams_riogrande_sf$elev)

r2 <- range(cams_riogrande_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_riogrande_sf$dist_casa))/sd(cams_riogrande_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)

##### conditional

tigrinus_perro_no <- predict(fit2, "state", 
                         species="tigrinus", 
                         cond='-perros',
                         newdata=nd)

tigrinus_perro_no$Species <- "perro ausente"

tigrinus_perro_si <- predict(fit2, "state", 
                         species="tigrinus", 
                         cond='perros',
                         newdata=nd)

tigrinus_perro_si$Species <- "perro presente"





# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_perro_si, tigrinus_perro_no)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x1,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Distancia Casa", y="tigrinus conditional occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

TODOS JUNTOS COLAPSADOS A 26

26 equivale a cada 15 dias para la fe y rio grande

codigo R

# colapsa a una semana take time

# Ucumari
tigrinus_ucu_26 <- f.shrink.matrix.to26(matrix = all_data_ucu[[8]])
perros_ucu_26 <- f.shrink.matrix.to26(matrix = all_data_ucu[[75]])
# Huila
tigrinus_huila_26 <- f.shrink.matrix.to26(matrix = all_data_huila[[6]])
perros_huila_26 <- f.shrink.matrix.to26(matrix = all_data_huila[[43]])
# Cocha1
tigrinus_cocha1_26 <- f.shrink.matrix.to26(matrix = all_data_cocha1[[11]])
perros_cocha1_26 <- f.shrink.matrix.to26(matrix = all_data_cocha1[[2]])
#La Fe
tigrinus_lafe_26<-f.shrink.matrix.to26(matrix = all_data_lafe[[3]])
perros_lafe_26<-f.shrink.matrix.to26(matrix = all_data_lafe[[3]])
# Rio Grande
tigrinus_riogrande_26<-f.shrink.matrix.to26(matrix = all_data_riogrande[[11]])
perros_riogrande_26<-f.shrink.matrix.to26(matrix = all_data_riogrande[[1]])


tigrinus_full_26 <- rbind(tigrinus_ucu_26, 
                          tigrinus_huila_26,
                          tigrinus_cocha1_26,
                          tigrinus_lafe_26,
                          tigrinus_riogrande_26)

perros_full_26 <- rbind(perros_ucu_26, 
                          perros_huila_26,
                          perros_cocha1_26,
                          perros_lafe_26,
                          perros_riogrande_26)

# join raw data
site_covs_join <- rbind(data.frame(cams_ucu_sf[,c('elev','dist_casa')])[,1:2],
                        data.frame(cams_huila_sf[,c('elev','dist_casa')])[,1:2],
                        data.frame(cams_cocha1_sf[,c('elev','dist_casa')])[,1:2],
                        data.frame(cams_lafe_sf[,c('elev','dist_casa')])[,1:2],
                        data.frame(cams_riogrande_sf[,c('elev','dist_casa')])[,1:2])

# standarize
site_covs_full <-as.data.frame(apply(site_covs_join,2,scale))
# put names
names(site_covs_full) <- c("elev", "dist_casa")

MODELOS DE CO-OCURRENCIA

codigo R

detformulas <- c( "~1", "~1")#, "~1")
#stateformulas <- c('~elev','~elev', '~elev', "~1", "~1", "~1", "~0")# 3 sp
stateformulas <- c('~elev','~elev', "~0")
y_full <- list(tigrinus_full_26, perros_full_26)# , ocelote_ucu)
names(y_full) <- c("tigrinus", "perros")#, "ocelote")


# obs_covs <-as.data.frame(scale(cams_ucu_sf$dist_casa))
# names(obs_covs) <- "dist_casa"


umf <-  unmarkedFrameOccuMulti(y=y_full, 
                              siteCovs=site_covs_full,
                              obsCovs = NULL)
plot(umf)

codigo R

#umf

# occFormulas Length should match number/order of columns in fDesign
umf@fDesign

    f1[tigrinus] f2[perros] f3[tigrinus:perros]

psi[11] 1 1 1 psi[10] 1 0 0 psi[01] 0 1 0 psi[00] 0 0 0

codigo R

fit1 <- occuMulti(detformulas, stateformulas, umf,    
        method="BFGS", se=TRUE, engine=c("C"), silent=FALSE)

fit1

Call: occuMulti(detformulas = detformulas, stateformulas = stateformulas, data = umf, method = “BFGS”, se = TRUE, engine = c(“C”), silent = FALSE, maxOrder = 2L)

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) 4.62 3.12 1.48054 0.139 [tigrinus] elev -3.49 2.19 -1.59028 0.112 [perros] (Intercept) 28.26 16777.22 0.00168 0.999 [perros] elev -1.54 3.17 -0.48613 0.627

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -3.72 0.164 -22.7 1.09e-113 [perros] (Intercept) -4.55 0.197 -23.1 1.32e-117

AIC: 775.9273

codigo R

# fit_opt <- optimizePenalty(fit1, penalties=c(0,1,2))

# update model
# occFormulas2 <- c('~dist_casa', '~dist_casa', '~dist_casa', "~1", "~1", "~1", "~0")
occFormulas2 <- c('~dist_casa', '~dist_casa', "~0")
fit2 <- update(fit1, stateformulas=occFormulas2)
fit2

Occupancy: Estimate SE z P(>|z|) [tigrinus] (Intercept) -0.7606 0.337 -2.254 2.42e-02 [tigrinus] dist_casa -0.2583 0.259 -0.998 3.18e-01 [perros] (Intercept) -1.7709 0.332 -5.335 9.55e-08 [perros] dist_casa -0.0737 0.239 -0.309 7.57e-01

Detection: Estimate SE z P(>|z|) [tigrinus] (Intercept) -2.65 0.279 -9.50 2.04e-21 [perros] (Intercept) -2.24 0.396 -5.65 1.63e-08

AIC: 752.2022

Model Selection

codigo R

#List of fitted models
fl <- fitList(elev=fit1, dist=fit2)
# coef(fl)


###################
# Model selection #
###################

modSel(fl)

 nPars    AIC delta AICwt cumltvWt

dist 6 752.20 0.00 1e+00 1.00 elev 6 775.93 23.73 7e-06 1.00

codigo R

#############
# Model fit #
#############

# bt <- parboot(fit1, nsim=50) # takes time best model
# plot(bt)

plot predicted marginal occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(site_covs_join$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(site_covs_join$elev))/sd(site_covs_join$elev)

r2 <- range(site_covs_join$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(site_covs_join$dist_casa))/sd(site_covs_join$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)



tigrinus_pred <- predict(fit2, "state", species="tigrinus", newdata=nd)
tigrinus_pred$Species <- "tigrinus"

perros_pred <- predict(fit2, "state", species="perros", newdata=nd)
perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x2,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Dist Casas", y="Marginal occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted co-occurrence occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_ucu_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_ucu_sf$elev))/sd(cams_ucu_sf$elev)

r2 <- range(cams_ucu_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_ucu_sf$dist_casa))/sd(cams_ucu_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)


tigrinus_perros_pred <- predict(fit2, "state", 
                         species=c("tigrinus", "perros"), newdata=nd)

# tigrinus_pred$Species <- c("tigrinus", "perros")

# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- tigrinus_perros_pred #rbind(tigrinus_pred, perros_pred)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill = "grey50"), alpha=0.3) +
  geom_line(aes(y=Predicted), col="blue") +
  labs(x="Dist Casa", y="co-occurence") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))

plot predicted conditional occupancy

codigo R

#Plot predicted marginal occupancy as a function of disturbance
r <- range(cams_ucu_sf$elev)
x1 <- seq(r[1],r[2],length.out=100)
x_scale <- (x1-mean(cams_ucu_sf$elev))/sd(cams_ucu_sf$elev)

r2 <- range(cams_ucu_sf$dist_casa)
x2 <- seq(r2[1],r2[2],length.out=100)
x2_scale <- (x2-mean(cams_ucu_sf$dist_casa))/sd(cams_ucu_sf$dist_casa)

nd <- matrix(NA, 100, 2)
nd <- data.frame(elev=x_scale, dist_casa= x2_scale)

##### conditional

tigrinus_perro_no <- predict(fit2, "state", 
                         species="tigrinus", 
                         cond='-perros',
                         newdata=nd)

tigrinus_perro_no$Species <- "perro ausente"

tigrinus_perro_si <- predict(fit2, "state", 
                         species="tigrinus", 
                         cond='perros',
                         newdata=nd)

tigrinus_perro_si$Species <- "perro presente"





# perros_pred <- predict(fit1, "state", species="perros", newdata=nd)
# perros_pred$Species <- "perros"

# ocelote_pred <- predict(fit2, "state", species="ocelote", newdata=nd)
# ocelote_pred$Species <- "ocelote"


plot_dat <- rbind(tigrinus_perro_si, tigrinus_perro_no)#, ocelote_pred)

ggplot(data=plot_dat, aes(x=rep(x2,2), y=Predicted)) + # change to 3 sp and x2 to distance
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=Species), alpha=0.3) +
  geom_line(aes(col=Species)) +
  labs(x="Dist Casa", y="tigrinus conditional occupancy") +
  theme_bw() +
  theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
        axis.text=element_text(size=12), axis.title=element_text(size=14),
        legend.text=element_text(size=12), legend.title=element_text(size=14))