DATOS: PALMER PENGUINS

Palmer Penguins

El conjunto de datos Palmer Penguins (Pingüinos del Archipiélago Palmer) es una base de datos de origen real que se ha convertido en el estándar moderno para la enseñanza de ciencia de datos, visualización y aprendizaje automático, diseñado como una alternativa excelente y éticamente superior al clásico dataset Iris de Fisher.

Origen y Recolección

Los datos fueron recolectados y puestos a disposición del público por la Dra. Kristen Gorman en colaboración con el programa de Investigación Ecológica a Largo Plazo (LTER, por sus siglas en inglés) de la Estación Palmer, en la Antártida. Las observaciones se realizaron durante un período de tres años (2007-2009).

¿Qué contiene la base de datos?

El dataset documenta las características físicas y demográficas de 344 pingüinos adultos en etapa reproductiva, pertenecientes a tres especies diferentes que habitan en el Archipiélago Palmer en la Antártida.

Los datos están segmentados por las siguientes categorías principales:

Especies observadas:

• Adélie (Pingüino de Adelia)

• Chinstrap (Pingüino Barbijo)

• Gentoo (Pingüino Papúa)

Islas de hábitat: Biscoe, Dream y Torgersen.

Variables Morfométricas (Físicas) El núcleo del análisis estadístico de esta base de datos radica en cuatro mediciones físicas fundamentales (generalmente medidas en milímetros y gramos):

• Bill Length (Longitud del pico/culmen): La medida desde la base del cráneo hasta la punta del pico superior.

• Bill Depth (Profundidad del pico/culmen): El grosor o altura del pico.

• Flipper Length (Longitud de la aleta): Medida clave que suele diferenciar drásticamente a ciertas especies (como los Gentoo, que son más grandes).

• Body Mass (Masa corporal): El peso del pingüino en gramos.

Adicionalmente, el dataset incluye el sexo (macho/hembra) y el año en que se registró la observación, lo que permite realizar análisis de dimorfismo sexual (diferencias físicas entre géneros) y estudios de clasificación.

CODE R

Subida de datos

library(palmerpenguins)
library(dplyr)
library(MVN)
library(MASS)
library(car)
library(ggplot2)
library(GGally)
library(corrplot)
suppressWarnings(suppressMessages(library(biotools)))
head(penguins)

# A tibble: 6 × 8
  species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
  <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
1 Adelie  Torgersen           39.1          18.7               181        3750
2 Adelie  Torgersen           39.5          17.4               186        3800
3 Adelie  Torgersen           40.3          18                 195        3250
4 Adelie  Torgersen           NA            NA                  NA          NA
5 Adelie  Torgersen           36.7          19.3               193        3450
6 Adelie  Torgersen           39.3          20.6               190        3650
# ℹ 2 more variables: sex <fct>, year <int>

De esta base de datos se puede observar las especies a las que pertenecen, la isla proveniente, longitud y profundidad del pico}, logitud de aleta, masa, sexo y año de medición, ademas de que la base de datos presentaba datos perdidos, pero para ejercicio de las actividades requeridas se ha hecho un filtrado de los datos completos quitando los NA.

set.seed(2026)
vars <- c("bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g")

datos_completos <- penguins %>%
  dplyr::select(species, sex, all_of(vars)) %>%
  na.omit()
summary(datos_completos$species)

   Adelie Chinstrap    Gentoo 
      146        68       119 

Dado que las especies de animales suelen ser distintas bajo diversos criterios, para el ejercicio se ha colocado 2 filtros:

1. Especies: Adelie, Chinstrap y Gentoo
2. Sexo: male / female

datos_gentoo <- datos_completos %>% filter(species == "Gentoo")
datos_gentoo_h <- datos_gentoo %>% filter(sex == "female") %>% dplyr::select(all_of(vars))
datos_gentoo_m <- datos_gentoo %>% filter(sex == "male") %>% dplyr::select(all_of(vars))

cat("Dimensiones Gentoo Hembras:", dim(datos_gentoo_h), "\n")

Dimensiones Gentoo Hembras: 58 4 

cat("Dimensiones Gentoo Machos:", dim(datos_gentoo_m), "\n")

Dimensiones Gentoo Machos: 61 4 

Normalidad

