1 Módulo 1 — Análisis de patrones de puntos y riqueza de especies

1.1 Contexto ecológico

Los patrones de puntos espaciales (point patterns) permiten analizar si las ocurrencias de una especie se distribuyen de forma aleatoria, agrupada o regular en el espacio. La intensidad \(\lambda(x,y)\) estima la densidad esperada de individuos por unidad de área en cada ubicación.

A partir de las intensidades individuales de múltiples especies podemos construir un mapa de riqueza estimada sumando rasters de presencia/ausencia, lo que es una aproximación espacial al concepto de acumulación de diversidad.

1.2 Paso 1 — Elección de especies y país

Instrucción: Complete el bloque siguiente con sus 10 especies y el código ISO-2 de su país. Use nombres científicos válidos (Género + epíteto específico). Puede verificar nombres en gbif.org/species/search o usando name_backbone() en el Paso 2.

# Código ISO-2 del país donde buscará los registros
# Ejemplos: "CO" Colombia | "BR" Brasil | "MX" México |
#           "PE" Perú     | "EC" Ecuador| "VE" Venezuela
MI_PAIS <- "CO"

# Sus 10 especies (nombre científico completo: Género especie)
mis_especies <- c(
  "Panthera onca",
  "Saguinus oedipus",
  "Tremarctos ornatus",
  "Bolitoglossa walkeri",
  "Phyllobates terribilis",
  "Minyobates bombetes",
  "Ensifera ensifera",
  "Chlorochrysa nitidissima",
  "Hyalinobatrachium fleischmanni",
  "Oophaga lehmanni"
)

# Ancho de banda del kernel (en grados). Orientación:
#   países pequeños (Ecuador, El Salvador): sigma = 0.3 – 0.5
#   países medianos (Colombia, Venezuela):  sigma = 0.8 – 1.2
#   países grandes (Brasil, México):        sigma = 1.5 – 2.0
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]))

Se definen los parámetros principales del análisis: el país de interés (Colombia) y las 10 especies seleccionadas con base en su representatividad ecológica y disponibilidad de registros en GBIF; Además del ancho de banda σ = 1.0° para la estimación de densidad kernel (KDE).

1.3 Paso 2 — Verificación taxonómica en GBIF

# ── Consultar el backbone taxonómico de GBIF ──────────────────────────────────
verificar_sp <- function(nombre) {
  res <- tryCatch(
    name_backbone(name = nombre, rank = "species"),
    error = function(e) NULL
  )
  if (is.null(res) || nrow(res) == 0)
    return(tibble(ingresado = nombre, aceptado = NA, rango = NA,
                  confianza = NA, estado = "NO ENCONTRADO"))
  tibble(
    ingresado  = nombre,
    aceptado   = res$canonicalName[1],
    rango      = res$rank[1],
    confianza  = res$confidence[1],
    estado     = ifelse(res$rank[1] == "SPECIES", "OK", "REVISAR")
  )
}

tabla_tax <- map_dfr(mis_especies, verificar_sp)

kable(tabla_tax,
      caption = "Verificación taxonómica — backbone GBIF") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  column_spec(5, color = ifelse(tabla_tax$estado == "OK",
                                 "darkgreen", "red"),
              bold = TRUE)

Verificación taxonómica — backbone GBIF
ingresado	aceptado	rango	confianza	estado
Panthera onca	Panthera onca	SPECIES	99	OK
Saguinus oedipus	Saguinus oedipus	SPECIES	99	OK
Tremarctos ornatus	Tremarctos ornatus	SPECIES	99	OK
Bolitoglossa walkeri	Bolitoglossa walkeri	SPECIES	99	OK
Phyllobates terribilis	Phyllobates terribilis	SPECIES	99	OK
Minyobates bombetes	Minyobates bombetes	SPECIES	98	OK
Ensifera ensifera	Ensifera ensifera	SPECIES	99	OK
Chlorochrysa nitidissima	Chlorochrysa nitidissima	SPECIES	99	OK
Hyalinobatrachium fleischmanni	Hyalinobatrachium fleischmanni	SPECIES	99	OK
Oophaga lehmanni	Oophaga lehmanni	SPECIES	99	OK

Las 10 especies fueron verificadas contra el backbone taxonómico de GBIF mediante la función name_backbone(). Todas obtuvieron estado OK con un nivel de confianza de 98–99%, confirmando que los nombres científicos ingresados son válidos y reconocidos por la base de datos global.

Acción: Si el estado es REVISAR o NO ENCONTRADO, corrija el nombre en el Paso 1 y vuelva a ejecutar este bloque. Continúe cuando todas las especies aparezcan en verde.

1.4 Paso 3 — Descarga de ocurrencias desde GBIF

# ── 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
  })
}

# Descarga con pausa cortés (1.2 s) entre peticiones
message("Iniciando descarga...")
lista_occ <- vector("list", length(mis_especies))
for (i in seq_along(mis_especies)) {
  message(sprintf("  [%d/%d] %s", i, length(mis_especies), mis_especies[i]))
  lista_occ[[i]] <- descargar_especie(mis_especies[i], MI_PAIS)
  Sys.sleep(1.2)
}

ocurrencias <- bind_rows(lista_occ)

# Resumen de registros
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) %>%
  column_spec(2, color = ifelse(resumen_desc$n_registros < 10,
                                 "red", "black"))

Registros de presencia descargados por especie
especie	n_registros
Chlorochrysa nitidissima	500
Ensifera ensifera	500
Minyobates bombetes	500
Oophaga lehmanni	500
Panthera onca	500
Saguinus oedipus	500
Tremarctos ornatus	500
Bolitoglossa walkeri	408
Phyllobates terribilis	312
Hyalinobatrachium fleischmanni	257

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

## 
## Total de registros: 4477

Se descargaron hasta 500 registros por especie desde GBIF con coordenadas válidas para Colombia. Siete de las diez especies alcanzaron el límite máximo de 500 registros, lo que indica una buena representación en la base de datos.

Todas las especies superan el mínimo de 10 registros requeridos para el análisis de patrones de puntos.

Atención: Las especies con menos de 10 registros (marcadas en rojo) serán excluidas del análisis de patrones de puntos. Si varias especies tienen pocos registros, considere elegir otras con mayor cobertura en GBIF para su país.

1.5 Paso 4 — Mapa de registros crudos

# ── Polígono del país ─────────────────────────────────────────────────────────
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,
             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 = "Fuente: GBIF (rgbif)",
    x = "Longitud", y = "Latitud"
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.text = element_text(face = "italic", size = 8))

Como e puede observar en el mapa, este muestra la distribución geográfica de los registros de presencia de las 10 especies en Colombia. Se observa una concentración notable en la región andina y el noroccidente del país, particularmente entre los 4°N y 8°N. Especies como Panthera onca presentan registros más dispersos a lo largo del territorio, mientras que anfibios como Phyllobates terribilis y Oophaga lehmanni muestran distribuciones más restringidas, coherentes con su condición de endemismos del Chocó biogeográfico.

1.6 Paso 5 — Construcción de la ventana espacial (`owin`)

El objeto owin (observation window) de spatstat define el área de estudio. La forma más robusta de construirlo es extraer directamente las coordenadas del polígono sf y pasarlas como lista de vértices, lo que evita el error “cannot interpret W as a window” que ocurre cuando as.owin() no puede interpretar el objeto Spatial de versiones recientes de sf.

# ── Construcción robusta de la ventana owin ───────────────────────────────────
# 1. Simplificar y disolver el polígono del país
pais_union  <- pais_sf %>% st_union() %>% st_simplify(dTolerance = 0.05)

# 2. Seleccionar el polígono más grande (continente, ignorar islas menores)
pais_partes <- st_cast(pais_union, "POLYGON")
areas       <- st_area(pais_partes)
pais_main   <- pais_partes[which.max(areas)]

# 3. Extraer coordenadas del anillo exterior
coords_ext  <- st_coordinates(pais_main)
# st_coordinates devuelve columnas X, Y, L1, L2 — tomamos solo el anillo 1
anillo_ext  <- coords_ext[coords_ext[, "L2"] == 1, c("X", "Y")]

# 4. Construir la ventana como polígono explícito
#    spatstat requiere que el polígono NO esté cerrado (sin repetir el primer punto)
n_pts <- nrow(anillo_ext)
if (all(anillo_ext[1, ] == anillo_ext[n_pts, ])) {
  anillo_ext <- anillo_ext[-n_pts, ]     # eliminar el punto de cierre
}

ventana_pais <- tryCatch({
  owin(poly = list(x = anillo_ext[, "X"],
                   y = anillo_ext[, "Y"]))
}, error = function(e) {
  # Respaldo: ventana rectangular (bbox del país)
  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")

Se construyó exitosamente la ventana espacial de análisis con forma poligonal, siguiendo el contorno real de Colombia. La ventana cubre un rango longitudinal de -79.03° a -66.88° (eje X) y latitudinal de -4.24° a 12.43° (eje Y), abarcando todo el territorio continental colombiano. Esta ventana poligonal es más precisa que una ventana rectangular, ya que restringe el análisis de densidad kernel únicamente al área terrestre del país, evitando estimaciones en zonas marítimas o territorios vecinos.

1.7 Paso 6 — Análisis de patrones de puntos por especie

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

  # Filtrar al bbox de la ventana
  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)
  }

  # Construir patrón de puntos
  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)

  # Eliminar duplicados exactos
  ppp_obj <- unique.ppp(ppp_obj)

  # Densidad de kernel Gaussiano
  # eps controla la resolución de la grilla de salida (en grados)
  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)
}

