Resampling: test de permutación, simulación de Monte Carlo y Bootstrapping

Más sobre ciencia de datos en cienciadedatos.net

• T-test
• Test Wilcoxon-Mann-Whitney
• Test de los rangos con signo de Wilcoxon
• Test estadístico para proporciones
• ANOVA

Versión PDF: Github

Introducción

Los métodos basados en remuestreo (resampling) se engloban dentro de los test no paramétricos ya que no requieren de ninguna asunción sobre la distribución de la o las poblaciones estudiadas. Son por lo tanto, una alternativa a los test paramétricos (t-test, anova…) cuando no se satisfacen las condiciones del Teorema del Límite Central o cuando se quiere hacer inferencia sobre un parámetro distinto a la media. A lo largo de este documento se describe la idea general de los métodos de resampling más utilizados.

Test de permutaciones

El test de permutaciones es un test de significancia estadística para el estudio de diferencias entre grupos. Fue desarrollado por Ronald Fisher y E.J.G. Pitman en 1930. La distribución del estadístico estudiado (media, mediana…) se obtiene calculando el valor de dicho estadístico para todas las posibles reorganizaciones de las observaciones en los distintos grupos. Dado que implica calcular todas las posibles situaciones, se trata de un test exacto.

Para ilustrar la idea de los test de permutación se plantea un experimento simple. Supóngase un conjunto de sujetos que se distribuye en dos grupos, A y B, de tamaños \(n_{A}\) y \(n_{B}\), cuyas medias muestrales tras realizar el experimento resultan ser \(\bar{x}_{A}\) y \(\bar{x}_{B}\). Se desea determinar si existe una diferencia significativa entre la media de los dos grupos, o lo que es lo mismo, comprobar si hay evidencias en contra de la hipótesis nula de que la diferencia observada es debida, únicamente, a la asignación aleatoria de los sujetos a los dos grupos y que ambas muestras proceden realmente de la misma población. Para solucionarlo se aplica un test de permutación:

---

• En primer lugar, se calcula la diferencia entre las medias de los dos grupos, lo que se conoce como diferencia observada \((Dif_{observada})\).

• Todas las observaciones se combinan juntas sin tener en cuenta el grupo al que pertenecían.

• Se calculan todas las posibles permutaciones en las que las observaciones pueden ser distribuidas en dos grupos de tamaño \(n_{A}\) y \(n_{B}\).

• Para cada permutación, se calcula la diferencia entre medias \((Dif_{calculada})\). El conjunto de valores calculados forman la distribución exacta de las posibles diferencias siendo cierta la hipótesis nula. A esta distribución se le conoce como permutation distribution of the mean difference.

• El p-value de dos colas se calcula como la proporción de permutaciones en las que, el valor absoluto de la diferencia calculada, es mayor o igual al valor absoluto de la diferencia observada.

---

En este ejemplo se emplea como estadístico la media pero podría ser cualquier otro (mediana, varianza…)

La condición necesaria para un test de permutaciones se conoce como exchangeability, según la cual, todas las posibles permutaciones tienen la misma probabilidad de ocurrir siendo cierta la hipótesis nula. Las conclusiones de un test de permutación solo son aplicables a diseños de tipo experimental, en los que, tras haber elegido los sujetos del estudio, se realiza una asignación aleatoria a los diferentes grupos.

Los test de permutación son test de significancia y por lo tanto se emplean para calcular p-values, no para intervalos de confianza.

Simulación de Monte Carlo

En los test de permutación, el p-value obtenido es exacto ya que se calculan todas las posibles permutaciones de las observaciones. Esto resulta complicado o imposible cuando el tamaño muestral es mediano o grande. La simulación de Monte Carlo o randomization test consiste en realizar el test empleando únicamente una muestra aleatoria de todas las posibles permutaciones, evitando así tener que calcularlas todas. Al no calcularse todas las posibles permutaciones, no es un test exacto sino aproximado.

El método de Monte Carlo no es insesgado, lo que hace necesaria una corrección (Davison and Hinkley) tal que: \[p_{value}=\frac{r+1}{k+1}\] siendo r el número de permutaciones igual o más extremas que el estadístico observado y k el número de permutaciones empleadas en la simulación.

A pesar de que la corrección es mínima cuando el número de simulaciones es alto, presenta la ventaja de que, si el valor observado es mayor que cualquiera de los calculados, el p-value es muy bajo pero no cero. Por facilidades de cálculo se suele emplear un número de permutaciones terminado en 9 (4999, 9999).

Bootstrapping

El escenario ideal para realizar inferencia estadística sobre una población es disponer de infinitas (o una gran cantidad) de muestras de dicha población, sin embargo, en la práctica, raramente es posible. Si solo se dispone de una muestra, y esta es representativa de la población, los valores de la variable aleatoria en la muestra aparecen aproximadamente con la misma proporción que en la población. El método de Bootstrapping se basa en generar nuevas pseudomuestras del mismo tamaño que la muestra real mediante muestreo con reemplazo (sampling with replacement) de las observaciones. Si la muestra original es representativa de la población, la distribución del estadístico calculada a partir de las pseudomuestras (bootstraping distribution) se asemejará a la distribución muestral que se obtendría si se pudiera acceder a la población para generar nuevas muestras. Esta idea fue desarrollado por Bradley Efron en 1979.

Dado que se podrían generar infinitas nuevas muestras mediante el sampling with replacement, el método de bootstrapping emplea únicamente una cantidad determinada por el usuario. Hace uso, por lo tanto, de la simulación de Monte Carlo.

El bootstrapping no asume una asignación aleatoria de los grupos, sino que las muestras han sido obtenidas aleatoriamente de la o las poblaciones. Se aplica por lo tanto en diseños muestrales, no experimentales. Esta es la diferencia clave respecto a los test de permutación.

El método bootstrapping se puede emplear para:

• Calcular intervalos de confianza de un parámetro poblacional: se emplea el sampling with replacement a partir de la muestra original.

• Calcular significancia estadística (p-value) para la diferencia entre dos poblaciones: Si bien este uso se asemeja al de los test de permutación, no es igual. Los test de permutación contrastan la hipótesis nula de que las muestras pertenecen a una misma población (distribución) mediante el estudio de las diferencias debidas a la asignación aleatoria de los grupos. El método de bootstrapping también contrasta la hipótesis de que ambas muestras proceden de la misma población (distribución), pero lo hace mediante el estudio de las diferencias debidas al muestreo aleatorio. Se aplica por lo tanto a estudios en los que no ha habido una asignación aleatoria a los grupos previa realización de los experimentos. Los pasos a seguir son: se mezclan las observaciones de ambas muestras (pool), se emplea el sampling with replacement sobre este pool para generar una nueva pseudomuestra del mismo tamaño, se separan las observaciones de la pseudomuestra en dos grupos de igual tamaño a los originales y se calcula la diferencia del estadístico entre ambas. El proceso se repite múltiples veces generando así la distribución de las diferencias esperadas debido al muestreo aleatorio. El p-value de dos colas se calcula como la proporción de pseudomuestras en las que el valor absoluto de la diferencia calculada es mayor o igual al valor absoluto de la diferencia observada.

• Calcular intervalos de confianza para la diferencia entre dos poblaciones: para esta finalidad, se considera como hipótesis nula que las observaciones proceden de dos poblaciones distintas. Se emplea el sampling with replacement con la observaciones de cada muestra (sin mezclarlas) para generar dos nuevas pseudomuestras independientes y se calcula la diferencia del estadístico. Este proceso se repite múltiples veces, generando así la distribución que se obtendría si se obtuviesen cada vez dos muestras, cada una de su respectiva población, y se calculara la diferencia. La distribución resultante está centrada en la verdadera diferencia entre las poblaciones.

El método de bootstrapping se subdivide en paramétrico o no paramétrico dependiendo si se considera que la distribución generada sigue un determinado modelo teórico.

Comparación de los test de permutación y métodos bootstrapping

Tanto los test de permutación como los test de bootstrap se pueden emplear para estudiar diferencias entre grupos. Existe una lista muy extensa de referencias en las que se debate cuál de los dos métodos es el más adecuado. En general, todas ellas concluyen en que el método más adecuado depende del objetivo de la inferencia, y a su vez, el objetivo de la inferencia determina que diseño de estudio a seguir. La siguiente tabla contiene los diferentes tipos de diseños que se pueden emplear para comparar dos grupos y el tipo de inferencia (conclusiones) que se puede realizar en cada uno:

Muestreo aleatorio	Asignación de grupos aleatoria	Objetivo de la inferencia	Permite determinar causalidad
Sí	No	Población	No
No	Sí	Muestra	Sí
Sí	Sí	Población y Muestra	Sí

La principal diferencia entre ambos métodos aparece cuando se emplean para calcular p-values. Los test de significancia (cálculo de p-value) se basan en la hipótesis nula de que todas las observaciones proceden de la misma población. El objetivo del test es determinar si la diferencia observada entre los grupos se debe a un determinado factor (tratamiento) o solo a la variabilidad esperada por la naturaleza de un proceso aleatorio. Cuando la aleatoriedad se debe a la asignación de los sujetos (supuestamente iguales) a los distintos grupos, se emplean los test de permutación. La estructura de un experimento que puede analizarse mediante test de permutación es: selección de sujetos del estudio, asignación aleatoria a diferentes grupos, aplicación de o los “tratamientos” y comparación de resultados. Los test de permutación responden a la pregunta ¿Cuánta variabilidad se espera en un determinado test estadístico debido únicamente a la aleatoriedad de las asignaciones, si todos los sujetos proceden realmente de una misma población? Cuando se compara la media entre dos grupos, la pregunta anterior equivale a ¿Qué diferencia entre medias cabe esperar dependiendo de cómo se distribuyan los sujetos en los dos grupos, si todos proceden de una misma población? (Aun siendo todos de una misma población, dado que no serán exactamente idénticos, habrá pequeñas diferencias dependiendo de cómo se agrupen).

