Capa de sesgo
Patron-Rivero, C.
2026
Insumos
occ <- read.csv("D:/2026_ENM/Sesgo/rep_PY.csv", sep = ",")Kernel Utilization Distribution Probability
Regla de Silverman
\(h = \sigma n^{-1/6}\)
donde:
\(\sigma^2 = 0.5(\text{var}(x) + \text{var}(y))\)
El ancho de banda (h) es proporcional a la desviación estándar (σ) de los datos y varÃa inversamente con la raÃz sexta del tamaño de la muestra (n).
- Cuanto mayor sea la desviación estándar de los datos mayor será el valor de h.
- Cuanto mayor sea el tamaño de la muestra, menor será el valor h.
## PhantomJS not found. You can install it with webshot::install_phantomjs(). If it is installed, please make sure the phantomjs executable can be found via the PATH variable.Shiny applications not supported in static R Markdown documents
library(raster)
library(terra)
library(sp)
library(adehabitatHR)
library(sf)
var <- raster("D:/2026_ENM/Sesgo/bio1.tif")
r <- rasterToPoints(var)
r1 <- SpatialPoints(r[, 1:2])
r2 <- SpatialPixels(r1)
den <- occ[, c(1, 2)]
den <- SpatialPoints(den)
kde <- kernelUD(den, grid = r2, h = "href")
k <- rast(kde)
var_terra <- rast(var)
kr <- resample(k, var_terra)
b <- crop(kr, var_terra, mask = TRUE)Gausian Kernel Density / Rule-of-thumb
\(f(x,y) = \frac {\sum_{s} \phi (({x - x_s}) /{h_x}) \phi (({y - y_s}) /{h_y})}{nh_xh_y}\)
Donde:
ϕ = función de densidad en una distribución gausiana normal, no es negativa e integra a 1.
\(H = (Q3-Q1) / 1.34\)
\(IQR = (2 * z) * σ\)
\((2*0.674) ≈ \frac{IQR} {σ}\)
\(1.34 ≈ \frac{IQR} {σ}\)
Para el rango intercuartilico (del 25% al 75%) z = 0.674, Entonces si se supone una distribución normal, es una constante de 1.34
library(ntbox)
bl <- biaslayer(occs_df = occ, longitude = "decimalLongitude", latitude = "decimalLatitude", raster_mold = var)Geographic sampling biases - Bayesian framework
library(terra)
#install.packages("devtools")
#library(devtools)
#install_github("azizka/sampbias")
library(sampbias)
library(dplyr)
library(tidyr)
library(forcats)
library(viridis)
library(ggplot2)
library(rnaturalearth)
library(patchwork)
occ$species <- "Reptiles"
occ_b <- as.data.frame(occ[, c(4, 1, 2)])
var2 <- rast("D:/2026_ENM/Sesgo/bio1.tif")
wgs84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
terra::crs(var2) <- wgs84
res <- res(var2)[1]
air <- terra::vect("D:/2026_ENM/Sesgo/gaz/ne_10m_airports.shp")
air <- crop(air, var2)
air <- vect(geom(air)[, c("x", "y")])
cit <- terra::vect("D:/2026_ENM/Sesgo/gaz/ne_10m_urban_areas.shp")
cit <- crop(cit, var2)
cit <- vect(geom(cit)[, c("x", "y")])
roa <- terra::vect("D:/2026_ENM/Sesgo/gaz/ne_10m_roads.shp")
roa <- crop(roa, var2)
roa <- vect(geom(roa))
gazetteers <- list(airports = air, cities = cit, roads = roa)
time_taken <- system.time({
bay <- calculate_bias(x = occ_b, gaz = gazetteers, res = res, inp_raster = var2, buffer = 0)
})## |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= time_taken[[3]]/60## [1] 6.197833head(bay$bias_estimate)## it likA priorA q w_airports w_cities w_roads
## 1 20100 -45511.69 6.333344 5.037781 0.006024717 2.273047e-05 0.02176653
## 2 20200 -45513.77 6.271517 5.057486 0.005872948 2.081016e-05 0.02259387
## 3 20300 -45516.25 6.367092 4.751719 0.005513123 2.220263e-05 0.02108195
## 4 20400 -45511.83 6.269544 4.976132 0.005962226 1.905330e-05 0.02182307
## 5 20500 -45514.37 6.264357 4.901690 0.005870358 1.879697e-05 0.02198907
## 6 20600 -45511.65 6.284722 4.975259 0.006061241 1.980831e-05 0.02153435
## hp_rate
## 1 88.73541
## 2 124.25852
## 3 80.28490
## 4 80.28490
## 5 80.28490
## 6 80.28490plo1 <- bay$bias_estimate %>%
pivot_longer(cols = contains("w_"), names_to = "bias",
values_to = "posterior_estimate") %>%
mutate(bias = gsub("w_", "", .data$bias)) %>%
mutate(bias = fct_reorder(.data$bias, .data$posterior_estimate, .fun = median, .desc = FALSE))
head(plo1)## # A tibble: 6 × 7
## it likA priorA q hp_rate bias posterior_estimate
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <dbl>
## 1 20100 -45512. 6.33 5.04 88.7 airports 0.00602
## 2 20100 -45512. 6.33 5.04 88.7 cities 0.0000227
## 3 20100 -45512. 6.33 5.04 88.7 roads 0.0218
## 4 20200 -45514. 6.27 5.06 124. airports 0.00587
## 5 20200 -45514. 6.27 5.06 124. cities 0.0000208
## 6 20200 -45514. 6.27 5.06 124. roads 0.0226p1 <- ggplot(plo1)+
geom_boxplot(aes(x = .data$bias, y = .data$posterior_estimate, fill = .data$bias))+
scale_fill_viridis(discrete = TRUE)+
xlab("Biasing factor")+
ylab("Posterior weight")+
coord_flip()+
theme_bw()+
theme(panel.grid.minor.x = element_blank(),
panel.grid.major.y = element_blank(),
legend.position = "none")
p1
proj <- project_bias(bay)
result <- map_bias(proj, type = "log_sampling_rate")
proj <- raster(proj[[2]])
var <- crop(var, proj)
proj <- crop(proj, var)
pr <- mask(proj, var)Generar y guardar puntos de background
n depende de la cantidad de datos que necesites y del area
https://www.sciencedirect.com/science/article/pii/S0304380023003344
#library(remotes)
# remotes::install_github("sjevelazco/flexsdm")
library(flexsdm)
library(dplyr)
# b - Kernel utilization
# bl = Gausian RT
# pr = Bayesian
ncell(b)## [1] 11520bg_ku <- sample_background(
data = occ,
x = "decimalLongitude",
y = "decimalLatitude",
n = 1000,
method = "biased",
rlayer = var,
rbias = b
)
bg_kr <- sample_background(
data = occ,
x = "decimalLongitude",
y = "decimalLatitude",
n = 1000,
method = "biased",
rlayer = var,
rbias = bl
)
bg_b <- sample_background(
data = occ,
x = "decimalLongitude",
y = "decimalLatitude",
n = 1000,
method = "biased",
rlayer = var,
rbias = pr
)
Guardar los resultados
writeRaster(b, "D:/2026_ENM/Sesgo/ku_bias.tiff", format = "GTiff", overwrite = TRUE)
writeRaster(bl, "D:/2026_ENM/Sesgo/kr_bias.tiff", format = "GTiff", overwrite = TRUE)
writeRaster(pr, "D:/2026_ENM/Sesgo/bay_bias.tiff", format = "GTiff", overwrite = TRUE)
write.csv(bg_ku, "D:/2026_ENM/Sesgo/bg_ku.csv", row.names = FALSE)
write.csv(bg_kr, "D:/2026_ENM/Sesgo/bg_kr.csv", row.names = FALSE)
write.csv(bg_b, "D:/2026_ENM/Sesgo/bg_bay.csv", row.names = FALSE)