Taller Bioestadística

options(repos = c(CRAN = "https://cloud.r-project.org"))
knitr::opts_chunk$set(
  echo       = TRUE,
  warning    = FALSE,
  message    = FALSE,
  fig.width  = 9,
  fig.height = 6,
  fig.align  = "center",
  cache      = FALSE
)

# ── Instalación automática ────────────────────────────────────────────────────
paquetes <- c(
  "rgbif", "spatstat", "spatstat.geom", "spatstat.explore",
  "terra", "sf", "stars",
  "rnaturalearth", "rnaturalearthdata", "countrycode",
  "dplyr", "purrr","furrr", "ggplot2", "viridis", "scales", "gridExtra",
  "tidymodels", "corrplot", "broom",
  "palmerpenguins", "ggpubr", "emmeans", "rstatix",
  "knitr", "kableExtra", "mapview"
)

instalar_si_falta <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) install.packages(pkg)
}
invisible(lapply(paquetes, instalar_si_falta))

suppressPackageStartupMessages({
  library(rgbif);          library(spatstat)
  library(spatstat.geom);  library(spatstat.explore)
  library(terra);          library(sf);           library(stars)
  library(rnaturalearth);  library(rnaturalearthdata)
  library(countrycode);    library(dplyr)
  library(purrr);          library(ggplot2)
  library(viridis);        library(scales)
  library(gridExtra);      library(tidymodels)
  library(corrplot);       library(broom)
  library(palmerpenguins); library(ggpubr)
  library(emmeans);        library(rstatix)
  library(knitr);          library(kableExtra)
  library(mapview)
})

# ═══════════════════════════════════════════════════════════════════════
# >>> ESTUDIANTE: modifique SOLO este bloque con sus especies y país
# ═══════════════════════════════════════════════════════════════════════

MI_PAIS <- "CO"

mis_especies <- c(
  "Bothrops asper",
  "Bothrops atrox",
  "Bothrops punctatus",
  "Lachesis muta",
  "Crotalus durissus",
  "Micrurus mipartitus",
  "Micrurus dumerilii",
  "Boa constrictor",
  "Leptodeira annulata",
  "Spilotes pullatus"
)

SIGMA_KDE <- 1.0

# ═══════════════════════════════════════════════════════════════════════

cat("País:", MI_PAIS, "\n")
cat("Especies ingresadas:\n")
for (i in seq_along(mis_especies)) cat(sprintf("  %2d. %s\n", i, mis_especies[i]))

# ── Función de descarga con manejo de errores ─────────────────────────────────
descargar_especie <- function(nombre, pais, n_max = 500) {
  tryCatch({
    res <- occ_search(
      scientificName     = nombre,
      country            = pais,
      hasCoordinate      = TRUE,
      hasGeospatialIssue = FALSE,
      limit              = n_max,
      fields             = c("decimalLongitude", "decimalLatitude",
                             "year", "basisOfRecord", "species")
    )
    df <- res$data
    if (is.null(df) || nrow(df) == 0) {
      message("  Sin registros: ", nombre); return(NULL)
    }
    df %>%
      filter(!is.na(decimalLongitude), !is.na(decimalLatitude)) %>%
      rename(lon = decimalLongitude, lat = decimalLatitude) %>%
      mutate(especie = nombre)
  }, error = function(e) {
    message("  Error en ", nombre, ": ", conditionMessage(e)); NULL
  })
}

lista_occ <- list()

for (sp in mis_especies) {
  message("Descargando: ", sp)
  df <- descargar_especie(sp, MI_PAIS)
  if (!is.null(df) && nrow(df) > 0) {
    lista_occ[[sp]] <- df
  } else {
    message("  Sin datos: ", sp)
  }
  Sys.sleep(1.2)
}

ocurrencias <- bind_rows(lista_occ) %>%
  filter(especie %in% mis_especies)

resumen_desc <- ocurrencias %>%
  count(especie, name = "n_registros") %>%
  arrange(desc(n_registros))

