PESO

ggbetweenstats(
  data = datos,
  x = Variedad,
  y = Peso,
  type = "p",
  var.equal = TRUE,
  title = "Comparación de Peso entre Variedades",
  xlab = "Variedad",
  ylab = "Peso (g)"
)

El gráfico muestra la distribución del peso de las manzanas según la variedad. Cada caja representa la mediana y el rango intercuartílico, mientras que los puntos indican los valores individuales. Visualmente, se observa que una variedad tiende a tener pesos mayores que la otra.

El efecto de la variedad sobre el peso, medido con Hedges’ g, es grande y negativo (-1.14), lo que indica que la primera variedad pesa significativamente menos que la segunda. El intervalo de confianza al 95% [-1.87, -0.40] refuerza que esta diferencia es estadísticamente significativa. El número de observaciones total es 32.

Conclusión: En conjunto, esto sugiere que la variedad es un factor importante que influye en el peso de las manzanas en esta localidad de Oxapampa.

DIÁMETRO

ggbetweenstats(
  data = datos,
  x = Variedad,
  y = Diametro,
  type = "p",
  var.equal = TRUE,
  title = "Comparación de Diámetro entre Variedades",
  xlab = "Variedad",
  ylab = "Diámetro (mm)")

El gráfico muestra la distribución del diámetro de las manzanas según la variedad. A diferencia del peso, las cajas y la dispersión de los datos indican que las diferencias en diámetro entre variedades son más pequeñas.

El efecto de la variedad sobre el diámetro, medido con Hedges’ g, es pequeño (0.37) y su intervalo de confianza al 95% [-0.32, 1.06] incluye cero. Esto indica que, aunque puede existir una ligera diferencia, no es estadísticamente significativa en esta muestra.

En conclusión, la variedad parece influir más sobre el peso que sobre el diámetro de las manzanas en Oxapampa.

Distancia de Mahalanobis

Calcular las distancias de Mahalanobis dentro de cada grupo usando las dos variables.

# Filtrar variables de interés (Peso y Diámetro)
grupoA <- datos %>% filter(Variedad == "A") %>% dplyr::select(Peso, Diametro)
grupoB <- datos %>% filter(Variedad == "B") %>% dplyr::select(Peso, Diametro)

# Calcular distancia de Mahalanobis para cada grupo
mahalA <- mahalanobis(grupoA,
                      colMeans(grupoA),
                      cov(grupoA))
mahalB <- mahalanobis(grupoB,
                      colMeans(grupoB),
                      cov(grupoB))

# Umbral para detectar outliers (p=2, alfa=0.975)
umbral <- qchisq(0.975, df = 2)

# Identificar posibles outliers
outliersA <- which(mahalA > umbral)
outliersB <- which(mahalB > umbral)

# Mostrar resultados
cat("Grupo A - posibles outliers:", outliersA, "\n")

## Grupo A - posibles outliers:

cat("Grupo B - posibles outliers:", outliersB, "\n")

## Grupo B - posibles outliers:

Objetivo: Detectar observaciones atípicas multivariadas que puedan influir en el análisis.

Método: Se calcularon las distancias de Mahalanobis de cada observación respecto a la media multivariada de su grupo, utilizando la matriz de covarianza muestral. Luego se compararon estas distancias con los percentiles de una distribución χ2 con𝑝=2 grados de libertad (por las dos variables analizadas).

Interpretación: Las observaciones que superan el umbral 𝜒2(2,0.975) se consideran posibles outliers. En el Q-Q Plot, estos puntos aparecen alejados de la línea teórica (roja). Si existen, podrían afectar la prueba de Hotelling T² y se recomienda analizarlas antes de continuar.

# Graficar Q-Q Plot de Chi-cuadrado para Grupo A
dfA <- data.frame(Distancia = sort(mahalA),
                  Cuantiles = qchisq(ppoints(length(mahalA)), df = 2))
