1 Módulo — 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 <- "EC"

# Sus 10 especies (nombre científico completo: Género especie)
mis_especies <- c(
"Tapirus terrestris" ,
"Ateles belzebuth" ,
"Cebus albifrons",
"Panthera onca",
"Puma concolor",
"Tremarctos ornatus",
"Chelonoidis nigra"
)

# 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:", "Ecuador", "\n")

## País: Ecuador

cat("Especies ingresadas:\n")

## Especies ingresadas:

for (i in seq_along(mis_especies)) cat(sprintf("  %2d. %s\n", i, mis_especies[i]))

##    1. Tapirus terrestris
##    2. Ateles belzebuth
##    3. Cebus albifrons
##    4. Panthera onca
##    5. Puma concolor
##    6. Tremarctos ornatus
##    7. Chelonoidis nigra

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
Tapirus terrestris	Tapirus terrestris	SPECIES	99	OK
Ateles belzebuth	Ateles belzebuth	SPECIES	99	OK
Cebus albifrons	Cebus albifrons	SPECIES	99	OK
Panthera onca	Panthera onca	SPECIES	99	OK
Puma concolor	Puma concolor	SPECIES	99	OK
Tremarctos ornatus	Tremarctos ornatus	SPECIES	99	OK
Chelonoidis nigra	Chelonoidis nigra	SPECIES	98	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

En este paso se descargan registros de ocurrencia desde GBIF para cada especie, filtrando datos con coordenadas válidas y sin errores. Luego, se limpian y unifican en un solo conjunto de datos, generando un resumen del número de registros por especie para su uso en los análisis espaciales posteriores.

# ── 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
Chelonoidis nigra	500
Panthera onca	500
Puma concolor	500
Tapirus terrestris	500
Tremarctos ornatus	500
Ateles belzebuth	321
Cebus albifrons	29

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

## 
## Total de registros: 2850

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.

library(ggplot2)

ggplot(ocurrencias, aes(x = lon, y = lat, color = species)) +
  geom_point(alpha = 0.6, size = 2) +

  labs(
    title = "Registros de presencia en Ecuador",
    subtitle = "Fuente: GBIF (rgbif)",
    x = "Longitud",
    y = "Latitud",
    color = "Especie"
  ) +

  theme_minimal(base_size = 12) +

  theme(
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 12),
    legend.position = "right",
    legend.text = element_text(size = 10)
  )

La gráfica muestra que la mayoría de los registros de las especies se concentran en la región continental de Ecuador, especialmente en la zona andina y amazónica (longitudes entre -80 y -75), lo que sugiere una mayor disponibilidad de datos o una mayor diversidad reportada en estas áreas. En contraste, se observa un grupo claramente aislado hacia el oeste (alrededor de -90 de longitud), correspondiente a especies propias de las islas Galápagos, como Chelonoidis nigra, lo que evidencia una distribución geográficamente separada del territorio continental. Además, varias especies como Panthera onca, Puma concolor y Tremarctos ornatus presentan solapamiento espacial, indicando que comparten hábitats similares o rangos de distribución cercanos, mientras que otras tienen registros más dispersos o limitados, reflejando diferencias ecológicas o en el esfuerzo de muestreo.

1.5 Paso 4 — Mapa de registros crudo

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.2, alpha = 0.6) +
  
  scale_color_viridis_d(option = "turbo", name = "Especie") +
  
  labs(
    title    = paste("Registros de presencia —", nombre_pais),
    subtitle = "Fuente: GBIF (rgbif)",
    x = "Longitud", y = "Latitud"
  ) +
  
  coord_sf() +  # 🔥 ajusta proporciones del mapa
  
  theme_minimal(base_size = 11) +
  
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    plot.subtitle = element_text(size = 10),
    
    legend.position = "bottom",            
    legend.title = element_text(size = 10),
    legend.text = element_text(face = "italic", size = 8),
    
    legend.box = "vertical",
    
    plot.margin = margin(10, 10, 10, 10)   
  ) +
  
  guides(color = guide_legend(nrow = 2))

