3.3 MANOVA con dos factores de variación
En esta sección suponemos que hay dos factores: A y B. Del factor A podemos distinguir \(a\) niveles, mientras que en el factor B podemos encontrar \(b\) niveles. En cada una de las \(a\cdot b\) posibilidades realizamos \(n\) observaciones de un vector aleatorio \(\bm{y}=(y_1,\ldots,y_d)^\prime{}\). El objetivo será estudiar la influencia de los factores A y B, o de su interacción, en la media del vector \(\bm{y}\). Así, planteamos el siguiente modelo:
\[\bm{y_{ijk}}=\bm{\mu}+\bm{\alpha_i}+\bm{\beta_j}+\bm{\gamma_{ij}}+{\bm{\epsilon_{ijk}}} \qquad i=1,\ldots,a, \ \ j=1,\ldots,b, \ \ k=1,\ldots,n.\] Se asume que \({\bm{\epsilon_{ijk}}}\in N_d(\bm{0},\bm\Sigma)\). El parámetro \(\bm{\mu}\) representa la media global, los parámetros \(\bm{\alpha_i}\) representan el efecto principal del factor A, los parámetros \(\bm{\beta_j}\) representan el efecto principal del factor B y los parámetros \(\bm{\gamma_{ij}}\) representan la interacción de los factores A y B. Además verifican
\[\sum_{i=1}^a \bm{\alpha_i}=\sum_{j=1}^b \bm{\beta_j}=\sum_{i=1}^a \bm{\gamma_{ij}} =\sum_{j=1}^b\bm{\gamma_{ij}}=\bm{0}\]
Este modelo también se puede ver como un caso particular del modelo lineal general multivariante. Las labores de inferencia sobre el modelo se realizan en base a la siguiente descomposición de la variabilidad, mediante lo que llamaremos tabla MANOVA II:
| Fuente de variación | Matriz de covarianzas | Grados de libertad |
|---|---|---|
| Factor A | \({ \bm{H_A} = nb\sum_{i=1}^a\left(\bm{\overline y_{i\bullet\bullet}}-\bm{\overline y_{\bullet\bullet\bullet}}\right) \left(\bm{\overline y_{i\bullet\bullet}}-\bm{\overline y_{\bullet\bullet\bullet}}\right)^\prime{}}\) | \(a-1\) |
| Factor B | \({ \bm{H_B} = na\sum_{j=1}^b\left(\bm{\overline y_{\bullet j \bullet}}-\bm{\overline y_{\bullet\bullet\bullet}}\right)\left(\bm{\overline y_{\bullet j \bullet}}-\bm{\overline y_{\bullet\bullet\bullet}}\right)^\prime{}}\) | \(b-1\) |
| Interacción | \({ \bm{H_{AB}} = n\sum_{i=1}^a\sum_{j=1}^b\left(\bm{\overline{y}_{ij\bullet}}-\bm{\overline y_{i\bullet\bullet}}-\bm{\overline y_{\bullet j\bullet}}+\bm{\overline y_{\bullet\bullet\bullet}}\right)}\) | \((a-1)(b-1)\) |
| \({\times \left(\bm{\overline{y}_{ij\bullet}}-\bm{\overline y_{i\bullet\bullet}} -\bm{\overline y_{\bullet j\bullet}}+\bm{\overline y_{\bullet\bullet\bullet}}\right)^\prime{}}\) | ||
| Error | \({ \bm{E} = \sum_{i=1}^a \sum_{j=1}^b \sum_{k=1}^n\left( \bm{y_{ijk}}-\bm{\overline y_{ij\bullet}}\right)\left( \bm{y_{ijk}}-\bm{\overline y_{ij\bullet}}\right)^\prime{}}\) | \(ab(n-1)\) |
| Total | \({\sum_{i=1}^a \sum_{j=1}^b \sum_{k=1}^n\left( \bm{y_{ijk}}-\bm{\overline y_{\bullet\bullet\bullet}}\right)\left( \bm{y_{ijk}}-\bm{\overline y_{\bullet\bullet\bullet}}\right)^\prime{}}\) | \(abn-1\) |
La hipótesis relativa a la interacción \(H_{AB}:\bm{\gamma_{ij}}=\bm{0}\ \forall i,j\) se contrasta en base al estadístico
\[\frac{|\bm{E}|}{|\bm{E}+\bm{H_{AB}}|}\in \Lambda\left(d,(a-1)(b-1),ab(n-1)\right).\]
Si se obtiene que la interacción no es significativa, se pasaría a un modelo sin interacción, sobre el cual se pueden cuestionar los efectos principales. Esto sería el proceder habitual dentro de un planteamiento jerárquico, cuestión que ya surgía en los modelos univariantes.
Ejemplo. En una especie de césped, denominada Paspalum, se está investigando el efecto que experimenta al ser infectada con un hongo. Al mismo tiempo se tiene en cuenta la temperatura, dentro de un diseño con cuatro valores diferentes: \(14, 18, 22, 26^oC\). En cada realización del experimento se miden tres variables:
\[\begin{aligned}
y_1&=\mbox{El peso fresco de las raíces (medido en gramos)} \\
y_2&=\mbox{La longitud máxima de las raíces (medida en milímetros)} \\
y_3&=\mbox{El peso fresco de las hojas (medido en gramos)}\end{aligned}\]
Los datos se encuentran en el fichero hongos.txt. Vamos a contrastar
el efecto del tratamiento con hongos, el efecto de la temperatura y la
posible interacción entre ambos efectos.
> hongos <- read.table(file = "data/hongos.txt", header = TRUE)
> head(hongos)
tratamiento temperatura y1 y2 y3
1 1 1 2.2 23.5 1.7
2 1 1 3.0 27.0 2.3
3 1 1 3.3 24.5 3.2
4 1 1 2.2 20.5 1.5
5 1 1 2.0 19.0 2.0
6 1 1 3.5 23.5 2.9> y <- cbind(hongos$y1, hongos$y2, hongos$y3)
> tratamiento <- factor(hongos$tratamiento)
> temperatura <- factor(hongos$temperatura)Empezando por el modelo más general, con ambos efectos principales y la interacción, si se contrasta la interacción resulta:
> m12i = lm(y ~ tratamiento * temperatura)
> m12 = lm(y ~ tratamiento + temperatura)
> anova(m12, m12i, test = "Wilks")
Analysis of Variance Table
Model 1: y ~ tratamiento + temperatura
Model 2: y ~ tratamiento * temperatura
Res.Df Df Gen.var. Wilks approx F num Df den Df Pr(>F)
1 43 12.852
2 40 -3 12.299 0.70548 1.5864 9 92.633 0.1307> anova(m12, m12i, test = "Roy")
Analysis of Variance Table
Model 1: y ~ tratamiento + temperatura
Model 2: y ~ tratamiento * temperatura
Res.Df Df Gen.var. Roy approx F num Df den Df Pr(>F)
1 43 12.852
2 40 -3 12.299 0.25155 3.3541 3 40 0.02819 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> anova(m12, m12i, test = "Hotelling-Lawley")
Analysis of Variance Table
Model 1: y ~ tratamiento + temperatura
Model 2: y ~ tratamiento * temperatura
Res.Df Df Gen.var. Hotelling-Lawley approx F num Df den Df Pr(>F)
1 43 12.852
2 40 -3 12.299 0.38126 1.5533 9 110 0.1384> anova(m12, m12i, test = "Pillai")
Analysis of Variance Table
Model 1: y ~ tratamiento + temperatura
Model 2: y ~ tratamiento * temperatura
Res.Df Df Gen.var. Pillai approx F num Df den Df Pr(>F)
1 43 12.852
2 40 -3 12.299 0.32058 1.5953 9 120 0.1241Únicamente el test de Roy produce un nivel crítico algo pequeño, mientras que los otros tests no muestran significación. En consecuencia, se suprime la interacción del modelo. En el modelo sin interacción y con los dos efectos principales, se contrasta la significación de cada uno de los efectos, con los resultados siguientes. Para el tratamiento:
> m1 = lm(y ~ tratamiento)
> anova(m1, m12, test = "Wilks")
Analysis of Variance Table
Model 1: y ~ tratamiento
Model 2: y ~ tratamiento + temperatura
Res.Df Df Gen.var. Wilks approx F num Df den Df Pr(>F)
1 46 29.940
2 43 -3 12.852 0.064601 23.12 9 99.934 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> anova(m1, m12, test = "Roy")
Analysis of Variance Table
Model 1: y ~ tratamiento
Model 2: y ~ tratamiento + temperatura
Res.Df Df Gen.var. Roy approx F num Df den Df Pr(>F)
1 46 29.940
2 43 -3 12.852 5.1828 74.287 3 43 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> anova(m1, m12, test = "Hotelling-Lawley")
Analysis of Variance Table
Model 1: y ~ tratamiento
Model 2: y ~ tratamiento + temperatura
Res.Df Df Gen.var. Hotelling-Lawley approx F num Df den Df Pr(>F)
1 46 29.940
2 43 -3 12.852 6.6858 29.467 9 119 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> anova(m1, m12, test = "Pillai")
Analysis of Variance Table
Model 1: y ~ tratamiento
Model 2: y ~ tratamiento + temperatura
Res.Df Df Gen.var. Pillai approx F num Df den Df Pr(>F)
1 46 29.940
2 43 -3 12.852 1.4391 13.215 9 129 8.45e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Para la temperatura:
> m2 = lm(y ~ temperatura)
> anova(m2, m12, test = "Wilks")
Analysis of Variance Table
Model 1: y ~ temperatura
Model 2: y ~ tratamiento + temperatura
Res.Df Df Gen.var. Wilks approx F num Df den Df Pr(>F)
1 44 13.439
2 43 -1 12.852 0.8163 3.0756 3 41 0.03807 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> anova(m2, m12, test = "Roy")
Analysis of Variance Table
Model 1: y ~ temperatura
Model 2: y ~ tratamiento + temperatura
Res.Df Df Gen.var. Roy approx F num Df den Df Pr(>F)
1 44 13.439
2 43 -1 12.852 0.22504 3.0756 3 41 0.03807 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> anova(m2, m12, test = "Hotelling-Lawley")
Analysis of Variance Table
Model 1: y ~ temperatura
Model 2: y ~ tratamiento + temperatura
Res.Df Df Gen.var. Hotelling-Lawley approx F num Df den Df Pr(>F)
1 44 13.439
2 43 -1 12.852 0.22504 3.0756 3 41 0.03807 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1> anova(m2, m12, test = "Pillai")
Analysis of Variance Table
Model 1: y ~ temperatura
Model 2: y ~ tratamiento + temperatura
Res.Df Df Gen.var. Pillai approx F num Df den Df Pr(>F)
1 44 13.439
2 43 -1 12.852 0.1837 3.0756 3 41 0.03807 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Por tanto, la temperatura tiene un efecto muy significativo sobre las variables que se están midiendo. El tratamiento también resulta significativo, aunque no es tan notorio. Los niveles críticos de los cuatro test coinciden para el contraste del tratamiento. Este hecho se debe a que hay únicamente dos tratamientos, por lo que la matriz \(H\) asociada al contraste tiene rango uno.