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(
  "Panthera onca",
  "Puma concolor",
  "Dicotyles tajacu",
  "Herpailurus yagouaroundi",
  "Leopardus pardalis",
  "Leopardus wiedii",
  "Leopardus tigrinus",
  "Tremarctos ornatus",
  "Odocoileus virginianus",
  "Myrmecophaga tridactyla"
)

# 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")
## País: CO
cat("Especies ingresadas:\n")
## Especies ingresadas:
for (i in seq_along(mis_especies)) cat(sprintf("  %2d. %s\n", i, mis_especies[i]))
##    1. Panthera onca
##    2. Puma concolor
##    3. Dicotyles tajacu
##    4. Herpailurus yagouaroundi
##    5. Leopardus pardalis
##    6. Leopardus wiedii
##    7. Leopardus tigrinus
##    8. Tremarctos ornatus
##    9. Odocoileus virginianus
##   10. Myrmecophaga tridactyla

Se escogió como país a Colombia, en cuanto a las (10) especies, estos son mamíferos de tamaño mediano-grande del neotropico, encontrados en bosques tropicales y andinos, en esta lista se encuentra una gran diversidad funcional; es decir especies arboricolas, terrestres, generalistas y especialistas con dietas diversas ya que hay carnívoros, herbívoros e insectívoros. representando así un ecosistema complejo que muestra diferentes nichos ecológicos.

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

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

tabla_tax <- map_dfr(mis_especies, verificar_sp)

kable(tabla_tax,
      caption = "Verificación taxonómica — backbone GBIF") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  column_spec(5, color = ifelse(tabla_tax$estado == "OK",
                                 "darkgreen", "red"),
              bold = TRUE)
Verificación taxonómica — backbone GBIF
ingresado aceptado rango confianza estado
Panthera onca Panthera onca SPECIES 99 OK
Puma concolor Puma concolor SPECIES 99 OK
Dicotyles tajacu Dicotyles tajacu SPECIES 98 OK
Herpailurus yagouaroundi Herpailurus yagouaroundi SPECIES 98 OK
Leopardus pardalis Leopardus pardalis SPECIES 99 OK
Leopardus wiedii Leopardus wiedii SPECIES 99 OK
Leopardus tigrinus Leopardus tigrinus SPECIES 99 OK
Tremarctos ornatus Tremarctos ornatus SPECIES 99 OK
Odocoileus virginianus Odocoileus virginianus SPECIES 99 OK
Myrmecophaga tridactyla Myrmecophaga tridactyla 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.

Aquí se validan los nombres de las especies, detecta errores taxonómicos, los estandariza y muestra los que necesitan ser revisados, ademas, muestra un nivel de confianza para saber que tan seguro esta el sistema de que el nombre ingresado es el correcto.

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
Dicotyles tajacu 500
Herpailurus yagouaroundi 500
Leopardus pardalis 500
Leopardus tigrinus 500
Myrmecophaga tridactyla 500
Odocoileus virginianus 500
Panthera onca 500
Puma concolor 500
Tremarctos ornatus 500
Leopardus wiedii 447 
cat(sprintf("\nTotal de registros: %d\n", nrow(ocurrencias)))
## 
## Total de registros: 4947

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.

Se tomaron 4947 registros en total en donde en casi todas las especies hubo un numero de 500 registros excepto por Leopardus wiedii el cual tuvo 447 registros.

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

En este mapa se puede observar el registro de presencia de varias especies en Colombia, se puede observar una concentración de datos en la zona andina (en el centro del pais), por otro lado, en la zona de la Amazonia, Orinoquia y algunas zonas del pacifico se observan menos datos, que no siempre significa menos animales, si no que puede estar mostrando una falta de muestreo en esa zona.

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")
## Ventana creada exitosamente
cat("Tipo:", ventana_pais$type, "\n")
## Tipo: polygonal
cat("Rango X (lon):", paste(round(ventana_pais$xrange, 2), collapse = " – "), "\n")
## Rango X (lon): -79.03 – -66.88
cat("Rango Y (lat):", paste(round(ventana_pais$yrange, 2), collapse = " – "), "\n")
## Rango Y (lat): -4.24 – 12.43

Aquí se pueden ver las áreas definidas para este estudio y se confirma que no hay ningun error.

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: 10 / 10

Se observa que todas las especies se procesaron de manera correcta.

1.8 Paso 7 — Visualización de densidades KDE

