40 puntos parcial con Ancova Robusto de Wilcox
2024-09-11
\[\text{Juan Felipe Paz Paiba; Paula Gabriela Riveros Herrera}\] \[\text{Diseño en el cual la covarianza tiene efecto en la respuesta}\]
set.seed(1033686870)
# Fertilización (factor)
fertilizacion <- gl(4, 12, 48, labels = c('Control', 'Baja', 'Media', 'Alta'))
# Peso del tuberculo (respuesta, en gramos)
peso_tuberculo <- c(
rnorm(12, 78.12, 7.39), # Control
rnorm(12, 135.15, 7.74), # Baja
rnorm(12, 184.39,7.58), # Media
rnorm(12, 230.78, 7.81)) # Alta
# Salinidad
# Generar datos de salinidad (covariable)
salinidad <- sort(rnorm(48, 6.2,2.0 )) # 48 observaciones para 12 repeticiones por tratamiento
datos=data.frame(fertilizacion,peso_tuberculo,salinidad)
head(datos)## fertilizacion peso_tuberculo salinidad
## 1 Control 93.74991 0.2828302
## 2 Control 75.90368 2.1738425
## 3 Control 76.70521 2.1748540
## 4 Control 87.53347 2.3483939
## 5 Control 87.83263 3.1299089
## 6 Control 70.54761 3.2088930\[\text{Analisis Descriptivo}\]
# Calcular las medias de los tratamientos
medias <- tapply(datos$peso_tuberculo, datos$fertilizacion, mean)
# Calcular la media global de las respuestas
media_global <- mean(datos$peso_tuberculo)
# Crear el boxplot
boxplot(datos$peso_tuberculo ~ datos$fertilizacion, col = c("purple", "pink", "yellow", "blue"), notch = FALSE,
main = "Diagrama de Caja", xlab = "Tratamiento", ylab = "Peso del Tuberculo (gramos)")
# Añadir la línea horizontal en la media de todas las respuestas
abline(h = mean(datos$peso_tuberculo), col = "green")
# Añadir puntos en las medias de cada tratamiento
points(x = 1, y = medias["Control"], col = "red", pch = 19)
points(x = 2, y = medias["Baja"], col = "red", pch = 19)
points(x = 3, y = medias["Media"], col = "red", pch = 19)
points(x = 4, y = medias["Alta"], col = "red", pch = 19)
# Calcular la media del peso del Tuberculo para cada tratamiento de fertilización
medias <- tapply(datos$peso_tuberculo, datos$fertilizacion, mean)
# Calcular la varianza del peso del Tuberculo para cada tratamiento de fertilización
varianzas <- tapply(datos$peso_tuberculo, datos$fertilizacion, var)
# Calcular el coeficiente de variación del peso del Tuberculo para cada tratamiento de fertilización
cv <- tapply(datos$peso_tuberculo, datos$fertilizacion, function(x) (sd(x) / mean(x)) * 100)
# Crear un data frame con los resultados
resultados <- data.frame(
Tratamiento = names(medias),
Media = medias,
Varianza = varianzas,
Coeficiente_de_Variacion = cv
)
# Mostrar la tabla de resultados
print(resultados)## Tratamiento Media Varianza Coeficiente_de_Variacion
## Control Control 80.12683 57.31884 9.448671
## Baja Baja 137.90316 78.27230 6.415490
## Media Media 185.27351 46.30660 3.672892
## Alta Alta 231.52492 33.23541 2.490020#Grafico mentiroso
#Analisis de regresión
# Crear un gráfico de dispersión con una línea de regresión
plot(datos$salinidad, datos$peso_tuberculo, col = datos$fertilizacion, pch = 19,
main = "Grafico Mentiroso", xlab = "Salinidad del Suelo", ylab = "Peso del tuberculo (gramos)")
# Añadir la línea de regresión
abline(lm(datos$peso_tuberculo ~ datos$salinidad), col = 'orange', lwd = 2)
#Analisis de regresión Grafico real
#Crear el gráfico de dispersión
plot(datos$salinidad, datos$peso_tuberculo, pch = 19, col = datos$fertilizacion, cex = 2, asp = 1,
main = "Peso del tuberculo vs. Salinidad", xlab = "Salinidad del Suelo", ylab = "Peso del tuberculo (gramos)")
# Añadir la línea de regresión
abline(lm(datos$peso_tuberculo ~ datos$salinidad), col = 'orange', lwd = 2)
# Análisis de regresión
reg <- lm(datos$peso_tuberculo ~ datos$salinidad)
summary(reg) ##
## Call:
## lm(formula = datos$peso_tuberculo ~ datos$salinidad)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.901 -19.127 0.374 17.163 86.763
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.145 13.453 -0.085 0.933
## datos$salinidad 28.752 2.311 12.443 2.51e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27.67 on 46 degrees of freedom
## Multiple R-squared: 0.771, Adjusted R-squared: 0.766
## F-statistic: 154.8 on 1 and 46 DF, p-value: 2.512e-16# Calcular el ángulo en grados
atan(36.89) * 180 / pi## [1] 88.44723#Correlación entre las variables
cor(datos$peso_tuberculo, datos$salinidad)## [1] 0.878044\[\text{Estandarizar con el fin de realizar un gráfico estandarizado}\]
# Estandarización de las variables
peso_tuberculo_est <- scale(datos$peso_tuberculo)
salinidad_est <- scale(datos$salinidad)
# Crear un data frame con las variables estandarizadas
datos_est <- data.frame(
fertilizacion = datos$fertilizacion,
peso_tuberculo_est = peso_tuberculo_est,
salinidad_est = salinidad_est)
# Crear el gráfico de dispersión
plot(datos_est$peso_tuberculo_est, datos_est$salinidad_est, pch = 19, col = as.factor(datos_est$fertilizacion), cex = 2,
main = "Peso del Tuberculo Estandarizado vs. Salinidad Estandarizada", xlab = "Peso del Tuberculo (Estandarizado)", ylab = "Salinidad del Suelo (Estandarizada)")
# Añadir la línea de regresión
abline(lm(datos_est$peso_tuberculo_est ~ datos_est$salinidad_est), col = 'orange', lwd = 2)
\[Hipotesis\] \[H_{0_1}= \mu_{Tc} = \mu_{T1} =
\mu_{T2}=\mu_{T3}\] \[H_{0_2}:
\delta=0\] \[H_{0_c3}:\gamma=0\] \[\text{ANCOVA No Parametrico robusto de Wilcox}
\]
# Ajustar el modelo de ANCOVA no paramétrico
# Primero, ajustamos la respuesta por la covariable
residuos <- residuals(lm(peso_tuberculo ~ salinidad, data = datos))
# Agregar los residuos al data frame
datos$residuos <- residuos
# Realizar comparaciones por pares utilizando la prueba de Wilcoxon
niveles <- levels(datos$fertilizacion)
combinaciones <- combn(niveles, 2, simplify = FALSE)
resultados <- list()
for (i in seq_along(combinaciones)) {
grupo1 <- combinaciones[[i]][1]
grupo2 <- combinaciones[[i]][2]
resultado <- wilcox.test(residuos ~ fertilizacion, data = subset(datos, fertilizacion %in% c(grupo1, grupo2)))
resultados[[paste(grupo1, grupo2, sep = "_vs_")]] <- resultado
}
# Mostrar los resultados
for (nombre in names(resultados)) {
print(paste("Comparación:", nombre))
print(resultados[[nombre]])
print("")
}## [1] "Comparación: Control_vs_Baja"
##
## Wilcoxon rank sum exact test
##
## data: residuos by fertilizacion
## W = 60, p-value = 0.5137
## alternative hypothesis: true location shift is not equal to 0
##
## [1] ""
## [1] "Comparación: Control_vs_Media"
##
## Wilcoxon rank sum exact test
##
## data: residuos by fertilizacion
## W = 46, p-value = 0.1432
## alternative hypothesis: true location shift is not equal to 0
##
## [1] ""
## [1] "Comparación: Control_vs_Alta"
##
## Wilcoxon rank sum exact test
##
## data: residuos by fertilizacion
## W = 31, p-value = 0.01727
## alternative hypothesis: true location shift is not equal to 0
##
## [1] ""
## [1] "Comparación: Baja_vs_Media"
##
## Wilcoxon rank sum exact test
##
## data: residuos by fertilizacion
## W = 18, p-value = 0.001115
## alternative hypothesis: true location shift is not equal to 0
##
## [1] ""
## [1] "Comparación: Baja_vs_Alta"
##
## Wilcoxon rank sum exact test
##
## data: residuos by fertilizacion
## W = 19, p-value = 0.001433
## alternative hypothesis: true location shift is not equal to 0
##
## [1] ""
## [1] "Comparación: Media_vs_Alta"
##
## Wilcoxon rank sum exact test
##
## data: residuos by fertilizacion
## W = 29, p-value = 0.01209
## alternative hypothesis: true location shift is not equal to 0
##
## [1] ""\[\text{Se encontraron diferencias significativas en las medias de los siguientes tratamientos:} \] \[\text{Control-Alta; Baja-Media; Baja-Alta y Media-Alta } \]
#Analisis de Covarianza
mod_cov <- aov(datos$peso_tuberculo ~ datos$salinidad + datos$fertilizacion)
summary(mod_cov)## Df Sum Sq Mean Sq F value Pr(>F)
## datos$salinidad 1 118541 118541 2278.3 <2e-16 ***
## datos$fertilizacion 3 32979 10993 211.3 <2e-16 ***
## Residuals 43 2237 52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\[\text{Los datos proveen evidencia en contra de la hipotesis nula, por lo tanto, al menos un tratamiento se comportó diferente} \] \[\text{Los datos proveen evidencia en contra de la hipotesis nula, por lo tanto, la salinidad, tubo algún efecto en la respuesta} \] \[\text{Los datos proveen evidencia en contra de la hipotesis nula, por lo tanto, la pendiente no fue 0} \]
#Supuestos para el análisis paramétrico
#Extraer los residuales del modelo ANOVA
residuales <- residuals(mod_cov)
# Calcular la varianza de los residuales
varianza_residuales <- var(residuales)
# Crear el histograma de los residuales
hist(residuales, main = "Histograma de Residuales", xlab = "Residuales", ylab = "Frecuencia", col = "purple2", border = "black")
# Añadir una línea vertical en la varianza de los residuales
abline(v = varianza_residuales, col = "black", lwd = 2)
# Añadir una leyenda con la varianza de los residuales
legend("topleft", legend = paste("Varianza = ", round(varianza_residuales, 4)), col = "black", lwd = 2)
#Prueba de normalidad de los residuales
# Realizar la prueba de Shapiro-Wilk en los residuales
shapiro_test_result = shapiro.test(residuales)
# Mostrar los resultados de la prueba
print(shapiro_test_result)##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.95287, p-value = 0.05208#Prueba igualdad de varianzas
# Realizar la prueba de Bartlett
bartlett_test_result <- bartlett.test(datos$peso_tuberculo ~ datos$fertilizacion)
# Mostrar los resultados de la prueba
print(bartlett_test_result)##
## Bartlett test of homogeneity of variances
##
## data: datos$peso_tuberculo by datos$fertilizacion
## Bartlett's K-squared = 2.0363, df = 3, p-value = 0.5649# Crear el diagrama de caja de los residuales (sospechas de datos atipicos)
boxplot(residuales, main = "Diagrama de Caja de Residuales", ylab = "Residuales", col = "purple", border = "black", outline = TRUE)
# Realizar la prueba de Tukey HSD
tukey_result <- TukeyHSD(mod_cov, "datos$fertilizacion")## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## datos$salinidad# Mostrar los resultados de la prueba de Tukey HSD
print(tukey_result)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = datos$peso_tuberculo ~ datos$salinidad + datos$fertilizacion)
##
## $`datos$fertilizacion`
## diff lwr upr p adj
## Baja-Control -1.875415 -9.745043 5.994214 0.9195058
## Media-Control 16.870951 9.001323 24.740580 0.0000053
## Alta-Control 28.983890 21.114262 36.853519 0.0000000
## Media-Baja 18.746366 10.876737 26.615994 0.0000006
## Alta-Baja 30.859305 22.989676 38.728933 0.0000000
## Alta-Media 12.112939 4.243310 19.982567 0.0009602\[\text{Hay diferencias entre los tratamientos: Media-Control, Alta-Control, Media-baja, Alta-Baja y Alta-Media} \] \[\text{Diseño donde la covariable no tiene efecto en la respuesta} \]
set.seed(1033686870)
# Generar datos aleatoreos de salinidad suelo (covariable) utilizando runif
salinidad_aleatorizada <- runif(48, 4.2, 8.2) # 48 observaciones para 12 repeticiones por tratamiento
# Crear el data frame con los datos
datos_2 <- data.frame(
fertilizacion = fertilizacion,
peso_tuberculo = peso_tuberculo,
salinidad_aleatorizada = salinidad_aleatorizada
)
# Redondear los datos a dos decimales
datos_2$peso_tuberculo <- round(datos_2$peso_tuberculo, 2)
datos_2$salinidad_aleatorizada <- round(datos_2$salinidad_aleatorizada, 2)
datos_2## fertilizacion peso_tuberculo salinidad_aleatorizada
## 1 Control 93.75 8.13
## 2 Control 75.90 6.12
## 3 Control 76.71 5.73
## 4 Control 87.53 5.88
## 5 Control 87.83 5.90
## 6 Control 70.55 4.83
## 7 Control 71.08 7.79
## 8 Control 76.65 8.09
## 9 Control 80.09 7.82
## 10 Control 71.76 5.79
## 11 Control 84.08 4.81
## 12 Control 85.59 5.09
## 13 Baja 137.76 4.88
## 14 Baja 133.11 7.38
## 15 Baja 147.76 5.89
## 16 Baja 150.42 4.47
## 17 Baja 134.93 6.62
## 18 Baja 131.07 5.88
## 19 Baja 151.76 4.98
## 20 Baja 123.90 6.35
## 21 Baja 135.09 7.36
## 22 Baja 144.43 4.95
## 23 Baja 128.60 7.58
## 24 Baja 136.01 7.77
## 25 Media 179.61 6.73
## 26 Media 189.60 5.39
## 27 Media 184.70 5.78
## 28 Media 180.35 7.83
## 29 Media 187.05 7.99
## 30 Media 181.39 6.51
## 31 Media 194.19 8.10
## 32 Media 179.51 7.03
## 33 Media 177.15 6.15
## 34 Media 179.20 6.13
## 35 Media 195.72 5.40
## 36 Media 194.81 4.23
## 37 Alta 228.32 8.14
## 38 Alta 235.21 7.73
## 39 Alta 236.01 4.49
## 40 Alta 225.00 5.07
## 41 Alta 231.66 6.19
## 42 Alta 245.56 4.70
## 43 Alta 230.74 7.74
## 44 Alta 229.99 6.89
## 45 Alta 229.14 5.00
## 46 Alta 234.76 6.60
## 47 Alta 226.61 6.38
## 48 Alta 225.30 6.06#Analisis de regresión
# Crear un gráfico de dispersión con una línea de regresión
plot(datos_2$salinidad_aleatorizada, datos$peso_tuberculo, col = datos_2$fertilizacion, pch = 19,
main = "Peso del Tuberculo vs. Salinidad", xlab = "Salinidad del Suelo", ylab = "Peso del tuberculo (gramos)")
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.3# Análisis de regresión
ggplot(data = datos_2, aes(x = salinidad_aleatorizada, y = peso_tuberculo, color = fertilizacion)) +
geom_point() +
geom_smooth(method = 'lm') +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
# Análisis de regresión
reg <- lm(datos_2$peso_tuberculo ~ datos_2$salinidad_aleatorizada)
summary(reg)##
## Call:
## lm(formula = datos_2$peso_tuberculo ~ datos_2$salinidad_aleatorizada)
##
## Residuals:
## Min 1Q Median 3Q Max
## -90.686 -41.491 4.484 44.982 84.100
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 169.553 45.617 3.717 0.000546 ***
## datos_2$salinidad_aleatorizada -1.722 7.120 -0.242 0.809990
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 57.78 on 46 degrees of freedom
## Multiple R-squared: 0.00127, Adjusted R-squared: -0.02044
## F-statistic: 0.05848 on 1 and 46 DF, p-value: 0.81# Calcular el ángulo en grados
atan(36.89) * 180 / pi## [1] 88.44723# Correlación entre las variables
cor(datos_2$peso_tuberculo, datos_2$salinidad_aleatorizada)## [1] -0.0356327\[Hipotesis\] \[H_{0_1}= \mu_{Tc} = \mu_{T1} = \mu_{T2}=\mu_{T3}\] \[H_{0_2}: \delta=0\] \[H_{0_c3}:\gamma=0\]
#Analisis de Covarianza
mod_cov_2 <- aov(datos_2$peso_tuberculo ~ datos_2$salinidad_aleatorizada + datos$fertilizacion)
summary(mod_cov_2)## Df Sum Sq Mean Sq F value Pr(>F)
## datos_2$salinidad_aleatorizada 1 195 195 3.876 0.0555 .
## datos$fertilizacion 3 151396 50465 1001.857 <2e-16 ***
## Residuals 43 2166 50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\[\text{Los datos proveen evidencia a favor de la hipotesis nula, por lo tanto, la salinidad no tubo algún efecto en la respuesta}\] \[\text{Los datos proveen evidencia a favor de que la pendiente de la salinidad fué nula, por lo cual, esta no afectó la respuesta}\] \[\text{Analisis Parametrico, sin tomar en cuenta la covariable puesto que esta no tuvo ningun efecto en la respuesta}\]
#Analisis de Varianza
mod_cov_3 <- aov(datos$peso_tuberculo ~ datos$fertilizacion)
summary(mod_cov_3)## Df Sum Sq Mean Sq F value Pr(>F)
## datos$fertilizacion 3 151390 50463 938.3 <2e-16 ***
## Residuals 44 2366 54
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\[\text{Los datos proveen evidencia a en contra de la hipotesis nula, por lo cual, algun tratamiento se comportó diferente}\]
#Supuestos para el análisis paramétrico
# Extraer los residuales del modelo ANOVA
residuales_2 <- residuals(mod_cov_3)
# Calcular la varianza de los residuales
varianza_residuales_2 <- var(residuales)
# Crear el histograma de los residuales
hist(residuales_2, main = "Histograma de Residuales", xlab = "Residuales", ylab = "Frecuencia", col = "purple2", border = "black")
# Añadir una línea vertical en la varianza de los residuales
abline(v = varianza_residuales_2, col = "black", lwd = 2)
# Añadir una leyenda con la varianza de los residuales
legend("topleft", legend = paste("Varianza = ", round(varianza_residuales_2, 4)), col = "black", lwd = 2)
#Prueba de normalidad de los residuales
# Realizar la prueba de Shapiro-Wilk en los residuales
shapiro_test_result = shapiro.test(residuales)
# Mostrar los resultados de la prueba
print(shapiro_test_result)##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.95287, p-value = 0.05208#Prueba igualdad de varianzas
# Realizar la prueba de Bartlett
bartlett_test_result <- bartlett.test(datos$peso_tuberculo ~ datos$fertilizacion)
# Mostrar los resultados de la prueba
print(bartlett_test_result)##
## Bartlett test of homogeneity of variances
##
## data: datos$peso_tuberculo by datos$fertilizacion
## Bartlett's K-squared = 2.0363, df = 3, p-value = 0.5649# Crear el diagrama de caja de los residuales (sospechas de datos atipicos)
boxplot(residuales, main = "Diagrama de Caja de Residuales", ylab = "Residuales", col = "purple", border = "black", outline = TRUE)
# Realizar la prueba de Tukey HSD
tukey_result <- TukeyHSD(mod_cov, "datos$fertilizacion")## Warning in replications(paste("~", xx), data = mf): non-factors ignored:
## datos$salinidad# Mostrar los resultados de la prueba de Tukey HSD
print(tukey_result)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = datos$peso_tuberculo ~ datos$salinidad + datos$fertilizacion)
##
## $`datos$fertilizacion`
## diff lwr upr p adj
## Baja-Control -1.875415 -9.745043 5.994214 0.9195058
## Media-Control 16.870951 9.001323 24.740580 0.0000053
## Alta-Control 28.983890 21.114262 36.853519 0.0000000
## Media-Baja 18.746366 10.876737 26.615994 0.0000006
## Alta-Baja 30.859305 22.989676 38.728933 0.0000000
## Alta-Media 12.112939 4.243310 19.982567 0.0009602\[\text{Hay diferencias entre los tratamientos: Media-Control, Alta-Control, Media-baja, Alta-Baja y Alta-Media}\] \[\text{Analisis no Parametrico}\]
#ANALISIS NO PARAMETRICO
# Realizar la prueba de Kruskal-Wallis
kruskal_result <- kruskal.test(datos$peso_tuberculo ~ datos$fertilizacion)
# Mostrar los resultados de la prueba
print(kruskal_result)##
## Kruskal-Wallis rank sum test
##
## data: datos$peso_tuberculo by datos$fertilizacion
## Kruskal-Wallis chi-squared = 44.082, df = 3, p-value = 1.45e-09#Prueba de comparaciónes multiples de Bonferroni
# Cargar el paquete dunn.test
library(dunn.test)## Warning: package 'dunn.test' was built under R version 4.2.3# Realizar la prueba de Dunn (COMPARACIONES MULTIPLES NO PARAMETRICO)
dunn_result <- dunn.test(datos$peso_tuberculo, datos$fertilizacion, method = "bonferroni")## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 44.0816, df = 3, p-value = 0
##
##
## Comparison of x by group
## (Bonferroni)
## Col Mean-|
## Row Mean | Alta Baja Control
## ---------+---------------------------------
## Baja | 4.199125
## | 0.0001*
## |
## Control | 6.298687 2.099562
## | 0.0000* 0.1073
## |
## Media | 2.099562 -2.099562 -4.199125
## | 0.1073 0.1073 0.0001*
##
## alpha = 0.05
## Reject Ho if p <= alpha/2\[\text{Hay diferencias entre los tratamientos: Alta-Baja, Alta-Control y Media-Control }\]