message("\nCalculando KDE por especie (sigma = ", SIGMA_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)))

Las 10 especies fueron procesadas exitosamente para el análisis de patrones de puntos. Los registros únicos (sin duplicados exactos) varían entre 21 para Phyllobates terribilis y 339 para Panthera onca. La baja cantidad de registros de Phyllobates terribilis (n = 21) refleja su distribución extremadamente restringida al sur del Chocó, mientras que Panthera onca y Saguinus oedipus presentan la mayor cobertura geográfica con 339 y 277 registros respectivamente. Todas las especies superaron el umbral mínimo de 10 puntos requerido para estimar la densidad kernel, por lo que ninguna fue excluida del análisis.

1.8 Paso 7 — Visualización de densidades KDE

n_sp  <- length(resultados_ppp)
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))

Los mapas de densidad kernel (KDE) muestran la intensidad estimada λ̂(x,y) de cada especie en Colombia con σ = 1.0°. Las zonas amarillas indican alta densidad de registros y las oscuras ausencia. La mayoría de especies presentan patrones agrupados concentrados en la región andina y el Chocó biogeográfico, lo que sugiere distribuciones no aleatorias asociadas a condiciones ambientales específicas.

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)

Estadísticos de intensidad estimada λ̂(x,y) por especie
Especie	n	λ media	λ máximo	λ p75
Bolitoglossa walkeri	84	0.901915	12.166735	0.320253
Chlorochrysa nitidissima	122	1.342173	19.092415	0.563710
Ensifera ensifera	162	1.814464	16.526584	1.654621
Hyalinobatrachium fleischmanni	121	1.337560	10.313733	1.873127
Minyobates bombetes	237	2.635576	35.746707	0.810842
Oophaga lehmanni	197	2.040807	47.168919	0.162982
Panthera onca	339	3.612844	56.459485	3.224255
Phyllobates terribilis	21	0.203606	6.154213	0.002327
Saguinus oedipus	277	2.856986	42.538398	0.417068
Tremarctos ornatus	205	2.315739	14.414257	4.027909

La tabla presenta tres estadísticos de la intensidad estimada λ̂(x,y) por especie. Panthera onca registra la mayor intensidad media (λ = 3.61) y máxima (λ = 56.46), coherente con sus 339 registros distribuidos ampliamente por el territorio. En contraste, Phyllobates terribilis presenta los valores más bajos en todos los estadísticos (λ media = 0.20, λ máximo = 6.15), reflejando su distribución extremadamente restringida al Chocó con solo 21 registros.

El percentil 75 (λ p75) indica que el 75% del área de Colombia tiene una intensidad menor a ese valor para cada especie. Valores muy bajos de λ p75 como los de Phyllobates terribilis (0.002) y Oophaga lehmanni (0.163) indican que la especie ocupa una fracción muy pequeña del territorio, con alta densidad concentrada en pocas celdas. Por el contrario, Tremarctos ornatus (λ p75 = 4.03) y Panthera onca (λ p75 = 3.22) muestran una distribución más homogénea a lo largo del país.

1.9 Paso 8 — Riqueza estimada por álgebra de rasters

La riqueza en cada celda se define como:

\[R(x,y) = \sum_{i=1}^{S} \mathbf{1}\left[\hat{\lambda}_i(x,y) > p_{50,i}\right]\]

donde \(p_{50,i}\) es la mediana de la intensidad estimada para la especie \(i\). Una especie “ocupa” una celda si su densidad local supera la mediana de su distribución de densidad.

# ── KDE → SpatRaster binario ──────────────────────────────────────────────────
kde_a_raster_bin <- function(r_obj) {
  kde <- r_obj$kde
  mat <- as.matrix(kde)          # filas = y, cols = x

  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))   # transponer para alinear ejes x/y

  # Umbral = mediana de la densidad positiva
  v <- values(rst)
  umbral <- median(v[!is.na(v) & v > 0])
  ifel(rst > umbral, 1, 0)
}

rasters_bin <- map(resultados_ppp, kde_a_raster_bin)

# Sumar rasters (álgebra): raster de referencia = primera especie
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
}

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")

library(leaflet)
library(terra)
library(sf)

# Restaurado para que los chunks siguientes lo encuentren
riqueza_df <- as.data.frame(riqueza, xy = TRUE) %>%
  setNames(c("x", "y", "riqueza")) %>%
  filter(!is.na(riqueza), riqueza > 0)

# Forzar carga de valores en memoria
riqueza_mem <- rast(riqueza)
values(riqueza_mem) <- values(riqueza)

# Recortar raster a la silueta de Colombia
riqueza_clip <- mask(riqueza_mem, vect(pais_sf))

# Dominio calculado desde el raster recortado
rango <- range(values(riqueza_clip), na.rm = TRUE)

pal <- colorNumeric(
  palette  = c("#2166ac", "#67a9cf", "#f7f7f7", "#f4a582", "#d6604d"),
  domain   = rango,
  na.color = "transparent"
)

leaflet() %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addRasterImage(
    riqueza_clip,
    colors  = pal,
    opacity = 0.85
  ) %>%
  addPolylines(
    data    = pais_sf,
    color   = "white",
    weight  = 1.5,
    opacity = 0.9
  ) %>%
  addLegend(
    position = "bottomright",
    pal      = pal,
    values   = rango,
    title    = "Riqueza<br>estimada",
    opacity  = 0.9
  ) %>%
  addControl(
    html = paste0(
      "<div style='background:rgba(0,0,0,0.6);color:white;padding:6px 10px;",
      "border-radius:4px;font-size:12px;max-width:280px'>",
      "<b>Riqueza estimada — ", nombre_pais, "</b><br>",
      length(resultados_ppp), " especies · KDE (σ = ", SIGMA_KDE, "°)",
      "</div>"
    ),
    position = "topleft"
  )

El mapa interactivo muestra la riqueza estimada de especies por celda, calculada como la suma de los rasters binarios de las 10 especies (escala 0–10). Las zonas en rojo intenso representan alta coocurrencia de especies (hasta 10), concentradas principalmente en la región andina occidental y el Chocó biogeográfico, que coincide con uno de los hotspots de biodiversidad más importantes del mundo.

ggplot(riqueza_df, aes(x = riqueza)) +
  geom_histogram(aes(fill = after_stat(x)), binwidth = 1,
                 color = "white", linewidth = 0.3) +
  scale_fill_viridis_c(option = "plasma", guide = "none") +
  scale_x_continuous(breaks = 0:length(resultados_ppp)) +
  labs(
    title   = "Distribución de riqueza estimada",
    x       = "Número de especies coocurrentes",
    y       = "Celdas del raster"
  ) +
  theme_minimal(base_size = 12)

La distribución de riqueza estimada revela un patrón interesante: la mayoría de las celdas del raster presentan either muy baja (1–3 especies) o muy alta coocurrencia (9–10 especies), con relativamente pocas celdas en valores intermedios. Las celdas con 9 y 10 especies coocurrentes son las más abundantes del territorio, lo que sugiere que las 10 especies seleccionadas comparten en gran medida las mismas zonas de alta idoneidad ambiental, particularmente en la región andina y el Chocó biogeográfico. Por otro lado, las celdas con 1 a 3 especies corresponden a zonas de transición o áreas donde solo las especies de distribución más amplia como Panthera onca logran estar presentes.

1.9.1 Preguntas de interpretación — Módulo 1

¿Qué especies presentaron mayor intensidad máxima \(\hat{\lambda}_{max}\)? ¿Podría deberse a la ecología de la especie o al esfuerzo de muestreo?

Las tres especies con mayor λ_max fueron Panthera onca (56.46), Oophaga lehmanni (47.17) y Saguinus oedipus(42.51). Como se puede observar en la tabla de intensidad estimada, Panthera onca y Saguinus oedipus tienen los dos n más altos de la tabla (338 y 272). Esto sugiere que su λ_max elevado está influenciado principalmente por el esfuerzo de muestreo, ya que los datos de GBIF presentan sesgos espaciales derivados de un muestreo desigual (Beck et al., 2014). Al haber más registros, aumenta la probabilidad de acumulación en ciertas zonas, generando picos artificialmente altos en el mapa de calor. Sin embargo, Oophaga lehmanni posee solo 197 registros y alcanza un λ_max de 47.17, pero su λ p75 es apenas 0.16, siendo el más bajo de toda la tabla junto con Phyllobates terribilis. Esto significa que el pico enorme ocurre en un área geográfica muy pequeña. Esto es consistente con su ecología real, es una rana endémica de distribución muy restringida en el Pacífico colombiano, por lo que todos sus registros naturalmente se concentran en una zona pequeña (WCS - Oophaga Lehmanni, 2025).

¿Los hot-spots del mapa de riqueza coinciden con áreas de alta biodiversidad conocidas en su país?

Los hotspots del mapa están representados en rojo intenso con riqueza estimada cercana a 10 especies, las cuales se concentran en el occidente colombiano y coinciden con el Chocó biogeográfico y la vertiente occidental de los Andes, dos de las regiones de mayor biodiversidad del planeta (Herrera et al., 2018; Bax & Francesconi, 2018). Esto es ecológicamente coherente dado que varias de las especies seleccionadas son endémicas o tienen distribución preferencial en esa región, como Oophaga lehmanni y Phyllobates terribilis (Casagranda & Goloboff, 2019; IUCN, 2020).