# ── 1. Estandarizar cada KDE individualmente ────────────────────────────────
resultados_ppp_std <- lapply(resultados_ppp, function(r) {
  v <- r$kde$v
  
  v_std <- (v - min(v, na.rm = TRUE)) / 
           (max(v, na.rm = TRUE) - min(v, na.rm = TRUE))
  
  r$kde$v <- v_std
  return(r)
})

# ── 2. Configurar layout ────────────────────────────────────────────────────
n_sp  <- length(resultados_ppp_std)
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, 2, 0.8),
    oma = c(0, 0, 3, 0))

# ── 3. Graficar con misma escala (0–1) ──────────────────────────────────────
for (r in resultados_ppp_std) {
  plot(r$kde,
       main = "",
       col  = hcl.colors(64, "Inferno"),
       las  = 1,
       ribside = "right",
       zlim = c(0, 1))   # ← clave
  
  plot(ventana_pais, add = TRUE, lwd = 0.9, border = "white")
  points(r$ppp, pch = ".", cex = 1.6, col = "white")
  
  title(main = paste0(r$nombre, "\n(n=", r$n_pts, ")"),
        cex.main = 0.8)
}

par(mfrow = c(1, 1))

Para este análisis se estandarizo individualmente cada escala de datos por especie, esto con el fin de poder hacer un análisis comparativo entre estas. De esta forma, cada figura muestra las zonas de menor densidad y el hotspot máximo de cada especie, aunque no se puede hacer un análisis de la densidad entre especies si se puede ver en cuales casos se presenta solapamiento en los hotspots de estas.

Teniendo esto en cuenta se puede observar como Hay una alta concentración de especies que están concentrados en la zona andina (Leopardus pardalis, Leopardus tigrinus , Puma concolor, Odocoileus virginianus) Esto muestra un alto solapamiento entre especies carnívoras que muestran una distribución amplia dentro del país en contraste con especies como Tremarctos ornatus el cual muestra un hotspot grande pero concentrado en la zona andina central, que detodas formas se solapa con otras especies, pero que muestra que su ecosistema es mayormente el de montaña, mostrando un poco mas de especificidad a comparación de los anteriores que son un poco mas generalistas. Otras especies como Herpailurus yagouaroundi, Odocoileus virginianus ,Dicotyles tajacu también coinciden en la región andina, sin embargo muestran distribuciones en otras zonas, o sea que el solapamiento es parcial y hay mayor flexibilidad ecológica.

Hay que tener en cuenta también que aunque muchas de las especies comparten espacios, los hotspots son mas compactos en algunos casos o en otros se extienden mas, lo cual sugiere un uso distinto del espacio como también puede atribuirse a las diferencias de detectabilidad pues tienen diferentes tamaños haciendo algunos mas fáciles de captar que otros.

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
Dicotyles tajacu 86 0.940129 6.582343 1.257257
Herpailurus yagouaroundi 304 3.378974 16.618950 5.401322
Leopardus pardalis 334 3.691518 16.931664 5.753296
Leopardus tigrinus 255 2.886254 15.106428 4.893465
Leopardus wiedii 279 3.022768 13.497663 4.736430
Myrmecophaga tridactyla 167 1.812694 16.818934 1.884075
Odocoileus virginianus 139 1.495462 10.983394 1.695232
Panthera onca 338 3.601360 56.459485 3.206528
Puma concolor 178 1.933688 12.396097 2.594158
Tremarctos ornatus 205 2.315739 14.414257 4.027909

Las especies que repesentan la mayor intensidad máxima son Panthera onca (56.46) y después Leopardus pardis (16.93) , Herpailurus yagouaroundi (16.62) , Myrmecophaga tridactyla (18.82), el numero tan alto de Panthera onca se puede atribuir a la existencia de hotspots o puntos en los que los registros están concentrados en pocas áreas. En general también se puede observar que los felinos son el grupo con mayor solapamiento pues sus medias son bastante extendidas.

En cuanto a la razón, probablemente ambos factores influyen en el resultado, si bien estos animales tiene hábitats específicos en especial en selvas altamente conservadas, los registros suelen ser de áreas accesibles, parques nacionales o sitios de investigación, por lo que ese sesgo debe ser considerado al analizar los datos y entender por que puede que se vena ciertas zonas tan concentradas.

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

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

# Convertir sf → SpatVector (IMPORTANTE)
pais_vect <- terra::vect(pais_sf)

# Crear rasters binarios y enmascarar desde el inicio
rasters_bin <- map(resultados_ppp, function(r) {
  rst <- kde_a_raster_bin(r)
  rst <- terra::crop(rst, pais_vect)
  rst <- terra::mask(rst, pais_vect)
  return(rst)
})

# Sumar rasters
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
}