Este mapa de registros en crudo muestra una clara concentración de ocurrencias en la región continental de Ecuador, especialmente en la zona norte y centro de la Amazonía y la región andina, lo que puede reflejar tanto una mayor diversidad biológica como un mayor esfuerzo de muestreo en estas áreas. Se observa un fuerte solapamiento espacial entre varias especies, como Panthera onca, Puma concolor y Tapirus terrestris, lo que sugiere que comparten hábitats similares o rangos ecológicos cercanos. Por otro lado, los registros de Chelonoidis nigra aparecen completamente aislados hacia el oeste, en las islas Galápagos, evidenciando su distribución geográficamente separada del continente. En general, el patrón indica una distribución desigual de los datos, donde ciertas zonas están más representadas que otras, lo cual es característico de datos en crudo provenientes de bases como GBIF.

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

Rango Y (lat): -4.99 – 1.46

El resultado del rango en longitud (X) y latitud (Y) corresponde a los límites espaciales de la ventana de estudio construida a partir del polígono principal del país. En este caso, el rango de latitud (-4.99 a 1.46) indica que el área analizada se extiende desde el sur hasta el norte del Ecuador continental, abarcando toda su extensión geográfica en el eje norte-sur. Estos valores se obtienen automáticamente a partir de las coordenadas del polígono simplificado y representan los extremos dentro de los cuales se ubican todos los puntos analizados. En conjunto, estos rangos delimitan el espacio donde se construirán los patrones de puntos y se calcularán las densidades, asegurando que el análisis se realice únicamente dentro del área geográfica definida.

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

Especies procesadas exitosamente: 6 / 7

El resultado “6/7” indica que, de las especies analizadas, solo seis cumplieron las condiciones para construir su patrón de puntos y estimar la intensidad. En el código, la función analizar_ppp crea un objeto ppp (patrón de puntos espaciales) a partir de las coordenadas filtradas dentro de la ventana del país, elimina duplicados y luego calcula la densidad mediante un kernel gaussiano (KDE). Sin embargo, si una especie tiene menos de 10 puntos válidos dentro de esa ventana, se omite del análisis. Por eso, aunque inicialmente tuviera muchos registros, una especie pudo quedar fuera tras el filtrado espacial o la depuración, resultando en 6 de 7 especies procesadas exitosamente.

1.8 Paso 7 — Visualización de densidades KDE

Para complementar el análisis de los registros de ocurrencia, se estimaron superficies de densidad espacial para cada especie utilizando los datos obtenidos de GBIF. Estas densidades permiten visualizar las zonas con mayor concentración de registros, proporcionando una aproximación a las áreas de mayor presencia o de mayor esfuerzo de muestreo. Cada mapa representa la intensidad de registros mediante un gradiente de color, donde los tonos más claros indican mayores densidades.

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

Las densidades muestran que la mayoría de las especies se concentran en la región amazónica del Ecuador, con núcleos bien definidos en el oriente. Especies como Tapirus terrestris, Tremarctos ornatus y Panthera onca presentan concentraciones claras, mientras que Puma concolor tiene una distribución más amplia y dispersa. En cambio, Cebus albifrons muestra baja densidad debido a pocos registros, y Ateles belzebuth una distribución más localizada. En general, los patrones reflejan tanto la distribución de las especies como el esfuerzo de muestreo.

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
Ateles belzebuth	306	14.87682	83.046203	21.357133
Cebus albifrons	17	0.88310	1.640981	1.173905
Panthera onca	388	18.98665	92.806057	25.779700
Puma concolor	367	19.59336	48.050113	30.365070
Tapirus terrestris	451	23.96087	80.393297	44.055678
Tremarctos ornatus	486	26.02114	113.765480	40.115110

Los resultados muestran que Tremarctos ornatus y Tapirus terrestris presentan las mayores intensidades promedio y valores altos de λ p75, lo que indica concentraciones importantes y relativamente consistentes de registros. Panthera onca y Puma concolor también presentan valores elevados, aunque con una distribución ligeramente más variable. En contraste, Cebus albifrons muestra valores muy bajos en todas las métricas, lo que se explica por su reducido número de registros, limitando la estimación de densidad. En general, las diferencias observadas reflejan tanto patrones reales de distribución como variaciones en el esfuerzo de muestreo entre especies.

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

## Riqueza mínima: 0

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

## Riqueza máxima: 6 especies

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

La riqueza estimada representa el número de especies que “ocupan” cada celda, entendida como una unidad espacial del raster donde se evalúa la presencia relativa. Una especie se considera presente si su intensidad local supera su mediana (p50), lo que indica que la celda corresponde a una zona de alta densidad dentro de su propia distribución. Por ello, la riqueza en cada celda es la suma de especies que cumplen esta condición.

