Occupancy five species. Parita Bay, Panama

spatial prediction

Read spatial layers

Code

# load function
source("C:/CodigoR/AudubonPanama/R/matrix_creator.R")
############### load Packages
library(readr) # readr 
library(readxl) # read excel
library(tidyverse) # 
library(unmarked) # occupancy by maximun likelihood
library(sf) # simple feature
library(mapview) # view map
library(ubms) # occupancy by Bayesian made easy
library(stocc) # spatial data for example
library(terra) 
library(raster)
library(RColorBrewer)
library(rasterVis)


############### load data
# old 
# parita_data_full <- read_delim("C:/CodigoR/AudubonPanama/data/all_records_audubon_panama_to_occu.csv", 
#                               delim = "\t", escape_double = FALSE, 
#                               col_types = cols(eventDate = col_date(format = "%Y-%m-%d")), 
#                               trim_ws = TRUE)

# new
parita_data_full <- read_excel("C:/CodigoR/AudubonPanama/data/all_records_audubon_panama_Darwin_Core_revision.xlsx", 
    col_types = c("text", "text", "text", 
        "date", "text", "text", "numeric", 
        "numeric", "text", "text", "numeric", 
        "text", "numeric", "numeric", "text", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "text", "text", 
        "numeric", "text", "numeric", "text", 
        "text", "text", "text", "text", "numeric", 
        "text", "numeric", "text", "numeric", 
        "text", "text"))

covs <- read_csv("C:/CodigoR/AudubonPanama/shp/sites_covs_parita_nona.csv")

## Nyctidromus albicollis
## Aramus guarauna
## Buteogallus anthracinus
## Dryocopus lineatus
## Piranga rubra

# filter by threshold and species
parita_data_5sp <-  parita_data_full %>% 
                      filter(confidence >= threshold) %>%
                      filter(scientificName == c("Nyctidromus albicollis", 
                                                 "Aramus guarauna",
                                                 "Buteogallus anthracinus",
                                                 "Dryocopus lineatus",
                                                 "Piranga rubra"
                                                 ))  

# insert Max and min to get matrix limits
# start_date
parita_data2 <- parita_data_5sp %>% 
  group_by(locationID) %>% 
  mutate(start_date = min(eventDate),
         end_date= max(eventDate))

# apply function to get matrix per species
parita <- f.matrix.creator2(data = parita_data2, year = 2023)


# centroide for terra
# centoide <- centroids(puntos, TRUE)
centroide <- c(mean(as.numeric(parita_data_5sp$decimalLongitude)), mean(as.numeric(parita_data_5sp$decimalLatitude)))
clip_window <- extent(-80.67105, -80.05822, 7.713744, 8.489888) # 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)
 #rast(srtm_crop) # convert to terra


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

# convierte covs a puntos terra
puntos <- vect(covs, geom=c("Longitude", "Latitude"), crs="EPSG:4326")
# convierte a sf
puntos_sf <- sf::st_as_sf(puntos)
# extract elev from raster
cam_covs <- raster::extract(srtm_crop, puntos_sf)


AGB_Spawn <- rast("C:/CodigoR/AudubonPanama/raster/AGB Spawn.tif")
NDVI <- rast("C:/CodigoR/AudubonPanama/raster/S2_NDVI_median_v2.tif")
roads <- rast("C:/CodigoR/AudubonPanama/raster/roads_final_v2.tif")
carbon_stock <- rast("C:/CodigoR/AudubonPanama/raster/soil_organic_carbon_stock_0-30m.tif")
canopy <- rast("C:/CodigoR/AudubonPanama/raster/canopy_height_jetz2.tif")
# make elevation equal
srtm_projected <- projectRaster(srtm_crop, crs=projection(NDVI), method="ngb") 
elev <- mask(resample(rast(srtm_projected), NDVI), NDVI)



# list of terras SpatRasters  
many_rasters <- list(AGB_Spawn, NDVI, roads, carbon_stock, elev, canopy)
# terra stack
covs_raster <- rast(many_rasters)
names(covs_raster) <- c("AGB_Spawn", "NDVI", "roads", "carbon_stock", "elevation", "canopy")

# change NA to 0 
covs_raster0 <- subst(covs_raster, NA, 0)
# extract covariates
covs_raster_val <- extract(covs_raster0, puntos_sf)

plot(covs_raster)

Code

##########
# chk names of species
# dete es hora, occur es presencia en el dia
# names(parita$dete[]) 

Buteogallus anthracinus

One species, single season

Code

# Make UMF object
umf_y_full_2<- unmarkedFrameOccu(y= as.data.frame(parita$occur[[2]])) #"Aramus guarauna"
siteCovs(umf_y_full_2) <- data.frame(elevation =  covs_raster_val$elevation,
                                   AGB_Spawn = covs_raster_val$AGB_Spawn,
                                   roads = covs_raster_val$roads,
                                   carbon = covs_raster_val$carbon_stock,
                                   NDVI = covs_raster_val$NDVI,
                                   canopy = covs_raster_val$canopy
) # data.frame(Elev=full_covs$Elev) # Full


plot(umf_y_full_2, main="Buteogallus anthracinus")  #""Buteogallus anthracinus 

Code

# obscanop<-matrix(rep(covs_raster_val$canopy_height_jetz2,27),ncol = 27) # make elevetion per observation

# detection
obsCovs(umf_y_full_2) <- list(hour = matrix(parita$dete[[2]]),
                              canopy = covs_raster_val$canopy ) # hora


# fit unmarked
fit_unm_sp2_0 <- unmarked::occu(~1~1, data=umf_y_full_2) # ok!
fit_unm_sp2_0h <- unmarked::occu(~scale(hour)~1, data=umf_y_full_2) # It work!
fit_unm_sp2_0c <- unmarked::occu(~scale(canopy)~1, data=umf_y_full_2) # It work!
fit_unm_sp2_1 <- unmarked::occu(~1~scale(elevation), data=umf_y_full_2)
fit_unm_sp2_2 <- unmarked::occu(~1~scale(AGB_Spawn), data=umf_y_full_2)
fit_unm_sp2_3 <- unmarked::occu(~1~scale(roads), data=umf_y_full_2) # ok!
fit_unm_sp2_4 <- unmarked::occu(~1~scale(carbon), data=umf_y_full_2)
fit_unm_sp2_5 <- unmarked::occu(~1~scale(NDVI), data=umf_y_full_2)
fit_unm_sp2_6 <- unmarked::occu(~1~scale(canopy), data=umf_y_full_2)