# 🔥 CORRECCIÓN FINAL CLAVE
riqueza <- terra::crop(riqueza, pais_vect)
riqueza <- terra::mask(riqueza, pais_vect)

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: 10 especies
library(mapview)

# Paleta DISCRETA (pocos colores, bien diferenciados)
pal_azul_disc <- c("#deebf7", "#9ecae1", "#6baed6", "#3182bd", "#08519c")

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

cortes <- seq(0, max(vals), by = 2)  # puedes ajustar (1 o 2)

mapview(
  riqueza,
  col.regions = pal_azul_disc,  # 👈 ahora sí correcto
  at = cortes,
  alpha.regions = 0.8,
  legend = TRUE,
  layer.name = paste("Riqueza —", nombre_pais)
) + 
mapview(
  pais_sf,
  color = "#4a6fa5",
  alpha.regions = 0
)

En su mayoría los hot-spots de riqueza coinciden con las áreas reconocidas de alta biodiversidad en Colombia, en donde coincide mas es en la región andina. Pero se debe resaltar que la zona del choco no esta tan representada como debería, esto puede deberse a el sesgo mencionado anteriormente, lo que indicaría que hay que aumentar los esfuerzos de muestreo en estas zonas que se conoce que son diversas.

Esto puede deberse a la limitación principal es el sesgo espacial que hay, ya que los datos tienden a ser tomados en lugares protegidos por lo que se dejan otros lugares que de pronto son de acceso mas difícil por fuera. Así que no se sabe muy bien de donde no se han tomado datos o en donde simplemente no hay presencia de la especie.

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

    sigma = 1.0 (Normal)

    sigma = 0.3

    sigma = 3.0

se puede ver que cuando se reduce el valor a 0.3 los hot-spots aumentan pero se vuelven mas pequeños, es decir que los patrones se ven mas a nivel local, mientras que cuando se aumenta a 3.0 se ve mucho mas generalizado lo que hace que se pierda el detalle local pero permite identificar tendencias espaciales mas amplias en la distribución de las especies.

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

library(rgbif);          library(spatstat)
  library(spatstat.geom);  library(spatstat.explore)
  library(terra);          library(sf)
  library(rnaturalearth);  library(rnaturalearthdata)
  library(countrycode);    library(dplyr)
  library(purrr);          library(ggplot2)
  library(viridis);        library(scales)
  library(gridExtra);      library(tidymodels)
  library(corrplot);       library(broom)
  library(palmerpenguins); library(ggpubr)
  library(emmeans);        library(rstatix)
  library(knitr);          library(kableExtra)
# ── 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,…

La base de datos contiene 9 diferentes variables pertenecientes a Chrysemys picta y su hábitat, estas incluyen: longitud de la concha (mm), sexo, peso (g), ancho concha (mm), altura concha (mm), largo cabeza (mm), temperatura del agua (Celsius), edad en años y profundidad (m). Todas las variables son cuantitativas continuas, a excepción del sexo y la edad que son cualitativas nominales. A partir de estas variables se busca hacer un análisis predictivo de la longitud de concha, determinando que variables son influyentes y cuáles no.

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)

Con este primer análisis exploratorio se puede observar diferente información sobre los datos. Con la gráfica de variable de respuesta: longitud de concha se observa como las hembras presentan una mayor frecuencia en la mayoría de las longitudes de concha y son las únicas que presentan valores cercanos a 200mm, además se observa un significativo mayor número de registros de hembras. Finalmente es importante destacar que se observa un comportamiento normal de los datos, donde la mayor frecuencia de longitudes está entre los 120 y 160 mm.

Continuando con la segunda grafica se tiene una comparación entre la longitud y el peso, adicionalmente se ve el comportamiento de cada sexo. Se observa un ajuste lineal del comportamiento de la longitud de la concha acorde al peso, donde solo algunos datos extremos no corresponden al comportamiento lineal. Se observa como la pendiente en machos es ligeramente mayor que en hembras lo que indica que a un mismo peso los machos presentan una mayor longitud de concha que las hembras, sin embargo, las hembras presentan pesos cercanos a los 3000 g que no están presentes en los machos. A pesar de esas diferencias en ambos se mantiene el comportamiento lineal entre las dos variables analizadas por lo que a mayor peso se espera mayor longitud de concha.

Continuando con la tercera gráfica, esta compara la longitud de concha con la temperatura mediante una gráfica de puntos. Los puntos se observan dispersos donde no se identifica un comportamiento claro, lo cual se ve reflejado en la línea de ajuste. Esto indica que para ninguno de los dos sexos se espera una correlación significativa entre temperatura y longitud de concha.

