## =============================================================================
## ANALISIS MULTIVARIANTE DE DEPENDENCIA: WISHART, HOTELLING T^2, WILKS LAMBDA
## Y MANOVA SOBRE EL CONJUNTO DE DATOS PALMER PENGUINS
##
## Curso : Analisis Multivariante - Maestria en Estadistica
## Docente : Alfredo Valencia Toledo
## Alumnos : Percy Elbis Colque Caillahua
## Oscar Chullo Puclla
## Fecha : Septiembre 2026
##
## Este script reproduce, con datos REALES del paquete {palmerpenguins},
## todo el desarrollo teorico revisado en la Unidad II del curso:
## 1) Distribucion de Wishart (verificacion empirica de E(M) = m*Sigma)
## 2) T^2 de Hotelling (contraste de un vector de medias)
## 3) Lambda de Wilks y MANOVA (comparacion de perfiles multivariados)
## 4) Diagnostico de supuestos: normalidad multivariante (Mardia, Royston,
## Henze-Zirkler) y homogeneidad de covarianzas (Box's M)
## =============================================================================
## 0. PAQUETES
paquetes <- c("palmerpenguins", "dplyr", "MVN", "MASS", "car", "biotools",
"ggplot2", "GGally", "corrplot")
instalar_faltantes <- paquetes[!(paquetes %in% installed.packages()[, "Package"])]
if (length(instalar_faltantes) > 0) install.packages(instalar_faltantes, dependencies = TRUE)
library(palmerpenguins)
##
## Adjuntando el paquete: 'palmerpenguins'
## The following objects are masked from 'package:datasets':
##
## penguins, penguins_raw
library(dplyr)
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(MVN)
## Registered S3 method overwritten by 'lme4':
## method from
## na.action.merMod car
library(MASS)
##
## Adjuntando el paquete: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(car)
## Cargando paquete requerido: carData
##
## Adjuntando el paquete: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(ggplot2)
library(GGally)
library(corrplot)
## corrplot 0.95 loaded
suppressWarnings(suppressMessages(library(biotools))) # Box's M
set.seed(2026)
## 1. PREPARACION DE DATOS
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()
# Subconjunto para las secciones de Wishart y Hotelling: Gentoo, hembras
datos_gentoo_h <- datos_completos %>%
filter(species == "Gentoo", sex == "female") %>%
dplyr::select(all_of(vars)) %>%
as.data.frame()
cat("Dimensiones Gentoo (hembras):", dim(datos_gentoo_h), "\n")
## Dimensiones Gentoo (hembras): 58 4
summary(datos_gentoo_h)
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## Min. :40.90 Min. :13.10 Min. :203.0 Min. :3950
## 1st Qu.:43.85 1st Qu.:13.80 1st Qu.:210.0 1st Qu.:4462
## Median :45.50 Median :14.25 Median :212.0 Median :4700
## Mean :45.56 Mean :14.24 Mean :212.7 Mean :4680
## 3rd Qu.:46.88 3rd Qu.:14.60 3rd Qu.:215.0 3rd Qu.:4875
## Max. :50.50 Max. :15.50 Max. :222.0 Max. :5200
## 2. DIAGNOSTICO DE NORMALIDAD MULTIVARIANTE (Gentoo, hembras)
## Mardia, Royston y Henze-Zirkler (H-Z), + univariante Shapiro-Wilk
pM <- mvn(datos_gentoo_h, mvn_test = "mardia", univariate_test = "SW")
pR <- mvn(datos_gentoo_h, mvn_test = "royston", univariate_test = "SW")
pHZ <- mvn(datos_gentoo_h, mvn_test = "hz", univariate_test = "SW")
cat("\n================= NORMALIDAD MULTIVARIANTE =================\n")
##
## ================= NORMALIDAD MULTIVARIANTE =================
cat("\n--- Mardia (asimetria y curtosis) ---\n"); print(pM$multivariate_normality)
##
## --- Mardia (asimetria y curtosis) ---
## Test Statistic p.value Method MVN
## 1 Mardia Skewness 14.961 0.779 asymptotic ✓ Normal
## 2 Mardia Kurtosis -1.444 0.149 asymptotic ✓ Normal
cat("\n--- Royston ---\n"); print(pR$multivariate_normality)
##
## --- Royston ---
## Test Statistic p.value Method MVN
## 1 Royston 1.935 0.753 asymptotic ✓ Normal
cat("\n--- Henze-Zirkler (H-Z) ---\n"); print(pHZ$multivariate_normality)
##
## --- Henze-Zirkler (H-Z) ---
## Test Statistic p.value Method MVN
## 1 Henze-Zirkler 0.729 0.646 asymptotic ✓ Normal
cat("\n--- Shapiro-Wilk univariante ---\n"); print(pM$univariate_normality)
##
## --- Shapiro-Wilk univariante ---
## 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
# Grafico Q-Q chi-cuadrado de distancias de Mahalanobis (diagnostico visual)
xbar <- colMeans(datos_gentoo_h)
S <- cov(datos_gentoo_h)
D2 <- mahalanobis(datos_gentoo_h, center = xbar, cov = S)
n_g <- nrow(datos_gentoo_h); p <- ncol(datos_gentoo_h)
qq_teorico <- qchisq(ppoints(n_g), df = p)
# GRAFICO 1: Q-Q chi-cuadrado (Se imprime directo en la consola/visor)
plot(qq_teorico, sort(D2), pch = 19, col = "steelblue",
xlab = expression(paste("Cuantiles teoricos ", chi[4]^2)),
ylab = expression(paste("Distancias de Mahalanobis ", D[i]^2, " ordenadas")),
main = "Q-Q chi-cuadrado: normalidad multivariante (Gentoo, hembras)")
abline(0, 1, col = "firebrick", lwd = 2, lty = 2)
## 3. DISTRIBUCION DE WISHART: verificacion empirica de E(M) = m*Sigma
Sigma_hat <- cov(datos_gentoo_h) # usada como "Sigma poblacional" de referencia
m_wishart <- 40
B <- 3000
M_prom <- array(0, dim = c(p, p))
for (b in 1:B) {
Xb <- mvrnorm(n = m_wishart, mu = rep(0, p), Sigma = Sigma_hat)
M_prom <- M_prom + t(Xb) %*% Xb
}
M_prom <- M_prom / B
cat("\n================= DISTRIBUCION DE WISHART =================\n")
##
## ================= DISTRIBUCION DE WISHART =================
cat("\nE(M) simulado (promedio de", B, "replicas):\n"); print(round(M_prom, 1))
##
## E(M) simulado (promedio de 3000 replicas):
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm 167.2 19.1 65.3 6270.6
## bill_depth_mm 19.1 11.7 25.7 2287.8
## flipper_length_mm 65.3 25.7 607.8 21367.2
## body_mass_g 6270.6 2287.8 21367.2 3168106.9
cat("\nm*Sigma teorico:\n"); print(round(m_wishart * Sigma_hat, 1))
##
## m*Sigma teorico:
## 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("\nError relativo maximo:",
round(max(abs(M_prom - m_wishart * Sigma_hat) / abs(m_wishart * Sigma_hat)), 4), "\n")
##
## Error relativo maximo: 0.0128
# (n-1)*S como realizacion Wishart de los propios datos (Teorema de Fisher multivariante)
W_muestral <- (n_g - 1) * S
cat("\n(n-1)*S obtenida de los datos reales de Gentoo (realizacion W_p(Sigma, n-1)):\n")
##
## (n-1)*S obtenida de los datos reales de Gentoo (realizacion W_p(Sigma, n-1)):
print(round(W_muestral, 1))
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## bill_length_mm 239.8 27.2 94.3 8952.5
## bill_depth_mm 27.2 16.6 36.9 3224.6
## flipper_length_mm 94.3 36.9 866.0 30505.6
## body_mass_g 8952.5 3224.6 30505.6 4519321.1
## 4. T^2 DE HOTELLING (una muestra): contraste sobre el vector de medias
mu0 <- c(bill_length_mm = 45.0, bill_depth_mm = 15.5,
flipper_length_mm = 210.0, body_mass_g = 4900.0)
T2 <- as.numeric(n_g * t(xbar - mu0) %*% solve(S) %*% (xbar - mu0))
Fobs <- ((n_g - p) / (p * (n_g - 1))) * T2
pval <- pf(Fobs, df1 = p, df2 = n_g - p, lower.tail = FALSE)
Fcrit <- qf(0.95, df1 = p, df2 = n_g - p)
T2crit <- (p * (n_g - 1) / (n_g - p)) * Fcrit
cat("\n================= T^2 DE HOTELLING (una muestra) =================\n")
##
## ================= T^2 DE HOTELLING (una muestra) =================
cat("T^2 =", round(T2, 3), "\n")
## T^2 = 572.907
cat("F_obs =", round(Fobs, 3), " (df1 =", p, ", df2 =", n_g - p, ")\n")
## F_obs = 135.689 (df1 = 4 , df2 = 54 )
cat("valor p =", format.pval(pval, digits = 4), "\n")
## valor p = < 2.2e-16
cat("F_crit(0.05) =", round(Fcrit, 3), " -> T2_crit =", round(T2crit, 3), "\n")
## F_crit(0.05) = 2.543 -> T2_crit = 10.737
if (pval < 0.05) cat("Decision: se rechaza H0: mu = mu0\n") else cat("Decision: no se rechaza H0\n")
## Decision: se rechaza H0: mu = mu0
# Elipse de confianza al 95% para dos variables (bill_length vs body_mass)
idx <- c("bill_length_mm", "body_mass_g")
S2 <- S[idx, idx]; xbar2 <- xbar[idx]; mu02 <- mu0[idx]
Fcrit2 <- qf(0.95, df1 = 2, df2 = n_g - 2)
T2crit2 <- (2 * (n_g - 1) / (n_g - 2)) * Fcrit2
c2 <- T2crit2 / n_g
eig <- eigen(S2)
Ssqrt <- eig$vectors %*% diag(sqrt(eig$values)) %*% t(eig$vectors)
th <- seq(0, 2 * pi, length.out = 300)
ell <- matrix(xbar2, 2, length(th)) + sqrt(c2) * Ssqrt %*% rbind(cos(th), sin(th))
ell_df <- data.frame(bill_length_mm = ell[1, ], body_mass_g = ell[2, ])
# GRAFICO 2: Elipse de confianza
p_ellipse <- ggplot(datos_gentoo_h, aes(bill_length_mm, body_mass_g)) +
geom_point(alpha = 0.6, color = "#1f6f8b") +
geom_path(data = ell_df, aes(bill_length_mm, body_mass_g), color = "#c0392b", linewidth = 1) +
geom_point(data = data.frame(bill_length_mm = xbar2[1], body_mass_g = xbar2[2]),
aes(bill_length_mm, body_mass_g), size = 4) +
geom_point(data = data.frame(bill_length_mm = mu02[1], body_mass_g = mu02[2]),
aes(bill_length_mm, body_mass_g), size = 4, shape = 17, color = "darkorange") +
labs(title = "Region de confianza T^2 de Hotelling (95%)",
subtitle = "Punto = media muestral; triangulo = mu0 hipotetico",
x = "Longitud del pico (mm)", y = "Masa corporal (g)") +
theme_minimal()
print(p_ellipse) # Imprime explÃcitamente el gráfico en el visor
## 5. LAMBDA DE WILKS Y MANOVA: comparacion de las 3 especies
modelo <- manova(cbind(bill_length_mm, bill_depth_mm,
flipper_length_mm, body_mass_g) ~ species,
data = datos_completos)
cat("\n================= MANOVA: TODAS LAS ESPECIES =================\n")
##
## ================= MANOVA: TODAS LAS ESPECIES =================
cat("\n--- Wilks ---\n"); print(summary(modelo, test = "Wilks"))
##
## --- Wilks ---
## Df Wilks approx F num Df den Df Pr(>F)
## species 2 0.018698 516.09 8 654 < 2.2e-16 ***
## Residuals 330
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n--- Pillai ---\n"); print(summary(modelo, test = "Pillai"))
##
## --- Pillai ---
## Df Pillai approx F num Df den Df Pr(>F)
## species 2 1.6379 370.89 8 656 < 2.2e-16 ***
## Residuals 330
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n--- Hotelling-Lawley ---\n"); print(summary(modelo, test = "Hotelling-Lawley"))
##
## --- Hotelling-Lawley ---
## Df Hotelling-Lawley approx F num Df den Df Pr(>F)
## species 2 17.367 707.69 8 652 < 2.2e-16 ***
## Residuals 330
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n--- Roy ---\n"); print(summary(modelo, test = "Roy"))
##
## --- Roy ---
## Df Roy approx F num Df den Df Pr(>F)
## species 2 15.03 1232.5 4 328 < 2.2e-16 ***
## Residuals 330
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n--- ANOVAs univariados de seguimiento ---\n")
##
## --- ANOVAs univariados de seguimiento ---
print(summary.aov(modelo))
## Response bill_length_mm :
## Df Sum Sq Mean Sq F value Pr(>F)
## species 2 7015.4 3507.7 397.3 < 2.2e-16 ***
## Residuals 330 2913.5 8.8
## ---
## 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 870.79 435.39 344.83 < 2.2e-16 ***
## Residuals 330 416.67 1.26
## ---
## 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 50526 25262.9 567.41 < 2.2e-16 ***
## Residuals 330 14693 44.5
## ---
## 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 145190219 72595110 341.89 < 2.2e-16 ***
## Residuals 330 70069447 212332
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculo manual de Lambda de Wilks a partir de H y E (fines pedagogicos)
Y <- as.matrix(datos_completos[, vars])
g <- datos_completos$species
grand_mean <- colMeans(Y)
H <- matrix(0, p, p); E <- matrix(0, p, p)
for (lev in levels(g)) {
Yg_ <- Y[g == lev, ]
ng_ <- nrow(Yg_)
d_ <- matrix(colMeans(Yg_) - grand_mean, ncol = 1)
H <- H + ng_ * d_ %*% t(d_)
cen <- scale(Yg_, center = TRUE, scale = FALSE)
E <- E + t(cen) %*% cen
}
Lambda_manual <- det(E) / det(E + H)
cat("\nLambda de Wilks (calculo manual H, E):", round(Lambda_manual, 5), "\n")
##
## Lambda de Wilks (calculo manual H, E): 0.0187
# Transformacion EXACTA a F para vh = k - 1 = 2 (numero de especies - 1)
vh <- nlevels(g) - 1
ve <- nrow(Y) - nlevels(g)
F_manual <- ((ve - p + 1) / p) * ((1 - sqrt(Lambda_manual)) / sqrt(Lambda_manual))
df1_manual <- 2 * p
df2_manual <- 2 * (ve - p + 1)
p_manual <- pf(F_manual, df1_manual, df2_manual, lower.tail = FALSE)
cat("F exacta (vh=2) =", round(F_manual, 3),
" df=(", df1_manual, ",", df2_manual, ") p =", format.pval(p_manual, digits = 4), "\n")
## F exacta (vh=2) = 516.095 df=( 8 , 654 ) p = < 2.2e-16
## 6. SUPUESTO DE HOMOGENEIDAD DE COVARIANZAS: TEST DE BOX (M)
cat("\n================= TEST DE BOX (M) =================\n")
##
## ================= TEST DE BOX (M) =================
print(boxM(Y, g))
##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: Y
## Chi-Sq (approx.) = 74.731, df = 20, p-value = 3.02e-08
# GRAFICO: Matriz de dispersion por especie
p_pairs <- ggpairs(datos_completos, columns = vars, ggplot2::aes(color = species, alpha = 0.6)) +
theme_minimal()
print(p_pairs) # Imprime explÃcitamente
# GRAFICO: Comparacion H vs E por variable (diagonal), escala log
barplot(rbind(diag(H), diag(E)), beside = TRUE, log = "y",
names.arg = vars, col = c("firebrick", "steelblue"),
legend.text = c("H (entre grupos)", "E (dentro de grupos)"),
main = "Comparacion H vs E por variable (escala log)",
ylab = "Suma de cuadrados")
## Warning in yinch(0.1): y log scale: yinch() is nonsense
# GRAFICO: Correlograma de los residuos del modelo MANOVA
res <- residuals(modelo)
corrplot(cor(res), method = "ellipse", type = "upper",
title = "Correlacion entre residuos multivariados", mar = c(0, 0, 2, 0))
