1.- INTRODUCCIÓN 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:

2.- EJERCICIO

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   H

2.1.- ANÁLISIS EXPLORATORIO

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
library(ggpubr)
## Warning: package 'ggpubr' was built under R version 4.2.3
ggboxplot(ratas, x = "sex", y = c("y1", "y2"),
          merge = T, palette="jco")

Y1 “El peso antes” es similar para hombres y mujeres, la media del peso antes es similar para ambos casos. Para el peso antes la dispersión de los hombres es mayor a la de las mujeres. Y2 “El peso después” es similar para hombres y mujeres, la media del peso antes es similar para ambos casos. En el peso de después, los hombres presentan menor dispersión que las mujeres.

Para los hombres el peso antes es similar al peso de después. La dispersión de los hombres antes es muy superior a la de después.

Para las mujeres el peso antes es similar al peso de después. La dispersión de las mujeres antes es muy inferior a la de después.

ggboxplot(ratas, x = "temp", y = c("y1", "y2"),
          merge = T, palette="aaas")

Y1 “El peso antes” para la duración de 4 se obitene un promedio que para 20 y 34. Comparando 20 y 34, tienen una media superior a 4, pero entre ellas son promedios similares. En cuanto a la dispersión 4, presenta la menor dispersión, seguido de 20 y la mayor dispersión es para 34. Por tanto, a medida que el tratamiento es más dilatado en el tiempo, la dispersión del peso antes va aumentando.

Y2 “El peso después” se obtiene que el promedio más bajo es para 34, después tendríamos a 4 y por último a 20. La dispersión más elevada la tenemos para 34, después para 20 y por último para 4. Por tanto, a medida que el tratamiento es más dilatado en el tiempo, la dispersión del peso va aumentando.

Temp = 4. Se observa que la media de antes es inferior a la media de después. La dispersión es mucho mayor para el peso de antes.

Temp = 20. Se observa que la media de antes es similar a la media de después (porque comparten valores). La dispersión es ligeramente mayor para el peso de antes.

Temp = 34. Se observa que la media de antes es similar a la media después (porque comparten valores). La dispersión es ligeramente mayor para el peso de antes.

aggregate(ratas[,1:2], list(ratas$temp, ratas$sex), mean)
##   Group.1 Group.2       y1       y2
## 1       4       H 19.49000 19.69667
## 2      20       H 20.79000 20.12333
## 3      34       H 18.81000 17.90667
## 4       4       M 18.66000 18.50000
## 5      20       M 19.90333 22.53333
## 6      34       M 20.32667 18.31667
aggregate(ratas[,1:2], list(ratas$temp, ratas$sex), var)
##   Group.1 Group.2        y1         y2
## 1       4       H 1.2451000 0.04963333
## 2      20       H 0.5097000 1.16943333
## 3      34       H 3.4693000 4.00013333
## 4       4       M 0.2503000 2.97010000
## 5      20       M 1.5233333 0.46333333
## 6      34       M 0.1381333 2.00143333
aggregate(ratas[,1:2], list(ratas$temp, ratas$sex), sd)
##   Group.1 Group.2        y1        y2
## 1       4       H 1.1158405 0.2227854
## 2      20       H 0.7139328 1.0814034
## 3      34       H 1.8626057 2.0000333
## 4       4       M 0.5002999 1.7233978
## 5      20       M 1.2342339 0.6806859
## 6      34       M 0.3716629 1.4147202
aggregate(ratas[,1:2], list(ratas$temp, ratas$sex), length)
##   Group.1 Group.2 y1 y2
## 1       4       H  3  3
## 2      20       H  3  3
## 3      34       H  3  3
## 4       4       M  3  3
## 5      20       M  3  3
## 6      34       M  3  3
# Normal univariada
normal <- function(vec){
  shapiro.test(vec)$p
}
aggregate(ratas[,1:2], list(ratas$sex, ratas$temp), normal)
##   Group.1 Group.2        y1          y2
## 1       H       4 0.9258674 0.128686974
## 2       M       4 0.9338670 0.005540978
## 3       H      20 0.6975632 0.666926507
## 4       M      20 0.5491264 0.424350926
## 5       H      34 0.6155710 0.988973545
## 6       M      34 0.5202800 0.175746733

En todos los casos, los p valores son mayores a 0.05. Todos los casos aceptan la hipótesis de normalidad. Solamente para m & 4 se rechaza la hipótesis de normalidad, pues su p valor es 0.005 < 0.05.

