Comparaciones entre grupos

ronny_hdez & Álvaro Vega

15 de enero de 2017

Taller de R Comparaciones entre grupos I

El material de este taller ha sido creado por R-luminatis para estudiantes de biología de la Universidad Nacional de Costa Rica (e interesados en aprender).

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 a los correos:

Contacto: ronny.hernandez@gmail.com / alvarovh95@gmail.com

Temas de la sesión:

Introducción
- ¿Para qué hacer análisis estadísticos?
- Occam´s razor
- ¿Qué quieren decir mis p-values?
- ¿Cuándo hacemos pruebas clásicas?
- Errores tipo I y tipo II
- ¿Qué prueba realizo?
Dos grupos:
- t.test (independiente/dependiente)
- wilcox.test

¿Para qué hacer análisis estadísticos?

En ciencia, la objetividad es importante.
Casi nunca tenemos la posibilidad como investigadores de estudiar toda la población, por ende usamos muestras.
La estadística es una serie de fundamentos matemáticos y probabilísticos que nos ayudan a hacer inferencias a partir de las muestras de una población.

Occam´s razor:

Lo simple es lo mejor

¿Qué quieren decir mis p-values?

Ayudan a rechazar H0

Si el p-value es bajo (más bajo que el nivel de significancia alfa, que usualmente es 5%) podemos decir que es muy improbable datos que se comporten de la manera que predice la H0, por ende, rechazamos H0.

Si el p-value es alto (más alto que el nivel de significancia alfa, que usualmente es 5%) podemos decir que es muy probable observar datos que se comporten como predice H0, y por ende no rechazamos H0.

¡p-value está basado en el supuesto de que H0 es verdadera!

Errores tipo I y tipo II

Tipo I: No se acepta H0 cuando en la población es verdadera

Tipo II: Se acepta H0 cuando en la población es falsa.

¿Qué prueba realizo?

Paquetes a utilizar en la sesión:

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","car","pgirmess","lsr","ggplot2"))

Recordar que hay que cargar los paquetes luego de instalados.

library("car")
library("pgirmess")
library("lsr")
library("ggplot2")
library("psych")

Inicio de sesión:

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/Documents/Rmarkdown I Taller"

Para establecer el directorio de trabajo:

setwd("~/Rmarkdown I Taller")

Traer los datos que tenemos en el directorio de trabajo. Los datos con los que trabajaremos poseen dos variables continuas que corresponden a índices de complejidad estructural como MIG y Aniso.

dt <- read.csv("fotos.csv", h=T, stringsAsFactors = F)

Estructura datos

Vamos a revisar la estructura del set de datos:

head(dt)

##   i   plot  MIG.grey Aniso.grey        Timestamp         DateAndTime sitio
## 1 1 LADERA 0.2220831  0.9363847 21/06/2016 11:47 2015:05:21 13:21:50  Piro
## 2 2 LADERA 0.2334211  0.9371384 21/06/2016 11:47 2015:05:21 14:21:50  Piro
## 3 3 LADERA 0.1685896  0.9400231 21/06/2016 11:48 2015:05:22 09:00:01  Piro
## 4 4 LADERA 0.1554022  0.9261470 21/06/2016 11:48 2015:05:22 10:00:01  Piro
## 5 5 LADERA 0.2006901  0.9033331 21/06/2016 11:48 2015:05:22 11:00:01  Piro
## 6 6 LADERA 0.2428779  0.9122833 21/06/2016 11:49 2015:05:22 12:00:01  Piro
##   camara
## 1      1
## 2      1
## 3      1
## 4      1
## 5      1
## 6      1

Estructura datos

str(dt)

## 'data.frame':    3652 obs. of  8 variables:
##  $ i          : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ plot       : chr  "LADERA" "LADERA" "LADERA" "LADERA" ...
##  $ MIG.grey   : num  0.222 0.233 0.169 0.155 0.201 ...
##  $ Aniso.grey : num  0.936 0.937 0.94 0.926 0.903 ...
##  $ Timestamp  : chr  "21/06/2016 11:47" "21/06/2016 11:47" "21/06/2016 11:48" "21/06/2016 11:48" ...
##  $ DateAndTime: chr  "2015:05:21 13:21:50" "2015:05:21 14:21:50" "2015:05:22 09:00:01" "2015:05:22 10:00:01" ...
##  $ sitio      : chr  "Piro" "Piro" "Piro" "Piro" ...
##  $ camara     : int  1 1 1 1 1 1 1 1 1 1 ...

Esta función nos muestra de manera resumida cómo están estructurados los datos. En primer lugar nos muestra el tipo de datos y la cantidad de variables y observaciones. Luego, nos indica por columna, la clase de cada una y una muestra de los primeros valores que contiene.

t.test

Vamos a crear subsets por cada uno de los plots que existen (CIMA & LADERA):

cima <- subset(dt,plot=="CIMA")
ladera <- subset(dt,plot=="LADERA")

Obtener estadísticas de resumen por cada uno de los subconjuntos generados:

describe(cima)

