3.2 MANOVA con un factor de variación
En general, cuando disponemos de una muestra, no medimos una única variable. Lo más habitual es realizar mediciones de diversas variables dependientes sobre cada unidad muestral.
Ejemplo. En el trabajo clásico de Andrews and Herzberg (1985), se
comparan manzanos sobre seis tipos de portainjertos. En concreto, se
analizan ocho manzanos sobre cada tipo de portainjerto. Para cada
manzano se determina \(\bm{y}=(y_1,\ldots,y_4)^\prime\), siendo \(y_1\) la
circunferencia del tronco (mm \(\times\) 100) a los 4 años, \(y_2\) el
crecimiento (m) a los 4 años, \(y_3\) la circunferencia del tronco (mm
\(\times\) 100) a los 15 años e \(y_4\) el peso del árbol (lb \(\times\) 1000)
a los 15 años. Los datos se encuentran en el fichero manzanos.txt.
> manzanos <- read.table("data/manzanos.txt", header = TRUE)
> head(manzanos)
portainjerto cir4 crec4 cir15 peso15
1 1 1.11 2.569 3.58 0.760
2 1 1.19 2.928 3.75 0.821
3 1 1.09 2.865 3.93 0.928
4 1 1.25 3.844 3.94 1.009
5 1 1.11 3.027 3.60 0.766
6 1 1.08 2.336 3.51 0.726Se trata de comparar las medias de \(k\) poblaciones normales multivariantes independientes y con matriz de dispersión común. Consideremos entonces \(k\) muestras independientes
\[\begin{array}{cccccc} \bm{y_{11}} & \bm{y_{12}} & \cdots & \bm{y_{1n_1}} & \mbox{de una población} & N_d(\bm{\mu_1},\bm\Sigma), \\ \bm{y_{21}} & \bm{y_{22}} & \cdots & \bm{y_{2n_2}} & \mbox{de una población} & N_d(\bm{\mu_2},\bm\Sigma), \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \bm{y_{k1}} & \bm{y_{k2}} & \cdots & \bm{y_{kn_k}} &\mbox{de una población} & N_d(\bm{\mu_k},\bm\Sigma). \end{array}\]
Cada una de las \(k\) muestras está formada por vectores independientes y con la misma distribución. Se trata, por tanto, de \(k\) muestras aleatorias simples. Además se supone que las \(k\) muestras son, entre sí, independientes. La hipótesis de homocedasticidad consiste ahora en que la matriz de covarianzas se supone que es la misma en todas las poblaciones. Recuerda que la homocedasticidad se puede contrastar mediante un test de igualdad de las matrices de covarianzas.
3.2.1 Descomposición de la variabilidad. Contraste de igualdad de medias
Igual que ocurría en el caso univariante, se puede expresar el modelo de análisis multivariante de la varianza como un caso particular del modelo lineal múltiple (la variable respuesta es un vector \(d\)-dimensional). Para cada observación se tiene,
\[\bm{y_{ij}}=\bm{\mu}+\bm{\alpha_i}+\bm{\epsilon_{ij}}=\bm{\mu_i}+\bm{\epsilon_{ij}}\] con \(i=1,\ldots,k\) y \(j=1,\ldots, n_i\). En el modelo, \(\bm{\mu_i}=\bm{\mu}+\bm{\alpha_i}\) es el vector de medias de la población \(i\)-ésima.
El objetivo del modelo de análisis de la varianza multivariante es contrastar la hipótesis nula de que los vectores de medias son iguales. Planteamos así el contraste:
\[H_0: \bm{\mu_1}=\bm{\mu_2}=\cdots=\bm{\mu_k}.\]
Es decir,
\[H_0: \left(\begin{array}{c}\mu_{11}\\ \mu_{12}\\\vdots\\ \mu_{1d}\end{array}\right)=\left(\begin{array}{c}\mu_{21}\\ \mu_{22}\\\vdots\\ \mu_{2d}\end{array}\right)=\cdots=\left(\begin{array}{c}\mu_{k1}\\ \mu_{k2}\\\vdots\\ \mu_{kd}\end{array}\right).\]
Al igual que en el caso univariante, podremos calcular la media global de la muestra y las medias muestrales por grupo. Así denotamos mediante \(\bm{\overline{y}_{i\bullet}}\) a la media muestral procedente de la población \(i\)-ésima, es decir,
\[\bm{\overline{y}_{i\bullet}}=\frac{1}{n_i}\sum_{j=1}^{n_i} \bm{y_{ij}}\] para \(i=1,\ldots,k\). Sea \(n=n_1+\ldots+n_k\) el número total de datos. Entonces, la media global de la muestra se calcula como: \[\bm{\overline{y}_{\bullet\bullet}} =\frac{1}{n}\sum_{i=1}^k \sum_{j=1}^{n_i} \bm{y_{ij}}\]
Para determinar si hay diferencias significativas entre las respuestas medias a distintos niveles del factor, el MANOVA descompone la variabilidad total de forma similar a como se hacía en el ANOVA. La diferencia está en que ahora la variabilidad total se medirá mediante la siguiente matriz \(d\times d\):
\[\sum_{i=1}^k \sum_{j=1}^{n_i} \left( \bm{y_{ij}}-\bm{\overline{y}_{\bullet\bullet}}\right) \left( \bm{y_{ij}}-\bm{\overline{y}_{\bullet\bullet}}\right)^\prime{}.\] Por analogía con el caso univariante, en el que se comparaban la variabilidad inter-grupos con la variabilidad intra-grupos, se calculan ahora las matrices
\[\bm{H}=\sum_{i=1}^kn_i(\bm{\overline{y}_{i\bullet}}-\bm{\overline{y}_{\bullet\bullet}})(\bm{\overline{y}_{i\bullet}}-\bm{\overline{y}_{\bullet\bullet}})^\prime{}\] y
\[\bm{E} = \sum_{i=1}^k \sum_{j=1}^{n_i} \left( \bm{y_{ij}}- \bm{\overline{y}_{i\bullet}}\right) \left( \bm{y_{ij}}- \bm{\overline{y}_{i\bullet}}\right)^\prime{}.\]
Se muestra a continuación la tabla del análisis multivariante de la varianza o tabla MANOVA:
| Fuente de variación | Matriz de covarianzas | Grados de libertad |
|---|---|---|
| Entre poblaciones | \({ \bm{H}=\sum_{i=1}^kn_i(\bm{\overline{y}_{i\bullet}}-\bm{\overline{y}_{\bullet\bullet}})(\bm{\overline{y}_{i\bullet}}-\bm{\overline{y}_{\bullet\bullet}})^\prime{}}\) | \(k-1\) |
| Error | \({\bm{E} = \sum_{i=1}^k \sum_{j=1}^{n_i}\left( \bm{y_{ij}}- \bm{\overline{y}_{i\bullet}}\right)\left( \bm{y_{ij}}- \bm{\overline{y}_{i\bullet}}\right)^\prime{}}\) | \({\sum_{i=1}^k \left(n_i-1\right)}\) |
| Total | \({\bm{E_H} = \sum_{i=1}^k \sum_{j=1}^{n_i}\left( \bm{y_{ij}}-\bm{\overline{y}_{\bullet\bullet}}\right)\left( \bm{y_{ij}}-\bm{\overline{y}_{\bullet\bullet}}\right)^\prime{}}\) | \({\sum_{i=1}^k n_i-1}\) |
Observamos que esta tabla no es más que una extensión de la tabla ANOVA al caso multivariante. El sentido común nos invita a rechazar la hipótesis nula cuando la variabilidad proveniente de las diferencias entre poblaciones (que medimos mediante la matriz \(\bm{H}\)) sea grande comparada con la proveniente del error (medida por \(\bm{E}\)).
3.2.2 El estadístico \(\Lambda\) de Wilks
El contraste de razón de verosimiltudes para el contraste \(H_0: \bm{\mu_1}=\bm{\mu_2}=\cdots=\bm{\mu_k}\) viene dado por
\[\Lambda=\frac{|\bm{E}|}{|\bm{E}+\bm{H}|}.\]
Este estadístico, conocido como estadístico de Wilks, compara la matriz \(\bm{E}\), que mide la variabilidad intra-grupos, con la matriz \(\bm{E_H}=\bm{E}+\bm{H}\), que mide la variabilidad total8. Se puede demostrar que, bajo la hipótesis nula,
\[\begin{aligned} \bm{H} &\in & \mbox{Wishart}_d(\bm\Sigma,k-1) \\ \bm{E} &\in & \mbox{Wishart}_d(\bm\Sigma,n-k)\end{aligned}\] y además son independientes. Entonces tenemos que9
\[\frac{|\bm{E}|}{|\bm{E}+\bm{H}|}\in \Lambda(d,k-1,n-k).\] Por tanto, rechazaremos la hipótesis nula cuando el estadístico \(|\bm{E}|/|\bm{E}+\bm{H}|\) tome un valor menor que el cuantil \(\alpha\) de la distribución \(\Lambda(d,k-1,n-k)\), siendo \(\alpha\) el nivel de significación fijado de antemano. El estadístico se puede aproximar por una distribución \(F\) de Snédecor.
Ejemplo. Volvemos al ejemplo sobre manzanos. Se trata de contrastar \(H_0: \bm{\mu_1}=\bm{\mu_2}=\cdots=\bm{\mu_6}\), donde para cada \(i\) el vector \(\bm{\mu_i}=(\mu_{i1},\ldots,\mu{i4})^\prime\) representa el vector de medias de \(\bm{y}=(y_1,\ldots,y_4)^\prime{}\) para el tipo \(i\) de portainjerto. Los datos son:
> x <- factor(manzanos[, 1])
> y <- as.matrix(manzanos[, 2:5])
> n <- dim(y)[1]
> d <- dim(y)[2]
> n
[1] 48Calculamos a continuación la media global de la muestra \(\bm{\overline{y}_{\bullet\bullet}}\) y las medias \(\bm{\overline{y}_{i\bullet}}\) correspondientes a cada grupo:
> medyi <- t(apply(y, 2, tapply, x, mean))
> medyi
1 2 3 4 5 6
cir4 1.137500 1.157500 1.107500 1.09750 1.08000 1.036250
crec4 2.977125 3.109125 2.815250 2.87975 2.55725 2.214625
cir15 3.738750 4.515000 4.455000 3.90625 4.31250 3.596250
peso15 0.871125 1.280500 1.391375 1.03900 1.18100 0.735000A continuación, realizamos el MANOVA:
> z <- lm(y ~ x)
> summary(manova(z), test = c("Wilks"))
Df Wilks approx F num Df den Df Pr(>F)
x 5 0.15401 4.9369 20 130.3 7.714e-09 ***
Residuals 42
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comprobamos que el valor del estadístico es \(\Lambda=0.15401\). Para ello, calculamos en primer lugar la matriz \(\bm{E}\)
> res <- z$residuals
> E <- crossprod(res) # Equivalente a (n - 1) * cov(res)
> E
cir4 crec4 cir15 peso15
cir4 0.3199875 1.696564 0.5540875 0.217140
crec4 1.6965637 12.142790 4.3636125 2.110214
cir15 0.5540875 4.363612 4.2908125 2.481656
peso15 0.2171400 2.110214 2.4816562 1.722525Calculamos ahora la matriz \(\bm{E_H}=\bm{E}+\bm{H}\) aplicando la definición.
> resh <- y - matrix(medy, nrow = n, ncol = d, byrow = TRUE)
> EH <- crossprod(resh) # Equivalente a (n - 1) * cov(resh)
> EH
cir4 crec4 cir15 peso15
cir4 0.3935479 2.233949 0.8863521 0.425610
crec4 2.2339490 16.342452 6.7190010 3.747322
cir15 0.8863521 6.719001 10.4047479 6.262700
peso15 0.4256100 3.747322 6.2627000 4.215616Efectivamente, se tiene
En base a la aproximación a la \(F\) de Snédecor proporcionada por R y al \(p\)-valor del estadístico, se rechaza la hipótesis nula de igualdad de los vectores de medias en los 6 tipos de portainjerto.
3.2.3 Test de Roy
Existen otros procedimientos para contrastar \(H_0: \bm{\mu_1}=\bm{\mu_2}=\cdots=\bm{\mu_k}\) que dan lugar a estadísticos de contraste diferentes al \(\Lambda\) de Wilks. Uno de esos procedimientos alternativos es el de unión-intersección.
Un test de unión-intersección para contrastar una hipótesis \(H_0=\cap_a H_{0a}\) es un test cuya región de rechazo es de la forma \(R=\cup_a R_a\), siendo \(R_a\) la región de rechazo para \(H_{0a}\).
Esta idea se aplica en problemas de inferencia multivariante, expresando la hipótesis multivariante como intersección de hipótesis univariantes, para las cuales ya hay tests disponibles. Siguiendo este procedimiento, el test de unión-intersección de Roy tiene como estadístico de contraste:
\[\lambda_1=\mbox{máximo autovalor de}\ \bm{H}\bm{E}^{-1}\] y rechazaremos la hipótesis nula para valores grandes de \(\lambda_1\). De nuevo, se utilizarán aproximaciones de este estadístico a la \(F\) de Snédecor.
Ejemplo. Calculamos, para el ejemplo sobre manzanos, el estadístico de Roy:
> summary(manova(z), test = c("Roy"))
Df Roy approx F num Df den Df Pr(>F)
x 5 1.8757 15.756 5 42 1.002e-08 ***
Residuals 42
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En base a este estadístico también se rechaza la hipótesis nula de igualdad de medias. Comprobamos que el valor del estadístico es \(\lambda_1=1.8757\). Para ello, calculamos los autovalores de \(\bm{H}\bm{E}^{-1}.\)
3.2.4 Tests de Pillai y Lawley-Hotelling
Existen otros dos tests adicionales para contrastar \(H_0: \bm{\mu_1}=\bm{\mu_2}=\cdots=\bm{\mu_k}\) basados en los autovalores de la matriz \(\bm{H}\bm{E}^{-1}\). El estadístico de Pillai viene dado por
\[V^{(^d)}=\textnormal{traza}\left[\bm{H}(\bm{E}+\bm{H})^{-1}\right]=\sum_{i=1}^d\theta_i=\sum_{i=1}^d\frac{\lambda_i}{1+\lambda_i}\] siendo \(\theta_1,\ldots,\theta_d\) los autovalores de \(\bm{H}(\bm{E}+\bm{H})^{-1}\) y \(\lambda_1,\ldots,\lambda_d\) los autovalores de \(\bm{H}\bm{E}^{-1}\).
El estadístico de Lawley-Hotelling se define como
\[U^{(^d)}=\textnormal{traza}\left[\bm{H}\bm{E}^{-1}\right]=\sum_{i=1}^d\lambda_i\] donde \(\lambda_1,\ldots,\lambda_d\) son los autovalores de \(\bm{H}\bm{E}^{-1}\).
Ejemplo. Calculamos, para el ejemplo sobre manzanos, los estadísticos de Pillai y Lawley-Hotelling y comprobamos los valores obtenidos:
> summary(manova(z), test = c("Pillai"))
Df Pillai approx F num Df den Df Pr(>F)
x 5 1.3055 4.0697 20 168 1.983e-07 ***
Residuals 42
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> summary(manova(z), test = c("Hotelling-Lawley"))
Df Hotelling-Lawley approx F num Df den Df Pr(>F)
x 5 2.9214 5.4776 20 150 2.568e-10 ***
Residuals 42
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En base a ambos estadísticos también se rechaza la hipótesis nula de igualdad de medias. Comprobamos el valor del estadístico de Pillai y de Lawley-Hotelling
References
La razón para no comparar directamente \(\bm{H}\) con \(\bm{E}\), sino con \(\bm{E}+\bm{H}\) es que \(\bm{H}\) podría tener algún autovalor cero y eso haría que su determinante fuera cero.↩︎
Ver Anexo A de distribuciones multidimensionales para la definición formal de la distribución \(\Lambda\) de Wilks↩︎