kable(resumen_desc,
      caption = "Registros de presencia descargados por especie") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

cat(sprintf("\nTotal de registros: %d\n", nrow(ocurrencias)))

nombre_pais <- countrycode(MI_PAIS, "iso2c", "country.name")
pais_sf     <- ne_countries(country = nombre_pais, scale = "medium",
                             returnclass = "sf")

ggplot() +
  geom_sf(data = pais_sf, fill = "#f0ede4", color = "gray35", linewidth = 0.4) +
  geom_point(
    data = ocurrencias %>%
      filter(especie %in% mis_especies) %>%
      filter(!is.na(lon), !is.na(lat)),
    aes(x = lon, y = lat, color = especie),
    size = 1.5, alpha = 0.65
  ) +
  scale_color_viridis_d(option = "turbo", name = "Especie") +
  labs(
    title    = paste("Registros de presencia —", nombre_pais),
    subtitle = paste("Total registros:", nrow(ocurrencias)),
    x = "Longitud", y = "Latitud"
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.text = element_text(face = "italic", size = 8))

# ── Construcción robusta de la ventana owin ───────────────────────────────────
pais_union  <- pais_sf %>% st_union() %>% st_simplify(dTolerance = 0.05)

pais_partes <- st_cast(pais_union, "POLYGON")
areas       <- st_area(pais_partes)
pais_main   <- pais_partes[which.max(areas)]

coords_ext  <- st_coordinates(pais_main)
anillo_ext  <- coords_ext[coords_ext[, "L2"] == 1, c("X", "Y")]

n_pts <- nrow(anillo_ext)
if (all(anillo_ext[1, ] == anillo_ext[n_pts, ])) {
  anillo_ext <- anillo_ext[-n_pts, ]
}

ventana_pais <- tryCatch({
  owin(poly = list(x = anillo_ext[, "X"],
                   y = anillo_ext[, "Y"]))
}, error = function(e) {
  message("Polígono complejo — usando ventana rectangular como respaldo")
  bb <- st_bbox(pais_sf)
  owin(xrange = c(bb["xmin"], bb["xmax"]),
       yrange = c(bb["ymin"], bb["ymax"]))
})

cat("Ventana creada exitosamente\n")
cat("Tipo:", ventana_pais$type, "\n")
cat("Rango X (lon):", paste(round(ventana_pais$xrange, 2), collapse = " – "), "\n")
cat("Rango Y (lat):", paste(round(ventana_pais$yrange, 2), collapse = " – "), "\n")

# ── Función: construir ppp y estimar intensidad KDE ───────────────────────────
analizar_ppp <- function(df_sp, ventana, sigma) {
  nombre <- unique(df_sp$especie)[1]

  pts <- df_sp %>%
    filter(lon >= ventana$xrange[1], lon <= ventana$xrange[2],
           lat >= ventana$yrange[1], lat <= ventana$yrange[2])

  if (nrow(pts) < 10) {
    message(sprintf("  Omitida (<%d pts): %s", 10, nombre)); return(NULL)
  }

  ppp_obj <- tryCatch(
    ppp(x = pts$lon, y = pts$lat, window = ventana),
    error = function(e) {
      message("  Error ppp en ", nombre, ": ", conditionMessage(e)); NULL
    }
  )
  if (is.null(ppp_obj)) return(NULL)

  ppp_obj <- unique.ppp(ppp_obj)

  kde <- density.ppp(ppp_obj, sigma = sigma, eps = 0.08)

  message(sprintf("  OK %-38s | n = %d", nombre, npoints(ppp_obj)))

  list(nombre = nombre, n_pts = npoints(ppp_obj), ppp = ppp_obj, kde = kde)
}

resultados_ppp <- ocurrencias %>%
  group_by(especie) %>%
  group_split() %>%
  map(~ analizar_ppp(.x, ventana_pais, sigma = SIGMA_KDE)) %>%
  compact()

cat(sprintf("\nEspecies procesadas exitosamente: %d / %d\n",
            length(resultados_ppp), length(mis_especies)))