##              vars    n   mean     sd median trimmed    mad  min     max
## i               1 1435 504.47 355.10 457.00  488.25 450.71 1.00 1174.00
## plot*           2 1435    NaN     NA     NA     NaN     NA  Inf    -Inf
## MIG.grey        3 1435   0.21   0.04   0.20    0.21   0.03 0.09    0.56
## Aniso.grey      4 1435   0.97   0.04   0.97    0.97   0.03 0.85    1.06
## Timestamp*      5 1435    NaN     NA     NA     NaN     NA  Inf    -Inf
## DateAndTime*    6 1435    NaN     NA     NA     NaN     NA  Inf    -Inf
## sitio*          7 1435    NaN     NA     NA     NaN     NA  Inf    -Inf
## camara          8 1435   7.18   0.39   7.00    7.10   0.00 7.00    8.00
##                range  skew kurtosis   se
## i            1173.00  0.30    -1.25 9.37
## plot*           -Inf    NA       NA   NA
## MIG.grey        0.47  1.69    12.53 0.00
## Aniso.grey      0.21 -0.92     0.18 0.00
## Timestamp*      -Inf    NA       NA   NA
## DateAndTime*    -Inf    NA       NA   NA
## sitio*          -Inf    NA       NA   NA
## camara          1.00  1.65     0.72 0.01

describe(ladera)

##              vars    n   mean     sd median trimmed    mad  min     max
## i               1 2217 556.51 323.10 555.00  554.75 410.68 1.00 1171.00
## plot*           2 2217    NaN     NA     NA     NaN     NA  Inf    -Inf
## MIG.grey        3 2217   0.18   0.03   0.18    0.18   0.03 0.07    0.30
## Aniso.grey      4 2217   0.96   0.02   0.96    0.96   0.02 0.89    1.02
## Timestamp*      5 2217    NaN     NA     NA     NaN     NA  Inf    -Inf
## DateAndTime*    6 2217    NaN     NA     NA     NaN     NA  Inf    -Inf
## sitio*          7 2217    NaN     NA     NA     NaN     NA  Inf    -Inf
## camara          8 2217   2.06   1.00   3.00    2.07   0.00 1.00    3.00
##                range  skew kurtosis   se
## i            1170.00  0.03    -1.16 6.86
## plot*           -Inf    NA       NA   NA
## MIG.grey        0.23  0.59     0.92 0.00
## Aniso.grey      0.13 -0.36    -0.22 0.00
## Timestamp*      -Inf    NA       NA   NA
## DateAndTime*    -Inf    NA       NA   NA
## sitio*          -Inf    NA       NA   NA
## camara          2.00 -0.11    -1.99 0.02

Para explorar datos vamos a generar gráficos y así visualizar alguna diferencia.

ggplot(dt,aes(x=plot,y=MIG.grey,fill=sitio))+geom_boxplot()

Lo que observamos es que hay varios outliers. Estos los eliminamos ya que son observaciones debidas a un mal funcionamiento de los aparatos de medición.

dt <- subset(dt,MIG.grey<=0.5)

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(dt$MIG.grey~dt$plot)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    1   25.37 4.961e-07 ***
##       3647                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ahora, para correr la prueba de comparación de medias, utilizaremos la función t.test(). Para esta debemos de especificar los argumentos:

x: Columna de los datos que contiene el índice del primer grupo.

y: Columna de los datos que contiene el índice del segundo grupo.

var.equal: indicamos si la varianza es igual o no en los grupos a comparar

paired: este argumento sirve para indicar si hay independencia o no.

Para medir el tamaño de efecto lo hacemos mediante Cohen´s D a través de la función cohensD.

#Comparación índice MIG:
t.test(cima$MIG.grey,ladera$MIG.grey,var.equal=F,paired = F)

## 
##  Welch Two Sample t-test
## 
## data:  cima$MIG.grey and ladera$MIG.grey
## t = 25.861, df = 2567.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.02778092 0.03233948
## sample estimates:
## mean of x mean of y 
## 0.2090410 0.1789808

#Comparación índice Anisotropía:
t.test(cima$Aniso.grey,ladera$Aniso.grey,var.equal = F,paired = F)

## 
##  Welch Two Sample t-test
## 
## data:  cima$Aniso.grey and ladera$Aniso.grey
## t = 5.6208, df = 1937.2, p-value = 2.176e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.004202733 0.008707204
## sample estimates:
## mean of x mean of y 
## 0.9650715 0.9586166

#Calcular Cohen´s D:
cohensD(cima$Aniso.grey,ladera$Aniso.grey,method="pooled")

## [1] 0.2156327

De acuerdo a la prueba que realizamos, vamos a crear un gráfico que nos ayude a visualizar los datos:

ggplot(dt,aes(x=plot,y=MIG.grey,color=plot))+geom_boxplot()

¿Qué determinamos con lo anterior?

No hay homocedasticidad
Hay independencia (por diseño experimental)
Supuesto de normalidad por CLT
t.test indica que hay diferencias significativas en cada uno de los dos índices (MIG y ANISO) entre los plots ubicados en CIMA y LADERA.