La última grafica presenta un complot que compara la longitud de concha de ambos sexos. En las hembras se observa una aproximación a una distribución normal, pero con desviaciones debido a los valores atípicos y colas muy largas; la mediana es un poco superior a la de los machos y se observan datos atípicos tanto mayores como menores, además se observa cómo hay una alta concentración de registros alrededor de la mediana. De este modo la longitud de concha de las hembras parece ser mayor que la de los machos en general, donde adicionalmente los machos tienen una asimetría a la izquierda donde hay un rango amplio entre el cuartil y el límite inferior.

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

La grafica anterior permite analizar la correlación que existe entre las diferentes variables utilizando una escala del -1 a 1 con colores representativos, de modo que entre más cercano el valor a 1 o -1 hay una correlación más fuerte. De este modo, la mayor correlación con la longitud de concha en orden seria: el peso (0.98), ancho de concha (0.91), altura de concha (0.86), largo de la cabeza (0.79) y la edad (0.66); Esto tiene sentido ya que al ser una variable física las condiciones corporales deberían ser bastante influyentes, especialmente las de la concha y el peso siendo una tortuga más pesada más grande y por lo tanto con un caparazón de mayor longitud. Mientras que las variables de cabeza y edad tienen una ligera menor correlación ya que la cabeza por más de que se puede ver limitada por un caparazón muy grande no es necesariamente el caso; y con la edad por más de que el caparazón tiene crecimiento indefinido este se ralentiza después de la madurez sexual lo cual puede hacer la relación no tan directa como las otras. En cuanto las variables ambientales, temperatura de agua y profundidad muestran una muy baja correlación con la longitud del caparazón, en primer lugar porque al ser organismos ectotermos por lo que si las condiciones de temperatura no son ideales buscaran alternativas para termo regularse como enterrarse en el barro para enfriarse o tomar el sol para aumentar temperatura, en segundo lugar sus hábitats son aguas someras por lo cual los cambios en profundidad son mínimos para afectar ese desarrollo del caparazón (Van dijk, 2011).

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)))
## Entrenamiento: 210 | Prueba: 70
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")
## Predictores tras preproceso: ancho_concha_mm, altura_concha_mm, largo_cabeza_mm, temp_agua_c, profundidad_m, edad_años, sexo_Macho

Se dividen los datos y se hace un preprocesamiento para poder hacer la predicción. De este modo el 75% de los datos se utilizarán como entrenamiento del modelo, mientras que los últimos 25% serán obtenidos utilizando el modelo de forma predictiva. De este modo quedan definidos los coeficientes del modelo.

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

Observamos como tras el ajuste y evaluación del modelo de entrenamiento las variables con mayor relación a la longitud de la concha se mantienen y se observa como para el ancho de la concha, altura de la concha y largo de la cabeza hay un nivel de significancia del 0.01 lo que indica que son las más influyentes en la longitud de la concha, mientras que la edad tiene una significancia del 0.01, de modo que esta última tiene menos relación que las otras variables, pero igual muestra un nivel de significancia que debe ser tomado en cuenta. Las variables temperatura del agua, profundidad y sexo del macho no muestran significancia. También es importante resaltar que el ajuste deja por fuera la variable peso (g) al estar esta muy correlacionada con otras variables como ancho y alto de la concha y por eso la detecta como redundante.

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

Posteriormente se obtuvo la bondad de ajuste para el modelo en base a los datos de entrenamiento, de esto se obtuvo que utilizando las variables entregadas el 91.74% de la variabilidad de los datos está presente; posteriormente, al compararlo con el R² ajustado se observan valores similares lo que indica que no existe un sobreajuste de los datos obtenidos y la mayoría de las variables aportan al modelo para predecir la longitud de la concha. Finalmente, vemos el error promedio de la predicción (RMSE), este arroja un valor de 6.5093 mm lo cual es un error bastante bajo teniendo en cuenta la variabilidad en la longitud de la concha (80.15mm - 205.44 mm).

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

Continuando con los datos de prueba se obtuvieron las métricas de predicción. En primer lugar, se tiene un RMSE de 6.3460 mm el cual es un error bastante bajo y similar al obtenido en el entrenamiento, lo cual es un buen indicador debido a que muestra que los errores más grandes no toman valores demasiado desajustados y al ser similar al RMSE de entrenamiento se demuestra que el modelo está captando correctamente el comportamiento de este fenómeno. Adicionalmente, el MAE (error absoluto promedio) obtenido es de 5.1580 mm también muestra un error bastante bajo en la predicción de la prueba. Pasando a la última métrica el R² muestra un valor de 89.35%, lo que representa una alta variabilidad de los datos con una disminución menor respecto al entrenamiento, lo cual no indica sobreajuste y que el modelo es capaz de predecir nuevos datos sin quedar limitado en los patrones específicos presentes en la base de datos. De esto modo, la prueba presenta un nivel de error bajo dando valores cercanos a la realidad y representa en buena medida la variabilidad de los datos, pero con una disminución menor respecto del entrenamiento, lo que indica una buena predicción de los datos.

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)

