MANOVA Y PERMANOVA

library(readxl)

df <- read_excel("~/ANDREA/Universidad/metodos multivariados/Clase 7/Analisis Estadistico Datos Laboratorio 21.08.24.xlsx", 
    sheet = "Consolidado proximal MANOVA")

View(df)

#### revisar que tan balanceado estan mis tratamientos
table(df$Tto)

T1 T2 T3 T4 T5 T6 T7 
 9  9  9  9  9  9  9 

Matriz de correlación de Pearson

library(corrplot)

# seleccionar solo variables numericas
df_num <- df[sapply(df, is.numeric)]

# matriz de correlacion
cor_matrix <- cor(df_num, use = "complete.obs", method = "spearman")### o pearson

# grafica con valores de correlacion
corrplot(cor_matrix,
         method = "color",
         type = "upper",
         tl.col = "black",
         tl.srt = 45,
         addCoef.col = "black",   # muestra el valor de r
         number.cex = 0.8)        # tamaño del numero

# convertir Tto en factor
df$Tto <- as.factor(df$Tto)

MANOVA

# 1. Asegurarnos de que el Tratamiento sea un factor (categoría)
df$Tto <- as.factor(df$Tto)

# 2. Unir las variables dependientes en una sola matriz
# Usamos cbind() para "pegar" las columnas que queremos analizar juntas
variables_respuesta <- cbind(df$Grasas, df$Humedad, df$Proteina, df$Cenizas)

# 3. Ajustar el modelo MANOVA
# La fórmula significa: "Analizar el conjunto de variables en función del Tratamiento"
modelo_manova <- manova(variables_respuesta ~ Tto, data = df)

# 4. Ver el resultado global del MANOVA
# Por defecto, R usa la prueba Traza de Pillai (Pillai's trace), que es la más robusta
summary(modelo_manova)

          Df  Pillai approx F num Df den Df   Pr(>F)   
Tto        6 0.68727   1.9363     24    224 0.007133 **
Residuals 56                                           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Existe evidencia de que al menos uno de los tratamientos presenta una combinación promedio de grasas, humedad, proteína y cenizas diferente de la de los demás tratamientos.

# Extraer la matriz de residuales del MANOVA
residuales_manova <- residuals(modelo_manova)

# Ver las primeras filas
head(residuales_manova)

       [,1]        [,2]       [,3]         [,4]
1 -1.117767 -1.70445140 -0.3144008  0.008991832
2  2.831233  0.04888266  0.2395623 -0.048837767
3  7.519033 -1.41858384  0.6670862  0.008084027
4  2.650233  0.40893784  0.5389723  0.017541847
5 -4.826367  0.58546408  0.7482227 -0.002550635
6 -2.346867  0.33654203  0.7293994 -0.007978019

# Revisar dimensiones
dim(residuales_manova)

[1] 63  4

colnames(residuales_manova) <- c(
  "Grasas",
  "Humedad",
  "Proteina",
  "Cenizas"
)

head(residuales_manova)

     Grasas     Humedad   Proteina      Cenizas
1 -1.117767 -1.70445140 -0.3144008  0.008991832
2  2.831233  0.04888266  0.2395623 -0.048837767
3  7.519033 -1.41858384  0.6670862  0.008084027
4  2.650233  0.40893784  0.5389723  0.017541847
5 -4.826367  0.58546408  0.7482227 -0.002550635
6 -2.346867  0.33654203  0.7293994 -0.007978019

Prueba de normalidad multivariada

library(mvnormtest)

# =========================================================
# NORMALIDAD MULTIVARIANTE DE LOS RESIDUALES DEL MANOVA
# =========================================================

# Extraer residuales
residuales_manova <- residuals(modelo_manova)

# Asignar nombres
colnames(residuales_manova) <- c(
  "Grasas",
  "Humedad",
  "Proteina",
  "Cenizas"
)

# Revisar la matriz
head(residuales_manova)

     Grasas     Humedad   Proteina      Cenizas
1 -1.117767 -1.70445140 -0.3144008  0.008991832
2  2.831233  0.04888266  0.2395623 -0.048837767
3  7.519033 -1.41858384  0.6670862  0.008084027
4  2.650233  0.40893784  0.5389723  0.017541847
5 -4.826367  0.58546408  0.7482227 -0.002550635
6 -2.346867  0.33654203  0.7293994 -0.007978019

dim(residuales_manova)

[1] 63  4

# Instalar una sola vez
# install.packages("mvnormtest")

library(mvnormtest)

# mshapiro.test requiere:
# filas = variables
# columnas = observaciones
resultado_normalidad <- mshapiro.test(
  t(residuales_manova)
)

resultado_normalidad

    Shapiro-Wilk normality test

data:  Z
W = 0.9164, p-value = 0.0003942