# fit list for detection
fms1<-fitList("p(.) Ocu(.)"=fit_unm_sp2_0,
              #"p(h) Ocu(.)"=fit_unm_sp2_0h, # NA problem!
              "p(c) Ocu(.)"=fit_unm_sp2_0c)#,
#              "p(.) Ocu(elev)"=fit_unm_sp2_1,
#              "p(.) Ocu(slope)"=fit_unm_sp2_2,
#              "p(.) Ocu(aspect)"=fit_unm_sp2_3,
#              "p(.) Ocu(roughness)"=fit_unm_sp2_4,
#              "p(elev^2) Ocu(elev^2)"=fit_unm_sp2_5 )

# model selection detection unmarked
modSel(fms1) #problem!.. null model ranks better

            nPars    AIC delta AICwt cumltvWt
p(.) Ocu(.)     2 422.80  0.00  0.71     0.71
p(c) Ocu(.)     3 424.60  1.80  0.29     1.00

Code

# fit list for detection
fms2<-fitList("p(.) Ocu(.)"=fit_unm_sp2_0,
              #"p(h) Ocu(.)"=fit_unm_sp2_0h, #,
              #"p(c) Ocu(.)"=fit_unm_sp2_0c,
              "p(.) Ocu(elev)"=fit_unm_sp2_1,
              "p(.) Ocu(AGB_Spawn)"=fit_unm_sp2_2,
              "p(.) Ocu(roads)"=fit_unm_sp2_3,
              "p(.) Ocu(carbon)"=fit_unm_sp2_4,
              "p(.) Ocu(NDVI)"=fit_unm_sp2_5,
              "p(.) Ocu(canopy)"=fit_unm_sp2_6)

# model selection detection unmarked
modSel(fms2)

                    nPars    AIC delta AICwt cumltvWt
p(.) Ocu(.)             2 422.80  0.00  0.31     0.31
p(.) Ocu(canopy)        3 424.80  1.99  0.11     0.43
p(.) Ocu(NDVI)          3 424.80  1.99  0.11     0.54
p(.) Ocu(AGB_Spawn)     3 424.80  1.99  0.11     0.66
p(.) Ocu(elev)          3 424.80  1.99  0.11     0.77
p(.) Ocu(carbon)        3 424.80  1.99  0.11     0.89
p(.) Ocu(roads)         3 424.80  1.99  0.11     1.00

Code

######## ubms, bayesian occupancy
# Double formula: first part is for detection, second for occupancy
#  form <- ~detCov1 + detCov2 ~habCov1 + habCov2 + RSR(x, y, threshold=1)

# fit stan
fit_stan_sp2_0 <- stan_occu(~1~1, data=umf_y_full_2, chains=3, iter=1000, cores=3)
fit_stan_sp2_0h <- stan_occu(~scale(hour)~1, data=umf_y_full_2, chains=3, iter=1000, cores=3)
fit_stan_sp2_0c <- stan_occu(~scale(canopy)~1, data=umf_y_full_2, chains=3, iter=1000, cores=3)
fit_stan_sp2_1 <- stan_occu(~scale(hour)~scale(elevation), data=umf_y_full_2, chains=3, iter=1000, cores=3)
fit_stan_sp2_2 <- stan_occu(~scale(hour)~scale(AGB_Spawn), data=umf_y_full_2, chains=3, iter=1000, cores=3)
fit_stan_sp2_3 <- stan_occu(~scale(hour)~scale(roads), data=umf_y_full_2, chains=3, iter=1000, cores=3)
fit_stan_sp2_4 <- stan_occu(~scale(hour)~scale(carbon), data=umf_y_full_2, chains=3, iter=1000, cores=3)
fit_stan_sp2_5 <- stan_occu(~scale(hour)~scale(NDVI), data=umf_y_full_2, chains=3, iter=1000, cores=3)
fit_stan_sp2_6 <- stan_occu(~scale(hour)~scale(canopy), data=umf_y_full_2, chains=3, iter=1000, cores=3)


mods2 <- fitList(fit_stan_sp2_0,
                fit_stan_sp2_0h,
                fit_stan_sp2_0c,
                fit_stan_sp2_1,
                fit_stan_sp2_2,
                fit_stan_sp2_3,
                fit_stan_sp2_4,
                fit_stan_sp2_5,
                fit_stan_sp2_6)

# model selection
round(modSel(mods2), 3)

                    elpd nparam elpd_diff se_diff weight
fit_stan_sp2_0h -166.235  4.188     0.000   0.000  0.323
fit_stan_sp2_1  -166.236  4.159     0.000   0.122  0.306
fit_stan_sp2_6  -166.290  4.227    -0.054   0.165  0.222
fit_stan_sp2_4  -166.366  4.305    -0.130   0.092  0.137
fit_stan_sp2_2  -166.464  4.469    -0.229   0.122  0.012
fit_stan_sp2_3  -166.672  4.656    -0.436   0.159  0.000
fit_stan_sp2_5  -166.843  4.807    -0.608   0.257  0.000
fit_stan_sp2_0  -212.648  3.487   -46.412  10.437  0.000
fit_stan_sp2_0c -213.615  4.264   -47.379  10.705  0.000

Code

# plot top model
print("Hour")

[1] "Hour"

Code

plot_effects(fit_stan_sp2_5, "det") # Detection

Code

print("NDVI")

[1] "NDVI"

Code

plot_effects(fit_stan_sp2_5, "state") # Occupancy

Spatial Prediction

Code

#### predict


# srtm_crop_s <- mask(elev, NDVI)#stack(raster(elev))#, 
#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"
nd2 <- NDVI # mask(elev, NDVI) #data.frame(srtm_crop_s, hour=0.5)
names(nd2) <- "NDVI"
# newd <- data.frame(srtm_crop_s, hour=0.5)
pred_psi_s_unm <-predict(fit_unm_sp2_5, type="state", newdata=nd2)
#pred_psi_s_stan <-predict(fit_stan_sp2_1, submodel="state", newdata=nd, na.rm=FALSE)
#pred_psi_ras <- raster(srtm_crop_s) # make raster = to elev 
#pred_psi_ras[] <- pred_psi_s_stan$Predicted # drape on top


