Instrucciones generales. Este taller tiene tres módulos. En el Módulo 1 usted elige 10 especies de un país a su elección y completa los bloques marcados con >>> ESTUDIANTE. Los módulos 2 y 3 tienen datos incluidos y son completamente reproducibles sin modificaciones. Ejecute cada bloque en orden, analice los resultados y responda las preguntas de interpretación al final de cada sección.

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.

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

# 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(
  "Bothrops asper",
  "Bothrops atrox",
  "Bothrops punctatus",
  "Micrurus mipartitus",
  "Micrurus dumerilii",
  "Micrurus dissoleucus",
  "Porthidium lansbergii",
  "Bothriechis schlegelii",
  "Crotalus durissus",
  "Bothrocophias hyoprora"
)

# 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]))

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
Bothrops asper Bothrops asper SPECIES 99 OK
Bothrops atrox Bothrops atrox SPECIES 99 OK
Bothrops punctatus Bothrops punctatus SPECIES 99 OK
Micrurus mipartitus Micrurus mipartitus SPECIES 99 OK
Micrurus dumerilii Micrurus dumerilii SPECIES 99 OK
Micrurus dissoleucus Micrurus dissoleucus SPECIES 99 OK
Porthidium lansbergii Porthidium lansbergii SPECIES 99 OK
Bothriechis schlegelii Bothriechis schlegelii SPECIES 99 OK
Crotalus durissus Crotalus durissus SPECIES 97 OK
Bothrocophias hyoprora Bothrocophias hyoprora SPECIES 99 OK

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
Bothrops asper 500
Bothrops atrox 500
Micrurus dumerilii 500
Micrurus mipartitus 500
Porthidium lansbergii 500
Crotalus durissus 352
Bothriechis schlegelii 295
Micrurus dissoleucus 126
Bothrops punctatus 94
Bothrocophias hyoprora 73 
cat(sprintf("\nTotal de registros: %d\n", nrow(ocurrencias)))
## 
## Total de registros: 3440

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

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

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

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

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
Bothriechis schlegelii 232 2.642460 20.194355 3.210638
Bothrocophias hyoprora 64 0.624741 16.638135 0.290547
Bothrops asper 449 4.915995 31.843007 8.731600
Bothrops atrox 448 4.732505 34.561332 5.965446
Bothrops punctatus 75 0.794289 7.909455 0.817628
Crotalus durissus 308 3.178813 44.799050 1.887858
Micrurus dissoleucus 92 0.906146 16.126499 0.271418
Micrurus dumerilii 404 4.451644 23.936863 7.409379
Micrurus mipartitus 429 4.817903 25.419250 7.778107
Porthidium lansbergii 449 4.524251 61.935366 5.175008

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)

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

  # 🔥 ARREGLA EL VOLTEO
  rst <- flip(rst, direction = "vertical")

  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")
riqueza_df <- as.data.frame(riqueza, xy = TRUE) %>%
  setNames(c("x", "y", "riqueza")) %>%
  filter(!is.na(riqueza), riqueza > 0)

ggplot() +
  geom_raster(data = riqueza_df,
              aes(x = x, y = y, fill = riqueza)) +
  geom_sf(data = pais_sf, fill = NA, color = "white", linewidth = 0.5) +
  scale_fill_gradientn(
    colors = c("#0d0221","#3a1c71","#d76d77","#ffaf7b","#fffde4"),
    name   = "Riqueza\nestimada",
    breaks = pretty(riqueza_df$riqueza, n = 5)
  ) +
  labs(
    title    = paste("Riqueza estimada de especies —", nombre_pais),
    subtitle = paste(length(resultados_ppp),
                     "especies · KDE (sigma =", SIGMA_KDE, "°) · spatstat + terra"),
    x = "Longitud", y = "Latitud",
    caption  = "Álgebra de rasters: presencia = λ̂(x,y) > mediana(λ̂) por especie"
  ) +
  theme_void(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 10, color = "gray40"),
    legend.position = "right"
  )

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)

1.9.1 Preguntas de interpretación — Módulo 1

  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?