La prueba multivariante de Shapiro-Wilk mostró que los residuales no siguen una distribución normal multivariante (p=0.0003942) a pesar de que W=0.9164 sugiere que la desviación de la normalidad no es extrema.

Comparación de matriz de varianza y covarianza de los residuales

library(biotools)

# Extraer los residuales del MANOVA
residuales_manova <- residuals(modelo_manova)

# Asignar nombres a las columnas
colnames(residuales_manova) <- c(
  "Grasas",
  "Humedad",
  "Proteina",
  "Cenizas"
)

# Comparar las matrices de varianzas y covarianzas
# de los residuales entre tratamientos
resultado_boxm_residuales <- boxM(
  residuales_manova,
  df$Tto
)

resultado_boxm_residuales

    Box's M-test for Homogeneity of Covariance Matrices

data:  residuales_manova
Chi-Sq (approx.) = 77.071, df = 60, p-value = 0.06803

No hay evidencia estadística suficiente para afirmar que las matrices de varianzas y covarianzas de los residuales sean diferentes entre los tratamientos.

Como conclusión, para estos datos no es recomendable realizar un MANOVA, dado que no se cumple el supuesto de normalidad multivariante, por lo cual se llevará a cabo la prueba de PERMANOVA.

PERMANOVA

De acuerdo con el articulo, para estos datos no todas las distancias tienen la misma pertinencia. La propuesta principal se realiza con la distancia euclídea sobre variables estandarizadas. Las demás servirán como análisis de sensibilidad para comprobar si la conclusión depende de la distancia elegida.

# =========================================================
# 1. PAQUETES
# =========================================================

library(vegan)
library(cluster)

# Tratamiento como factor
df$Tto <- as.factor(df$Tto)

# Seleccionar solamente las variables respuesta
Y <- df[, c(
  "Grasas",
  "Humedad",
  "Proteina",
  "Cenizas"
)]

# Estandarizar las variables
Y_z <- scale(Y)

# Comprobar que se creó correctamente
head(Y_z)

         Grasas    Humedad   Proteina    Cenizas
[1,] -0.2672462 -1.8028943 -1.3113283 -0.9431330
[2,]  0.4188379  0.1580559 -0.3687332 -1.8778252
[3,]  1.2332782 -1.4831765  0.3587196 -0.9578057
[4,]  0.3873916  0.5607460  0.1407275 -0.8049402
[5,] -0.9115641  0.7581751  0.4967772 -1.1296923
[6,] -0.4807853  0.4797777  0.4647485 -1.2174144

# Revisar medias y desviaciones estándar
round(colMeans(Y_z), 6)

  Grasas  Humedad Proteina  Cenizas 
       0        0        0        0 

apply(Y_z, 2, sd)

  Grasas  Humedad Proteina  Cenizas 
       1        1        1        1 

# Comprobar estructura
str(Y)

tibble [63 × 4] (S3: tbl_df/tbl/data.frame)
 $ Grasas  : num [1:63] 115 119 123 119 111 ...
 $ Humedad : num [1:63] 75.1 76.9 75.4 77.2 77.4 ...
 $ Proteina: num [1:63] 19 19.5 20 19.8 20 ...
 $ Cenizas : num [1:63] 1.08 1.03 1.08 1.09 1.07 ...

summary(Y)

     Grasas         Humedad         Proteina        Cenizas     
 Min.   :106.4   Min.   :74.33   Min.   :17.89   Min.   :1.025  
 1st Qu.:113.3   1st Qu.:76.13   1st Qu.:19.51   1st Qu.:1.099  
 Median :115.1   Median :76.77   Median :19.83   Median :1.135  
 Mean   :116.4   Mean   :76.73   Mean   :19.75   Mean   :1.141  
 3rd Qu.:119.1   3rd Qu.:77.29   3rd Qu.:20.10   3rd Qu.:1.172  
 Max.   :127.6   Max.   :78.73   Max.   :20.96   Max.   :1.348  

# Comprobar valores faltantes
colSums(is.na(Y))

  Grasas  Humedad Proteina  Cenizas 
       0        0        0        0 

# Conservar solamente filas completas
filas_completas <- complete.cases(Y, df$Tto)

Y <- Y[filas_completas, ]
grupo <- droplevels(df$Tto[filas_completas])

# Revisar nuevamente el balance
table(grupo)

grupo
T1 T2 T3 T4 T5 T6 T7 
 9  9  9  9  9  9  9 

La matriz multivariada contiene 63 observaciones completas y cuatro variables continuas. Presenta diferencias importantes de escala entre las variables, se recomienda usar distancia euclídea sobre datos estandarizados como análisis principal.

Construcción de distancias

# =========================================================
# 3. MATRICES DE DISTANCIA
# =========================================================