# pred_psi_r <- pred_psi_s # * sd(full_covs$elevation) + mean(full_covs$elevation)
# crs(pred_psi_ras) <- "+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_s_unm[[1]], par.settings = YlOrRdTheme(), margin=FALSE, main="Ocupancy")

Code

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

Dryocopus lineatus

One species, single season

Code

# Make UMF object
umf_y_full_3<- unmarkedFrameOccu(y= as.data.frame(parita$occur[[3]])) #Dryocopus lineatus
siteCovs(umf_y_full_3) <- data.frame(elevation =  covs_raster_val$elevation,
                                   AGB_Spawn = covs_raster_val$AGB_Spawn,
                                   roads = covs_raster_val$roads,
                                   carbon = covs_raster_val$carbon_stock,
                                   NDVI = covs_raster_val$NDVI,
                                   canopy = covs_raster_val$canopy
) # data.frame(Elev=full_covs$Elev) # Full


plot(umf_y_full_3, main="Dryocopus lineatus")  #""Buteogallus anthracinus no converge

Code

# obscanop<-matrix(rep(covs_raster_val$canopy_height_jetz2,27),ncol = 27) # make elevetion per observation

# detection
obsCovs(umf_y_full_3) <- list(hour = matrix(parita$dete[[3]]),
                              canopy = covs_raster_val$canopy ) # hora


# fit unmarked
fit_unm_sp3_0 <- unmarked::occu(~1~1, data=umf_y_full_3) # ok!
fit_unm_sp3_0h <- unmarked::occu(~scale(hour)~1, data=umf_y_full_3) # It work!
fit_unm_sp3_0c <- unmarked::occu(~scale(canopy)~1, data=umf_y_full_3) # It work!
fit_unm_sp3_1 <- unmarked::occu(~1~scale(elevation), data=umf_y_full_3)
fit_unm_sp3_2 <- unmarked::occu(~1~scale(AGB_Spawn), data=umf_y_full_3)
fit_unm_sp3_3 <- unmarked::occu(~1~scale(roads), data=umf_y_full_3) # ok!
fit_unm_sp3_4 <- unmarked::occu(~1~scale(carbon), data=umf_y_full_3)
fit_unm_sp3_5 <- unmarked::occu(~1~scale(NDVI), data=umf_y_full_3)
fit_unm_sp3_6 <- unmarked::occu(~1~scale(canopy), data=umf_y_full_3)


# fit list for detection
fms1<-fitList("p(.) Ocu(.)"=fit_unm_sp3_0,
              # "p(h) Ocu(.)"=fit_unm_sp3_0h, # NA problem!
              "p(c) Ocu(.)"=fit_unm_sp3_0c)#,
#              "p(.) Ocu(elev)"=fit_unm_sp2_1,
#              "p(.) Ocu(slope)"=fit_unm_sp2_2,
#              "p(.) Ocu(aspect)"=fit_unm_sp2_3,
#              "p(.) Ocu(roughness)"=fit_unm_sp2_4,
#              "p(elev^2) Ocu(elev^2)"=fit_unm_sp2_5 )

# model selection detection unmarked
modSel(fms1) #problem!.. null model ranks better

            nPars   AIC delta AICwt cumltvWt
p(.) Ocu(.)     2 55.03  0.00  0.73     0.73
p(c) Ocu(.)     3 56.97  1.94  0.27     1.00

Code

# fit list for detection
fms2<-fitList("p(.) Ocu(.)"=fit_unm_sp3_0,
              #"p(h) Ocu(.)"=fit_unm_sp2_0h, #,
              #"p(c) Ocu(.)"=fit_unm_sp3_0c,
              "p(.) Ocu(elev)"=fit_unm_sp3_1,
              "p(.) Ocu(AGB_Spawn)"=fit_unm_sp3_2,
              "p(.) Ocu(roads)"=fit_unm_sp3_3,
              "p(.) Ocu(carbon)"=fit_unm_sp3_4,
              "p(.) Ocu(NDVI)"=fit_unm_sp3_5,
              "p(.) Ocu(canopy)"=fit_unm_sp3_6)

# model selection detection unmarked
modSel(fms2)

                    nPars   AIC delta AICwt cumltvWt
p(.) Ocu(AGB_Spawn)     3 52.42  0.00 0.441     0.44
p(.) Ocu(.)             2 55.03  2.61 0.120     0.56
p(.) Ocu(canopy)        3 55.09  2.67 0.116     0.68
p(.) Ocu(NDVI)          3 55.21  2.79 0.109     0.78
p(.) Ocu(elev)          3 55.52  3.10 0.094     0.88
p(.) Ocu(carbon)        3 55.89  3.47 0.078     0.96
p(.) Ocu(roads)         3 57.03  4.61 0.044     1.00

Code

# fit stan
fit_stan_sp3_0 <- stan_occu(~1~1, data=umf_y_full_3, chains=3, iter=1000, cores=3)
fit_stan_sp3_0h <- stan_occu(~scale(hour)~1, data=umf_y_full_3, chains=3, iter=1000, cores=3)
fit_stan_sp3_1 <- stan_occu(~scale(hour)~scale(elevation), data=umf_y_full_3, chains=3, iter=1000, cores=3)
fit_stan_sp3_2 <- stan_occu(~scale(hour)~scale(AGB_Spawn), data=umf_y_full_3, chains=3, iter=1000, cores=3)
fit_stan_sp3_3 <- stan_occu(~scale(hour)~scale(roads), data=umf_y_full_3, chains=3, iter=1000, cores=3)
fit_stan_sp3_4 <- stan_occu(~scale(hour)~scale(carbon), data=umf_y_full_3, chains=3, iter=1000, cores=3)
fit_stan_sp3_5 <- stan_occu(~scale(hour)~scale(NDVI), data=umf_y_full_3, chains=3, iter=1000, cores=3)
fit_stan_sp3_6 <- stan_occu(~scale(hour)~scale(canopy), data=umf_y_full_3, chains=3, iter=1000, cores=3)


mods3 <- fitList(fit_stan_sp3_0,
                fit_stan_sp3_0h,
                fit_stan_sp3_1,
                fit_stan_sp3_2,
                fit_stan_sp3_3,
                fit_stan_sp3_4,
                fit_stan_sp3_5,
                fit_stan_sp3_6)