El bootstraping como test de significancia se emplea cuando la aleatoriedad es debida al proceso de obtención de las muestras y no a la asignación en grupos. Responden a la pregunta ¿Cuánta variabilidad se espera en un determinado estadístico debido únicamente al muestreo aleatorio, si todos los sujetos proceden realmente de una misma población? Debido a las pequeñas diferencias entre los individuos de una población, si se extraen dos muestras aleatorias de ella y se comparan, no van a ser exactamente iguales, además, esta diferencia será distinta para cada par de muestras aleatorias extraídas. La estructura de un experimento que puede analizarse mediante bootstrapping es: Se obtienen dos muestras aleatorias de dos poblaciones y se comparan.

En conclusión, aunque ambos test pueden emplearse para calcular p-values, sus aplicaciones no se solapan. Los test de permutación y randomization se emplean para diseños experimentales, mientras que el bootstraping se emplea para diseños muestrales.

El método de bootstrapping puede emplearse, además, para generar intervalos de confianza de la verdadera diferencia de un parámetro entre dos poblaciones, mientras que los test de permutación no. El porqué reside en que la simulación que se realiza en el bootstrapping. Cuando se emplea por separado para cada muestra, se genera una aproximación de la verdadera distribución del estadístico que se está estudiando, lo que permite generar intervalos de confianza que acoten su valor con una determinada seguridad. Cuando se quiere calcular un p-value, se tiene que comparar el valor observado del estadístico con la distribución que cabría esperar de él bajo la hipótesis nula, no respecto a su verdadera distribución. Los test de permutación y los test de bootstrapping (aplicado al conjunto de observaciones mezcladas), generan la distribución esperada si se cumple la \(H_0\).

Es importante tener en cuenta que ninguno de estos métodos está al margen de los problemas que implica tener muestras pequeñas.

Test de permutación exacto

Este tipo de test solo es posible aplicarlo cuando el tamaño de las muestras es limitado, ya que implica generar todas las posibles permutaciones de los datos y calcular para cada una de ellas el estadístico de interés. No he encontrado en R ningún paquete que lo realice

Test de permutación con simulación de Monte Carlo

Ejemplo1

Supóngase un estudio que pretende determinar si la participación en actividades extraescolares aumenta la capacidad empática de los estudiantes. Para ello, el colegio ofrece un programa voluntario en el que cada participante se designa de forma aleatoria a un grupo “control” que no recibe clases extraescolares o a un grupo “tratamiento” que sí las recibe. A final del año, todos los sujetos del estudio realizan un examen que determina su capacidad empática. En vista de los resultados ¿Se puede considerar que las clases extraescolares tienen un impacto en cómo se relacionan socialmente (en promedio) los estudiantes?

Ejemplo obtenido del libro Comparing Groups Randomization and Bootstrap Methods Using R. Los datos ya no están disponibles en la página web del libro, pero puede encontrarse en Github.

library(tidyverse)
datos <- read_csv("./datos/AfterSchool.csv")
datos <- datos %>% select(Treatment, Delinq)
datos <- datos %>% rename(grupo = Treatment, puntuacion = Delinq)
datos <- datos %>% mutate(grupo =  recode(grupo, `0` = "control",
                                                 `1` = "tratamiento"))
head(datos)

El diseño experimental del estudio emplea una asignación aleatoria de los sujetos a dos grupos (tratamiento y control) para llevar a cabo el experimento. Esta aleatorización en la asignación implica que, en promedio, los dos grupos son iguales para todas las características, por lo que la única diferencia entre ellos es si reciben o no el tratamiento. Determinar si la diferencia observada es significativa, equivale a preguntarse cómo de probable es obtener esta misma diferencia si el tratamiento no tiene efecto, o lo que es lo mismo, determinar si la diferencia observada es mayor de lo que cabría esperar debido únicamente a la variabilidad producida por la formación aleatoria de los grupos.

El diseño experimental de asignación aleatoria a grupos permite obtener conclusiones de tipo causa efecto. Sin embargo, dado que la selección de los sujetos no ha sido aleatoria, sino por voluntarios, no son extrapolables a toda la población.

Analisis exploratorio de los datos

ggplot(data = datos, aes(x = grupo, y = puntuacion, colour = grupo)) +
    geom_boxplot(outlier.shape = NA) +
    geom_jitter(width = 0.25) +
    scale_color_manual(values = c("gray50", "orangered2")) +
    coord_flip() +
    theme_bw() +
    theme(legend.position = "none")

ggplot(data = datos, aes(x = puntuacion, fill = grupo)) +
    geom_density(alpha = 0.5) +
    scale_fill_manual(values = c("gray50", "orangered2")) +
    theme_bw() +
    theme(legend.position = "bottom")

La representación gráfica muestra una clara distribución asimétrica de los valores, con una cola pronunciada hacia la derecha en ambos grupos. No se observa una diferencia evidente entre grupos.

datos %>% group_by(grupo) %>% summarise(media = mean(puntuacion),
                                        sd = sd(puntuacion),
                                        mediana = median(puntuacion))

Las medias de localización (media y mediana) de ambos grupos son similares. Lo mismo ocurre con su dispersión.

Test de permutaciones

Dado que las observaciones no se distribuyen de forma normal, los test de hipótesis paramétricos no son adecuados. Como alternativa, se recurre a un test no paramétrico basado en resampling. Al tratarse de un diseño experimental en el que se ha partiendo de un conjunto de sujetos y se han asignado aleatoriamente a cada grupo, el test adecuado es el de permutaciones.

En primer lugar se calcula la diferencia entre las medias de ambos grupos (diferencia observada).

media_control <- datos %>% filter(grupo == "control") %>%
                           pull(puntuacion) %>%
                           mean()
media_tratamiento <- datos %>% filter(grupo == "tratamiento") %>%
                           pull(puntuacion) %>%
                           mean()
dif_obs <- media_control - media_tratamiento
dif_obs

## [1] 1.706636

Determinar si la diferencia observada es significativa equivale a preguntarse cómo de probable es obtener esta diferencia si el tratamiento no tiene efecto y los estudiantes se han asignado de forma aleatoria en cada grupo. Para poder obtener la probabilidad exacta, se necesita generar todas las posibles permutaciones en las que, 356 sujetos, pueden repartirse en dos grupos, y calcular la diferencia de medias para cada una. El número de permutaciones posibles es muy elevado, \((3.93x10^{105})\), por lo que se recurre a una simulación de Monte Carlo. El algoritmo seguido es:

---

1. Se reasignan aleatoriamente los sujetos a cada uno de los grupos manteniendo el tamaño original de cada uno.

2. Se calcula la diferencia de sus medias y se almacena el valor. Existen múltiples formas de realizar las permutaciones, lo importante es que mimeticen una asignación aleatoria de los grupos manteniendo el tamaño original.

3. Se repiten los pasos 1 y 2 n veces.

4. Se calcula el p_value como la proporción de permutaciones en las que, la diferencia absoluta obtenida, es igual o mayor a la observada.

---

distribucion_permut <- rep(NA, 9999)
n_control <- datos %>% filter(grupo == "control") %>% nrow()
n_tratamiento <- datos %>% filter(grupo == "tratamiento") %>% nrow()

for (i in 1:9999) {
  # mezclado aleatorio de las observaciones
  datos_aleatorizados <- sample(datos$puntuacion) 
  mean_control <- mean(datos_aleatorizados[1:n_control])
  mean_tratamiento <- mean(datos_aleatorizados[n_control + 1:n_tratamiento])
  distribucion_permut[i] <- mean_control - mean_tratamiento
}

head(distribucion_permut)

## [1] -0.2612364 -1.0629620 -0.3341206  0.7591416 -0.2612364 -0.1154681

Los datos simulados forman lo que se conoce como distribución de permutaciones o de Monte Carlo, que representa la variación esperada en la diferencia de medias debido únicamente a la asignación aleatoria de grupos.

ggplot(data = data.frame(permutaciones = distribucion_permut),
      aes(permutaciones)) +
  geom_histogram(aes(y = ..density..), alpha = 0.4, colour = "white") + 
  geom_line(aes(y = ..density..), stat = 'density', colour = "red") +
  geom_vline(xintercept = mean(distribucion_permut)) +
  geom_vline(xintercept = dif_obs, colour = "blue") +
  geom_vline(xintercept = -dif_obs, colour = "blue") +
  labs(title = "Distribución de permutaciones", x = "diferencia de medias") +
  theme_bw()

summary(distribucion_permut)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -4.342749 -0.698541 -0.042584  0.005852  0.759142  4.330465

sd(distribucion_permut)

## [1] 1.043649

Como era de esperar, si el tratamiento no es efectivo, la diferencia media entre grupos es muy próxima a cero (línea vertical negra). La desviación típica (1.044) de la distribución de permutaciones indica que la variabilidad debida a la asignación aleatoria de los sujetos es muy pequeña.

Finalmente, se calcula la probabilidad (p-value) de obtener diferencias igual o más extremas que la observada (líneas verticales azules) con y sin corrección de continuidad.

# Sin corrección
p_value = (sum(abs(distribucion_permut) > abs(dif_obs)))/9999
p_value

## [1] 0.1064106

# Con corrección
p_value_corregido = ((sum(abs(distribucion_permut) > abs(dif_obs))) + 1)/(9999 + 1)
p_value_corregido

## [1] 0.1065

Conclusión

