MANOVA IRIS
2023-09-02
#1.- PREVIEW DATA SETS
library(MASS)## Warning: package 'MASS' was built under R version 4.2.3data("iris")library(DT)## Warning: package 'DT' was built under R version 4.2.3DT::datatable(iris, rownames=TRUE, options = list(autowidth=TRUE,sClass="alignRight", className = 'dt-center', pageLength=10, digit=4))2.BOXPLOTS
str(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.3library(gridExtra)## Warning: package 'gridExtra' was built under R version 4.2.3p1 <- ggplot(iris, aes(x=Species, y=Sepal.Length), aes(fill = factor(Species))) +
ggtitle("Sepal.Length") + geom_boxplot(aes(fill = factor(Species))) +
guides(fill=guide_legend(title="Species")) + geom_jitter()
p2 <- ggplot(iris, aes(x=Species, y=Sepal.Width), aes(fill = factor(Species))) +
ggtitle("Sepal.Width") + geom_boxplot(aes(fill = factor(Species))) +
guides(fill=guide_legend(title="Species")) + geom_jitter()
p3 <- ggplot(iris, aes(x=Species, y=Petal.Length), aes(fill = factor(Species))) +
ggtitle("Petal.Length") + geom_boxplot(aes(fill = factor(Species))) +
guides(fill=guide_legend(title="Species")) + geom_jitter()
p4 <- ggplot(iris, aes(x=Species, y=Petal.Width), aes(fill = factor(Species))) +
ggtitle("Petal.Width") + geom_boxplot(aes(fill = factor(Species))) +
guides(fill=guide_legend(title="Species")) + geom_jitter()
grid.arrange(p1, p2, p3, p4, ncol=2, nrow=2)
grid.arrange(p1, p2, ncol=2, nrow=1)
En Sepal Length, las medias son diferentes, y en Sepal Width la media de Setosa es distinta a la de las otras variables.
grid.arrange(p3, p4, ncol=2, nrow=1)
En Petal.Length y en Petal.Width, las medias son diferentes.
str(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...3.- Histogramas
ggplot(data=iris, mapping = aes(x=Sepal.Length, colour = Species))+
geom_histogram()+
facet_grid(.~Species) +
theme(legend.position = "bottom")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data=iris, mapping = aes(x=Sepal.Width, colour = Species))+
geom_histogram()+
facet_grid(.~Species) +
theme(legend.position = "bottom")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data=iris, mapping = aes(x=Petal.Length, colour = Species))+
geom_histogram()+
facet_grid(.~Species) +
theme(legend.position = "bottom")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data=iris, mapping = aes(x=Petal.Width, colour = Species))+
geom_histogram()+
facet_grid(.~Species) +
theme(legend.position = "bottom")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
4.- Análisis descriptivo (EDA)
iris.Setosa <- iris[iris$Species == "setosa",]
str(iris.Setosa)## 'data.frame': 50 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...summary(iris.Setosa)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.300 Min. :1.000 Min. :0.100
## 1st Qu.:4.800 1st Qu.:3.200 1st Qu.:1.400 1st Qu.:0.200
## Median :5.000 Median :3.400 Median :1.500 Median :0.200
## Mean :5.006 Mean :3.428 Mean :1.462 Mean :0.246
## 3rd Qu.:5.200 3rd Qu.:3.675 3rd Qu.:1.575 3rd Qu.:0.300
## Max. :5.800 Max. :4.400 Max. :1.900 Max. :0.600
## Species
## setosa :50
## versicolor: 0
## virginica : 0
##
##
## iris.Versicolor <- iris[iris$Species == "versicolor",]
str(iris.Versicolor)## 'data.frame': 50 obs. of 5 variables:
## $ Sepal.Length: num 7 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 ...
## $ Sepal.Width : num 3.2 3.2 3.1 2.3 2.8 2.8 3.3 2.4 2.9 2.7 ...
## $ Petal.Length: num 4.7 4.5 4.9 4 4.6 4.5 4.7 3.3 4.6 3.9 ...
## $ Petal.Width : num 1.4 1.5 1.5 1.3 1.5 1.3 1.6 1 1.3 1.4 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 2 2 2 2 2 2 2 2 2 2 ...summary(iris.Versicolor)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## Min. :4.900 Min. :2.000 Min. :3.00 Min. :1.000 setosa : 0
## 1st Qu.:5.600 1st Qu.:2.525 1st Qu.:4.00 1st Qu.:1.200 versicolor:50
## Median :5.900 Median :2.800 Median :4.35 Median :1.300 virginica : 0
## Mean :5.936 Mean :2.770 Mean :4.26 Mean :1.326
## 3rd Qu.:6.300 3rd Qu.:3.000 3rd Qu.:4.60 3rd Qu.:1.500
## Max. :7.000 Max. :3.400 Max. :5.10 Max. :1.800iris.Virginica <- iris[iris$Species == "virginica",]
str(iris.Virginica)## 'data.frame': 50 obs. of 5 variables:
## $ Sepal.Length: num 6.3 5.8 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 ...
## $ Sepal.Width : num 3.3 2.7 3 2.9 3 3 2.5 2.9 2.5 3.6 ...
## $ Petal.Length: num 6 5.1 5.9 5.6 5.8 6.6 4.5 6.3 5.8 6.1 ...
## $ Petal.Width : num 2.5 1.9 2.1 1.8 2.2 2.1 1.7 1.8 1.8 2.5 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 3 3 3 3 3 3 3 3 3 3 ...summary(iris.Virginica)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.900 Min. :2.200 Min. :4.500 Min. :1.400
## 1st Qu.:6.225 1st Qu.:2.800 1st Qu.:5.100 1st Qu.:1.800
## Median :6.500 Median :3.000 Median :5.550 Median :2.000
## Mean :6.588 Mean :2.974 Mean :5.552 Mean :2.026
## 3rd Qu.:6.900 3rd Qu.:3.175 3rd Qu.:5.875 3rd Qu.:2.300
## Max. :7.900 Max. :3.800 Max. :6.900 Max. :2.500
## Species
## setosa : 0
## versicolor: 0
## virginica :50
##
##
## summary(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## # Create Dataframe of Summary Statistics for each Sepal/Pedal metric by Specie
summary.dataframe <- data.frame(c(summary(iris.Setosa$Sepal.Length)), c(summary(iris.Versicolor$Sepal.Length)), c(summary(iris.Virginica$Sepal.Length)), c(summary(iris.Setosa$Sepal.Width)), c(summary(iris.Versicolor$Sepal.Width)), c(summary(iris.Virginica$Sepal.Width)), c(summary(iris.Setosa$Petal.Length)), c(summary(iris.Versicolor$Petal.Length)), c(summary(iris.Virginica$Petal.Length)), c(summary(iris.Setosa$Petal.Width)), c(summary(iris.Versicolor$Petal.Width)), c(summary(iris.Virginica$Petal.Width)))
# Transpose dataframe
summary.dataframe1 <- t(summary.dataframe)
measures<- c(rep("Sepal Length",3) , rep("Sepal Width",3), rep("Petal Length",3), rep("Petal Width",3))
species.type<- c(rep(c("Setosa", "Versicolor", "Virginica"),4))
summary.dataframe2<- data.frame(measures, species.type, summary.dataframe1)
colnames(summary.dataframe2) <- c("Pedal/Sepal", "Species", "Minimum", "1st Quartile", "Median", "Mean", "3rd Quartile", "Maximum")
# Output Table
DT::datatable(summary.dataframe2, rownames=FALSE, options = list(autowidth=TRUE,sClass="alignRight", className = 'dt-center', dom='tips', pageLength=12, digit=3))5.- MANOVA
str(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...modelo1a <- manova(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, iris)
summary(modelo1a)## Df Pillai approx F num Df den Df Pr(>F)
## Species 2 1.1919 53.466 8 290 < 2.2e-16 ***
## Residuals 147
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1El estadístico 53.466 y el p valor es 0.000 < 0.05. Tenemos evidencia empírica suficiente para rechazaer H0. En alguna/s variables, las medias no son iguales para cada planta.
modelo1a <- manova(cbind(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) ~ Species, iris)
summary.aov(modelo1a)## Response Sepal.Length :
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 63.212 31.606 119.26 < 2.2e-16 ***
## Residuals 147 38.956 0.265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Sepal.Width :
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.345 5.6725 49.16 < 2.2e-16 ***
## Residuals 147 16.962 0.1154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Petal.Length :
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 437.10 218.551 1180.2 < 2.2e-16 ***
## Residuals 147 27.22 0.185
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Petal.Width :
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 80.413 40.207 960.01 < 2.2e-16 ***
## Residuals 147 6.157 0.042
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sepal.Length
Tenemos evidencia empírica suficiente para rechazar H0 al 5% e indicar que la media de la variable es diferente en cada una de las plantas.
Sepal.Width
Tenemos evidencia empírica suficiente para rechazar H0 al 5% e indicar que la media de la variable es diferente en cada una de las plantas.
Petal.Length
Tenemos evidencia empírica suficiente para rechazar H0 al 5% e indicar que la media de la variable es diferente en cada una de las plantas.
Petal.Width
Tenemos evidencia empírica suficiente para rechazar H0 al 5% e indicar que la media de la variable es diferente en cada una de las plantas.
m1 <- aov(Sepal.Length~Species, data = iris)
summary(m1)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 63.21 31.606 119.3 <2e-16 ***
## Residuals 147 38.96 0.265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m1 <- aov(Sepal.Length~Species, data = iris)
plot(m1)
Se aceptan las hipótesis de independencia, media cero, ausencia de
extremos, normalidad y homocedasticidad. Los resultados son
concluyentes.
m2 <- aov(Sepal.Width~Species, data = iris)
plot(m2)
Se aceptan las hipótesis de independencia, media cero, ausencia de
extremos, normalidad y homocedasticidad. Los resultados son
concluyentes.
m3 <- aov(Petal.Length~Species, data = iris)
plot(m3)
Se aceptan las hipótesis de independencia, media cero, ausencia de
extremos, normalidad y homocedasticidad. Los resultados son
concluyentes.
m4 <- aov(Petal.Width~Species, data = iris)
plot(m4)
Se aceptan las hipótesis de independencia, media cero, ausencia de
extremos, normalidad y homocedasticidad. Los resultados son
concluyentes.
#M1
shapiro.test(m1$residuals)##
## Shapiro-Wilk normality test
##
## data: m1$residuals
## W = 0.9879, p-value = 0.2189Se acepta la normalidad.
#M2
shapiro.test(m2$residuals)##
## Shapiro-Wilk normality test
##
## data: m2$residuals
## W = 0.98948, p-value = 0.323Se acepta la normalidad.
#M3
shapiro.test(m3$residuals)##
## Shapiro-Wilk normality test
##
## data: m3$residuals
## W = 0.98108, p-value = 0.03676Al 1% sí se acepta la hipótesis de normalidad.
#M4
shapiro.test(m4$residuals)##
## Shapiro-Wilk normality test
##
## data: m4$residuals
## W = 0.97217, p-value = 0.003866No se acepta la normalidad.
#M1
t1<-TukeyHSD(m1)
t1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Length ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa 0.930 0.6862273 1.1737727 0
## virginica-setosa 1.582 1.3382273 1.8257727 0
## virginica-versicolor 0.652 0.4082273 0.8957727 0plot(t1)
Virginica-Setosa muestra una diferencia de medias positiva y mayor al
resto; por lo tanto, su media es especialmente mayor al resto.
$Species diff lwr upr p adj versicolor-setosa 0.930 0.6862273 1.1737727 0 virginica-setosa 1.582 1.3382273 1.8257727 0 virginica-versicolor 0.652 0.4082273 0.8957727 0
Versicolor y Setosa tienen un IC al 95% que oscila de 0.68 a 1.17 (+;+); nos indica que la media de Sepal.Length es superior en Versicolor que en Setosa.
Virginica y Setosa tienen un IC al 95% que oscila de 1.33 a 1.82 (+;+); nos indica que la media de Sepal.Length es superior en Virginica que en Setosa.
Virginica y Versicolor tienen un IC al 95% que oscila de 0.40 a 0.89 (+;+); nos indica que la media de Sepal.Length es superior en Virginica que en versicolor.
#M2
t2<-TukeyHSD(m2)
t2## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802plot(t2)
Virginica y Versicolor muestran una media muy superior al resto de
comparaciones para la variable Sepal.Width.
$Species diff lwr upr p adj versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000 virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000 virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802
Versicolor y Setosa tienen un IC al 95% que oscila de -0.81 a -0.49 (-;-); nos indica que la media de Sepal.Width es inferior en Versicolor que en Setosa.
Virginica y Setosa tienen un IC al 95% que oscila de -0.61 a -0.29 (-;-); nos indica que la media de Sepal.Width es inferior en Virginica que en Setosa.
Virginica y Versicolor tienen un IC al 95% que oscila de 0.43 a 0.36 (+;+); nos indica que la media de Sepal.Width es superior en Virginica que en versicolor.
#M3
t3<-TukeyHSD(m3)
t3## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Petal.Length ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa 2.798 2.59422 3.00178 0
## virginica-setosa 4.090 3.88622 4.29378 0
## virginica-versicolor 1.292 1.08822 1.49578 0plot(t3)
$Species diff lwr upr p adj versicolor-setosa 2.798 2.59422 3.00178 0 virginica-setosa 4.090 3.88622 4.29378 0 virginica-versicolor 1.292 1.08822 1.49578 0
Versicolor y Setosa tienen un IC al 95% que oscila de 2.59 a 3-00 (+;+); nos indica que la media de Petal.Length es superior en Versicolor que en Setosa.
Virginica y Setosa tienen un IC al 95% que oscila de 3.88 a 4.29 (+;+); nos indica que la media de Petal.Length es superior en Virginica que en Setosa.
Virginica y Versicolor tienen un IC al 95% que oscila de 1.08 a 1.49 (+;+); nos indica que la media de Petal.Length es superior en Virginica que en versicolor.
# m4 <- aov(Petal.Width~Species, data = iris) esto cuando sea normal
m4b <- kruskal.test(Petal.Width~Species, data=iris)
m4b##
## Kruskal-Wallis rank sum test
##
## data: Petal.Width by Species
## Kruskal-Wallis chi-squared = 131.19, df = 2, p-value < 2.2e-16Tenemos evidencia empírica para rechazar la hipótesis de igualdad de medias para cada especie en la variable Petal.Width.
pairwise.wilcox.test(x=iris$Petal.Width, g=iris$Species, p.adjust.method = "holm")##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: iris$Petal.Width and iris$Species
##
## setosa versicolor
## versicolor <2e-16 -
## virginica <2e-16 <2e-16
##
## P value adjustment method: holmVersicolor con Setosa tiene p valor 0.000 < 0.05. Las medias son diferentes en el Petal.Width.
Virginica con Setosa tiene p valor 0.000 < 0.05. Las medias son diferentes en el Petal.Width.
Virginica con Setosa tiene p valor 0.000 < 0.05. Las medias son diferentes en el Petal.Width.
6.- Informe de la Sesión
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
)## # A tibble: 70 × 3
## package loadedversion source
## <chr> <chr> <chr>
## 1 bslib 0.4.2 CRAN (R 4.2.2)
## 2 cachem 1.0.6 CRAN (R 4.2.2)
## 3 callr 3.7.3 CRAN (R 4.2.3)
## 4 cli 3.6.0 CRAN (R 4.2.2)
## 5 colorspace 2.1-0 CRAN (R 4.2.3)
## 6 crayon 1.5.2 CRAN (R 4.2.3)
## 7 crosstalk 1.2.0 CRAN (R 4.2.3)
## 8 devtools 2.4.5 CRAN (R 4.2.3)
## 9 digest 0.6.31 CRAN (R 4.2.2)
## 10 dplyr 1.1.2 CRAN (R 4.2.3)
## # ℹ 60 more rows