En contraste, el oriente del país el cual corresponde a los Llanos y la Amazonia, aparece en azul con riqueza cercana a cero, esto responde a dos factores combinados. Por un lado, un sesgo de muestreo en GBIF hacia regiones más accesibles; por otro, que las especies elegidas son predominantemente andinas y del Pacífico, por lo que el mapa refleja la distribución de ese grupo específico y no la biodiversidad total del país.

¿Qué limitaciones tiene usar registros de GBIF para análisis de patrones de puntos?

El uso de registros de GBIF para analizar patrones de puntos presenta varias limitaciones. Una de las principales es el sesgo espacial de muestreo, ya que los datos se concentran en zonas accesibles o ampliamente estudiadas; esto explica por qué regiones como el oriente colombiano aparecen casi vacías en el mapa de riqueza, no por falta de biodiversidad, sino por escasez de registros. Aunque GBIF reúne millones de datos de colecciones y observaciones que permiten analizar la distribución de las especies (López Tobar, 2025), su cobertura no es uniforme. Además, existe un sesgo temporal, pues muchos registros son antiguos y pueden no reflejar la distribución actual, y un sesgo de detección por especie, donde especies más visibles como Panthera onca tienen muchos más registros que especies crípticas como Phyllobates terribilis, lo que dificulta comparaciones directas.

¿Cómo cambiaría el mapa si usara sigma = 0.3 en lugar de su valor actual? ¿Y si usara sigma = 3.0?

Con sigma = 0.3 cada especie muestra manchas pequeñas y concentradas justo donde están los registros reales, lo que produce un mapa más preciso pero con grandes zonas vacías. Con sigma = 3.0 las manchas se expanden considerablemente cubriendo regiones más amplias, lo que suaviza el mapa pero puede sobreestimar el área de distribución real de especies de rango restringido como Oophaga lehmanni o Phyllobates terribilis. El valor de sigma = 1.0 representa un balance adecuado para Colombia, evitando tanto la fragmentación excesiva del patrón como la sobreestimación del área ocupada.

2 Módulo 2 — Regresión lineal múltiple con enfoque predictivo

2.1 Contexto biológico

Usaremos datos morfométricos de tortugas pintadas (Chrysemys picta), un modelo clásico en ecología de reptiles. El objetivo es predecir la longitud de la concha a partir de variables morfológicas y ambientales.

La diferencia entre enfoque explicativo y predictivo es central en bioestadística moderna: en el explicativo nos interesan los coeficientes e inferencia; en el predictivo nos interesa el desempeño en datos nuevos no vistos.

2.2 Datos

# ── Datos morfométricos basados en relaciones alométricas reales ──────────────
# Referencia: Janzen et al. (1992) Ecology
set.seed(42)
n <- 280

datos_tortugas <- tibble(
  longitud_concha_mm = rnorm(n, mean = 138, sd = 22),
  sexo               = sample(c("Hembra", "Macho"), n, replace = TRUE,
                               prob = c(0.55, 0.45))
) %>%
  mutate(
    longitud_concha_mm = longitud_concha_mm + ifelse(sexo == "Hembra", 8, 0),
    peso_g          = 0.0018 * longitud_concha_mm^2.7 +
                      ifelse(sexo == "Hembra", 60, 0) + rnorm(n, 0, 35),
    ancho_concha_mm = 0.82 * longitud_concha_mm + rnorm(n, 0, 8),
    altura_concha_mm= 0.45 * longitud_concha_mm + rnorm(n, 0, 6),
    largo_cabeza_mm = 0.28 * longitud_concha_mm + rnorm(n, 0, 5),
    temp_agua_c     = runif(n, 12, 28),
    profundidad_m   = runif(n, 0.2, 3.5),
    edad_años       = pmax(round(longitud_concha_mm/18 + rnorm(n, 0, 1.5)), 1)
  )

glimpse(datos_tortugas)

## Rows: 280
## Columns: 9
## $ longitud_concha_mm <dbl> 168.16109, 133.57664, 153.98883, 159.92298, 154.893…
## $ sexo               <chr> "Macho", "Hembra", "Hembra", "Hembra", "Hembra", "H…
## $ peso_g             <dbl> 1867.4922, 997.0991, 1513.8624, 1645.5372, 1564.640…
## $ ancho_concha_mm    <dbl> 140.24963, 112.67478, 118.26409, 128.53102, 118.946…
## $ altura_concha_mm   <dbl> 80.23671, 57.76213, 75.35438, 68.20350, 59.89479, 7…
## $ largo_cabeza_mm    <dbl> 56.97275, 34.32343, 50.71398, 42.61396, 38.21520, 3…
## $ temp_agua_c        <dbl> 13.69241, 24.35050, 22.17509, 21.47685, 12.18190, 1…
## $ profundidad_m      <dbl> 2.9048596, 0.7727506, 0.5567349, 1.5877060, 0.42013…
## $ edad_años          <dbl> 10, 6, 7, 8, 9, 8, 8, 10, 11, 9, 9, 11, 8, 9, 9, 6,…

2.3 Análisis exploratorio

p1 <- ggplot(datos_tortugas, aes(x = longitud_concha_mm, fill = sexo)) +
  geom_histogram(bins = 25, color = "white", alpha = 0.85) +
  scale_fill_manual(values = c("#3F51B5","#E91E63")) +
  labs(title = "Variable respuesta: longitud de concha",
       x = "Longitud (mm)", y = "Frecuencia") + theme_minimal()

p2 <- ggplot(datos_tortugas, aes(x = peso_g, y = longitud_concha_mm,
                                   color = sexo)) +
  geom_point(alpha = 0.45, size = 1.6) +
  geom_smooth(method = "lm", se = TRUE, linewidth = 1) +
  scale_color_manual(values = c("#3F51B5","#E91E63")) +
  labs(title = "Longitud vs Peso", x = "Peso (g)",
       y = "Longitud concha (mm)") + theme_minimal()

p3 <- ggplot(datos_tortugas, aes(x = temp_agua_c, y = longitud_concha_mm,
                                   color = sexo)) +
  geom_point(alpha = 0.45, size = 1.6) +
  geom_smooth(method = "loess", se = TRUE, color = "gray30", linewidth = 1) +
  scale_color_manual(values = c("#3F51B5","#E91E63")) +
  labs(title = "Longitud vs Temperatura", x = "Temperatura (°C)",
       y = "Longitud concha (mm)") + theme_minimal()

p4 <- ggplot(datos_tortugas, aes(x = sexo, y = longitud_concha_mm,
                                   fill = sexo)) +
  geom_violin(alpha = 0.35) +
  geom_boxplot(width = 0.18, alpha = 0.85) +
  scale_fill_manual(values = c("#3F51B5","#E91E63"), guide = "none") +
  labs(title = "Longitud por sexo", x = NULL,
       y = "Longitud concha (mm)") + theme_minimal()

grid.arrange(p1, p2, p3, p4, ncol = 2)

vars_num <- datos_tortugas %>%
  select(longitud_concha_mm, peso_g, ancho_concha_mm, altura_concha_mm,
         largo_cabeza_mm, temp_agua_c, profundidad_m, edad_años)

corrplot(cor(vars_num, use = "complete.obs"),
         method = "color", type = "upper",
         addCoef.col = "black", number.cex = 0.72,
         tl.cex = 0.9, tl.col = "gray20",
         col = COL2("PuOr", 10),
         title = "Correlaciones entre variables morfométricas",
         mar = c(0, 0, 2, 0))

2.4 División train / test y preprocesamiento

set.seed(123)

split <- initial_split(datos_tortugas, prop = 0.75, strata = sexo)
train <- training(split);  test <- testing(split)

cat(sprintf("Entrenamiento: %d | Prueba: %d\n", nrow(train), nrow(test)))

receta <- recipe(longitud_concha_mm ~ ., data = train) %>%
  step_dummy(all_nominal_predictors()) %>%
  step_normalize(all_numeric_predictors()) %>%
  step_corr(all_numeric_predictors(), threshold = 0.90)

receta_prep <- prep(receta, training = train)
train_proc  <- bake(receta_prep, new_data = NULL)
test_proc   <- bake(receta_prep, new_data = test)

cat("Predictores tras preproceso:",
    paste(setdiff(names(train_proc), "longitud_concha_mm"), collapse = ", "), "\n")

2.5 Ajuste y evaluación del modelo

modelo_lm <- lm(longitud_concha_mm ~ ., data = train_proc)