El valor máximo de 6 especies indica que no existe coincidencia total entre las 10 especies analizadas, evidenciando una distribución espacial heterogénea en Ecuador. En la gráfica, los colores claros representan celdas con mayor riqueza (mayor coincidencia de especies en zonas de alta densidad), mientras que los tonos oscuros indican baja riqueza, asociada a menor solapamiento o distribuciones más restringidas. Este patrón sugiere diferencias en los nichos ecológicos y en la respuesta de las especies a las condiciones ambientales.

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 gráfica muestra la frecuencia de celdas según el número de especies coocurrentes. El valor más alto corresponde a 6 especies (~790 celdas), lo que indica que existen muchas áreas donde varias especies coinciden en zonas de alta densidad relativa. También destaca el valor de 3 especies (~600 celdas), sugiriendo un nivel intermedio de solapamiento bastante común.

En contraste, las celdas con 2 (~190) y 4 (~210 especies) son menos frecuentes, lo que indica que estos niveles de coincidencia son menos representativos en el espacio. Los valores de 1 (~270) y 5 (~280) muestran frecuencias moderadas.

En conjunto, la distribución sugiere que el paisaje está dominado por zonas con alta (6) y media (3) coocurrencia, lo que refuerza la idea de agrupamientos espaciales de especies en ciertas áreas, más que una distribución uniforme.

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?

R// Las especies con mayor intensidad máxima λmax son Tremarctos ornatus y Panthera onca, seguidas por Tapirus terrestris. Esto puede explicarse tanto por su ecología (preferencia por hábitats específicos donde se concentran) como por el esfuerzo de muestreo, ya que ciertas zonas han sido más estudiadas que otras, generando mayores acumulaciones de registros.

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

R// Sí, los hot-spots del mapa de riqueza coinciden en gran medida con áreas reconocidas de alta biodiversidad en Ecuador, especialmente la región amazónica y zonas de transición andina, que son ampliamente conocidas por su alta diversidad de especies y complejidad ecológica.

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

R// El uso de registros de GBIF tiene limitaciones importantes, como el sesgo de muestreo (más datos en zonas accesibles), errores en georreferenciación, identificaciones taxonómicas incorrectas y ausencia de datos en áreas poco estudiadas, lo que puede distorsionar los patrones espaciales reales.

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

R// Si se usa un sigma = 0.3, el mapa mostraría patrones más fragmentados y detallados, con muchos picos locales (más “ruido”). En cambio, con sigma = 3.0, la densidad se suavizaría mucho más, generando patrones más generales y continuos, pero perdiendo detalle fino en la distribución espacial.

2 Módulo — 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

Se realizó un análisis exploratorio de la variable respuesta longitud_concha_mm mediante cuatro gráficos.

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)

El histograma muestra una distribución aproximadamente normal centrada entre 140-150 mm, con las hembras presentando valores ligeramente superiores a los machos. El gráfico de dispersión longitud vs peso revela una relación lineal casi perfecta entre ambas variables, anticipando problemas de multicolinealidad. El gráfico longitud vs temperatura muestra una línea loess prácticamente horizontal, indicando que la temperatura del agua no tiene relación con el tamaño corporal. Finalmente, el violin plot confirma que las hembras tienen una mediana de longitud ligeramente mayor que los machos, con distribuciones simétricas en ambos sexos.

2.3.1 Correlaciones

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("RdYlBu", 10),
         title = "Correlaciones entre variables morfométricas",
         mar = c(0, 0, 2, 0))

La matriz de correlaciones de Pearson reveló dos pares de variables con multicolinealidad crítica (|r| > 0.90): longitud_concha_mm peso_g (r = 0.98) y longitud_concha_mm ancho_concha_mm (r = 0.91). Estas variables fueron eliminadas automáticamente en el paso de preprocesamiento mediante step_corr(threshold = 0.90). Las variables temp_agua_c y profundidad_m mostraron correlaciones cercanas a cero con todas las demás variables, confirmando su independencia como predictores ambientales.

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

## Entrenamiento: 210 | Prueba: 70

Predictores tras pre proceso:

ancho_concha_mm, altura_concha_mm, largo_cabeza_mm, temp_agua_c, profundidad_m, edad_años, sexo_Macho