Los 356 sujetos del estudio fueron asignados de forma aleatoria a un grupo control (n=187) o a un grupo tratamiento (n=169) que asistió a clases extra escolares. Un test de permutación se empleó para determinar si existía una diferencia significativa en la capacidad empática promedio entre ambos grupos. El p-value fue calculado mediante una simulación de Monte Carlo con 9999 permutaciones usando la corrección de continuidad sugerida por Davison and Hikley(1997). El p-value obtenido muestra una evidencia muy débil en contra de la hipótesis nula de que el tratamiento no tiene efecto, sugiriendo que asistir a clases extra escolares no mejora la capacidad empática para los estudiantes que formaron parte del experimento. Siendo estadísticamente estrictos, no se puede extrapolar a la población de estudiantes ya que la selección de sujetos no fue aleatoria.

A modo de comprobación, siendo los tamaños muestrales de más de 30 observaciones, el test paramétrico t-test debería dar un resultado similar.

t.test(puntuacion ~ grupo, data = datos, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  puntuacion by grupo
## t = 1.6379, df = 354, p-value = 0.1023
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.3425852  3.7558562
## sample estimates:
##     mean in group control mean in group tratamiento 
##                  50.72559                  49.01896

Los p-values son similares: Bootstrapping = 0.1065, t-test = 0.1023

Ejemplo2

Una investigación quiere determinar si las clases de repaso tienen un efecto significativo en el resultado de los exámenes de los estudiantes. Para ello, un conjunto de estudiantes se reparte al azar en dos grupos, uno que asiste a clases de repaso y otro que no, y se evalúa su conocimiento con un examen. ¿Existe una diferencia significativa en el promedio entre ambos grupos? Solucionarlo con t-test y test de permutación si se cumplen sus respectivas condiciones. Ejemplo del libro Bootstrap Methods ans Permutation Test by Tim Hestemberg.

grupo <- c("clases repaso", "clases repaso", "clases repaso","clases repaso",
           "clases repaso", "clases repaso", "clases repaso", "clases repaso",
           "clases repaso", "clases repaso", "clases repaso", "clases repaso",
           "clases repaso", "clases repaso", "clases repaso", "clases repaso",
           "clases repaso", "clases repaso", "clases repaso", "clases repaso",
           "clases repaso", "control", "control", "control", "control", "control",
           "control", "control", "control", "control", "control", "control",
           "control", "control", "control", "control", "control", "control",
           "control", "control", "control", "control", "control", "control")
resultado <- c(24, 43, 58, 71, 43, 49, 61, 44, 67, 49,  53, 56, 59, 52, 62,  54,
               57, 33, 46, 43, 57, 26, 62, 37, 42, 43, 55, 54, 20, 85, 33, 41, 19,
               60,  53, 42, 46, 10, 17, 28, 48, 37, 42, 55)
datos <- data_frame(grupo, resultado)
head(datos)

Exploración de los datos

El t-test es un test estadístico válido para comparar medias poblacionales siempre y cuando se cumpla que los datos son independientes y que las observaciones proceden de distribuciones normales (son robustos frente a cierta asimetría si el tamaño de ambas muestras es mayor de 30).

datos %>% group_by(grupo) %>%  
          summarise(media = mean(resultado), mediana = median(resultado),
                    sd = sd(resultado))

ggplot(data = datos, aes(x = grupo, y = resultado, colour = grupo)) +
    geom_boxplot(outlier.shape = NA) +
    geom_jitter(width = 0.2) +
    scale_color_manual(values = c("gray50", "orangered2")) +
    theme_bw() +
    theme(legend.position = "none")

ggplot(data = datos, aes(x = resultado, fill = grupo)) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  geom_vline(xintercept = mean(datos%>%filter(grupo=="clases repaso")%>%pull(resultado)),
             color = "gray50", linetype = "dashed", size = 1) +
  geom_vline(xintercept = mean(datos%>%filter(grupo=="control")%>%pull(resultado)),
             color = "orangered2", linetype = "dashed", size = 1) +
  theme_bw() +
  theme(legend.position = "bottom")

cutom_qqnorm <- function(grupo, df){
  qqnorm(df$resultado, main = grupo, pch = 19)
  qqline(df$resultado)
}

datos %>% group_by(grupo) %>% nest() %>% 
          mutate(plot = walk2(.x = grupo, .y = data, .f = cutom_qqnorm))

custom_shapiro_test <- function(df){
  return(shapiro.test(df$resultado)$p.value)
}

datos %>% group_by(grupo) %>% nest() %>% 
          mutate(pvalue_shapiro = map_dbl(.x = data, .f = custom_shapiro_test))

Ambas muestras se distribuyen aproximadamente de forma normal con varianzas similares, por lo que se puede aplicar un t-test.

t.test(resultado ~ grupo, data = datos, altermative = "two.sided", mu = 0,
       var.equal = TRUE, paired = FALSE , conf.level = 0.95)

## 
##  Two Sample t-test
## 
## data:  resultado by grupo
## t = 2.2666, df = 42, p-value = 0.02863
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   1.091253 18.817650
## sample estimates:
## mean in group clases repaso       mean in group control 
##                    51.47619                    41.52174

Solución mediante test de permutación

distribucion_permut <- rep(NA, 9999)
n_control <- datos %>% filter(grupo == "control") %>% nrow()
n_tratamiento <- datos %>% filter(grupo == "clases repaso") %>% nrow()

for (i in 1:9999) {
  # mezclado aleatorio de las observaciones
  datos_aleatorizados <- sample(datos$resultado) 
  mean_control <- mean(datos_aleatorizados[1:n_control])
  mean_tratamiento <- mean(datos_aleatorizados[n_control + 1:n_tratamiento])
  distribucion_permut[i] <- mean_control - mean_tratamiento
}

media_control <- datos %>% filter(grupo == "control") %>%
                           pull(resultado) %>%
                           mean()
media_tratamiento <- datos %>% filter(grupo == "clases repaso") %>%
                           pull(resultado) %>%
                           mean()
dif_obs <- media_control - media_tratamiento
dif_obs

## [1] -9.954451

ggplot(data = data.frame(permutaciones = distribucion_permut),
      aes(permutaciones)) +
  geom_histogram(aes(y = ..density..), alpha = 0.4, colour = "white") + 
  geom_line(aes(y = ..density..), stat = 'density', colour = "red") +
  geom_vline(xintercept = mean(distribucion_permut)) +
  geom_vline(xintercept = dif_obs, colour = "blue") +
  geom_vline(xintercept = -dif_obs, colour = "blue") +
  labs(title = "Distribución de permutaciones", x = "diferencia de medias") +
  theme_bw()

summary(distribucion_permut)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -17.51553  -3.21325  -0.11594  -0.05698   3.16356  15.46170

sd(distribucion_permut)

## [1] 4.619516

Cálculo de p-values

# Sin corrección
p_value = (sum(abs(distribucion_permut) > abs(dif_obs)))/9999
p_value

## [1] 0.02770277

# Con corrección
p_value_corregido = ((sum(abs(distribucion_permut) > abs(dif_obs))) + 1)/(9999 + 1)
p_value_corregido

## [1] 0.0278

Conclusión

Ambos p-values, el obtenido por t-test y por test de permutación, se asemejan mucho. Esto es debido a que, como muestra la distribución de permutaciones, la diferencia de medias se distribuye de forma normal, por lo que los t-test son válidos. En caso de que no fuese así, el test de permutación es una mejor opción. Existen evidencias suficientes para rechazar la hipótesis nula y concluir que sí existen evidencias significativas de que la media de los dos grupos es distinta. Las clases de repaso sí influyen en los resultados de los alumnos.

En caso de querer saber si el t-test es adecuado, una forma de comprobarlo es generar la permutation distribution y verificar si se distribuye de forma normal.

Test de permutación para comparar varianzas

En la mayoría de casos, la comparación entre grupos se centra en estudiar diferencias en la posición de las distribuciones (media, mediana…), sin embargo, también puede ser de interés comparar si las varianzas de dos grupos son iguales. Para este tipo de estudios, se puede emplear un test de permutación en el que el estadístico empleado es la varianza.

Ejemplo

Supóngase el estudio de un departamento médico que quiere evaluar la si un fármaco aumenta el razonamiento matemático de los estudiantes. Se sabe que no existe diferencia en el promedio de las puntuaciones obtenidas en el examen, pero se quiere evaluar si hay diferencia en la variabilidad. Esto es interesante porque, aunque el fármaco no sea capaz de incrementar el promedio de capacidad de razonamiento, podría disminuir la diferencia entre sujetos

La hipótesis nula considera que ambos grupos son iguales por lo que: \[\sigma^2_{control} = \sigma^2_{tratamiento}\] \[\sigma^2_{control} - \sigma^2_{tratamiento} = 0\]

datos <- read_csv("./datos/AfterSchool.csv")
datos <- datos %>% select(Treatment, Delinq)
datos <- datos %>% rename(grupo = Treatment, puntuacion = Delinq)
datos <- datos %>% mutate(grupo =  recode(grupo, `0` = "control",
                                                 `1` = "tratamiento"))

La diferencia observada en las varianzas es:

datos %>% group_by(grupo) %>% summarise(varianza = var(puntuacion))

var_control <- datos %>% filter(grupo == "control") %>%
                         pull(puntuacion) %>%
                         var()
var_tratamiento <- datos %>% filter(grupo == "tratamiento") %>%
                             pull(puntuacion) %>%
                             var()
dif_obs <- var_control - var_tratamiento
dif_obs

## [1] 30.15232

Mediante permutaciones se obtiene la distribución de la diferencia de varianza esperada únicamente debido a la asignación aleatoria de las observaciones a cada grupo, siendo cierta la hipótesis nula.

distribucion_permut <- rep(NA, 9999)
n_control <- datos %>% filter(grupo == "control") %>% nrow()
n_tratamiento <- datos %>% filter(grupo == "tratamiento") %>% nrow()

for (i in 1:9999) {
  datos_aleatorizados <- sample(datos$puntuacion) 
  var_control <- var(datos_aleatorizados[1:n_control])
  var_tratamiento <- var(datos_aleatorizados[n_control + 1:n_tratamiento])
  distribucion_permut[i] <- var_control - var_tratamiento
}

ggplot(data = data.frame(permutaciones = distribucion_permut),
      aes(permutaciones)) +
  geom_histogram(aes(y = ..density..), alpha = 0.4, colour = "white") + 
  geom_line(aes(y = ..density..), stat = 'density', colour = "red") +
  geom_vline(xintercept = mean(distribucion_permut)) +
  geom_vline(xintercept = dif_obs, colour = "blue") +
  geom_vline(xintercept = -dif_obs, colour = "blue") +
  labs(title = "Distribución de permutaciones", x = "diferencia de medias") +
  theme_bw()

summary(distribucion_permut)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -97.7295 -17.8909  -0.1047   0.1294  17.8412  88.7256

sd(distribucion_permut)

## [1] 25.88688

La distribución de las permutaciones muestra que la diferencia media entre varianzas, si el fármaco no tiene efecto, es muy próxima a cero (0.1293638, línea vertical negra). La desviación típica de la distribución de permutaciones indica la variabilidad de la diferencia de varianzas debida a la asignación aleatoria de los sujetos a los diferentes grupos.

Finalmente, se calcula la probabilidad (p-value) de obtener diferencias igual o más extremas que la observada (líneas verticales azules) con y sin corrección de continuidad:

# Sin corrección
p_value = (sum(abs(distribucion_permut) > dif_obs))/9999
p_value

## [1] 0.2488249

# Con corrección
p_value_corregido = ((sum(abs(distribucion_permut) > dif_obs)) + 1)/(9999 + 1)
p_value_corregido

## [1] 0.2489

Conclusión

Los 356 sujetos del estudio fueron asignados de forma aleatoria a un grupo control (n=187) o a un grupo tratamiento (n=169). Un test de permutación se empleó para determinar si existía una diferencia significativa en la varianza de la capacidad de razonamiento matemático entre ambos grupos. El p-value fue calculado mediante una simulación de Monte Carlo con 9999 permutaciones usando la corrección de continuidad sugerida por Davison and Hikley(1997). El p-value obtenido muestra una evidencia muy débil en contra de la hipótesis nula de que el tratamiento no tiene efecto, sugiriendo que tomar el medicamento no reduce la varianza en la capacidad de razonamiento matemático de los estudiantes que formaron parte del experimento.

Boostraping para intervalos de confianza de una población

Ejemplo1

Se dispone de una muestra formada por 30 observaciones de una variable aleatoria continua. Se quiere inferir el valor medio de dicha variable en la población mediante el cálculo de un intervalo de confianza del 95%.

Exploración de los datos

library(ggpubr)

datos <- c(81.372918, 25.700971, 4.942646, 43.020853, 81.690589, 51.195236,
           55.659909, 15.153155, 38.745780, 12.610385, 22.415094, 18.355721,
           38.081501, 48.171135, 18.462725, 44.642251, 25.391082, 20.410874,
           15.778187, 19.351485, 20.189991, 27.795406, 25.268600, 20.177459,
           15.196887, 26.206537, 19.190966, 35.481161, 28.094252, 30.305922)

p1 <- ggplot(data = data_frame(datos), aes(x = datos)) +
      geom_density()  +
      theme_bw()

p2 <- ggplot(data = data_frame(datos), aes(x = datos)) +
      geom_histogram()  +
      theme_bw()

ggarrange(p1, p2)

qqnorm(datos)
qqline(datos)

shapiro.test(datos)

## 
##  Shapiro-Wilk normality test
## 
## data:  datos
## W = 0.85765, p-value = 0.0009003

La representación gráfica y el test de normalidad Shapiro-Wilk muestran que los datos no se distribuyen de forma normal. Esto implica que la aproximación basada en el teorema del límite central para estimar el error estándar \(SE = \frac{sd}{\sqrt{n}}\) deja de ser buena y, por lo tanto, también los intervalos paramétricos basados en la estructura \([parámetro \ estimado \pm t_{\alpha}SE]\) que se generen con esa estimación del \(SE\). Una alternativa para poder calcular intervalos de confianza es emplear bootstrapping.

El método de bootstrapping se basa en generar nuevas pseudomuestras del mismo tamaño que la muestra real, realizando muestreo repetido de las observaciones. Si la muestra original es representativa de la población, la distribución del estadístico calculada a partir de las pseudomuestras (bootstraping distribution) se asemejará a la distribución muestral que se obtendría si se pudiera acceder a la población para generar nuevas muestras. De tal forma que:

• La sd de la bootstraping distribution es un estimador del SE.

• La media de la bootstraping distribution es un estimador del verdadero parámetro poblacional.

Así pues, a partir de la bootstraping distribution se pueden obtener los valores necesarios para crear el intervalo de confianza sin recurrir al teorema del límite central.

Generar la bootstrapping distribution

boot_distribution <- rep(NA, 9999)
for (i in 1:9999) {
    boot_distribution[i] <- mean(sample(x = datos, size = 30, replace = TRUE))
}

Estudio la bootstrapping distribution

Dependiendo de la distribución que tengan los estadísticos obtenidos por bootstraping, se emplea un método distinto para estimar los límites del intervalo de confianza. Cuando la distribución obtenida no es de tipo normal, no se pueden emplear intervalos de confianza basados en la t-distribution ni en percentiles.

Distribución:

ggplot(data = data.frame(boot_distribution), aes(x = boot_distribution)) +
  geom_histogram(aes(y = ..density.. , fill = ..count..)) +
  scale_fill_gradient(low = "#DCDCDC", high = "#7C7C7C") +
  stat_function(fun = dnorm, colour = "firebrick",
                args = list(mean = mean(boot_distribution),
                            sd = sd(boot_distribution))) +
  ggtitle("Histograma + curva normal teórica") +
  theme_bw()

qqnorm(boot_distribution,  xlab = "Quantiles teóricos")
qqline(boot_distribution)

La distribución de los estadísticos obtenidos por bootstraping se distribuyen de forma normal.

Centro:
La media de la bootstraping distribution debe de ser cercana a la media de la muestra inicial a partir de la cual se está generando el bootstrapping. A esta diferencia se le llama bias.

mean(datos) - mean(boot_distribution)

## [1] 0.106298

En este caso el bias es muy pequeño comparado al valor de la media calculada.

Intervalo de confianza

Existen varios métodos para generar intervalos de confianza a partir de una bootstrapping distribution. Los intervalos de tipo t y los basados en percentiles emplean la desviación estándar de la bootstrapping distribution como estimación del SE.

Intervalo basado en t-distribution:

\[[parámetro \ estimado \pm t_{\alpha,df}SE]\]

• media = mean(datos) = 30.9686559

• \(t_{\alpha = 0.05,df=29}\) = qt(p=1-0.05/2, df = 29) = 2.045

• SE = sd(boot_distribution) = 3.309775

• IC 95% = [24.2, 37.74]

Intervalo de confianza basado en percentiles:

# Un IC del 95% debe abarcar desde el percentil 0.025 al 0.975
quantile(x = boot_distribution,probs = c(0.025, 0.975))

##     2.5%    97.5% 
## 24.88947 37.87902

Tanto el intervalo basado en el valor t como el basado en percentiles son únicamente válidos cuando la bootstrapping distribution se asemeja a una normal y el bias es pequeño. El intervalo basado en percentiles no ignora la asimetría de los datos por lo que se suele considerar más adecuado. En caso de no cumplirse estas condiciones, o de que ambos intervalos no se asemejen entre sí, ninguno es fiable.

El intervalo de confianza más adecuado para bootstrapping se conoce como Intervalo Bootstrapping Bias Corrected Acelerated, BCa. En R existen varios paquetes que permiten obtener la distribución de bootstrapping así como el intervalo BCa.

Solución mediante R

La función boot() del paquete boot recibe como argumentos un vector con las observaciones, una función que calcule el estadístico de interés (tiene que ser definida previamente por el usuario) y el número de simulaciones.

library(boot)
media <- function(valores, i) {
    mean(valores[i])
}

boot_distribution <- boot(datos, media, R = 9999)
boot_distribution

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = datos, statistic = media, R = 9999)
## 
## 
## Bootstrap Statistics :
##     original     bias    std. error
## t1* 30.96866 0.01633382    3.287712

Una vez generado el objeto boot, se puede mostrar gráficamente la distribución y calcular los diferentes tipos de intervalos.

plot(boot_distribution)

boot.ci(boot_distribution, conf = 0.95)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 9999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_distribution, conf = 0.95)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (24.51, 37.40 )   (24.12, 36.98 )  
## 
## Level     Percentile            BCa          
## 95%   (24.96, 37.81 )   (25.59, 38.86 )  
## Calculations and Intervals on Original Scale

Otra opción disponible en R es la función bootES() del paquete bootES, que internamente hace una llamada a la función boot() pero sin necesidad de definir la función del estadístico.

library(bootES)
bootES(data = datos, R = 9999, mean)

Ejemplo 2

Supóngase que se quiere estimar el precio de las viviendas en una ciudad a partir de una muestra de 50 propiedades. ¿Cuál es el intervalo de confianza del 95% del precio de las viviendas de la ciudad?

precios <- c(142, 175, 198, 150, 705, 232, 50, 146, 155, 1850, 132, 215,117, 245,
             290, 200, 260, 450, 66, 165, 362, 307, 266, 166, 375,245, 211, 265,
             296, 335, 335, 1370, 256, 148, 988, 324, 216, 684,270, 330, 222, 180,
             257, 253, 150, 225, 217, 570, 507, 190)

Exploración de los datos

ggplot(data = data.frame(precios), aes(x = precios)) +
  geom_histogram(aes(y = ..density.. , fill = ..count..)) +
  scale_fill_gradient(low = "#DCDCDC", high = "#7C7C7C") +
  stat_function(fun = dnorm, colour = "firebrick",
                args = list(mean = mean(precios),
                            sd = sd(precios))) +
  ggtitle("Histograma + curva normal teórica") +
  theme_bw()

qqnorm(precios)
qqline(precios)

La representación gráfica de los datos refleja que la muestra es muy asimétrica.

Si los datos no se distribuyen de forma normal, no se pueden aplicar métodos paramétricos basados en el teorema del límite central. Además, dado que hay asimetría, es recomendable emplear estadísticos más robustos que la media como la mediana o la media truncada (trimmed mean). Por lo general, para los métodos de bootstrapping, se recomienda la media truncada antes que la mediana.

• 25% media truncada: Se define como la media del 50% de los valores centrales. Es decir, la media de los valores que caen entre el primer y el tercer percentil (lo que forma la caja en los diagramas box-plot). En R se calcula mediante mean(x=, trim= 0.25).

La idea de Bootstrapping es que, la bootstrapping distribution, se debe asemejar a la sampling distribution del parámetro estimado (en este caso la trimmed mean). Y que la sd de la bootstrapping distribution se aproxima al SE de la sampling distribution. Es decir, mediante bootstrapping podemos obtener los valores necesarios para calcular el intervalo de confianza sin recurrir a los métodos teóricos basados en teorema del límite central.

Cálculo de la bootstraping distribution

library(boot)
media_truncada <- function(valores, i) {
    mean(valores[i], trim = 0.25)
}
boot_distribution <- boot(precios, media_truncada, R = 9999)
boot_distribution

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = precios, statistic = media_truncada, R = 9999)
## 
## 
## Bootstrap Statistics :
##     original    bias    std. error
## t1* 244.0385 0.5635179    17.46588

Estudio de la bootstrapping distribution

Distribución:

plot(boot_distribution)

La distribución generada no cumple las condiciones de normalidad. Esto implica que los intervalos de confianza obtenidos por percentiles o basados en la t-distribution no son fiables, es más precisa la estimación obtenida por los intervalos BCa.

Centro:

La media de la bootstrapping distribution debe de ser cercana al valor del estadístico estudiado en la muestra inicial a partir de la cual se está generando el bootstrapping. A esta diferencia se le llama bias.

mean(precios, trim = 0.25) - mean(boot_distribution$t)

## [1] -0.5635179

Intervalo de confianza

boot.ci(boot_distribution)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 9999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_distribution)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (209.2, 277.7 )   (207.0, 275.3 )  
## 
## Level     Percentile            BCa          
## 95%   (212.8, 281.1 )   (213.1, 281.7 )  
## Calculations and Intervals on Original Scale

Bootstrapping para cálculo de p-value.

El bootstrapping, como test de significancia para la diferencia entre grupos, se emplea cuando se quiere estudiar si la diferencia entre dos poblaciones es significativa, empleando muestras aleatorias y separadas de cada una de las poblaciones.

Ejemplo1

Un estudio pretende determinar si existe diferencia entre el promedio de los logros académicos de los estudiantes inmigrantes procedentes de México y los procedentes de otros países latinoamericanos. Para ello, se obtiene una muestra aleatoria de estudiantes que llegaron al país entre 2005 y 2006 entre los que hay mexicanos y no mexicanos (otros paises).

Ejemplo obtenido del libro Comparing Groups Randomization and Bootstrap Methods Using R. Los datos ya no están disponibles en la página web del libro, pero puede encontrarse en Github.

Dado que se trata de un diseño muestral, en el que se ha obtenido una muestra aleatoria con individuos de cada uno de los grupos que se van a comparar, el método adecuado es bootstrapping.

Exploración de datos

datos <- read_csv("./datos/LatinoEd.csv")
datos <- datos %>% select(Achieve, Mex)
datos <- datos %>% rename(nota = Achieve, nacionalidad = Mex)
datos <- datos %>% mutate(nacionalidad = recode(nacionalidad,
                                                `1` = "mexicana",
                                                `0` = "otras"))
head(datos)

El promedio del rendimiento académico es de 5.92 puntos menor en estudiantes de origen mexicano que el de otras nacionalidades latinas.

datos %>% group_by(nacionalidad) %>% summarise(media = mean(nota),
                                               mediana = median(nota),
                                               sd = sd(nota),
                                               n = n())

mean_mexicana <- datos %>% filter(nacionalidad == "mexicana") %>%
                           pull(nota) %>%
                           mean()
mean_otras <- datos %>% filter(nacionalidad == "otras") %>%
                           pull(nota) %>%
                           mean()
dif_obs <- mean_mexicana - mean_otras
dif_obs

## [1] -5.921602

La exploración gráfica de los datos indica posibles diferencias en la posición de la distribución.

ggplot(data = datos, aes(x = nacionalidad, y = nota, colour = nacionalidad)) +
    geom_boxplot(outlier.shape = NA) +
    geom_jitter(width = 0.25) +
    scale_color_manual(values = c("gray50", "orangered2")) +
    coord_flip() +
    theme_bw() +
    theme(legend.position = "none")

ggplot(data = datos, aes(x = nota, fill = nacionalidad)) +
    geom_density(alpha = 0.5) +
    scale_fill_manual(values = c("gray50", "orangered2")) +
    theme_bw() +
    theme(legend.position = "bottom")

par(mfrow = c(1, 2))
qqnorm(filter(datos, nacionalidad == "mexicana") %>% pull(nota),
       main = "mexicanos", col = "gray50", pch = 19)
qqline(filter(datos, nacionalidad == "mexicana") %>% pull(nota))

qqnorm(filter(datos, nacionalidad == "otras") %>% pull(nota),
       main = "otras", col = "orangered2", pch = 19)
qqline(filter(datos, nacionalidad == "otras") %>% pull(nota))

par(mfrow =  c(1,1))

Los datos no se distribuyen de forma normal, siendo más marcada la falta de normalidad en las observaciones de nacionalidad no mexicana, que es la muestra de menor tamaño.

Distribución bootstrapping

El objetivo del estudio es determinar si la diferencia observada de 5.9 unidades está dentro de lo que cabría esperar solo por azar debido al muestreo aleatorio, si no existiera diferencia entre las poblaciones (\(H_0\): todas las observaciones proceden de la misma población/distribución). Si se obtuviera una nueva muestra de 2000 estudiantes, la diferencia promedio entre estudiantes mexicanos y no mexicanos sería ligeramente distinta aunque no existiera una diferencia real entre nacionalidades. Por lo tanto, dar respuesta a este problema pasa por determinar cuanta diferencia se espera por el simple hecho de repetir el muestreo.

Mediante bootstrapping se generan nuevas pseudomuestras de individuos empleando las dos muestras originales (mexicanos y no mexicanos) combinadas. Con cada una de las pseudomuestras, se generan dos grupos de tamaños iguales a los grupos de las muestras originales (34, 116) y se calcula la diferencia del estadístico, en este caso la media. El proceso se asemeja mucho al empleado en los test de permutación. La diferencia radica en que en los test de permutación, se emplean siempre las observaciones de la muestra original ordenadas de forma distinta en cada iteración. En este caso, en cada iteración y antes del reparto en grupos, se genera una nueva pseudomuestra que tiene el mismo tamaño que la muestra original pero formada por distintas observaciones.

boot_distribution <- rep(NA, 9999)
n_mexicana <- datos %>% filter(nacionalidad == "mexicana") %>% nrow()
n_otras <- datos %>% filter(nacionalidad == "otras") %>% nrow()

for (i in 1:9999) {
  # Se crea un paseudomuestra
  pseudomuestra <- sample(datos$nota, size = length(datos$nota), replace = TRUE)
  # Se reparten las observaciones en dos grupos y se calculan las medias
  mean_mexicana <- mean(pseudomuestra[1:n_mexicana])
  mean_otras    <- mean(pseudomuestra[(n_mexicana + 1):length(pseudomuestra)])
  # Se calcula la diferencia de medias
  boot_distribution[i] <- mean_mexicana - mean_otras
}

Los datos simulados forman lo que se conoce como bootstrapping distribution y representa la variación esperada en la diferencia de medias debida únicamente al muestreo aleatorio.

ggplot(data = data.frame(boot_distribution), aes(boot_distribution)) +
  geom_histogram(aes(y = ..density..), alpha = 0.4, colour = "white") + 
  geom_line(aes(y = ..density..), stat = 'density', colour = "red") +
  geom_vline(xintercept = mean(boot_distribution)) +
  geom_vline(xintercept = dif_obs, colour = "blue") +
  geom_vline(xintercept = -dif_obs, colour = "blue") +
  labs(title = "Distribución de Bootstrapping",
       x = "diferencia de medias en la pseudomuestra") +
  theme_bw()

summary(boot_distribution)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -12.48606  -2.02264  -0.01684   0.01073   1.98433  11.80502

sd(boot_distribution)

## [1] 3.001633

Como era de esperar, la diferencia entre la media de los grupos considerando que todos proceden de una misma población está centrada en cero. La desviación típica de la bootsrapping distribution indica que la variabilidad en la diferencia de medias debida al muestreo, en este caso es de 3.002 unidades.

Cálculo de p-value

Finalmente, se calcula la probabilidad (p-value) de obtener diferencias igual o más extremas que la observada (líneas verticales azules) con y sin corrección de continuidad.

# Sin corrección
p_value = (sum(abs(boot_distribution) > abs(dif_obs)))/9999
p_value

## [1] 0.04740474

# Con corrección
p_value_corregido = ((sum(abs(boot_distribution) > abs(dif_obs))) + 1)/(9999 + 1)
p_value_corregido

## [1] 0.0475

Conclusión

Un test de bootstrapping no paramétrico se empleó para determinar si existe diferencia en el promedio de logros académicos entre inmigrantes mexicanos y no mexicanos. Los datos se obtuvieron mediante muestreo aleatorio asegurando así la independencia de los datos. El p-value se calculó mediante 9999 simulaciones de Monte Carlo, empleando la corrección de continuidad sugerida por Davison y Hinkley (1997) resultando en 0.048. Esta es una evidencia moderada en contra de la hipótesis nula de que no hay diferencia entre poblaciones y sugiere que los logros académicos entre inmigrantes de nacionalidad mexicana y no mexicana son distintos.

Boostraping para intervalos de confianza de la diferencia entre poblaciones

Los test de hipótesis que contrastan la diferencia entre dos poblaciones generan un p-value basado en la hipótesis nula de que no existe diferencia del estadístico (media, mediana…) entre las poblaciones. Por esta razón, cuando el bootstrapping se emplea como test de hipótesis, las observaciones se mezclan dotas juntas y, posteriormente, se generan las pseudomuestras. Cuando la finalidad del estudio es generar intervalos de confianza para la verdadera diferencia de un estadístico entre dos poblaciones, la hipótesis nula considera que las muestras sí proceden de dos poblaciones diferentes, cada una con un valor distinto para el estadístico estudiado. Para simular esta hipótesis nula, las pseudomuestras obtenidas por sampling with replacement se tienen que generar de forma separada para cada grupo. Esta es la diferencia entre el empleo de bootstrapping para calcular p-values y para calcular intervalos de confianza.

El algoritmo del bootstrapping para calcular intervalos de confianza de la diferencia entre dos poblaciones es:

---

1. Generar una nueva pseudomuestra del grupo A del mismo tamaño que la muestra original \(n_A\) y empleando únicamente las observaciones pertenecientes a dicho grupo.

2. Generar una nueva pseudomuestra del grupo B del mismo tamaño que la muestra original \(n_B\) y empleando únicamente las observaciones pertenecientes a dicho grupo.

3. Calcular la diferencia del estadístico entre las dos nuevas pseudomuestras.

4. Repetir los pasos 1,2 y 3 múltiples veces, almacenando la diferencia calculada en cada iteración. El conjunto de valores generado forma lo que se conoce como la bootstap distribution de la diferencia del estadístico si ambas muestras proceden de dos poblaciones distintas. Esta distribución tiende a estar centrada en el verdadero valor de la diferencia entre ambas poblaciones.

5. A partir de la bootstap distribution, generar un intervalo de confianza para el parámetro poblacional. Las mismas consideraciones explicadas para intervalos de confianza en una única población se aplican también a la diferencia entre poblaciones.

---

Tamaño del efecto (effect size)

El cálculo de p-value y los intervalos de confianza permiten responder a dos de las preguntas fundamentales cuando se comparan poblaciones:

• ¿Existen evidencias que apunten a la existencia de una diferencia real entre las dos poblaciones? o ¿Es la diferencia observada la esperada por simple variabilidad?

• ¿Si el la diferencia es real, cómo de grande es esta diferencia?

A pesar de la importancia de estas dos preguntas, existe una tercera que es fundamental a la hora de tomar decisiones en base a los resultados. ¿Es la diferencia suficientemente grande como para ser práctica? Cuando las unidades en la que se expresa la diferencia son familiares para el analista, se pueden interpretar fácilmente. Sin embargo, en muchas ocasiones, la escala en la que se mide la diferencia no es intuitiva, lo que complica su interpretación así como su divulgación. Por esta razón siempre se debe reportar la diferencia obtenida en unidades estandarizadas, lo que se conoce como tamaño de efecto o effect size. Uno de los tamaños de efecto más empleados es la \(Cohen's \ d\), que consiste en expresar la diferencia entre las medias en términos de desviación típica (pooled sd).

Ejemplo 1

Empleando los datos del estudio de logros académicos entre inmigrantes de nacionalidad mexicana y no mexicana, generar un intervalo de confianza de la diferencia en la media de rendimiento académico entre ambas poblaciones.

datos <- read_csv("./datos/LatinoEd.csv")
datos <- as.data.frame(datos)
datos <- datos %>% select(Achieve, Mex)
datos <- datos %>% rename(nota = Achieve, nacionalidad = Mex)
datos <- datos %>% mutate(nacionalidad = recode(nacionalidad,
                                                `1` = "mexicana",
                                                `0` = "otras"))
head(datos)

datos %>% group_by(nacionalidad) %>% summarise(media = mean(nota),
                                               mediana = median(nota),
                                               sd = sd(nota),
                                               n = n())

mean_mexicana <- datos %>% filter(nacionalidad == "mexicana") %>%
                           pull(nota) %>%
                           mean()
mean_otras <- datos %>% filter(nacionalidad == "otras") %>%
                           pull(nota) %>%
                           mean()
dif_obs <- mean_mexicana - mean_otras
dif_obs

## [1] -5.921602

Se emplean las dos funciones de R que automatizan el proceso de resampling teniendo en cuenta cada grupo:

Función boot()

Esta función requiere que se defina previamente la función que calcula el estadístico de interés, en este caso la media.

library(boot)

#Función que calcula la diferencia entre medias de los dos grupos
diferencia_medias <- function(data, i) {
  pseudomuestra <- data[i,]
  mean_mexicana <- mean(pseudomuestra[pseudomuestra$nacionalidad == "mexicana", "nota"])
  mean_otras <- mean(pseudomuestra[pseudomuestra$nacionalidad == "otras", "nota"])
  return(mean_mexicana - mean_otras)
}


#El argumento strata de la funcion boot() permite identificar los grupos
set.seed(5290)
boot_distribution <- boot(data = datos, statistic = diferencia_medias, R = 9999,
                          strata = as.factor(datos$nacionalidad))
plot(boot_distribution)

boot_distribution

## 
## STRATIFIED BOOTSTRAP
## 
## 
## Call:
## boot(data = datos, statistic = diferencia_medias, R = 9999, strata = as.factor(datos$nacionalidad))
## 
## 
## Bootstrap Statistics :
##      original     bias    std. error
## t1* -5.921602 0.04400832    2.640973

boot.ci(boot_distribution)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 9999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_distribution)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   (-11.142,  -0.789 )   (-11.222,  -0.830 )  
## 
## Level     Percentile            BCa          
## 95%   (-11.013,  -0.621 )   (-11.045,  -0.642 )  
## Calculations and Intervals on Original Scale

Funcion bootES()

La función bootES() hace una llamada interna a boot() pero facilita la sintaxis al no necesitar definir la función del estadístico cando ya existen en R (mean, median…).

library(bootES)
set.seed(5290)
bootES(data = datos, data.col = "nota", group.col = "nacionalidad",
       contrast = c("mexicana", "otras"), R = 9999, plot = TRUE)

## 
## User-specified lambdas: (mexicana, otras)
## Scaled lambdas: (-1, 1)
## 95.00% bca Confidence Interval, 9999 replicates
## Stat        CI (Low)    CI (High)   bias        SE          
## 5.922       0.642       11.045      -0.044      2.641

La distribución generada está centrada en el valor de la diferencia observada (-5.9), ya que el resampling se está haciendo acorde a la hipótesis nula de que las muestras proceden de poblaciones distintas. Además, como la distribución obtenida cumple la normalidad, todos los diferentes tipos de intervalos (percentiles, t, BCa) se consideran válidos.

Tamaño del efecto

La función bootES() permite, además de generar los intervalos de confianza para la diferencia de dos medias en las unidades en las que se está midiendo la variable, calcular intervalos para múltiples tipos de tamaño de efecto. Se puede encontrar una descripción detallada de todos ellos en BootES: An R package for bootstrap confidence intervals on effect sizes (Kris N.Kirby & Daniel Gerlanc).

library(bootES)
bootES(data = datos, data.col = "nota", group.col = "nacionalidad",
       contrast = c("mexicana", "otras"), R = 9999,
       effect.type = "cohens.d")

## 
## User-specified lambdas: (mexicana, otras)
## Scaled lambdas: (-1, 1)
## 95.00% bca Confidence Interval, 9999 replicates
## Stat        CI (Low)    CI (High)   bias        SE          
## 0.392       0.023       0.727       0.003       0.179

Acorde al resultado devuelto por bootES(), el tamaño de efecto es de 0.392 y su CI del 95% (0.03, 0.722). Este valor sugiere que la diferencia observada en las notas aunque significativa, es pequeña.

Ejemplo 2

Se dispone de un set de datos con información sobre el tiempo que una empresa de telefonía tarda en reparar los problemas que tienen dos grupos de consumidores (CLEC y ILEC). Se desea conocer cuál es la diferencia en la media los tiempos de reparación de ambos grupos empleando un intervalo de confianza del 95%.

library("resample")
data("Verizon")
head(Verizon)

Exploración de los datos

Verizon %>% group_by(Group) %>% summarise(media = mean(Time),
                                          median = median(Time),
                                          sd = sd(Time),
                                          n = n())

ggplot(data = Verizon, aes(x = Group, y = Time, colour = Group)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.25) +
  scale_color_manual(values = c("gray50", "orangered2")) +
  coord_flip() +
  theme_bw() +
  theme(legend.position = "none")