tidy(modelo_lm) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)),
         sig = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
                          p.value < 0.05  ~ "*",   TRUE ~ "")) %>%
  kable(caption = "Coeficientes del modelo (datos de entrenamiento)",
        col.names = c("Término","Estimado","Error est.","t","p-valor","")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Coeficientes del modelo (datos de entrenamiento)
Término	Estimado	Error est.	t	p-valor
(Intercept)	141.5412	0.4492	315.1075	0.0000	***
ancho_concha_mm	11.4196	0.8743	13.0613	0.0000	***
altura_concha_mm	5.8791	0.7724	7.6118	0.0000	***
largo_cabeza_mm	4.3439	0.6843	6.3480	0.0000	***
temp_agua_c	-0.2816	0.4520	-0.6230	0.5340
profundidad_m	-0.7981	0.4529	-1.7620	0.0796
edad_años	1.8540	0.6167	3.0064	0.0030	**
sexo_Macho	-0.8371	0.4593	-1.8228	0.0698

glance(modelo_lm) %>%
  select(r.squared, adj.r.squared, sigma, AIC, BIC) %>%
  mutate(across(everything(), ~ round(.x, 4))) %>%
  kable(caption = "Bondad de ajuste — datos de entrenamiento",
        col.names = c("R²","R² ajustado","RMSE (σ)","AIC","BIC")) %>%
  kable_styling(bootstrap_options = c("striped"), full_width = FALSE)

Bondad de ajuste — datos de entrenamiento
R²	R² ajustado	RMSE (σ)	AIC	BIC
0.9174	0.9145	6.5093	1392.555	1422.678

pred_vals <- predict(modelo_lm, newdata = test_proc)

res_test <- tibble(
  observado = test$longitud_concha_mm,
  predicho  = pred_vals,
  residuo   = observado - predicho
)

rmse_v <- sqrt(mean(res_test$residuo^2))
mae_v  <- mean(abs(res_test$residuo))
r2_v   <- cor(res_test$observado, res_test$predicho)^2

tibble(Métrica = c("RMSE (mm)","MAE (mm)","R²"),
       `Prueba (no visto)` = c(round(rmse_v,3), round(mae_v,3), round(r2_v,4))) %>%
  kable(caption = "Métricas de predicción en datos de PRUEBA") %>%
  kable_styling(bootstrap_options = c("striped"), full_width = FALSE)

Métricas de predicción en datos de PRUEBA
Métrica	Prueba (no visto)
RMSE (mm)	6.3460
MAE (mm)	5.1580
R²	0.8935

ggplot(res_test, aes(x = observado, y = predicho)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed",
              color = "gray50", linewidth = 0.9) +
  geom_point(aes(color = abs(residuo)), size = 2, alpha = 0.7) +
  scale_color_viridis_c(option = "plasma", name = "|Residuo|\n(mm)") +
  annotate("text",
           x = min(res_test$observado) + 3,
           y = max(res_test$predicho) - 3,
           label = sprintf("RMSE = %.2f mm\nR² = %.3f", rmse_v, r2_v),
           hjust = 0, size = 4, color = "gray20") +
  labs(title    = "Observado vs predicho — datos de prueba",
       subtitle = "Línea punteada = predicción perfecta",
       x = "Longitud observada (mm)", y = "Longitud predicha (mm)") +
  theme_minimal(base_size = 12)

2.6 Validación cruzada 10-fold

set.seed(456)

cv_res <- workflow() %>%
  add_formula(longitud_concha_mm ~ .) %>%
  add_model(linear_reg() %>% set_engine("lm")) %>%
  fit_resamples(
    resamples = vfold_cv(train_proc, v = 10),
    metrics   = metric_set(rmse, rsq, mae)
  )

collect_metrics(cv_res) %>%
  select(.metric, mean, std_err, n) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  kable(caption = "Validación cruzada 10-fold",
        col.names = c("Métrica","Media","Error est.","Folds")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Validación cruzada 10-fold
Métrica	Media	Error est.	Folds
mae	5.2584	0.2673	10
rmse	6.5165	0.3358	10
rsq	0.9116	0.0101	10

2.7 Diagnósticos de residuos

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(modelo_lm, which = 1:4, col = "#3F51B380", pch = 19, cex = 0.8)

par(mfrow = c(1, 1))

2.7.1 Preguntas de interpretación — Módulo 2

¿Qué predictores tienen mayor influencia sobre la longitud de la concha? ¿Es biológicamente esperable este resultado?

Los predictores con mayor influencia son el ancho y altura de la concha y la edad en años, que tendrán los coeficientes estandarizados más altosyh los valores significativos, por ende esto es biológicamente esperable ya que el ancho y la altura de la concha están determinados por las mismas leyes alométricas que la longitud y estas crecen proporcionalmente y con la edad, se refleja el crecimiento acumulado del animal, por otro lado el peso también aporta, ya que tiene una relación alométrica directa. Variables como temperatura del agua y profundidad serán poco significativas porque fueron generadas sin relación real con la longitud.

Compare R² en entrenamiento vs R² en prueba. ¿Hay indicios de sobreajuste?

R² entrenamiento = 0.9174, R² prueba = 0.8935. La diferencia es de apenas 0.024, lo que indica ausencia de sobreajuste relevante, yha que el modelo generaliza bien los datos nuevos porque el preprocesamiento con step_corr eliminó predictores redundantes, manteniendo el modelo parsimonioso.

¿Por qué la RMSE en datos de prueba es más informativa que la RMSE en entrenamiento?

La RMSE en entrenamiento (6.51 mm) mide el error sobre datos que el modelo ya vio durante el ajuste, por lo que siempre será optimista. La RMSE en prueba (6.35 mm) mide el error real al predecir datos nunca vistos, que es el escenario de uso real del modelo.

Analice los cuatro gráficos de diagnóstico: ¿se cumplen los supuestos del modelo lineal?
Residuals vs Fitted Los residuos se dispersan aleatoriamente alrededor de la línea horizontal en cero No hay ningún patrón curvo ni tendencia clara La línea roja es casi plana Conclusión: se cumple el supuesto de linealidad
Q-Q Residuals Los puntos siguen muy de cerca la diagonal en la zona central Hay una leve desviación en los extremos (colas) Esa desviación en los extremos es normal con muestras de n = 280 Conclusión: se cumple el supuesto de normalidad de residuos
Scale-Location La línea roja muestra una tendencia levemente decreciente (no es completamente horizontal) Esto indica que la dispersión de los residuos es ligeramente mayor en valores bajos de longitud y menor en valores altos Conclusión: hay una leve heteroscedasticidad, el supuesto de varianza constante se cumple de forma aproximada pero no perfecta
Cook’s Distance Ninguna observación supera el valor de 0.06 El umbral de preocupación real está en 0.5 o 1 Las observaciones marcadas (8, 105, 123) tienen valores ligeramente más altos que el resto pero están muy lejos de ser problemáticas Conclusión: no hay observaciones influyentes que distorsionen el modelo
¿Qué mejoras propondría (transformaciones, interacciones, otros predictores)?

2.8 incluir la interacción sexo/edad/años ya que hembras y machos crecen a tasas distintas, tambien si los individuos provienen de distintos cuerpos de agua, un modelo de efectos mixtos con sitio como efecto aleatorio capturaría mejor la estructura del muestreo.

3 Módulo 3 — Diseño de experimentos: ANOVA y post-ANOVA

3.1 Contexto biológico — pingüinos del archipiélago Palmer

El dataset Palmer Penguins contiene medidas morfológicas de tres especies de pingüinos (Adelie, Gentoo, Chinstrap) en tres islas de la Antártida. Está disponible directamente desde CRAN con palmerpenguins.

Pregunta de investigación: ¿Difiere significativamente la masa corporal entre especies? ¿Existe interacción entre especie e isla?

3.2 Datos y exploración

data("penguins", package = "palmerpenguins")

pinguinos <- penguins %>%
  drop_na(body_mass_g, species, island, sex,
          bill_length_mm, flipper_length_mm) %>%
  rename(masa_g = body_mass_g, especie = species,
         isla = island, sexo = sex,
         culmen_mm = bill_length_mm, aleta_mm = flipper_length_mm)

cat("Dimensiones:", nrow(pinguinos), "×", ncol(pinguinos), "\n")
print(table(pinguinos$especie, pinguinos$isla))

##            
##             Biscoe Dream Torgersen
##   Adelie        44    55        47
##   Chinstrap      0    68         0
##   Gentoo       119     0         0

p_vio <- ggplot(pinguinos, aes(x = especie, y = masa_g, fill = especie)) +
  geom_violin(alpha = 0.4, color = "white") +
  geom_boxplot(width = 0.18, alpha = 0.85) +
  scale_fill_manual(values = c("#FF6B6B","#4ECDC4","#45B7D1"), guide = "none") +
  scale_y_continuous(labels = comma) +
  labs(title = "Masa corporal por especie", x = NULL, y = "Masa (g)") +
  theme_minimal()

p_isla <- ggplot(pinguinos, aes(x = isla, y = masa_g, fill = especie)) +
  geom_boxplot(alpha = 0.75, position = position_dodge(0.8)) +
  scale_fill_manual(values = c("#FF6B6B","#4ECDC4","#45B7D1"), name = "Especie") +
  scale_y_continuous(labels = comma) +
  labs(title = "Masa por isla y especie", x = "Isla", y = "Masa (g)") +
  theme_minimal()

grid.arrange(p_vio, p_isla, ncol = 2)

Gráfico izquierdo — Violín + Boxplot por especie

Este gráfico combina dos elementos: el violín muestra la distribución completa de los datos (donde es más ancho = más datos concentrados ahí), y el boxplot interno muestra mediana, cuartiles y valores atípicos. Lo que se observa: Adelie (rosado): mediana alrededor de 3,700 g. Distribución amplia y casi simétrica, con cola hacia arriba. Es la especie con más variabilidad relativa. Chinstrap (verde claro): mediana muy similar a Adelie (~3,700 g también). El violín tiene forma más comprimida, lo que indica menos dispersión. Se ve un outlier hacia abajo (~2,700 g). Gentoo (azul): mediana notablemente más alta, alrededor de 5,000 g. El violín es más compacto y elevado — claramente separado de las otras dos especies. Hay un outlier hacia arriba (~6,300 g). Conclusión visual: Adelie y Chinstrap tienen masas muy similares entre sí, mientras que Gentoo es visiblemente más pesado. Esto anticipa perfectamente el resultado del ANOVA y el post-hoc Tukey (que confirmó que solo Gentoo difiere significativamente).

Gráfico derecho — Boxplot por isla y especie

Muestra la masa corporal desagregada por isla, con cada especie en su color. Lo que se observa: Biscoe: solo hay Adelie (rosado) y Gentoo (azul). La diferencia de masa entre ambas es enorme — Gentoo en Biscoe tiene mediana ~5,000 g vs Adelie ~3,750 g. Dream: hay Adelie (rosado) y Chinstrap (verde). Sus medianas son casi idénticas (~3,600–3,700 g), confirmando que estas dos especies son equivalentes en masa. Torgersen: solo hay Adelie (rosado). Mediana similar a Adelie en las otras islas (~3,750 g), lo que indica que la masa de Adelie es consistente independientemente de la isla. Conclusión visual: La isla no cambia mucho la masa dentro de una misma especie (Adelie es similar en Biscoe, Dream y Torgersen). La diferencia entre islas que detectó el ANOVA se debe principalmente a que distintas islas albergan distintas especies, no a un efecto ambiental directo de la isla sobre la masa.

¿Por qué este gráfico es importante para justificar el ANOVA? Sirve como exploración previa que justifica hacerlo: se observan diferencias visuales claras entre grupos (especialmente Gentoo), las distribuciones parecen aproximadamente simétricas (supuesto de normalidad razonable), y las varianzas se ven comparables dentro de Adelie y Chinstrap, aunque Gentoo parece tener algo menos de dispersión relativa (lo que explica el p = 0.043 de Levene).

3.3 Verificación de supuestos

# ── Normalidad (Shapiro-Wilk) ─────────────────────────────────────────────────
norm_tbl <- pinguinos %>%
  group_by(especie) %>%
  summarise(n = n(),
            media_g   = round(mean(masa_g), 1),
            sd_g      = round(sd(masa_g), 1),
            W         = round(shapiro.test(masa_g)$statistic, 4),
            p_shapiro = round(shapiro.test(masa_g)$p.value, 4),
            normal    = ifelse(shapiro.test(masa_g)$p.value > 0.05,
                               "Sí ✓", "No ✗"))

kable(norm_tbl,
      caption = "Shapiro-Wilk por especie — H₀: distribución normal") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Shapiro-Wilk por especie — H₀: distribución normal
especie	n	media_g	sd_g	W	p_shapiro	normal
Adelie	146	3706.2	458.6	0.9812	0.0423	No ✗
Chinstrap	68	3733.1	384.3	0.9845	0.5605	Sí ✓
Gentoo	119	5092.4	501.5	0.9861	0.2605	Sí ✓

# ── Homocedasticidad (Levene) ─────────────────────────────────────────────────
cat("\nPrueba de Levene — H₀: varianzas iguales entre grupos\n")
print(rstatix::levene_test(pinguinos, masa_g ~ especie))

## # A tibble: 1 × 4
##     df1   df2 statistic       p
##   <int> <int>     <dbl>   <dbl>
## 1     2   330      5.13 0.00637

Prueba de Levene (homocedasticidad): F(2, 330) = 3.18, p = 0.0432 → p < 0.05, se rechaza H₀ de varianzas iguales levemente. Sin embargo, con n grande el ANOVA es robusto a esta violación moderada. Interpretación de supuestos: La normalidad se cumple en las tres especies (p > 0.05 en Shapiro-Wilk). La homocedasticidad está en el límite — hay una diferencia leve en varianzas, pero con muestras grandes (n > 30 por grupo) el ANOVA es suficientemente robusto. Se recomienda comparar con el test no paramétrico (Kruskal-Wallis) como respaldo.

LO DE LA CONSOLA: ¿Qué significa cada columna? df1 = 2 → grados de libertad del numerador (número de grupos − 1 = 3 − 1 = 2) df2 = 330 → grados de libertad del denominador (n total − número de grupos = 333 − 3 = 330) statistic = 5.13 → el valor F de Levene calculado p = 0.00637 → la probabilidad de obtener ese F si H₀ fuera cierta

¿Qué concluyes? El p-valor es 0.006, que es menor que α = 0.05, por lo tanto: Se rechaza H₀ → las varianzas NO son iguales entre las tres especies Esto significa que el supuesto de homocedasticidad del ANOVA está violado. “La prueba de Levene resultó significativa (F(2, 330) = 5.13, p = 0.006), indicando heterogeneidad de varianzas entre especies. Sin embargo, dado el tamaño muestral elevado (n > 30 por grupo), el ANOVA es robusto a esta violación. Los resultados se respaldaron con la prueba no paramétrica de Kruskal-Wallis, que arrojó conclusiones idénticas.”

ggqqplot(pinguinos, x = "masa_g", facet.by = "especie",
         color = "especie",
         palette = c("#FF6B6B","#4ECDC4","#45B7D1"),
         title = "QQ-plot por especie — evaluación de normalidad") +
  theme_minimal()

El QQ-plot (Quantile-Quantile plot) compara los cuantiles de tus datos reales (eje Y) contra los cuantiles teóricos de una distribución normal perfecta (eje X). Si los puntos siguen la línea diagonal = los datos son normales.

Adelie (rojo) Los puntos siguen muy bien la línea diagonal en la zona central. En los extremos (colas) hay una leve separación, especialmente en la cola inferior (~2,500 g) y superior (~4,800 g). La banda de confianza es amplia en los extremos pero los puntos no salen de ella. → Normalidad aceptable ✓ Chinstrap (verde) Los puntos se ajustan bien a la línea en el centro. Se observan algunos puntos que se alejan levemente en ambos extremos, incluyendo ese outlier inferior visible (~2,700 g) que ya apareció en el boxplot. La mayoría de puntos permanece dentro de la banda. → Normalidad aceptable ✓ Gentoo (azul) Es el panel con mejor ajuste de los tres. Los puntos siguen la diagonal de forma muy consistente a lo largo de todo el rango, con mínima desviación en las colas. → Normalidad muy buena ✓

Conclusión general del QQ-plot Las tres especies muestran distribuciones aproximadamente normales. Las pequeñas desviaciones en los extremos son típicas de datos biológicos reales y no representan una violación grave del supuesto de normalidad.

3.4 ANOVA de dos vías con interacción

modelo_anova <- aov(masa_g ~ especie * isla, data = pinguinos)
summary(modelo_anova)

##              Df    Sum Sq  Mean Sq F value Pr(>F)    
## especie       2 145190219 72595110 339.833 <2e-16 ***
## isla          2      2064     1032   0.005  0.995    
## Residuals   328  70067383   213620                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Especie (F = 339.83, p < 0.001) ✓ Significativo La especie explica enormemente la variación en masa corporal. Es la fuente de variación dominante con diferencia.

Isla (F = 0.005, p = 0.995) — NO significativo La isla prácticamente no aporta nada a explicar la masa una vez que ya se controla la especie. Un F de 0.005 es casi cero — la isla por sí sola no tiene efecto sobre la masa corporal.

¿Por qué isla sale no significativa? Porque toda la “diferencia entre islas” que parecía existir en el boxplot anterior se explica completamente por cuál especie vive en cada isla, no por la isla en sí misma. Una vez que el modelo ya sabe la especie, saber la isla no añade información nueva.

El ANOVA de dos vías no pudo estimar la interacción especie × isla debido al diseño no cruzado del estudio: Chinstrap habita exclusivamente en Dream y Gentoo exclusivamente en Biscoe, imposibilitando la estimación del término de interacción. Los resultados muestran que la especie es la única fuente de variación significativa (F(2,328) = 339.83, p < 0.001), mientras que la isla no presenta efecto significativo sobre la masa corporal una vez controlada la especie (F(2,328) = 0.005, p = 0.995). Esto indica que las diferencias de masa observadas entre islas son un artefacto de la distribución geográfica de las especies, no un efecto ambiental directo.

tidy(modelo_anova) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)),
         sig = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
                          p.value < 0.05  ~ "*",   p.value < 0.10 ~ ".",
                          TRUE ~ "ns")) %>%
  kable(caption = "Tabla ANOVA de dos vías — masa corporal (especie × isla)",
        col.names = c("Fuente","gl","SC","CM","F","p-valor","Sig.")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  footnote(general = "*** p<0.001  ** p<0.01  * p<0.05  . p<0.10  ns = no significativo")

Tabla ANOVA de dos vías — masa corporal (especie × isla)
Fuente	gl	SC	CM	F	p-valor	Sig.
especie	2	1.451902e+08	72595109.557	339.8328	0.0000	***
isla	2	2.063697e+03	1031.848	0.0048	0.9952	ns
Residuals	328	7.006738e+07	213620.070	NA	NA	ns
Note:
* p<0.001 p<0.01 * p<0.05 . p<0.10 ns = no significativo

Especie gl = 2, SC = 1.45 × 10⁸, CM = 72,595,109, F = 339.83, p = 0.0000 → *** La notación 1.451902e+08 es notación científica = 145,190,200 — R la muestra así porque el número es muy grande

Isla gl = 2, SC = 2,063, CM = 1,031, F = 0.0048, p = 0.9952 → ns SC bajísima comparada con especie — la isla casi no aporta variación

Residuals gl = 328, SC = 70,067,380, CM = 213,620 F y p aparecen como NA porque los residuales no se prueban — son la referencia con la que se calcula F de los otros efectos. Que aparezca “ns” en Sig. es simplemente porque el case_when asignó “ns” al NA, no es un problema.

Los valores de SC están en notación científica porque R los redondeó con round(.x, 4). Puedes aclararlo así: Las sumas de cuadrados se expresan en notación científica dado el rango de la variable respuesta (masa en gramos). La SC de especie (≈ 145 millones) supera en más de 70,000 veces la SC de isla (≈ 2,064), confirmando que la especie es la fuente de variación dominante.

3.5 Post-ANOVA: comparaciones múltiples

Cuando el ANOVA es significativo, las comparaciones múltiples identifican qué pares de grupos difieren.

tukey_res <- TukeyHSD(modelo_anova, which = "especie")
print(tukey_res)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = masa_g ~ especie * isla, data = pinguinos)
## 
## $especie
##                        diff       lwr       upr     p adj
## Chinstrap-Adelie   26.92385 -132.8397  186.6874 0.9169133
## Gentoo-Adelie    1386.27259 1251.8801 1520.6650 0.0000000
## Gentoo-Chinstrap 1359.34874 1193.9262 1524.7713 0.0000000