Práctica:

Podemos realizar otra prueba t.test, pero en lugar de comparar CIMA y LADERA, podríamos comparar los sitios Piro y Rincon:

-Realice la prueba t.test entre los sitios Piro y Rincon y presente los resultados con un gráfico de cajas.

Solución práctica:

#Subconjuntos
piro <- subset(dt,sitio=="Piro")
rincon <- subset(dt,sitio=="Rincon")

leveneTest(dt$MIG.grey~dt$sitio)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    1  13.445 0.0002491 ***
##       3647                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

t.test(rincon$MIG.grey,piro$MIG.grey,var.equal = F,paired = F)

## 
##  Welch Two Sample t-test
## 
## data:  rincon$MIG.grey and piro$MIG.grey
## t = 1.4108, df = 3556.4, p-value = 0.1584
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0006271554  0.0038456240
## sample estimates:
## mean of x mean of y 
## 0.1913150 0.1897057

cohensD(rincon$MIG.grey,piro$MIG.grey,method="pooled")

## [1] 0.04670805

Datos no independientes

¿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

El set de datos contiene mediciones de estrés de hojas de una planta obtenidas a través de su fluorescencia. Este parámetro se midió en las mismas hojas antes y después de someterlas al tratamiento de estrés.

fv <- read.csv("fv.csv",h=T, stringsAsFactors = F)
head(fv)

##     X IND POS ID TRT Medicion      fv_a
## 1 151   A   D 11 34C  Inicial 0.7543683
## 2 152   A   D 12 34C  Inicial 0.7220108
## 3 153   A   D 13 34C  Inicial 0.8236920
## 4 154   A   D 14 34C  Inicial 0.7837575
## 5 155   A   D 15 34C  Inicial 0.7754600
## 6 156   A   D 16 34C  Inicial 0.7753925

Vamos a comprobar si hay diferencias significativas entre la primera medición de fluorescencia y luego del tratamiento por estrés térmico a 34 C.

tc_1 <- subset(fv,Medicion=="Inicial")
tc_2 <- subset(fv,Medicion=="F")

#Homocedasticidad:
leveneTest(fv$fv_a~fv$Medicion)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  1.1677 0.2821
##       118

#normalidad
shapiro.test(tc_1$fv_a)

## 
##  Shapiro-Wilk normality test
## 
## data:  tc_1$fv_a
## W = 0.98438, p-value = 0.6377

shapiro.test(tc_2$fv_a)

## 
##  Shapiro-Wilk normality test
## 
## data:  tc_2$fv_a
## W = 0.9391, p-value = 0.004935

#Prueba pareada (paired t.test):
t.test(tc_1$fv_a,tc_2$fv_a,var.equal = T,paired = T)

## 
##  Paired t-test
## 
## data:  tc_1$fv_a and tc_2$fv_a
## t = 82.448, df = 59, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.2485877 0.2609542
## sample estimates:
## mean of the differences 
##               0.2547709

ggplot(fv,aes(x=Medicion,y=fv_a,color=Medicion))+geom_boxplot()

Datos sin normalidad:

Si no tenemos normalidad, debemos de aplicar una prueba más robusta que lidie con esto. Para esto podemos utilizar la prueba de wilcox

Para este ejemplo utilizaremos el set de datos FV. El set de datos contiene mediciones de estrés de hojas de una planta obtenidas a través de su fluorescencia. Este parámetro se midió en las mismas hojas antes y después de someterlas al tratamiento de estrés.

fv_2 <- read.csv("FVI.csv",h=T,stringsAsFactors = F)
head(fv_2)

##   IND POS ID TRT Medicion       FVI
## 1   A   D  1 27C        F 0.6913333
## 2   A   D  2 27C        F 0.6773333
## 3   A   D  3 27C        F 0.6450000
## 4   A   D  4 27C        F 0.6703333
## 5   A   D  5 27C        F 0.6673333
## 6   A   D  6 27C        F 0.6910000

Creamos un subset con las variables a comparar. En este caso compararemos dos grupos de mediciones finales (que son independientes entre sí ya que cada medición corresponde a una muestra distinta). Lo haremos con el TRT 34C y el TRT 45C y comprobaremos si hay normalidad:

shapiro.test(fv_2$FVI)

## 
##  Shapiro-Wilk normality test
## 
## data:  fv_2$FVI
## W = 0.89699, p-value < 2.2e-16

qqnorm(fv_2$FVI);qqline(fv_2$FVI)

Como observamos en la prueba y en el gráfico anterior, no tenemos normalidad, y por lo tanto es preferible hacer una prueba robusta que nos ayude a no incurrir en errores tipo I o tipo II.

tc <- subset(fv_2,TRT=="34C")
cc <- subset(fv_2,TRT== "45C")
wilcox.test(tc$FVI,cc$FVI)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  tc$FVI and cc$FVI
## W = 13379, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

El gráfico para visualizar las diferencias:

a <- rbind(tc,cc)
ggplot(a, aes(y=FVI,x=TRT,fill=TRT))+geom_boxplot()