CIENCIA DE DATOS/Estadística paramétrica y no paramétrica/ Programación R/ ANOVA versus Kruskal Wallis
Dr. Roberto Chang López/ rchang@unah.edu.hn
13/12/2020
##Haga Click Aqui para ver Certificado Machine Learning MIT https://www.credential.net/4dd365ea-ea5a-46a2-a72e-539e70545c6e
##Haga Click Aqui para ver Certificado Columbia Python for Managers https://certificates.emeritus.org/0a2e1de7-add2-4710-ad49-417d1dadfb61#gs.4a92hv ##Contacto: rchang@unah.edu.hn
Algunos Dashboards elaborados son:
Para Bolsa de Valores https://rchang.shinyapps.io/rchang-stock-exchange/
Para el Estado del Clima https://rchang.shinyapps.io/rchang-app_clima_ho/
Para Machine Learning https://rchang.shinyapps.io/rchang-app/
Para Empresariales e Industriales https://rchang.shinyapps.io/rchang-app_final_emp/
Para Dashboards con log in https://rchang.shinyapps.io/clase_3-shiny-2/_w_ae4e775f/_w_f249a9a1/?page=sign_in
y para Sistemas de Información Geográfica
Este documento presenta algunos ejemplos de ciencia de datos aplicados en ANOVA y Kruskal Wallis y en programación R.
Ejercicios Para estudiar ¿cuál de los dos tratamientos contra la artrosis es más eficaz se eligen aleatoriamente dos muestras de 10 y 22 pacientes a los cuales se les somete a los tratamientos 1 y 2, respectivamente. Pasados tres meses se valoran ambos tratamientos de manera que el que tenga mayor puntuación será el mas eficaz.
Use un test adecuado para evaluar si existen diferencias entre los dos tratamientos
H0: Los tratamientos tienen efectos iguales
H1: Existe una diferencia entre las dos tratamientos
Tratamiento1<-c(12,15,21,17,38,42,10,23,35,28)
Tratamiento2<-c(21,18,42,25,14,52,65,40,43,35,18,56,29,32,44,15,68,41,37,43,58,42)
Puntuaciones<-c(Tratamiento1, Tratamiento2)
Tratamiento<-factor(x=rep(c(1,2),c(10,22)), levels = c(1,2),labels = c("Tratamiento1","Tratamiento2"))
data<-data.frame(Tratamiento,Puntuaciones)
data## Tratamiento Puntuaciones
## 1 Tratamiento1 12
## 2 Tratamiento1 15
## 3 Tratamiento1 21
## 4 Tratamiento1 17
## 5 Tratamiento1 38
## 6 Tratamiento1 42
## 7 Tratamiento1 10
## 8 Tratamiento1 23
## 9 Tratamiento1 35
## 10 Tratamiento1 28
## 11 Tratamiento2 21
## 12 Tratamiento2 18
## 13 Tratamiento2 42
## 14 Tratamiento2 25
## 15 Tratamiento2 14
## 16 Tratamiento2 52
## 17 Tratamiento2 65
## 18 Tratamiento2 40
## 19 Tratamiento2 43
## 20 Tratamiento2 35
## 21 Tratamiento2 18
## 22 Tratamiento2 56
## 23 Tratamiento2 29
## 24 Tratamiento2 32
## 25 Tratamiento2 44
## 26 Tratamiento2 15
## 27 Tratamiento2 68
## 28 Tratamiento2 41
## 29 Tratamiento2 37
## 30 Tratamiento2 43
## 31 Tratamiento2 58
## 32 Tratamiento2 42#Tratamiento1<-c(12,15,21,17,38,42,10,23,35,28) #Tratamiento2<-c(21,18,42,25,14,52,65,40,43,35,18,56,29,32,44,15,68,41,37,43,58,42) #Puntuaciones<-c(Tratamiento1, Tratamiento2)
#Tratamiento<-factor(x=rep(c(1,2),c(10,22)), levels = c(1,2),labels = #c(“Tratamiento1”,“Tratamiento2”)) #data<-data.frame(Tratamiento,Puntuaciones)
#data ##El test de wilcoxon con su p avlor se presenta a continuación:
wilcox.test(data$Puntuaciones~ data$Tratamiento, paired=FALSE, alternative="two.sided")## Warning in wilcox.test.default(x = c(12, 15, 21, 17, 38, 42, 10, 23, 35, :
## cannot compute exact p-value with ties##
## Wilcoxon rank sum test with continuity correction
##
## data: data$Puntuaciones by data$Tratamiento
## W = 50.5, p-value = 0.01637
## alternative hypothesis: true location shift is not equal to 0###Con respecto a la media y a la mediana el tratamiento 2 es más efectivo
tapply(data$Puntuaciones, data$Tratamiento, mean)## Tratamiento1 Tratamiento2
## 24.10000 38.09091tapply(data$Puntuaciones, data$Tratamiento, median)## Tratamiento1 Tratamiento2
## 22.0 40.5Ejercicio 1: Suponga que se desea comparar el rendimiento de 5 tipos de neumáticos,
A, B, C, D y E, para lo que se decide probarlos en distintos coches
de similares características. Sus vidas medias en rodaje, medidas en miles
de kilómetros, están en el archivo neumaticos
library(readr)## Warning: package 'readr' was built under R version 4.0.5neumaticos_ <- read.csv("neumaticos_")
neumaticos_## Km Neum
## 1 68 A
## 2 72 A
## 3 77 A
## 4 42 A
## 5 53 A
## 6 72 B
## 7 53 B
## 8 63 B
## 9 53 B
## 10 48 B
## 11 60 C
## 12 82 C
## 13 64 C
## 14 75 C
## 15 72 C
## 16 48 D
## 17 61 D
## 18 57 D
## 19 64 D
## 20 50 D
## 21 64 E
## 22 65 E
## 23 70 E
## 24 68 E
## 25 53 EH0: µ1=µ2=…=µk No existe diferencia en el rendimiento de los 5 neumáticos
H1: Al menos dos medias son diferentes. Existe diferencia en el rendimiento de los 5 neumáticos
boxplot(Km~Neum,data = neumaticos_,col=c("yellow","blue","white","green", "red"))
calcular la media por rendimiento de neumático
tapply(neumaticos_$Km, neumaticos_$Neum, mean)## A B C D E
## 62.4 57.8 70.6 56.0 64.0calcular la mediana por rendimiento de neumático
tapply(neumaticos_$Km, neumaticos_$Neum, median)## A B C D E
## 68 53 72 57 65H0: µ1=µ2=…=µk No existe diferencia de insectos capturados por color de trampa
H1: Al menos dos medias son diferentes
#ANOVA
modelo4<-aov (lm(Km~Neum,data = neumaticos_) )
summary(modelo4)## Df Sum Sq Mean Sq F value Pr(>F)
## Neum 4 658.2 164.54 1.747 0.179
## Residuals 20 1883.2 94.16#Como el p-valor es muy pequeño se concluye que hay diferencia significativas #entre los 3 colores
##Los intervalos de confianza con TukeyHsD
intervalos<-TukeyHSD(modelo4)
intervalos## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = lm(Km ~ Neum, data = neumaticos_))
##
## $Neum
## diff lwr upr p adj
## B-A -4.6 -22.964503 13.764503 0.9418752
## C-A 8.2 -10.164503 26.564503 0.6728994
## D-A -6.4 -24.764503 11.964503 0.8325385
## E-A 1.6 -16.764503 19.964503 0.9988915
## C-B 12.8 -5.564503 31.164503 0.2644283
## D-B -1.8 -20.164503 16.564503 0.9982420
## E-B 6.2 -12.164503 24.564503 0.8476487
## D-C -14.6 -32.964503 3.764503 0.1620755
## E-C -6.6 -24.964503 11.764503 0.8167814
## E-D 8.0 -10.364503 26.364503 0.6921456plot(intervalos)
# validación del modelo
normalidad
hist(modelo4$residuals)
qqnorm(modelo4$residuals)
qqline(modelo4$residuals)
shapiro.test(modelo4$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo4$residuals
## W = 0.96747, p-value = 0.5817homocedasticidad
Test de bartlett
#H0: las varianzas en los grupos son iguales #H1: al menos un par de grupos tienen varianzas diferentes
bartlett.test(modelo4$residuals ~ neumaticos_$Neum)##
## Bartlett test of homogeneity of variances
##
## data: modelo4$residuals by neumaticos_$Neum
## Bartlett's K-squared = 3.1723, df = 4, p-value = 0.5294Ejercicio 2: Para evaluar el índice de alfabetización de cuatro minicipios
de una determinada región, se ha pasado un test a varios habitantes de cada
una de ellas con los siguientes resultados
Pueblo1: 78 85 90 77 69
Pueblo2: 52 48 60 35 51 47
Pueblo3: 82 91 85 74 70
Pueblo4: 57 61 45 46
Si suponemos que los datos son normales y que las varianzas son iguales
Determine si existen diferencias en el indice de alfabetización entre los pueblos
indices<-c(78,85,90,77,69,
52,48,60,35,51,47,
82,91,85,74,70,
57,61,45,46)
pueblos<-as.factor(rep(c("pueblo1","pueblo2","pueblo3","pueblo4"),c(5,6,5,4)))
datos<-data.frame(indices, pueblos)
datos## indices pueblos
## 1 78 pueblo1
## 2 85 pueblo1
## 3 90 pueblo1
## 4 77 pueblo1
## 5 69 pueblo1
## 6 52 pueblo2
## 7 48 pueblo2
## 8 60 pueblo2
## 9 35 pueblo2
## 10 51 pueblo2
## 11 47 pueblo2
## 12 82 pueblo3
## 13 91 pueblo3
## 14 85 pueblo3
## 15 74 pueblo3
## 16 70 pueblo3
## 17 57 pueblo4
## 18 61 pueblo4
## 19 45 pueblo4
## 20 46 pueblo4#indices<-c(78,85,90,77,69, # 52,48,60,35,51,47, # 82,91,85,74,70, # 57,61,45,46) #pueblos<-as.factor(rep(c(“pueblo1”,“pueblo2”,“pueblo3”,“pueblo4”),c(5,6,5,4)))
#datos<-data.frame(indices, pueblos) #datos
#ANOVA
modelo5<-aov (lm(indices~pueblos,data = datos) )
summary(modelo5)## Df Sum Sq Mean Sq F value Pr(>F)
## pueblos 3 4499 1499.7 22.43 5.63e-06 ***
## Residuals 16 1070 66.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Si existe diferencia en la alfabetización en los 4 pueblos
Prueba de Kruskal-Wallis
#Cuando no se cumplen las hipótesis exigidas por el modelo ANOVA, #es posible utilizar la prueba no paramétrica Kruskal-Wallis: #¿hay diferencias significativas entre las poblaciones?
#H0: la variable respuesta es la misma en las poblaciones valoradas #H1: la variable respuesta es mayor en alguna de las poblaciones
Ejemplo: Un estudio compara el número de huevos que pone un determinado insecto bajo
3 condiciones distintas. ¿Existen diferencias significativas dependiendo de las condiciones?
datos <- data.frame(
condicion = c(rep("condicion1", 18), rep("condicion2", 18), rep("condicion3", 18)),
n_huevos = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 16, 27, 28, 29, 30, 51, 52, 53, 342, 40,
41, 42, 43, 44, 45, 46, 47, 48, 67, 88, 89, 90, 91,92, 93, 94, 293,
19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 25, 36, 37, 58, 59, 60, 71, 72)
)tapply(datos$n_huevos,datos$condicion,median)## condicion1 condicion2 condicion3
## 12.5 57.5 26.5tapply(datos$n_huevos,datos$condicion,sd)## condicion1 condicion2 condicion3
## 78.10637 58.41750 18.59097boxplot(n_huevos~condicion,data = datos,col=c("yellow","blue","green"))
library("ggplot2")## Warning: package 'ggplot2' was built under R version 4.0.5ggplot(data = datos, mapping = aes(x = n_huevos, colour = condicion)) +
geom_histogram() +
theme_bw() +
facet_grid(. ~ condicion) +
theme(legend.position = "none")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
kruskal.test(n_huevos~condicion, data=datos)##
## Kruskal-Wallis rank sum test
##
## data: n_huevos by condicion
## Kruskal-Wallis chi-squared = 19.964, df = 2, p-value = 4.623e-05pairwise.wilcox.test(x = datos$n_huevos, g = datos$condicion, p.adjust.method = "holm" )## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: datos$n_huevos and datos$condicion
##
## condicion1 condicion2
## condicion2 0.00029 -
## condicion3 0.04795 0.00058
##
## P value adjustment method: holmejercicio
#Se dispone de un registro que contiene cientos de emails con informaciónn de cada uno de ellos. #El objetivo de estudio es intentar crear un modelo que permita filtrar ¿cuáles emails son un #spam? y ¿cuáles no, en función de determinadas características.
#En este caso se van a emplear únicamente como posibles predictores variables categóricas. #En particular, las variables que se van a estudiar como posibles predictores son:
#spam: si el email es spam (1) si no lo es (0). #to_multiple: si hay más de una persona en la lista de distribución. #format: si está en formato HTLM. #cc: si hay otras direcciones en copia. #attach: si hay archivos adjuntos. #dollar: si el email contiene la palabra dollar o el símbolo $. #inherit: si contiene la palabra inherit. #winner: si el email contiene la palabra winner. #password: si el email contiene la palabra password. #re_subj: si la palabra âRe:â está escrita en el asunto del email. #exclaim_subj: si se incluye algún signo de exclamación en el email.
email<-openintro::email
email## # A tibble: 3,921 x 21
## spam to_multiple from cc sent_email time image attach
## <fct> <fct> <fct> <int> <fct> <dttm> <dbl> <dbl>
## 1 0 0 1 0 0 2012-01-01 00:16:41 0 0
## 2 0 0 1 0 0 2012-01-01 01:03:59 0 0
## 3 0 0 1 0 0 2012-01-01 10:00:32 0 0
## 4 0 0 1 0 0 2012-01-01 03:09:49 0 0
## 5 0 0 1 0 0 2012-01-01 04:00:01 0 0
## 6 0 0 1 0 0 2012-01-01 04:04:46 0 0
## 7 0 1 1 0 1 2012-01-01 11:55:06 0 0
## 8 0 1 1 1 1 2012-01-01 12:45:21 1 1
## 9 0 0 1 0 0 2012-01-01 15:08:59 0 0
## 10 0 0 1 0 0 2012-01-01 12:12:00 0 0
## # ... with 3,911 more rows, and 13 more variables: dollar <dbl>, winner <fct>,
## # inherit <dbl>, viagra <dbl>, password <dbl>, num_char <dbl>,
## # line_breaks <int>, format <fct>, re_subj <fct>, exclaim_subj <dbl>,
## # urgent_subj <fct>, exclaim_mess <dbl>, number <fct>modelo_spam <- glm(spam ~ to_multiple + format+cc+attach + dollar+inherit+winner+ password+re_subj+exclaim_subj, data = email,
family = "binomial")
summary(modelo_spam)##
## Call:
## glm(formula = spam ~ to_multiple + format + cc + attach + dollar +
## inherit + winner + password + re_subj + exclaim_subj, family = "binomial",
## data = email)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6348 -0.4325 -0.2566 -0.0945 3.8846
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.79976 0.08935 -8.950 < 2e-16 ***
## to_multiple1 -2.84097 0.31158 -9.118 < 2e-16 ***
## format1 -1.52284 0.12270 -12.411 < 2e-16 ***
## cc 0.03134 0.01895 1.654 0.098058 .
## attach 0.20351 0.05851 3.478 0.000505 ***
## dollar -0.07304 0.02306 -3.168 0.001535 **
## inherit 0.32999 0.15223 2.168 0.030184 *
## winneryes 1.83103 0.33641 5.443 5.24e-08 ***
## password -0.75953 0.29597 -2.566 0.010280 *
## re_subj1 -3.11857 0.36522 -8.539 < 2e-16 ***
## exclaim_subj 0.24399 0.22502 1.084 0.278221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2437.2 on 3920 degrees of freedom
## Residual deviance: 1936.2 on 3910 degrees of freedom
## AIC: 1958.2
##
## Number of Fisher Scoring iterations: 7modelo_spam1 <- glm(spam ~ to_multiple + format+attach + dollar+inherit+winner+ password+re_subj+exclaim_subj, data = email,
family = "binomial")
summary(modelo_spam1)##
## Call:
## glm(formula = spam ~ to_multiple + format + attach + dollar +
## inherit + winner + password + re_subj + exclaim_subj, family = "binomial",
## data = email)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6445 -0.4337 -0.2587 -0.0938 3.8838
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.79079 0.08910 -8.875 < 2e-16 ***
## to_multiple1 -2.79514 0.30808 -9.073 < 2e-16 ***
## format1 -1.52592 0.12263 -12.443 < 2e-16 ***
## attach 0.20509 0.05798 3.537 0.000404 ***
## dollar -0.07378 0.02310 -3.194 0.001404 **
## inherit 0.33499 0.15258 2.196 0.028122 *
## winneryes 1.84356 0.33755 5.462 4.72e-08 ***
## password -0.76351 0.29641 -2.576 0.009999 **
## re_subj1 -3.10746 0.36537 -8.505 < 2e-16 ***
## exclaim_subj 0.24100 0.22486 1.072 0.283830
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2437.2 on 3920 degrees of freedom
## Residual deviance: 1938.5 on 3911 degrees of freedom
## AIC: 1958.5
##
## Number of Fisher Scoring iterations: 7modelo_spam2 <- glm(spam ~ to_multiple + format+attach + dollar+inherit+winner+ password+re_subj, data = email,
family = "binomial")
summary(modelo_spam2)##
## Call:
## glm(formula = spam ~ to_multiple + format + attach + dollar +
## inherit + winner + password + re_subj, family = "binomial",
## data = email)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6591 -0.4373 -0.2544 -0.0944 3.8707
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.78138 0.08860 -8.820 < 2e-16 ***
## to_multiple1 -2.77682 0.30752 -9.030 < 2e-16 ***
## format1 -1.51770 0.12226 -12.414 < 2e-16 ***
## attach 0.20419 0.05789 3.527 0.00042 ***
## dollar -0.06970 0.02239 -3.113 0.00185 **
## inherit 0.33614 0.15073 2.230 0.02575 *
## winneryes 1.86675 0.33652 5.547 2.9e-08 ***
## password -0.76035 0.29680 -2.562 0.01041 *
## re_subj1 -3.11329 0.36519 -8.525 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2437.2 on 3920 degrees of freedom
## Residual deviance: 1939.6 on 3912 degrees of freedom
## AIC: 1957.6
##
## Number of Fisher Scoring iterations: 7tabla1 <- table(email$spam, email$winner,
dnn = c("Spam","Ganador"))
addmargins(tabla1)## Ganador
## Spam no yes Sum
## 0 3510 44 3554
## 1 347 20 367
## Sum 3857 64 3921addmargins(prop.table(tabla1))## Ganador
## Spam no yes Sum
## 0 0.89517980 0.01122163 0.90640143
## 1 0.08849783 0.00510074 0.09359857
## Sum 0.98367763 0.01632237 1.00000000mosaicplot(tabla1)
library(ggplot2)
ggplot(email, aes(x=factor(spam),y= num_char, fill=factor(winner))) +
geom_boxplot()
Comparación de las predicciones con las observaciones
#Para este estudio se va a emplear un umbral de 0.5. Si la probabilidad predicha de que recibamos un spam #es superior a 0.5 se asigna al nivel 1 (Alerta de Spam), si es menor se asigna al #nivel 0 (No es Spam).
predicciones1 <- ifelse(modelo_spam2$fitted.values > 0.5, 1, 0)
matriz_confusion1 <- table(modelo_spam2$model$spam, predicciones1,
dnn = c("observaciones", "predicciones"))
addmargins(matriz_confusion1)## predicciones
## observaciones 0 1 Sum
## 0 3551 3 3554
## 1 355 12 367
## Sum 3906 15 3921mosaicplot(matriz_confusion1)
##¿Cómo explica el modelo en que porcentaje con la matriz de confusión?
Ccorrecta<-(3551+12)/3921*100
Ccorrecta## [1] 90.86968