Porthidium lansbergii tuvo la intensidad máxima más alta (61.93), concentrada en el norte de Colombia. Esto puede reflejar tanto su ecología real (especie de tierras bajas y zonas áridas del Caribe) como un sesgo de muestreo: las regiones norte y occidente del país tienen mayor accesibilidad y tradición de colecta herpetológica. Crotalus durissus (λ = 44.8) también es llamativo porque habita zonas abiertas y áridas (Llanos, Caribe) que suelen estar mejor prospectadas. La combinación de ecología de hábitat abierto + esfuerzo acumulado en GBIF inflaciona localmente la intensidad estimada.

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

Los hot-spots estimados se ubican en el noroccidente y la zona Andina central, áreas que efectivamente coinciden con la alta biodiversidad de la región del Chocó biogeográfico y los Andes colombianos, reconocidos mundialmente como puntos calientes de biodiversidad. Sin embargo, regiones igualmente megadiversas como la Amazonía y la Orinoquía aparecen con baja riqueza estimada, lo que refleja subregistro en GBIF más que ausencia real de especies. El mapa captura bien las zonas mejor muestreadas, no necesariamente las más diversas.

  1. ¿Qué limitaciones tiene usar registros de GBIF para análisis de patrones de puntos?
  1. Sesgo de detección: los registros se acumulan donde hay más colectores, instituciones o ciudadanos científicos, no donde hay más individuos. (2) Registros duplicados de la misma observación reportada por múltiples fuentes inflaciona artificialmente la intensidad local. (3) Resolución temporal variable: datos históricos y recientes se mezclan sin indicar cambios poblacionales. (4) Errores de georeferenciación pueden desplazar puntos al mar o a coordenadas imposibles. (5) La intensidad KDE estima densidad de registros, no de individuos reales, por lo que refleja el esfuerzo de muestreo tanto como la distribución real.
  1. ¿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 el kernel sería muy estrecho: cada punto generaría un pico muy concentrado, el mapa mostraría “manchas” aisladas donde se concentran registros y zonas de cero entre ellas. Capturaría bien la heterogeneidad local pero sería sensible a registros erróneos. Con sigma = 3.0 el suavizado sería excesivo: las densidades se difuminarían sobre grandes áreas, borrando patrones de distribución reales y haciendo que especies de distribuciones muy distintas aparezcan como si tuvieran rangos similares. El valor actual (1°) representa un balance razonable para el territorio colombiano.

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)

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

library(caret)

set.seed(123)

# ── Split estratificado por sexo ─────────────────────────────────────
train <- datos_tortugas %>%
  group_by(sexo) %>%
  slice_sample(prop = 0.75) %>%
  ungroup()

test <- anti_join(datos_tortugas, train)

cat(sprintf("Entrenamiento: %d | Prueba: %d\n", nrow(train), nrow(test)))
# ── One-hot encoding (equivalente a step_dummy) ─────────────────────
train_proc <- model.matrix(longitud_concha_mm ~ ., data = train) %>%
  as.data.frame()

test_proc <- model.matrix(longitud_concha_mm ~ ., data = test) %>%
  as.data.frame()

# Quitar intercepto
train_proc <- train_proc %>% select(-`(Intercept)`)
test_proc  <- test_proc  %>% select(-`(Intercept)`)

# ── Separar respuesta
y_train <- train$longitud_concha_mm
y_test  <- test$longitud_concha_mm

# ── Normalización (equivalente a step_normalize) ─────────────────────
num_cols <- names(train_proc)

means <- sapply(train_proc, mean)
sds   <- sapply(train_proc, sd)

train_proc <- as.data.frame(scale(train_proc, center = means, scale = sds))
test_proc  <- as.data.frame(scale(test_proc, center = means, scale = sds))

# ── Eliminar alta correlación (equivalente a step_corr) ──────────────
cor_mat <- cor(train_proc)
high_corr <- findCorrelation(cor_mat, cutoff = 0.90)  # de caret

if (length(high_corr) > 0) {
  train_proc <- train_proc[, -high_corr]
  test_proc  <- test_proc[, -high_corr]
}

# ── Volver a añadir la respuesta
train_proc$longitud_concha_mm <- y_train
test_proc$longitud_concha_mm  <- y_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 ***
sexoMacho -0.8371 0.4593 -1.8228 0.0698
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 ** 
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)