ggplot(dfA, aes(x = Cuantiles, y = Distancia)) +
  geom_point(color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  geom_hline(yintercept = umbral, color = "tomato", linetype = "dotted") +
  labs(title = "Q-Q Plot Chi-cuadrado - Grupo A",
       x = expression(paste(chi^2, " teórico")),
       y = "Distancia de Mahalanobis") +
  theme_minimal()

# Graficar Q-Q Plot para Grupo B (igual procedimiento)
dfB <- data.frame(Distancia = sort(mahalB),
                  Cuantiles = qchisq(ppoints(length(mahalB)), df = 2))
ggplot(dfB, aes(x = Cuantiles, y = Distancia)) +
  geom_point(color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  geom_hline(yintercept = umbral, color = "tomato", linetype = "dotted") +
  labs(title = "Q-Q Plot Chi-cuadrado - Grupo B",
       x = expression(paste(chi^2, " teórico")),
       y = "Distancia de Mahalanobis") +
  theme_minimal()

Evaluación de Outliers Multivariados:

• Se calcularon las distancias de Mahalanobis para cada observación en los grupos A y B considerando las variables Peso y Diámetro. Estas distancias se compararon con el valor crítico de la distribución 𝜒2(2,0.975)=7.3778.

• Ninguna observación superó este umbral en ninguno de los grupos, por lo que no se detectaron outliers multivariados. Esto permite continuar con el análisis de Hotelling T² sin necesidad de depuración de datos.

Prueba de Hotelling T² para muestras independientes

La prueba de Hotelling T² compara de manera multivariante los vectores de medias entre dos grupos independientes.

• Hipótesis nula (H₀): 𝜇𝑋=𝜇𝑌

Es decir, los vectores de medias de Peso y Diámetro son iguales en las dos variedades. No hay diferencias multivariadas entre A y B.

• Hipótesis alternativa (H₁): 𝜇𝑋≠𝜇𝑌

Al menos una de las medias (Peso o Diámetro) es diferente o la combinación de ambas medias es distinta entre variedades.

# Definir X como la matriz de variables (Peso y Diámetro) para la variedad A.
X <- datos %>% filter(Variedad=="A") %>% dplyr::select(Peso, Diametro)

# Definir Y como la matriz para la variedad B.
Y <- datos %>% filter(Variedad=="B") %>% dplyr::select(Peso, Diametro)

hot <- HotellingsT2(X, Y)
hot

## 
##  Hotelling's two sample T2-test
## 
## data:  X and Y
## T.2 = 5.3374, df1 = 2, df2 = 29, p-value = 0.01062
## alternative hypothesis: true location difference is not equal to c(0,0)

• Se aplicó la prueba de Hotelling T² para comparar simultáneamente las medias de Peso y Diámetro entre las variedades A y B.

• El estadístico obtenido fue T² = 5.34, con df1=2 y df2=29, y un p-valor = 0.0106.

  • Como el p-valor es menor a 0.05, se rechaza la hipótesis nula de igualdad de medias multivariadas, concluyendo que existe diferencia significativa en el conjunto de variables Peso y Diámetro entre ambas variedades.

Conclusión: Este resultado confirma que las variedades presentan diferencias no solo en el peso, sino también en el patrón combinado de peso y diámetro.

Comparación con vector de diferencias esperado

Se tiene un vector de diferencias esperado (5,−2), podemos evaluar si las diferencias observadas coinciden con ese vector.

• Hipótesis nula (H₀): 𝜇𝑌−𝜇𝑋=𝜇0

La diferencia real de medias entre las variedades es exactamente igual a (5, −2). Es decir, se espera que B tenga 5 g más de peso y 2 mm menos de diámetro que A.

• Hipótesis alternativa (H₁): 𝜇𝑌−𝜇𝑋≠𝜇0

La diferencia observada no coincide con el vector hipotético, al menos en una de las variables.

mu0 <- c(5, -2)  # B - A

hot_mu <- HotellingsT2(X, Y, mu = mu0)
hot_mu

## 
##  Hotelling's two sample T2-test
## 
## data:  X and Y
## T.2 = 22.127, df1 = 2, df2 = 29, p-value = 1.461e-06
## alternative hypothesis: true location difference is not equal to c(5,-2)

Prueba de Hotelling T² con vector de medias hipotético:

Además de verificar la igualdad de medias, se evaluó si la diferencia de medias entre las variedades era igual a mu_o = (5,-2), donde se esperaba que la variedad B presentara 5 g más de peso y 2 mm menos de diámetro en comparación con la variedad A.

El estadístico obtenido fue T²= 22.13, con un p-valor = 1.46 × 10⁻⁶, por lo que se rechaza la hipótesis nula.

SIMULACIÓN

set.seed(123)

# Marcas y repeticiones
marcas <- rep(c("Toyota", "Ford", "Hyundai", "Chevrolet", "Volkswagen"), each = 70)

# Variables dependientes simuladas
velocidad   <- rnorm(350, mean = rep(c(190, 200, 185, 195, 188), each = 70), sd = 10)
consumo     <- rnorm(350, mean = rep(c(16, 14, 17, 15, 16.5), each = 70), sd = 1.5)
aceleracion <- rnorm(350, mean = rep(c(9.5, 8.8, 10.2, 9.3, 9.0), each = 70), sd = 0.6)
precio      <- rnorm(350, mean = rep(c(30, 28, 25, 27, 29), each = 70), sd = 3)
emisiones   <- rnorm(350, mean = rep(c(120, 130, 115, 125, 118), each = 70), sd = 10)

# Data frame final
autos <- data.frame(
  Marca = marcas,
  Velocidad_max = round(velocidad, 1),
  Consumo = round(consumo, 2),
  Aceleracion = round(aceleracion, 2),
  Precio = round(precio, 1),
  Emisiones = round(emisiones, 1)
)

3.1 Supuesto de normalidad multivariada

Hipótesis de Shapiro-Wilk multivariado

• H₀: Los datos provienen de una distribución normal multivariada.

• H₁: Los datos NO provienen de una normal multivariada.

library(tidyverse)
toyota = autos %>% filter(Marca == "Toyota") %>%
  dplyr::select(Velocidad_max, Consumo, Aceleracion, Precio, Emisiones)
ford = autos %>% filter(Marca == "Ford") %>%
  dplyr::select(Velocidad_max, Consumo, Aceleracion, Precio, Emisiones)
hyundai = autos %>% filter(Marca == "Hyundai") %>%
  dplyr::select(Velocidad_max, Consumo, Aceleracion, Precio, Emisiones)
chevrolet = autos %>% filter(Marca == "Chevrolet") %>%
  dplyr::select(Velocidad_max, Consumo, Aceleracion, Precio, Emisiones)
volkswagen = autos %>% filter(Marca == "Volkswagen") %>%
  dplyr::select(Velocidad_max, Consumo, Aceleracion, Precio, Emisiones)

#Ejecutamos el test de normalidad multivariada

library(mvnormtest)
mshapiro.test(t(toyota))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.96508, p-value = 0.04803

mshapiro.test(t(ford))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.92552, p-value = 0.0004645

mshapiro.test(t(hyundai))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.94714, p-value = 0.005153

mshapiro.test(t(chevrolet))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.96199, p-value = 0.03229

mshapiro.test(t(volkswagen))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.9522, p-value = 0.009468

conclusion:

Como en todas las marcas p< 0.05, rechazamos H₀. A un nivel de significación del 0.5%, existe suficiente evidencia estadística para afirmar que los datos no cumplen el supuesto de normalidad multivariada.

Sin embargo, debido a que contamos con un tamaño de muestra grande (n = 350, 70 por grupo), el Teorema Central del Límite (TCL) nos respalda: al aumentar el tamaño muestral.

3.2 Supuesto de homogeneidad de matrices variancia covariancia

• H₀: Las matrices de varianzas-covarianzas son iguales entre los grupos (marcas).

• H₁: Al menos una matriz de varianzas-covarianzas difiere entre los grupos.

En este caso se evaluaron 3 pruebas para evaluar este supuesto:

# MANOVA con 5 respuestas: Velocidad, Consumo, Aceleración, Precio, Emisiones
res <- boxM(autos[, -1], autos$Marca)
res

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  autos[, -1]
## Chi-Sq (approx.) = 36.548, df = 60, p-value = 0.9928

summary(res)

## Summary for Box's M-test of Equality of Covariance Matrices
## 
## Chi-Sq:   36.54757 
## df:   60 
## p-value: 0.9928 
## 
## log of Covariance determinants:
##  Chevrolet       Ford    Hyundai     Toyota Volkswagen 
##   11.13877   10.72263   10.93171   10.75061   11.48488 
## 
## Eigenvalues:
##     Chevrolet       Ford     Hyundai     Toyota  Volkswagen      pooled
## 1 123.5055709 97.9369713 109.3741906 91.5285303 121.6573435 102.8216093
## 2  80.6349734 79.4242280  89.7329031 77.2978375 105.1141321  92.2590158
## 3   8.7326293  7.7449957   9.2115855  8.5522379   8.1838240   8.6009040
## 4   2.4113255  2.2995390   1.6776331  2.1567286   2.2290557   2.1821441
## 5   0.3280159  0.3275034   0.3687094  0.3575441   0.4168155   0.3772915
## 
## Statistics based on eigenvalues:
##              Chevrolet         Ford      Hyundai       Toyota   Volkswagen
## product   6.878691e+04 4.537090e+04 55921.925473 4.665829e+04 9.723453e+04
## sum       2.156125e+02 1.877332e+02   210.365022 1.798929e+02 2.376012e+02
## precision 2.779049e-01 2.747109e-01     0.290944 2.940041e-01 3.347068e-01
## max       1.235056e+02 9.793697e+01   109.374191 9.152853e+01 1.216573e+02
##                 pooled
## product   6.717346e+04
## sum       2.062410e+02
## precision 3.081126e-01
## max       1.028216e+02

Conclución: Como el p-valor=0.9928 > 0.05, no se rechaza Ho, existe suficiente evidencia estadística para afirmar que se cumple el supuesto de homogeneidad de las matrices de varianzas-covarianzas entre las cinco marcas de vehículos.

library(biotools)

boxM(autos[, -1], autos[, 1])

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  autos[, -1]
## Chi-Sq (approx.) = 36.548, df = 60, p-value = 0.9928

Con la anterior prueba tambien cumple el supuestro de homogeneidad de varianzas.

Y con esta prueba de WRAPER

library(DFA.CANCOR)

## **************************************************************************************************
## DFA.CANCOR 0.3.6
## 
## Please contact Brian O'Connor at brian.oconnor@ubc.ca if you have questions or suggestions.
## **************************************************************************************************

HOMOGENEITY(data = autos,
            groups = 'Marca',
            variables = c('Velocidad_max','Consumo','Aceleracion','Precio','Emisiones'))

## 
## Covariance matrix for GroupToyota

##               Velocidad_max Consumo Aceleracion Precio Emisiones
## Velocidad_max         84.58   -0.17       -0.04  -3.02    -12.56
## Consumo               -0.17    2.42       -0.37   0.31      2.53
## Aceleracion           -0.04   -0.37        0.40  -0.01     -1.20
## Precio                -3.02    0.31       -0.01   8.84      0.27
## Emisiones            -12.56    2.53       -1.20   0.27    119.37

## 
## Covariance matrix for GroupFord

##               Velocidad_max Consumo Aceleracion Precio Emisiones
## Velocidad_max         87.25   -0.08        0.29   4.52      9.23
## Consumo               -0.08    2.36       -0.06   0.44     -1.74
## Aceleracion            0.29   -0.06        0.34   0.24      0.09
## Precio                 4.52    0.44        0.24   7.98     -0.58
## Emisiones              9.23   -1.74        0.09  -0.58     89.81

## 
## Covariance matrix for GroupHyundai

##               Velocidad_max Consumo Aceleracion Precio Emisiones
## Velocidad_max        102.05   -0.74       -0.73  -0.86     -9.52
## Consumo               -0.74    1.81        0.06   0.71     -2.29
## Aceleracion           -0.73    0.06        0.39   0.25     -0.76
## Precio                -0.86    0.71        0.25   9.15     -0.42
## Emisiones             -9.52   -2.29       -0.76  -0.42     96.97

## 
## Covariance matrix for GroupChevrolet

##               Velocidad_max Consumo Aceleracion Precio Emisiones
## Velocidad_max         82.52   -1.91        1.10  -2.82      6.78
## Consumo               -1.91    2.24       -0.11   0.47     -1.47
## Aceleracion            1.10   -0.11        0.38   0.00     -0.30
## Precio                -2.82    0.47        0.00   8.71     -2.58
## Emisiones              6.78   -1.47       -0.30  -2.58     86.05

## 
## Covariance matrix for GroupVolkswagen

##               Velocidad_max Consumo Aceleracion Precio Emisiones
## Velocidad_max        105.23   -0.81        1.42  -0.51      1.53
## Consumo               -0.81    2.26       -0.04   0.34     -0.88
## Aceleracion            1.42   -0.04        0.44  -0.11     -0.03
## Precio                -0.51    0.34       -0.11   8.19      1.72
## Emisiones              1.53   -0.88       -0.03   1.72    121.48

## 
## 
## Bartlett test of HOMOGENEITY of variances (parametric):

## 
## Bartlett's K-squared =2446.808  df =4  p value =0

## 
## 
## Fligner-Killeen test of HOMOGENEITY of variances (non parametric):

## 
## Fligner-Killeen chi-squared =869.084  df =4  p value =0

## 
## 
## Pooled within groups covariance matrix from SPSS:

##               Velocidad_max Consumo Aceleracion Precio Emisiones
## Velocidad_max        92.325  -0.741       0.409 -0.540    -0.909
## Consumo              -0.741   2.219      -0.104  0.457    -0.768
## Aceleracion           0.409  -0.104       0.388  0.075    -0.441
## Precio               -0.540   0.457       0.075  8.573    -0.317
## Emisiones            -0.909  -0.768      -0.441 -0.317   102.735

## 
## 
## Pooled within groups correlation matrix from SPSS:

##               Velocidad_max Consumo Aceleracion Precio Emisiones
## Velocidad_max         1.000  -0.052       0.068 -0.019    -0.009
## Consumo              -0.052   1.000      -0.112  0.105    -0.051
## Aceleracion           0.068  -0.112       1.000  0.041    -0.070
## Precio               -0.019   0.105       0.041  1.000    -0.011
## Emisiones            -0.009  -0.051      -0.070 -0.011     1.000

## 
## 
## Box Test of equality of covariance matrices:

## 
## Log determinants:

##            Log Determinant
## Toyota              11.139
## Ford                10.723
## Hyundai             10.932
## Chevrolet           10.751
## Volkswagen          11.485
## Pooled              11.115

## 
## 
## M = 37.714  F = 0.609  df1 = 60  df2 = 238906.29  p = 0.99278

3.3 Supuesto de variables dependientes correlacionadas (Bartlett)

• H₀: Las variables dependientes no están correlacionadas (matriz identidad)

• H₁: Las variables dependientes están correlacionadas

library(psych)

## 
## Adjuntando el paquete: 'psych'

## The following object is masked from 'package:car':
## 
##     logit

## The following object is masked from 'package:effectsize':
## 
##     phi

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

options(scipen = 0)

# Bartlett sobre la matriz de correlación de las variables dependientes
cortest.bartlett(cor(autos[, -1]), n = nrow(autos[, -1]))

## $chisq
## [1] 127.8793
## 
## $p.value
## [1] 1.264207e-22
## 
## $df
## [1] 10

Conclución:Dado que el p-valor < 0.05, se rechaza la hipótesis nula. Por lo tanto, existe evidencia estadística significativa para afirmar que las variables dependientes están correlacionadas entre sí, lo que valida el supuesto de correlación necesario para aplicar un análisis MANOVA.

library(MVTests)

res1 <- Bsper(autos[, -1])  # Variables dependientes
summary(res1)

##        Bartlett's Sphericity Test 
## 
## Chi-Squared Value = 127.8793 , df = 10  and p-value: <2e-16 
## 
##        Correlation Matrix 
## 
##               Velocidad_max      Consumo Aceleracion       Precio   Emisiones
## Velocidad_max    1.00000000 -0.293988633  -0.2101076  0.047949543  0.22523062
## Consumo         -0.29398863  1.000000000   0.2257921  0.001177129 -0.31740846
## Aceleracion     -0.21010757  0.225792120   1.0000000 -0.136911411 -0.29881962
## Precio           0.04794954  0.001177129  -0.1369114  1.000000000  0.04294112
## Emisiones        0.22523062 -0.317408461  -0.2988196  0.042941115  1.00000000

Además,la matriz de correlación muestra relaciones moderadas entre algunas variables, como por ejemplo:

• Velocidad máxima y Consumo: correlación negativa (-0.29)

• Consumo y Emisiones: correlación negativa (-0.32)

• Velocidad máxima y Emisiones: correlación positiva (0.23)

Estas correlaciones confirman que las variables no son independientes

4.1 Determinación de la matriz residual y la matriz factorial del MANOVA.

A continuación en la salida se muestran amabas matrices

Matrices = summary(modelo)$SS
F = Matrices$Marca       # Variabilidad explicada por el factor Marca
W = Matrices$Residuals   # Variabilidad residual (dentro de cada marca)
T = F + W                # Variabilidad total
F

##               Velocidad_max    Consumo Aceleracion    Precio  Emisiones
## Velocidad_max     8611.6367 -1719.9534   -794.1877  795.6666 10219.0869
## Consumo          -1719.9534   350.4071    152.3376 -155.0666 -2053.2561
## Aceleracion       -794.1877   152.3376    104.9375 -159.4726  -857.6446
## Precio             795.6666  -155.0666   -159.4726 1034.4562   702.6053
## Emisiones        10219.0869 -2053.2561   -857.6446  702.6053 12357.0519

W

##               Velocidad_max    Consumo Aceleracion     Precio  Emisiones
## Velocidad_max    31852.2023 -255.68800   141.01897 -186.23114  -313.5566
## Consumo           -255.6880  765.65421   -35.76323  157.55137  -265.0955
## Aceleracion        141.0190  -35.76323   133.89877   25.78241  -152.0174
## Precio            -186.2311  157.55137    25.78241 2957.80957  -109.4071
## Emisiones         -313.5566 -265.09546  -152.01743 -109.40714 35443.5680

5.1 Criterio de Razón de verosimilitud (Lambda de Wilk’s)

eta2 = 1 - det(W)/det(T)
eta2

## [1] 0.7692905

Conclusión: El 76.9% de la variabilidad total está explicada por las diferencias entre los grupos (Marca),indica efecto fuerte del factor Marca sobre las variables dependientes.

Lo que se quiere probar es lo siguiente:

H₀:μ1=μ2=μ3=μ4=μ5

H₁: al menos un μi es diferente a los demás

summary(modelo, test = "Pillai")

##            Df Pillai approx F num Df den Df    Pr(>F)    
## Marca       4 1.0182   23.492     20   1376 < 2.2e-16 ***
## Residuals 345                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo, test = "Hotelling-Lawley")

##            Df Hotelling-Lawley approx F num Df den Df    Pr(>F)    
## Marca       4           2.3386   39.699     20   1358 < 2.2e-16 ***
## Residuals 345                                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo, test = "Roy")

##            Df    Roy approx F num Df den Df    Pr(>F)    
## Marca       4 1.8748   128.99      5    344 < 2.2e-16 ***
## Residuals 345                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusión: Dado que todas las pruebas (Pillai, Hotelling, Roy) son significativas con p < 0.05, podemos afirmar que:Existe suficiente evidencia estadística para rechazar la hipótesis nula Ho, los vectores de medias de las variables dependientes difieren entre las marcas de autos.

Por lo tanto, el factor Marca tiene un efecto significativo sobre el conjunto de variables dependientes.

summary.aov(modelo)

##  Response Velocidad_max :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Marca         4   8612 2152.91  23.319 < 2.2e-16 ***
## Residuals   345  31852   92.33                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Consumo :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Marca         4 350.41  87.602  39.473 < 2.2e-16 ***
## Residuals   345 765.65   2.219                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Aceleracion :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Marca         4 104.94 26.2344  67.595 < 2.2e-16 ***
## Residuals   345 133.90  0.3881                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Precio :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Marca         4 1034.5 258.614  30.165 < 2.2e-16 ***
## Residuals   345 2957.8   8.573                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Emisiones :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Marca         4  12357 3089.26   30.07 < 2.2e-16 ***
## Residuals   345  35444  102.73                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusión: De la salida se puede apreciar que todas las variables aportan para que se rechace la hipotesis nula H₀, es decir que contribuyen para que se esta hipótesis se rechace, esto se debe que cada uno de los p-valores son menores al nivel de significación alfa.

5.2 Comparamos dos Marcas de automóviles

Primero generamos un modelo lineal

modelo1 = manova(cbind(Velocidad_max, Consumo, Aceleracion,Precio, Emisiones) ~ Marca, data = autos,
               subset = Marca %in% c("Toyota","Ford"))

summary(modelo1,test="Pillai")

##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.59459   39.305      5    134 < 2.2e-16 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo1,test="Wilks")

##            Df   Wilks approx F num Df den Df    Pr(>F)    
## Marca       1 0.40541   39.305      5    134 < 2.2e-16 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo1,test="Hotelling-Lawley")

##            Df Hotelling-Lawley approx F num Df den Df    Pr(>F)    
## Marca       1           1.4666   39.305      5    134 < 2.2e-16 ***
## Residuals 138                                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo1,test="Roy")

##            Df    Roy approx F num Df den Df    Pr(>F)    
## Marca       1 1.4666   39.305      5    134 < 2.2e-16 ***
## Residuals 138                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como en todos los casos el p-valor < 0.05, rechazamos H₀. Se puede afirmar que con las distintas pruebas de pillai,willks,hotelling y roy, en todas ellas rechazamos H₀, es decir podemos afirmar que el vector de promedios de las marcas “Toyota” y “Ford” no son iguales. Es decir las 5 variables dicen que hay diferencias entre los promedios de las marcas “Toyota” y “Ford”.

5.3 Comparación General por pares

# H0:los 2 vectores de promedios son iguales
# H1:los 2 vectores de promedios difieren

# Hipótesis:
# H0: los 2 vectores de promedios (de las variables respuesta) son iguales
# H1: los 2 vectores de promedios difieren

# Definir los grupos de comparación
Marcas <- c("Toyota", "Ford", "Hyundai", "Chevrolet", "Volkswagen")

# Crear todas las combinaciones posibles de pares de marcas
comb <- t(combn(length(Marcas), 2))

# Comparación general por pares
for(i in 1:nrow(comb)){
  modelo.comp = manova(cbind(Velocidad_max, Consumo, Aceleracion, Precio, Emisiones) ~ Marca, 
                       data = autos,
                       subset = Marca %in% Marcas[comb[i,]])
  
  cat("Comparación entre:", Marcas[comb[i,]][1], "y", Marcas[comb[i,]][2], "\n")
  print(summary(modelo.comp, test = "Pillai"))
  cat("\n-------------------------------------\n\n")
}

## Comparación entre: Toyota y Ford 
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.59459   39.305      5    134 < 2.2e-16 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------
## 
## Comparación entre: Toyota y Hyundai 
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.62557   44.775      5    134 < 2.2e-16 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------
## 
## Comparación entre: Toyota y Chevrolet 
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.37394   16.008      5    134 2.376e-12 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------
## 
## Comparación entre: Toyota y Volkswagen 
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.23155   8.0756      5    134 1.101e-06 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------
## 
## Comparación entre: Ford y Hyundai 
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.80654   111.73      5    134 < 2.2e-16 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------
## 
## Comparación entre: Ford y Chevrolet 
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.27808   10.323      5    134 2.173e-08 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------
## 
## Comparación entre: Ford y Volkswagen 
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.55347   33.219      5    134 < 2.2e-16 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------
## 
## Comparación entre: Hyundai y Chevrolet 
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.70518   64.104      5    134 < 2.2e-16 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------
## 
## Comparación entre: Hyundai y Volkswagen 
##            Df Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.6167   43.118      5    134 < 2.2e-16 ***
## Residuals 138                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------
## 
## Comparación entre: Chevrolet y Volkswagen 
##            Df  Pillai approx F num Df den Df    Pr(>F)    
## Marca       1 0.37657   16.188      5    134 1.811e-12 ***
## Residuals 138                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -------------------------------------

Interpretación: El test de Pillai muestra que todas las comparaciones arrojan p < 0.05, lo que implica diferencias estadísticamente significativas entre cada par de marcas.

Por lo tanto, podemos concluir:

Existe suficiente evidencia estadística para afirmar que los vectores de medias de Velocidad máxima, Consumo, Aceleración, Precio y Emisiones difieren entre todas las marcas comparadas. Ninguna pareja de marcas tiene vectores de medias iguales.

Contexto del estudio

La productividad de los trabajadores en entornos de oficina está determinada por diversos factores, como la modalidad de trabajo y las herramientas digitales empleadas para la gestión de tareas. En un contexto organizacional, resulta relevante identificar qué combinaciones de estas condiciones favorecen un mejor rendimiento laboral, optimizan el tiempo de respuesta y disminuyen la carga asociada al uso intensivo de pantallas.

Para responder a esta problemática, se implementó un diseño factorial 3×4 en bloques completos al azar (DBCA)

Los factores principales fueron:

• Estrategia de trabajo: Híbrido, Teletrabajo y Presencial.

• Herramienta digital: Google Drive, Microsoft Excel, Microsoft Teams y Microsoft Outlook.

Bloque:

• Género: Masculino y Femenino

Con el propósito de controlar la variabilidad atribuible a esta característica y obtener estimaciones más precisas de los efectos de los factores de interés.

Unidad experimental:

• Cada trabajador individual

Las variables de respuesta registradas fueron:

• Tareas completadas : número de actividades finalizadas por trabajador en la jornada laboral.

• Horas de pantalla (horas/día) : tiempo de exposición frente al computador.

• Tiempo de respuesta a tareas (minutos) : rapidez promedio en atender y resolver solicitudes laborales.

    #Cargar librería para leer archivos Excel
    library(readxl)
    # Importar los datos desde el archivo Excel "productividad_oficinas.xlsx"
    datos <- read_excel("productividad_oficinas_ajustada_v10.xlsx")
    datos

## # A tibble: 1,000 × 6
##    Estrategia_Trabajo Herramienta_Digital Tareas_Completadas Horas_Pantalla
##    <chr>              <chr>                            <dbl>          <dbl>
##  1 Híbrido            Google Drive                        11           8.07
##  2 Híbrido            Microsoft Excel                     10           9.19
##  3 Híbrido            Microsoft Teams                      9           9.64
##  4 Presencial         Microsoft Excel                      9           7.48
##  5 Híbrido            Google Drive                        12           7.91
##  6 Presencial         Google Drive                        10           7.27
##  7 Presencial         Google Drive                        11           7.68
##  8 Presencial         Microsoft Outlook                   10           7.40
##  9 Híbrido            Microsoft Teams                      9           8.00
## 10 Teletrabajo        Microsoft Outlook                    8           5.44
## # ℹ 990 more rows
## # ℹ 2 more variables: Tiempo_Respuesta_Tareas <dbl>, Genero <chr>

# Ver la estructura inicial de los datos
str(datos)

## tibble [1,000 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Estrategia_Trabajo     : chr [1:1000] "Híbrido" "Híbrido" "Híbrido" "Presencial" ...
##  $ Herramienta_Digital    : chr [1:1000] "Google Drive" "Microsoft Excel" "Microsoft Teams" "Microsoft Excel" ...
##  $ Tareas_Completadas     : num [1:1000] 11 10 9 9 12 10 11 10 9 8 ...
##  $ Horas_Pantalla         : num [1:1000] 8.07 9.19 9.64 7.48 7.91 ...
##  $ Tiempo_Respuesta_Tareas: num [1:1000] 41.8 43.1 40 39.5 42.2 ...
##  $ Genero                 : chr [1:1000] "Masculino" "Femenino" "Masculino" "Femenino" ...

# Convertir variables categóricas en factores
datos$Estrategia_Trabajo = as.factor(datos$Estrategia_Trabajo)
datos$Herramienta_Digital = as.factor(datos$Herramienta_Digital)
datos$Genero = as.factor(datos$Genero)
# Verificar nuevamente la estructura para confirmar los cambios
str(datos)

## tibble [1,000 × 6] (S3: tbl_df/tbl/data.frame)
##  $ Estrategia_Trabajo     : Factor w/ 3 levels "Híbrido","Presencial",..: 1 1 1 2 1 2 2 2 1 3 ...
##  $ Herramienta_Digital    : Factor w/ 4 levels "Google Drive",..: 1 2 4 2 1 1 1 3 4 3 ...
##  $ Tareas_Completadas     : num [1:1000] 11 10 9 9 12 10 11 10 9 8 ...
##  $ Horas_Pantalla         : num [1:1000] 8.07 9.19 9.64 7.48 7.91 ...
##  $ Tiempo_Respuesta_Tareas: num [1:1000] 41.8 43.1 40 39.5 42.2 ...
##  $ Genero                 : Factor w/ 2 levels "Femenino","Masculino": 2 1 2 1 2 2 2 1 2 1 ...

# Mostrar las primeras filas de la base de datos para inspeccionar su contenido
head(datos)

## # A tibble: 6 × 6
##   Estrategia_Trabajo Herramienta_Digital Tareas_Completadas Horas_Pantalla
##   <fct>              <fct>                            <dbl>          <dbl>
## 1 Híbrido            Google Drive                        11           8.07
## 2 Híbrido            Microsoft Excel                     10           9.19
## 3 Híbrido            Microsoft Teams                      9           9.64
## 4 Presencial         Microsoft Excel                      9           7.48
## 5 Híbrido            Google Drive                        12           7.91
## 6 Presencial         Google Drive                        10           7.27
## # ℹ 2 more variables: Tiempo_Respuesta_Tareas <dbl>, Genero <fct>

# Supuesto de normalidad multivariada

# Desagregamos la base por tratamientos
datosc = datos
datosc

## # A tibble: 1,000 × 6
##    Estrategia_Trabajo Herramienta_Digital Tareas_Completadas Horas_Pantalla
##    <fct>              <fct>                            <dbl>          <dbl>
##  1 Híbrido            Google Drive                        11           8.07
##  2 Híbrido            Microsoft Excel                     10           9.19
##  3 Híbrido            Microsoft Teams                      9           9.64
##  4 Presencial         Microsoft Excel                      9           7.48
##  5 Híbrido            Google Drive                        12           7.91
##  6 Presencial         Google Drive                        10           7.27
##  7 Presencial         Google Drive                        11           7.68
##  8 Presencial         Microsoft Outlook                   10           7.40
##  9 Híbrido            Microsoft Teams                      9           8.00
## 10 Teletrabajo        Microsoft Outlook                    8           5.44
## # ℹ 990 more rows
## # ℹ 2 more variables: Tiempo_Respuesta_Tareas <dbl>, Genero <fct>

Hipótesis:

Hipótesis de Shapiro-Wilk multivariado

• H₀: Los datos provienen de una distribución normal multivariada.

• H₁: Los datos NO provienen de una normal multivariada.

# Cargamos tidyverse para manipulación de datos
library(tidyverse)

# Tratamiento 1: Teletrabajo con Microsoft Outlook
trat1 = datosc %>% 
  filter(Estrategia_Trabajo == "Teletrabajo", Herramienta_Digital == "Microsoft Outlook") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 2: Teletrabajo con Microsoft Excel
trat2 = datosc %>% 
  filter(Estrategia_Trabajo == "Teletrabajo", Herramienta_Digital == "Microsoft Excel") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 3: Presencial con Microsoft Outlook
trat3 = datosc %>% 
  filter(Estrategia_Trabajo == "Presencial", Herramienta_Digital == "Microsoft Outlook") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 4: Presencial con Microsoft Excel
trat4 = datosc %>% 
  filter(Estrategia_Trabajo == "Presencial", Herramienta_Digital == "Microsoft Excel") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 5: Teletrabajo con Google Drive
trat5 = datosc %>% 
  filter(Estrategia_Trabajo == "Teletrabajo", Herramienta_Digital == "Google Drive") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 6: Teletrabajo con Microsoft Teams
trat6 = datosc %>% 
  filter(Estrategia_Trabajo == "Teletrabajo", Herramienta_Digital == "Microsoft Teams") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 7: Presencial con Google Drive
trat7 = datosc %>% 
  filter(Estrategia_Trabajo == "Presencial", Herramienta_Digital == "Google Drive") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 8: Presencial con Microsoft Teams
trat8 = datosc %>% 
  filter(Estrategia_Trabajo == "Presencial", Herramienta_Digital == "Microsoft Teams") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 9: Híbrido con Microsoft Outlook
trat9 = datosc %>% 
  filter(Estrategia_Trabajo == "Híbrido", Herramienta_Digital == "Microsoft Outlook") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 10: Híbrido con Microsoft Excel
trat10 = datosc %>% 
  filter(Estrategia_Trabajo == "Híbrido", Herramienta_Digital == "Microsoft Excel") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 11: Híbrido con Google Drive
trat11 = datosc %>% 
  filter(Estrategia_Trabajo == "Híbrido", Herramienta_Digital == "Google Drive") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Tratamiento 12: Híbrido con Microsoft Teams
trat12 = datosc %>% 
  filter(Estrategia_Trabajo == "Híbrido", Herramienta_Digital == "Microsoft Teams") %>%
  dplyr::select(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas)

# Usamos la prueba de Shapiro-Wilk multivariada (mshapiro.test)
# Se necesita transponer (t) porque la función espera variables en filas

library(mvnormtest)
mshapiro.test(t(trat1))  # normalidad para el tratamiento 1

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.97972, p-value = 0.2075

mshapiro.test(t(trat2))  # normalidad para el tratamiento 2

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.98594, p-value = 0.479

mshapiro.test(t(trat3))  # normalidad para el tratamiento 3

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.97575, p-value = 0.08065

mshapiro.test(t(trat4))  # normalidad para el tratamiento 4

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.98305, p-value = 0.4941

mshapiro.test(t(trat5))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.95902, p-value = 0.01713

mshapiro.test(t(trat6))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.98755, p-value = 0.6183

mshapiro.test(t(trat7))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.982, p-value = 0.3467

mshapiro.test(t(trat8))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.98784, p-value = 0.5417

mshapiro.test(t(trat9))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.97889, p-value = 0.2438

mshapiro.test(t(trat10))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.9527, p-value = 0.001342

mshapiro.test(t(trat11))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.99043, p-value = 0.7981

mshapiro.test(t(trat12))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.98435, p-value = 0.3945

Conclusión:

De los 12 tratamientos analizados, 10 presentan p > 0.05, por lo que no se rechaza H₀ de normalidad multivariada en la mayoría de los casos. Sin embargo, en 2 tratamientos el supuesto no se cumple (p ≤ 0.05).

Conclusión: A un nivel de significación del 5%, se puede afirmar que la mayor parte de los datos cumplen con el supuesto de normalidad multivariada, aunque existen algunas excepciones que deben considerarse.

NOTA: Esto sugiere que el uso de MANOVA es en general válido, pero se recomienda cautela en la interpretación de los resultados o considerar métodos robustos que toleren leves desviaciones de normalidad.

Hipótesis:

• H₀: Las matrices de varianzas-covarianzas son iguales entre los grupos (tratamientos).

• H₁: Al menos una matriz de varianzas-covarianzas difiere entre los grupos.

# Supuesto de homogeneidad de matrices varianza-covarianza
# Usamos la prueba de Box M (heplots::boxM)
# Compara si las matrices de varianza-covarianza son iguales entre grupos

library(heplots)
heplots::boxM(cbind(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas) ~ 
                Estrategia_Trabajo * Herramienta_Digital, data = datos)

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Y
## Chi-Sq (approx.) = 57.323, df = 66, p-value = 0.768

Conclusión:

Como el p-valor = 0.768 > 0.05, no se rechaza H₀. Existe suficiente evidencia estadística para afirmar que se cumple el supuesto de homogeneidad de las matrices de varianzas-covarianzas entre los 12 tratamientos (combinaciones de Estrategia de trabajo y Herramienta digital).

Hipótesis:

• H₀: Las variables dependientes no están correlacionadas(matriz identidad)

• H₁: Las variables dependientes están correlacionadas

# Supuesto de variables dependientes correlacionadas. Prueba de esfericidad de Bartlett
# Usamos la prueba de esfericidad de Bartlett
# Verifica si la matriz de correlaciones es distinta de la identidad
# (es decir, si las variables dependientes están correlacionadas)
datos1 = datos[3:5]  # seleccionamos las 3 variables dependientes
head(datos1)

## # A tibble: 6 × 3
##   Tareas_Completadas Horas_Pantalla Tiempo_Respuesta_Tareas
##                <dbl>          <dbl>                   <dbl>
## 1                 11           8.07                    41.8
## 2                 10           9.19                    43.1
## 3                  9           9.64                    40.0
## 4                  9           7.48                    39.5
## 5                 12           7.91                    42.2
## 6                 10           7.27                    39.2

library(psych)
options(scipen=0) #permitirá que aparezcan números en notación científica cuando lo considere conveniente
cortest.bartlett(cor(datos1), n = nrow(datos1))

## $chisq
## [1] 696.6452
## 
## $p.value
## [1] 1.120669e-150
## 
## $df
## [1] 3

Conclusión: Dado que el p-valor< 0.05, se rechaza la hipótesis nula. Por lo tanto, existe evidencia estadística significativa para afirmar que las variables dependientes están correlacionadas entre sí, lo que valida el supuesto de correlación necesario para aplicar un análisis MANOVA.

MANOVA en ANOVAs individuales

# Descompone el MANOVA en ANOVAs individuales, mostrando los efectos de los factores
# y sus interacciones sobre cada variable dependiente (Tareas_Completadas,
# Horas_Pantalla y Tiempo_Respuesta_Tareas).
summary.aov(modelo)

##  Response Tareas_Completadas :
##                                         Df Sum Sq Mean Sq F value Pr(>F)    
## Estrategia_Trabajo                       2 210.67 105.337 153.439 <2e-16 ***
## Herramienta_Digital                      3 198.69  66.230  96.473 <2e-16 ***
## Genero                                   1   1.25   1.249   1.820 0.1776    
## Estrategia_Trabajo:Herramienta_Digital   6 140.20  23.367  34.038 <2e-16 ***
## Residuals                              987 677.58   0.687                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Horas_Pantalla :
##                                         Df Sum Sq Mean Sq  F value  Pr(>F)    
## Estrategia_Trabajo                       2 138.24  69.119 106.6602 < 2e-16 ***
## Herramienta_Digital                      3 129.55  43.184  66.6387 < 2e-16 ***
## Genero                                   1   2.01   2.011   3.1025 0.07848 .  
## Estrategia_Trabajo:Herramienta_Digital   6 175.15  29.192  45.0478 < 2e-16 ***
## Residuals                              987 639.60   0.648                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Tiempo_Respuesta_Tareas :
##                                         Df Sum Sq Mean Sq  F value  Pr(>F)    
## Estrategia_Trabajo                       2 834.64  417.32 659.3301 < 2e-16 ***
## Herramienta_Digital                      3 747.87  249.29 393.8572 < 2e-16 ***
## Genero                                   1   3.58    3.58   5.6584 0.01756 *  
## Estrategia_Trabajo:Herramienta_Digital   6 707.15  117.86 186.2045 < 2e-16 ***
## Residuals                              987 624.72    0.63                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusión:

De la salida se puede apreciar que todas las variables aportan para que se rechace la hipotesis nula H₀, es decir que contribuyen para que esta hipótesis se rechace, esto se debe que cada uno de los p-valores son menores al nivel de significación alfa igual a 0.05

7.1 Efectos simples

Cuando la interacción SÍ es significativa)

# Modelo MANOVA completo con interacción Estrategia_Trabajo * Herramienta_Digital
# Prueba de hipótesis usando el estadístico de Roy 
modelo_manova <- manova(cbind(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas) ~ Estrategia_Trabajo * Herramienta_Digital, data = datos)
summary(modelo_manova, test = "Pillai")

##                                         Df  Pillai approx F num Df den Df
## Estrategia_Trabajo                       2 0.69572  175.492      6   1974
## Herramienta_Digital                      3 0.66095   93.060      9   2964
## Estrategia_Trabajo:Herramienta_Digital   6 0.68325   48.563     18   2964
## Residuals                              988                               
##                                           Pr(>F)    
## Estrategia_Trabajo                     < 2.2e-16 ***
## Herramienta_Digital                    < 2.2e-16 ***
## Estrategia_Trabajo:Herramienta_Digital < 2.2e-16 ***
## Residuals                                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusión:

Tanto la Estrategia de trabajo, la Herramienta digital y, sobre todo, su Interacción influyen significativamente en el desempeño de los trabajadores medido a través de las tres variables de respuesta.

library(dplyr)

# Para cada herramienta (Excel, Outlook, Teams, etc.):
#   1. Filtramos el dataset por la herramienta en cuestión.
#   2. Ajustamos un MANOVA considerando únicamente Estrategia_Trabajo como factor.
#   3. Evaluamos la significancia mediante el test de Roy.

levels(datos$Herramienta_Digital) %>% lapply(function(nivel) {
  cat("\n--- Análisis para Herramienta =", nivel, "---\n")
  subdatos <- filter(datos, Herramienta_Digital == nivel)
  summary(manova(cbind(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas) ~ Estrategia_Trabajo, data = subdatos), test = "Pillai")
})

## 
## --- Análisis para Herramienta = Google Drive ---
## 
## --- Análisis para Herramienta = Microsoft Excel ---
## 
## --- Análisis para Herramienta = Microsoft Outlook ---
## 
## --- Análisis para Herramienta = Microsoft Teams ---

## [[1]]
##                     Df  Pillai approx F num Df den Df    Pr(>F)    
## Estrategia_Trabajo   2 0.71972   43.286      6    462 < 2.2e-16 ***
## Residuals          232                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[2]]
##                     Df  Pillai approx F num Df den Df    Pr(>F)    
## Estrategia_Trabajo   2 0.78707   53.642      6    496 < 2.2e-16 ***
## Residuals          249                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[3]]
##                     Df  Pillai approx F num Df den Df    Pr(>F)    
## Estrategia_Trabajo   2 0.85167    61.31      6    496 < 2.2e-16 ***
## Residuals          249                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[4]]
##                     Df  Pillai approx F num Df den Df    Pr(>F)    
## Estrategia_Trabajo   2 0.89728   69.707      6    514 < 2.2e-16 ***
## Residuals          258                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Para cada estrategia (Presencial, Híbrido, Teletrabajo):
#   1. Filtramos el dataset por la estrategia en cuestión.
#   2. Ajustamos un MANOVA considerando únicamente Herramienta_Digital como factor.
#   3. Evaluamos la significancia mediante el test de Roy.
levels(datos$Estrategia_Trabajo) %>% lapply(function(nivel) {
  cat("\n--- Análisis para Estrategia =", nivel, "---\n")
  subdatos <- filter(datos, Estrategia_Trabajo == nivel)
  summary(manova(cbind(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas) ~ 
                   Herramienta_Digital, data = subdatos), test = "Pillai")
})

## 
## --- Análisis para Estrategia = Híbrido ---
## 
## --- Análisis para Estrategia = Presencial ---
## 
## --- Análisis para Estrategia = Teletrabajo ---

## [[1]]
##                      Df  Pillai approx F num Df den Df    Pr(>F)    
## Herramienta_Digital   3 0.72794   36.204      9   1017 < 2.2e-16 ***
## Residuals           339                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[2]]
##                      Df  Pillai approx F num Df den Df    Pr(>F)    
## Herramienta_Digital   3 0.89611   46.426      9    981 < 2.2e-16 ***
## Residuals           327                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[3]]
##                      Df  Pillai approx F num Df den Df    Pr(>F)    
## Herramienta_Digital   3 0.78466   38.017      9    966 < 2.2e-16 ***
## Residuals           322                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Visualización
promedios <- datos %>%
  dplyr::group_by(Estrategia_Trabajo, Herramienta_Digital) %>%
  dplyr::summarise(Tareas_Completadas = mean(Tareas_Completadas), Horas_Pantalla = mean(Horas_Pantalla), Tiempo_Respuesta_Tareas = mean(Tiempo_Respuesta_Tareas), .groups = "drop")

ggplot(promedios, aes(x = Herramienta_Digital, y = Tareas_Completadas, group = Estrategia_Trabajo, color = Estrategia_Trabajo)) +
  geom_line() + geom_point(size = 2) +
  labs(title = "Interacción Estrategia*Herramienta sobre Tareas_Completadas")

ggplot(promedios, aes(x = Herramienta_Digital, y = Horas_Pantalla, group = Estrategia_Trabajo, color = Estrategia_Trabajo)) +
  geom_line() + geom_point(size = 2) +
  labs(title = "Interacción Estrategia*Herramienta sobre Horas_Pantalla")

library(dplyr)
library(tibble)

# Inicializamos listas para guardar resultados
efecto_Estrategia_en_Herramienta <- list()
efecto_Herramienta_en_Estrategia <- list()

# ---------------------------------------------------------
# Efecto de la Estrategia de Trabajo dentro de cada nivel de Herramienta Digital
# ---------------------------------------------------------
for (nivel_herramienta in levels(datos$Herramienta_Digital)) {
  subdatos <- filter(datos, Herramienta_Digital == nivel_herramienta)
  manova_mod <- manova(cbind(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas) ~ 
                         Estrategia_Trabajo, data = subdatos)
  pval <- summary(manova_mod, test = "Pillai")$stats["Estrategia_Trabajo", "Pr(>F)"]
  efecto_Estrategia_en_Herramienta[[nivel_herramienta]] <- pval
}

# ---------------------------------------------------------
# Efecto de la Herramienta Digital dentro de cada nivel de Estrategia de Trabajo
# ---------------------------------------------------------
for (nivel_estrategia in levels(datos$Estrategia_Trabajo)) {
  subdatos <- filter(datos, Estrategia_Trabajo == nivel_estrategia)
  manova_mod <- manova(cbind(Tareas_Completadas, Horas_Pantalla, Tiempo_Respuesta_Tareas) ~ 
                         Herramienta_Digital, data = subdatos)
  pval <- summary(manova_mod, test = "Pillai")$stats["Herramienta_Digital", "Pr(>F)"]
  efecto_Herramienta_en_Estrategia[[nivel_estrategia]] <- pval
}

# ---------------------------------------------------------
# Convertir resultados en tablas
# ---------------------------------------------------------
tabla_efecto_Estrategia_en_Herramienta <- tibble(
  `Herramienta Digital` = names(efecto_Estrategia_en_Herramienta),
  `p-valor efecto de Estrategia` = unlist(efecto_Estrategia_en_Herramienta)
)

tabla_efecto_Herramienta_en_Estrategia <- tibble(
  `Estrategia de Trabajo` = names(efecto_Herramienta_en_Estrategia),
  `p-valor efecto de Herramienta` = unlist(efecto_Herramienta_en_Estrategia)
)

# Mostrar tablas de efectos simples
print(tabla_efecto_Estrategia_en_Herramienta)

## # A tibble: 4 × 2
##   `Herramienta Digital` `p-valor efecto de Estrategia`
##   <chr>                                          <dbl>
## 1 Google Drive                                6.32e-42
## 2 Microsoft Excel                             6.68e-51
## 3 Microsoft Outlook                           9.95e-57
## 4 Microsoft Teams                             2.40e-63

print(tabla_efecto_Herramienta_en_Estrategia)

## # A tibble: 3 × 2
##   `Estrategia de Trabajo` `p-valor efecto de Herramienta`
##   <chr>                                             <dbl>
## 1 Híbrido                                        7.81e-56
## 2 Presencial                                     8.79e-70
## 3 Teletrabajo                                    5.06e-58

library(car)
linearHypothesis(modelo, "Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Excel") #Contraste individual de interacción

## 
## Sum of squares and products for the hypothesis:
##                         Tareas_Completadas Horas_Pantalla
## Tareas_Completadas                5.985609       5.600459
## Horas_Pantalla                    5.600459       5.240091
## Tiempo_Respuesta_Tareas          17.681303      16.543580
##                         Tiempo_Respuesta_Tareas
## Tareas_Completadas                     17.68130
## Horas_Pantalla                         16.54358
## Tiempo_Respuesta_Tareas                52.23002
## 
## Sum of squares and products for error:
##                         Tareas_Completadas Horas_Pantalla
## Tareas_Completadas             677.5838227     -0.6715254
## Horas_Pantalla                  -0.6715254    639.6013548
## Tiempo_Respuesta_Tareas          1.8676226    -17.9051263
##                         Tiempo_Respuesta_Tareas
## Tareas_Completadas                     1.867623
## Horas_Pantalla                       -17.905126
## Tiempo_Respuesta_Tareas              624.719901
## 
## Multivariate Tests: 
##                  Df test stat approx F num Df den Df     Pr(>F)    
## Pillai            1 0.0926015 33.50694      3    985 < 2.22e-16 ***
## Wilks             1 0.9073985 33.50694      3    985 < 2.22e-16 ***
## Hotelling-Lawley  1 0.1020516 33.50694      3    985 < 2.22e-16 ***
## Roy               1 0.1020516 33.50694      3    985 < 2.22e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coef(modelo)

##                                                                    Tareas_Completadas
## (Intercept)                                                               10.98699671
## Estrategia_TrabajoPresencial                                              -0.79123551
## Estrategia_TrabajoTeletrabajo                                             -1.37166599
## Herramienta_DigitalMicrosoft Excel                                        -0.07250708
## Herramienta_DigitalMicrosoft Outlook                                      -1.27265077
## Herramienta_DigitalMicrosoft Teams                                        -0.39301630
## GeneroMasculino                                                           -0.04147226
## Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Excel           -0.54687475
## Estrategia_TrabajoTeletrabajo:Herramienta_DigitalMicrosoft Excel           0.89651714
## Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Outlook          0.04586897
## Estrategia_TrabajoTeletrabajo:Herramienta_DigitalMicrosoft Outlook         0.66516398
## Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Teams            0.74725895
## Estrategia_TrabajoTeletrabajo:Herramienta_DigitalMicrosoft Teams          -0.31184016
##                                                                    Horas_Pantalla
## (Intercept)                                                            8.79330324
## Estrategia_TrabajoPresencial                                          -0.63521255
## Estrategia_TrabajoTeletrabajo                                         -0.73113511
## Herramienta_DigitalMicrosoft Excel                                    -0.29009112
## Herramienta_DigitalMicrosoft Outlook                                  -0.85983564
## Herramienta_DigitalMicrosoft Teams                                    -0.58782685
## GeneroMasculino                                                       -0.02685689
## Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Excel       -0.51168550
## Estrategia_TrabajoTeletrabajo:Herramienta_DigitalMicrosoft Excel       0.30886404
## Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Outlook     -0.08386567
## Estrategia_TrabajoTeletrabajo:Herramienta_DigitalMicrosoft Outlook    -0.28753587
## Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Teams        1.29928716
## Estrategia_TrabajoTeletrabajo:Herramienta_DigitalMicrosoft Teams      -0.62874025
##                                                                    Tiempo_Respuesta_Tareas
## (Intercept)                                                                   41.749185978
## Estrategia_TrabajoPresencial                                                  -1.874431035
## Estrategia_TrabajoTeletrabajo                                                 -2.486308704
## Herramienta_DigitalMicrosoft Excel                                             0.382629408
## Herramienta_DigitalMicrosoft Outlook                                          -2.352888107
## Herramienta_DigitalMicrosoft Teams                                            -0.745433221
## GeneroMasculino                                                                0.009831585
## Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Excel               -1.615450960
## Estrategia_TrabajoTeletrabajo:Herramienta_DigitalMicrosoft Excel               0.843336101
## Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Outlook             -0.348959011
## Estrategia_TrabajoTeletrabajo:Herramienta_DigitalMicrosoft Outlook             1.761059213
## Estrategia_TrabajoPresencial:Herramienta_DigitalMicrosoft Teams                2.214663986
## Estrategia_TrabajoTeletrabajo:Herramienta_DigitalMicrosoft Teams              -0.219016740

lh.out <- linearHypothesis(modelo, hypothesis.matrix =
                             c("Estrategia_TrabajoPresencial = 0", 
                               "Herramienta_DigitalMicrosoft Excel = 0")) #Contraste conjunto (bloque) de hipótesis
lh.out

## 
## Sum of squares and products for the hypothesis:
##                         Tareas_Completadas Horas_Pantalla
## Tareas_Completadas                30.66418       20.36450
## Horas_Pantalla                    20.36450       16.18991
## Tiempo_Respuesta_Tareas           82.81166       48.62335
##                         Tiempo_Respuesta_Tareas
## Tareas_Completadas                     82.81166
## Horas_Pantalla                         48.62335
## Tiempo_Respuesta_Tareas               238.87814
## 
## Sum of squares and products for error:
##                         Tareas_Completadas Horas_Pantalla
## Tareas_Completadas             677.5838227     -0.6715254
## Horas_Pantalla                  -0.6715254    639.6013548
## Tiempo_Respuesta_Tareas          1.8676226    -17.9051263
##                         Tiempo_Respuesta_Tareas
## Tareas_Completadas                     1.867623
## Horas_Pantalla                       -17.905126
## Tiempo_Respuesta_Tareas              624.719901
## 
## Multivariate Tests: 
##                  Df test stat  approx F num Df den Df     Pr(>F)    
## Pillai            2 0.3194481  62.47469      6   1972 < 2.22e-16 ***
## Wilks             2 0.6839937  68.66523      6   1970 < 2.22e-16 ***
## Hotelling-Lawley  2 0.4569699  74.94306      6   1968 < 2.22e-16 ***
## Roy               2 0.4456794 146.47998      3    986 < 2.22e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Simulación de los datos

set.seed(123)

# Factores
Fertilizante <- factor(rep(c("F1","F2","F3"), each=12))  # 12 plantas por fertilizante
Bloque <- factor(rep(1:4, times=9))  # 4 bloques

# Covariables simuladas
HUMEDAD <- round(rnorm(36, mean=60, sd=5), 1)
PESOINI <- round(rnorm(36, mean=15, sd=2), 1)

# Efecto del fertilizante en variables respuesta
# (F3 es el mejor, F1 el peor)
media_ALTURA <- c(F1=50, F2=55, F3=62)
media_FRUTOS <- c(F1=10, F2=13, F3=16)

ALTURA <- rnorm(36, mean=media_ALTURA[Fertilizante] + 0.1*(HUMEDAD-60) + 0.5*(PESOINI-15), sd=3)
FRUTOS <- rpois(36, lambda=media_FRUTOS[Fertilizante] + 0.05*(HUMEDAD-60))

# Data frame final
datos <- data.frame(Fertilizante, Bloque, HUMEDAD, PESOINI, ALTURA, FRUTOS)
head(datos)

##   Fertilizante Bloque HUMEDAD PESOINI   ALTURA FRUTOS
## 1           F1      1    57.2    16.1 53.28722      9
## 2           F1      2    58.8    14.9 47.70240      9
## 3           F1      3    67.8    14.4 48.41597     13
## 4           F1      4    60.4    14.2 52.71671      8
## 5           F1      1    60.6    13.6 48.50568     11
## 6           F1      2    68.6    14.6 46.99785      5

# Verificación
summary(datos)

##  Fertilizante Bloque    HUMEDAD         PESOINI          ALTURA     
##  F1:12        1:9    Min.   :50.20   Min.   :10.40   Min.   :47.00  
##  F2:12        2:9    1st Qu.:57.12   1st Qu.:14.07   1st Qu.:51.52  
##  F3:12        3:9    Median :60.60   Median :14.90   Median :55.65  
##               4:9    Mean   :60.28   Mean   :15.02   Mean   :55.88  
##                      3rd Qu.:63.65   3rd Qu.:15.95   3rd Qu.:59.00  
##                      Max.   :68.90   Max.   :19.30   Max.   :68.65  
##      FRUTOS     
##  Min.   : 5.00  
##  1st Qu.: 9.75  
##  Median :13.50  
##  Mean   :13.19  
##  3rd Qu.:16.00  
##  Max.   :25.00

Hipótesis:

• H₀: Las variables dependientes siguen una distribución normal multivariada dentro de cada combinación de tratamiento (fertilizante) y bloque.

• H₁: Al menos en un grupo (combinación de fertilizante y bloque), las variables dependientes no siguen una distribución normal multivariada.

La siguiente funcion (mshapiro.test) se aplica grupo a grupo, por lo tanto primero es necesario dividir la base de datos en los grupos, en nuestro caso 3.

library(tidyverse)
library(mvnormtest)
library(dplyr)

datosc=datos
datosc

##    Fertilizante Bloque HUMEDAD PESOINI   ALTURA FRUTOS
## 1            F1      1    57.2    16.1 53.28722      9
## 2            F1      2    58.8    14.9 47.70240      9
## 3            F1      3    67.8    14.4 48.41597     13
## 4            F1      4    60.4    14.2 52.71671      8
## 5            F1      1    60.6    13.6 48.50568     11
## 6            F1      2    68.6    14.6 46.99785      5
## 7            F1      3    62.3    12.5 49.52391     10
## 8            F1      4    53.7    19.3 51.10333     14
## 9            F1      1    56.6    17.4 50.87729     11
## 10           F1      2    57.8    12.8 49.83584      9
## 11           F1      3    66.1    14.2 49.09802     10
## 12           F1      4    61.8    14.1 51.66313      8
## 13           F2      1    62.0    16.6 55.33854     10
## 14           F2      2    60.6    14.8 55.95535     16
## 15           F2      3    57.2    15.5 58.26052     14
## 16           F2      4    68.9    14.9 57.14554     14
## 17           F2      1    62.5    14.9 54.22221     15
## 18           F2      2    50.2    17.7 58.81642     19
## 19           F2      3    63.5    14.5 58.08051     15
## 20           F2      4    57.6    18.0 57.90519     10
## 21           F2      1    54.7    11.9 53.63620     13
## 22           F2      2    58.9    16.2 53.60628      9
## 23           F2      3    54.9    15.2 58.67196     16
## 24           F2      4    56.4    15.4 53.03922     17
## 25           F3      1    56.9    15.8 68.65200     22
## 26           F3      2    51.6    14.0 65.25783     16
## 27           F3      3    64.2    14.3 61.36290     14
## 28           F3      4    60.8    13.0 58.00074     20
## 29           F3      1    54.3    12.9 58.24878      9
## 30           F3      2    66.3    15.6 63.70065     18
## 31           F3      3    62.1    15.9 61.91992     15
## 32           F3      4    58.5    15.1 60.85737     17
## 33           F3      1    64.5    16.8 60.49514     11
## 34           F3      2    64.4    19.1 64.35492     25
## 35           F3      3    64.1    14.0 59.55529      6
## 36           F3      4    63.4    10.4 55.03617     17

library(tidyverse)
trat1 = datosc %>% filter(Fertilizante == "F1") %>%
  dplyr::select(ALTURA, FRUTOS)
trat2 = datosc %>% filter(Fertilizante == "F2") %>%
  dplyr::select(ALTURA, FRUTOS)
trat3 = datosc %>% filter(Fertilizante == "F3") %>%
  dplyr::select(ALTURA, FRUTOS)


# Ejecutamos el test
library(mvnormtest)
mshapiro.test(t(trat1)) 

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.90802, p-value = 0.2012

mshapiro.test(t(trat2))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.96417, p-value = 0.8413

mshapiro.test(t(trat3))

## 
##  Shapiro-Wilk normality test
## 
## data:  Z
## W = 0.92249, p-value = 0.3072

Hipótesis:

• H₀: Las matrices de varianza-covarianza de las variables dependientes son iguales en todos los tratamientos.

• H₁: Al menos una de las matrices de varianza-covarianza difiere entre los tratamientos.

library(heplots)
library(biotools)

datosc=datos
str(datosc)

## 'data.frame':    36 obs. of  6 variables:
##  $ Fertilizante: Factor w/ 3 levels "F1","F2","F3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Bloque      : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
##  $ HUMEDAD     : num  57.2 58.8 67.8 60.4 60.6 68.6 62.3 53.7 56.6 57.8 ...
##  $ PESOINI     : num  16.1 14.9 14.4 14.2 13.6 14.6 12.5 19.3 17.4 12.8 ...
##  $ ALTURA      : num  53.3 47.7 48.4 52.7 48.5 ...
##  $ FRUTOS      : int  9 9 13 8 11 5 10 14 11 9 ...

head(datosc)

##   Fertilizante Bloque HUMEDAD PESOINI   ALTURA FRUTOS
## 1           F1      1    57.2    16.1 53.28722      9
## 2           F1      2    58.8    14.9 47.70240      9
## 3           F1      3    67.8    14.4 48.41597     13
## 4           F1      4    60.4    14.2 52.71671      8
## 5           F1      1    60.6    13.6 48.50568     11
## 6           F1      2    68.6    14.6 46.99785      5

res <- boxM(datosc[, 5:6], datosc[, "Fertilizante"])
res

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  datosc[, 5:6]
## Chi-Sq (approx.) = 10.974, df = 6, p-value = 0.08917

summary(res)

## Summary for Box's M-test of Equality of Covariance Matrices
## 
## Chi-Sq:   10.97435 
## df:   6 
## p-value: 0.08917 
## 
## log of Covariance determinants:
##       F1       F2       F3 
## 3.090107 3.712221 5.731208 
## 
## Eigenvalues:
##         F1       F2        F3   pooled
## 1 5.746389 9.922473 33.015471 16.19210
## 2 3.824911 4.126456  9.339307  5.79957
## 
## Statistics based on eigenvalues:
##                  F1        F2         F3    pooled
## product   21.979427 40.944649 308.341627 93.907210
## sum        9.571300 14.048929  42.354778 21.991669
## precision  2.296389  2.914432   7.279973  4.270127
## max        5.746389  9.922473  33.015471 16.192099

heplots::boxM(cbind(ALTURA, FRUTOS) ~ Fertilizante, data = datosc)

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Y
## Chi-Sq (approx.) = 10.974, df = 6, p-value = 0.08917

boxM(datos[, 5:6], datos[, 1])

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  datos[, 5:6]
## Chi-Sq (approx.) = 10.974, df = 6, p-value = 0.08917

library(covTestR)

tomates <- unique(datosc$Fertilizante)
tomates1 <- lapply(tomates,
                    function(x){as.matrix(datosc[datosc$Fertilizante == x, 5:6])}
)

names(tomates1) <- tomates
Ahmad2017(tomates1)

## 
##  Ahmad 2017 Homogeneity of Covariance Matrices Test
## 
## data:  F1, F2 and F3
## Standard Normal = 2187.1, Mean = 0, Variance = 1, p-value < 2.2e-16
## alternative hypothesis: true difference in covariance matrices is not equal to 0

data = datosc[, c("Fertilizante", "ALTURA", "FRUTOS")]
homogeneityCovariances(data, group = Fertilizante, covTest = BoxesM)

## Warning: `invoke()` is deprecated as of rlang 0.4.0.
## Please use `exec()` or `inject()` instead.
## This warning is displayed once every 8 hours.

## 
##  Boxes' M Homogeneity of Covariance Matrices Test
## 
## data:  F1, F2 and F3
## Chi-Squared = 12.027, df = 156, p-value = 1
## alternative hypothesis: true difference in covariance matrices is not equal to 0

# Bartlett test of HOMOGENEITY of variances (parametric):

# Bartlett's K-squared =1.175  df =1  p value =0.27847


# Fligner-Killeen test of HOMOGENEITY of variances (non parametric):

# Fligner-Killeen chi-squared =1.646  df =1  p value =0.19954

library(DFA.CANCOR)
HOMOGENEITY(data = datosc[c(1,5,6)],groups = 'Fertilizante', 
            variables = c('ALTURA', 'FRUTOS'))

## 
## Covariance matrix for GroupF1

##        ALTURA FRUTOS
## ALTURA   3.91   0.40
## FRUTOS   0.40   5.66

## 
## Covariance matrix for GroupF2

##        ALTURA FRUTOS
## ALTURA   4.78   1.83
## FRUTOS   1.83   9.27

## 
## Covariance matrix for GroupF3

##        ALTURA FRUTOS
## ALTURA  13.48   8.99
## FRUTOS   8.99  28.88

## 
## 
## Bartlett test of HOMOGENEITY of variances (parametric):

## 
## Bartlett's K-squared =1.175  df =1  p value =0.27847

## 
## 
## Fligner-Killeen test of HOMOGENEITY of variances (non parametric):

## 
## Fligner-Killeen chi-squared =1.646  df =1  p value =0.19954

## 
## 
## Pooled within groups covariance matrix from SPSS:

##        ALTURA FRUTOS
## ALTURA  7.388  3.740
## FRUTOS  3.740 14.604

## 
## 
## Pooled within groups correlation matrix from SPSS:

##        ALTURA FRUTOS
## ALTURA   1.00   0.36
## FRUTOS   0.36   1.00

## 
## 
## Box Test of equality of covariance matrices:

## 
## Log determinants:

##        Log Determinant
## F1               3.090
## F2               3.712
## F3               5.731
## Pooled           4.542

## 
## 
## M = 12.027  F = 1.829  df1 = 6  df2 = 27141.23  p = 0.0893

Hipótesis:

• H₀:La matriz de correlación es una matriz identidad (las variables no están correlacionadas)

• H₁: Al menos una correlación es diferente de 0 (hay correlación entre las variables dependientes).

head(datosc)

##   Fertilizante Bloque HUMEDAD PESOINI   ALTURA FRUTOS
## 1           F1      1    57.2    16.1 53.28722      9
## 2           F1      2    58.8    14.9 47.70240      9
## 3           F1      3    67.8    14.4 48.41597     13
## 4           F1      4    60.4    14.2 52.71671      8
## 5           F1      1    60.6    13.6 48.50568     11
## 6           F1      2    68.6    14.6 46.99785      5

datos1=datosc[, 5:6]
datos1

##      ALTURA FRUTOS
## 1  53.28722      9
## 2  47.70240      9
## 3  48.41597     13
## 4  52.71671      8
## 5  48.50568     11
## 6  46.99785      5
## 7  49.52391     10
## 8  51.10333     14
## 9  50.87729     11
## 10 49.83584      9
## 11 49.09802     10
## 12 51.66313      8
## 13 55.33854     10
## 14 55.95535     16
## 15 58.26052     14
## 16 57.14554     14
## 17 54.22221     15
## 18 58.81642     19
## 19 58.08051     15
## 20 57.90519     10
## 21 53.63620     13
## 22 53.60628      9
## 23 58.67196     16
## 24 53.03922     17
## 25 68.65200     22
## 26 65.25783     16
## 27 61.36290     14
## 28 58.00074     20
## 29 58.24878      9
## 30 63.70065     18
## 31 61.91992     15
## 32 60.85737     17
## 33 60.49514     11
## 34 64.35492     25
## 35 59.55529      6
## 36 55.03617     17

library(psych)
#library(rela)
options(scipen=0)
cortest.bartlett(cor(datos1), n = nrow(datos1))

## $chisq
## [1] 17.32741
## 
## $p.value
## [1] 3.146165e-05
## 
## $df
## [1] 1

Conclusión: A un nivel de significancia del 5%, existe suficiente evidencia estadística para afirmar que si hay correlación entre las variables dependientes.

Dado esto podemos continuar con el MANCOVA en DBCA

Para el factor Fertilizante (tu factor principal):

• H₀: Una vez controlado el efecto de las covariables (HUMEDAD, PESOINI) y del bloque, los vectores de medias de ALTURA y FRUTOS son iguales en los tres fertilizantes.

• H₁: Una vez controlado el efecto de las covariables y del bloque, al menos un vector de medias es diferente entre los fertilizantes.

Para el bloque (DBCA):

• H₀: No existen diferencias en los vectores de medias de ALTURA y FRUTOS entre los bloques.

• H₁: Al menos un bloque difiere en su vector de medias.

Para las covariables:

summary(modelo, test = "Pillai")

##              Df  Pillai approx F num Df den Df    Pr(>F)    
## Fertilizante  2 0.84854  10.3170      4     56 2.489e-06 ***
## Bloque        3 0.08691   0.4240      6     56   0.85994    
## HUMEDAD       1 0.11241   1.7097      2     27   0.19993    
## PESOINI       1 0.26470   4.8598      2     27   0.01575 *  
## Residuals    28                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo, test = "Wilks")

##              Df   Wilks approx F num Df den Df    Pr(>F)    
## Fertilizante  2 0.16606  19.6287      4     54 5.057e-10 ***
## Bloque        3 0.91446   0.4115      6     54   0.86819    
## HUMEDAD       1 0.88759   1.7097      2     27   0.19993    
## PESOINI       1 0.73530   4.8598      2     27   0.01575 *  
## Residuals    28                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo, test = "Hotelling-Lawley")

##              Df Hotelling-Lawley approx F num Df den Df    Pr(>F)    
## Fertilizante  2           4.9341   32.072      4     52 1.783e-13 ***
## Bloque        3           0.0920    0.399      6     52   0.87639    
## HUMEDAD       1           0.1266    1.710      2     27   0.19993    
## PESOINI       1           0.3600    4.860      2     27   0.01575 *  
## Residuals    28                                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo, test = "Roy")

##              Df    Roy approx F num Df den Df    Pr(>F)    
## Fertilizante  2 4.9162   68.827      2     28 1.554e-11 ***
## Bloque        3 0.0710    0.663      3     28   0.58179    
## HUMEDAD       1 0.1266    1.710      2     27   0.19993    
## PESOINI       1 0.3600    4.860      2     27   0.01575 *  
## Residuals    28                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(modelo)

##  Response ALTURA :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## Fertilizante  2 792.28  396.14 68.6885 1.591e-11 ***
## Bloque        3   4.96    1.65  0.2868  0.834514    
## HUMEDAD       1  20.01   20.01  3.4699  0.073016 .  
## PESOINI       1  57.35   57.35  9.9446  0.003829 ** 
## Residuals    28 161.48    5.77                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response FRUTOS :
##              Df Sum Sq Mean Sq F value   Pr(>F)   
## Fertilizante  2 233.72 116.861  7.6501 0.002235 **
## Bloque        3  20.53   6.843  0.4479 0.720704   
## HUMEDAD       1   9.69   9.691  0.6344 0.432452   
## PESOINI       1  23.98  23.979  1.5697 0.220607   
## Residuals    28 427.72  15.276                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7.1 Comparar dos Grupos (F1 vs F3)

modelo1 = manova(cbind(ALTURA, FRUTOS) ~ 
                   Fertilizante + Bloque + HUMEDAD + PESOINI, data = datos,
                 subset = Fertilizante %in% c("F1", "F3"))

summary(modelo1, test = "Pillai")

##              Df  Pillai approx F num Df den Df    Pr(>F)    
## Fertilizante  1 0.87225   54.621      2     16 7.095e-08 ***
## Bloque        3 0.21571    0.685      6     34   0.66285    
## HUMEDAD       1 0.13755    1.276      2     16   0.30612    
## PESOINI       1 0.30127    3.449      2     16   0.05682 .  
## Residuals    17                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo1, test = "Wilks")

##              Df   Wilks approx F num Df den Df    Pr(>F)    
## Fertilizante  1 0.12775   54.621      2     16 7.095e-08 ***
## Bloque        3 0.79492    0.649      6     32   0.69087    
## HUMEDAD       1 0.86245    1.276      2     16   0.30612    
## PESOINI       1 0.69873    3.449      2     16   0.05682 .  
## Residuals    17                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo1, test = "Hotelling-Lawley")

##              Df Hotelling-Lawley approx F num Df den Df    Pr(>F)    
## Fertilizante  1           6.8276   54.621      2     16 7.095e-08 ***
## Bloque        3           0.2446    0.612      6     30   0.71917    
## HUMEDAD       1           0.1595    1.276      2     16   0.30612    
## PESOINI       1           0.4312    3.449      2     16   0.05682 .  
## Residuals    17                                                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(modelo1, test = "Roy")

##              Df    Roy approx F num Df den Df    Pr(>F)    
## Fertilizante  1 6.8276   54.621      2     16 7.095e-08 ***
## Bloque        3 0.1621    0.919      3     17   0.45283    
## HUMEDAD       1 0.1595    1.276      2     16   0.30612    
## PESOINI       1 0.4312    3.449      2     16   0.05682 .  
## Residuals    17                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En la comparación específica entre los fertilizantes F1 y F3

• Se observó que existen diferencias significativas en la combinación de ALTURA y FRUTOS, reforzando la evidencia de que el tipo de fertilizante influye en el desempeño de las plantas.

summary.aov(modelo1)

##  Response ALTURA :
##              Df Sum Sq Mean Sq  F value  Pr(>F)    
## Fertilizante  1 790.22  790.22 115.8864 5.2e-09 ***
## Bloque        3  14.94    4.98   0.7304 0.54795    
## HUMEDAD       1  17.54   17.54   2.5721 0.12718    
## PESOINI       1  42.87   42.87   6.2867 0.02261 *  
## Residuals    17 115.92    6.82                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response FRUTOS :
##              Df  Sum Sq Mean Sq F value  Pr(>F)   
## Fertilizante  1 222.042 222.042 12.5916 0.00247 **
## Bloque        3  28.458   9.486  0.5379 0.66264   
## HUMEDAD       1   0.217   0.217  0.0123 0.91292   
## PESOINI       1  51.460  51.460  2.9182 0.10578   
## Residuals    17 299.781  17.634                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7.2 Comparación General por pares

• H₀: los 2 vectores son iguales.

• H₁: los 2 vectores difieren.

Ferti <- c("F1", "F2", "F3")
comb<-t(combn(length(Ferti), 2))

for(i in 1:nrow(comb)){
  modelo.comp = manova(cbind(ALTURA, FRUTOS) ~ 
                         Fertilizante + Bloque + HUMEDAD + PESOINI,
                       data=datos,
                       subset=Fertilizante %in% Ferti[comb[i,]])
  print(paste("Trat: ",Ferti[comb[i,]][1], "y",Ferti[comb[i,]][2]))
  print(summary(modelo.comp, test = "Pillai"))
  cat("\n")
  
}

## [1] "Trat:  F1 y F2"
##              Df  Pillai approx F num Df den Df    Pr(>F)    
## Fertilizante  1 0.82867   38.693      2     16 7.425e-07 ***
## Bloque        3 0.24740    0.800      6     34   0.57679    
## HUMEDAD       1 0.32506    3.853      2     16   0.04307 *  
## PESOINI       1 0.08082    0.703      2     16   0.50956    
## Residuals    17                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [1] "Trat:  F1 y F3"
##              Df  Pillai approx F num Df den Df    Pr(>F)    
## Fertilizante  1 0.87225   54.621      2     16 7.095e-08 ***
## Bloque        3 0.21571    0.685      6     34   0.66285    
## HUMEDAD       1 0.13755    1.276      2     16   0.30612    
## PESOINI       1 0.30127    3.449      2     16   0.05682 .  
## Residuals    17                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [1] "Trat:  F2 y F3"
##              Df  Pillai approx F num Df den Df    Pr(>F)    
## Fertilizante  1 0.62366  13.2575      2     16 0.0004024 ***
## Bloque        3 0.47280   1.7543      6     34 0.1384472    
## HUMEDAD       1 0.03192   0.2638      2     16 0.7714198    
## PESOINI       1 0.33673   4.0614      2     16 0.0374573 *  
## Residuals    17                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretación General:

• El efecto significativo encontrado en la MANCOVA global (para Fertilizante) se debe principalmente a la diferencia entre (F1 y F2), (F1 y F3), y (F2 y F3).

Conclusión práctica: si tu objetivo es encontrar qué tratamientos difieren, el contraste entre (F1 y F2), (F1 y F3), y (F2 y F3) contribuyen de manera significativa al rechazo de H₀; aunque (F2 y F3) con menor magnitud que las comparaciones que involucran a F1.

Esto es coherente con tu resultado global, donde Fertilizante fue significativa en la MANCOVA: la diferencia detectada está impulsada por F1 vs F2, F1 vs F3, y F2 vs F3

CONCLUSIÓN GENERAL:

De manera práctica, esto confirma que el factor Fertilizante influye de manera clara y consistente en las variables respuesta, especialmente en ALTURA, tal como lo sugería la MANCOVA global y los ANOVA univariados. Por lo tanto, los tres fertilizantes generan respuestas multivariadas distintas en las plantas, siendo F1 el que más contribuye a las diferencias observadas.