# Normal multivariada
library(MVN)
## Warning: package 'MVN' was built under R version 4.2.3
mvn(ratas[,c(1,2)], mvnTest = "hz")
## $multivariateNormality
##            Test        HZ   p value MVN
## 1 Henze-Zirkler 0.3144069 0.6573042 YES
## 
## $univariateNormality
##               Test  Variable Statistic   p value Normality
## 1 Anderson-Darling    y1        0.2429    0.7291    YES   
## 2 Anderson-Darling    y2        0.4135    0.3021    YES   
## 
## $Descriptives
##     n     Mean  Std.Dev Median  Min   Max    25th  75th        Skew   Kurtosis
## y1 18 19.66333 1.209458 19.795 17.2 21.56 18.7275 20.64 -0.26369339 -1.0800902
## y2 18 19.51278 1.955492 19.495 15.9 23.30 18.9625 20.23 -0.03861505 -0.5963898

El HZ tiene un estadístico de 0.31 con un p valor de 0.65 > 0.05. Tenemos evidencia empírica suficiente para aceptar la hipótesis de normalidad multivariada.

El Anderson Darling de y1 tiene un estadístico de 0.24 con p valor de 0.729 > 0.05. Tenemos evidencia empírica suficiente para aceptar la hipótesis de normalidad de y1.

El Anderson Darling de y2 tiene un estadístico de 0.41 con p valor de 0.302 > 0.05. Tenemos evidencia empírica suficiente para aceptar la hipótesis de normalidad de y2.

Se acepta la hipótesis de normalidad.

library(biotools)
## Warning: package 'biotools' was built under R version 4.2.3
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 4.2.3
## ---
## biotools version 4.2
library(car)
## Warning: package 'car' was built under R version 4.2.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.2.3
library(carData)
boxM(ratas[, c(1,2)], grouping = ratas [,c(4)])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  ratas[, c(1, 2)]
## Chi-Sq (approx.) = 2.6213, df = 6, p-value = 0.8547

“Y1” El estadístico es 2.62 con p valor de 0.85. Se acepta la hipótesis de que las matrices de varianzas y covarianzas son similares (son esféricas).

boxM(ratas[, c(1,2)], grouping = ratas [,c(5)])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  ratas[, c(1, 2)]
## Chi-Sq (approx.) = 5.6997, df = 3, p-value = 0.1272

“Y2” El estadístico es 5.69 con p valor de 0.12. Se acepta la hipótesis de que las matrices de varianzas y covarianzas son similares (son esféricas).

leveneTest(ratas$y1~ratas$temp)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.3282 0.7253
##       15

Se acepta la hipótesis de homogeneidad de varianzas.

leveneTest(ratas$y1~ratas$sex)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.7524 0.3985
##       16

Se acepta la hipótesis de homogeneidad de varianzas.

leveneTest(ratas$y2~ratas$temp)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.8605 0.4428
##       15

Se acepta la hipótesis de homogeneidad de varianzas.

leveneTest(ratas$y2~ratas$sex)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.2716 0.2761
##       16

Se acepta la hipótesis de homogeneidad de varianzas.

Por tanto, se acepta la hipótesis de homocedasticidad o esfericidad.

##2.2.- MANOVA

modelo1<- manova(cbind(y1, y2)~sex+temp, data =ratas)
summary(modelo1)
##           Df  Pillai approx F num Df den Df  Pr(>F)  
## sex        1 0.05292  0.36322      2     13 0.70228  
## temp       2 0.61789  3.12945      4     28 0.03017 *
## Residuals 14                                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se estima que diferenciando el sexo el estadístico es 0.36 con p valor de 0.70>0.05. Tenemos evidencia de que el peso antes y después es similar independientemente del sexo.

Se estima que diferenciando el tiempo el estadístico es 3.12 con p valor de 0.03>0.05. Tenemos evidencia de que el peso antes y después es diferente dependiendo del tiempo del tratamiento.

# Modelo con interacción
modelo2 <- manova(cbind(y1,y2)~sex*temp, data=ratas)
summary(modelo2)
##           Df  Pillai approx F num Df den Df  Pr(>F)  
## sex        1 0.09768   0.5954      2     11 0.56819  
## temp       2 0.78457   3.8730      4     24 0.01448 *
## sex:temp   2 0.71498   3.3384      4     24 0.02618 *
## Residuals 12                                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se estima que diferenciando el sexo el estadístico es 0.59 con p valor de 0.56>0.05. Tenemos evidencia de que el peso antes y después es similar independientemente del sexo.

Se estima que diferenciando el tiempo el estadístico es 3.87 con p valor de 0.01>0.05. Tenemos evidencia de que el peso antes y después es diferente dependiendo del tiempo del tratamiento.

Para la interacción entre el sexo y el tiempo de tratamiento, el estadístico asciende a 3.33 con un p valor de 0.026 < 0.05. Tenemos evidencia empírica para rechazar la hipótesis nula; por lo tanto, la interacción de estos factores proporciona media del peso antes y después distintas.

##2.3.- Modelos univariados ANOVA

m1 <- aov(y1~temp)
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.329