La tabla de coeficientes del modelo de regresión lineal ajustado sobre los datos de entrenamiento muestra que los predictores morfológicos ancho_concha_mm (β = 11.42), altura_concha_mm (β = 5.88) y largo_cabeza_mm (β = 4.34) son los más influyentes y altamente significativos (p < 0.001). La variable edad_años también resultó significativa (p = 0.003). En contraste, temp_agua_c, profundidad_m y sexo_Macho no alcanzaron significancia estadística (p > 0.05). El intercepto de 141.54 mm corresponde a la longitud media predicha cuando todos los predictores están en su valor promedio.

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

Las métricas de bondad de ajuste en los datos de entrenamiento indican un modelo con excelente desempeño: R² = 0.9174, lo que significa que el modelo explica el 91.74% de la variabilidad en la longitud de la concha. El R² ajustado (0.9145) es prácticamente idéntico al R², confirmando que todos los predictores incluidos aportan información relevante y no inflan el R² artificialmente. El RMSE de 6.51 mm representa un error promedio del 4.7% respecto al rango total de longitudes observadas

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

La evaluación sobre los 70 datos de prueba, nunca vistos por el modelo, arrojó un RMSE de 6.35 mm, MAE de 5.16 mm y R² de 0.8935. La mínima diferencia entre el R² de entrenamiento (0.9174) y el de prueba (0.8935) — apenas 0.024 — confirma que el modelo generaliza correctamente a datos nuevos y no presenta sobreajuste.

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)

El gráfico de valores observados versus predichos muestra que la mayoría de los puntos se ubican sobre la línea de predicción perfecta (línea punteada), con errores predominantemente pequeños (puntos en color morado oscuro). Se identifican 2-3 casos con errores mayores a 15 mm (puntos amarillos), pero son minoría. La distribución homogénea de los errores a lo largo de todo el rango de longitudes (90-180 mm) sugiere que el modelo predice igual de bien para tortugas pequeñas que para tortugas grandes.

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

La validación cruzada 10-fold, que divide los datos de entrenamiento en 10 particiones y evalúa el modelo en cada una, arrojó un RMSE medio de 6.52 mm (±0.34), MAE de 5.26 mm (±0.27) y R² de 0.9116 (±0.01). Los errores estándar pequeños en todas las métricas indican que el desempeño del modelo es consistente independientemente de qué datos se usen para entrenar, confirmando su robustez y estabilidad.

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

Los cuatro gráficos de diagnóstico evalúan el cumplimiento de los supuestos del modelo lineal. El gráfico Residuals vs Fitted confirma linealidad con residuos distribuidos aleatoriamente alrededor de cero. El Q-Q plot muestra que los residuos siguen una distribución normal con leves desviaciones solo en las colas. El gráfico Scale-Location revela una leve heterocedasticidad, con residuos algo mayores en tortugas de menor tamaño, aunque no compromete la validez del modelo. Finalmente, la distancia de Cook muestra que los casos más influyentes (8, 105, 123) están muy por debajo del umbral crítico de 0.5, descartando observaciones problemáticas.

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?

R//Los predictores más influyentes fueron el ancho y la altura de la concha, el largo de la cabeza y la edad, todos altamente significativos. Esto es biológicamente esperable porque el crecimiento de la concha es proporcional en sus dimensiones y las tortugas crecen con la edad, mientras que variables ambientales como temperatura y profundidad no influyen directamente en el tamaño corporal.

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

R// No hay indicios de sobreajuste, ya que el R² de entrenamiento (0.9174) y de prueba (0.8935) son muy cercanos, con una diferencia pequeña. Además, la validación cruzada confirma la estabilidad del modelo, lo que indica que generaliza bien y no está memorizando los datos.

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

R// La RMSE en prueba es más informativa porque evalúa el modelo con datos que no fueron usados en el entrenamiento, reflejando su desempeño real. En cambio, la RMSE de entrenamiento suele ser optimista. En este caso, la RMSE de prueba incluso es ligeramente menor, lo que sugiere un buen ajuste y capacidad predictiva.

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

R// Los supuestos del modelo se cumplen en general: hay linealidad (residuos distribuidos aleatoriamente), normalidad (Q-Q alineado), no hay puntos influyentes problemáticos (Cook bajo), y solo se observa una leve heterocedasticidad, que no compromete significativamente la validez del modelo.

¿Qué mejoras propondría (transformaciones, interacciones, otros predictores)?