La grafica relaciona el modelo en entrenamiento y el de prueba. Presenta el RMSE y R² obtenidos con anterioridad junto a una línea punteada representando la predicción perfecta junto a los datos obtenidos mostrando el valor residual que tienen. Se observa como los residuos toman valores bajos en su mayoría presentando entre 10 a 5 mm entre los 120 y 160 mm, sin embargo en las longitudes mayores a los 160 mm se observa como los valores de la predicción se encuentran significativamente más alejados a la predicción perfecta.

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

Posteriormente se realizó una validación cruzada 10-fold donde se obtuvo el MAE, RMSE y RSQ; con esto se busca verificar el comportamiento predictivo del modelo. El MAE obtenido presenta un valor promedio de 5.2584 mm el cual es muy cercano al valor de la prueba que realizamos anteriormente lo cual muestra que el modelo está teniendo predicciones bastante correctas. Este análisis se ve reforzado por el RMSE que muestra una media de 6.5165 mm el cual es nuevamente bastante cercano al valor de nuestra prueba individual y aún más al entrenamiento, lo cual indica que el modelo está haciendo predicciones acertadas. Por otro lado, el R² presenta un valor de 0.9116 lo que indica que un 91.16% de la variabilidad de los datos ha sido representada en los 10 folds, este valor es incluso más cercano al entrenamiento que la prueba que se había realizado, lo cual muestra que en definitiva no existe un sobreajuste y que el modelo está siendo capaz de realizar una predicción precisa y variada de las longitudes de concha. Finalmente, los errores estándar bajos para cada una de las métricas muestran como en ninguno de los 10 folds realizados hubo un desajuste o un comportamiento inusual, lo cual da aún más seguridad sobre el correcto comportamiento del modelo.

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

Para terminar el análisis se hizo un diagnóstico de residuos para el modelo, de estas graficas la importancia esta sobre el Residuals vs Fitted y el Q-Q residuals. La grafica de residuals vs fitted muestra los errores distribuidos de forma aleatoria donde se observa un poco de no linealidad en el comportamiento de los errores sin embargo esta no es lo suficientemente pronunciada para determinar que la varianza cuadrada de los errores presenta un patrón y que por lo tanto el modelo este desajustado.

Con la gráfica de Q-Q residuals se observa como la mayoría de los errores presentan una distribución normal la cual solo se pierde en valores extremos inferiores y superiores que pueden asociarse a los valores atípicos encontrados en la exploración de los datos; de modo que esta grafica nuevamente confirma que el modelo está teniendo una buena capacidad predictiva con errores aceptables.

La grafica scale-location permite analizar si hay homocedasticidad donde la varianza de residuos es constante. En el grafico se observa que la línea promedio local de los puntos presenta una curvatura contraria al comportamiento constante esperado, así que se observa heterocedasticidad donde la varianza de residuos cambia en este caso siendo mas alta en valores ajustado bajos , un comportamiento similar de 120-160 y baja de 160 a 200. Este comportamiento puede causar que los valores p e intervalos de confianza no sean totalmente confiables al haber un mal cálculo del error estándar de cada coeficiente.

Terminando con la gráfica de distancia de Cook, esta permite ver la influencia de cada observación en el entrenamiento y diseño del modelo. Se observa como 3 valores muestran una influencia significativamente mas alta que el resto, sin embargo no están muy lejanos al umbral de 0.5 por lo que no se debe considerar directamente su exclusión, pero sí podrían revisarse para mejorar el modelo.

De acuerdo con lo analizado en estas graficas se puede determinar que el modelo cumple con los supuestos de linealidad donde los errores tienen una media de 0, la varianza es constante, hay una distribución normal y son independientes, esta última no hubo necesidad de validarla al no estar trabajando con varios periodos de tiempo distintos. Sin embargo, la heterocedasticidad puede ser revisada y corregida para mejorar el modelo y el comportamiento de los errores.