ggplot(data = Verizon, aes(x = Time, fill = Group)) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c("gray50", "orangered2")) +
  theme_bw() +
  theme(legend.position = "bottom")

qqnorm(Verizon[Verizon$Group == "CLEC", "Time"], col = "gray50", pch = 19,
       main = "CLEC" )
qqline(Verizon[Verizon$Group == "CLEC", "Time"])

qqnorm(Verizon[Verizon$Group == "ILEC", "Time"], col = "orangered2", pch = 19,
       main = "ILEC" )
qqline(Verizon[Verizon$Group == "ILEC", "Time"])

La representación gráfica muestra una marcada asimetría de los datos. Además, el tamaño de uno de los grupos de solo 23 observaciones.

Los métodos paramétricos basados en el teorema del límite central, tales como los t-test, emplean como Error Estándar la ecuación \(SE = \sqrt{\frac{sd_1^2}{n_1} + \frac{sd_2^2}{n_2}}\). Cuando las condiciones de normalidad no se cumplen, está aproximación no es válida. El bootstrapping permite obtener el valor de SE sin tener que recurrir al teorema del límite central, ya que, la desviación estándar de la bootstrapping distribution se aproxima al SE de la sampling distribution.

Cálculo de la bootstrapping distribution

Se obtienen nuevas pseudomuestras por separado de cada una de las muestras iniciales manteniendo el tamaño original y se calcula la diferencia.

boot_distribution <- rep(NA, 9999)
ILEC <- Verizon %>% filter(Group == "ILEC") %>% pull(Time)
CLEC <- Verizon %>% filter(Group == "CLEC") %>% pull(Time)

for (i in 1:9999) {
    boot_ILEC <- sample(x = ILEC, size = length(ILEC), replace = TRUE)
    boot_CLEC <- sample(x = ILEC, size = length(ILEC), replace = TRUE)
    boot_distribution[i] <- mean(boot_ILEC) - mean(boot_CLEC)
}

Con la función bootES() se puede realizar el mismo proceso de forma más eficiente.

library(bootES)
set.seed(5290)
boot_distribution <- bootES(data = Verizon, data.col = "Time", group.col = "Group",
                            contrast = c("CLEC", "ILEC"), R = 9999)

Estudio la bootstrapping distribution

mean(boot_distribution$t)

## [1] -8.062663

sd(boot_distribution$t)

## [1] 3.916538

plot(boot_distribution)

Dado que la distribución obtenida por bootstrapping no es normal, los intervalos basados en t-distribution o percentiles no son adecuados, solo el de tipo BCa.

Intervalo de confianza

# Intervalo de confianza obtenido por bca
boot_distribution

## 
## User-specified lambdas: (CLEC, ILEC)
## Scaled lambdas: (-1, 1)
## 95.00% bca Confidence Interval, 9999 replicates
## Stat        CI (Low)     CI (High)   bias        SE          
## -8.098      -17.607      -2.134      0.035       3.917

Tamaño del efecto

set.seed(5290)
bootES(data = Verizon, data.col = "Time", group.col = "Group",
       contrast = c("CLEC", "ILEC"), effect.type = "cohens.d",
       R = 9999)

## 
## User-specified lambdas: (CLEC, ILEC)
## Scaled lambdas: (-1, 1)
## 95.00% bca Confidence Interval, 9999 replicates
## Stat        CI (Low)    CI (High)   bias        SE          
## -0.549      -1.151      -0.115      -0.001      0.268

Conclusión

Sí existen evidencias significativas para afirmar que el tiempo promedio de reparación entre los dos grupos es distinto. El intervalo de confianza del 95% para la diferencia de medias obtenido por bootstrapping indica que, en promedio, las reparaciones del grupo CLEC tardan de 17.923 a 1.996 horas más.

Bootstrapping para datos pareados

Un determinado estudio pretende determinar si existe relación entre la luna llena y la agresividad de las personas dementes. Para ello, se mide con un test de comportamiento el nivel de agresividad de 15 pacientes los días de luna llena y los días sin luna llena. Obtener el intervalo de confianza del 95% para la verdadera diferencia entre las dos condiciones.

Dado que se trata de un diseño muestral, el método adecuado es el bootstrapping.

datos <- data.frame(paciente = c(1:15),
                    luna_llena = c(3.33, 3.67, 2.67, 3.33, 3.33, 3.67, 4.67, 2.67,
                                   6, 4.33, 3.33, 0.67, 1.33, 0.33, 2),
                    sin_luna = c(0.27, 0.59, 0.32, 0.19, 1.26, 0.11, 0.30, 0.4,
                                 1.59, 0.6, 0.65, 0.69, 1.26, 0.23,0.38))
head(datos)

mean(datos$luna_llena)

## [1] 3.022

mean(datos$sin_luna)

## [1] 0.5893333

dif_obs <- mean(datos$luna_llena) - mean(datos$sin_luna)
dif_obs

## [1] 2.432667

Cuando se quiere emplear bootstrapping para estudiar la diferencia en el promedio de una variable bajo dos condiciones y los datos son pareados, al igual que en un t-test, se calcula una nueva variable a partir de la diferencia entre las dos condiciones para cada sujeto.

datos$diff <- datos$luna_llena - datos$sin_luna
head(datos)

Una vez se dispone de esta nueva variable, el proceso es igual al seguido cuando se quiere realizar inferencia sobre el parámetro de una población (Boostraping para intervalos de confianza de una población).

Obtención de la distribución de bootstrapping

# Utilizando la función boot()
library(boot)
media <- function(valores,i){
    mean(valores[i])
}

set.seed(5290)
boot_distribution <- boot(datos$diff, media, R = 9999)
boot_distribution

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = datos$diff, statistic = media, R = 9999)
## 
## 
## Bootstrap Statistics :
##     original       bias    std. error
## t1* 2.432667 0.0009820982   0.3650497

plot(boot_distribution)

boot.ci(boot_distribution)

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 9999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_distribution)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 1.716,  3.147 )   ( 1.730,  3.175 )  
## 
## Level     Percentile            BCa          
## 95%   ( 1.690,  3.135 )   ( 1.660,  3.105 )  
## Calculations and Intervals on Original Scale

# Utilizando la función bootES()
set.seed(5290)
bootES(datos$diff, R = 9999, plot = TRUE)

## 
## 95.00% bca Confidence Interval, 9999 replicates
## Stat        CI (Low)    CI (High)   bias        SE          
## 2.433       1.660       3.105       0.001       0.365

Dado que el intervalo de confianza no contiene el valor 0, sí existe evidencia para un alpha de 0.05 de que el promedio de agresividad no es igual con y sin luna llena.

Tamaño del efecto

set.seed(5290)
bootES(datos$diff, R = 9999, effect.type = "cohens.d")

## 
## 95.00% bca Confidence Interval, 9999 replicates
## Stat        CI (Low)    CI (High)   bias        SE          
## 1.666       0.953       2.696       0.160       0.574

Bootstrapping para coeficientes de correlación

Cuando se estudia la correlación entre dos variables continuas, junto con el valor R de la correlación, hay que calcular su intervalo de confianza o su significancia. Por muy alto que sea un coeficiente de correlación, si no es significativo, se tiene que considerar como inexistente puesto que podría deberse simplemente al azar.

Supongase que se quiere saber si existe una relación entre el salario de los jugadores y su calidad (éxito al batear la pelota). Se dispone de una muestra de 50 jugadores seleccionados al azar de la liga. ¿Se puede considerar real la relación, en caso de que exista?

jugador = c("Matt_Williams", "Jim_Thome", "Jim_Edmonds", "Fred_McGriff",
            "Jermaine_Dye", "Edgar_Martinez", "Jeff_Cirillo", "Rey_Ordonez",
            "Edgardo_AlfonzoMoises_Alou", "Travis_Fryman", "Kevin_Young",
            "M._GrudzielanekTony_Batista", "Fernando_Tatis", "Doug_Glanville",
            "Miguel_Tejada", "Bill_Mueller", "Mark_McLemore", "Vinny_Castilla",
            "Brook_Fordyce", "Torii_Hunter", "Michael_Tucker", "Eric_Chavez",
            "Aaron_Boone", "Greg_Colbrunn", "Dave_Martinez", "Einar_Diaz",
            "Brian_L._HunterDavid_Ortiz", "Luis_Alicea", "Ron_Coomer", 
            "Enrique_Wilson", "Dave_Hansen", "Alfonso_SorianoKeith_Lockhart",
            "Mike_Mordecai", "Julio_Lugo", "Mark_L._JohnsonJason_LaRue",
            "Doug_MientkiewiczJay_Gibbons", "Corey_PattersonFelipe_Lopez",
            "Nick_Johnson", "Thomas_Wilson", "Dave_Roberts", "Pablo_Ozuna",
            "Alexis_Sanchez", "Abraham_Nunez","Corey_PattersonFelipe_Lopez",
            "Nick_Johnson", "Thomas_Wilson", "Dave_Roberts", "Pablo_Ozuna",
            "Alexis_Sanchez", "Abraham_Nunez")
salario = c(9500000, 8000000, 7333333, 7250000, 7166667, 7086668, 6375000, 6250000,
            6200000, 6000000, 5825000, 5625000, 5000000, 4900000, 4500000, 4000000,
            3625000, 3450000, 3150000, 3000000, 2500000, 2400000, 2250000, 2125000,
            2100000, 1800000, 1500000, 1087500, 1000000, 950000, 800000, 750000,
            720000, 675000, 630000, 600000, 500000, 325000, 320000, 305000, 285000,
            232500, 227500, 221000, 220650, 220000, 217500, 202000, 202000, 200000)
