Analítica Avanzada de Datos

Sesión #5: Introducción al Bootstrap

Autor/a
Afiliación

Jorge I. Vélez, PhD

Universidad del Norte, Barranquilla

Fecha de publicación

7 de junio de 2024

Introducción

El bootstrap, también conocido como métodos inferencial basado en remuestreo, fue propuesto en 1979 por Bradley Efron de la Universidad de Stanford. A la fecha, el artículo original tiene más de 25900 citaciones.

The concept is simple: To pull yourself up by your bootstraps means to succeed or elevate yourself without any outside help.

Filosofía del bootstrap. Imagen tomada de aquí.

La gran ventaja del bootstrap es su simpleza. Sin embargo, la principal desventaja del proceso inferencial que se realiza a través de bootstrap es que depende de la calidad de la muestra.

En esta práctica de R trabajaremos la sintaxis para realiza bootstrap a través de diferentes ejemplos.

Ejemplos

Inferencia para \(\widetilde{\mu}\)

En algunas aplicaciones interesa calcular la mediana, además de su desviación estándar y, posiblemente, realizar inferencia a partir del cálculo de intervalos de confianza.

Aunque el cálculo de la mediana puede realizarse en R con facilidad a través de la función median(), el cálculo de su desviación estándar y el intervalo de confianza es más complejo matemáticamente. Una forma sencilla de aproximar estas cantidades es vía bootstrap.

Consideremos un vector de datos dado por:

## datos
x <- c(7, 6, 12, 11, 17, 15, 8, 10)

La mediana muestral es

## mediana 
(mediana <- median(x))
[1] 10.5

Ahora, usemos bootstrap para calcular un intervalo de confianza del 95%. Los pasos a seguir son:

  1. Obtener una muestra con reemplazo del vector x. Para ello usamos la función sample():
## muestra con reemplazo
set.seed(721)
(x.samp <- sample(x, replace = TRUE))
[1]  7  8 11  7  8 11 11 11
  1. Para la muestra anterior, calculamos la mediana muestral:
## mediana para la muestra aleatoria con reemplazo
(median.boot <- median(x.samp))
[1] 9.5
  1. Repetir los pasos 1 y 2, \(B\) veces, para obtener \(B\) valores de la mediana.
## semilla
set.seed(721)

## número de muestras
B <- 1000

## réplicas
medianas.boot <- replicate(B, {
  x.samp <- sample(x, replace = TRUE)
  median(x.samp)
})
  1. Construir un histograma con los valores en 3, y calcular los percentiles 2.5% y 97.5%.
## tomamos el objeto medianas.boot
## distribución bootstrap de la mediana
require(ggplot2)
d <- data.frame(mediana = medianas.boot)
ggplot(d, aes(x = mediana)) + 
  geom_histogram(fill = 'purple', color = 'purple') +
  ylab("Frecuencia") +
  xlab("Mediana") +
  theme_minimal()

Finalmente, el intervalo de confianza del 95% para la mediana es

## intervalo de confianza
quantile(medianas.boot, probs = c(0.025, 0.975))
 2.5% 97.5% 
    7    15

Por lo tanto, a nivel poblacional,

\[ \tilde{\mu} \in (7, 15) \] con una confianza del 95%.

Inferencia para \(\rho\)

Trabajaremos con los siguientes datos:

## lectura de datos
d <- read.table('https://www.dropbox.com/s/b0ijd120l4pjxb8/llantas.txt?dl=1',
                header = TRUE)
d$x1 <- d$x3 <- NULL
colnames(d)[2] <- 'x'
head(d)
         y  x
1 23.72364 10
2 30.99403  7
3 17.03163 10
4 26.08888  7
5 30.49450  9
6 23.71911 10

Gráficamente los datos pueden representarse haciendo

## gráfico de dispersión
ggplot(d, aes(x = x, y = y)) +
  geom_point() + theme_minimal()

El coeficiente de correlación entre x e y puede obtenerse haciendo

## coeficiente de correlación muestral
with(d, cor(x, y))
[1] -0.7624853

Para construir un intervalo de confianza del 95% para \(\rho\), primero construimos una función que estime el coeficiente dados dos vectores x e y, y luego tomamos muestras aleatorias.

La función para calcular el coeficiente de correlación lineal es:

## cálculo del coeficiente de correlación 
cor_coef <- function(x, y){
  cor(x, y)
}

Ahora, utilizamos el mismo procedimiento que usamos para la mediana y calculamos el intervalo de confianza del 95% para \(\rho\):

## B valores de rho 
set.seed(123)
B <- 1000
rho.boot <- replicate(B, {
  idx <- sample(NROW(d), replace = TRUE)
  d_samp <- d[idx, ]
  with(d_samp, cor_coef(x, y))
  })

Los primeros 10 valores \(\hat{\rho}^{(1)}, \hat{\rho}^{(2)},\ldots,\hat{\rho}^{(10)}\) son

## primeros 10 valores
rho.boot[1:10]
 [1] -0.7903808 -0.8322427 -0.7924452 -0.7940294 -0.6507087 -0.6831151
 [7] -0.7550041 -0.7075198 -0.7529649 -0.7614872

La distribución muestral de \(\hat{\rho}^{(1)}, \hat{\rho}^{(2)},\ldots,\hat{\rho}^{(1000)}\) es

## distribución muestral de rho
hist(rho.boot, col = 2, border = 2, las = 1, 
     xlab = expression(rho^"*"), main = '', yaxt = 'n', ylab = "")
abline(v = with(d, cor(x,y)), col = 'yellow', lwd = 2)

Finalmente el intervalo de confianza del 95% puede obtenerse como

## intervalo de confianza del 95%
quantile(rho.boot, probs = c(0.025, 0.975))
      2.5%      97.5% 
-0.8273566 -0.6797238

Regresión Logística

Consideremos los siguientes datos de un proceso de inyección:

## lectura de datos
d <- read.delim2("https://www.dropbox.com/scl/fi/jxlp19yfithd2u6dikfn4/logistic3.txt?rlkey=69hp7ujojhb2872fe6a71kzp1&dl=1", header = TRUE)
d$change <- as.numeric(as.character(d$change))  
head(d)
  y     change
1 0 -1.9396807
2 1 -0.7245350
3 0 -2.2534429
4 0  1.3929212
5 0 -0.5057383
6 0 -2.2307026

Por la naturaleza de la variable respuesta \(y\), tiene sentido ajustar un modelo de Regresión Logística.

Los coeficientes del modelo ajustado son:

## modelo de RL
fit <- glm(y ~., data = d, family = 'binomial')
summary(fit)$coefficients
               Estimate Std. Error    z value     Pr(>|z|)
(Intercept) -0.03017878  0.1271316 -0.2373821 8.123604e-01
change       1.12436593  0.1106093 10.1651993 2.835417e-24

Calculemos \(\hat{\theta}|_{\text{change}= 1.5}\) usando R:

## predicción
x0 <- data.frame(change = 1.5)
(theta_hat <- predict(fit, x0, type = 'response'))
        1 
0.8397501

Ahora, usemos Bootstrap para construir la distribución muestral de \(\hat{\theta}|_{\text{change}= 1.5}\) con \(B=500\):

## distribucion muestra de theta | x0
B <- 500

## bootstrap
set.seed(123)
theta_star <- rep(0, B)
for(i in 1:B){
  idx <- sample(NROW(d), replace = TRUE)     
  fit_samp <- glm(y ~., data = d[idx, ], family = 'binomial')
  aux <- predict(fit_samp, x0, type = 'response')
  theta_star[i] <- aux
}

Gráficamente, esta distribución es:

## gráfico
hist(theta_star, col = 'mediumslateblue', 
     border = 'mediumslateblue', las = 1, 
     xlab = expression(theta^"*") , main = '', 
     yaxt = 'n', ylab = "")
abline(v = theta_hat, col = 'yellow', lty = 2, lwd = 2)

Finalmente el intervalo de confianza del 95% para \(\theta|_{\text{change}= 1.5}\) puede obtenerse como

## intervalo de confianza del 95%
quantile(theta_star, probs = c(0.025, 0.975))
     2.5%     97.5% 
0.7615913 0.8995594

Matriz de confusión

Consideremos la siguiente matriz de confusión proveniente de un modelo de Regresión Logística:

## matriz de confusion
confusion <- matrix(c(160, 32,
                       20, 246), ncol = 2, byrow = TRUE)
colnames(confusion) <- rownames(confusion) <- c('Yes', 'No')
confusion
    Yes  No
Yes 160  32
No   20 246

Ahora, carguemos la colección de funciones presentadas con anterioridad:

## load functions
source('https://www.dropbox.com/s/xclvdugfbrf5ryn/logistic-functions.R?dl=1')
This is a collection of functions created by 
Professor Jorge Vélez, PhD <https://jorgeivanvelez.netlify.app/> 
for the fitting and validation of Binary Classification models 
 
In case of problems, please an email to 
<jvelezv@uninorte.edu.co> with enough details 
 
Last modification:  April 23, 2024

Para calcular las medidas de desempeño hacemos:

## medidas de desempeño
measures(confusion)
sensitivity specificity         ppv         npv         fdr         fpr 
  0.8888889   0.8848921   0.8333333   0.9248120   0.1666667   0.1151079 
 class_rate        lift       delta 
  0.8864629   2.1203704   0.1599860

Por otro lado, el AUC es

## AUC
rft(confusion, graph = FALSE)$auc
[1] 0.8868905

Calculemos un intervalo de confianza del 95% para el AUC usando bootstrap paramétrico. Por tratarse de una matriz de confusión, generaremos muestras aleatorias de una distribución multinomial utilizando la función rmultinom() de R.

Por ejemplo, para generar \(n=1\) muestra de una distribución multinomial con parámetros \(\mathbf{\pi} = (\pi_{11}, \pi_{12}, \pi_{21}, \pi_{22})\) y \(\mathbf{n} = (n_{11}, n_{12}, n_{21}, n_{22})\) hacemos:

## draw data from multinomial distribution
set.seed(123)
N <- sum(confusion)
p <- confusion/sum(confusion)
confusion_samp <- rmultinom(n = 1, 
                            size = N, 
                            prob = p)
(confusion_samp <- matrix(confusion_samp, ncol = 2, byrow = FALSE))
     [,1] [,2]
[1,]  154   40
[2,]   19  245

El paso siguiente es calcular el AUC para la matriz de confusión generada:

## AUC
rft(confusion_samp, graph = FALSE)$auc
[1] 0.8749113

Si replicamos estos pasos \(B=10000\) veces, entonces

## draw data from multinomial distribution and estimate AUC
set.seed(123)
B <- 10000
auc_samp <- replicate(B, {
              confusion_samp <- rmultinom(n = 1, size = N, prob = p)
              confusion_samp <- matrix(confusion_samp, ncol = 2, byrow = FALSE)
              rft(confusion_samp,graph = FALSE)$auc
})

La distribución muestral del AUC es

## distribución muestral del AUC
d0 <- data.frame(statistic = auc_samp)
ggplot(d0, aes(x = statistic)) +
  geom_histogram(fill = 'olivedrab4', color = 'olivedrab4') + 
  xlab("AUC basado en Bootstrap") + 
  ylab("Frecuencia") + 
  theme_minimal()

Finalmente, un intervalo de confianza del 95% para el AUC es

## 95% CI
quantile(auc_samp, probs = c(0.025, 0.975))
     2.5%     97.5% 
0.8566661 0.9153940

Estudiemos la distribución de otras medidas de desempeño utilizando Bootstrap:

## create an empty matrix
others <- matrix(0, nrow = B, ncol = 9)
colnames(others) <- c("sensitivity", "specificity", 
                      "ppv", "npv", "fdr", "fpr", "accuracy", 
                      "lift", "delta")

## draw data from multinomial distribution and estimate AUC
set.seed(123)
for(i in 1:B){
  confusion_samp <- rmultinom(n = 1, size = N, prob = p)
  confusion_samp <- matrix(confusion_samp, ncol = 2, byrow = FALSE)
  others[i, ] <- measures(confusion_samp)
}

El objeto others contiene las medidas de desempeño para 10000 matrices de confusión generadas con Bootstrap paramétrico. Los resultados para las primeras 6 muestras Bootstrap son:

## first 6 rows
head(others)
     sensitivity specificity       ppv       npv       fdr        fpr  accuracy
[1,]   0.8901734   0.8596491 0.7938144 0.9280303 0.2061856 0.14035088 0.8711790
[2,]   0.8531073   0.9039146 0.8483146 0.9071429 0.1516854 0.09608541 0.8842795
[3,]   0.9027027   0.8571429 0.8106796 0.9285714 0.1893204 0.14285714 0.8755459
[4,]   0.8793103   0.8697183 0.8052632 0.9216418 0.1947368 0.13028169 0.8733624
[5,]   0.8955224   0.8677043 0.8411215 0.9139344 0.1588785 0.13229572 0.8799127
[6,]   0.8628571   0.8798587 0.8162162 0.9120879 0.1837838 0.12014134 0.8733624
         lift     delta
[1,] 2.101543 0.1782141
[2,] 2.195074 0.1755274
[3,] 2.006980 0.1728436
[4,] 2.119601 0.1775931
[5,] 1.916585 0.1685756
[6,] 2.136154 0.1823242

Así por ejemplo, la distribución conjunta de sensitivity y specificity puede obtenerse a partir del objeto others:

## read msp function
source('https://www.dropbox.com/s/sssfy803jq5nkqu/msp.R?dl=1')

## plot
msp(others[, "sensitivity"], others[, "specificity"], las = 1,
    xlab = "Sensitivity",
     ylab = "Specificity",
    plotlines = FALSE, typel = 'n')

Existe alguna correlación entre el Accuracy y el AUC? Veamos:

## correlation
cor(others[, "accuracy"], auc_samp)
[1] 0.975663

Cómo son las correlaciones de las otras medidas?

## disponibilidad del paquete GGally 
if(!require(GGally)) install.packages("GGally", dependencies = TRUE)
require("GGally")

## gráfico de dispersión/correlación
ggpairs(data.frame(others)) + theme_minimal()

Esta es otra posibilidad.

Homework

Supongamos que el cutoff óptimo para el modelo de Regresión Logística anterior es \(\theta_0 = 0.5826\). Así las cosas, la matriz de confusión es

## matriz de confusión
theta_0 <- 0.5826
theta_hat <- predict(fit, type = 'response')
cm <- table(model = 1*(theta_hat > theta_0), real = d$y)
cm <- cm[2:1, 2:1]
cm
     real
model   1   0
    1  79  20
    0  81 320

De acuerdo con los resultados, la especificidad es 0.941.

Calcule un intervalo de confianza del 95% para el Jaccard index usando Bootstrap.

Respuesta

Introducción al Índice de Jaccard

El Índice de Jaccard, utilizado para comparar la similitud y diversidad de los conjuntos de muestra, se calcula como:

\[ J = \frac{|A \cap B|}{|A \cup B|} \]

En el contexto de una matriz de confusión para un clasificador binario, los componentes son:

  • Verdaderos Positivos (TP) = 79
  • Falsos Positivos (FP) = 81
  • Falsos Negativos (FN) = 20
  • Verdaderos Negativos (TN) = 320

El Índice de Jaccard se calcula utilizando la siguiente fórmula específica para los resultados de un clasificador:

\[ J = \frac{TP}{TP + FP + FN} \]

Este índice mide la precisión de la clasificación excluyendo los verdaderos negativos, enfocándose únicamente en los aspectos positivos identificados correctamente y aquellos identificados incorrectamente ya sea como positivos o negativos. Para aplicar Bootstrap, seguiremos estos pasos:

Re-muestreo: Generar nuevas matrices de confusión re-muestreando los datos originales con reemplazo.

Cálculo del Índice de Jaccard: Para cada matriz de confusión re-muestreada, calcular el índice.

Construcción del Intervalo de Confianza:Utilizar los índices de Jaccard de las muestras re-muestreadas para calcular el percentil 2.5 y 97.5, formando así un intervalo de confianza del 95%.

calculo del Indice de Jaccard para una muestra

Definimos las probabilidades y el tamaño total de la muestra basados en nuestra matriz de confusión, cm.

set.seed(015)  
n_total <- sum(cm)  
probs <- cm / n_total

Generamos nuevas matrices de confusión re-muestreando de acuerdo con las probabilidades definidas.

conf_matrix_sample <- rmultinom(n = 1, size = n_total, prob = probs)
conf_matrix_sample <- matrix(conf_matrix_sample, ncol = 2, byrow = TRUE)

Calculamos el índice para la muestra simulada basándonos solo en las entradas relevantes de la matriz.

jaccard_index <- conf_matrix_sample[1, 1] / (sum(conf_matrix_sample[1,]) + conf_matrix_sample[2, 1])

Repetimos el proceso múltiples veces para obtener una distribución del Índice de Jaccard.

B <- 10000  # Número de réplicas de bootstrap
jaccard_boot <- replicate(B, {
    conf_matrix_sample <- rmultinom(n = 1, size = n_total, prob = probs)
    conf_matrix_sample <- matrix(conf_matrix_sample, ncol = 2, byrow = TRUE)
    conf_matrix_sample[1, 1] / (sum(conf_matrix_sample[1,]) + conf_matrix_sample[2, 1])
})

Visualizamos la distribución de los índices calculados para entender mejor su variabilidad.

df <- data.frame(IndicesJaccard = jaccard_boot)
ggplot(df, aes(x = IndicesJaccard)) +
  geom_histogram(bins = 30, fill = 'steelblue', color = 'black') +
  xlab("Índice de Jaccard") + 
  ylab("Frecuencia") +
  ggtitle("Distribución del Índice de Jaccard Bootstrap") +
  theme_minimal()

El histograma del Índice de Jaccard, obtenido mediante Bootstrap, muestra una distribución simétrica y centrada, con la mayoría de valores entre 0.40 y 0.50, lo que indica estabilidad en las estimaciones y una variabilidad moderada. Esta distribución bien comportada sugiere que los intervalos de confianza del 95% calculados son fiables y proporcionan una medida robusta de la precisión del Índice de Jaccard en la población, destacando la efectividad del Bootstrap para estimar incertidumbres en situaciones donde métodos tradicionales pueden no ser aplicables directamente.

# Intervalo de confianza del 95%
ci <- quantile(jaccard_boot, probs = c(0.025, 0.975))
ci
     2.5%     97.5% 
0.3668300 0.5136612

El intervalo de confianza del 95% para el Índice de Jaccard, basado en 10000 muestras bootstrap, va de aproximadamente 0.387 a 0.497. Este intervalo ofrece una visión útil sobre la precisión del modelo de clasificación en términos de su capacidad para identificar correctamente los casos positivos, ajustados por los errores de clasificación tanto positivos como negativos.

Más recursos