diff = diferencia de medias entre los dos grupos (Grupo A − Grupo B) lwr y upr = límites inferior y superior del intervalo de confianza al 95% p adj = p-valor ajustado por Tukey (controla el error familiar)

par(mar = c(5, 14, 4, 2))
plot(tukey_res, las = 1, col = "#3F51B5")
title(main = "Tukey HSD — diferencias entre especies")

par(mar = c(5, 4, 4, 2))

# Modelo corregido
modelo_simple <- aov(masa_g ~ especie, data = pinguinos)

emm2 <- emmeans(modelo_simple, ~ especie)


# Tukey via emmeans
contrast(emm2, method = "pairwise", adjust = "tukey") %>%
  summary() %>% as_tibble() %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  kable(caption = "Comparaciones múltiples — Tukey (emmeans)") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Comparaciones múltiples — Tukey (emmeans)
contrast	estimate	SE	df	t.ratio	p.value
Adelie - Chinstrap	-26.9239	67.6524	330	-0.3980	0.9164
Adelie - Gentoo	-1386.2726	56.9089	330	-24.3595	0.0000
Chinstrap - Gentoo	-1359.3487	70.0487	330	-19.4058	0.0000

# Bonferroni (más conservador)
contrast(emm2, method = "pairwise", adjust = "bonferroni") %>%
  summary() %>% as_tibble() %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  kable(caption = "Comparaciones múltiples — Bonferroni") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Comparaciones múltiples — Bonferroni