Se concluye que el modelo es capaz de hacer predicciones bastante precisas del comportamiento de la longitud de concha (mm) en Chrysemys picta en base a las variables de longitud de concha, ancho de concha, altura de concha, largo de la cabeza, temperatura del agua, profundidad y edad en años. Sin embargo, según lo visto en el diagnostico de residuos el modelo puede tener ciertas mejoras, en primer lugar transformaría la variable de respuesta a log para intentar reducir la diferencia entre valores grandes y pequeños, esto para poder arreglar la heterocedasticidad pasando a homocedasticidad que permite tener más confiabilidad en los valores p del modelo e intentar mejorar a su vez el ligero comportamiento no lineal observado en el diagnostico de errores.

.En segundo lugar añadir el peso como predictor podría mejorar significativamente el modelo especialmente porque se observó una fuerte correlación con el peso y la longitud de concha, que es incluso la más alta de todas; en esa misma línea se podría ver como este y las otras variables morfométricas interactúan con el sexo ya que se está asumiendo que en ambos sexos la correlación entre la longitud y las variables morfológicas son iguales, cuando esto podría no ser así. Otra correlación existente a tener en cuenta es entre las variables morfológicas las cuales pueden estar siendo redundantes entre si según lo observado en la gráfica de matriz de correlación, por lo que probar una variable que agrupe todas o algunas de ellas (ancho y altura de concha por ejemplo) podría mejorar el modelo al disminuir redundancia.

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")
## Dimensiones: 333 × 8
print(table(pinguinos$especie, pinguinos$isla))
##            
##             Biscoe Dream Torgersen
##   Adelie        44    55        47
##   Chinstrap      0    68         0
##   Gentoo       119     0         0

Primero, se extrajeron las dos variables que vamos a comparar donde se observan las especies en la columna de la izquierda y las islas en la parte superior. De aquí se observa cómo solo los pingüinos Adelie habitan las tres islas, mientras que Chinstrap y Gentoo se encuentran en una sola (Biscoe y Dream respectivamente)

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)

Con las gráficas de boxplot que se observan anteriormente se puede hacer una primera aproximación a la interacción entre la masa corporal, las especies y las islas. En primer lugar, se logra observar como la especie Gentoo tiene una mediana bastante mayor a la de las otras especies, esta se encuentra ligeramente más cercana a Q1 lo que puede indicar un ligero sesgo a la derecha, pero el violín muestra una distribución bastante simétrica. Además, se tiene un rango amplio en las masas. En Chinstrap se observa un rango más reducido y una distancia intercuartil más uniforme, donde adicionalmente el violín parece mostrar un leve sesgo hacia la derecha y esta presenta valores atípicos tanto superiores como inferiores. Finalmente, en Adelie, se presenta una mediana que se ve bastante cercana a la especie anterior, pero esta se encuentra más cercana al Q3 y el violín muestra ser más ancho en la parte inferior indicando un sesgo a la izquierda, también tiene mayor rango que Chinstrap. De este modo parece que, si existe una diferencia entre la masa corporal entre especies, especialmente entre Gentoo y las otras dos.

En cuanto a la comparación entre islas, se observa como en Biscoe se encuentra la especie Gentoo que es la que presenta mayor masa corporal de las especies, junto a la mayor mediana de los Adelie, pero el comportamiento de estos últimos parece tener un fuerte sesgo a la izquierda. Continuando con la isla Dream, en este caso los Adelie presentan su menor mediana, aunque esta vez los datos aparentan un significativo sesgo a la derecha con valores muy bajos predominando sobre los demás, adicionalmente en esta isla se observa la especie Chinstrap con una mediana cercana a los Adelie de las otras dos listas. Para la isla Torgersen solo se observa la especie Adelie, esta presenta un sesgo ligero a la derecha, esta mediana se encuentra menor a los individuos de Biscoe, pero mayor a los de Dream. Según esta exploración la relación entre masa e isla no se ve con total claridad.

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")
## 
## Prueba de Levene — H₀: varianzas iguales entre grupos
print(rstatix::levene_test(pinguinos, masa_g ~ especie))
## # A tibble: 1 × 4
##     df1   df2 statistic       p
##   <int> <int>     <dbl>   <dbl>
## 1     2   330      5.13 0.00637

Se realizo una prueba Shapiro-Wilk para determinar si los datos presentan una distribución normal. En estos se observa como Chinstrap y Gentoo presentan una distribución normal de los datos, mientras que Adelie no la presenta.

Seguido de esa se realizó una prueba de homocedasticidad de Levene para determinar si la varianza es igual entre especies. Podemos observar como el valor p es de 0.00637 el cual es menor al 0.05, por lo cual se rechaza la hipótesis nula y la varianza de cada especie es distinta, esto debe tenerse en cuenta ya que puede afectar la ANOVA y los valores p.

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