# model selection
round(modSel(mods3), 3) #AGB_Spawn

                   elpd nparam elpd_diff se_diff weight
fit_stan_sp3_0h -22.435  1.776     0.000   0.000  0.722
fit_stan_sp3_5  -22.562  2.123    -0.127   0.660  0.221
fit_stan_sp3_6  -22.800  2.240    -0.365   0.469  0.000
fit_stan_sp3_2  -22.963  2.440    -0.529   0.364  0.000
fit_stan_sp3_4  -22.985  2.768    -0.550   1.168  0.057
fit_stan_sp3_1  -23.093  2.665    -0.658   0.753  0.000
fit_stan_sp3_3  -23.384  2.932    -0.949   0.754  0.000
fit_stan_sp3_0  -27.272  1.237    -4.837   4.850  0.000

Code

# plot top model
print("Hour")

[1] "Hour"

Code

plot_effects(fit_stan_sp3_5, "det") # Detectiom

Code

print("NDVI")

[1] "NDVI"

Code

plot_effects(fit_stan_sp3_5, "state") # Ocupancy

Spatial Prediction

Code

#### predict

# srtm_crop_s <- mask(elev, NDVI)#stack(raster(elev))#, 
# 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"
nd3 <- NDVI #mask(AGB_Spawn, NDVI) #data.frame(srtm_crop_s, hour=0.5)
names(nd3) <- "NDVI"
# newd <- data.frame(srtm_crop_s, hour=0.5)
pred_psi_s_unm <-predict(fit_unm_sp3_5, type="state", newdata=nd3)
#pred_psi_s_stan <-predict(fit_stan_sp2_1, submodel="state", newdata=nd, na.rm=FALSE)
#pred_psi_ras <- raster(srtm_crop_s) # make raster = to elev 
#pred_psi_ras[] <- pred_psi_s_stan$Predicted # drape on top


# pred_psi_r <- pred_psi_s # * sd(full_covs$elevation) + mean(full_covs$elevation)
# crs(pred_psi_ras) <- "+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_s_unm[[1]], par.settings = YlOrRdTheme(), margin=FALSE, main="Ocupancy")

Code

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

Nyctidromus albicollis

One species, single season

Code

# Make UMF object
umf_y_full_4<- unmarkedFrameOccu(y= as.data.frame(parita$occur[[4]])) #Dryocopus lineatus
siteCovs(umf_y_full_4) <- data.frame(elevation =  covs_raster_val$elevation,
                                   AGB_Spawn = covs_raster_val$AGB_Spawn,
                                   roads = covs_raster_val$roads,
                                   carbon = covs_raster_val$carbon_stock,
                                   NDVI = covs_raster_val$NDVI,
                                   canopy = covs_raster_val$canopy
) # data.frame(Elev=full_covs$Elev) # Full


plot(umf_y_full_4, main="Nyctidromus albicollis")  #""Buteogallus anthracinus no converge

Code

# obscanop<-matrix(rep(covs_raster_val$canopy_height_jetz2,27),ncol = 27) # make elevetion per observation

# detection
obsCovs(umf_y_full_4) <- list(hour = matrix(parita$dete[[4]]),
                              canopy = covs_raster_val$canopy ) # hora


# fit unmarked
fit_unm_sp4_0 <- unmarked::occu(~1~1, data=umf_y_full_4) # ok!
fit_unm_sp4_0h <- unmarked::occu(~scale(hour)~1, data=umf_y_full_4) # It work!
fit_unm_sp4_0c <- unmarked::occu(~scale(canopy)~1, data=umf_y_full_4) # It work!
fit_unm_sp4_1 <- unmarked::occu(~scale(canopy)~scale(elevation), data=umf_y_full_4)
fit_unm_sp4_2 <- unmarked::occu(~scale(canopy)~scale(AGB_Spawn), data=umf_y_full_4)
fit_unm_sp4_3 <- unmarked::occu(~scale(canopy)~scale(roads), data=umf_y_full_4) # ok!
fit_unm_sp4_4 <- unmarked::occu(~scale(canopy)~scale(carbon), data=umf_y_full_4)
fit_unm_sp4_5 <- unmarked::occu(~scale(canopy)~scale(NDVI), data=umf_y_full_4)
fit_unm_sp4_6 <- unmarked::occu(~scale(canopy)~scale(canopy), data=umf_y_full_4)
fit_unm_sp4_2b <- unmarked::occu(~1~scale(AGB_Spawn), data=umf_y_full_4)



# fit list for detection
fms1<-fitList("p(.) Ocu(.)"=fit_unm_sp4_0,
              # "p(h) Ocu(.)"=fit_unm_sp3_0h, # NA problem!
              "p(c) Ocu(.)"=fit_unm_sp4_0c)#,
#              "p(.) Ocu(elev)"=fit_unm_sp2_1,
#              "p(.) Ocu(slope)"=fit_unm_sp2_2,
#              "p(.) Ocu(aspect)"=fit_unm_sp2_3,
#              "p(.) Ocu(roughness)"=fit_unm_sp2_4,
#              "p(elev^2) Ocu(elev^2)"=fit_unm_sp2_5 )

# model selection detection unmarked
modSel(fms1) #good!.. canopy  model ranks better

            nPars    AIC delta AICwt cumltvWt
p(c) Ocu(.)     3 165.54  0.00  0.53     0.53
p(.) Ocu(.)     2 165.80  0.26  0.47     1.00

Code

# fit list for detection
fms2<-fitList("p(.) Ocu(.)"=fit_unm_sp4_0,
              #"p(h) Ocu(.)"=fit_unm_sp2_0h, #,
              "p(c) Ocu(.)"=fit_unm_sp4_0c,
              "p(c) Ocu(elev)"=fit_unm_sp4_1,
              "p(c) Ocu(AGB_Spawn)"=fit_unm_sp4_2,
              "p(c) Ocu(roads)"=fit_unm_sp4_3,
              "p(c) Ocu(carbon)"=fit_unm_sp4_4,
              "p(c) Ocu(NDVI)"=fit_unm_sp4_5,
              "p(c) Ocu(canopy)"=fit_unm_sp4_6)

# model selection detection unmarked
modSel(fms2)

                    nPars    AIC delta AICwt cumltvWt