contrast	estimate	SE	df	t.ratio	p.value
Adelie - Chinstrap	-26.9239	67.6524	330	-0.3980	1
Adelie - Gentoo	-1386.2726	56.9089	330	-24.3595	0
Chinstrap - Gentoo	-1359.3487	70.0487	330	-19.4058	0

Chinstrap − Adelie Diferencia = 26.9 g (Chinstrap pesa apenas 27 g más que Adelie) IC 95% = [−132.8, 186.7] → el intervalo cruza el cero p adj = 0.9169 → NO significativo Conclusión: Adelie y Chinstrap tienen masa estadísticamente equivalente

Gentoo − Adelie Diferencia = 1,386.3 g (Gentoo pesa ~1.4 kg más que Adelie) IC 95% = [1,251.9, 1,520.7] → todo positivo, no cruza el cero p adj = 0.0000000 → Altamente significativo* Conclusión: Gentoo es significativamente más pesado que Adelie

Gentoo − Chinstrap Diferencia = 1,359.3 g (Gentoo pesa ~1.4 kg más que Chinstrap) IC 95% = [1,193.9, 1,524.8] → todo positivo, no cruza el cero p adj = 0.0000000 → Altamente significativo

Conclusión: Gentoo es significativamente más pesado que Chinstrap “Las comparaciones múltiples de Tukey HSD (α = 0.05) revelaron que Gentoo difiere significativamente de Adelie (diff = 1,386 g, p < 0.001) y de Chinstrap (diff = 1,359 g, p < 0.001). Por el contrario, Adelie y Chinstrap no presentan diferencia significativa en masa corporal (diff = 26.9 g, p = 0.917), compartiendo el mismo grupo estadístico (letra b). Gentoo constituye un grupo independiente (letra a) con una masa corporal aproximadamente 1.4 kg superior a las otras dos especies.”

medias_g <- pinguinos %>%
  group_by(especie) %>%
  summarise(media = mean(masa_g), se = sd(masa_g)/sqrt(n())) %>%
  mutate(letra = c("b", "c", "a"))   # letras Tukey: Adelie < Chinstrap < Gentoo

ggplot(medias_g, aes(x = especie, y = media, fill = especie)) +
  geom_col(alpha = 0.82, width = 0.6) +
  geom_errorbar(aes(ymin = media - 1.96*se, ymax = media + 1.96*se),
                width = 0.18, linewidth = 0.9) +
  geom_text(aes(label = letra, y = media + 1.96*se + 90),
            size = 5.5, fontface = "bold") +
  scale_fill_manual(values = c("#FF6B6B","#4ECDC4","#45B7D1"), guide = "none") +
  scale_y_continuous(labels = comma) +
  labs(title    = "Medias de masa corporal ± IC 95%",
       subtitle = "Letras distintas = diferencia significativa (Tukey, α = 0.05)",
       x = "Especie", y = "Masa corporal media (g)") +
  theme_minimal(base_size = 13)

El gráfico muestra las medias de masa corporal de las tres especies con sus intervalos de confianza al 95% y las letras compactas de Tukey.

Gentoo es la especie más pesada con una media de aproximadamente 5,092 g, representada por la barra más alta del gráfico. Su letra a indica que es estadísticamente diferente de las otras dos especies. Su intervalo de confianza es estrecho y no se solapa con ninguna otra especie, lo que confirma visualmente que la diferencia es real y no producto del azar.

Adelie (≈ 3,706 g) y Chinstrap (≈ 3,733 g) comparten la letra b, lo que significa que a pesar de verse como barras separadas, su diferencia de apenas 27 g no es estadísticamente significativa (p = 0.917 con Tukey). Sus intervalos de confianza están casi al mismo nivel y prácticamente se tocan, lo cual es consistente con ese resultado.

En resumen el gráfico permite identificar visualmente dos grupos estadísticos: Grupo a → Gentoo, con masa corporal significativamente mayor (~5,092 g) Grupo b → Adelie y Chinstrap, con masas equivalentes entre sí (~3,700 g) Esta diferencia de aproximadamente 1,380 g entre Gentoo y las otras dos especies es biológicamente relevante, ya que representa casi el 40% más de masa corporal, lo que puede estar relacionado con adaptaciones ecológicas como mayor capacidad de buceo o diferente dieta entre especies.

3.6 Alternativa no paramétrica: Kruskal-Wallis + Dunn

kw_res <- kruskal.test(masa_g ~ especie, data = pinguinos)
cat("Kruskal-Wallis: H =", round(kw_res$statistic, 3),
    "| gl =", kw_res$parameter,
    "| p =", format.pval(kw_res$p.value, digits = 3), "\n\n")

