EJemplo de MANOVA
Luis Bolaños
9/5/2020
Manova
El análisis de varianza multivariante ( MANOVA ) es un ANOVA con dos o más variables de resultado (o respuesta) continuas.
El MANOVA unidireccional prueba simultáneamente las diferencias estadísticas para las variables de respuesta múltiple mediante una agrupación de variables.
Por ejemplo, podemos realizar un experimento en el que le damos dos tratamientos (A y B) a dos grupos de ratones, y estamos interesados en el peso y la altura de los ratones. En ese caso, el peso y la altura de los ratones son nuestras variables de resultado (o dependientes), y nuestra hipótesis es que ambos juntos se ven afectados por la diferencia en el tratamiento. Se podría utilizar un análisis de varianza multivariado para probar esta hipótesis.
El procedimiento de MANOVA se puede resumir de la siguiente manera:
- Cree una nueva variable compuesta que sea una combinación lineal de todas las variables de respuesta.
- Compare los valores medios de esta nueva variable entre grupos.
Ejemplo
En un experimento para inhibir un tumor, se quiere investigar el efecto del sexo (S) y de la temperatura ambiental (T). Se consideran las variables: Y1 =peso inicial, Y2 =peso final, Y3 =peso del tumor.
pesos <-c (18.15,16.51, 0.24, 19.15, 19.49, 0.16, 18.68,
19.50, 0.32, 18.35, 19.81, 0.17, 19.54, 19.84,
0.20, 20.58, 19.44, 0.22,21.27, 23.30, 0.33,
18.87, 22.00, 0.25, 19.57, 22.30, 0.45, 20.66,
21.08, 0.20, 20.15, 18.95, 0.35, 21.56, 20.34,
0.20,20.74, 16.69, 0.31, 20.22, 19.00, 0.18,
20.02, 19.26, 0.41, 18.38, 17.92, 0.30, 17.20,
15.90, 0.28, 20.85, 19.90, 0.17)
y1<-pesos[c(T,F,F)]
y2<-pesos[c(F,T,F)]
y3<-pesos[c(F,F,T)]
temp <-factor(c(rep(4,6),rep(20,6),rep(34,6)))
sex<-factor(rep(c(rep("M",3),rep("H",3)),3))
ratas<-data.frame(y1,y2,y3,temp,sex)
head(ratas)## y1 y2 y3 temp sex
## 1 18.15 16.51 0.24 4 M
## 2 19.15 19.49 0.16 4 M
## 3 18.68 19.50 0.32 4 M
## 4 18.35 19.81 0.17 4 H
## 5 19.54 19.84 0.20 4 H
## 6 20.58 19.44 0.22 4 Hvisualización
ggboxplot(ratas,x="sex",y=c("y1","y2"),
merge=T,palette = "jco")
ggboxplot(ratas,x="temp",y=c("y1","y2"),
merge=T,palette = "jco")
Cálculo de estadísticas
aggregate(ratas[,1:3],list(ratas$sex,ratas$temp),mean)## Group.1 Group.2 y1 y2 y3
## 1 H 4 19.49000 19.69667 0.1966667
## 2 M 4 18.66000 18.50000 0.2400000
## 3 H 20 20.79000 20.12333 0.2500000
## 4 M 20 19.90333 22.53333 0.3433333
## 5 H 34 18.81000 17.90667 0.2500000
## 6 M 34 20.32667 18.31667 0.3000000aggregate(ratas[,1:3],list(ratas$sex,ratas$temp),length)## Group.1 Group.2 y1 y2 y3
## 1 H 4 3 3 3
## 2 M 4 3 3 3
## 3 H 20 3 3 3
## 4 M 20 3 3 3
## 5 H 34 3 3 3
## 6 M 34 3 3 3aggregate(ratas[,1:3],list(ratas$sex,ratas$temp),sd)## Group.1 Group.2 y1 y2 y3
## 1 H 4 1.1158405 0.2227854 0.02516611
## 2 M 4 0.5002999 1.7233978 0.08000000
## 3 H 20 0.7139328 1.0814034 0.08660254
## 4 M 20 1.2342339 0.6806859 0.10066446
## 5 H 34 1.8626057 2.0000333 0.07000000
## 6 M 34 0.3716629 1.4147202 0.11532563Normalidad univariada
normal<-function(vec){
shapiro.test(vec)$p
}
aggregate(ratas[,1:3],list(ratas$sex,ratas$temp),normal)## Group.1 Group.2 y1 y2 y3
## 1 H 4 0.9258674 0.128686974 0.7804408
## 2 M 4 0.9338670 0.005540978 1.0000000
## 3 H 20 0.6975632 0.666926507 0.0000000
## 4 M 20 0.5491264 0.424350926 0.7804408
## 5 H 34 0.6155710 0.988973545 0.2737737
## 6 M 34 0.5202800 0.175746733 0.8564460Normalidad Multiple
library(MVN)
mvn(ratas[,c(1,2,3)],mvnTest = "hz")## $multivariateNormality
## Test HZ p value MVN
## 1 Henze-Zirkler 0.7029199 0.151332 YES
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Shapiro-Wilk y1 0.9686 0.7711 YES
## 2 Shapiro-Wilk y2 0.9574 0.5518 YES
## 3 Shapiro-Wilk y3 0.9240 0.1519 YES
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th Skew
## y1 18 19.6633333 1.20945783 19.795 17.20 21.56 18.7275 20.6400 -0.26369339
## y2 18 19.5127778 1.95549243 19.495 15.90 23.30 18.9625 20.2300 -0.03861505
## y3 18 0.2633333 0.08574929 0.245 0.16 0.45 0.2000 0.3175 0.59716549
## Kurtosis
## y1 -1.0800902
## y2 -0.5963898
## y3 -0.8050386Pruebas de varianzas y covarianzas
library(biotools)## ---
## biotools version 3.1library(car)
boxM(ratas[,c(1,2,3)],grouping = ratas[,c(4)])##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: ratas[, c(1, 2, 3)]
## Chi-Sq (approx.) = 3.7453, df = 12, p-value = 0.9876boxM(ratas[,c(1,2,3)],grouping = ratas[,c(5)])##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: ratas[, c(1, 2, 3)]
## Chi-Sq (approx.) = 11.034, df = 6, p-value = 0.08734leveneTest(ratas$y1~ratas$sex)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.7524 0.3985
## 16leveneTest(ratas$y1~ratas$temp)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.3282 0.7253
## 15leveneTest(ratas$y2~ratas$sex)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.2716 0.2761
## 16leveneTest(ratas$y2~ratas$temp)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.8605 0.4428
## 15leveneTest(ratas$y3~ratas$sex)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.3213 0.2673
## 16leveneTest(ratas$y3~ratas$temp)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.0105 0.3875
## 15Modelo Manova
mod1<-manova(cbind(y1,y2,y3)~sex+temp,data=ratas)
summary(mod1)## Df Pillai approx F num Df den Df Pr(>F)
## sex 1 0.21483 1.0945 3 12 0.38905
## temp 2 0.80234 2.9030 6 26 0.02656 *
## Residuals 14
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#modelo con interacción
mod2<-manova(cbind(y1,y2,y3)~sex*temp,data=ratas)
summary(mod2)## Df Pillai approx F num Df den Df Pr(>F)
## sex 1 0.27510 1.2650 3 10 0.33838
## temp 2 0.97019 3.4544 6 22 0.01475 *
## sex:temp 2 0.75271 2.2128 6 22 0.08041 .
## Residuals 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Modelos univariados
m1<-aov(y1~temp,data=ratas)
m2<-aov(y2~temp,data=ratas)
m3<-aov(y3~temp,data=ratas)
m4<-aov(y1~temp*sex,data=ratas)
m5<-aov(y2~temp*sex,data=ratas)
m6<-aov(y3~temp*sex,data=ratas)
summary(m1)## Df Sum Sq Mean Sq F value Pr(>F)
## temp 2 4.933 2.466 1.856 0.19
## Residuals 15 19.935 1.329summary(m2)## Df Sum Sq Mean Sq F value Pr(>F)
## temp 2 32.59 16.293 7.538 0.00542 **
## Residuals 15 32.42 2.161
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m3)## Df Sum Sq Mean Sq F value Pr(>F)
## temp 2 0.01963 0.009817 1.398 0.278
## Residuals 15 0.10537 0.007024summary(m4)## Df Sum Sq Mean Sq F value Pr(>F)
## temp 2 4.933 2.466 2.074 0.168
## sex 1 0.020 0.020 0.017 0.899
## temp:sex 2 5.643 2.821 2.372 0.135
## Residuals 12 14.272 1.189summary(m5)## Df Sum Sq Mean Sq F value Pr(>F)
## temp 2 32.59 16.293 9.176 0.00382 **
## sex 1 1.32 1.318 0.742 0.40590
## temp:sex 2 9.79 4.897 2.758 0.10339
## Residuals 12 21.31 1.776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m6)## Df Sum Sq Mean Sq F value Pr(>F)
## temp 2 0.01963 0.009817 1.374 0.290
## sex 1 0.01742 0.017422 2.439 0.144
## temp:sex 2 0.00221 0.001106 0.155 0.858
## Residuals 12 0.08573 0.007144# análisis post hoc
etaSquared(m2)## eta.sq eta.sq.part
## temp 0.5012788 0.5012788etaSquared(m5)## eta.sq eta.sq.part
## temp 0.50127879 0.60463503
## sex 0.02026862 0.05823481
## temp:sex 0.15067126 0.31491368post<-TukeyHSD(m2)
post## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = y2 ~ temp, data = ratas)
##
## $temp
## diff lwr upr p adj
## 20-4 2.2300000 0.02527797 4.434722 0.0472590
## 34-4 -0.9866667 -3.19138869 1.218055 0.4924782
## 34-20 -3.2166667 -5.42138869 -1.011945 0.0047632plot(post)