p(c) Ocu(AGB_Spawn)     4 164.85  0.00 0.226     0.23
p(c) Ocu(canopy)        4 165.42  0.57 0.170     0.40
p(c) Ocu(.)             3 165.54  0.69 0.160     0.56
p(.) Ocu(.)             2 165.80  0.95 0.141     0.70
p(c) Ocu(carbon)        4 166.25  1.40 0.112     0.81
p(c) Ocu(elev)          4 167.14  2.29 0.072     0.88
p(c) Ocu(roads)         4 167.49  2.64 0.060     0.94
p(c) Ocu(NDVI)          4 167.54  2.69 0.059     1.00

Code

# fit stan
fit_stan_sp4_0 <- stan_occu(~1~1, data=umf_y_full_4, chains=3, iter=1000, cores=3)
fit_stan_sp4_0h <- stan_occu(~scale(hour)~1, data=umf_y_full_4, chains=3, iter=1000, cores=3)
fit_stan_sp4_0c <- stan_occu(~scale(canopy)~1, data=umf_y_full_4, chains=3, iter=1000, cores=3)
fit_stan_sp4_1 <- stan_occu(~1~scale(elevation), data=umf_y_full_4, chains=3, iter=1000, cores=3)
fit_stan_sp4_2 <- stan_occu(~scale(hour)~scale(AGB_Spawn), data=umf_y_full_4, chains=3, iter=1000, cores=3)
fit_stan_sp4_3 <- stan_occu(~scale(hour)~scale(roads), data=umf_y_full_4, chains=3, iter=1000, cores=3)
fit_stan_sp4_4 <- stan_occu(~scale(hour)~scale(carbon), data=umf_y_full_4, chains=3, iter=1000, cores=3)
fit_stan_sp4_5 <- stan_occu(~scale(hour)~scale(NDVI), data=umf_y_full_4, chains=3, iter=1000, cores=3)


mods4 <- fitList(fit_stan_sp4_0,
                fit_stan_sp4_0h,
                fit_stan_sp4_0c,
                fit_stan_sp4_1,
                fit_stan_sp4_2,
                fit_stan_sp4_3,
                fit_stan_sp4_4,
                fit_stan_sp4_5)

# model selection stan
round(modSel(mods4), 3)

                   elpd nparam elpd_diff se_diff weight
fit_stan_sp4_2  -64.898  9.173     0.000   0.000  0.678
fit_stan_sp4_0h -65.492  8.621    -0.595   1.685  0.321
fit_stan_sp4_5  -65.824  9.388    -0.926   1.754  0.000
fit_stan_sp4_4  -66.306  9.799    -1.409   2.110  0.000
fit_stan_sp4_3  -66.364  9.529    -1.467   1.881  0.000
fit_stan_sp4_0c -86.170  8.784   -21.272  10.101  0.000
fit_stan_sp4_0  -86.632  8.468   -21.735  10.231  0.000
fit_stan_sp4_1  -88.611 10.676   -23.713  10.551  0.000

Code

# plot top model
print("Hour")

[1] "Hour"

Code

plot_effects(fit_stan_sp4_2, "det") # Detection

Code

print("AGB_Spawn")

[1] "AGB_Spawn"

Code

plot_effects(fit_stan_sp4_2, "state") # Occupancy

Spatial Prediction

Code

#### predict

# srtm_crop_s <- mask(elev, NDVI)#stack(raster(elev))#, 
# 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"
nd4 <- mask(AGB_Spawn, NDVI) #data.frame(srtm_crop_s, hour=0.5)
names(nd4) <- "AGB_Spawn"
# newd <- data.frame(srtm_crop_s, hour=0.5)
pred_psi_s_unm <-predict(fit_unm_sp4_2b, type="state", newdata=nd4)
#pred_psi_s_stan <-predict(fit_stan_sp2_1, submodel="state", newdata=nd, na.rm=FALSE)
#pred_psi_ras <- raster(srtm_crop_s) # make raster = to elev 
#pred_psi_ras[] <- pred_psi_s_stan$Predicted # drape on top


# pred_psi_r <- pred_psi_s # * sd(full_covs$elevation) + mean(full_covs$elevation)
# crs(pred_psi_ras) <- "+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_s_unm[[1]], par.settings = YlOrRdTheme(), margin=FALSE, main="Ocupancy")

Code

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

Piranga rubra

One species, single season

Code

# Make UMF object
umf_y_full_5<- unmarkedFrameOccu(y= as.data.frame(parita$occur[[5]])) #Dryocopus lineatus
siteCovs(umf_y_full_5) <- data.frame(elevation =  covs_raster_val$elevation,
                                   AGB_Spawn = covs_raster_val$AGB_Spawn,
                                   roads = covs_raster_val$roads,
                                   carbon = covs_raster_val$carbon_stock,
                                   NDVI = covs_raster_val$NDVI,
                                   canopy = covs_raster_val$canopy
) # data.frame(Elev=full_covs$Elev) # Full


plot(umf_y_full_5, main="Piranga rubra")  #""Buteogallus anthracinus no converge

Code

# detection
obsCovs(umf_y_full_5) <- list(hour = matrix(parita$dete[[5]]),
                              canopy = covs_raster_val$canopy ) # hora

# Double formula: first part is for detection, second for occupancy
#  form <- ~detCov1 + detCov2 ~habCov1 + habCov2 + RSR(x, y, threshold=1)

# fit unmarked
fit_unm_sp5_0 <- unmarked::occu(~1~1, data=umf_y_full_5) # ok!
# fit_unm_sp4_0h <- unmarked::occu(~hour~1, data=umf_y_full_4) # Does not work!
fit_unm_sp5_0c <- unmarked::occu(~scale(canopy)~1, data=umf_y_full_5) # Does not work!
fit_unm_sp5_1 <- unmarked::occu(~scale(canopy)~scale(elevation), data=umf_y_full_5)
fit_unm_sp5_2 <- unmarked::occu(~scale(canopy)~scale(AGB_Spawn), data=umf_y_full_5)
fit_unm_sp5_3 <- unmarked::occu(~scale(canopy)~scale(roads), data=umf_y_full_5) # ok!
fit_unm_sp5_4 <- unmarked::occu(~scale(canopy)~scale(carbon), data=umf_y_full_5)
fit_unm_sp5_5 <- unmarked::occu(~scale(canopy)~scale(NDVI), data=umf_y_full_5)
fit_unm_sp5_2b <- unmarked::occu(~1~scale(AGB_Spawn), data=umf_y_full_5)