k <- 10
n <- nrow(train_proc)

folds <- sample(rep(1:k, length.out = n))

rmse_vec <- c()
rsq_vec  <- c()
mae_vec  <- c()

for (i in 1:k) {
  
  train_fold <- train_proc[folds != i, ]
  val_fold   <- train_proc[folds == i, ]
  
  # Modelo
  modelo <- lm(longitud_concha_mm ~ ., data = train_fold)
  
  # Predicción
  pred <- predict(modelo, newdata = val_fold)
  real <- val_fold$longitud_concha_mm
  
  # Métricas
  rmse <- sqrt(mean((real - pred)^2))
  mae  <- mean(abs(real - pred))
  rsq  <- cor(real, pred)^2
  
  rmse_vec <- c(rmse_vec, rmse)
  mae_vec  <- c(mae_vec, mae)
  rsq_vec  <- c(rsq_vec, rsq)
}

# Resumen
cv_resumen <- data.frame(
  Métrica = c("RMSE", "R²", "MAE"),
  Media   = c(mean(rmse_vec), mean(rsq_vec), mean(mae_vec)),
  Error_est = c(sd(rmse_vec)/sqrt(k),
                sd(rsq_vec)/sqrt(k),
                sd(mae_vec)/sqrt(k)),
  Folds = k
)

cv_resumen %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  kable(caption = "Validación cruzada 10-fold") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)
Validación cruzada 10-fold
Métrica Media Error_est Folds
RMSE 6.5774 0.4297 10
R² 0.9077 0.0137 10
MAE 5.3429 0.3723 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

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

Los predictores más influyentes son el ancho de concha (β = 11.42), la altura de concha (β = 5.88) y el largo de cabeza (β = 4.34), todos con p < 0.001. Esto es biológicamente muy esperable: en tortugas, las dimensiones de la concha están fuertemente correlacionadas entre sí (alometría positiva), por lo que conocer el ancho o la altura permite predecir con alta precisión la longitud. La edad también contribuye significativamente (β = 1.85, p = 0.003), consistente con el crecimiento continuo de tortugas. En contraste, la temperatura del agua (p = 0.53) y la profundidad (p = 0.08) no son predictores significativos, lo que tiene sentido ya que son variables ambientales externas con menos relación directa con la morfología individual.

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

El R² en entrenamiento (0.917) y en datos de prueba (0.893) son muy similares, con una diferencia de apenas 0.024. Esta pequeña caída es normal y esperada al evaluar en datos no vistos. No hay indicios de sobreajuste severo. La validación cruzada 10-fold confirma esto con R² medio = 0.908 y error estándar muy bajo (0.014), lo que indica que el modelo generaliza de forma estable y robusta a nuevas observaciones.

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

La RMSE en entrenamiento mide qué tan bien el modelo se ajusta a los datos con los que fue construido; siempre será optimista porque el modelo “conoce” esos datos. La RMSE en prueba (6.35 mm) mide el error en datos completamente nuevos, lo que refleja la capacidad de generalización real del modelo. Un modelo sobreajustado puede tener RMSE de entrenamiento muy baja pero alta en prueba. En este caso ambas son similares (6.51 vs 6.35 mm), confirmando que el modelo tiene buena capacidad predictiva en situaciones reales.

  1. Analice los cuatro gráficos de diagnóstico: ¿se cumplen los supuestos del modelo lineal?

Residuals vs Fitted: los puntos se distribuyen aleatoriamente alrededor de cero, cumpliendo linealidad. Q-Q Residuals: los puntos siguen la línea diagonal, confirmando normalidad de los residuos. Scale-Location: la línea es aproximadamente horizontal aunque con ligera tendencia, sugiriendo homocedasticidad aceptable pero con posible heterocedasticidad leve en los extremos. Cook’s distance: la observación 81 tiene la mayor influencia, seguida de 71 y 77, aunque ninguna supera valores críticos que comprometan el modelo.

  1. ¿Qué mejoras propondría (transformaciones, interacciones, otros predictores)?
  1. Transformación logarítmica de la variable respuesta y/o predictores de peso/longitud para linealizar relaciones alométricas y mejorar homocedasticidad. (2) Incluir interacciones como sexo × edad o ancho × altura, ya que el dimorfismo sexual puede afectar la relación entre predictores y longitud de forma diferencial. (3) Eliminar predictores no significativos (temperatura, profundidad) con selección de variables (AIC stepwise o LASSO) para obtener un modelo más parsimonioso. (4) Investigar la observación 81 para verificar si es un error de medición o un individuo biológicamente atípico.

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

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

