2026-05-29
Esta presentación desarrolla un ejemplo completo de MANOVA aplicado a un diseño experimental simulado.
El proceso se organiza en tres momentos:
El MANOVA, o análisis multivariado de varianza, es una extensión del ANOVA que permite comparar grupos o tratamientos cuando existen varias variables respuesta cuantitativas medidas simultáneamente.
En lugar de analizar una sola variable dependiente, el MANOVA analiza un vector de respuestas.
\[ \mathbf{Y}_{ij}=\begin{pmatrix} Y_{ij1}\\ Y_{ij2}\\ \vdots\\ Y_{ijp} \end{pmatrix} \]
Donde \(p\) representa el número de variables respuesta.
El MANOVA es útil cuando las variables respuesta están relacionadas entre sí.
Si se aplicaran varios ANOVA independientes, uno por cada variable respuesta, aumentaría el riesgo de cometer error Tipo I.
El MANOVA evalúa el efecto del tratamiento sobre el comportamiento conjunto de las variables.
Por esta razón, es una técnica adecuada en diseños experimentales donde el fenómeno estudiado no puede representarse mediante una sola respuesta.
| Criterio | ANOVA | MANOVA |
|---|---|---|
| Número de respuestas | Una variable dependiente | Dos o más variables dependientes |
| Pregunta central | ¿Cambian las medias entre grupos? | ¿Cambian los vectores de medias entre grupos? |
| Hipótesis | Igualdad de medias | Igualdad de vectores de medias |
| Error Tipo I | Puede aumentar si se aplican varios ANOVA | Se controla mejor al evaluar conjuntamente |
| Relación entre respuestas | No se modela | Se considera mediante covarianzas |
Se plantea un experimento agrícola con tres tratamientos de fertilización:
Se evalúan tres variables respuesta:
El objetivo es determinar si los tratamientos producen diferencias conjuntas en las respuestas agronómicas.
El diseño tiene la siguiente estructura:
\[ g=3 \]
tratamientos experimentales.
\[ n_i=8 \]
réplicas por tratamiento.
\[ p=3 \]
variables respuesta.
\[ N=24 \]
observaciones totales.
datos <- data.frame(
Tratamiento = rep(c("Control", "Biofertilizante", "Integral"), each = 8),
Replica = rep(1:8, times = 3),
Rendimiento = c(
48, 51, 49, 52, 50, 47, 53, 49,
58, 61, 60, 63, 59, 62, 64, 60,
68, 70, 67, 72, 69, 71, 73, 68
),
Altura = c(
24, 25, 23, 26, 25, 24, 27, 24,
28, 29, 30, 31, 28, 30, 31, 29,
34, 35, 33, 36, 34, 35, 37, 34
),
Biomasa = c(
9.5, 10.1, 9.8, 10.5, 10.0, 9.2, 10.8, 9.7,
12.5, 13.0, 13.2, 13.8, 12.7, 13.5, 14.0, 13.1,
16.0, 16.5, 15.8, 17.0, 16.2, 16.7, 17.4, 16.1
)
)
datos$Tratamiento <- factor(datos$Tratamiento)
datos
## Tratamiento Replica Rendimiento Altura Biomasa ## 1 Control 1 48 24 9.5 ## 2 Control 2 51 25 10.1 ## 3 Control 3 49 23 9.8 ## 4 Control 4 52 26 10.5 ## 5 Control 5 50 25 10.0 ## 6 Control 6 47 24 9.2 ## 7 Control 7 53 27 10.8 ## 8 Control 8 49 24 9.7 ## 9 Biofertilizante 1 58 28 12.5 ## 10 Biofertilizante 2 61 29 13.0 ## 11 Biofertilizante 3 60 30 13.2 ## 12 Biofertilizante 4 63 31 13.8 ## 13 Biofertilizante 5 59 28 12.7 ## 14 Biofertilizante 6 62 30 13.5 ## 15 Biofertilizante 7 64 31 14.0 ## 16 Biofertilizante 8 60 29 13.1 ## 17 Integral 1 68 34 16.0 ## 18 Integral 2 70 35 16.5 ## 19 Integral 3 67 33 15.8 ## 20 Integral 4 72 36 17.0 ## 21 Integral 5 69 34 16.2 ## 22 Integral 6 71 35 16.7 ## 23 Integral 7 73 37 17.4 ## 24 Integral 8 68 34 16.1
El modelo multivariado del diseño experimental es:
\[ \mathbf{Y}_{ij}=\boldsymbol{\mu}+\boldsymbol{\tau}_i+\boldsymbol{\varepsilon}_{ij} \]
Donde:
\[ \mathbf{Y}_{ij}=\begin{pmatrix} \text{Rendimiento}_{ij}\\ \text{Altura}_{ij}\\ \text{Biomasa}_{ij} \end{pmatrix} \]
\[ \boldsymbol{\mu} \]
es el vector de medias generales.
\[ \boldsymbol{\tau}_i \]
es el efecto multivariado del tratamiento \(i\).
La hipótesis nula establece que los tres tratamientos tienen el mismo vector de medias.
\[ H_0: \boldsymbol{\mu}_{Control}=\boldsymbol{\mu}_{Biofertilizante}=\boldsymbol{\mu}_{Integral} \]
La hipótesis alternativa establece que al menos un tratamiento tiene un vector de medias diferente.
\[ H_1: \text{al menos un vector de medias difiere entre tratamientos} \]
aggregate( cbind(Rendimiento, Altura, Biomasa) ~ Tratamiento, data = datos, FUN = mean )
## Tratamiento Rendimiento Altura Biomasa ## 1 Biofertilizante 60.875 29.50 13.2250 ## 2 Control 49.875 24.75 9.9500 ## 3 Integral 69.750 34.75 16.4625
La media general multivariada es:
colMeans(datos[, c("Rendimiento", "Altura", "Biomasa")])
## Rendimiento Altura Biomasa ## 60.16667 29.66667 13.21250
Las medias utilizadas para el desarrollo manual son:
| Tratamiento | Rendimiento | Altura | Biomasa |
|---|---|---|---|
| Control | 49.875 | 24.750 | 9.950 |
| Biofertilizante | 60.875 | 29.500 | 13.225 |
| Integral | 69.750 | 34.750 | 16.463 |
La media general es:
\[ \bar{\mathbf{Y}}=\begin{pmatrix} 60.167\\ 29.667\\ 13.213 \end{pmatrix} \]
El MANOVA descompone la variabilidad total en dos partes:
\[ \mathbf{T}=\mathbf{H}+\mathbf{E} \]
Donde:
\[ \mathbf{T} \]
es la matriz de variabilidad total.
\[ \mathbf{H} \]
es la matriz de variabilidad entre tratamientos.
\[ \mathbf{E} \]
es la matriz de variabilidad dentro de los tratamientos o matriz de error.
La matriz de error se calcula como:
\[ \mathbf{E}=\sum_{i=1}^{g}\sum_{j=1}^{n_i} (\mathbf{Y}_{ij}-\bar{\mathbf{Y}}_i) (\mathbf{Y}_{ij}-\bar{\mathbf{Y}}_i)' \]
Para este ejemplo:
\[ \mathbf{E}=\begin{pmatrix} 89.250 & 49.750 & 22.500\\ 49.750 & 33.000 & 13.025\\ 22.500 & 13.025 & 5.854 \end{pmatrix} \]
Esta matriz representa la variación interna de las unidades experimentales dentro de cada tratamiento.
La matriz de hipótesis se calcula como:
\[ \mathbf{H}=\sum_{i=1}^{g}n_i (\bar{\mathbf{Y}}_i-\bar{\mathbf{Y}}) (\bar{\mathbf{Y}}_i-\bar{\mathbf{Y}})' \]
Para este ejemplo:
\[ \mathbf{H}=\begin{pmatrix} 1586.083 & 793.583 & 517.850\\ 793.583 & 400.333 & 260.475\\ 517.850 & 260.475 & 169.653 \end{pmatrix} \]
Esta matriz representa la separación multivariada entre los tratamientos.
Y <- as.matrix(datos[, c("Rendimiento", "Altura", "Biomasa")])
grupo <- datos$Tratamiento
media_general <- colMeans(Y)
E <- matrix(0, nrow = 3, ncol = 3)
H <- matrix(0, nrow = 3, ncol = 3)
for(nivel in levels(grupo)){
Yi <- Y[grupo == nivel, ]
ni <- nrow(Yi)
media_i <- colMeans(Yi)
residuos_i <- sweep(Yi, 2, media_i, "-")
E <- E + t(residuos_i) %*% residuos_i
diferencia <- matrix(media_i - media_general, ncol = 1)
H <- H + ni * diferencia %*% t(diferencia)
}
T <- t(sweep(Y, 2, media_general, "-")) %*% sweep(Y, 2, media_general, "-")
E
## Rendimiento Altura Biomasa ## Rendimiento 89.25 49.750 22.50000 ## Altura 49.75 33.000 13.02500 ## Biomasa 22.50 13.025 5.85375
H
## [,1] [,2] [,3] ## [1,] 1586.0833 793.5833 517.8500 ## [2,] 793.5833 400.3333 260.4750 ## [3,] 517.8500 260.4750 169.6525
T
## Rendimiento Altura Biomasa ## Rendimiento 1675.3333 843.3333 540.3500 ## Altura 843.3333 433.3333 273.5000 ## Biomasa 540.3500 273.5000 175.5063
H + E
## Rendimiento Altura Biomasa ## Rendimiento 1675.3333 843.3333 540.3500 ## Altura 843.3333 433.3333 273.5000 ## Biomasa 540.3500 273.5000 175.5062
El estadístico Lambda de Wilks se define como:
\[ \Lambda=\frac{|\mathbf{E}|}{|\mathbf{E}+\mathbf{H}|} \]
Para este ejemplo:
\[ |\mathbf{E}|=64.511 \]
\[ |\mathbf{E}+\mathbf{H}|=14342.542 \]
Por tanto:
\[ \Lambda=\frac{64.511}{14342.542}=0.00450 \]
Un valor de Lambda de Wilks cercano a cero indica que la variabilidad explicada por los tratamientos es grande frente a la variabilidad interna del error.
En este ejemplo:
\[ \Lambda=0.00450 \]
Este valor sugiere que los tratamientos de fertilización generan diferencias multivariadas importantes en rendimiento, altura y biomasa.
La decisión formal se realiza con el valor p reportado por RStudio.
det_E <- det(E) det_EH <- det(E + H) lambda_wilks <- det_E / det_EH det_E
## [1] 64.51055
det_EH
## [1] 14342.54
lambda_wilks
## [1] 0.004497846
Además de Lambda de Wilks, se utilizan otros estadísticos:
\[ V=tr\left[\mathbf{H}(\mathbf{H}+\mathbf{E})^{-1}\right] \]
corresponde al criterio de Pillai.
\[ T=tr(\mathbf{E}^{-1}\mathbf{H}) \]
corresponde al criterio de Hotelling-Lawley.
\[ \theta=\lambda_{max}(\mathbf{E}^{-1}\mathbf{H}) \]
corresponde al mayor valor característico de Roy.
Los supuestos principales son:
En un diseño experimental, la independencia depende principalmente de la asignación de tratamientos y de la forma en que se recolectan las unidades experimentales.
cor(datos[, c("Rendimiento", "Altura", "Biomasa")])
## Rendimiento Altura Biomasa ## Rendimiento 1.0000000 0.9897776 0.9965022 ## Altura 0.9897776 1.0000000 0.9917449 ## Biomasa 0.9965022 0.9917449 1.0000000
El MANOVA es más apropiado cuando las variables respuesta tienen relación estadística, pero no son redundantes de manera perfecta.
library(ggplot2)
library(tidyr)
datos_largos <- pivot_longer(
datos,
cols = c(Rendimiento, Altura, Biomasa),
names_to = "Variable",
values_to = "Valor"
)
ggplot(datos_largos, aes(x = Tratamiento, y = Valor, fill = Tratamiento)) +
geom_boxplot(alpha = 0.75) +
facet_wrap(~ Variable, scales = "free_y") +
labs(
title = "Distribución de respuestas por tratamiento",
x = "Tratamiento",
y = "Valor observado"
) +
theme_minimal() +
theme(legend.position = "none")
# Instalar una sola vez si no se tiene el paquete:
# install.packages("GGally")
library(GGally)
ggpairs(
datos,
columns = c("Rendimiento", "Altura", "Biomasa"),
aes(color = Tratamiento, alpha = 0.8)
)
modelo_manova <- manova( cbind(Rendimiento, Altura, Biomasa) ~ Tratamiento, data = datos ) modelo_manova
## Call: ## manova(cbind(Rendimiento, Altura, Biomasa) ~ Tratamiento, data = datos) ## ## Terms: ## Tratamiento Residuals ## Rendimiento 1586.083 89.250 ## Altura 400.3333 33.0000 ## Biomasa 169.6525 5.8538 ## Deg. of Freedom 2 21 ## ## Residual standard errors: 2.061553 1.253566 0.5279678 ## Estimated effects may be unbalanced
summary(modelo_manova, test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F) ## Tratamiento 2 0.0044978 88.101 6 38 < 2.2e-16 *** ## Residuals 21 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
La interpretación se realiza observando el valor p.
Si:
\[ p\text{-value}<0.05 \]
se rechaza la hipótesis nula y se concluye que existen diferencias multivariadas entre tratamientos.
summary(modelo_manova, test = "Pillai")
## Df Pillai approx F num Df den Df Pr(>F) ## Tratamiento 2 1.5626 23.819 6 40 9.062e-12 *** ## Residuals 21 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El criterio de Pillai suele considerarse robusto cuando existen desviaciones moderadas de algunos supuestos.
summary(modelo_manova, test = "Hotelling-Lawley")
## Df Hotelling-Lawley approx F num Df den Df Pr(>F) ## Tratamiento 2 95.24 285.72 6 36 < 2.2e-16 *** ## Residuals 21 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Este criterio se basa en la suma de las raíces características de la matriz:
\[ \mathbf{E}^{-1}\mathbf{H} \]
summary(modelo_manova, test = "Roy")
## Df Roy approx F num Df den Df Pr(>F) ## Tratamiento 2 93.897 625.98 3 20 < 2.2e-16 *** ## Residuals 21 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El criterio de Roy se concentra en la mayor raíz característica.
Es útil cuando la separación entre grupos se explica principalmente por una dimensión multivariada dominante.
Si el MANOVA resulta significativo, se revisan los ANOVA univariados para identificar en qué variables aparecen diferencias.
summary.aov(modelo_manova)
## Response Rendimiento : ## Df Sum Sq Mean Sq F value Pr(>F) ## Tratamiento 2 1586.08 793.04 186.6 4.249e-14 *** ## Residuals 21 89.25 4.25 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Response Altura : ## Df Sum Sq Mean Sq F value Pr(>F) ## Tratamiento 2 400.33 200.167 127.38 1.81e-12 *** ## Residuals 21 33.00 1.571 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Response Biomasa : ## Df Sum Sq Mean Sq F value Pr(>F) ## Tratamiento 2 169.653 84.826 304.31 3.112e-16 *** ## Residuals 21 5.854 0.279 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Estos ANOVA no reemplazan el MANOVA, sino que ayudan a interpretar cuáles variables contribuyen a la diferencia global.
modelo_rendimiento <- aov(Rendimiento ~ Tratamiento, data = datos) modelo_altura <- aov(Altura ~ Tratamiento, data = datos) modelo_biomasa <- aov(Biomasa ~ Tratamiento, data = datos) TukeyHSD(modelo_rendimiento)
## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = Rendimiento ~ Tratamiento, data = datos) ## ## $Tratamiento ## diff lwr upr p adj ## Control-Biofertilizante -11.000 -13.598144 -8.401856 0e+00 ## Integral-Biofertilizante 8.875 6.276856 11.473144 1e-07 ## Integral-Control 19.875 17.276856 22.473144 0e+00
TukeyHSD(modelo_altura)
## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = Altura ~ Tratamiento, data = datos) ## ## $Tratamiento ## diff lwr upr p adj ## Control-Biofertilizante -4.75 -6.329851 -3.170149 6e-07 ## Integral-Biofertilizante 5.25 3.670149 6.829851 1e-07 ## Integral-Control 10.00 8.420149 11.579851 0e+00
TukeyHSD(modelo_biomasa)
## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = Biomasa ~ Tratamiento, data = datos) ## ## $Tratamiento ## diff lwr upr p adj ## Control-Biofertilizante -3.2750 -3.94039 -2.60961 0 ## Integral-Biofertilizante 3.2375 2.57211 3.90289 0 ## Integral-Control 6.5125 5.84711 7.17789 0
Estas comparaciones permiten identificar qué tratamientos difieren entre sí para cada variable respuesta.
residuos <- residuals(modelo_manova) shapiro.test(residuos[, "Rendimiento"])
## ## Shapiro-Wilk normality test ## ## data: residuos[, "Rendimiento"] ## W = 0.93545, p-value = 0.1291
shapiro.test(residuos[, "Altura"])
## ## Shapiro-Wilk normality test ## ## data: residuos[, "Altura"] ## W = 0.93552, p-value = 0.1295
shapiro.test(residuos[, "Biomasa"])
## ## Shapiro-Wilk normality test ## ## data: residuos[, "Biomasa"] ## W = 0.95494, p-value = 0.3453
La normalidad univariada de residuos es una revisión inicial.
Para el enfoque multivariado, se puede usar la distancia de Mahalanobis.
centro_res <- colMeans(residuos) cov_res <- cov(residuos) dist_maha <- mahalanobis( residuos, center = centro_res, cov = cov_res ) qqplot( qchisq(ppoints(nrow(residuos)), df = 3), dist_maha, main = "QQ plot multivariado basado en Mahalanobis", xlab = "Cuantiles teóricos Chi-cuadrado", ylab = "Distancias de Mahalanobis" ) abline(0, 1, col = "red", lwd = 2)
La homogeneidad de matrices de covarianzas se evalúa con la prueba de Box-M.
\[ H_0: \boldsymbol{\Sigma}_{Control}=\boldsymbol{\Sigma}_{Biofertilizante}=\boldsymbol{\Sigma}_{Integral} \]
# Instalar una sola vez si no se tiene el paquete:
# install.packages("biotools")
library(biotools)
boxM(
datos[, c("Rendimiento", "Altura", "Biomasa")],
datos$Tratamiento
)
by(
datos[, c("Rendimiento", "Altura", "Biomasa")],
datos$Tratamiento,
cov
)
## datos$Tratamiento: Biofertilizante ## Rendimiento Altura Biomasa ## Rendimiento 4.125000 2.214286 1.0178571 ## Altura 2.214286 1.428571 0.6000000 ## Biomasa 1.017857 0.600000 0.2678571 ## ------------------------------------------------------------ ## datos$Tratamiento: Control ## Rendimiento Altura Biomasa ## Rendimiento 4.125 2.2500000 1.0500000 ## Altura 2.250 1.6428571 0.5714286 ## Biomasa 1.050 0.5714286 0.2714286 ## ------------------------------------------------------------ ## datos$Tratamiento: Integral ## Rendimiento Altura Biomasa ## Rendimiento 4.500000 2.6428571 1.1464286 ## Altura 2.642857 1.6428571 0.6892857 ## Biomasa 1.146429 0.6892857 0.2969643
Si las matrices son similares, el supuesto de homogeneidad de covarianzas es razonable.
p <- 3 limite <- qchisq(0.975, df = p) atipicos <- data.frame( datos, Distancia_Mahalanobis = dist_maha, Atipico = dist_maha > limite ) atipicos
## Tratamiento Replica Rendimiento Altura Biomasa Distancia_Mahalanobis ## 1 Control 1 48 24 9.5 1.2895588 ## 2 Control 2 51 25 10.1 2.5913700 ## 3 Control 3 49 23 9.8 12.9722113 ## 4 Control 4 52 26 10.5 1.1939083 ## 5 Control 5 50 25 10.0 0.1466268 ## 6 Control 6 47 24 9.2 7.0950509 ## 7 Control 7 53 27 10.8 3.6843598 ## 8 Control 8 49 24 9.7 0.5024056 ## 9 Biofertilizante 1 58 28 12.5 2.1915262 ## 10 Biofertilizante 2 61 29 13.0 8.4167498 ## 11 Biofertilizante 3 60 30 13.2 6.3061150 ## 12 Biofertilizante 4 63 31 13.8 1.6162024 ## 13 Biofertilizante 5 59 28 12.7 1.8280318 ## 14 Biofertilizante 6 62 30 13.5 0.3982317 ## 15 Biofertilizante 7 64 31 14.0 2.7869048 ## 16 Biofertilizante 8 60 29 13.1 1.7658896 ## 17 Integral 1 68 34 16.0 1.3068562 ## 18 Integral 2 70 35 16.5 0.2828327 ## 19 Integral 3 67 33 15.8 2.5851864 ## 20 Integral 4 72 36 17.0 1.4490228 ## 21 Integral 5 69 34 16.2 0.9358413 ## 22 Integral 6 71 35 16.7 1.4994942 ## 23 Integral 7 73 37 17.4 4.5804221 ## 24 Integral 8 68 34 16.1 1.5752013 ## Atipico ## 1 FALSE ## 2 FALSE ## 3 TRUE ## 4 FALSE ## 5 FALSE ## 6 FALSE ## 7 FALSE ## 8 FALSE ## 9 FALSE ## 10 FALSE ## 11 FALSE ## 12 FALSE ## 13 FALSE ## 14 FALSE ## 15 FALSE ## 16 FALSE ## 17 FALSE ## 18 FALSE ## 19 FALSE ## 20 FALSE ## 21 FALSE ## 22 FALSE ## 23 FALSE ## 24 FALSE
Una observación se considera potencialmente atípica si su distancia de Mahalanobis supera el cuantil crítico de la distribución Chi-cuadrado.
# Instalar una sola vez si no se tiene el paquete:
# install.packages("heplots")
library(heplots)
heplot(modelo_manova)
El gráfico HE compara visualmente la variabilidad asociada al tratamiento con la variabilidad residual.
# Instalar una sola vez si no se tiene el paquete:
# install.packages("candisc")
library(candisc)
modelo_can <- candisc(modelo_manova)
summary(modelo_can)
plot(modelo_can)
El análisis canónico discriminante permite visualizar las combinaciones lineales de variables que separan mejor los tratamientos.
Si RStudio reporta un valor p menor que 0.05 en la prueba MANOVA, se concluye que los tratamientos tienen un efecto estadísticamente significativo sobre el vector conjunto de respuestas.
En términos del experimento:
Los tratamientos de fertilización no solo modifican una variable aislada, sino el perfil agronómico conjunto compuesto por rendimiento, altura y biomasa.
Desde el diseño experimental, el factor tratamiento explica una parte importante de la variabilidad multivariada observada.
El tratamiento integral muestra valores promedio superiores en las tres respuestas, seguido por el biofertilizante y luego por el control.
Por tanto, el MANOVA permite afirmar que la elección del tratamiento modifica de manera conjunta el desempeño productivo de las plantas.
El MANOVA permite evaluar simultáneamente varias respuestas cuantitativas en un diseño experimental.
En este ejemplo, la descomposición matricial permitió identificar la variabilidad entre tratamientos y dentro de tratamientos.
El valor de Lambda de Wilks cercano a cero sugiere una separación multivariada fuerte entre los tratamientos.
El análisis en RStudio complementa el desarrollo matemático mediante pruebas formales, visualizaciones, diagnóstico de supuestos y análisis posteriores.
Después del MANOVA se recomienda avanzar hacia:
Johnson, R. A., & Wichern, D. W. (2007). Applied Multivariate Statistical Analysis. Pearson.
Rencher, A. C., & Christensen, W. F. (2012). Methods of Multivariate Analysis. Wiley.
Hair, J. F., Black, W. C., Babin, B. J., & Anderson, R. E. (2019). Multivariate Data Analysis. Cengage.
R Core Team. (2026). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing.