1: T de student para comparasión de dos medias

1.1 Introducción:

La teoria de distribución que llevo a la solución del problema de estimar la media de la población \(\mu\), cuando la varianza \(\sigma^{2}\) es desconocida fue abordado por William S. Gossett bajo el seudónimo de “Student”. El determino que cuando el tamaño de la muestra es grande, su distribución se acercaba a la normalidad. De esta manera, la prueba para comparar las medias de dos muestras pequeñas se le conoce como prueba “t de Student”.

1.2 Supuestos de la prueba

Para esta prueba es necesario considerar el nivel de medición de las variables, para lo cual se requiere de una variable cuantitativa y una variable nominal (dos grupos).

Las condiciones para calcular intervalos de confianza o aplicar un test basados en la distribución T-student son:

  • Independencia: Las observaciones tienen que ser independientes unas de las otras. Para ello el muestreo tiene que ser aleatoria y el tamaño de la muestra inferior al 10% de la población.
  • Normalidad: Las poblaciones que se comparan tienn que distribuirse de forma normal. En caso de cierta asimetría, los T-test son considerablemente robustos cuando el tamaño de muestra es mayor o igual a 30.
  • Homogeneidad de varianzas: La varianza de ambas poblaciones comparadas debe ser igual. En caso de no cumplirse con esta condición se puede usar una prueba Welch Two Sample t-test. Esta corrección se incorpora a traves de los grados de libertad permitiendo compensar la diferencia de varianzas (con el inconveniente de que se pierde precisión).

NOTA: La prueba de Welch se conisdera en general la mas segura. Usualmente, ambos procedimientos proporcionan resultados muy similares a menos que el tamaño de ambos grupos y las desviaciones estandar sean muy diferentes. Para saber mas sobre la prueba Welch entra aqui

1.3: Establecer las hipótesis

La hopótesis nula (\(H_0\)) considera que no hay diferencia o cambio. Suele contener en su definición el simbolo \(=\) de manera que al comparar dos medias independientes, suponemos que \(\mu_1 = \mu_2\). La hipotesis alternativa (\(H_A\)) considera que el valor real de las medias de ambas poblaciones son diferentes, de manera que \(\mu_1 \ne \mu_2\)

Box: Como interpretar la prueba

El concepto básico en esta prueba es el error estandar de la media (SEM). Este describe la variación en el promedio de n valores aleatorios with una media \(\mu\) y una varianza \(\sigma^{2}\) de manera:

\(SEM = \sigma / \sqrt{n}\)

Este valor implica que si se repitiera el experimento varias veces y se calculara el promedio en cada experimento, entonces estos promedios seguirian una distribucición mas estrecha que la distribución original.

De esta manera, podemos calcular un valor de distribución t para una muestra:

\(t =\frac{\hat{x} - \mu_0} {SEM}\)

O para comparar la media de dos muestras

\(t = \frac{\hat{x}_2 - \hat{x}_1}{SEDM}\)

Entonces, los pasos para determinar una prueba, son:

  • PASO 1: Planteamiento de la hipótesis nula y alternativa