R// Se podrían proponer mejoras como aplicar una transformación logarítmica para reducir la heterocedasticidad, incluir interacciones como sexo × edad, construir un índice de condición corporal en lugar de eliminar variables correlacionadas, y explorar modelos no lineales que capturen mejor el crecimiento biológico.

3 Módulo — 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")

## Dimensiones: 333 × 8

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)

La tabla muestra que las especies no se distribuyen de manera uniforme entre las islas: Adelie está presente en Biscoe, Dream y Torgersen (con mayor abundancia en Dream), Chinstrap solo aparece en Dream y Gentoo exclusivamente en Biscoe. La gráfica complementa esto al evidenciar que también existen diferencias claras en la masa corporal entre especies, donde Gentoo presenta los valores más altos, mientras que Adelie y Chinstrap tienen masas menores y similares. Además, al considerar la isla, se observa que estas diferencias se mantienen dentro de cada ubicación, destacando nuevamente Gentoo en Biscoe. En conjunto, esto indica que tanto la distribución de individuos como la masa corporal varían según la especie y la isla, sugiriendo un efecto conjunto de estos factores, lo cual es clave para el análisis ANOVA.

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)

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

La prueba de Shapiro-Wilk reveló normalidad en Chinstrap (p = 0.56) y Gentoo (p = 0.26), mientras que Adelie mostró un p = 0.042 borderline. Sin embargo, dado que el estadístico W = 0.98 está muy cercano a 1 y el QQ-plot confirma normalidad visual, este rechazo se atribuye a la sensibilidad del test con muestras grandes (n = 146). La prueba de Levene resultó significativa (p = 0.006), indicando heterocedasticidad entre grupos, lo que refuerza el uso del test no paramétrico Kruskal-Wallis como análisis complementario.

3.3.1 QQ-plot por 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()

Los tres QQ-plots muestran puntos siguiendo la diagonal teórica con gran precisión, confirmando normalidad visual en las tres especies. Las leves desviaciones en las colas extremas de Adelie son menores y no comprometen el supuesto de normalidad. El QQ-plot es más informativo que el p-valor de Shapiro-Wilk cuando el tamaño muestral es grande.

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

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

El ANOVA de dos vías reveló que la especie es la única fuente de variación significativa (F = 339.83, p < 0.001), explicando la mayor parte de la variabilidad en masa corporal. La isla no resultó significativa (F = 0.005, p = 0.995), lo cual es esperable dado el diseño desbalanceado donde cada especie habita islas distintas. El término de interacción especie × isla fue eliminado automáticamente por R debido a las celdas vacías en el diseño.

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

tukey_df <- as.data.frame(tukey_res$especie) %>%
  tibble::rownames_to_column("comparacion")

ggplot(tukey_df, aes(x = diff, y = comparacion)) +
  geom_vline(xintercept = 0, linetype = "dashed",
             color = "red", linewidth = 0.8) +
  geom_errorbarh(aes(xmin = lwr, xmax = upr),
                 height = 0.2, linewidth = 0.9, color = "#3F51B5") +
  geom_point(size = 3.5, color = "#3F51B5") +
  labs(title = "Tukey HSD — diferencias entre especies",
       x = "Diferencia de medias (g)", y = "Comparación") +
  theme_minimal(base_size = 13)

El gráfico muestra que Gentoo difiere significativamente de Adelie (+1386g) y Chinstrap (+1359g), cuyos intervalos de confianza no cruzan la línea roja de cero. El par Chinstrap-Adelie sí cruza el cero, confirmando que no difieren significativamente.

modelo_simple <- aov(masa_g ~ especie, data = pinguinos)
emm <- emmeans(modelo_simple, ~ especie)

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

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

Tanto Tukey como Bonferroni confirmaron que Gentoo difiere significativamente de las otras dos especies (p = 0.000), mientras que Adelie y Chinstrap no difieren entre sí (p Tukey = 0.917, p Bonferroni = 1.000).

3.5.1 Gráfico de medias

medias_g <- pinguinos %>%
  group_by(especie) %>%
  summarise(media = mean(masa_g), se = sd(masa_g)/sqrt(n())) %>%
  mutate(letra = c("b", "b", "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 con sus intervalos de confianza al 95%. Adelie y Chinstrap comparten la letra b (no difieren entre sí), mientras que Gentoo tiene la letra a (difiere significativamente de ambas).

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)

El test no paramétrico de Kruskal-Wallis confirmó diferencias significativas entre especies en masa corporal. El post-hoc de Dunn con corrección Bonferroni reprodujo exactamente las mismas conclusiones que el ANOVA paramétrico: Gentoo difiere significativamente de Adelie (Z = 13.61, p = 0.000) y Chinstrap (Z = 10.73, p = 0.000), mientras que Adelie y Chinstrap no difieren (Z = 0.34, p = 1.000). La concordancia entre ambos métodos refuerza la robustez de los resultados ante la violación 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)