n_sp  <- length(resultados_ppp)

if (n_sp > 0) {
  n_col <- min(3, n_sp)
  n_row <- ceiling(n_sp / n_col)

  par(mfrow = c(n_row, n_col), mar = c(1.2, 1.2, 3.2, 0.8))

  for (r in resultados_ppp) {
    plot(r$kde,
         main     = paste0(r$nombre, "  (n=", r$n_pts, ")"),
         col      = hcl.colors(64, "Inferno"),
         las      = 1,
         ribside  = "right",
         cex.main = 0.82)
    plot(ventana_pais, add = TRUE, lwd = 0.9, border = "white")
    points(r$ppp, pch = ".", cex = 1.6, col = "white")
  }

  par(mfrow = c(1, 1))
}

if (length(resultados_ppp) > 0) {
  tabla_int <- map_dfr(resultados_ppp, function(r) {
    vals <- as.vector(as.matrix(r$kde))
    vals <- vals[!is.na(vals) & vals > 0]
    tibble(
      Especie      = r$nombre,
      N_registros  = r$n_pts,
      lambda_media = round(mean(vals), 6),
      lambda_max   = round(max(vals), 6),
      lambda_p75   = round(quantile(vals, 0.75), 6)
    )
  })

  kable(tabla_int,
        caption = "Estadísticos de intensidad estimada λ̂(x,y) por especie",
        col.names = c("Especie", "n", "λ media", "λ máximo", "λ p75")) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                  full_width = FALSE)
}

kde_a_raster_bin <- function(r_obj) {
  kde <- r_obj$kde
  mat <- as.matrix(kde)

  rst <- rast(
    nrows = nrow(mat), ncols = ncol(mat),
    xmin  = min(kde$xcol), xmax = max(kde$xcol),
    ymin  = min(kde$yrow), ymax = max(kde$yrow),
    crs   = "EPSG:4326"
  )
  values(rst) <- as.vector(t(mat))

  v <- values(rst)
  umbral <- median(v[!is.na(v) & v > 0])
  ifel(rst > umbral, 1, 0)
}

pais_vect <- terra::vect(pais_sf)

rasters_bin <- map(resultados_ppp, function(r) {
  rst <- kde_a_raster_bin(r)
  rst <- terra::crop(rst, pais_vect)
  rst <- terra::mask(rst, pais_vect)
  return(rst)
})

if (length(rasters_bin) > 0) {
  rst_ref <- rasters_bin[[1]]
  riqueza <- rst_ref

  for (i in 2:length(rasters_bin)) {
    rst_i   <- resample(rasters_bin[[i]], rst_ref, method = "near")
    riqueza <- riqueza + rst_i
  }

  riqueza <- terra::crop(riqueza, pais_vect)
  riqueza <- terra::mask(riqueza, pais_vect)

  cat("Riqueza mínima:", global(riqueza, "min", na.rm = TRUE)[[1]], "\n")
  cat("Riqueza máxima:", global(riqueza, "max", na.rm = TRUE)[[1]], "especies\n")
}

## Riqueza mínima: 0 
## Riqueza máxima: 10 especies

if (exists("riqueza")) {
  pal_azul_disc <- c("#deebf7", "#9ecae1", "#6baed6", "#3182bd", "#08519c")

  vals <- terra::values(riqueza)
  vals <- vals[!is.na(vals)]

  cortes <- seq(0, max(vals), by = 2)

  mapview(
    riqueza,
    col.regions = pal_azul_disc,
    at = cortes,
    alpha.regions = 0.8,
    legend = TRUE,
    layer.name = paste("Riqueza -", nombre_pais)
  ) + 
  mapview(
    pais_sf,
    color = "#4a6fa5",
    alpha.regions = 0
  )
}

Taller Bioestadística

Patrones puntuales · Regresión predictiva · Diseño de experimentos

Sebastián Hidalgo y María Camila Chaves

23 de abril de 2026