Indicación:
Comparación de todas las variables de las especies versicolor y setosa.
Decidir que parmetros usar: Mean(SD) o Median (IQR).
Incluir p value ( decidir que test usar para comparar cada variable).
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(knitr)
library(dplyr)
df<- iris
Para realizar ello, tenemos que saber si las variables tienen una distribución normal o no. Para ello realizaremos qqplots y realizaremos pruebas de normalidad
for (col in names(df)) {
if (is.numeric(df[[col]])) {
test <- shapiro.test(df[[col]])
print(sprintf("Columna %s: p-value = %.5f W=%.5f" , col, test$p.value, test$statistic))
qqnorm(df[[col]], main = sprintf('QQplot de %s', col), col = 'blue')
qqline(df[[col]], col = 'red')
}
}
## [1] "Columna Sepal.Length: p-value = 0.01018 W=0.97609"
## [1] "Columna Sepal.Width: p-value = 0.10115 W=0.98492"
## [1] "Columna Petal.Length: p-value = 0.00000 W=0.87627"
## [1] "Columna Petal.Width: p-value = 0.00000 W=0.90183"
Podemos observar que a excepción de la columna Sepal Width, las variables no siguen una distribución normal a pesar de que el QQ-plot podría sugerir que sí. La columna Sepal Width tiene un estadísitico W=0.98, muy cercano a 1 y un p-value mayor a 0.05, observando el QQ-plot podemos asegurar que sí sigue una distribución normal.
Para realizar las comparaciones entre las variables usaremos entonces la media y desviación estandar para Sepal Width y para las columnas Sepal.Length, Petal.Length y Petal.Width la mediana y el rango intercuartil.
Hagamos la separación de variables
versicolor<-df[df$Species=='versicolor',]
setosa<-df[df$Species=='setosa',]
datos_filtrados <- filter(df, Species == "versicolor" | Species == "setosa")
Ya que tenemos los datos, haremos las comparaciones entre las variables, empezando por Sepal Width.
Sepal Width
Dado que sigue una distribución normal, la prueba que se hace es la t-test la cual compara las medias de dos grupos en unidades de desviación estándar.
La interpretación del p-value es la siguiente: Si p < 0.05: Diferencias significativas (rechazas H₀). Si p ≥ 0.05: No hay evidencia suficiente para rechazar H₀.
print(sprintf("Media versicolor= %.2f. Desviación estándar= %.2f",
mean(versicolor$Sepal.Width),sd(versicolor$Sepal.Width) ))
## [1] "Media versicolor= 2.77. Desviación estándar= 0.31"
print(sprintf("Media setosa= %.2f, Desviación estándar= %.2f",
mean(setosa$Sepal.Width), sd(setosa$Sepal.Width)))
## [1] "Media setosa= 3.43, Desviación estándar= 0.38"
t.test(versicolor$Sepal.Width, setosa$Sepal.Width)
##
## Welch Two Sample t-test
##
## data: versicolor$Sepal.Width and setosa$Sepal.Width
## t = -9.455, df = 94.698, p-value = 2.484e-15
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.7961652 -0.5198348
## sample estimates:
## mean of x mean of y
## 2.770 3.428
ggplot(datos_filtrados, aes(x = Species, y = Sepal.Width, fill = Species)) +
geom_boxplot() +
labs(title = "Sepal Width",
y = "Length", x = "")
Dado que el p-value<0.05, entonces se puede asegurar que las diferencias entre las medias no es cero y que además hay una diferencia de 9.455 desviaciones estándar entre ellas y el signo negativo nos indica que la media de versicolor es mayor que la de setosa.
Ahora comparemos las medianas y IQR del resto de variables usando la prueba U de Mann-Whitney o también llamada Wilcoxon.
Sepal Length
print(sprintf("Mediana versicolor= %.2f. IQR= %.2f",
median(versicolor$Sepal.Length),IQR(versicolor$Sepal.Length) ))
## [1] "Mediana versicolor= 5.90. IQR= 0.70"
print(sprintf("Mediana setosa= %.2f, IQR= %.2f",
median(setosa$Sepal.Length), IQR(setosa$Sepal.Length)))
## [1] "Mediana setosa= 5.00, IQR= 0.40"
wilcox.test(versicolor$Sepal.Length,setosa$Sepal.Length)
##
## Wilcoxon rank sum test with continuity correction
##
## data: versicolor$Sepal.Length and setosa$Sepal.Length
## W = 2331.5, p-value = 8.346e-14
## alternative hypothesis: true location shift is not equal to 0
ggplot(datos_filtrados, aes(x = Species, y = Sepal.Length, fill = Species)) +
geom_boxplot() +
labs(title = "Sepal Length",
y = "Length", x = "")
EL valor obtenido del estadístico de la prueba de Wilcoxon (W = 2331.5) nos dice que los valores de un grupo son significativamente mayores que en el otro, del boxplot se puede apreciar que la mediana del Sepal Length de la especie versicolor es mucho mayor a la obtenida de setosa.
El p-value = 8.346e-14 mucho menor a 0.05,nos dice que hay una diferencia estadísticamente significativa entre las longitudes de sépalos de versicolor y setosa.
Petal length
print(sprintf("Mediana versicolor= %.2f. IQR= %.2f",
median(versicolor$Petal.Length),IQR(versicolor$Petal.Length) ))
## [1] "Mediana versicolor= 4.35. IQR= 0.60"
print(sprintf("Mediana setosa= %.2f, IQR= %.2f",
median(setosa$Petal.Length), IQR(setosa$Petal.Length)))
## [1] "Mediana setosa= 1.50, IQR= 0.18"
wilcox.test(versicolor$Petal.Length,setosa$Petal.Length)
##
## Wilcoxon rank sum test with continuity correction
##
## data: versicolor$Petal.Length and setosa$Petal.Length
## W = 2500, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
ggplot(datos_filtrados, aes(x = Species, y = Petal.Length, fill = Species)) +
geom_boxplot() +
labs(title = "Petal Length",
y = "Length", x = "")
El valor del estadístico fue W = 2500 en conjunto con el p-value, nos dice que los valores de un grupo son significativamente mayores que en el otro, del boxplot se puede apreciar que la mediana del Petal Length de la especie versicolor es mucho mayor a la obtenida de setosa.
Petal width
print("Petal Width")
## [1] "Petal Width"
print(sprintf("Mediana versicolor= %.2f. IQR= %.2f",
median(versicolor$Petal.Width),IQR(versicolor$Petal.Width) ))
## [1] "Mediana versicolor= 1.30. IQR= 0.30"
print(sprintf("Mediana setosa= %.2f, IQR= %.2f",
median(setosa$Petal.Width), IQR(setosa$Petal.Width)))
## [1] "Mediana setosa= 0.20, IQR= 0.10"
wilcox.test(versicolor$Petal.Width,setosa$Petal.Width)
##
## Wilcoxon rank sum test with continuity correction
##
## data: versicolor$Petal.Width and setosa$Petal.Width
## W = 2500, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
datos_filtrados <- filter(iris, Species == "versicolor" | Species == "setosa")
ggplot(datos_filtrados, aes(x = Species, y = Petal.Width, fill = Species)) +
geom_boxplot() +
labs(title = "Petal.Width",
y = "Length", x = "")
Al igual que en último caso, el valor sel estadístico (W = 2500) y el p-value confirman que son diferentes las medianas, observano el boxplot, la mediana del Petal Width de la especie versicolor es mucho mayor a la obtenida de setosa.
Finalmente, automatizamos la obtención de los parámetros y los presentamos en una tabla:
library(knitr)
library(dplyr)
# Filtrar datos (versicolor y setosa)
datos <- iris %>% filter(Species %in% c("versicolor", "setosa"))
#Vemos qué variables están distribuidas normalmente
analizar_normalidad <- function(var) {
p_valor <- shapiro.test(datos[[var]])$p.value
if (p_valor>0.05){
return("t-test")
} else{
return("Mann-Whitney")
}
}
#Aplicamos a cada variable
resultados <- lapply(names(datos)[1:4], function(col) {
metodo <- analizar_normalidad(col)
p_valor <- if (metodo == "t-test") {
t.test(datos[[col]] ~ Species, data=datos)$p.value
} else {
wilcox.test(datos[[col]] ~ Species, data=datos)$p.value
}
#Creamos la tabla para presentar los datos
data.frame(
#Columna Variable
Variable = col,
#Columna setosa
setosa = if (metodo == "t-test") {
sprintf("%.2f (%.1f)",
mean(datos[[col]][datos$Species == "setosa"]),
sd(datos[[col]][datos$Species == "setosa"]))
} else {
sprintf("%.2f [%.1f]",
median(datos[[col]][datos$Species == "setosa"]),
IQR(datos[[col]][datos$Species == "setosa"]))
},
#Columna versicolor
versicolor = if (metodo == "t-test") {
sprintf("%.2f (%.1f)",
mean(datos[[col]][datos$Species == "versicolor"]),
sd(datos[[col]][datos$Species == "versicolor"]))
} else {
sprintf("%.2f [%.1f]",
median(datos[[col]][datos$Species == "versicolor"]),
IQR(datos[[col]][datos$Species == "versicolor"]))
},
#Columna p-value
`p-value` = ifelse(p_valor < 0.001, "<0.001", sprintf("%.3f", p_valor)),
#Columna Metodo
Método = metodo
)
}) %>% bind_rows()#Juntamos todas los dataframes en uno solo
# 4. Mostrar tabla
kable(
resultados,
align = "c", # Alineación centrada
caption = "Comparación entre Iris setosa e Iris versicolor"
)
| Variable | setosa | versicolor | p.value | Método |
|---|---|---|---|---|
| Sepal.Length | 5.00 [0.4] | 5.90 [0.7] | <0.001 | Mann-Whitney |
| Sepal.Width | 3.43 (0.4) | 2.77 (0.3) | <0.001 | t-test |
| Petal.Length | 1.50 [0.2] | 4.35 [0.6] | <0.001 | Mann-Whitney |
| Petal.Width | 0.20 [0.1] | 1.30 [0.3] | <0.001 | Mann-Whitney |
Nota: Los valores representan Media (SD) para t-test o Mediana [IQR] para Mann-Whitney.