La grafica Q-Q por especie permite observar visualmente lo analizado con la prueba Shapiro-Wilk. Se observa como los Adelie en ambos extremos presentan desviaciones de la normalidad ligeras que podrían haber pasado desapercibidas sin la prueba previa. Continuando con los Chinstrap estos muestran un ajuste mayor a la normalidad excepto por los datos atípicos y unos pocos datos en la parte superior e inferior. Por último, los Gentoo presentan como sus datos centrales tienen muy buen ajuste a la normalidad y nuevamente unos pocos datos en los extremos se ven ligeramente desajustados.

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

El modelo ANOVA generado tiene en cuenta el efecto de la especie e isla en la masa y la interacción de esas dos variables sobre la masa, de modo que evalúa si la diferencia entre especies cambia según la isla.

Con la tabla podemos observar como la diferencia de especie tiene un nivel de significancia de 0.001 lo cual la hace clave que se use en el post-ANOVA. Mientras que la isla no muestra significancia estadística y por lo cual no se le realiza análisis post-ANOVA

La interacción entre isla y especie no pudo evaluarse debido a que las tres especies no están en las tres islas y por lo tanto no se puede estimar su interacción y R la elimino del modelo, esta interacción era bastante significativa ya que permitía evaluar si dependiendo de la isla la diferencia de masa entre las especies cambiaba, lo cual podía deberse a diferencia de presas, condiciones climáticas, competencia y/o otras características específicas de cada hábitat, en caso de encontrarse esta interacción podía ser evaluado que condiciones de las mencionadas eran específicamente las que causaban el efecto.

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

Con este código se hace la identificacion de los pares que tienen masas corporales distintas, esto se hizo aplicando el modelo de Tukey HSD sobre el modelo ANOVA, de esta manera se realizan comparaciones múltiples entre las combinaciones posibles de los grupos, controlando el error “family-wise error rate”.

En este caso se ignoran los factores de la isla, solo se enfoca en las especies, en el texto se observan la diferencia de medidas (diff), el intervalo de confianza al 95% (lwr, upr) y el valor ajustado (p adj).

INTERPRETACION DE COMPARACIONES:

1- Chinstrap - Adelie

Se observa que hay una diferencia de medias de aproximadamente 26.9g, tienen un intervalo de confianza de -132.8g a 186.7g, finalmente el valor p ajustado es de 0.9169.

Como el intervalo de confianza incluye el 0 y tiene un valor p mayor a 0.05, no hay evidencia de una diferencia significativa en la comparación de la masa corporal de estas dos especies.

2- Gentoo - Adelie

Se observa que hay una diferencia de medias de aproximadamente 1386.3g , tienen un intervalo de confianza de 1251.9g a 1520.7g, finalmente el valor p ajustado es 0.

Lo que indica que esta comparación de especies tiene una diferencia de masa corporal significativa.

3- Gentoo - Chinstrap

Se observa que hay una diferencia de medias de aproximadamente 1359.3g , tienen un intervalo de confianza de 1193.9g a 1524.8g, finalmente el valor p ajustado es 0.

Lo que indica que esta comparación de especies tiene una diferencia de masa corporal significativa.

Ademas estos resultados muestran que la especie Gentoo presenta una masa corporal mayor a comparación de las otras dos especies, así como la similitud entre Adelie y Chinstrap que hace que la diferencia de masa corporal no sea significativa. Dando así una precisión a las diferencias marcadas por ANOVA.

par(mar = c(5, 18, 4, 2))   # más espacio a la izquierda
plot(tukey_res,
     las = 1,
     col = "#3F51B5",
     cex.axis = 0.9,        # tamaño de texto eje
     cex.lab = 1.1)         # tamaño etiquetas
title(main = "Tukey HSD — diferencias entre especies", cex.main = 1.2)

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

Este gráfico permite visualizar las diferencias de las medidas entre las especies así como los intervalos de confianza al 95%. Las lineas horizontales representan una comparación entre dos especies, mientras que el punto central indica la diferencia estimada de sus medidas. Las barras horizontales indican los intervalos de confianza y la linea vertical que esta en el cero muestra la ausencia de diferencia entre los grupos. Así que hay que fijarse en que grupos tienen intervalos que cruzan o no esta linea del valor cero.

Es así como se puede observar que este gráfico confirma de una manera mas visual los resultados de la tabla de esta forma se facilita la interpretación ya que se puede ver mejor la significancia y magnitud de diferencias entre los grupos. Consolidando también que la especie que posee la mayor masa corporal es la de los pinguinos Gentoo a comparación de los otros dos.

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