# Euclidea sin estandarizar
D_euclid_raw <- dist(
  Y,
  method = "euclidean"
)

# Euclidea estandarizada
D_euclid_z <- dist(
  Y_z,
  method = "euclidean"
)

# Manhattan estandarizada
D_manhattan_z <- dist(
  Y_z,
  method = "manhattan"
)

# Gower
D_gower <- cluster::daisy(
  Y,
  metric = "gower"
)

# Bray-Curtis
D_bray <- vegan::vegdist(
  Y,
  method = "bray"
)

# Hellinger + euclidea
Y_hellinger <- vegan::decostand(
  Y,
  method = "hellinger"
)

D_hellinger <- dist(
  Y_hellinger,
  method = "euclidean"
)

# =========================================================
# 4. PERMANOVA GLOBAL
# =========================================================

set.seed(123)

permanova_euclid_raw <- adonis2(
  D_euclid_raw ~ grupo,
  permutations = 9999
)

set.seed(123)

permanova_euclid_z <- adonis2(
  D_euclid_z ~ grupo,
  permutations = 9999
)

set.seed(123)

permanova_manhattan_z <- adonis2(
  D_manhattan_z ~ grupo,
  permutations = 9999
)

set.seed(123)

permanova_gower <- adonis2(
  D_gower ~ grupo,
  permutations = 9999
)

set.seed(123)

permanova_bray <- adonis2(
  D_bray ~ grupo,
  permutations = 9999
)

set.seed(123)

permanova_hellinger <- adonis2(
  D_hellinger ~ grupo,
  permutations = 9999
)

permanova_euclid_raw

Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 9999

adonis2(formula = D_euclid_raw ~ grupo, permutations = 9999)
         Df SumOfSqs     R2      F Pr(>F)
Model     6    63.33 0.0298 0.2867 0.9478
Residual 56  2061.94 0.9702              
Total    62  2125.27 1.0000              

permanova_euclid_z

Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 9999

adonis2(formula = D_euclid_z ~ grupo, permutations = 9999)
         Df SumOfSqs      R2      F Pr(>F)    
Model     6   48.093 0.19392 2.2454  7e-04 ***
Residual 56  199.907 0.80608                  
Total    62  248.000 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

permanova_manhattan_z

Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 9999

adonis2(formula = D_manhattan_z ~ grupo, permutations = 9999)
         Df SumOfSqs      R2      F Pr(>F)   
Model     6   139.48 0.19672 2.2858 0.0013 **
Residual 56   569.54 0.80328                 
Total    62   709.02 1.00000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

permanova_gower

Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 9999

adonis2(formula = D_gower ~ grupo, permutations = 9999)
         Df SumOfSqs      R2      F Pr(>F)   
Model     6  0.36244 0.17779 2.0182 0.0077 **
Residual 56  1.67610 0.82221                 
Total    62  2.03854 1.00000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

permanova_bray

Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 9999

adonis2(formula = D_bray ~ grupo, permutations = 9999)
         Df  SumOfSqs     R2     F Pr(>F)
Model     6 0.0008041 0.0521 0.513 0.8528
Residual 56 0.0146306 0.9479             
Total    62 0.0154348 1.0000             

permanova_hellinger

Permutation test for adonis under reduced model
Permutation: free
Number of permutations: 9999

adonis2(formula = D_hellinger ~ grupo, permutations = 9999)
         Df  SumOfSqs      R2      F Pr(>F)
Model     6 0.0007867 0.06636 0.6634 0.7149
Residual 56 0.0110682 0.93364              
Total    62 0.0118549 1.00000              

extraer_permanova <- function(modelo, nombre_distancia) {

  data.frame(
    Distancia = nombre_distancia,
    Df = modelo$Df[1],
    SumSq = modelo$SumOfSqs[1],
    R2 = modelo$R2[1],
    Pseudo_F = modelo$F[1],
    p_value = modelo$`Pr(>F)`[1]
  )
}

tabla_permanova <- rbind(
  extraer_permanova(
    permanova_euclid_raw,
    "Euclidea sin estandarizar"
  ),
  extraer_permanova(
    permanova_euclid_z,
    "Euclidea estandarizada"
  ),
  extraer_permanova(
    permanova_manhattan_z,
    "Manhattan estandarizada"
  ),
  extraer_permanova(
    permanova_gower,
    "Gower"
  ),
  extraer_permanova(
    permanova_bray,
    "Bray-Curtis"
  ),
  extraer_permanova(
    permanova_hellinger,
    "Hellinger + Euclidea"
  )
)

tabla_permanova[
  order(tabla_permanova$p_value),
]

Evaluación de la dispersión multivariada PERMDISP

# =========================================================
# 6. PERMDISP
# =========================================================