El gráfico de interacción muestra que Adelie mantiene una masa corporal prácticamente constante (~3650-3700g) en las tres islas, indicando ausencia de efecto de isla sobre esta especie. Chinstrap y Gentoo aparecen como puntos únicos al habitar exclusivamente en Dream y Biscoe respectivamente, lo que impide evaluar la interacción formalmente. Este diseño confirma que las diferencias en masa corporal son atribuibles exclusivamente a la especie y no a la isla de habitación.

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?

R// Los supuestos se cumplen parcialmente. La normalidad se confirma visualmente mediante QQ-plots para las tres especies, y el rechazo borderline de Shapiro-Wilk en Adelie (p = 0.042) se atribuye a la sensibilidad del test con muestras grandes (n = 146) y no a una desviación real. Sin embargo, la prueba de Levene sí rechaza la homocedasticidad (p = 0.006), indicando que las varianzas difieren entre especies. A pesar de esto, con n = 333 el ANOVA es robusto a violaciones moderadas de homocedasticidad, y la concordancia con el test no paramétrico Kruskal-Wallis valida las conclusiones obtenidas.

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

R// La especie es la única fuente de variación significativa (F = 339.83, p < 0.001), mientras que la isla no lo es (F = 0.005, p = 0.995). Biológicamente esto indica que las diferencias en masa corporal son una característica intrínseca de cada especie, independiente del hábitat geográfico.

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

R// Ambos métodos controlan el error tipo I familiar (FWER), pero difieren en su estrategia. Tukey HSD está diseñado específicamente para comparaciones por pares y es más potente que Bonferroni cuando se comparan todos los pares posibles, reduciendo el error tipo II. Bonferroni es más conservador porque divide α entre el número de comparaciones, aumentando la probabilidad de falsos negativos. En este caso ambos coincidieron en sus conclusiones. Se prefiere Tukey cuando se comparan todos los pares de grupos, y Bonferroni cuando se realizan pocas comparaciones planeadas de antemano o cuando se requiere mayor control del error tipo I.

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

R// Sí, ambos métodos coinciden exactamente en sus conclusiones. Gentoo difiere significativamente de Adelie y Chinstrap, mientras que Adelie y Chinstrap no difieren entre sí. Esta concordancia ocurre porque las diferencias entre especies son muy grandes, tanto con los valores originales (ANOVA) o sus rangos (Kruskal-Wallis). Además, con muestras grandes como las de este dataset, el ANOVA es robusto a las violaciones moderadas de homocedasticidad detectadas por Levene.

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

R// El gráfico de interacción revela que Adelie es la única especie con presencia en las tres islas, manteniendo una masa corporal prácticamente constante (~3650-3700g) independientemente de la isla, lo que sugiere que las condiciones ambientales entre islas no afectan su tamaño corporal. Chinstrap aparece exclusivamente en Dream y Gentoo exclusivamente en Biscoe, como puntos únicos sin posibilidad de trazar líneas de interacción. Esto refleja una distribución geográfica especializada por especie que genera un diseño estadístico desbalanceado, limitando la capacidad de evaluar formalmente si el efecto de la especie sobre la masa varía según la isla.

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.