Para estas tablas el análisis no es posible debido a que las comparaciones que se quieren hacer no son posibles por dos razones, la primera siendo que el análisis de ANOVA no encuentra significancia entre islas lo cual hace que este factor se descarte y la segunda razón es que no todas las especies de pinguinos están en las tres islas, esto genera que hayan vacíos en la información y al hacer el análisis genera celdas vacías que te traducen en lo que se esta viendo que es que todas las casillas estén nulas pues no hay información que permita hacer las comparaciones. Por lo que se recomienda escoger otros factores mas estables para poder compararlos.

Sin embargo, al hacer una búsqueda encontramos que estos ( Tukey y Bonferrini) son metodos de comparaciones múltiples que buscan controlar el error de los falsos positivos en los análisis que realizan pruebas simultáneamente, la diferencia esta en el enfoque y nivel de conservadurismo.

TUKEY: Este esta diseñado para comparar los pares de medias, es mas potente por lo que tiene una mayor capacidad de detectar diferencias reales. Por lo que detecta mas diferencias significativas.

BONFERRONI: Ajusta el nivel de significancia dividiendo el nivel de significancia entre el numero de comparaciones y es un método mas conservador lo que reduce bastante el riesgo de falsos positivos pero aumente el riesgo de no detectar diferencias reales. Lo que hace que tienda a mostrar menos diferencias significativas.

Si hubieran suficientes datos, la mejor opción hubiera sido Tukey ya que con este se tendría el análisis mas especifico en búsqueda de la mayor cantidad de diferencias, puede que Bonferroni no captara alguna de las diferencias.

library(dplyr)
library(ggplot2)
library(emmeans)
library(multcomp)
library(scales)

# Modelo (recomendado sin interacción para estas comparaciones)
modelo_anova <- aov(masa_g ~ especie, data = pinguinos)

# Obtener letras automáticamente desde Tukey
emm <- emmeans(modelo_anova, ~ especie)

letras_df <- multcomp::cld(emm, Letters = letters) %>%
  as.data.frame() %>%
  mutate(letra = trimws(.group)) %>%
  dplyr::select(especie, letra)

# Tus medias (igual que antes)
medias_g <- pinguinos %>%
  group_by(especie) %>%
  summarise(
    media = mean(masa_g),
    se = sd(masa_g)/sqrt(n()),
    .groups = "drop"
  ) %>%
  left_join(letras_df, by = "especie")

# Gráfico (igual que el tuyo)
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)

Aquí se puede observar esa misma diferencia pero ahora con la ayuda visual en la que las letras iguales son las que no tienen diferencia significativa mientras que las de letras distintas si tienen una diferencia significativa.

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")
## Kruskal-Wallis: H = 212.085 | gl = 2 | p = <2e-16
rstatix::dunn_test(pinguinos, masa_g ~ especie,
                   p.adjust.method = "bonferroni") %>%
  dplyr::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 ****

Aunque aquí se esta usando una alternativa, el análisis muestra la misma diferencia entre las masas de los pinguinos.

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)


Con esta alternativa se logra hacer un gráfico de interacción donde se muestra la relación entre la masa corporal media de las especies de los pinguinos y la isla en donde se han encontrado, las lineas representan cada especie y el eje X es correspondiente a las islas, mientras que el eje Y corresponde a la masa corporal.

Las lineas paralelas muestran que posiblemente no hay una interacción entre estos dos factores, pues se ve constante la especie y su peso en las islas. En términos visuales, no se observa evidencia clara de una interacción fuerte entre especie e isla, por que no se observan cruces evidentes ni cambios drásticos en las tendencias.

Sin embargo, esta conclusión esta basada en una estructura incompleta del diseño experimental, lo cual es de suma importancia para tener en cuenta. Por lo tanto, aunque el efecto de la especie sobre la masa corporal parece consistente, no tener suficientes datos en ciertas combinaciones limita la interpretación y análisis que se intente hacer sobre la interacción.

De este modo, se concluye que, si se observa una diferencia entre las especies donde Gentoo se diferencia de diferencia de Chinstrap y Adelie, mientras que entre Chinstrap y Adelie no se encuentra una diferencia significativa, mientras que la isla no presenta una interacción con la masa de los pingüinos. Por último, se recomienda recolectar datos en las tres islas para las tres especies para poder evaluar de mejor manera el efecto isla y especies, e isla y masa que podría ser más robusto si se hiciera para las tres especies y no solo para Adelie.

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