Luego de tener la separación de nuestros datos según el sexo para la especie Gentoo (Pingüino Papúa) comprobamos que los datos cumplan con los supuestos de normalidad

Macho

pM_m <- mvn(datos_gentoo_m, mvn_test = "mardia", univariate_test = "SW")
print(pM_m$multivariate_normality)

             Test Statistic p.value     Method      MVN
1 Mardia Skewness    23.145   0.282 asymptotic ✓ Normal
2 Mardia Kurtosis     0.038   0.970 asymptotic ✓ Normal

pM_m$univariate_normality

          Test          Variable Statistic p.value    Normality
1 Shapiro-Wilk    bill_length_mm     0.940   0.005 ✗ Not normal
2 Shapiro-Wilk     bill_depth_mm     0.980   0.401     ✓ Normal
3 Shapiro-Wilk flipper_length_mm     0.962   0.055     ✓ Normal
4 Shapiro-Wilk       body_mass_g     0.993   0.985     ✓ Normal

xbar_m <- colMeans(datos_gentoo_m); S_m <- cov(datos_gentoo_m)
D2_m   <- mahalanobis(datos_gentoo_m, center = xbar_m, cov = S_m)
qq_teorico_m <- qchisq(ppoints(nrow(datos_gentoo_m)), df = ncol(datos_gentoo_m))

plot(qq_teorico_m, sort(D2_m), pch = 19, col = "steelblue",
     xlab = expression(paste("Cuantiles teoricos ", chi[4]^2)),
     ylab = expression(paste("Distancias de Mahalanobis ", D[i]^2)),
     main = "Q-Q chi-cuadrado (Gentoo, machos)")
abline(0, 1, col = "firebrick", lwd = 2, lty = 2)

Hembra

pM_h <- mvn(datos_gentoo_h, mvn_test = "mardia", univariate_test = "SW")
print(pM_h$multivariate_normality)

             Test Statistic p.value     Method      MVN
1 Mardia Skewness    14.961   0.779 asymptotic ✓ Normal
2 Mardia Kurtosis    -1.444   0.149 asymptotic ✓ Normal

pM_h$univariate_normality

          Test          Variable Statistic p.value Normality
1 Shapiro-Wilk    bill_length_mm     0.989   0.895  ✓ Normal
2 Shapiro-Wilk     bill_depth_mm     0.986   0.736  ✓ Normal
3 Shapiro-Wilk flipper_length_mm     0.974   0.245  ✓ Normal
4 Shapiro-Wilk       body_mass_g     0.981   0.511  ✓ Normal

xbar_h <- colMeans(datos_gentoo_h); S_h <- cov(datos_gentoo_h)
D2_h   <- mahalanobis(datos_gentoo_h, center = xbar_h, cov = S_h)
qq_teorico_h <- qchisq(ppoints(nrow(datos_gentoo_h)), df = ncol(datos_gentoo_h))

plot(qq_teorico_h, sort(D2_h), pch = 19, col = "steelblue",
     xlab = expression(paste("Cuantiles teoricos ", chi[4]^2)),
     ylab = expression(paste("Distancias de Mahalanobis ", D[i]^2)),
     main = "Q-Q chi-cuadrado (Gentoo, hembras)")
abline(0, 1, col = "firebrick", lwd = 2, lty = 2)

A pesar de que para bill_length_mm en machos no se cumpla la normalidad univariante, el test de mardia tanto para machos y hembras, nos indica que los datos en conjunto cumplen con el supuesto de normalidad multivariante, por lo que dado el cumplimiento de la normalidad multivariante se procedió a realizar una prueba mas:

cat("\n================= HOMOGENEIDAD DE COVARIANZAS (BOX'S M - SEXO EN GENTOO) =================\n")