## R version 4.5.2 Patched (2026-01-31 r89366 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.4.0        knitr_1.51              rstatix_0.7.3          
##  [4] emmeans_2.0.3           ggpubr_0.6.3            palmerpenguins_0.1.1   
##  [7] corrplot_0.95           yardstick_1.3.2         workflowsets_1.1.1     
## [10] workflows_1.3.0         tune_2.0.1              tidyr_1.3.2            
## [13] tailor_0.1.0            rsample_1.3.2           recipes_1.3.2          
## [16] parsnip_1.5.0           modeldata_1.5.1         infer_1.1.0            
## [19] dials_1.4.2             broom_1.0.12            tidymodels_1.4.1       
## [22] gridExtra_2.3           scales_1.4.0            viridis_0.6.5          
## [25] viridisLite_0.4.2       ggplot2_4.0.2           purrr_1.2.1            
## [28] dplyr_1.2.0             countrycode_1.7.0       rnaturalearthdata_1.0.0
## [31] rnaturalearth_1.2.0     sf_1.0-24               terra_1.8-93           
## [34] spatstat_3.5-1          spatstat.linnet_3.4-1   spatstat.model_3.6-1   
## [37] rpart_4.1.24            spatstat.explore_3.7-0  nlme_3.1-168           
## [40] spatstat.random_3.4-4   spatstat.geom_3.7-0     spatstat.univar_3.1-6  
## [43] spatstat.data_3.1-9     rgbif_3.8.4            
## 
## 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-1  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          plyr_1.8.9            lubridate_1.9.5      
##  [22] cachem_1.1.0          whisker_0.4.1         lifecycle_1.0.5      
##  [25] pkgconfig_2.0.3       Matrix_1.7-4          R6_2.6.1             
##  [28] fastmap_1.2.0         future_1.70.0         digest_0.6.39        
##  [31] furrr_0.4.0           tensor_1.5.1          textshaping_1.0.4    
##  [34] labeling_0.4.3        urltools_1.7.3.1      spatstat.sparse_3.1-0
##  [37] timechange_0.4.0      httr_1.4.7            polyclip_1.10-7      
##  [40] abind_1.4-8           mgcv_1.9-4            compiler_4.5.2       
##  [43] proxy_0.4-29          withr_3.0.2           S7_0.2.1             
##  [46] backports_1.5.0       carData_3.0-6         DBI_1.2.3            
##  [49] ggsignif_0.6.4        MASS_7.3-65           lava_1.9.0           
##  [52] classInt_0.4-11       oai_0.4.0             tools_4.5.2          
##  [55] units_1.0-0           otel_0.2.0            future.apply_1.20.2  
##  [58] nnet_7.3-20           goftest_1.2-3         glue_1.8.0           
##  [61] modelenv_0.2.0        grid_4.5.2            generics_0.1.4       
##  [64] gtable_0.3.6          class_7.3-23          data.table_1.18.2.1  
##  [67] xml2_1.5.2            car_3.1-5             pillar_1.11.1        
##  [70] stringr_1.6.0         splines_4.5.2         lhs_1.2.1            
##  [73] lattice_0.22-7        survival_3.8-6        deldir_2.0-4         
##  [76] tidyselect_1.2.1      svglite_2.2.2         crul_1.6.0           
##  [79] xfun_0.56             hardhat_1.4.3         timeDate_4052.112    
##  [82] stringi_1.8.7         DiceDesign_1.10       lazyeval_0.2.2       
##  [85] yaml_2.3.12           httpcode_0.3.0        evaluate_1.0.5       
##  [88] codetools_0.2-20      tibble_3.3.1          cli_3.6.5            
##  [91] systemfonts_1.3.1     jquerylib_0.1.4       Rcpp_1.1.1           
##  [94] globals_0.19.1        triebeard_0.4.1       parallel_4.5.2       
##  [97] gower_1.0.2           GPfit_1.0-9           listenv_0.10.1       
## [100] mvtnorm_1.3-6         ipred_0.9-15          prodlim_2026.03.11   
## [103] e1071_1.7-17          rlang_1.1.7

Taller generado con R Markdown — 2026-04-23 21:18

Taller de Bioestadística Aplicada en R

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

Caroline Muñoz y Mariana Simmons

23 de abril de 2026

1 Módulo — 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 crudo

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 — Regresión lineal múltiple con enfoque predictivo

2.1 Contexto biológico

2.2 Datos

2.3 Análisis exploratorio

2.3.1 Correlaciones

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

3 Módulo — 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.3.1 QQ-plot por especie

3.4 ANOVA de dos vías con interacción

3.5 Post-ANOVA: comparaciones múltiples

3.5.1 Gráfico de medias

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

Caroline Muñoz y Mariana Simmons

23 de abril de 2026

1 Módulo — 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 crudo

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 — Regresión lineal múltiple con enfoque predictivo

2.1 Contexto biológico

2.2 Datos

2.3 Análisis exploratorio

2.3.1 Correlaciones

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

3 Módulo — 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.3.1 QQ-plot por especie

3.4 ANOVA de dos vías con interacción

3.5 Post-ANOVA: comparaciones múltiples

3.5.1 Gráfico de medias

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