Occupancy Panama

One species, one season

Author

Diego J. Lizcano, Jorge Velásquez-Tibata

Published

April 21, 2023

Code 
library(tidyverse)
library(ggplot2)
library(raster)
library(knitr)
library(terra)

Selected species

Camptostoma obsoletum

Threshold

Read data

Code 
library(readxl)
library(raster)

Camptostoma_obsoletum_occu_dia <- read_excel("C:/Audubon_Panama/result/occu/Camptostoma_obsoletum/Camptostoma_obsoletum.xlsx", 
    sheet = "dia", col_types = c("numeric", 
        "text", "numeric", "numeric", "text", "text", 
        "numeric", "date", "date", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric"))

Camptostoma_obsoletum_occu_hora <- read_excel("C:/Audubon_Panama/result/occu/Camptostoma_obsoletum/Camptostoma_obsoletum.xlsx", 
    sheet = "hora", col_types = c("numeric", 
        "text", "numeric", "numeric", "text", "text", 
        "numeric", "date", "date", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric"))


# centroide for terra
# centoide <- centroids(puntos, TRUE)
centroide <- c(mean(as.numeric(Camptostoma_obsoletum_occu_dia$Longitude)), mean(as.numeric(Camptostoma_obsoletum_occu_dia$Latitude)))
clip_window <- extent(-80.612233, -80.180109, 7.852311, 8.390073) # 7.932311, -80.612233  8.360073, -80.180109
bb <- c(-80.612233, -80.180109, 7.932311, 8.360073)
srtm <- raster::getData('SRTM', lon=centroide[1], lat=centroide[2])

# crop the  raster using the vector extent
srtm_crop <- raster::crop(srtm, clip_window)


#altitud <- elevation_3s(-72.893262, 7.664081007, path="data")
altitud <- as.numeric(Camptostoma_obsoletum_occu_dia$Altitude)

# convierte a terra
puntos <- vect(Camptostoma_obsoletum_occu_dia, geom=c("Longitude", "Latitude"), crs="EPSG:4326")
# convierte a sf
puntos_sf <- sf::st_as_sf(puntos)

cam_covs <- raster::extract(srtm_crop, puntos_sf)

Fit Occupancy models

Make unmarked frame

Code 
library(unmarked)

umf_y_full<- unmarkedFrameOccu(y= Camptostoma_obsoletum_occu_dia[,10:36])
siteCovs(umf_y_full) <- data.frame(elevation =  as.numeric(Camptostoma_obsoletum_occu_dia$Altitude)) # data.frame(Elev=full_covs$Elev) # Full
obsCovs(umf_y_full) <- list(hour= as.data.frame(Camptostoma_obsoletum_occu_hora[,10:36]))

#######Graficar umf
plot(umf_y_full)

Build Models

Null model

Code 
mf0<-occu(~1 ~ 1, umf_y_full)
backTransform(mf0, type="det")  
Backtransformed linear combination(s) of Detection estimate(s)

 Estimate     SE LinComb (Intercept)
    0.173 0.0382   -1.57           1

Transformation: logistic
Code 
backTransform(mf0, type="state")  
Backtransformed linear combination(s) of Occupancy estimate(s)

 Estimate    SE LinComb (Intercept)
    0.315 0.117  -0.775           1

Transformation: logistic

Model selection

Code 
mf0<-occu(~1 ~ 1, umf_y_full)
mf1<-occu(~1 ~ scale(elevation), umf_y_full)
mf2<-occu(~1 ~ scale(elevation) +I(scale(elevation)^2), umf_y_full)
# mf3<-occu(~hour ~ 1, umf_y_full)
# mf4<-occu(~1 ~ aspect, umf_y_full)
# mf5<-occu(~1 ~ roughness, umf_y_full, starts = c(0.6, -3, 0))
# mf6<-occu(~elevation +I(elevation^2) ~ elevation +I(elevation^2), umf_y_full)
# mf7<-occu(~roughness ~ elevation +I(elevation^2), umf_y_full)
# mf8<-occu(~slope ~ elevation +I(elevation^2), umf_y_full)

# fit list
fms1<-fitList("p(.) Ocu(.)"=mf0,
              "p(.) Ocu(elev)"=mf1,
              "p(.) Ocu(elev^2)"=mf2 )

modSel(fms1)
                 nPars    AIC delta  AICwt cumltvWt
p(.) Ocu(elev)       3 101.46  0.00 0.7300     0.73
p(.) Ocu(elev^2)     4 103.45  1.99 0.2699     1.00
p(.) Ocu(.)          2 119.23 17.77 0.0001     1.00

Best model

Code 
summary(mf1)

Call:
occu(formula = ~1 ~ scale(elevation), data = umf_y_full)

Occupancy (logit-scale):
                 Estimate   SE      z P(>|z|)
(Intercept)         -4.21 14.0 -0.300   0.764
scale(elevation)    13.29 30.1  0.441   0.659

Detection (logit-scale):
 Estimate   SE     z  P(>|z|)
    -1.55 0.26 -5.98 2.18e-09

AIC: 101.461 
Number of sites: 16
optim convergence code: 0
optim iterations: 45 
Bootstrap iterations: 0

Model fit assesment

Code 
pb_f <- parboot(mf1, nsim=500, report=10) 
t0 = 14.85526 
Running parametric bootstrap in parallel on 3 cores.
Bootstrapped statistics not reported during parallel processing.
Code 
## t0 = 36.31634
plot (pb_f)

Prediction using “best” model

Code 
newdat_range<-data.frame(elevation=seq(-1,
                                       max(as.vector(srtm_crop), na.rm = T),length=100))#, 
#                         roughness=seq(min(full_covs_s$roughness))#,
#                                       max(full_covs_s$roughness), length=100))
# ## plot Detection en escala original
# pred_det <-predict(mf6, type="det", newdata=newdat_range, appendData=TRUE)
# plot(Predicted~elevation, pred_det,type="l",col="blue", 
#      xlab="Elevation",
#      ylab="Detection Probability",
#      xaxt="n")
# xticks <- c(-1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3) # -1:2
# xlabs <- xticks*sd(full_covs$roughness) + mean(full_covs$roughness) #Use the mean and sd of original value to change label name
# axis(1, at=xticks, labels=round(xlabs, 1))
# lines(lower~elevation, pred_det,type="l",col=gray(0.5))
# lines(upper~elevation, pred_det,type="l",col=gray(0.5))
# ###  Plot occupancy en escala original
pred_psi <-predict(mf1, type="state", newdata=newdat_range, appendData=TRUE) 
plot(Predicted ~ elevation, pred_psi, type="l", ylim=c(0,1), col="blue",
     xlab="Elevation",
     ylab="Occupancy Probability",
     xaxt="n")
xticks <- c(-1, -0.5, 0, 0.5, 1, 1.5, 2)  # -1:2
xlabs <- xticks*sd(cam_covs) + mean(cam_covs) #Use the mean and sd of original value to change label name
axis(1, at=xticks, labels=round(xlabs, 1))
lines(lower ~ elevation, pred_psi, type="l", col=gray(0.5))
lines(upper ~ elevation, pred_psi, type="l", col=gray(0.5))

Spatially Explicit Occupancy Model

Code 
library(RColorBrewer)
library(rasterVis)
Warning: package 'rasterVis' was built under R version 4.2.3
Loading required package: lattice
Code 
library(mapview)

 srtm_crop_s <- stack((srtm_crop))#, 
                      #scale(roughness)) # scale altitud
 names(srtm_crop_s) <- c("elevation")#, "roughness")
 crs(srtm_crop_s) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
pred_psi_s <-predict(mf1, type="state", newdata=srtm_crop_s) 
pred_psi_r <- pred_psi_s # * sd(full_covs$elevation) + mean(full_covs$elevation)
crs(pred_psi_r) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
clr <- colorRampPalette(brewer.pal(9, "YlGn"))
# mapview (pred_psi_r[[1]], col.regions = clr,  legend = TRUE, alpha=0.7)
# plot(pred_psi_s[[1]], main="Occupancy")
levelplot(pred_psi_r[[1]], par.settings = YlOrRdTheme(), margin=FALSE, main="Ocupancy")

Code 
pal = mapviewPalette("mapviewSpectralColors")
mapview(pred_psi_r[[1]], map.types = c("Esri.WorldImagery"), alpha=0.5, col.regions = pal(100)) + mapview(puntos_sf)