bateo = c(0.269, 0.282, 0.327, 0.259, 0.240, 0.270, 0.253, 0.238, 0.300, 0.247, 0.213,
          0.238, 0.245, 0.276, 0.268, 0.221, 0.301, 0.242, 0.273, 0.250, 0.208, 0.306,
          0.235, 0.277, 0.227, 0.307, 0.276, 0.216, 0.289, 0.237, 0.202, 0.344, 0.185,
          0.234, 0.324, 0.200, 0.214, 0.262, 0.207, 0.233, 0.259, 0.250, 0.278, 0.237,
          0.235, 0.243, 0.297, 0.333, 0.301, 0.224)
datos <- data.frame(jugador, salario, bateo)
head(datos)

Estudio de la correlación

ggplot(data = datos, mapping = aes(x = bateo, y = salario)) +
    geom_point(colour = "darkred") +
    geom_smooth(method = "lm", colour = "black") +
    theme_bw()

cor(x = datos$bateo, y = datos$salario, method = "pearson")

## [1] 0.1067575

El cálculo de la correlación lineal de Pearson muestran un r = 0.107. ¿Es significativa está correlación?

Distribución bootstrappig del coeficiente de correlación

# Por defecto, si a la función bootES() se le pasan dos columnas numéricas, calcula
# el intervalo de confianza para el coeficiente de correlación de Pearson.
bootES(data = datos[c("bateo", "salario")], R = 9999, ci.type = "bca", plot = TRUE)

## 
## 95.00% bca Confidence Interval, 9999 replicates
## Stat        CI (Low)    CI (High)   bias        SE          
## 0.107       -0.153      0.353       0.002       0.130

Dado que el intervalo de confianza del 95% contiene el valor 0, no se puede considerar significativa la correlación para un \(\alpha\) de 0.05.

Bootstrapping para variables cualitativas (proporciones)

A la hora de comparar variables cualitativas (proporciones), si se cumplen una serie de condiciones, se pueden emplear métodos basados en el teorema del límite central.

• Independencia: El muestreo es aleatorio y el tamaño de la muestra es menor que el 10% de la población.

• Tamaño y distribución: Debe haber al menos 10 eventos verdaderos y 10 eventos falsos esperados dada la hipótesis nula. Esto implica que en el cálculo hay que emplear el null value de p establecido en \(H_0\). \(np_0 >= 10\) , \(n(1-p_0) >= 10\).

Si estas condiciones no se satisfacen, la aproximación mediante el teorema del límite central no es válida. En estos casos es posible calcular el p-value mediante bootstrapping.

Ejemplo

El porcentaje de complicaciones que ocurren en una prueba clínica con el instrumental que tiene un hospital es del 10%. Una compañía considera que su instrumental consigue minimizar el riesgo de complicaciones argumentando que solo 3 intervenciones de las 62 realizadas han tenido complicaciones. ¿Hay evidencias de que sea cierto que su instrumental es mejor, utilizando un límite de significancia del 5%?

Hipótesis

Las hipótesis que se pretenden contrastar son:

• \(H_0\): No existe una asociación entre el nuevo instrumental y el indice de complicaciones. El valor p del nuevo instrumental es igual al del instrumento actual, p = 0.1.

• \(H_a\): Con el nuevo instrumental se consigue reducir el porcentaje de complicaciones p < 0.1.
  
  

Estadístico

Se calcula el valor de la proporción en la muestra \(\hat{p} = 3/62 = 0.0483871\)

Condiciones para solucionarlo mediante el teorema del límite central

• Independencia: El muestreo es aleatorio y el tamaño de la muestra es menor que el 10% de la población.

• Tamaño y distribución: Debe haber en la muestra al menos 10 eventos verdaderos y 10 eventos falsos dada la hipótesis nula.

  • \(np_0\) = 62 x 0.1 = 6.2 no se cumple la condición

  • \(n(1-p_0)\) = 62 x 0.9 = 55.8

Al no cumplirse las condiciones, no se puede considerar que \(\hat{p}\) se distribuye de forma normal por lo que no se pueden aplicar test paramétricos. Para poder saber cual es la probabilidad de que siendo \(H_0\) cierta (p = 0.1) de 62 eventos solo 3 o menos sean positivos, es decir, de que \(\hat{p}\) = 0.0483871, se pueden recurrir a bootstrapping.

Distribución de bootstrapping

boot_distribution <- rep(NA, 9999)
# La muestra tiene un 10% de true. Esto equivale a la Ho
muestra_original <- c(rep(TRUE, 10), rep(FALSE, 90)) 

for (i in 1:9999) {
    boot_muestra <- sample(muestra_original, 62, replace = TRUE)
    # Hacer la media de un vector lógico es como calcular la proporción de True
    boot_distribution[i] <- mean(boot_muestra)
}

hist(boot_distribution, breaks = 30,
     main = "distribución bootstrapping",
     xlab = "proporción de complicaciones")
abline(v = 0.0483871, col = "red")

Cálculo de p-value

El p-value es igual a la proporción de simulaciones en las que el valor \(\hat{p}\) es menor que el \(p\) observado.

sum(boot_distribution <= 0.0483871) / 9999

## [1] 0.1227123

Conclusión

Dado que p-value es mayor que \(\alpha = 0.05\), no hay evidencias suficientes para considerar que los resultados que muestra la compañía no se puedan obtener por fluctuaciones aleatorias siendo la verdadera probabilidad de complicaciones del 10%.

Test de permutaciones para variables cualitativas (proporciones)

Un estudio quiere determinar si un fármaco reduce el riesgo de muerte tras una intervención del corazón. Se diseña un experimento en el que pacientes que tienen que ser operados se distribuyen aleatoriamente en dos grupos (control y tratamiento). En vista de los resultados, ¿Se puede decir que el fármaco es efectivo para un nivel de significancia del 5%?.

grupo	vivo	muerto	total
control	11	39	50
tratamiento	14	26	40
———–	—-	——	—–
total	25	65	90

Hipótesis

• \(H_0\): El porcentaje de supervivencia es independiente del tratamiento, la proporción de vivos es la misma en ambos grupos, p(control) - p(tratamiento) = 0.

• \(H_a\): El porcentaje de supervivencia es dependiente del tratamiento, p(control) - p(tratamiento) != 0.
  
  

Estadístico

\(\hat{p}(control)\) - \(\hat{p}(tratamiento)\) = 11/50 - 14/40 = -0.13

Cálculo de \(p_{pool}\)

Siendo \(p_{pool} = \frac{total \ positivos}{total \ eventos} = \frac{positivos1 + positivos2}{n1 + n2}\)

\(p_{pool}\) = (11 + 14)/(50+40) = 0.278

Condiciones para solucionarlo mediante el teorema del límite central

Independencia:

• Dentro de cada muestra: El muestreo debe ser aleatorio y el tamaño de la muestra menor que el 10% de la población.

• Entre grupos/muestras: los datos no pueden ser parados.

Tamaño y distribución:

• \(n_1p_{pool}\) & \(n_1(1-p_{pool})\) = 13.9, 36.1 >= 10 Al menos 10 eventos positivos y negativos esperados en la muestra 1.

• \(n_2p_{pool}\) & \(n_(1-p_{pool})\)= 11.12, 28.88 >= 10 Al menos 10 eventos positivos y negativos esperados en la muestra 2.

A pesar de que las condiciones se cumplen, están muy cerca del límite, por lo que la aproximación por teorema del límite central puede no ser del todo precisa.

*Cálculo de p-value mediante test de permutaciones**

Mediante permutaciones se puede estimar la distribución muestral de p(control) - p(tratamiento). Acorde a la \(H_0\), la probabilidad de supervivencia en la misma en ambos grupos, para simular esto, se redistribuyen aleatoriamente (manteniendo el tamaño de cada grupo) los sujetos y se calcula la diferencia de proporciones. De los 90 individuos hay 50 control y 40 tratados.

total_individuos <- c(rep(TRUE, 25), rep(FALSE, 65))
posicion_individuos <- c(1:90)
permutation_distribution <- rep(NA,9999)

for (i in 1:9999) {
    a <- sample(posicion_individuos, 50)
    control <- total_individuos[a]
    # El resto de individuos va al grupo tratamiento
    b <- posicion_individuos[!(posicion_individuos %in% a)]
    tratamiento <- total_individuos[b]
    
    # La media de un vector lógico es la proporcion de TRUEs
    permutation_distribution[i] <- mean(control) - mean(tratamiento)
}

hist(permutation_distribution,breaks = 50, freq = FALSE ,
     main = "distribución permutaciones",
     xlab = "diferencia de proporcion entre grupos")
abline(v = -0.13, col = "red")
abline(v = 0.13, col = "red")

El p-value es igual al número de permutaciones en las que el valor de p calculado es menor igual a -0.13 o mayor igual a 0.13, dividido por el número total permutaciones.

mean(permutation_distribution >=0.13 | permutation_distribution <= -0.13)

## [1] 0.1671167

Conclusión

Dado que p-value es mayor que el nivel de significancia \((\alpha = 0.05)\), no hay evidencias suficientes para considerar que los resultados muestren una relación entre el tratamiento y la proporción de supervivientes.

Información sesión

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)

Bibliografía

Comparing groups Randomization and Bootstrap Methods using R. Andrew S.Zieffler

Bootstrap Methods and Permutation Test by Tim Hestemberg

Introduction to Statistical Learning, Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani

The RBook Michael J Crawley

¿Cómo citar este documento?

Resampling: test de permutación, simulación de Monte Carlo y Bootstrapping por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/23_resampling_test_permutacion_simulacion_de_monte_carlo_bootstrapping

¿Te ha gustado el artículo? Tu ayuda es importante

Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊

Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.