# fit list for detection
fms1<-fitList("p(.) Ocu(.)"=fit_unm_sp5_0,
              # "p(h) Ocu(.)"=fit_unm_sp3_0h, # NA problem!
              "p(c) Ocu(.)"=fit_unm_sp5_0c)#,
#              "p(.) Ocu(elev)"=fit_unm_sp2_1,
#              "p(.) Ocu(slope)"=fit_unm_sp2_2,
#              "p(.) Ocu(aspect)"=fit_unm_sp2_3,
#              "p(.) Ocu(roughness)"=fit_unm_sp2_4,
#              "p(elev^2) Ocu(elev^2)"=fit_unm_sp2_5 )

# model selection detection unmarked
modSel(fms1) #good!.. canopy  model ranks better

            nPars    AIC delta AICwt cumltvWt
p(.) Ocu(.)     2 240.99  0.00  0.67     0.67
p(c) Ocu(.)     3 242.41  1.42  0.33     1.00

Code

# fit list for detection
fms2<-fitList("p(.) Ocu(.)"=fit_unm_sp4_0,
              #"p(h) Ocu(.)"=fit_unm_sp2_0h, #,
              "p(c) Ocu(.)"=fit_unm_sp4_0c,
              "p(c) Ocu(elev)"=fit_unm_sp4_1,
              "p(c) Ocu(AGB_Spawn)"=fit_unm_sp4_2,
              "p(c) Ocu(roads)"=fit_unm_sp4_3,
              "p(c) Ocu(carbon)"=fit_unm_sp4_4,
              "p(c) Ocu(NDVI)"=fit_unm_sp4_5,
              "p(c) Ocu(canopy)"=fit_unm_sp4_6)

# model selection detection unmarked
modSel(fms2)

                    nPars    AIC delta AICwt cumltvWt
p(c) Ocu(AGB_Spawn)     4 164.85  0.00 0.226     0.23
p(c) Ocu(canopy)        4 165.42  0.57 0.170     0.40
p(c) Ocu(.)             3 165.54  0.69 0.160     0.56
p(.) Ocu(.)             2 165.80  0.95 0.141     0.70
p(c) Ocu(carbon)        4 166.25  1.40 0.112     0.81
p(c) Ocu(elev)          4 167.14  2.29 0.072     0.88
p(c) Ocu(roads)         4 167.49  2.64 0.060     0.94
p(c) Ocu(NDVI)          4 167.54  2.69 0.059     1.00

Code

# fit stan
fit_stan_sp5_0 <- stan_occu(~1~1, data=umf_y_full_5, chains=3, iter=1000, cores=3)
fit_stan_sp5_0h <- stan_occu(~scale(hour)~1, data=umf_y_full_5, chains=3, iter=1000, cores=3)
fit_stan_sp5_0c <- stan_occu(~scale(canopy)~1, data=umf_y_full_5, chains=3, iter=1000, cores=3)
fit_stan_sp5_1 <- stan_occu(~scale(hour)~scale(elevation), data=umf_y_full_5, chains=3, iter=1000, cores=3)
fit_stan_sp5_2 <- stan_occu(~scale(hour)~scale(AGB_Spawn), data=umf_y_full_5, chains=3, iter=1000, cores=3)
fit_stan_sp5_3 <- stan_occu(~scale(hour)~scale(roads), data=umf_y_full_5, chains=3, iter=1000, cores=3)
fit_stan_sp5_4 <- stan_occu(~scale(hour)~scale(carbon), data=umf_y_full_5, chains=3, iter=1000, cores=3)
fit_stan_sp5_5 <- stan_occu(~scale(hour)~scale(NDVI), data=umf_y_full_5, chains=3, iter=1000, cores=3)
fit_stan_sp5_6 <- stan_occu(~scale(hour)~scale(canopy), data=umf_y_full_5, chains=3, iter=1000, cores=3)


mods5 <- fitList(fit_stan_sp5_0,
                fit_stan_sp5_0h,
                fit_stan_sp5_0c,
                fit_stan_sp5_1,
                fit_stan_sp5_2,
                fit_stan_sp5_3,
                fit_stan_sp5_4,
                fit_stan_sp5_5)

# model selection
round(modSel(mods5), 3)

                    elpd nparam elpd_diff se_diff weight
fit_stan_sp5_4   -75.001 12.425     0.000   0.000  0.844
fit_stan_sp5_0h  -75.533 12.399    -0.531   1.512  0.001
fit_stan_sp5_2   -75.837 12.875    -0.835   2.134  0.000
fit_stan_sp5_1   -75.847 13.529    -0.845   1.347  0.000
fit_stan_sp5_3   -75.925 12.502    -0.923   1.497  0.000
fit_stan_sp5_5   -76.710 13.744    -1.708   1.925  0.000
fit_stan_sp5_0  -128.069 15.628   -53.067  23.952  0.155
fit_stan_sp5_0c -128.271 15.503   -53.270  23.799  0.000

Code

# plot top model
print("Hour")

[1] "Hour"

Code

plot_effects(fit_stan_sp5_2, "det") # Detection

Code

print("AGB_Spawn")

[1] "AGB_Spawn"

Code

plot_effects(fit_stan_sp5_2, "state") # Occupancy

Spatial Prediction

Code

#### predict

# srtm_crop_s <- mask(elev, NDVI)#stack(raster(elev))#, 
# 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"
nd5 <- mask(AGB_Spawn, NDVI) #data.frame(srtm_crop_s, hour=0.5)
names(nd5) <- "AGB_Spawn"
# newd <- data.frame(srtm_crop_s, hour=0.5)
pred_psi_s_unm <-predict(fit_unm_sp5_2b, type="state", newdata=nd5)
#pred_psi_s_stan <-predict(fit_stan_sp2_1, submodel="state", newdata=nd, na.rm=FALSE)
#pred_psi_ras <- raster(srtm_crop_s) # make raster = to elev 
#pred_psi_ras[] <- pred_psi_s_stan$Predicted # drape on top


# pred_psi_r <- pred_psi_s # * sd(full_covs$elevation) + mean(full_covs$elevation)
# crs(pred_psi_ras) <- "+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_s_unm[[1]], par.settings = YlOrRdTheme(), margin=FALSE, main="Ocupancy")

Code

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

