11 de abril de 2017
El material de este taller ha sido creado como parte del grupo de bioestadística de estudiantes de biología de la Universidad Nacional (https://www.facebook.com/groups/109321476192304/)
Si usted quiere ser partícipe de los talleres, o quiere organizar este mismo taller en su centro de estudio o grupo de trabajo, puede contactarnos al correo:
Contacto: ronny.hernandez@gmail.com | alvarovh@gmail.com
Necesitaremos los paquetes psych, pgirmess, lsr & ggplot2. En caso de no tenerlos instalados tendremos que hacer uso de la siguiente función:
install.packages(c("psych","pgirmess","lsr","ggplot2","car"))
Recordar que hay que cargar los paquetes luego de instalados.
library("car")
library("pgirmess")
library("lsr")
library("ggplot2")
library("psych")
library("car")
library("MASS")
Es importante que en cada proyecto o trabajo que realicemos, configuremos el directorio de trabajo, de manera tal que podamos traer los datos, guardar figuras y demás objetos que se generen durante nuestro análisis de datos.
Para revisar en qué directorio de trabajo hemos iniciado sesión:
getwd()
## [1] "C:/Users/Mechas/Desktop/SC_RUG"
Para establecer el directorio de trabajo:
setwd("C:/Users/Mechas/Desktop/SC_RUG")
Vamos a utilizar el set de datos UScrime que se encuentra en el paquete MASS
str(UScrime)
## 'data.frame': 47 obs. of 16 variables: ## $ M : int 151 143 142 136 141 121 127 131 157 140 ... ## $ So : int 1 0 1 0 0 0 1 1 1 0 ... ## $ Ed : int 91 113 89 121 121 110 111 109 90 118 ... ## $ Po1 : int 58 103 45 149 109 118 82 115 65 71 ... ## $ Po2 : int 56 95 44 141 101 115 79 109 62 68 ... ## $ LF : int 510 583 533 577 591 547 519 542 553 632 ... ## $ M.F : int 950 1012 969 994 985 964 982 969 955 1029 ... ## $ Pop : int 33 13 18 157 18 25 4 50 39 7 ... ## $ NW : int 301 102 219 80 30 44 139 179 286 15 ... ## $ U1 : int 108 96 94 102 91 84 97 79 81 100 ... ## $ U2 : int 41 36 33 39 20 29 38 35 28 24 ... ## $ GDP : int 394 557 318 673 578 689 620 472 421 526 ... ## $ Ineq: int 261 194 250 167 174 126 168 206 239 174 ... ## $ Prob: num 0.0846 0.0296 0.0834 0.0158 0.0414 ... ## $ Time: num 26.2 25.3 24.3 29.9 21.3 ... ## $ y : int 791 1635 578 1969 1234 682 963 1555 856 705 ...
El dataset contiene datos de criminalidad de Estados Unidos, en los cuales existe una columna con la probabilidad de que una persona sea encarcelada de acuerdo a su origen: Southern (grupo 1) Non-Southern (grupo 0)
Para revisar los datos rápidamente podríamos generar un gráfico:
ggplot(UScrime,aes(group = So,x=So,y= Prob))+
geom_boxplot()
Recordar que uno de los supuestos a cumplir para no incurrir en errores tipo I o tipo II, es la homocedasticidad.
Con leveneTest() podremos revisar si todos los grupos contienen una cantidad equivalente de varianza, que es una precondición necesaria para el uso de la desviación estándar agrupada en el estadístico t para grupos independientes.
leveneTest(UScrime$Prob ~ UScrime$So)
La columna con la que estamos trabajando, R la ha leído como de tipo "numeric", por lo que debemos de coercionarla a factor:
UScrime$So <- as.factor(UScrime$So)
Volvemos a correr la prueba de levene:
leveneTest(UScrime$Prob ~ UScrime$So)
## Levene's Test for Homogeneity of Variance (center = median) ## Df F value Pr(>F) ## group 1 1.15 0.2893 ## 45
Repaso: ¿Qué nos dice en esta prueba el p-value?
Ahora, para correr la prueba de comparación de medias, utilizaremos la función t.test(). Para esta debemos de especificar los argumentos:
t.test(Prob ~ So,data = UScrime, var.equal = T, paired = F)
## ## Two Sample t-test ## ## data: Prob by So ## t = -4.2021, df = 45, p-value = 0.0001236 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -0.03727856 -0.01312153 ## sample estimates: ## mean in group 0 mean in group 1 ## 0.03851265 0.06371269
Podemos crear un gráfico para visualizar el comportamiento de los datos:
ggplot(UScrime,aes(x=So,y=Prob, fill=So))+
geom_boxplot()
Por defecto, la función t.test en sus argumentos supone heterocedasticidad. Por lo que si no indicamos lo contrario, la operación incluirá la variante del t por la de Welch:
t.test(Prob ~ So,data = UScrime)
## ## Welch Two Sample t-test ## ## data: Prob by So ## t = -3.8954, df = 24.925, p-value = 0.0006506 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -0.03852569 -0.01187439 ## sample estimates: ## mean in group 0 mean in group 1 ## 0.03851265 0.06371269
¿Cómo sé que no tengo independencia?
Simple, si un mismo individuo fue medido en dos ocasiones, hay falta de independencia. Esto porque la respuesta va a estar ligada a factores propios del individuo sin importar el tratamiento o experimento que se le realice. Por lo tanto es necesario especificar en la prueba que no tenemos independencia lo que se logra con una prueba pareada
Por ejemplo, usaremos el set de datos CO2 que se encuentra en el paquete datasets.
p <- CO2 head(CO2)
## Grouped Data: uptake ~ conc | Plant ## Plant Type Treatment conc uptake ## 1 Qn1 Quebec nonchilled 95 16.0 ## 2 Qn1 Quebec nonchilled 175 30.4 ## 3 Qn1 Quebec nonchilled 250 34.8 ## 4 Qn1 Quebec nonchilled 350 37.2 ## 5 Qn1 Quebec nonchilled 500 35.3 ## 6 Qn1 Quebec nonchilled 675 39.2
Este set de datos consiste en 84 filas con 5 columnas de un experimento de tolerancia a al frío en una especie de planta herbácea.
En ocasiones cuando se tienen más de 30 observaciones, se asume que se cumple con la normalidad, sin embargo, es recomendable revisar para no cometer errores:
#Se puede hacer de manera visual: qqnorm(p$uptake) qqline(p$uptake)
#Otra forma: hist(p$uptake)
O si nos sentimos inseguros con la parte visual podemos realizar una prueba de normalidad:
shapiro.test(p$uptake)
## ## Shapiro-Wilk normality test ## ## data: p$uptake ## W = 0.94105, p-value = 0.0007908
Ante la falta de normalidad, podemos preferir hacer una prueba no paramétrica como lo es wilcox.test:
wilcox.test(p$uptake~p$Treatment)
## Warning in wilcox.test.default(x = c(16, 30.4, 34.8, 37.2, 35.3, 39.2, ## 39.7, : cannot compute exact p-value with ties
## ## Wilcoxon rank sum test with continuity correction ## ## data: p$uptake by p$Treatment ## W = 1187.5, p-value = 0.006358 ## alternative hypothesis: true location shift is not equal to 0
Un ANOVA no compara medias.
Nosotros al hacer un ANOVA, estamos interesados en la variación que se debe a nuestro diseño experimental, a las categorías que creamos debido a que pensamos que hay un efecto.
Hay dos tipos de variación:
El estadístico F divide la variación entre grupos / variación dentro de los grupos.
Existen diferentes tipos de ANOVA, de acuerdo a las categorías que se creen, sin embargo en este tutorial vamos a repasar el Análisis de varianza de una vía
Utilizaremos el set de datos amis que viene en el boot. Este contiene una columna con la velocidad de automóviles muestreadas en diferentes periodos:
La pregunta que podríamos formular es si la señal tuvo algún efecto en la velocidad de los automóviles que pasaban por cierto punto y si esos periodos diferían en su efecto en la velocidad
Podemos revisar los datos un poco por periodo:
velocidad <- boot::amis library(psych)
describeBy(velocidad$speed,velocidad$period)
## $`1` ## vars n mean sd median trimmed mad min max range skew kurtosis ## X1 1 2800 37.36 6.35 37 37.15 5.93 20 65 45 0.35 0.16 ## se ## X1 0.12 ## ## $`2` ## vars n mean sd median trimmed mad min max range skew kurtosis ## X1 1 2800 37.46 6.67 37 37.08 5.93 19 67 48 0.61 0.79 ## se ## X1 0.13 ## ## $`3` ## vars n mean sd median trimmed mad min max range skew kurtosis ## X1 1 2837 38.64 6.42 38 38.44 5.93 20 61 41 0.3 -0.02 ## se ## X1 0.12 ## ## attr(,"call") ## by.default(data = x, INDICES = group, FUN = describe, type = type)
Para iniciar y ver los supuestos, podemos comprobar si tenemos homocedasticidad:
velocidad$period <- as.factor(velocidad$period) leveneTest(velocidad$speed~velocidad$period)
## Levene's Test for Homogeneity of Variance (center = median) ## Df F value Pr(>F) ## group 2 1.1395 0.32 ## 8434
Podemos proceder a realizar el ANOVA
velocidad$period <- as.factor(velocidad$period) aov_model <- aov(velocidad$speed ~ velocidad$period) summary(aov_model)
## Df Sum Sq Mean Sq F value Pr(>F) ## velocidad$period 2 2836 1418 33.75 2.53e-15 *** ## Residuals 8434 354412 42 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Hay diferencias significativas entre las categorías en este caso, los valores de las velocidades difieren entre sí de acuerdo a los periodos en los cuales se levantó la señal.
Peeeero… no nos dice entre cuáles periodos hay diferencias.
Para encontrar entre cuáles periodos hay diferencias utilizamos lo que se llama "post-hoc tests"
TukeyHSD(aov_model)
## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = velocidad$speed ~ velocidad$period) ## ## $`velocidad$period` ## diff lwr upr p adj ## 2-1 0.09642857 -0.3096896 0.5025468 0.8430912 ## 3-1 1.27265472 0.8678628 1.6774466 0.0000000 ## 3-2 1.17622614 0.7714343 1.5810180 0.0000000
Lo anterior nos demuestra que en donde hay diferencias significativas es entre el grupo 3 con 1 y 3 con 2. Entre el 2 con el 1 no hay diferencias significativas.
Un gráfico nos ayuda a visualizar dichas diferencias.
ggplot(velocidad, aes(x = period, y = speed )) +
geom_boxplot()
Hacemos una prueba no paramétrica que se conoce como Kruskal-wallis. Si no estamos seguros de que nuestros datos cumplan con alguno de los supuestos, lo mejor es utilizar una prueba más robusta.
Vamos a utilizar el dataset "yields" que consiste en un experimento realizado en un campo agrícola en donde 10 unidades experimentales con diferentes tipo de suelo fueron muestreadas en cuanto a su rendimiento de cultivos.
cultivo <- read.table("yields.txt",h=T)
head(cultivo)
## sand clay loam ## 1 6 17 13 ## 2 10 15 16 ## 3 8 3 9 ## 4 6 11 12 ## 5 14 14 15 ## 6 17 12 16
str(cultivo)
## 'data.frame': 10 obs. of 3 variables: ## $ sand: int 6 10 8 6 14 17 9 11 7 11 ## $ clay: int 17 15 3 11 14 12 12 8 10 13 ## $ loam: int 13 16 9 12 15 16 17 13 18 14
Al observar el dataset anterior nos damos cuenta que los datos no se encuentran de manera ordenada de acuerdo a tidy data:
Para no cambiar los datos "a mano", podemos hacer uso de este paquete:
library(tidyr) cultivo_a <- gather(cultivo,soil,yield) head(cultivo_a)
## soil yield ## 1 sand 6 ## 2 sand 10 ## 3 sand 8 ## 4 sand 6 ## 5 sand 14 ## 6 sand 17
str(cultivo_a)
## 'data.frame': 30 obs. of 2 variables: ## $ soil : chr "sand" "sand" "sand" "sand" ... ## $ yield: int 6 10 8 6 14 17 9 11 7 11 ...
Con este dataset tan pequeño, podemos preferir realizar una prueba no parámetrica, que en R es bastante sencillo:
library(pgirmess) kruskal.test(cultivo_a$yield~cultivo_a$soil)
Nos dice que todos los niveles de grupo deben de ser finitos, y es porque la columna que muestra la variable "soils" no se encuentra como un factor, por lo que R no reconoce cada tipo de suelo como grupos.
Una simple coerción nos ayuda:
cultivo_a$soil <- as.factor(cultivo_a$soil)
kruskal.test(cultivo_a$yield~cultivo_a$soil)
## ## Kruskal-Wallis rank sum test ## ## data: cultivo_a$yield by cultivo_a$soil ## Kruskal-Wallis chi-squared = 7.5813, df = 2, p-value = 0.02258
Sí hay diferencias significativas en el rendimiento entre los tres tipos de suelo existentes. Sin embargo, no sabemos entre cuáles hay diferencias.
La función kruskalmc del paquete pgirmess nos ayudará.
kruskalmc(cultivo_a$yield~cultivo_a$soil)
## Multiple comparison test after Kruskal-Wallis ## p.value: 0.05 ## Comparisons ## obs.dif critical.dif difference ## clay-loam 6.35 9.425108 FALSE ## clay-sand 4.40 9.425108 FALSE ## loam-sand 10.75 9.425108 TRUE
No obtenemos un p-value, sin embargo nos indica entre cuáles grupos hay diferencias, las cuales se encuentran entre loan y sand.
Podemos visualizar los datos con el error estándar. Existe una función para esto y se pueden encontrar más opciones con ejemplos en:
http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/
#Para hacer el gráfico debo de hacer un summarySE
plot <- summarySE(cultivo_a,measurevar="yield",groupvars=c("soil"))
#Establecer distancia horizontal entre puntos para evitar solapamiento:
pd <- position_dodge(0.15)
#Establecer estilo de gráfico:
theme_set(theme_bw(base_size = 26))
#Crear gráfico como objeto:
b <- ggplot(plot, aes(x=soil, y=yield, colour=soil)) +
geom_errorbar(aes(ymin=yield-se, ymax=yield+se),colour="black",width=0.1,
position=pd) + geom_line(position=pd,size=1) +
geom_point(position=pd, size=1.5, shape=21, fill="white")+
xlab("Tipo de suelo") + ylab("Rendimiento")
b
## geom_path: Each group consists of only one observation. Do you need to ## adjust the group aesthetic?
Se pueden hacer otros gráficos más aptos para diferentes audiencias.
c <- ggplot(cultivo_a, aes(x = soil, y = yield)) +
stat_summary(fun.y = mean, geom = "bar", fill = c("skyblue")) +
stat_summary(fun.data = mean_sdl, fun.args = list(mult = 1),
geom = "errorbar", width = 0.1,size=1)+
ylab("Yield")
#coord_cartesian(ylim=c(8,16)) #Con este delimito el axis
c
Cualquier información, comentario, corrección y demás pueden contactarme: