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.93x10105), 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.

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: σ2control=σ2tratamiento

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.

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=sd√n deja de ser buena y, por lo tanto, también los intervalos paramétricos basados en la estructura [parámetro estimado±tα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±tα,dfSE]

• media = mean(datos) = 30.9686559

• tα=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

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 (H0: 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.

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=√sd21n1+sd22n2. 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.

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:

• H0: 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.

• Ha: 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 ˆ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.

  • np0 = 62 x 0.1 = 6.2 no se cumple la condición

  • n(1−p0) = 62 x 0.9 = 55.8

Al no cumplirse las condiciones, no se puede considerar que ˆ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 H0 cierta (p = 0.1) de 62 eventos solo 3 o menos sean positivos, es decir, de que ˆ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 ˆ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 α=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%.