[original: https://towardsdatascience.com/statistical-tests-when-to-use-which-704557554740]

  • PASO 2: Elegir el nivel de significancia (\(\alpha=0.05\) o en su forma 5%)

  • PASO 3: Determinación de la zona de aceptación y rechazo de la hipóthesis nula (\(H_0\)). Este paso se puede obtener utilizando tablas de distribución o un calculador, como este calculador. Este valor dependerá de los grados de libertad y el nivel de significancia \(\alpha\)

  • PASO 4: Determinación del valor t. Este es el que se obtiene usando la funcion en R

  • PASO 5: Ubicar el valor obtenido en el calculo dentro de la región de rechazo o de aceptación.

  • Finalmente, puedes obtener manualmente el valor P en esta zona usando esta herramienta

2: Comparacion de medias poblacionales independientes (no pareadas)

2.1: Datos

Para estos ejercicios, utilizaremos la base de datos de pingünos del Archipielago de Palmer.Puedes ver la descripción de la base aqui

El archivo penguins_size.csv contiene las variables:

  • species: Especies de pingüino (Chnstrap, Adelie o Gentoo)
  • culmen_length_mm: longitud del culmen (mm)
  • culmen_depth_mm: profundidad del culmen
  • flipper_length_mm: longitud de la aleta (mm)
  • body_mass_g: masa corporal (g)
  • island: nombre de la isla (Dream, Torgersen o Biscoe)
  • sex: sexo

Los pinguinos:

Descripción de las variables:

imagenes originales de:https://github.com/allisonhorst/palmerpenguins

Para esta sesión, necesitaremos tres paquetes:

  • Tidyverse: Para usar dplyr y ggplot2
  • car: contiene algunas funciones estadisticas
  • rstatixs: Contiene funciones estadisticas compatibles con el ambiente de tidyverse (pipe-friendly). Puedes conocer mas sobre este paquete aqui
library(tidyverse)
library(car)
library(rstatix)

Abrimos la tabla penguins_size.csv

pin <- read_csv("data/penguins_size.csv")
glimpse(pin)
## Rows: 344
## Columns: 7
## $ species           <chr> "Adelie", "Adelie", "Adelie", "Adelie", "Adelie",...
## $ island            <chr> "Torgersen", "Torgersen", "Torgersen", "Torgersen...
## $ culmen_length_mm  <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34....
## $ culmen_depth_mm   <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18....
## $ flipper_length_mm <dbl> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, ...
## $ body_mass_g       <dbl> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 347...
## $ sex               <chr> "MALE", "FEMALE", "FEMALE", NA, "FEMALE", "MALE",...

2.2: ¿Hay diferencias en la longitud del culmen entre entre las poblaciones de Dream y Biscoe de la especie Adelie?

pin_Ade <- pin %>%
  filter(species == "Adelie", island %in% c("Dream", "Biscoe"))

El primer paso es realizar una exploración de los datos y verificar los supuestos de la prueba; la prueba t requiere que los datos tengan una distribución normal y con varianzas iguales.

Grafico de histogramas para visualizar la distribución de los datos entre ambas islas

ggplot(pin_Ade, aes(x = culmen_length_mm, fill = island, col = island))+
  geom_histogram(alpha = 0.3)+
  labs(y = "frecuencia",
       x = "longitud del culmen (mm)",
       title="Grafico de densidad", 
       fill = "isla", 
       col = "isla")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Tambien es posible realizar un gráfico qq (quantil-quantil) para evaluar visualmente si la distribución de los datos de ambas islas se asemejan a una distribución normal

ggplot(pin_Ade, aes(sample = culmen_length_mm, col = island))+
  stat_qq()+
  stat_qq_line()+
  facet_grid(.~ island)

¿Que es un QQPLOT?

Método gráfico para el diagnostico de diferencias entre la distribución de probabilidad de una muestra de una población mediante el trazado de sus cuantiles uno contra el otro. De manera general, un grafico QQ permite obvservar cuan cerca está la distribución de un conjunto de datos a alguna distribución idal o comparar la distribución de dos conjuntos de datos. La linea que por defecto cruza los puntos del primer (0.25) y el tercer (0.75) cuartil es una aproximación robusta de los valores esperados de nuestros datos si siguieran una distribución normal. Si los datos se alejan de esta linea, especialmente cerca del centro, nos sugeriría que nuestros datos no se distribuyen normalmente. Para saber mas de este gráfico y como generarlo con comandos base, visita aqui

De acuerdo a los gráficos, es evidente que ambas poblaciones tienen distribuciones muy parecida a la normal. Para validar estos resutlados, realizaremos una prueba de Shapiro_wilks para evaluar la significancia de la normalidad.

El método de Shapiro-Wilks es ampliamente recomendado y tiene mayor poder en comparación de K-S. Para conocer mas sobre pruebas de normalidad visita aqui

pin_Ade_shap <- pin_Ade %>%
  group_by(island) %>%
  shapiro_test(culmen_length_mm) %>%
  print()
## # A tibble: 2 x 4
##   island variable         statistic     p
##   <chr>  <chr>                <dbl> <dbl>
## 1 Biscoe culmen_length_mm     0.977 0.520
## 2 Dream  culmen_length_mm     0.984 0.640
La hipótesis nula de esta prueba es que la muestra sigue una distribución normal. Si la prueba es significativa, la distribución es no-normal.

Tanto la gráfica de densidad como el qqplot dan soporte a que los datos tienen una distribucion normal. Para poder dalre validez, hace la prueba de shapiro. La prueba de shapiro nos indica que los datos tienen una distribución normal

El siguiente paso es evaluar el supuesto de homogeneidad de varianzas. Las gráficas de densidades nuevamente nos dan un indicativo visual de la disperción de los datos, pero es posible evaluar la significancia con una prueba de Levene implementada en el paquete car:

leveneTest(culmen_length_mm ~ island, data = pin_Ade)
## 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  0.0054 0.9415
##       98

La prueba nos dice que las varianzas son homogeneas entre los grupos

El siguiente paso para la exploración de los datos es obtener el promedio de los valores de la longitud del culmen:

pin_Ade_sum <- pin_Ade %>%
  group_by(island) %>%
  summarise(promedio = mean(culmen_length_mm),
            desviacion = sd(culmen_length_mm),
            N = n())
pin_Ade_sum
## # A tibble: 2 x 4
##   island promedio desviacion     N
##   <chr>     <dbl>      <dbl> <int>
## 1 Biscoe     39.0       2.48    44
## 2 Dream      38.5       2.47    56

De manera alternativa, este procedimiento se puede realizar con ayuda del paquete rstatix

pin_Ade_sum_r <- pin_Ade %>%
  group_by(island) %>%
  get_summary_stats(culmen_length_mm)
pin_Ade_sum_r
## # A tibble: 2 x 14
##   island variable     n   min   max median    q1    q3   iqr   mad  mean    sd
##   <chr>  <chr>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Biscoe culmen_~    44  34.5  45.6   38.7  37.7  40.7  3.02  2.74  39.0  2.48
## 2 Dream  culmen_~    56  32.1  44.1   38.6  36.8  40.4  3.62  2.74  38.5  2.46
## # ... with 2 more variables: se <dbl>, ci <dbl>

Posteriormente, hacemos un boxplot para visualizar los datos:

ggplot(pin_Ade, aes(x = island, y = culmen_length_mm, fill = island))+
  geom_boxplot(outlier.shape = "")+
  geom_point(position = position_jitter(0.1), col = "grey50")+
  theme_light()+
  labs(x = "Isla",
       y = "longitud del culmen (mm)",
       title = "Longitud del culmen por isla en el pingüino Adelie",
       fill = "Isla", col = "Isla")

Finalmente, hacemos la prueba t de student. Esta puede realizarse siguiendo el comando base en R t.test()

t.test(culmen_length_mm ~ island, data = pin_Ade, var.equal = TRUE )
## 
##  Two Sample t-test
## 
## data:  culmen_length_mm by island
## t = 0.95016, df = 98, p-value = 0.3444
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.5151265  1.4615551
## sample estimates:
## mean in group Biscoe  mean in group Dream 
##             38.97500             38.50179

En esta prueba, la hipótesis nula es que la media de las dos poblaciones es igual, lo que impricla que la hipótesis alternativa es que la diferencia de las medias no es igual a cero.

En la prueba de hiótesis, se desea aceptar o rechazar la hipotesis nula con ciertos intervalos de confianza. Dado que probamos la diferencias entre ambas medias, el intervalo de confianza especifica el intervalo de los valores de esta diferencia.

Alternativametne podemos usar la prueba implementada en rstatixs:

pin_Ade_ttest <- pin_Ade %>%
  t_test(culmen_length_mm ~ island, var.equal = TRUE) %>%
  print()
## # A tibble: 1 x 8
##   .y.              group1 group2    n1    n2 statistic    df     p
## * <chr>            <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
## 1 culmen_length_mm Biscoe Dream     44    56     0.950    98 0.344

Ahora con la información generada con la prueba, es posible agregar la anotacion a la gráfica

p.val = pin_Ade_ttest$p

ggplot(pin_Ade, aes(x = island, y = culmen_length_mm, fill = island))+
  geom_boxplot(outlier.shape = "")+
  geom_point(position = position_jitter(0.1), col = "grey50")+
  theme_light()+
  labs(x = "Isla",
       y = "longitud del culmen (mm)",
       title = "Longitud del culmen por isla en el pingüino Adelie",
       fill = "Isla", col = "Isla",
       caption = paste0("t-student; p = ", p.val))

2.3: ¿ Hay diferencia en la longitud de la aleta entre pingünio Adelie y Gentoo?

pin_Aleta <- pin %>%
  filter(species %in% c("Adelie", "Gentoo"))

Histograma de frecuencias

ggplot(pin_Aleta, aes(x = flipper_length_mm, fill = species, col = species))+
  geom_histogram(alpha = 0.3)+
  labs(y = "Frecuencia",
       x = "longitud de la aleta (mm)",
       title="Grafico de densidad", 
       fill = "especie", 
       col = "especie")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing non-finite values (stat_bin).

qqplot

ggplot(pin_Aleta, aes(sample = flipper_length_mm, col = species))+
  stat_qq()+
  stat_qq_line()+
  facet_grid(.~ species)
## Warning: Removed 2 rows containing non-finite values (stat_qq).
## Warning: Removed 2 rows containing non-finite values (stat_qq_line).

Prueba de shapiro utilizando el paquete RSTATIX

pin_Aleta_shap <- pin_Aleta %>%
  group_by(species) %>%
  shapiro_test(flipper_length_mm) %>%
  print()
## # A tibble: 2 x 4
##   species variable          statistic       p
##   <chr>   <chr>                 <dbl>   <dbl>
## 1 Adelie  flipper_length_mm     0.993 0.720  
## 2 Gentoo  flipper_length_mm     0.962 0.00162

Distribucion de las varianzas

leveneTest(flipper_length_mm ~ species, data = pin_Aleta)
## 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   0.157 0.6922
##       272

Boxplot

ggplot(pin_Aleta, aes(x = species, y = flipper_length_mm, fill = species))+
  geom_boxplot()+
  geom_point(position = position_jitter(0.1), col = "grey50")+
  theme_light()+
  labs(x = "Especie",
       y = "longitud de la aleta (mm)",
       title = "Longitud de la aleta por especie",
       fill = "Especie", col = "Especie")
## Warning: Removed 2 rows containing non-finite values (stat_boxplot).
## Warning: Removed 2 rows containing missing values (geom_point).

Hacemos la prueba T de student

t.test(flipper_length_mm ~ species, data = pin_Aleta, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  flipper_length_mm by species
## t = -34.415, df = 272, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -28.79125 -25.67545
## sample estimates:
## mean in group Adelie mean in group Gentoo 
##             189.9536             217.1870

t.student con rstatixs

pin_Aleta_ttest <- pin_Aleta %>%
  t_test(flipper_length_mm ~ species, var.equal = TRUE)
pin_Aleta_ttest
## # A tibble: 1 x 8
##   .y.               group1 group2    n1    n2 statistic    df         p
## * <chr>             <chr>  <chr>  <int> <int>     <dbl> <dbl>     <dbl>
## 1 flipper_length_mm Adelie Gentoo   152   124     -34.4   272 4.21e-101

La prueba nos dice que las medias de ambas poblaciones son diferentes. ***

3: Comparacion de medias poblacionales dependientes (pareadas)

Dos medias son dependientes o pareadas cuando proceden de grupos o muestras dependientes, esto es, cuando existe una relación entre las observaciones de las muestras. Este escenario ocurre a menudo cuando los resultados se generan a partir de los mismos individuos bajo condiciones distintas.

Las pruebas pareadas tienen la ventaja frente a los independientesde que se puede controlar mejor la variación no sistematica, ya que se bloquean al estar examinando los mismos individuos dos veces, no dos grupos de individuos distintos.

3.1: Datos

Se utilizaran los datos de consumo de oxígeno (\(MO_2\)) del abulon Haliotis fulgens aclimatado a 18°C y expuesto a un incremento agudo en la temperatura a 30°C.

mo2 <- read_csv("data/abalone_MO2.csv")

mo2 <- mo2 %>%
  mutate(Temperatura = factor(Temperatura))

3.2: ¿Hay diferencia en el consumo de oxígeno (\(MO_2\)) tras la exposición a la temperatura?

Exploración de los datos con geom_histogram

ggplot(mo2, aes(x = MO2,fill=Temperatura, col = Temperatura))+
  geom_histogram()+
  labs(x = "consumo de oxígeno (MO2)")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Verificar la normalidad con gráfico de qq seguido con una prueba de Shapiro

ggplot(mo2, aes(sample= MO2, col = Temperatura)) +
  stat_qq()+
  stat_qq_line()+
  facet_grid(.~ Temperatura)

mo2_shapiro <- mo2 %>%
  group_by(Temperatura) %>%
  shapiro_test(MO2)
mo2_shapiro
## # A tibble: 2 x 4
##   Temperatura variable statistic      p
##   <fct>       <chr>        <dbl>  <dbl>
## 1 18          MO2          0.829 0.0329
## 2 32          MO2          0.879 0.128

Ambas gráficas y la prueba indica que uno de los grupos no cumple con el supuesto de normalidad. Dado que el tamaño de la muestra es pequeño, hay que evaluar cual es el mejor procedimiento.

Una opción es analizar cada una de las muestras y evaluar si alguno puede ser considerado como outlier. Otra opcion es realizar una transformación de los datos.

A continuación, se tranformaran los datos de MO2 a logaritmo

mo2_log <- mo2 %>%
mutate(MO2 = log(MO2))

Verificar la normalidad de los datos transformados con gráfico de qq seguido con una prueba de Shapiro

ggplot(mo2_log, aes(sample= MO2, col = Temperatura)) +
  stat_qq()+
  stat_qq_line()+
  facet_grid(.~ Temperatura)+
  labs(title = "QQplot de datos transformados")

Prueba de Shapiro para los datos transformados

mo2_log_shapiro <- mo2_log %>%
  group_by(Temperatura) %>%
  shapiro_test(MO2)
mo2_log_shapiro
## # A tibble: 2 x 4
##   Temperatura variable statistic     p
##   <fct>       <chr>        <dbl> <dbl>
## 1 18          MO2          0.952 0.696
## 2 32          MO2          0.929 0.441

Prueba de levene para varianzas

mo2_log_leven <- mo2_log %>%
  levene_test(MO2 ~ Temperatura)
mo2_log_leven
## # A tibble: 1 x 4
##     df1   df2 statistic      p
##   <int> <int>     <dbl>  <dbl>
## 1     1    18      4.42 0.0499

Prueba t de student pareada

mo2_log_ttest <- mo2_log %>%
  t_test(MO2 ~ Temperatura, paired = TRUE, var.equal = TRUE)
mo2_log_ttest
## # A tibble: 1 x 8
##   .y.   group1 group2    n1    n2 statistic    df         p
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>     <dbl>
## 1 MO2   18     32        10    10     -8.17     9 0.0000187

Y finalmente graficamos los datos

ggplot(mo2, aes(x = Temperatura, y = MO2, fill = Temperatura))+
  geom_boxplot(outlier.shape = "")+
  geom_point(col = "grey35")+
  geom_line(aes(group = Individuo), col = "grey35")+
  labs(title = "Cambio en el consumo de oxígeno tras un incremento en la temperatura")