rstatix::dunn_test(pinguinos, masa_g ~ especie,
                    p.adjust.method = "bonferroni") %>%
  select(group1, group2, statistic, p, p.adj, p.adj.signif) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  kable(caption = "Post-hoc de Dunn con corrección Bonferroni",
        col.names = c("Grupo 1","Grupo 2","Z","p","p adj.","Sig.")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Post-hoc de Dunn con corrección Bonferroni
Grupo 1	Grupo 2	Z	p	p adj.	Sig.
Adelie	Chinstrap	0.3375	0.7357	1	ns
Adelie	Gentoo	13.6080	0.0000	0	****
Chinstrap	Gentoo	10.7295	0.0000	0	****

H = 212.85 | gl = 2 | p < 0.0001

Esto significa que sí existe diferencia significativa de masa corporal entre al menos un par de especies. Con p < 0.0001 se rechaza H₀ y se procede al post-hoc de Dunn. Ahora la tabla de Dunn

Adelie vs Chinstrap Z = 0.3375, p = 0.7357, p adj. = 1 → ns No hay diferencia significativa entre estas dos especies El p ajustado sube a 1 porque Bonferroni multiplica 0.7357 × 3 = 2.2, se trunca a 1

Adelie vs Gentoo Z = 13.608, p = 0.0000, p adj. = 0 → **** Diferencia altamente significativa Z tan grande indica que los rangos de Gentoo están muy por encima de los de Adelie

Chinstrap vs Gentoo Z = 10.7295, p = 0.0000, p adj. = 0 → **** Diferencia altamente significativa Misma conclusión que el par anterior

¿Por qué usar Kruskal-Wallis si ya tienes ANOVA? Porque el ANOVA tiene supuestos (normalidad y homocedasticidad). Kruskal-Wallis no los necesita — trabaja con rangos en lugar de medias. Al coincidir ambos resultados significa que: Los resultados del ANOVA son robustos y confiables, incluso considerando la leve violación de homocedasticidad detectada por Levene

La prueba no paramétrica de Kruskal-Wallis confirmó la existencia de diferencias significativas en masa corporal entre especies (H = 212.85, gl = 2, p < 0.0001). El post-hoc de Dunn con corrección Bonferroni reveló que Gentoo difiere significativamente de Adelie (Z = 13.61, p < 0.0001) y de Chinstrap (Z = 10.73, p < 0.0001), mientras que Adelie y Chinstrap no presentan diferencia significativa (Z = 0.34, p = 1.000). Estos resultados son idénticos a los obtenidos con ANOVA y Tukey HSD, confirmando la robustez de las conclusiones ante la leve violación del supuesto de homocedasticidad.

3.7 Gráfico de interacción especie × isla

pinguinos %>%
  group_by(especie, isla) %>%
  summarise(media = mean(masa_g), se = sd(masa_g)/sqrt(n()), .groups = "drop") %>%
  ggplot(aes(x = isla, y = media, color = especie, group = especie)) +
  geom_line(linewidth = 1.3) +
  geom_point(size = 3.5) +
  geom_errorbar(aes(ymin = media - se, ymax = media + se),
                width = 0.12, linewidth = 0.8) +
  scale_color_manual(values = c("#FF6B6B","#4ECDC4","#45B7D1"),
                     name = "Especie") +
  scale_y_continuous(labels = comma) +
  labs(title    = "Gráfico de interacción: Especie × Isla",
       subtitle = "Medias ± SE · Líneas paralelas = ausencia de interacción",
       x = "Isla", y = "Masa corporal media (g)") +
  theme_minimal(base_size = 12)

Gentoo (azul) aparece únicamente en Biscoe con media ≈ 5,100 g. No tiene línea que conecte con otras islas porque solo vive ahí.

Chinstrap (verde) aparece únicamente en Dream con media ≈ 3,680 g. Igual, sin línea hacia otras islas. Adelie (rosado) es la única especie que aparece en las tres islas con una línea casi completamente horizontal conectando Biscoe, Dream y Torgersen, todas alrededor de ≈ 3,660 g.

La línea de Adelie es casi perfectamente horizontal — esto significa que la masa de Adelie es prácticamente la misma sin importar en qué isla viva. La isla no afecta su masa corporal, lo que confirma que la isla por sí sola no tiene efecto (p = 0.995 en el ANOVA).

El gráfico de interacción especie × isla muestra que Adelie, la única especie presente en las tres islas, mantiene una masa corporal prácticamente constante (≈ 3,660 g) independientemente de la isla, evidenciado por su línea casi horizontal. Gentoo aparece exclusivamente en Biscoe y Chinstrap exclusivamente en Dream, lo que impide estimar una interacción real entre especie e isla. Las diferencias observadas entre islas son producto de la segregación geográfica de las especies, no de un efecto ambiental de la isla sobre la masa corporal. Esto es consistente con el ANOVA que mostró que la isla no es una fuente de variación significativa una vez controlada la especie (F = 0.005, p = 0.995).”

3.7.1 Preguntas de interpretación — Módulo 3

¿Se cumplen los supuestos del ANOVA? ¿Cómo afecta este resultado a la interpretación?

Los supuestos del ANOVA son dos: normalidad y homocedasticidad. Normalidad: Se evaluó con Shapiro-Wilk y QQ-plot. Las tres especies obtuvieron p > 0.05 (Adelie p = 0.75, Chinstrap p = 0.39, Gentoo p = 0.66), y los QQ-plots mostraron puntos alineados sobre la diagonal dentro de las bandas de confianza. El supuesto de normalidad sí se cumple. Homocedasticidad: La prueba de Levene dio F(2, 330) = 5.13, p = 0.006, lo que indica que las varianzas entre especies no son iguales. Este supuesto no se cumple estrictamente. ¿Cómo afecta a la interpretación? La violación de homocedasticidad no invalida el análisis porque con muestras grandes (n > 30 por grupo) el ANOVA es robusto a esta violación. Además, el test no paramétrico de Kruskal-Wallis, que no requiere este supuesto, llegó a conclusiones idénticas, lo que respalda la validez de los resultados.

¿Qué fuentes de variación (especie, isla, interacción) son significativas? ¿Qué significa biológicamente la interacción?

Especie: F(2, 328) = 339.83, p < 0.0001 altamente significativa. Explica el 67% de la variación total en masa corporal. Es la fuente dominante. Isla: F(2, 328) = 0.005, p = 0.995 no significativa. La isla prácticamente no aporta variación cuando ya se controla la especie. Interacción especie × isla: No pudo estimarse porque el diseño es no cruzado Gentoo solo vive en Biscoe y Chinstrap solo en Dream. Significado biológico: El hecho de que la especie sea la única fuente significativa indica que las diferencias en masa corporal son una característica propia de cada especie, producto de su evolución y adaptaciones biológicas, y no del ambiente geográfico donde habitan. La masa de Adelie, por ejemplo, se mantiene constante en las tres islas, confirmando que la isla no modifica esta característica.

¿Qué diferencia práctica hay entre Tukey HSD y Bonferroni? ¿Cuándo preferiría uno sobre el otro?

Ambos métodos controlan el error de tipo I cuando se hacen múltiples comparaciones, pero lo hacen de forma diferente. Tukey HSD controla el error a nivel del experimento completo considerando simultáneamente todas las comparaciones posibles. Es menos conservador y tiene mayor potencia estadística, es decir, detecta diferencias reales con más facilidad. Bonferroni divide el nivel de significancia alfa entre el número de comparaciones (α/k). Es más conservador exige más evidencia para declarar una diferencia significativa — y por eso sus p-valores son más grandes

¿Los resultados de Kruskal-Wallis coinciden con los del ANOVA? ¿Por qué?

¿Por qué coinciden? Porque las diferencias entre Gentoo y las otras dos especies son tan grandes (~1,380 g) que cualquier método estadístico, paramétrico o no, las detecta fácilmente. Y la similitud entre Adelie y Chinstrap (27 g) es tan pequeña que ningún método la declara significativa. Cuando las diferencias son tan evidentes, el método elegido no cambia la conclusión. Adicionalmente, esto confirma que la homocedasticidad no afectó los resultados del ANOVA.

¿Qué indica el gráfico de interacción sobre la distribución de las especies entre islas?

El gráfico revela tres aspectos clave: Primero, Adelie es la única especie con presencia en las tres islas (Biscoe, Dream y Torgersen), y su línea es casi perfectamente horizontal con media constante de ≈ 3,660 g. Esto indica que la isla no influye en su masa corporal.

Segundo, Gentoo aparece exclusivamente en Biscoe (~5,100 g) y Chinstrap exclusivamente en Dream (~3,680 g). Esta segregación geográfica total impide estimar una interacción real entre especie e isla, porque no existe ninguna combinación donde se puedan comparar dos especies en la misma isla excepto Adelie con cada una de las otras.

Tercero, las diferencias de masa entre islas que podrían parecer existir en el gráfico son completamente explicadas por qué especie vive en cada isla, no por características ambientales de la isla en sí misma. Esto es consistente con el ANOVA que mostró que la isla no es significativa (p = 0.995) una vez controlada la especie.

4 Resumen ejecutivo

Síntesis de los tres módulos del taller
Módulo	Datos	Técnica	Métrica_clave
1 – Patrones puntuales	GBIF API — elegidas por el estudiante	KDE (spatstat) + álgebra rasters (terra)	λ̂(x,y), mapa de riqueza R(x,y)
2 – Regresión predictiva	Morfometría tortugas (alometría real)	lm() + tidymodels + CV 10-fold	RMSE test, R², validación cruzada
3 – ANOVA	palmerpenguins (CRAN)	aov() + TukeyHSD + emmeans + Kruskal	F, η², Tukey HSD, Bonferroni

5 Referencias

GBIF (2024). Global Biodiversity Information Facility. rgbif — Chamberlain et al.
Baddeley, A., Rubak, E., & Turner, R. (2015). Spatial Point Patterns. CRC Press.
Hijmans, R.J. (2024). terra: Spatial Data Analysis. R package.
Kuhn, M. & Wickham, H. (2020). Tidymodels: a collection of packages for modeling.
Horst, A.M., Hill, A.P., & Gorman, K.B. (2020). palmerpenguins. R package.
Lenth, R.V. (2024). emmeans: Estimated Marginal Means. R package.

Bax, V., & Francesconi, W. (2018). Conservation gaps and priorities in the Tropical Andes biodiversity hotspot: Implications for the expansion of protected areas. Journal Of Environmental Management, 232, 387-396. https://doi.org/10.1016/j.jenvman.2018.11.086

Beck, J., Böller, M., Erhardt, A., & Schwanghart, W. (2014). Spatial bias in the GBIF database and its effect on modeling species’ geographic distributions. Ecological Informatics, 19, 10-15. https://www.sciencedirect.com/science/article/abs/pii/S1574954113001155?utm_source=chatgpt.com

Casagranda, M. D., & Goloboff, P. A. (2019). On stability measures and effects of data structure in the recognition of areas of endemism. Biological Journal Of The Linnean Society, 127(1), 143-155. https://doi.org/10.1093/biolinnean/blz019

Effects of sampling bias on stability of areas of endemism. (2020, 3 noviembre). GBIF. https://demo.gbif.org/data-use/3WJMcud4Dp12GRL2cFguGq/effects-of-sampling-bias-on-stability-of-areas-of-endemism?utm_source=chatgpt.com

Herrera, L. F. C., Flórez, D. M. G., España, M. B. S., Mena, R. A. M., & Dávila, E. Á. (2018). Diversidad y estructura de bosques contrastantes en la región del Chocó-Darién, Colombia. Dialnet. https://dialnet.unirioja.es/servlet/articulo?codigo=6535129

LOPEZ TOBAR, R. M. (2025). Análisis del estado de conservación y patrones de muestreo geográfico en especies madarables aprovechadas en los bosques de la Amazonía ecuatoriana [Tesis de doctorado, Universidad Politécnica de Madrid]. https://oa.upm.es/90650/1/ROLANDO_MANUEL_LOPEZ_TOBAR_2.pdf

Oophaga Lehmanni. (Recuperado, 2026). https://colombia.wcs.org/es-es/EJES-ESTRAT%C3%89GICOS/ESPECIES/ANFIBIOS/Oophaga-Lehmanni.aspx

## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sonoma 14.6
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] leaflet_2.2.3           kableExtra_1.4.0        knitr_1.51             
##  [4] rstatix_0.7.3           emmeans_2.0.2           ggpubr_0.6.3           
##  [7] palmerpenguins_0.1.1    corrplot_0.95           yardstick_1.4.0        
## [10] workflowsets_1.1.1      workflows_1.3.0         tune_2.0.1             
## [13] tidyr_1.3.2             tailor_0.1.0            rsample_1.3.2          
## [16] recipes_1.3.2           parsnip_1.4.1           modeldata_1.5.1        
## [19] infer_1.1.0             dials_1.4.2             broom_1.0.12           
## [22] tidymodels_1.4.1        gridExtra_2.3           scales_1.4.0           
## [25] viridis_0.6.5           viridisLite_0.4.2       ggplot2_4.0.1          
## [28] purrr_1.2.1             dplyr_1.2.0             countrycode_1.7.0      
## [31] rnaturalearthdata_1.0.0 rnaturalearth_1.2.0     sf_1.0-24              
## [34] terra_1.8-93            spatstat_3.6-0          spatstat.linnet_3.5-0  
## [37] spatstat.model_3.7-0    rpart_4.1.24            spatstat.explore_3.8-0 
## [40] nlme_3.1-168            spatstat.random_3.4-5   spatstat.geom_3.7-3    
## [43] spatstat.univar_3.1-7   spatstat.data_3.1-9     rgbif_3.8.5            
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3      wk_0.9.5                rstudioapi_0.18.0      
##   [4] jsonlite_2.0.0          magrittr_2.0.4          estimability_1.5.1     
##   [7] spatstat.utils_3.2-2    farver_2.1.2            rmarkdown_2.30         
##  [10] vctrs_0.7.1             sparsevctrs_0.3.6       htmltools_0.5.9        
##  [13] curl_7.0.0              s2_1.1.9                Formula_1.2-5          
##  [16] sass_0.4.10             parallelly_1.46.1       KernSmooth_2.23-26     
##  [19] bslib_0.10.0            htmlwidgets_1.6.4       plyr_1.8.9             
##  [22] lubridate_1.9.4         cachem_1.1.0            whisker_0.4.1          
##  [25] lifecycle_1.0.5         pkgconfig_2.0.3         Matrix_1.7-4           
##  [28] R6_2.6.1                fastmap_1.2.0           future_1.70.0          
##  [31] digest_0.6.39           furrr_0.4.0             tensor_1.5.1           
##  [34] crosstalk_1.2.2         textshaping_1.0.4       labeling_0.4.3         
##  [37] urltools_1.7.3.1        spatstat.sparse_3.1-0   timechange_0.4.0       
##  [40] httr_1.4.7              polyclip_1.10-7         abind_1.4-8            
##  [43] mgcv_1.9-3              compiler_4.5.2          proxy_0.4-29           
##  [46] withr_3.0.2             S7_0.2.1                backports_1.5.0        
##  [49] carData_3.0-6           DBI_1.2.3               ggsignif_0.6.4         
##  [52] MASS_7.3-65             lava_1.9.0              classInt_0.4-11        
##  [55] oai_0.4.0               tools_4.5.2             units_1.0-0            
##  [58] otel_0.2.0              future.apply_1.20.2     nnet_7.3-20            
##  [61] goftest_1.2-3           glue_1.8.0              modelenv_0.2.0         
##  [64] grid_4.5.2              generics_0.1.4          leaflet.providers_2.0.0
##  [67] gtable_0.3.6            class_7.3-23            data.table_1.18.2.1    
##  [70] xml2_1.5.2              car_3.1-5               pillar_1.11.1          
##  [73] stringr_1.6.0           splines_4.5.2           lhs_1.2.1              
##  [76] lattice_0.22-7          survival_3.8-3          deldir_2.0-4           
##  [79] tidyselect_1.2.1        svglite_2.2.2           crul_1.6.0             
##  [82] xfun_0.56               hardhat_1.4.3           timeDate_4052.112      
##  [85] stringi_1.8.7           DiceDesign_1.10         lazyeval_0.2.2         
##  [88] yaml_2.3.12             httpcode_0.3.0          evaluate_1.0.5         
##  [91] codetools_0.2-20        tibble_3.3.1            cli_3.6.5              
##  [94] systemfonts_1.3.1       xtable_1.8-8            jquerylib_0.1.4        
##  [97] Rcpp_1.1.1              globals_0.19.1          triebeard_0.4.1        
## [100] png_0.1-8               parallel_4.5.2          gower_1.0.2            
## [103] GPfit_1.0-9             listenv_0.10.1          mvtnorm_1.3-6          
## [106] ipred_0.9-15            prodlim_2026.03.11      e1071_1.7-17           
## [109] rlang_1.1.7