Se obtiene que el peso de antes es similar independientemente del tiempo del tratamiento.

m2 <- aov(y2~temp)
summary(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 ' ' 1

Se obtiene que el peso después es diferente dependiendo del tiempo de tratamiento.

m3 <- aov(y1~sex)
summary(m3)
##             Df Sum Sq Mean Sq F value Pr(>F)
## sex          1   0.02   0.020   0.013  0.911
## Residuals   16  24.85   1.553

Se obtiene que el peso de antes es similar, independientemente del sexo del individuo.

m4 <- aov(y2~sex)
summary(m4)
##             Df Sum Sq Mean Sq F value Pr(>F)
## sex          1   1.32   1.318   0.331  0.573
## Residuals   16  63.69   3.981

Se obtiene que el peso de después es similar, independientemente del sexo del individuo.

m5 <- aov(y1~temp*sex)
summary(m5)
##             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.189

Las medias son iguales para el peso de antes analizando el tiempo, el sexo y su interacción.

m6 <- aov(y2~temp*sex)
summary(m6)
##             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 ' ' 1

Las medias para el peso de después son diferentes para el tiempo de tratamiento, no es así para el sexo o la interacción entre el sexo y el tiempo de tratamiento.

library(lsr)
## Warning: package 'lsr' was built under R version 4.2.3
etaSquared(m1)
##         eta.sq eta.sq.part
## temp 0.1983574   0.1983574

El efecto es reducido (<25%)

etaSquared(m2)
##         eta.sq eta.sq.part
## temp 0.5012788   0.5012788

El efecto es grande (>25%).

etaSquared(m3)
##           eta.sq  eta.sq.part
## sex 0.0008042658 0.0008042658

El efecto es reducido (<25%).

etaSquared(m4)
##         eta.sq eta.sq.part
## sex 0.02026862  0.02026862

El efecto es reducido (<25%)

etaSquared(m5)
##                eta.sq eta.sq.part
## temp     0.1983574211  0.25684957
## sex      0.0008042658  0.00139941
## temp:sex 0.2269249432  0.28335925

Para el tiempo sería grande (>25%). Para sexo sería reducido (<25%). Para la interacción sería grande (>25%).

etaSquared(m6)
##              eta.sq eta.sq.part
## temp     0.50127879  0.60463503
## sex      0.02026862  0.05823481
## temp:sex 0.15067126  0.31491368

Para el tiempo sería grande (>25%). Para sexo sería reducido (<25%). Para la interacción sería grande (>25%).

##2.4.- Análisis Post-Hoc

ph1 <- TukeyHSD(m1)
ph1
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = y1 ~ temp)
## 
## $temp
##             diff        lwr       upr     p adj
## 20-4   1.2716667 -0.4571538 3.0004872 0.1699245
## 34-4   0.4933333 -1.2354872 2.2221538 0.7433607
## 34-20 -0.7783333 -2.5071538 0.9504872 0.4884990
plot(ph1)

En este caso, el valor de cero es compatible con los tres intervalos. Se obtienen las mismas medias en cada caso. El peso antes es similar para cualquier tiempo de tratamiento.

ph2 <- TukeyHSD(m2)
ph2
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = y2 ~ temp)
## 
## $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.0047632
plot(ph2)

La comparación entre 34-20 proporciona una media del peso después significativamente inferior al resto.

El IC al 95%, comparando la media de 34 con 20 toma valores entre [-5.42 y -1.011], ambos extremos son negativos, indicando que la media del peso de después del tiempo de 34 es inferior a la media del peso de después de tiempo 20.

ph3 <- TukeyHSD(m3)
ph3
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = y1 ~ sex)
## 
## $sex
##            diff       lwr      upr     p adj
## M-H -0.06666667 -1.312015 1.178681 0.9110586
plot(ph3)

No hay diferencias en el peso de antes entre hombres y mujeres.

ph4 <- TukeyHSD(m4)
ph4
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = y2 ~ sex)
## 
## $sex
##          diff       lwr      upr     p adj
## M-H 0.5411111 -1.452701 2.534924 0.5730696
plot(ph4)

No hay diferencias en el peso de después entre hombres y mujeres.