================= HOMOGENEIDAD DE COVARIANZAS (BOX'S M - SEXO EN GENTOO) =================

box_gentoo <- boxM(datos_gentoo[, vars], datos_gentoo$sex)
print(box_gentoo)

    Box's M-test for Homogeneity of Covariance Matrices

data:  datos_gentoo[, vars]
Chi-Sq (approx.) = 26.223, df = 10, p-value = 0.003452

La variabilidad de los datos no es igual entre los pingüinos Gentoo machos y hembras. Es decir, las características que se están midiendo están más dispersas o se relacionan de forma distinta en un sexo en comparación con el otro.

Con esta prueba se corrobora que la separación de los datos por sexo fue un procedimiento adecuado para agrupar características diferenciadas.

Distribución Wishart

m_wishart <- 40
B <- 3000
p <- length(vars)

S_h <- cov(datos_gentoo_h)
S_m <- cov(datos_gentoo_m)

macho

M_prom_m <- array(0, dim = c(p, p))
for (b in 1:B) {
  Xb_m <- mvrnorm(n = m_wishart, mu = rep(0, p), Sigma = S_m)
  M_prom_m <- M_prom_m + t(Xb_m) %*% Xb_m
}
M_prom_m <- M_prom_m / B

cat("E(M) simulado (promedio de", B, "réplicas):\n"); print(round(M_prom_m, 1))

E(M) simulado (promedio de 3000 réplicas):

                  bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
bill_length_mm             294.2          24.8             320.2     13311.7
bill_depth_mm               24.8          22.0              78.9      2392.3
flipper_length_mm          320.2          78.9            1288.3     23426.7
body_mass_g              13311.7        2392.3           23426.7   3918593.1

cat("m*Sigma teórico (Machos):\n");                   print(round(m_wishart * S_m, 1))

m*Sigma teórico (Machos):

                  bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
bill_length_mm             296.1          24.7             321.3     13335.5
bill_depth_mm               24.7          22.0              79.2      2352.8
flipper_length_mm          321.3          79.2            1287.4     23483.6
body_mass_g              13335.5        2352.8           23483.6   3922732.2

cat("Error relativo máximo (Machos):",
    round(max(abs(M_prom_m - m_wishart * S_m) / abs(m_wishart * S_m)), 4), "\n")

Error relativo máximo (Machos): 0.0168 

hembra

M_prom_h <- array(0, dim = c(p, p))
for (b in 1:B) {
  Xb_h <- mvrnorm(n = m_wishart, mu = rep(0, p), Sigma = S_h)
  M_prom_h <- M_prom_h + t(Xb_h) %*% Xb_h
}
M_prom_h <- M_prom_h / B

cat("E(M) simulado (promedio de", B, "réplicas):\n"); print(round(M_prom_h, 1))

E(M) simulado (promedio de 3000 réplicas):

                  bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
bill_length_mm             167.8          19.3              65.1      6236.0
bill_depth_mm               19.3          11.7              25.6      2291.1
flipper_length_mm           65.1          25.6             605.2     21264.8
body_mass_g               6236.0        2291.1           21264.8   3168137.7

cat("m*Sigma teórico (Hembras):\n");                  print(round(m_wishart * S_h, 1))

m*Sigma teórico (Hembras):

                  bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
bill_length_mm             168.3          19.1              66.2      6282.4
bill_depth_mm               19.1          11.7              25.9      2262.9
flipper_length_mm           66.2          25.9             607.7     21407.4
body_mass_g               6282.4        2262.9           21407.4   3171453.4

cat("Error relativo máximo (Hembras):",
    round(max(abs(M_prom_h - m_wishart * S_h) / abs(m_wishart * S_h)), 4), "\n")

Error relativo máximo (Hembras): 0.0168 

\(T^2\) de Hotelling

p <- length(vars)

#Definir vectores de medias hipoteticos (mu0) para cada sexo
mu0_h <- c(bill_length_mm = 45.0, bill_depth_mm = 14.0, flipper_length_mm = 212.0, body_mass_g = 4600.0)
mu0_m <- c(bill_length_mm = 49.5, bill_depth_mm = 15.5, flipper_length_mm = 221.0, body_mass_g = 5400.0)


n_h <- nrow(datos_gentoo_h); xbar_h <- colMeans(datos_gentoo_h); S_h <- cov(datos_gentoo_h)
n_m <- nrow(datos_gentoo_m); xbar_m <- colMeans(datos_gentoo_m); S_m <- cov(datos_gentoo_m)

Machos hotelling

T2_m   <- as.numeric(n_m * t(xbar_m - mu0_m) %*% solve(S_m) %*% (xbar_m - mu0_m))
Fobs_m <- ((n_m - p) / (p * (n_m - 1))) * T2_m
pval_m <- pf(Fobs_m, df1 = p, df2 = n_m - p, lower.tail = FALSE)
Fcrit_m  <- qf(0.95, df1 = p, df2 = n_m - p)
T2crit_m <- (p * (n_m - 1) / (n_m - p)) * Fcrit_m

cat("\n--- MACHOS GENTOO (H0: mu = mu0_m) ---\n")

--- MACHOS GENTOO (H0: mu = mu0_m) ---

cat("T^2      =", round(T2_m, 3), "\n")

T^2      = 9.943 

cat("F_obs    =", round(Fobs_m, 3), "  (df1 =", p, ", df2 =", n_m - p, ")\n")

F_obs    = 2.361   (df1 = 4 , df2 = 57 )

cat("valor p  =", format.pval(pval_m, digits = 4), "\n")

valor p  = 0.0639 

# Construccion de las elipses de confianza (95%) exactas de Hotelling
idx <- c("bill_length_mm", "body_mass_g")
th <- seq(0, 2 * pi, length.out = 300)

# Funcion auxiliar para calcular los puntos de la elipse
get_ellipse <- function(S, xbar, T2crit, n, idx) {
  S2 <- S[idx, idx]
  c2 <- T2crit / n
  eig <- eigen(S2)
  Ssqrt <- eig$vectors %*% diag(sqrt(eig$values)) %*% t(eig$vectors)
  ell <- matrix(xbar[idx], 2, length(th)) + sqrt(c2) * Ssqrt %*% rbind(cos(th), sin(th))
  data.frame(bill_length_mm = ell[1, ], body_mass_g = ell[2, ])
}

Hembras hotelling

T2_h   <- as.numeric(n_h * t(xbar_h - mu0_h) %*% solve(S_h) %*% (xbar_h - mu0_h))
Fobs_h <- ((n_h - p) / (p * (n_h - 1))) * T2_h
pval_h <- pf(Fobs_h, df1 = p, df2 = n_h - p, lower.tail = FALSE)
Fcrit_h  <- qf(0.95, df1 = p, df2 = n_h - p)
T2crit_h <- (p * (n_h - 1) / (n_h - p)) * Fcrit_h

cat("\n--- HEMBRAS GENTOO (H0: mu = mu0_h) ---\n")

--- HEMBRAS GENTOO (H0: mu = mu0_h) ---

cat("T^2      =", round(T2_h, 3), "\n")

T^2      = 12.568 

cat("F_obs    =", round(Fobs_h, 3), "  (df1 =", p, ", df2 =", n_h - p, ")\n")

F_obs    = 2.977   (df1 = 4 , df2 = 54 )

cat("valor p  =", format.pval(pval_h, digits = 4), "\n")

valor p  = 0.02709 

Gráfica

# Calculamos las elipses para hembras (usamos T2_crit con df=2 para gráfico bivariante)
Fcrit2_h <- qf(0.95, df1 = 2, df2 = n_h - 2)
T2crit2_h <- (2 * (n_h - 1) / (n_h - 2)) * Fcrit2_h
ell_df_h <- get_ellipse(S_h, xbar_h, T2crit2_h, n_h, idx)
ell_df_h$sex <- "female"

# Calculamos las elipses para machos
Fcrit2_m <- qf(0.95, df1 = 2, df2 = n_m - 2)
T2crit2_m <- (2 * (n_m - 1) / (n_m - 2)) * Fcrit2_m
ell_df_m <- get_ellipse(S_m, xbar_m, T2crit2_m, n_m, idx)
ell_df_m$sex <- "male"

# Unimos los datos para graficar
ell_df_both <- rbind(ell_df_h, ell_df_m)

# Puntos de las medias muestrales e hipoteticas
puntos_medias <- data.frame(
  bill_length_mm = c(xbar_h[1], xbar_m[1], mu0_h[1], mu0_m[1]),
  body_mass_g = c(xbar_h[4], xbar_m[4], mu0_h[4], mu0_m[4]),
  sex = c("female", "male", "female", "male"),
  tipo = c("Muestral", "Muestral", "Hipotetica (mu0)", "Hipotetica (mu0)")
)

# GRAFICO 2: Elipses de confianza de Hotelling separadas
p_ellipse <- ggplot() +
  # Puntos de dispersion originales
  geom_point(data = datos_gentoo, aes(x = bill_length_mm, y = body_mass_g, color = sex), alpha = 0.4) +
  # Contornos de las elipses exactas de Hotelling
  geom_path(data = ell_df_both, aes(x = bill_length_mm, y = body_mass_g, color = sex), linewidth = 1.2) +
  # Medias (muestral vs hipotetica)
  geom_point(data = puntos_medias, aes(x = bill_length_mm, y = body_mass_g, color = sex, shape = tipo), size = 4) +
  scale_shape_manual(values = c("Muestral" = 16, "Hipotetica (mu0)" = 17)) +
  scale_color_manual(values = c("female" = "#c0392b", "male" = "#2980b9")) +
  labs(title = "Regiones de Confianza T^2 de Hotelling (95%) por Sexo",
       subtitle = "Análisis separado: Cada sexo contrastado con su propio vector mu_0",
       x = "Longitud del pico (mm)", y = "Masa corporal (g)",
       shape = "Tipo de Media") +
  theme_minimal()

print(p_ellipse)

WILKS Y MANOVA

## -----------------------------------------------------------------------
## 1. EVALUACION DE NORMALIDAD MULTIVARIANTE Y FILTRADO DE OUTLIERS (POR GRUPO)
## -----------------------------------------------------------------------
library(palmerpenguins)
library(dplyr)
library(MVN)

vars <- c("bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g")

# Preparación de datos y creación de subgrupos (3 especies x 2 sexos = 6 grupos)
datos_completos <- penguins %>%
  dplyr::select(species, sex, all_of(vars)) %>%
  na.omit() %>%
  mutate(grupo = paste(species, sex, sep = "_"))

grupos <- unique(datos_completos$grupo)

# 1.1 EVALUACION INICIAL DE NORMALIDAD MULTIVARIANTE (Mardia) POR GRUPO
cat("\n================= 1.1 NORMALIDAD MULTIVARIANTE INICIAL (Mardia) =================\n")

================= 1.1 NORMALIDAD MULTIVARIANTE INICIAL (Mardia) =================

for (g in grupos) {
  df_sub <- datos_completos %>% filter(grupo == g) %>% dplyr::select(all_of(vars))
  mvn_res <- mvn(df_sub, mvn_test = "mardia")
  cat("\n--- Grupo:", g, "---\n")
  print(mvn_res$multivariate_normality)
}

--- Grupo: Adelie_male ---
             Test Statistic p.value     Method          MVN
1 Mardia Skewness    32.048   0.043 asymptotic ✗ Not normal
2 Mardia Kurtosis    -0.258   0.797 asymptotic     ✓ Normal

--- Grupo: Adelie_female ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    14.164   0.822 asymptotic ✓ Normal
2 Mardia Kurtosis    -0.788   0.431 asymptotic ✓ Normal

--- Grupo: Gentoo_female ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    14.961   0.779 asymptotic ✓ Normal
2 Mardia Kurtosis    -1.444   0.149 asymptotic ✓ Normal

--- Grupo: Gentoo_male ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    23.145   0.282 asymptotic ✓ Normal
2 Mardia Kurtosis     0.038   0.970 asymptotic ✓ Normal

--- Grupo: Chinstrap_female ---
             Test Statistic p.value     Method          MVN
1 Mardia Skewness    38.823   0.007 asymptotic ✗ Not normal
2 Mardia Kurtosis     2.218   0.027 asymptotic ✗ Not normal

--- Grupo: Chinstrap_male ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    21.092   0.392 asymptotic ✓ Normal
2 Mardia Kurtosis     0.040   0.968 asymptotic ✓ Normal

# 1.2 DETECCION Y REMOCION DE OUTLIERS MULTIVARIANTES CON MAHALANOBIS
p <- length(vars)
umbral_chi2 <- qchisq(0.99, df = p)  # Umbral alfa = 0.01 para chi-cuadrado con 4 grados de libertad (13.277)

datos_analizados <- datos_completos %>%
  group_by(grupo) %>%
  group_modify(~ {
    Y_g <- as.matrix(.x[, vars])
    S_g <- cov(Y_g)
    xbar_g <- colMeans(Y_g)
    D2 <- mahalanobis(Y_g, center = xbar_g, cov = S_g)
    .x %>% mutate(D2 = D2, es_outlier = D2 > umbral_chi2)
  }) %>%
  ungroup()

# Resumen de observaciones atipicas
outliers <- datos_analizados %>% filter(es_outlier)
cat("\n================= 1.2 OUTLIERS MULTIVARIANTES DETECTADOS (D^2 >", round(umbral_chi2, 3), ") =================\n")

================= 1.2 OUTLIERS MULTIVARIANTES DETECTADOS (D^2 > 13.277 ) =================

cat("Total de observaciones atipicas eliminadas:", nrow(outliers), "de", nrow(datos_completos), "\n")

Total de observaciones atipicas eliminadas: 3 de 333 

if (nrow(outliers) > 0) {
  print(outliers %>% dplyr::select(grupo, all_of(vars), D2))
}

# A tibble: 3 × 6
  grupo         bill_length_mm bill_depth_mm flipper_length_mm body_mass_g    D2
  <chr>                  <dbl>         <dbl>             <int>       <int> <dbl>
1 Adelie_male             46            21.5               194        4200  13.5
2 Chinstrap_fe…           58            17.8               181        3700  19.2
3 Gentoo_male             59.6          17                 230        6050  15.1

# Base de datos limpia
datos_limpios <- datos_analizados %>% filter(!es_outlier) %>% dplyr::select(-D2, -es_outlier)

# 1.3 RE-EVALUACION DE NORMALIDAD MULTIVARIANTE POST-LIMPIEZA
cat("\n================= 1.3 NORMALIDAD MULTIVARIANTE POST-LIMPIEZA =================\n")

================= 1.3 NORMALIDAD MULTIVARIANTE POST-LIMPIEZA =================

for (g in grupos) {
  df_sub_limpio <- datos_limpios %>% filter(grupo == g) %>% dplyr::select(all_of(vars))
  mvn_res_limpio <- mvn(df_sub_limpio, mvn_test = "mardia")
  cat("\n--- Grupo (Limpio):", g, "---\n")
  print(mvn_res_limpio$multivariate_normality)
}

--- Grupo (Limpio): Adelie_male ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    24.423   0.224 asymptotic ✓ Normal
2 Mardia Kurtosis    -0.623   0.534 asymptotic ✓ Normal

--- Grupo (Limpio): Adelie_female ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    14.164   0.822 asymptotic ✓ Normal
2 Mardia Kurtosis    -0.788   0.431 asymptotic ✓ Normal

--- Grupo (Limpio): Gentoo_female ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    14.961   0.779 asymptotic ✓ Normal
2 Mardia Kurtosis    -1.444   0.149 asymptotic ✓ Normal

--- Grupo (Limpio): Gentoo_male ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    15.457   0.750 asymptotic ✓ Normal
2 Mardia Kurtosis    -0.831   0.406 asymptotic ✓ Normal

--- Grupo (Limpio): Chinstrap_female ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    10.771   0.952 asymptotic ✓ Normal
2 Mardia Kurtosis    -0.217   0.828 asymptotic ✓ Normal

--- Grupo (Limpio): Chinstrap_male ---
             Test Statistic p.value     Method      MVN
1 Mardia Skewness    21.092   0.392 asymptotic ✓ Normal
2 Mardia Kurtosis     0.040   0.968 asymptotic ✓ Normal

## -----------------------------------------------------------------------
## 2. HOMOGENEIDAD DE COVARIANZAS (BOX'S M) SOBRE DATOS LIMPIOS
## -----------------------------------------------------------------------
# Separamos los datos limpios por sexo para prepararlos para el MANOVA
datos_limpios_h <- datos_limpios %>% filter(sex == "female")
datos_limpios_m <- datos_limpios %>% filter(sex == "male")

cat("\n================= TEST DE BOX (M): HEMBRAS =================\n")

================= TEST DE BOX (M): HEMBRAS =================

# H0: Sigma_Adelie_H = Sigma_Chinstrap_H = Sigma_Gentoo_H
box_h <- boxM(datos_limpios_h[, vars], datos_limpios_h$species)
print(box_h)

    Box's M-test for Homogeneity of Covariance Matrices

data:  datos_limpios_h[, vars]
Chi-Sq (approx.) = 45.78, df = 20, p-value = 0.0008639

cat("\n================= TEST DE BOX (M): MACHOS =================\n")

================= TEST DE BOX (M): MACHOS =================

# H0: Sigma_Adelie_M = Sigma_Chinstrap_M = Sigma_Gentoo_M
box_m <- boxM(datos_limpios_m[, vars], datos_limpios_m$species)
print(box_m)

    Box's M-test for Homogeneity of Covariance Matrices

data:  datos_limpios_m[, vars]
Chi-Sq (approx.) = 41.319, df = 20, p-value = 0.003389

## -----------------------------------------------------------------------
## 3. MANOVA: COMPARACION DE LAS 3 ESPECIES (SOBRE DATOS LIMPIOS)
## -----------------------------------------------------------------------
cat("\n================= MANOVA: DIFERENCIAS ENTRE ESPECIES (HEMBRAS LIMPIAS) =================\n")

================= MANOVA: DIFERENCIAS ENTRE ESPECIES (HEMBRAS LIMPIAS) =================

modelo_h_limpio <- manova(cbind(bill_length_mm, bill_depth_mm,
                                flipper_length_mm, body_mass_g) ~ species,
                          data = datos_limpios_h)

# Imprimimos Pillai como estadistico principal debido a Box's M
cat("\n--- Traza de Pillai (Recomendado por heterogeneidad) ---\n")

--- Traza de Pillai (Recomendado por heterogeneidad) ---

print(summary(modelo_h_limpio, test = "Pillai"))

           Df Pillai approx F num Df den Df    Pr(>F)    
species     2 1.6306   175.46      8    318 < 2.2e-16 ***
Residuals 161                                            
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n--- Wilks (Para comparacion) ---\n")

--- Wilks (Para comparacion) ---

print(summary(modelo_h_limpio, test = "Wilks"))

           Df    Wilks approx F num Df den Df    Pr(>F)    
species     2 0.015728   275.46      8    316 < 2.2e-16 ***
Residuals 161                                              
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n================= MANOVA: DIFERENCIAS ENTRE ESPECIES (MACHOS LIMPIOS) =================\n")

================= MANOVA: DIFERENCIAS ENTRE ESPECIES (MACHOS LIMPIOS) =================

modelo_m_limpio <- manova(cbind(bill_length_mm, bill_depth_mm,
                                flipper_length_mm, body_mass_g) ~ species,
                          data = datos_limpios_m)

cat("\n--- Traza de Pillai (Recomendado por heterogeneidad) ---\n")

--- Traza de Pillai (Recomendado por heterogeneidad) ---

print(summary(modelo_m_limpio, test = "Pillai"))

           Df Pillai approx F num Df den Df    Pr(>F)    
species     2 1.7279   255.59      8    322 < 2.2e-16 ***
Residuals 163                                            
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n--- Wilks (Para comparacion) ---\n")

--- Wilks (Para comparacion) ---

print(summary(modelo_m_limpio, test = "Wilks"))

           Df    Wilks approx F num Df den Df    Pr(>F)    
species     2 0.013441   305.02      8    320 < 2.2e-16 ***
Residuals 163                                              
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n================= ANOVAS UNIVARIADOS DE SEGUIMIENTO =================\n")

================= ANOVAS UNIVARIADOS DE SEGUIMIENTO =================

cat("\n--- Hembras ---\n"); print(summary.aov(modelo_h_limpio))

--- Hembras ---

 Response bill_length_mm :
             Df Sum Sq Mean Sq F value    Pr(>F)    
species       2 2968.2 1484.09  331.58 < 2.2e-16 ***
Residuals   161  720.6    4.48                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response bill_depth_mm :
             Df Sum Sq Mean Sq F value    Pr(>F)    
species       2 426.16 213.081   340.5 < 2.2e-16 ***
Residuals   161 100.75   0.626                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response flipper_length_mm :
             Df  Sum Sq Mean Sq F value    Pr(>F)    
species       2 21265.0 10632.5  418.15 < 2.2e-16 ***
Residuals   161  4093.8    25.4                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response body_mass_g :
             Df   Sum Sq  Mean Sq F value    Pr(>F)    
species       2 60354286 30177143  391.82 < 2.2e-16 ***
Residuals   161 12399994    77019                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\n--- Machos ---\n");  print(summary.aov(modelo_m_limpio))

--- Machos ---

 Response bill_length_mm :
             Df Sum Sq Mean Sq F value    Pr(>F)    
species       2 3858.0 1929.01   412.6 < 2.2e-16 ***
Residuals   163  762.1    4.68                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response bill_depth_mm :
             Df Sum Sq Mean Sq F value    Pr(>F)    
species       2 446.80 223.399  305.55 < 2.2e-16 ***
Residuals   163 119.18   0.731                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response flipper_length_mm :
             Df  Sum Sq Mean Sq F value    Pr(>F)    
species       2 28412.1 14206.0  375.28 < 2.2e-16 ***
Residuals   163  6170.2    37.9                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response body_mass_g :
             Df   Sum Sq  Mean Sq F value    Pr(>F)    
species       2 82687143 41343572  363.83 < 2.2e-16 ***
Residuals   163 18522303   113634                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El análisis MANOVA revela la existencia de diferencias morfológicas multivariadas altamente significativas entre las tres especies de pingüinos, un patrón que se mantiene constante tanto en la población de hembras como en la de machos. Utilizando la Traza de Pillai —el estadístico recomendado y más robusto frente a la heterogeneidad de covarianzas previamente detectada— confirmamos que el factor especie tiene un efecto determinante sobre el conjunto de características físicas. Tanto en las hembras (Pillai = 1.63, F = 175.46, p < 0.001) como en los machos (Pillai = 1.73, F = 255.59, p < 0.001), los valores de la traza se acercan fuertemente al máximo teórico de 2. Esto indica que la especie explica una inmensa proporción de la varianza total del modelo, logrando una separación casi perfecta de los grupos, lo cual es reafirmado por valores de Wilks extremadamente cercanos a cero en ambos sexos.

Por su parte, los ANOVAs univariados de seguimiento demuestran que esta contundente separación multivariada no es impulsada por una sola característica aislada, sino que las cuatro variables biométricas difieren significativamente entre las especies (p < 0.001 en todos los casos). No obstante, el poder de discriminación de cada rasgo varía ligeramente según el dimorfismo de cada sexo. En las hembras, la longitud de la aleta (F = 418.15) y la masa corporal (F = 391.82) son los atributos que marcan la mayor diferencia entre las especies. En contraste, en la población de machos, es la longitud del pico (F = 412.60) la característica que lidera la separación anatómica, seguida de cerca por la longitud de la aleta (F = 375.28).

library(MASS)
library(ggplot2)
library(dplyr)

# 1. Calculo de funciones discriminantes canonicas por sexo
lda_h <- lda(species ~ bill_length_mm + bill_depth_mm + flipper_length_mm + body_mass_g, data = datos_limpios_h)
lda_m <- lda(species ~ bill_length_mm + bill_depth_mm + flipper_length_mm + body_mass_g, data = datos_limpios_m)

# 2. Proyeccion de las observaciones en los ejes canonicos (LD1 y LD2)
df_cda_h <- data.frame(species = datos_limpios_h$species, sex = "Hembras", predict(lda_h)$x)
df_cda_m <- data.frame(species = datos_limpios_m$species, sex = "Machos", predict(lda_m)$x)
df_cda_total <- rbind(df_cda_h, df_cda_m)

# 3. Grafico de dispersion canonica con elipses de confianza al 95%
ggplot(df_cda_total, aes(x = LD1, y = LD2, color = species, fill = species)) +
  geom_point(alpha = 0.5, size = 2.5) +
  stat_ellipse(geom = "polygon", alpha = 0.15, level = 0.95, linewidth = 0.8) +
  facet_wrap(~ sex) +
  scale_color_manual(values = c("Adelie" = "#FF8C00", "Chinstrap" = "#A020F0", "Gentoo" = "#008B8B")) +
  scale_fill_manual(values = c("Adelie" = "#FF8C00", "Chinstrap" = "#A020F0", "Gentoo" = "#008B8B")) +
  labs(title = "Proyección de Especies en el Espacio Canónico del MANOVA",
       subtitle = "Separación clara entre centroides multivariados (Elipses de confianza al 95%)",
       x = "Primera Dimensión Canónica (LD1)",
       y = "Segunda Dimensión Canónica (LD2)",
       color = "Especie", fill = "Especie") +
  theme_minimal() +
  theme(legend.position = "bottom",
        strip.text = element_text(face = "bold", size = 11))