Aramus guarauna

One species, single season

Code

# Make UMF object
umf_y_full_1<- unmarkedFrameOccu(y= as.data.frame(parita$occur[[1]])) #Aramus guarauna
siteCovs(umf_y_full_1) <- data.frame(elevation =  covs_raster_val$elevation,
                                   AGB_Spawn = covs_raster_val$AGB_Spawn,
                                   roads = covs_raster_val$roads,
                                   carbon = covs_raster_val$carbon_stock,
                                   NDVI = covs_raster_val$NDVI,
                                   canopy = covs_raster_val$canopy
) # data.frame(Elev=full_covs$Elev) # Full


plot(umf_y_full_1, main="Aramus guarauna")  #""Buteogallus anthracinus no converge

Code

# detection
obsCovs(umf_y_full_1) <- list(hour = matrix(parita$dete[[1]]),
                              canopy = covs_raster_val$canopy ) # hora

# Double formula: first part is for detection, second for occupancy
#  form <- ~detCov1 + detCov2 ~habCov1 + habCov2 + RSR(x, y, threshold=1)

# fit unmarked
fit_unm_sp1_0 <- unmarked::occu(~1~1, data=umf_y_full_1) # Model did not converge!
fit_unm_sp1_0h <- unmarked::occu(~hour~1, data=umf_y_full_1) # Does work but not comparable!
fit_unm_sp1_0c <- unmarked::occu(~scale(canopy)~1, data=umf_y_full_1) # Does work!
fit_unm_sp1_1 <- unmarked::occu(~1~scale(elevation), data=umf_y_full_1)
fit_unm_sp1_2 <- unmarked::occu(~1~scale(AGB_Spawn), data=umf_y_full_1)
fit_unm_sp1_3 <- unmarked::occu(~1~scale(roads), data=umf_y_full_1) # ok!
fit_unm_sp1_4 <- unmarked::occu(~1~scale(carbon), data=umf_y_full_1)
fit_unm_sp1_5 <- unmarked::occu(~1~scale(NDVI), data=umf_y_full_1)
fit_unm_sp1_6 <- unmarked::occu(~1~scale(canopy), data=umf_y_full_1)
fit_unm_sp1_2b <- unmarked::occu(~1~scale(AGB_Spawn), data=umf_y_full_1)


# fit list for detection
fms1<-fitList("p(.) Ocu(.)"=fit_unm_sp1_0,
              #"p(h) Ocu(.)"=fit_unm_sp1_0h, # NA problem!
              "p(c) Ocu(.)"=fit_unm_sp1_0c)#,
#              "p(.) Ocu(elev)"=fit_unm_sp2_1,
#              "p(.) Ocu(slope)"=fit_unm_sp2_2,
#              "p(.) Ocu(aspect)"=fit_unm_sp2_3,
#              "p(.) Ocu(roughness)"=fit_unm_sp2_4,
#              "p(elev^2) Ocu(elev^2)"=fit_unm_sp2_5 )

# model selection detection unmarked
modSel(fms1) #bad!.. null  model ranks better

            nPars   AIC delta AICwt cumltvWt
p(.) Ocu(.)     2 37.70  0.00  0.61     0.61
p(c) Ocu(.)     3 38.57  0.87  0.39     1.00

Code

# fit list for detection
fms2<-fitList("p(.) Ocu(.)"=fit_unm_sp1_0,
              #"p(h) Ocu(.)"=fit_unm_sp2_0h, #,
              "p(c) Ocu(.)"=fit_unm_sp1_0c,
              "p(.) Ocu(elev)"=fit_unm_sp1_1,
              "p(.) Ocu(AGB_Spawn)"=fit_unm_sp1_2,
              "p(.) Ocu(roads)"=fit_unm_sp1_3,
              "p(.) Ocu(carbon)"=fit_unm_sp1_4,
              "p(.) Ocu(NDVI)"=fit_unm_sp1_5,
              "p(.) Ocu(canopy)"=fit_unm_sp1_6)

# model selection detection unmarked
modSel(fms2)

                    nPars   AIC delta AICwt cumltvWt
p(.) Ocu(roads)         3 36.91  0.00 0.278     0.28
p(.) Ocu(.)             2 37.70  0.79 0.187     0.46
p(c) Ocu(.)             3 38.57  1.66 0.121     0.59
p(.) Ocu(AGB_Spawn)     3 38.81  1.90 0.107     0.69
p(.) Ocu(NDVI)          3 39.11  2.20 0.093     0.79
p(.) Ocu(carbon)        3 39.49  2.58 0.076     0.86
p(.) Ocu(elev)          3 39.70  2.79 0.069     0.93
p(.) Ocu(canopy)        3 39.70  2.79 0.069     1.00

Code

# fit stan
fit_stan_sp1_0 <- stan_occu(~1~1, data=umf_y_full_1, chains=3, iter=1000, cores=3)
fit_stan_sp1_0h <- stan_occu(~scale(hour)~1, data=umf_y_full_1, chains=3, iter=1000, cores=3)
fit_stan_sp1_0c <- stan_occu(~scale(canopy)~1, data=umf_y_full_1, chains=3, iter=1000, cores=3)
fit_stan_sp1_1 <- stan_occu(~scale(hour)~scale(elevation), data=umf_y_full_1, chains=3, iter=1000, cores=3)
fit_stan_sp1_2 <- stan_occu(~scale(hour)~scale(AGB_Spawn), data=umf_y_full_1, chains=3, iter=1000, cores=3)
fit_stan_sp1_3 <- stan_occu(~scale(hour)~scale(roads), data=umf_y_full_1, chains=3, iter=1000, cores=3)
fit_stan_sp1_4 <- stan_occu(~scale(hour)~scale(carbon), data=umf_y_full_1, chains=3, iter=1000, cores=3)
fit_stan_sp1_5 <- stan_occu(~scale(hour)~scale(NDVI), data=umf_y_full_1, chains=3, iter=1000, cores=3)
fit_stan_sp1_6 <- stan_occu(~scale(hour)~scale(canopy), data=umf_y_full_1, chains=3, iter=1000, cores=3)


mods5 <- fitList(fit_stan_sp1_0,
                fit_stan_sp1_0h,
                fit_stan_sp1_0c,
                fit_stan_sp1_1,
                fit_stan_sp1_2,
                fit_stan_sp1_3,
                fit_stan_sp1_4,
                fit_stan_sp1_5)