evaluar_permdisp <- function(distancia, grupo, semillas = 123) {

  dispersion <- betadisper(
    distancia,
    grupo,
    type = "centroid"
  )

  set.seed(semillas)

  prueba <- permutest(
    dispersion,
    permutations = 9999
  )

  list(
    objeto = dispersion,
    resultado = prueba
  )
}

disp_euclid_raw <- evaluar_permdisp(
  D_euclid_raw,
  grupo
)

disp_euclid_z <- evaluar_permdisp(
  D_euclid_z,
  grupo
)

disp_manhattan_z <- evaluar_permdisp(
  D_manhattan_z,
  grupo
)

disp_gower <- evaluar_permdisp(
  D_gower,
  grupo
)

disp_bray <- evaluar_permdisp(
  D_bray,
  grupo
)

disp_hellinger <- evaluar_permdisp(
  D_hellinger,
  grupo
)

disp_euclid_raw$resultado

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999

Response: Distances
          Df Sum Sq Mean Sq      F N.Perm Pr(>F)
Groups     6  37.42  6.2365 0.5993   9999 0.7321
Residuals 56 582.80 10.4071                     

disp_euclid_z$resultado

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)   
Groups     6  6.2047 1.03412 3.1765   9999 0.0095 **
Residuals 56 18.2307 0.32555                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

disp_manhattan_z$resultado

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999

Response: Distances
          Df Sum Sq Mean Sq      F N.Perm Pr(>F)  
Groups     6 18.760  3.1266 2.8329   9999 0.0166 *
Residuals 56 61.807  1.1037                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

disp_gower$resultado

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     6 0.048486 0.008081 2.5081   9999 0.0311 *
Residuals 56 0.180431 0.003222                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

disp_bray$resultado

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999

Response: Distances
          Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     6 0.0002356 3.9275e-05 0.6729   9999  0.676
Residuals 56 0.0032685 5.8366e-05                     

disp_hellinger$resultado

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 9999

Response: Distances
          Df     Sum Sq    Mean Sq      F N.Perm Pr(>F)
Groups     6 0.00033987 5.6646e-05 1.2224   9999 0.3104
Residuals 56 0.00259510 4.6341e-05                     

extraer_permdisp <- function(objeto, nombre_distancia) {

  tabla <- objeto$resultado$tab

  data.frame(
    Distancia = nombre_distancia,
    F_dispersion = tabla$F[1],
    p_dispersion = tabla$`Pr(>F)`[1]
  )
}

tabla_permdisp <- rbind(
  extraer_permdisp(
    disp_euclid_raw,
    "Euclidea sin estandarizar"
  ),
  extraer_permdisp(
    disp_euclid_z,
    "Euclidea estandarizada"
  ),
  extraer_permdisp(
    disp_manhattan_z,
    "Manhattan estandarizada"
  ),
  extraer_permdisp(
    disp_gower,
    "Gower"
  ),
  extraer_permdisp(
    disp_bray,
    "Bray-Curtis"
  ),
  extraer_permdisp(
    disp_hellinger,
    "Hellinger + Euclidea"
  )
)

tabla_permdisp

tabla_completa <- merge(
  tabla_permanova,
  tabla_permdisp,
  by = "Distancia"
)

tabla_completa[
  order(tabla_completa$p_value),
]

• Con base en estos resultados, se detectaron diferencias multivariadas entre los tratamientos con distancia Euclídea estandarizada (PERMANOVA), pero también se encontraron diferencias en la dispersión interna de los grupos (PERMDISP). Por ello, la significancia del PERMANOVA no puede atribuirse exclusivamente a diferencias entre centroides. Estos resultados se observaron en la distancia de Manhattan estandarizada y Gower, siendo mas sensibles en la detección de las diferencias en comparación con el resto.

• La Euclidea sin estandarizar mostró que las direncias de escala de las variables evaluadas fueron determinantes en los resultados obtenidos, ya que al elevar las diferencias al cuadrado, da mayor peso a unas variables con respecto a otras.

• Manhattan estandarizada suma el valor absoluto de las diferencias multivariadas siendo altamente sensible a esta prueba a pesar de que se emplea principalmente en datos de clasificación, agrupamiento y análisis de datos con posibles valores extremos.

• Gower es la distancia característica cuando se mezclan variables numéricas, ordinales y binarias, fue funcional para estos datos, debido a su ponderación el rango, lo que hace que la distancia maneje una escala similar entre variables.

• Hellinger Euclidea no detectó diferencias porque transformó los datos en perfiles relativos por muestra, eliminando una gran parte de las diferencias absolutas que sí eran relevantes en el análisis estandarizado.

• Bray-Curtis no detectó diferencias porque responde a diferencias relativas en perfiles de abundancia y no se ajusta bien a cuatro variables químicas con escalas y naturalezas diferentes.