Párcial 2 Bioestadística
Patrones puntuales · Regresión predictiva · Diseño de experimentos
Salomé Lema y Lis Saavedra
24 de Abril
Introducción La bioestadística es una herramienta esencial para el análisis e interpretación de fenómenos biológicos y ecológicos, ya que permite identificar patrones, establecer relaciones entre variables y evaluar hipótesis a partir de datos cuantitativos. En el presente documento se integran tres enfoques estadísticos complementarios aplicados a diferentes contextos biológicos. En primer lugar, se analizan los patrones espaciales de distribución de especies, mediante mapas de riqueza, con el fin de identificar zonas de alta concentración y discutir la influencia de factores ecológicos y del esfuerzo de muestreo sobre la intensidad espacial observada. Este tipo de análisis es especialmente útil en estudios de biodiversidad y conservación.
También se emplean herramientas de regresión lineal múltiple y análisis de varianza (ANOVA) para abordar problemas de modelamiento predictivo y comparación entre grupos. La regresión permite evaluar qué variables explican mejor la variación de un rasgo biológico, así como la capacidad predictiva del modelo mediante métricas de ajuste y error. Por su parte, el ANOVA y las pruebas post-hoc facilitan la identificación de diferencias significativas entre tratamientos o grupos experimentales, incluyendo posibles interacciones entre factores. En conjunto, estos métodos permiten una interpretación integral de los datos, destacando la importancia de validar supuestos estadísticos y contextualizar los resultados desde una perspectiva biológica.
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 — Especies y país
# ═══════════════════════════════════════════════════════════════════════
# >>> 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(
"Vultur gryphus",
"Puma concolor",
"Lycalopex culpaeus",
"Dasypus novemcinctus",
"Caracara plancus",
"Ardea alba",
"Iguana iguana",
"Boa constrictor",
"Lantana camara",
"Danaus plexippus"
)
# 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## Especies ingresadas:## 1. Vultur gryphus
## 2. Puma concolor
## 3. Lycalopex culpaeus
## 4. Dasypus novemcinctus
## 5. Caracara plancus
## 6. Ardea alba
## 7. Iguana iguana
## 8. Boa constrictor
## 9. Lantana camara
## 10. Danaus plexippus1.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(mis_especies, verificar_sp) |> list_rbind()
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)| ingresado | aceptado | rango | confianza | estado |
|---|---|---|---|---|
| Vultur gryphus | Vultur gryphus | SPECIES | 99 | OK |
| Puma concolor | Puma concolor | SPECIES | 99 | OK |
| Lycalopex culpaeus | Lycalopex culpaeus | SPECIES | 99 | OK |
| Dasypus novemcinctus | Dasypus novemcinctus | SPECIES | 99 | OK |
| Caracara plancus | Caracara plancus | SPECIES | 99 | OK |
| Ardea alba | Ardea alba | SPECIES | 97 | OK |
| Iguana iguana | Iguana iguana | SPECIES | 99 | OK |
| Boa constrictor | Boa constrictor | SPECIES | 99 | OK |
| Lantana camara | Lantana camara | SPECIES | 99 | OK |
| Danaus plexippus | Danaus plexippus | SPECIES | 99 | OK |
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"))| especie | n_registros |
|---|---|
| Ardea alba | 500 |
| Boa constrictor | 500 |
| Caracara plancus | 500 |
| Danaus plexippus | 500 |
| Dasypus novemcinctus | 500 |
| Iguana iguana | 500 |
| Lantana camara | 500 |
| Puma concolor | 500 |
| Vultur gryphus | 500 |
| Lycalopex culpaeus | 12 |
##
## Total de registros: 45121.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))
Análisis — Mapa de registros crudos: Este mapa muestra la distribución geográfica de los registros de presencia descargados desde GBIF para las 10 especies seleccionadas en Colombia. Cada color representa una especie diferente. Se observa que la mayoría de los registros se concentran en la región andina y en los departamentos del centro y nororiente del país, lo cual refleja en parte el sesgo de muestreo propio de GBIF: las áreas con mayor accesibilidad, infraestructura científica y densidad poblacional humana tienden a tener más registros documentados. Las zonas del Pacífico, Amazonía y Llanos Orientales presentan una cobertura más escasa, no necesariamente porque las especies estén ausentes, sino porque históricamente han recibido menor esfuerzo de muestreo. Este patrón es una limitación importante a considerar en el análisis posterior: los puntos no representan abundancia real, sino probabilidad de detección condicionada al esfuerzo de observación.
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## Tipo: polygonal## Rango X (lon): -79.03 – -66.88## Rango Y (lat): -4.24 – 12.431.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 / 101.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))
Análisis — Mapas de densidad KDE por especie: Los diez mapas de densidad KDE con una escala de 0–1, se evidencian dos grandes grupos de patrones espaciales. Por un lado, las especies generalistas y de amplia valencia ecológica, como Ardea alba, Caracara plancus e Iguana iguana muestran distribuciones más extendidas y con múltiples focos de densidad moderada repartidos entre la región Caribe, los valles interandinos y las tierras bajas, lo cual es coherente con su capacidad de explotar diversos hábitats. En el extremo opuesto, Vultur gryphus es la especie con el patrón espacialmente más restringido y ecológicamente más consistente de todo el conjunto: su densidad queda confinada al corredor andino sin ningún registro en tierras bajas, reflejando una distribución real y no un artefacto. Lycalopex culpaeus también aparece restringida al sur andino, pero con apenas 12 registros su superficie KDE es prácticamente no interpretable, lo que la convierte en la especie más problemática del análisis.
Por otro lado, Boa constrictor, Danaus plexippus y Lantana camara comparten un rasgo preocupante: a pesar de tener conteos altos (n=470, 461 y 415 respectivamente), presentan focos de densidad máxima extremadamente concentrados en una sola zona del país, que es el centro y norte para la boa, el corredor andino occidental para la monarca y el Eje Cafetero para la lantana, con una caída abrupta hacia la periferia. Este patrón, más que reflejar la distribución real de las especies, delata un fuerte sesgo de muestreo: muchos registros acumulados en pocas localidades accesibles o con alta actividad de naturalistas. Dasypus novemcinctus y Puma concolor muestran el problema inverso: distribuciones potencialmente amplias pero con registros escasos (n=117 y 178) que subestiman severamente su presencia real debido a que son especies crípticas y de difícil detección. En síntesis, solo Ardea alba, Caracara plancus y Vultur gryphus producen mapas en los que la señal ecológica supera claramente al ruido del esfuerzo de muestreo desigual.
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## Riqueza máxima: 10 especiesriqueza_df <- as.data.frame(riqueza, xy = TRUE) %>%
setNames(c("x", "y", "riqueza")) %>%
filter(!is.na(riqueza), riqueza > 0)
# Límites exactos de Colombia para que no se recorte el mapa
lim_col <- st_bbox(pais_sf)
ggplot() +
geom_raster(data = riqueza_df,
aes(x = x, y = y, fill = riqueza)) +
geom_sf(data = pais_sf, fill = NA, color = "white", linewidth = 0.6) +
coord_sf(
xlim = c(lim_col["xmin"] - 0.5, lim_col["xmax"] + 0.5),
ylim = c(lim_col["ymin"] - 0.5, lim_col["ymax"] + 0.5),
expand = FALSE
) +
scale_fill_gradientn(
colors = c("#0d0221", "#2d1b69", "#7b2d8b", "#c0392b",
"#e67e22", "#f1c40f", "#ffffcc"),
values = scales::rescale(c(0, 1, 2, 4, 6, 8, 10)),
name = "Riqueza\nestimada",
breaks = pretty(riqueza_df$riqueza, n = 6),
guide = guide_colorbar(barheight = 12, barwidth = 1.2,
ticks.colour = "white",
frame.colour = "white")
) +
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_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 15),
plot.subtitle = element_text(size = 10, color = "gray40"),
axis.text = element_text(size = 8),
legend.position = "right",
panel.grid = element_line(color = "gray85", linewidth = 0.3)
)
Análisis — Mapa de riqueza estimada: El mapa de riqueza integra la información espacial de todas las especies mediante álgebra de rasters: cada celda recibe un valor entre 0 y 10 según cuántas especies presentan una intensidad estimada superior a su propia mediana en esa ubicación. La paleta de colores va desde tonos oscuros (baja riqueza) hasta tonos claros amarillentos (alta riqueza), pasando por rosas y naranjas.
Los sectores con mayor riqueza acumulada tienden a coincidir con la región andina central y nororiente colombiano, donde se superponen las distribuciones de varias especies de amplio rango (como Ardea alba, Caracara plancus y Dasypus novemcinctus) junto con especies más especializadas en hábitats andinos (Vultur gryphus, Puma concolor). Las áreas periféricas de la Amazonía y el Chocó muestran menor riqueza acumulada, lo cual, como se señaló antes, puede deberse tanto a distribuciones reales como al sesgo del esfuerzo de muestreo en GBIF.
Es fundamental recordar que este mapa no es un inventario de riqueza real: es una aproximación basada en la densidad estimada de registros. Especies con más puntos disponibles en GBIF influyen más en el resultado final que especies con pocos registros, independientemente de su distribución ecológica verdadera.
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)
Análisis — Histograma de distribución de riqueza: El histograma muestra cuántas celdas del raster alcanzan cada nivel de riqueza (número de especies coocurrentes). La distribución típicamente presenta una forma sesgada hacia la derecha: la mayoría de las celdas del territorio tienen una riqueza baja o media (pocas especies coocurrentes), mientras que un número reducido de celdas concentra la máxima riqueza. Este patrón es coherente con la teoría biogeográfica, que predice que los hotspots de diversidad ocupan proporciones pequeñas del territorio total.
La paleta plasma (de morado oscuro a amarillo) refuerza visualmente el gradiente de riqueza: las barras más a la derecha, con colores más cálidos, representan las celdas de mayor diversidad estimada. Si el histograma muestra un pico muy pronunciado en valores bajos con una cola larga, indica que la riqueza está fuertemente concentrada en pocas zonas. Si la distribución es más uniforme, sugiere que las especies tienen distribuciones más homogéneas o solapadas en todo el territorio.
- ¿Qué especies presentaron mayor intensidad máxima \(\hat{\lambda}_{max}\)? ¿Podría deberse a la ecología de la especie o al esfuerzo de muestreo? Las especies que presentaron mayor intensidad máxima (λ^max) son Dasypus novemcinctus, Lantana camara y Danaus plexippus, ya que muestran los hotspots más intensos y concentrados en los mapas de densidad KDE, lo que indica una alta acumulación de registros en áreas específicas.
Este patrón puede deberse tanto a factores ecológicos, como la preferencia por ciertos hábitats o una distribución naturalmente agregada, como al esfuerzo de muestreo, ya que las zonas más accesibles o más estudiadas tienden a concentrar un mayor número de registros, generando una sobreestimación de la intensidad.
¿Los hot-spots del mapa de riqueza coinciden con áreas de alta biodiversidad conocidas en su país? Sí, los hotspots de riqueza coinciden en gran medida con regiones reconocidas por su alta biodiversidad, lo que sugiere que el patrón observado es ecológicamente coherente. Sin embargo, esta coincidencia también puede estar influenciada por un mayor esfuerzo de muestreo en estas zonas, por lo que no necesariamente refleja únicamente la riqueza real sino también la intensidad de registro.
¿Qué limitaciones tiene usar registros de GBIF para análisis de patrones de puntos? El uso de datos de GBIF presenta varias limitaciones: sesgo de muestreo espacial, ya que los registros no se distribuyen de manera uniforme; errores de georreferenciación; datos incompletos o inconsistentes; y la ausencia de información sobre esfuerzo de muestreo. Esto puede generar sobreestimaciones de densidad en zonas más accesibles o estudiadas y subestimaciones en áreas poco muestreadas.
¿Cómo cambiaría el mapa si usara
sigma = 0.3en lugar de su valor actual? ¿Y si usarasigma = 3.0?
Al modificar el valor de sigma en el KDE se altera el nivel de suavizado de la distribución: un sigma bajo produce mapas con mayor detalle y variabilidad local, resaltando agrupaciones pequeñas, mientras que un sigma alto genera una superficie más suavizada, donde se diluyen los detalles finos y se destacan patrones generales. Por tanto, la elección de sigma influye directamente en la interpretación de los patrones espaciales.
2 Módulo 2 — Regresión lineal múltiple con enfoque predictivo
2.1 Contexto biológico
Usaremos datos morfométricos de tortugas pintadas (Chrysemys picta), un modelo clásico en ecología de reptiles. El objetivo es predecir la longitud de la concha a partir de variables morfológicas y ambientales.
La diferencia entre enfoque explicativo y predictivo es central en bioestadística moderna: en el explicativo nos interesan los coeficientes e inferencia; en el predictivo nos interesa el desempeño en datos nuevos no vistos.
2.2 Datos
# ── Datos morfométricos basados en relaciones alométricas reales ──────────────
# Referencia: Janzen et al. (1992) Ecology
set.seed(42)
n <- 280
datos_tortugas <- tibble(
longitud_concha_mm = rnorm(n, mean = 138, sd = 22),
sexo = sample(c("Hembra", "Macho"), n, replace = TRUE,
prob = c(0.55, 0.45))
) %>%
mutate(
longitud_concha_mm = longitud_concha_mm + ifelse(sexo == "Hembra", 8, 0),
peso_g = 0.0018 * longitud_concha_mm^2.7 +
ifelse(sexo == "Hembra", 60, 0) + rnorm(n, 0, 35),
ancho_concha_mm = 0.82 * longitud_concha_mm + rnorm(n, 0, 8),
altura_concha_mm= 0.45 * longitud_concha_mm + rnorm(n, 0, 6),
largo_cabeza_mm = 0.28 * longitud_concha_mm + rnorm(n, 0, 5),
temp_agua_c = runif(n, 12, 28),
profundidad_m = runif(n, 0.2, 3.5),
edad_años = pmax(round(longitud_concha_mm/18 + rnorm(n, 0, 1.5)), 1)
)
glimpse(datos_tortugas)## Rows: 280
## Columns: 9
## $ longitud_concha_mm <dbl> 168.16109, 133.57664, 153.98883, 159.92298, 154.893…
## $ sexo <chr> "Macho", "Hembra", "Hembra", "Hembra", "Hembra", "H…
## $ peso_g <dbl> 1867.4922, 997.0991, 1513.8624, 1645.5372, 1564.640…
## $ ancho_concha_mm <dbl> 140.24963, 112.67478, 118.26409, 128.53102, 118.946…
## $ altura_concha_mm <dbl> 80.23671, 57.76213, 75.35438, 68.20350, 59.89479, 7…
## $ largo_cabeza_mm <dbl> 56.97275, 34.32343, 50.71398, 42.61396, 38.21520, 3…
## $ temp_agua_c <dbl> 13.69241, 24.35050, 22.17509, 21.47685, 12.18190, 1…
## $ profundidad_m <dbl> 2.9048596, 0.7727506, 0.5567349, 1.5877060, 0.42013…
## $ edad_años <dbl> 10, 6, 7, 8, 9, 8, 8, 10, 11, 9, 9, 11, 8, 9, 9, 6,…2.3 Análisis exploratorio
p1 <- ggplot(datos_tortugas, aes(x = longitud_concha_mm, fill = sexo)) +
geom_histogram(bins = 25, color = "white", alpha = 0.85) +
scale_fill_manual(values = c("#3F51B5","#E91E63")) +
labs(title = "Variable respuesta: longitud de concha",
x = "Longitud (mm)", y = "Frecuencia") + theme_minimal()
p2 <- ggplot(datos_tortugas, aes(x = peso_g, y = longitud_concha_mm,
color = sexo)) +
geom_point(alpha = 0.45, size = 1.6) +
geom_smooth(method = "lm", se = TRUE, linewidth = 1) +
scale_color_manual(values = c("#3F51B5","#E91E63")) +
labs(title = "Longitud vs Peso", x = "Peso (g)",
y = "Longitud concha (mm)") + theme_minimal()
p3 <- ggplot(datos_tortugas, aes(x = temp_agua_c, y = longitud_concha_mm,
color = sexo)) +
geom_point(alpha = 0.45, size = 1.6) +
geom_smooth(method = "loess", se = TRUE, color = "gray30", linewidth = 1) +
scale_color_manual(values = c("#3F51B5","#E91E63")) +
labs(title = "Longitud vs Temperatura", x = "Temperatura (°C)",
y = "Longitud concha (mm)") + theme_minimal()
p4 <- ggplot(datos_tortugas, aes(x = sexo, y = longitud_concha_mm,
fill = sexo)) +
geom_violin(alpha = 0.35) +
geom_boxplot(width = 0.18, alpha = 0.85) +
scale_fill_manual(values = c("#3F51B5","#E91E63"), guide = "none") +
labs(title = "Longitud por sexo", x = NULL,
y = "Longitud concha (mm)") + theme_minimal()
grid.arrange(p1, p2, p3, p4, ncol = 2)
Sobre el histograma, se puede observar que la concha tiene una distribución aproximadamente normal, ya que esta tiende a tener una forma de “campana”. La mayor concentración está entre los 120 - 160 mm aproximadamente. se nota una ligera diferencia entre sexos, donde las hembras tienden a tomar valores un poco mayores que los machos.
En cuanto a la gráfica de Longitud vs Peso, se puede observar una fuerte relación lineal positiva entre el peso y la longitud, indicando que a mayor peso, mayor longitud de la concha. También se puede observar una diferencia significativa del peso entre machos y hembras.
Sobre la gráfica inferior izquierda (Longitud vs Temperatura), se observa que no existe una relación directa entre ambas variables, ya que esta es casi inexistente. Los puntos se encuetran distribuídos aleatoriamente por la gráfica y la línea de tendencia es casi horizontal indicando su bajo relacionamiento; por lo que se puede decir que la temperatura no tiene una relación significativa con la longitud de las tortugas.
La última gráfica (Longitud por sexo) muestra que en promedio, las hembras tienen mayor longitud que los machos, aunque la diferencia no es extremadamente marcada. Además, se pueden observar algunos valores atípicos en la longitud de las hembras.
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))
Se observa que las variables peso, ancho de la concha y altura de la concha presentan valores de correlación muy cercanos a 1, lo que indica una relación positiva fuerte con relación a la longitud de la concha. En otras palabras, indica que a medida que cualquiera de estas 3 variables aumenta, la longitud de la concha también lo hará, evidenciando una relación directamente proporcional.
Por otro lado, variables como la temperatura del agua y la profundidad presentan valores muy cercanos a cero, lo cual indica que estos no tienen una relación significativa con la longitud de la concha.
Finalmente, la edad muestra una correlación central (0.66), lo cual sugiere que si influye en la longitud de la concha, aunque a menor medida que las variables morfométricas.
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: 70receta <- 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_Macho2.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)| 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 |
En la tabla se pueden observar los p-valores de distintas variables. El modelo muestra que variables morfométricas como el ancho, altura de la concha, y largo de la cabeza presentan coeficientes positivos y altamente significativos (p < 0.05), lo cual indica que están altamente relacionados con respecto al aumento de la longitud de la concha; esto se puede explicar de modo que todas forman parte del tamaño corporal del organismo, por lo que crecen de manera conjunta. La edad también resulta significativa, aunque con un efecto menor, ya que también maneja un valor p menor a 0.05 pero mayor a 0.00 (como lo tienen el ancho, altura de la concha y el largo de la cabeza), lo cual sugiere que esta si influye en el crecimiento de una forma no tan directa como las variables morfométricas mencionadas.
Por otro lado, las variables de la temperatura del agua, la profundidad y el sexo no muestran efectos estadísticamente significativos (p > 0.05), indicando que no aportan de manera relevante a la predicción de la longitud de la concha en este modelo.
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)| R² | R² ajustado | RMSE (σ) | AIC | BIC |
|---|---|---|---|---|
| 0.9174 | 0.9145 | 6.5093 | 1392.555 | 1422.678 |
El modelo presenta un alto nivel de ajuste, el cual se evidencia por un valor de R2 de 0.9174 y un R2 ajustado de 0.9145, lo cual indica que aproximadamente un 91% de la variabilidad de la concha es explicada por los predictores incluidos. La cercanía entre ambos valores sugiere que el modelo no está sobreajustado y que las variables incorporadas aportan información relevante sin redundancia excesiva.
El RMSE (6.5093) indica que el error promedio de predicción es relativamente bajo en relación con la escala de la variable respuesta, lo cual refleja una buena capacidad predictiva.
Por último tenemos los valores AIC y BIC, quienes son útiles para la comparación con otros modelos, penalizando la complejidad; en este caso, apoyan la adecuación del modelo ajustado.
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étrica | Prueba (no visto) |
|---|---|
| RMSE (mm) | 6.3460 |
| MAE (mm) | 5.1580 |
| R² | 0.8935 |
El modelo mantiene un buen desempeño en los datos de prueba, con un RMSE de 6.35 mm, un MAE de 5.16 aproximadamente, indicando que el error de predicción se mantiene bajo incluso en datos no vistos. El valor de R2 muestra que el modelo sigue explicando cerca del 89% de la variabilidad de la longitud de la concha fuera del conjunto de entrenamiento.
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 mayoría de los puntos se distribuyen cercanos a la línea diagonal de predicción perfecta, lo cual indica que el modelo presenta un buen ajuste y predice adecuadamente la longitud de la concha en los datos de prueba. La dispersión alrededor de la línea es relativamente baja, lo cual coincide con los valores de error (RMSE = 6.35 mm), evidenciando una buena precisión en las predicciones. Además, aunque existan unos puntos más alejados.
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)| Métrica | Media | Error est. | Folds |
|---|---|---|---|
| mae | 5.2584 | 0.2673 | 10 |
| rmse | 6.5165 | 0.3358 | 10 |
| rsq | 0.9116 | 0.0101 | 10 |
Esta validación cruzada muestra que el modelo tiene un desempeño consistente, con un RMSE promedio de 6.52 mm y un MAE de 5.26 mm, valores muy similares a los obtenidos en el conjunto de prueba. Además, el R2 promedio de 0.91 indica que el modelo mantiene su alta capacidad explicativa en diferentes particiones de los datos.
El bajo error estándar en todos los datos indica que los resultados son estables y no dependen fuertemente de una única división de datos. Esto indica que el modelo generaliza adecuadamente y presenta un buen equilibrio entre ajuste y capacidad predictiva.
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)
El gráfico de Residuals vs Fitted muestra como los residuos se distribuyen de forma aproximadamente aleatoria alrededor de la línea horizontal en cero, sin mostrar un patrón sistemático claro. La línea suavizada en rojo presenta una ligera tendencia decreciente, especialmente en los extremos, lo cual sugiere una leve desviación de la varianza constante. Las observaciones de 8, 123 y 188 se destacan como valores atípicos por su mayor distancia respecto a los demás.
El gráfico Q-Q Residuals evalúa el supuesto de normalidad de los residuos. Los puntos siguen de forma muy cercana a la línea diagonal teórica, indicando que los residuos se distribuyen de forma aproximadamente normal. Las desviaciones se encuentran únicamente en los extremos (8, 123 y 188) comportamiento esperado dado el tamaño muestral de 280 individuos.
El Scale-Location complementa el de Residuals vs Fitted evaluando la homocedasticidad a partir de los residuos estandarizados. Se observa una tendencia decreciente en la línea roja, lo que indica que la variabilidad de los residuos disminuye ligeramente a medida que aumentan los valores ajustados; esto sugiere una leve heterocedasticidad, aunque no lo suficientemente marcada como para invalidar el modelo. Se siguen destacando como puntos con mayor incersión los de los valores 8, 123 y 188.
Por último, el gráfico Cook’s distance identifica las observaciones con mayor influencia sobre las estimaciones del modelo. Los valores de distancia Cook para el modelo se mantienen por debajo de 0.06; las observaciones más influyentes son las de los datos 8, 123 y 105, esto debido a que son los más altos sin comprometer la estabilidad del modelo.
2.7.1 Preguntas de interpretación — Módulo 2
¿Qué predictores tienen mayor influencia sobre la longitud de la concha? ¿Es biológicamente esperable este resultado?
Los predictores con mayor influencia estadística sobre la longitud de la concha fueron el ancho de la concha, la altura de la concha, largo de la cabeza y la edad, todos con valores p signficativos (p < 0.05). Este resultado es biológicamente esperable, dado que en reptiles, el crecimiento corporal es un proceso integrado: todas las estructuras del cuerpo (concha y cabeza) aumentan de forma coordinada y proporcional a medida que el individuo se desarrolla, lo cual se refleja en las fuertes correlaciones morfométricas observadas en el análisis exploratorio. La edad también resultó significativa, aunque con un efecto de menor magnitud, consistente con el hecho de que el crecimiento en tortugas es asintótico y se desacelere con la madurez. Por su parte, la temperatura del agua, la profundidad y el sexo no resultaron significativos, lo que indica que, controladas las variables morfológicas, las condiciones ambientales puntuales no aportan poder predictivo adicional al modelo.
Compare R² en entrenamiento vs R² en prueba. ¿Hay indicios de sobreajuste?
El modelo obtuvo un R² de 0.9174 en el conjunto de entrenamiento y de 0.8935 en el conjunto de prueba, con un R² ajustado de 0.9145 en entrenamiento. La diferencia entre ambos valores es aproximadamente de 0.024, lo cual no constituye una señal de sobreajuste relevante. La cercanía entre el R² y el R² ajustado en entrenamiento indica que las variables incorporadas aportan información genuina sin redundancia excesiva. Estos resultados evidencian que el modelo captura patrones reales en los datos y generaliza adecuadamente a observaciones nuevas.
¿Por qué la RMSE en datos de prueba es más informativa que la RMSE en entrenamiento?
La RMSE en entrenamiento mide el error sobre los mismos datos con los que se construyó el modelo, por lo que tiende a subestimar el error real al no reflejar cómo se comporta el modelo frente a datos desconocidos. La RMSE en datos de prueba (6.35 mm) cuantifica el error promedio de predicción sobre observaciones que el modelo nunca vio durante su ajuste, constituyendo una estimación honesta de su capacidad de generalización. Dado que el objetivo de este módulo es predictivo, esta métrica es la que determina la utilidad práctica del model; una RMSE baja en entrenamiento acompañada de una alta en prueba revelaría sobreajuste y haría al modelo inaplicable en contextos reales. La consistencia entre ambos valores (RMSE entrenamiento = 6.51 vs prueba = 6.35) confirma que el modelo generaliza correctamente.
Analice los cuatro gráficos de diagnóstico: ¿se cumplen los supuestos del modelo lineal?
Teniendo en cuenta la información suministrada en el apartado Diagnósticos de Residuos, donde se analiza detalladamente cada gráfica, se puede concluir que si, los supuestos del modelo lineal se cumplen de manera satisfactoria.
¿Qué mejoras propondría (transformaciones, interacciones, otros predictores)?
Dado que la leve heterocedasticidad detectada en el gráfico Scale-Location podría estar relacionada con la escala de la variable respuesta, sería recomendable evaluar una transformación logarítmica de la logitud de la concha para estabilizar la varianza de los residuos.
También se puede explorar una interacción entre sexo y las variables morfométricas, dado que el dimorfismo sexual en Chrysemys picta sugiere que el efecto de estas variables sobre la longitud podría diferir entre machos y hembras.
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##
## Biscoe Dream Torgersen
## Adelie 44 55 47
## Chinstrap 0 68 0
## Gentoo 119 0 0p_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)
Análisis En la figura izquierda (masa corporal por especie), se observa una diferenciación clara entre especies. Gentoo presenta las mayores masas corporales, con una mediana cercana a 5000 g y una distribución relativamente amplia, lo que indica mayor tamaño corporal en comparación con las otras especies. Por otro lado, Adelie y Chinstrap muestran valores más bajos y bastante similares entre sí, con medianas alrededor de 3500–3800 g, aunque Chinstrap presenta ligeramente menor variabilidad. Esto sugiere, desde un enfoque exploratorio, que la especie es un factor importante en la variación de la masa corporal.
En la figura derecha (masa por isla y especie), se observa que la distribución de la masa varía no solo entre especies sino también entre islas. Por ejemplo, Gentoo (presente principalmente en Biscoe) mantiene valores altos, mientras que Adelie muestra variabilidad entre islas como Dream y Torgersen. Además, la dispersión y la presencia de algunos valores atípicos indican que la variabilidad no es completamente homogénea entre grupos. Este patrón sugiere una posible interacción entre especie e isla, donde el efecto de la especie sobre la masa corporal podría depender del contexto geográfico.
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)| 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## # A tibble: 1 × 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 2 330 5.13 0.00637ggqqplot(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()
Análisis El QQ-plot muestra la relación entre los cuantiles teóricos de una distribución normal y los cuantiles observados de la masa corporal para cada especie. En general, para Adelie y Chinstrap, los puntos siguen de manera bastante cercana la línea de referencia, lo que indica que sus distribuciones se aproximan razonablemente a la normalidad. Sin embargo, se observan ligeras desviaciones en los extremos, lo que sugiere pequeñas colas o posibles valores atípicos, aunque no lo suficientemente marcados como para invalidar el supuesto.
En el caso de Gentoo, aunque la mayoría de los puntos se alinean con la recta, se evidencia una mayor desviación en los extremos superiores, indicando una posible asimetría o presencia de valores extremos altos. Esto sugiere que la normalidad no se cumple perfectamente en esta especie, aunque el patrón general sigue siendo aceptable.
3.4 ANOVA de dos vías con interacción
## 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 ' ' 1tidy(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")| Fuente | gl | SC | CM | F | p-valor | Sig. |
|---|---|---|---|---|---|---|
| especie | 2 | 1.451902e+08 | 72595109.557 | 339.8328 | 0.0000 | *** |
| isla | 2 | 2.063697e+03 | 1031.848 | 0.0048 | 0.9952 | ns |
| Residuals | 328 | 7.006738e+07 | 213620.070 | NA | NA | ns |
| Note: | ||||||
| *** p<0.001 ** p<0.01 * p<0.05 . p<0.10 ns = no significativo |
3.5 Post-ANOVA: comparaciones múltiples
Cuando el ANOVA es significativo, las comparaciones múltiples identifican qué pares de grupos difieren.
## Tukey 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.0000000par(mar = c(5, 14, 4, 2))
plot(tukey_res, las = 1, col = "#3F51B5")
title(main = "Tukey HSD — diferencias entre especies")
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)| 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)| 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 |
Análisis El gráfico muestra los intervalos de confianza al 95% para las diferencias en masa corporal media entre cada par de especies. La línea vertical punteada en cero es clave, ya que indica que, si el intervalo de un par cruza el cero, la diferencia no es estadísticamente significativa; y que si el intervalo está completamente a un lado del cero, la diferencia sí es significativa. Se observa que los pingüinos Gentoo son significativamente más pesados que tanto Adelie como Chinstrap, mientras que estas dos últimas especies no difieren de forma significativa entre sí en términos de masa corporal. Esto es consistente con las letras de agrupación Tukey que se observan en el gráfico de barras de medias del mismo módulo.
Medias de masa corporal Este gráfico muestra la masa corporal media de cada una de las tres especies de pingüinos, acompañada de barras de error que representan el intervalo de confianza al 95%. Las letras sobre cada barra indican los grupos formados por el test de Tukey, el cual indica que letras distintas significan diferencias estadísticamente significativas entre ese par de especies (α = 0.05). Las tres especies difieren significativamente entre sí en masa corporal, cada una con una letra distinta (b, c, a). Gentoo es la especie más pesada por un margen amplio, mientras que Adelie y Chinstrap son más similares entre sí, aunque estadísticamente distinguibles. Este resultado es consistente con lo observado en el gráfico analizado anteriormente (Tukey HSD), en el cual se observa que la diferencia Gentoo−Adelie y Gentoo−Chinstrap fue grande y significativa, y Chinstrap−Adelie mostró una diferencia más pequeña pero igualmente detectable con este tamaño muestral.
medias_g <- pinguinos %>%
group_by(especie) %>%
summarise(media = mean(masa_g), se = sd(masa_g)/sqrt(n())) %>%
mutate(letra = c("b", "c", "a")) # letras Tukey: Adelie < Chinstrap < Gentoo
ggplot(medias_g, aes(x = especie, y = media, fill = especie)) +
geom_col(alpha = 0.82, width = 0.6) +
geom_errorbar(aes(ymin = media - 1.96*se, ymax = media + 1.96*se),
width = 0.18, linewidth = 0.9) +
geom_text(aes(label = letra, y = media + 1.96*se + 90),
size = 5.5, fontface = "bold") +
scale_fill_manual(values = c("#FF6B6B","#4ECDC4","#45B7D1"), guide = "none") +
scale_y_continuous(labels = comma) +
labs(title = "Medias de masa corporal ± IC 95%",
subtitle = "Letras distintas = diferencia significativa (Tukey, α = 0.05)",
x = "Especie", y = "Masa corporal media (g)") +
theme_minimal(base_size = 13)
Análisis La prueba arrojó un valor p muy cercano a cero, indicando que existe una diferencia estadísticamente significativa en la masa corporal entre al menos una de las 3 especies de pingüinos. Los resultados de Kruskal-Wallis y Dunn son completamente consistentes con los del ANOVA paramétrico, mostrando que Gentoo difiere significativamente de las otras dos especies, mientras que Adelie y Chinstrap no presentan diferencias significativas entre sí. Esta concordancia refuerza la robustez de las conclusiones, ya que dos enfoques metodológicos distintos (uno paramétrico y uno no paramétrico) llegan a las mismas conclusiones sobre la masa corporal de los pingüinos del archipiélago Palmer.
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-16rstatix::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)| Grupo 1 | Grupo 2 | Z | p | p adj. | Sig. |
|---|---|---|---|---|---|
| Adelie | Chinstrap | 0.3375 | 0.7357 | 1 | ns |
| Adelie | Gentoo | 13.6080 | 0.0000 | 0 | **** |
| Chinstrap | Gentoo | 10.7295 | 0.0000 | 0 | **** |
3.7 Gráfico de interacción especie × isla
pinguinos %>%
group_by(especie, isla) %>%
summarise(media = mean(masa_g), se = sd(masa_g)/sqrt(n()), .groups = "drop") %>%
ggplot(aes(x = isla, y = media, color = especie, group = especie)) +
geom_line(linewidth = 1.3) +
geom_point(size = 3.5) +
geom_errorbar(aes(ymin = media - se, ymax = media + se),
width = 0.12, linewidth = 0.8) +
scale_color_manual(values = c("#FF6B6B","#4ECDC4","#45B7D1"),
name = "Especie") +
scale_y_continuous(labels = comma) +
labs(title = "Gráfico de interacción: Especie × Isla",
subtitle = "Medias ± SE · Líneas paralelas = ausencia de interacción",
x = "Isla", y = "Masa corporal media (g)") +
theme_minimal(base_size = 12)
Análisis El gráfico muestra las medias de masa corporal de cada especie en cada isla. Las líneas conectan los puntos de una misma especie a través de las islas, y su comportamiento permite evaluar si existe interacción entre ambos factores. Este revela que la distribución de las especies entre islas no es uniforme, sino que refleja la segregación geográfica real del archipiélago Palmer. Gentoo se distingue claramente por su mayor masa en Biscoe, mientras que Adelie mantiene una masa consistente en todas las islas donde habita. La aparente ausencia de líneas paralelas completas no indica necesariamente una interacción biológica real, sino que es consecuencia de la distribución restringida de cada especie.
3.7.1 Preguntas de interpretación — Módulo 3
1. ¿Se cumplen los supuestos del ANOVA? ¿Cómo afecta este resultado a la interpretación? Los supuestos del ANOVA se cumplen de manera aproximada. El QQ-plot muestra que los residuos siguen una distribución cercana a la normal, aunque con ligeras desviaciones en los extremos, especialmente en la especie Gentoo. Asimismo, las gráficas exploratorias sugieren que la varianza entre grupos no es completamente homogénea, pero las diferencias no son lo suficientemente marcadas como para invalidar el análisis.
Debido a que el ANOVA es relativamente robusto a desviaciones moderadas de normalidad y homocedasticidad, estos resultados permiten continuar con la interpretación del modelo, aunque se recomienda hacerlo con cautela y, de ser necesario, contrastar con métodos no paramétricos.
2. ¿Qué fuentes de variación (especie, isla, interacción) son significativas? ¿Qué significa biológicamente la interacción? La especie es claramente una fuente de variación significativa, ya que las diferencias en masa corporal entre Gentoo y las otras especies son marcadas. La isla también contribuye a la variabilidad observada, aunque en menor medida. Además, las gráficas sugieren la presencia de una interacción entre especie e isla, evidenciada por cambios en la distribución de la masa dependiendo de la combinación de ambos factores.
Biológicamente, esta interacción implica que el efecto de la especie sobre la masa corporal no es constante en todas las islas, lo que podría reflejar diferencias en condiciones ambientales, disponibilidad de alimento o adaptación local que modulan el tamaño corporal de cada especie.
3. ¿Qué diferencia práctica hay entre Tukey HSD y Bonferroni? ¿Cuándo preferiría uno sobre el otro? La diferencia principal radica en el control del error tipo I. Tukey HSD está diseñado específicamente para comparaciones múltiples entre todos los pares de medias, controlando el error familiar de manera eficiente y manteniendo mayor poder estadístico. Por otro lado, Bonferroni es un método más conservador que ajusta el nivel de significancia dividiéndolo entre el número de comparaciones, lo que reduce la probabilidad de falsos positivos, pero también disminuye la capacidad de detectar diferencias reales.
En la práctica, Tukey HSD se prefiere cuando se realizan comparaciones pareadas entre todos los grupos, como en este caso con las especies. Bonferroni es más adecuado cuando el número de comparaciones es limitado o cuando se requiere un control más estricto del error, aunque a costa de menor sensibilidad.
4. ¿Los resultados de Kruskal-Wallis coinciden con los del ANOVA? ¿Por qué? Si, los resultados de Kruskal-Wallis coinciden con los del ANOVA debido a que ambos indican diferencias significativas en la masa corporal entre especies. Esto es a causa de que las diferencias entre grupos son marcadas, y los supuestos del ANOVA se cumplen razonablemente, permitiendo así que ambas pruebas detecten el mismo patrón. Otro dato que refuerza esta validez de los resultados es la concordancia entre los métodos paramétricos y no paramétricos.
5. ¿Qué indica el gráfico de interacción sobre la distribución de las especies entre islas? El gráfico de interacción muestra líneas aproximadamente paralelas, indicando una ausencia de una interacción fuerte entre especie e isla. Esto sugiere que las diferencias en masa corporal entre especies son consistentes en todas las islas, manteniéndose el mismo patrón relativo. Biológicamente implica que la masa corporal está principalmente determinada por la especie y no varía significativamente en función del ambiente insular.
4 Resumen ejecutivo
| 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.