# model selection
round(modSel(mods5), 3)

                   elpd nparam elpd_diff se_diff weight
fit_stan_sp1_5  -13.582  2.247     0.000   0.000  0.648
fit_stan_sp1_0h -13.583  1.449    -0.001   0.801  0.243
fit_stan_sp1_3  -13.820  1.930    -0.238   0.999  0.000
fit_stan_sp1_2  -13.923  2.004    -0.340   1.091  0.000
fit_stan_sp1_4  -14.460  2.444    -0.878   0.942  0.000
fit_stan_sp1_1  -14.597  2.619    -1.015   0.717  0.000
fit_stan_sp1_0  -18.495  1.190    -4.913   5.316  0.000
fit_stan_sp1_0c -18.671  1.682    -5.088   5.894  0.108

Code

# plot top model
print("Hour")

[1] "Hour"

Code

plot_effects(fit_stan_sp1_5, "det") # Detection

Code

print("NDVI")

[1] "NDVI"

Code

plot_effects(fit_stan_sp1_5, "state") # Occupancy

Spatial Prediction

Code

#### predict

# srtm_crop_s <- mask(elev, NDVI)#stack(raster(elev))#, 
# 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"
nd1 <- NDVI #data.frame(srtm_crop_s, hour=0.5)
names(nd1) <- "NDVI"
# newd <- data.frame(srtm_crop_s, hour=0.5)
pred_psi_s_unm <-predict(fit_unm_sp1_5, type="state", newdata=nd1)
#pred_psi_s_stan <-predict(fit_stan_sp2_1, submodel="state", newdata=nd, na.rm=FALSE)
#pred_psi_ras <- raster(srtm_crop_s) # make raster = to elev 
#pred_psi_ras[] <- pred_psi_s_stan$Predicted # drape on top


# pred_psi_r <- pred_psi_s # * sd(full_covs$elevation) + mean(full_covs$elevation)
# crs(pred_psi_ras) <- "+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_s_unm[[1]], par.settings = YlOrRdTheme(), margin=FALSE, main="Ocupancy")

Code

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

Información de la sesión en R.

Code

print(sessionInfo(), locale = FALSE)

R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] rasterVis_0.51.5   lattice_0.20-45    RColorBrewer_1.1-3 raster_3.6-20     
 [5] sp_1.6-0           terra_1.7-3        stocc_1.31         ubms_1.2.6        
 [9] mapview_2.11.0     sf_1.0-9           unmarked_1.4.1     lubridate_1.9.2   
[13] forcats_1.0.0      stringr_1.5.0      dplyr_1.1.0        purrr_1.0.1       
[17] tidyr_1.3.0        tibble_3.2.0       ggplot2_3.4.1      tidyverse_2.0.0   
[21] readxl_1.4.2       readr_2.1.4       

loaded via a namespace (and not attached):
  [1] minqa_1.2.5             leafem_0.2.0            colorspace_2.1-0       
  [4] deldir_1.0-6            ellipsis_0.3.2          class_7.3-20           
  [7] rgdal_1.6-4             leaflet_2.1.1           satellite_1.0.4        
 [10] QuickJSR_1.1.3          base64enc_0.1-3         rstudioapi_0.14        
 [13] proxy_0.4-27            farver_2.1.1            hexbin_1.28.3          
 [16] rstan_2.32.5            bit64_4.0.5             RSpectra_0.16-1        
 [19] fansi_1.0.4             splines_4.2.2           codetools_0.2-18       
 [22] knitr_1.42              spam_2.10-0             jsonlite_1.8.4         
 [25] nloptr_2.0.3            png_0.1-8               compiler_4.2.2         
 [28] backports_1.4.1         Matrix_1.5-3            fastmap_1.1.1          
 [31] cli_3.6.0               leaflet.providers_1.9.0 htmltools_0.5.4        
 [34] prettyunits_1.1.1       tools_4.2.2             dotCall64_1.1-1        
 [37] coda_0.19-4             gtable_0.3.2            glue_1.6.2             
 [40] posterior_1.4.1         maps_3.4.2              Rcpp_1.0.10            
 [43] cellranger_1.1.0        vctrs_0.5.2             svglite_2.1.1          
 [46] nlme_3.1-160            crosstalk_1.2.0         tensorA_0.36.2         
 [49] xfun_0.37               ps_1.7.2                lme4_1.1-31            
 [52] timechange_0.2.0        lifecycle_1.0.3         MASS_7.3-58.1          
 [55] zoo_1.8-12              scales_1.2.1            vroom_1.6.1            
 [58] hms_1.1.2               parallel_4.2.2          inline_0.3.19          
 [61] fields_15.2             yaml_2.3.7              pbapply_1.7-0          
 [64] gridExtra_2.3           loo_2.6.0               StanHeaders_2.32.5     
 [67] latticeExtra_0.6-30     stringi_1.7.12          leafpop_0.1.0          
 [70] checkmate_2.1.0         e1071_1.7-13            boot_1.3-28            
 [73] pkgbuild_1.4.0          truncnorm_1.0-9         systemfonts_1.0.4      
 [76] rlang_1.1.0             pkgconfig_2.0.3         matrixStats_1.0.0      
 [79] distributional_0.3.2    evaluate_0.20           labeling_0.4.2         
 [82] rstantools_2.4.0        htmlwidgets_1.6.2       bit_4.0.5              
 [85] processx_3.8.0          tidyselect_1.2.0        magrittr_2.0.3         
 [88] R6_2.5.1                generics_0.1.3          DBI_1.1.3              
 [91] pillar_1.8.1            withr_2.5.0             units_0.8-1            
 [94] abind_1.4-5             crayon_1.5.2            rARPACK_0.11-0         
 [97] uuid_1.1-0              interp_1.1-3            KernSmooth_2.23-20     
[100] utf8_1.2.3              tzdb_0.3.0              rmarkdown_2.21         
[103] jpeg_0.1-10             grid_4.2.2              callr_3.7.3            
[106] digest_0.6.31           classInt_0.4-8          webshot_0.5.4          
[109] brew_1.0-8              RcppParallel_5.1.6      stats4_4.2.2           
[112] munsell_0.5.0           viridisLite_0.4.1      

Reuse

https://creativecommons.org/licenses/by-sa/4.0/