ph5 <- TukeyHSD(m5)
ph5
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = y1 ~ temp * sex)
## 
## $temp
##             diff        lwr       upr     p adj
## 20-4   1.2716667 -0.4081069 2.9514402 0.1497835
## 34-4   0.4933333 -1.1864402 2.1731069 0.7197791
## 34-20 -0.7783333 -2.4581069 0.9014402 0.4556772
## 
## $sex
##            diff       lwr      upr     p adj
## M-H -0.06666667 -1.186779 1.053445 0.8989693
## 
## $`temp:sex`
##                 diff       lwr       upr     p adj
## 20:H-4:H   1.3000000 -1.690900 4.2908999 0.6934084
## 34:H-4:H  -0.6800000 -3.670900 2.3108999 0.9687400
## 4:M-4:H   -0.8300000 -3.820900 2.1608999 0.9303733
## 20:M-4:H   0.4133333 -2.577567 3.4042333 0.9966097
## 34:M-4:H   0.8366667 -2.154233 3.8275666 0.9282065
## 34:H-20:H -1.9800000 -4.970900 1.0108999 0.2951324
## 4:M-20:H  -2.1300000 -5.120900 0.8608999 0.2326552
## 20:M-20:H -0.8866667 -3.877567 2.1042333 0.9106784
## 34:M-20:H -0.4633333 -3.454233 2.5275666 0.9942430
## 4:M-34:H  -0.1500000 -3.140900 2.8408999 0.9999754
## 20:M-34:H  1.0933333 -1.897567 4.0842333 0.8157131
## 34:M-34:H  1.5166667 -1.474233 4.5075666 0.5542845
## 20:M-4:M   1.2433333 -1.747567 4.2342333 0.7287434
## 34:M-4:M   1.6666667 -1.324233 4.6575666 0.4612261
## 34:M-20:M  0.4233333 -2.567567 3.4142333 0.9962095
plot(ph5)

La interacción no es significativa, ya que todos los intervalos son compatibles con el valor de cero.

ph6 <- TukeyHSD(m6)
ph6
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = y2 ~ temp * sex)
## 
## $temp
##             diff        lwr       upr     p adj
## 20-4   2.2300000  0.1774909  4.282509 0.0332276
## 34-4  -0.9866667 -3.0391757  1.065842 0.4309489
## 34-20 -3.2166667 -5.2691757 -1.164158 0.0033709
## 
## $sex
##          diff        lwr      upr     p adj
## M-H 0.5411111 -0.8275494 1.909772 0.4058955
## 
## $`temp:sex`
##                 diff        lwr        upr     p adj
## 20:H-4:H   0.4266667 -3.2279028  4.0812361 0.9984713
## 34:H-4:H  -1.7900000 -5.4445695  1.8645695 0.5874694
## 4:M-4:H   -1.1966667 -4.8512361  2.4579028 0.8723511
## 20:M-4:H   2.8366667 -0.8179028  6.4912361 0.1687646
## 34:M-4:H  -1.3800000 -5.0345695  2.2745695 0.7959119
## 34:H-20:H -2.2166667 -5.8712361  1.4379028 0.3774021
## 4:M-20:H  -1.6233333 -5.2779028  2.0312361 0.6752967
## 20:M-20:H  2.4100000 -1.2445695  6.0645695 0.2986266
## 34:M-20:H -1.8066667 -5.4612361  1.8479028 0.5786896
## 4:M-34:H   0.5933333 -3.0612361  4.2479028 0.9928641
## 20:M-34:H  4.6266667  0.9720972  8.2812361 0.0110897
## 34:M-34:H  0.4100000 -3.2445695  4.0645695 0.9987361
## 20:M-4:M   4.0333333  0.3787639  7.6879028 0.0278169
## 34:M-4:M  -0.1833333 -3.8379028  3.4712361 0.9999753
## 34:M-20:M -4.2166667 -7.8712361 -0.5620972 0.0209219
plot(ph6)

Existen intervalos incompatibles con el valor de cero. Indicando que la interacción es estadísticamente significativa.

20:M-34:H 4.6266667 0.9720972 8.2812361 0.0110897 Positivos

34:M-20:M -4.2166667 -7.8712361 -0.5620972 0.0209219 Negativos

La media de los hombres con 20 es significativamente superior a la media de los hombres de 34.

El peso de después sólo se incrementa cuando un hombre reduce su tratamiento de 20 a 34.

La media de las mujeres de 34 es inferior a la media del peso dd después de las mujeres de 20.

El peso de las mujeres después del tratamiento cuando pasan de 34 a 20 disminuye.

#3.- INFORME DE LA SESIÓN

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)
## # A tibble: 90 × 3
##    package   loadedversion source        
##    <chr>     <chr>         <chr>         
##  1 abind     1.4-5         CRAN (R 4.2.0)
##  2 backports 1.4.1         CRAN (R 4.2.0)
##  3 biotools  4.2           CRAN (R 4.2.3)
##  4 boot      1.3-28        CRAN (R 4.2.2)
##  5 broom     1.0.5         CRAN (R 4.2.3)
##  6 bslib     0.4.2         CRAN (R 4.2.2)
##  7 cachem    1.0.6         CRAN (R 4.2.2)
##  8 callr     3.7.3         CRAN (R 4.2.3)
##  9 car       3.1-2         CRAN (R 4.2.3)
## 10 carData   3.0-5         CRAN (R 4.2.3)
## # ℹ 80 more rows