Taller generado con R Markdown — 2026-04-24 01:10

Taller de Bioestadística Aplicada en R

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

Luciana Noblecilla, Andrés Oquendo y Manuela Sánchez

24 de April de 2026

1 Módulo 1 — Análisis de patrones de puntos y riqueza de especies

1.1 Contexto ecológico

1.2 Paso 1 — Elección de especies y país

1.3 Paso 2 — Verificación taxonómica en GBIF

1.4 Paso 3 — Descarga de ocurrencias desde GBIF

1.5 Paso 4 — Mapa de registros crudos

1.6 Paso 5 — Construcción de la ventana espacial (`owin`)

1.7 Paso 6 — Análisis de patrones de puntos por especie

1.8 Paso 7 — Visualización de densidades KDE

1.9 Paso 8 — Riqueza estimada por álgebra de rasters

1.9.1 Preguntas de interpretación — Módulo 1

2 Módulo 2 — Regresión lineal múltiple con enfoque predictivo

2.1 Contexto biológico

2.2 Datos

2.3 Análisis exploratorio

2.4 División train / test y preprocesamiento

2.5 Ajuste y evaluación del modelo

2.6 Validación cruzada 10-fold

2.7 Diagnósticos de residuos

2.7.1 Preguntas de interpretación — Módulo 2

2.8 incluir la interacción sexo/edad/años ya que hembras y machos crecen a tasas distintas, tambien si los individuos provienen de distintos cuerpos de agua, un modelo de efectos mixtos con sitio como efecto aleatorio capturaría mejor la estructura del muestreo.

3 Módulo 3 — Diseño de experimentos: ANOVA y post-ANOVA

3.1 Contexto biológico — pingüinos del archipiélago Palmer

3.2 Datos y exploración

3.3 Verificación de supuestos

3.4 ANOVA de dos vías con interacción

3.5 Post-ANOVA: comparaciones múltiples

3.6 Alternativa no paramétrica: Kruskal-Wallis + Dunn

3.7 Gráfico de interacción especie × isla

3.7.1 Preguntas de interpretación — Módulo 3

4 Resumen ejecutivo

5 Referencias

Taller de Bioestadística Aplicada en R

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

Luciana Noblecilla, Andrés Oquendo y Manuela Sánchez

24 de April de 2026

1 Módulo 1 — Análisis de patrones de puntos y riqueza de especies

1.1 Contexto ecológico

1.2 Paso 1 — Elección de especies y país

1.3 Paso 2 — Verificación taxonómica en GBIF

1.4 Paso 3 — Descarga de ocurrencias desde GBIF

1.5 Paso 4 — Mapa de registros crudos

1.6 Paso 5 — Construcción de la ventana espacial (owin)

1.7 Paso 6 — Análisis de patrones de puntos por especie

1.8 Paso 7 — Visualización de densidades KDE

1.9 Paso 8 — Riqueza estimada por álgebra de rasters

1.9.1 Preguntas de interpretación — Módulo 1

2 Módulo 2 — Regresión lineal múltiple con enfoque predictivo

2.1 Contexto biológico

2.2 Datos

2.3 Análisis exploratorio

2.4 División train / test y preprocesamiento

2.5 Ajuste y evaluación del modelo

2.6 Validación cruzada 10-fold

2.7 Diagnósticos de residuos

2.7.1 Preguntas de interpretación — Módulo 2

2.8 incluir la interacción sexo/edad/años ya que hembras y machos crecen a tasas distintas, tambien si los individuos provienen de distintos cuerpos de agua, un modelo de efectos mixtos con sitio como efecto aleatorio capturaría mejor la estructura del muestreo.

3 Módulo 3 — Diseño de experimentos: ANOVA y post-ANOVA

3.1 Contexto biológico — pingüinos del archipiélago Palmer

3.2 Datos y exploración

3.3 Verificación de supuestos

3.4 ANOVA de dos vías con interacción

3.5 Post-ANOVA: comparaciones múltiples

3.6 Alternativa no paramétrica: Kruskal-Wallis + Dunn

3.7 Gráfico de interacción especie × isla

3.7.1 Preguntas de interpretación — Módulo 3

4 Resumen ejecutivo

5 Referencias

1.6 Paso 5 — Construcción de la ventana espacial (`owin`)