3.4 ANOVA de dos vías con interacción

modelo_anova <- aov(masa_g ~ especie * isla, data = pinguinos)
summary(modelo_anova)
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

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)
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))
emm <- emmeans(modelo_anova, ~ especie)

# Tukey via emmeans
contrast(emm, 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 z.ratio p.value
Adelie - Chinstrap NA NA NA NA NA
Adelie - Gentoo NA NA NA NA NA
Chinstrap - Gentoo NA NA NA NA NA 
# Bonferroni (más conservador)
contrast(emm, 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 z.ratio p.value
Adelie - Chinstrap NA NA NA NA NA
Adelie - Gentoo NA NA NA NA NA
Chinstrap - Gentoo NA NA NA NA NA 
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)

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 ****

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)

3.7.1 Preguntas de interpretación — Módulo 3

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

Los supuestos se cumplen de forma parcial. Chinstrap (p = 0.56) y Gentoo (p = 0.26) pasan el test de Shapiro-Wilk sin problemas. Adelie rechaza la normalidad marginalmente (p = 0.042), lo que puede atribuirse a su mayor tamaño muestral (n = 146) que hace al test más sensible. El ANOVA es robusto a desviaciones moderadas de normalidad cuando n es grande, por lo que los resultados siguen siendo válidos. Esto se confirma con el análisis no paramétrico de Kruskal-Wallis que llega a las mismas conclusiones.

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

Solo la especie es una fuente de variación significativa (F = 339.83, p < 0.001). La isla no tiene efecto significativo (p = 0.995), y la interacción especie × isla tampoco es significativa, lo que biológicamente indica que las diferencias de masa corporal entre especies son consistentes independientemente de la isla en que habitan. El peso de Gentoo es mayor que Adelie y Chinstrap independientemente del sitio, lo que sugiere que estas diferencias son intrínsecas de la especie (genética, morfología, dieta) y no están moduladas por el ambiente local de cada isla.

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

Tukey HSD está diseñado específicamente para comparaciones entre todos los pares de medias y controla la tasa de error familiar de forma más eficiente; es preferible cuando se quieren comparar todos los grupos entre sí con mayor potencia estadística. Bonferroni divide el α entre el número de comparaciones; es más conservador (mayor riesgo de falsos negativos) pero más flexible, aplicable a cualquier conjunto de hipótesis a priori, no solo comparaciones de medias. En este caso con 3 grupos y pocas comparaciones (3 pares), la diferencia práctica es mínima y ambos deberían coincidir en las conclusiones.

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

Sí coinciden. El post-hoc de Dunn con corrección Bonferroni indica que Adelie vs Gentoo (Z = 13.61, p < 0.0001) y Chinstrap vs Gentoo (Z = 10.73, p < 0.0001) son significativamente diferentes, mientras que Adelie vs Chinstrap no lo es (p = 1.0). Esto replica exactamente los resultados de Tukey HSD. La concordancia se explica porque las distribuciones no se apartan demasiado de la normalidad y los tamaños muestrales son grandes, condiciones bajo las cuales ANOVA y Kruskal-Wallis tienen potencia similar.

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

El gráfico muestra que las líneas de Adelie y Chinstrap son prácticamente horizontales y paralelas entre las islas donde coexisten (Biscoe y Dream), confirmando la ausencia de interacción: el efecto de la especie sobre el peso no cambia según la isla. Gentoo solo aparece en la isla Biscoe en este estudio, lo que imposibilita comparar su comportamiento entre islas pero no afecta las conclusiones generales. Biológicamente, esto sugiere que las diferencias de masa corporal entre especies son un rasgo estable de cada especie y no están moduladas por condiciones locales de cada isla del archipiélago Palmer.

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.

Taller generado con R Markdown — 2026-04-28 17:03