Estas notas se desarrollaron con base en el libro An Introduction to the Bootstrap de Effron y Tibshirani y abordan los siguientes temas:

Muestras aleatorias

Supongamos que tenemos una población finita o universo \(U\), conformado por unidades individuales \(U_1,U_2,...,U_N\) cada una tiene la misma probabilidad de ser seleccionada en una extracción aleatoria. Las unidades individuales \(U_i\) tienen propiedades que nos gustaría aprender (opinión política,…). Debido a que es muy difícil y caro examinar cada unidad en \(U\) seleccionamos una muestra aleatoria.

Una muestra aleatoria de tamaño \(n\) se define como una colección de \(n\) unidades \(u_1,...,u_n\) seleccionadas aleatoriamente de un universo \(U\).

En principio el proceso de muestreo es como sigue:

  1. Seleccionamos \(n\) enteros de manera independiente \(j_1,...,j_n\) (con probabilidad \(1/N\)), cada uno de ellos asociado a un número entre \(1\) y \(N\).

  2. Los enteros determinan las unidades que seleccionamos: \(u_1=U_{j_1},u_2=U_{j_2},...,u_n=U_{j_n}\).

En la práctica el proceso de selección suele ser más complicado y la definición de la población \(U\) suele ser deficiente; sin embargo el marco conceptual sigue siendo útil para entender la inferencia estadística.

Observación: Nuestra definición de muestra aleatoria permite que una unidad particular \(U_i\) aparezca más de una vez, podríamos evitar esto si realizamos un muestreo sin remplazo; sin embargo, es un poco más sencillo permitir repeticiones y si el tamaño de la muestra \(n\) es mucho más chico que la población \(N\), la probabilidad de muestrear la misma unidad más de una vez es chica.

Una vez que se selecciona una muestra aleatoria \(u_1,...,u_n\) obtenemos una o más medidas de interés para cada unidad. Los datos observados son la colección de medidas \(x_1,...,x_n\), que también denotaremos \(\textbf{x} = (x_1,...,x_n)\).

También podemos obtener las medidas de interés de cada unidad en la población \(U_1,U_2,...,U_N\), obteniendo así los valores \(X_1,...,X_N\), esto sería un censo, y denotamos al conjunto de mediciones de la población por \(\mathcal{X}\). El objetivo de la inferencia estadística es expresar lo que hemos aprendido de la población \(\mathcal{X}\) a partir de los datos observados \(\textbf{x}\). En particular, vamos a usar el bootstrap para determinar la precisión con la que una estadística (e.g. media o mediana) calculada de la muestra \(x_1,...,x_n\) estima la cantidad correspondiente en la población.

Veamos un ejemplo artificial donde tenemos una muestra de 500 escuelas primarias de la Ciudad de México, tomada de un universo de nrow(prim) escuelas,

set.seed(16021)
n <- 500
prim_muestra <- sample_n(prim, n, replace = TRUE)
glimpse(prim_muestra)
#> Observations: 500
#> Variables: 6
#> $ clave <fctr> 09DPR3043G, 09PPR1369C, 09DBN0026X, 09DPR2982T, 09DPR29...
#> $ turno <fctr> MATUTINO, MATUTINO, NOCTURNO, VESPERTINO, MATUTINO, MAT...
#> $ tipo  <chr> "GENERAL", "PARTICULAR", "GENERAL", "GENERAL", "GENERAL"...
#> $ mun   <int> 16, 18, 26, 16, 26, 18, 25, 16, 13, 16, 26, 20, 16, 13, ...
#> $ esp3  <dbl> 539, 670, 630, 536, 539, 543, 689, 587, 579, 627, 615, 5...
#> $ esp6  <dbl> 543, 592, 0, 563, 588, 618, 711, 567, 557, 545, 608, 558...

para cada escuela en la muestra tenemos la medida \(x_i\), conformada por el promedio de las calificaciones en español de los alumnos de tercero y sexto de primaria (prueba ENLACE 2010): \[x_i=(esp3_i, esp6_i)\]

Este ejemplo es artificial pues contamos con un censo de las escuelas, sin embargo es común contar únicamente con la muestra, esta tiene una media de 571.575, con un error estándar estimado de 4.015. Debido a que nuestro ejemplo es artificial podemos comparar con la población, la media de las 3311 escuela es 574.775.

El principio del plug-in

Recordemos la definición de la distribución empírica.

Dada una muestra aleatoria de tamaño \(n\) de una distribución de probabilidad \(P\), la función de distribución empírica \(P_n\) se define como la distribución que asigna probabilidad \(1/n\) a cada valor \(x_i\) con \(i=1,2,..., n\). En otras palabras, \(P_n\) asigna a un conjunto \(A\) en el espacio muestral de \(x\) la probabilidad empírica:

\[P_n(A)=\#\{x_i \in A \}/n\]

Ahora, muchos problemas de inferencia estadística involucran la estimación de algún aspecto de una distribución de de probabilidad \(P\) en base a una muestra aleatoria obtenida de \(P\). La función de distribución empírica \(P_n\) es una estimación de la distribución completa \(P\), por lo que una manera inmediata de estimar aspectos de \(P\) (e.g media o mediana) es calcular el aspecto correspondiente de \(P_n\).

Podemos comparar el histograma de la distribución completa con el histograma de la distribución empírica para el ejemplo de las calificaciones de la prueba ENLACE.

claves_muestra <- prim_muestra$clave
prim_long <- prim %>% 
    gather(grado, calif, esp3, esp6) %>%
    mutate(muestra = ifelse(clave %in% claves_muestra, "muestra", "población"))

ggplot(prim_long, aes(x = calif)) +
  geom_histogram(aes(y = ..density..), binwidth = 20, fill = "darkgray") +
  facet_grid(grado ~ muestra)

Cuando la variable de interés toma pocos valores es fácil ver la distribución empírica, supongamos que la medición de las unidades que nos interesa es la variable tipo de escuela, entonces la distribución empírica en la muestra es

table(prim_muestra$tipo) / n
#> 
#>    GENERAL PARTICULAR 
#>      0.654      0.346

Vale la pena notar que pasar de la muestra desagregada a la distribución empírica (lista de valores y la proporción que ocurre cada una en la muestra) no conlleva ninguna pérdida de información: el vector de frecuencias observadas es un estadístico suficiente para la verdadera distribución. Esto quiere decir que toda la información de \(P\) contenida en el vector de observaciones \(\textbf{x}\) está también contenida en \(P_n\).

Nota: el teorema de suficiencia asume que las observaciones \(\textbf{x}\) son una muestra aleatoria de la distribución \(P\), este no es siempre el caso (e.g. si tenemos una serie de tiempo).

Cuando aplicamos teoría estadística a problemas reales, es común que las respuestas estén dadas en términos de distribuciones de probabilidad. Por ejemplo, podemos preguntarnos que tan correlacionados están los resultados de las pruebas de español correspondientes a 3^o y 6^o. Si conocemos la distribución de probabilidad \(P\) contestar esta pregunta es simplemente cuestión de aritmética, el coeficiente de correlación poblacional esta dado por:

\[corr(y,z) = \frac{\sum_{j=1}^{N}(Y_j - \mu_y)(Z_j-\mu_z)} {[\sum_{j=1}^{N}(Y_j - \mu_y)^2\sum_{j=1}^{N}(Z_j - \mu_z)^2]^{1/2}}\]

en nuestro ejemplo \((Y_j,Z_j)\) son el j-ésimo punto en la población de escuelas primarias \(\mathcal{X}\), \(\mu_y=\sum Y_j/3311\) y \(\mu_z=\sum Z_j/3311\).

ggplot(prim, aes(x = esp3, y = esp6)) +
  geom_point(alpha = 0.5)

cor(prim$esp3, prim$esp6)
#> [1] 0.24

Si no tenemos un censo debemos inferir, podríamos estimar la correlación \(corr(y,z)\) a través del coeficiente de correlación muestral: \[\hat{corr}(y,z) = \frac{\sum_{j=1}^{n}(y_j - \hat{\mu}_y)(z_j-\hat{\mu}_z)} {[\sum_{j=1}^{n}(y_j - \hat{\mu}_y)^2\sum_{j=1}^{n}(z_j - \hat{\mu}_z)^2]^{1/2}}\]

cor(prim_muestra$esp3, prim_muestra$esp6)
#> [1] 0.28

Otros ejemplos de estimaciones plug-in:

median(prim_muestra$esp3)
#> [1] 562

\[\theta=\frac{1}{N}\sum_{j=1}^N I_{\{Y_i>700\}}\]

donde \(I_{\{\cdot\}}\) es la función indicadora.

Hacemos la estimación plug-in \(\hat{\theta}\):

sum(prim_muestra$esp3 > 700) / n
#> [1] 0.046

Volvamos ahora al ejemplo de los 100 lanzamientos de un dado, en este ejemplo obteníamos la siguiente distribución empírica:

dado <- read.table("../04-Probabilidad/data/dado.csv", header=TRUE, quote="\"")
table(dado$observado) / n
#> 
#>     1     2     3     4     5     6 
#> 0.026 0.038 0.020 0.034 0.028 0.054

En este caso no tenemos un censo, solo contamos con la muestra. Una pregunta de inferencia que surge de manera natural es si el dado es justo, esto es, si la distribución que generó esta muestra tiene una distribución \(P = (1/6, 1/6, 1/6,1/6, 1/6, 1/6)\). Para resolver esta pregunta, debemos hacer inferencia de la distribución empírica.

Antes de proseguir repasemos dos conceptos importantes: parámetros y estadísticos:

Un parámetro es una función de la distribución de probabilidad \(\theta=t(P)\), mientras que un estadístico es una función de la muestra \(\textbf{x}\).

Por ejemplo, la \(corr(x,y)\) es un parámetro de \(P\) y \(\hat{corr}(x,y)\) es un estadístico con base en \(\textbf{x}\) y \(\textbf{y}\).

Entonces:

El principio del plug-in es un método para estimar parámetros a partir de muestras; la estimación plug-in de un parámetro \(\theta=t(P)\) se define como: \[\hat{\theta}=t(P_n).\]

Es decir, estimamos la función \(\theta = t(P)\) de la distribución de probabilidad \(P\) con la misma función de la distribución empírica \(\hat{\theta}=t(P_n)\).

Una pregunta natural es: ¿qué tan bien funciona el principio del plug-in? Suele ser muy bueno cuando la única información disponible de \(P\) es la muestra \(\textbf{x}\), bajo esta circunstancia \(\hat{\theta}=t(P_n)\) no puede ser superado como estimador de \(\theta=t(P)\), al menos no en el sentido asintótico de teoría estadística \((n\to\infty)\).

El principio del plug-in provee de una estimación más no habla de precisión, por ello usaremos el bootstrap para estudiar el sesgo y el error estándar del estimador plug-in \(\hat{\theta}=t(P_n)\), la maravilla del bootstrap es que produce errores estándar y sesgos de manera automática, sin importar que tan complicada es la función \(t(P)\).

Errores estándar y sus estimaciones

Los estadísticos como \(\hat{\theta}=t(P_n)\) suelen ser el primer paso en el análisis de datos, el siguiente paso es investigar la precisión de las estimaciones; el bootstrap es un método para calcular precisión de estimaciones que utiliza el principio del plug-in para estimar el error estándar de una estadística.

Ejemplo: el error estándar de una media

Supongamos que \(x\) es una variable aleatoria que toma valores en los reales con distribución de probabilidad P. Denotamos por \(\mu_P\) y \(\sigma_P^2\) la media y varianza de P,

\[\mu_P = E_P(x),\] \[\sigma_P^2=var_P(x)=E_P[(x-\mu_P)^2]\]

en la notación enfatizamos la dependencia de la media y varianza en la distribución \(P\).

Ahora, sea \((x_1,...,x_n)\) una muestra aleatoria de \(P\), de tamaño \(n\), la media de la muestra \(\bar{x}=\sum_{i=1}^nx_i/n\) tiene esperanza \(\mu_P\) y varianza \(\sigma_P^2/n\).

En palabras: la esperanza de \(\bar{x}\) es la misma que la esperanza de \(x\), pero la varianza de \(\bar{x}\) es \(1/n\) veces la varianza de \(x\), así que entre mayor es la \(n\) tenemos una mejor estimación de \(\mu_P\).

El error estándar es la desviación estándar de una estadística.

En el caso de la media \(\bar{x}\), el error estándar, que denotamos \(se_P(\bar{x})\), es la raíz de la varianza de \(\bar{x}\), \[se_P(\bar{x}) = [var_P(\bar{x})]^{1/2}= \sigma_P/ \sqrt{n}.\]

En este punto podemos usar el principio del plug-in, simplemente sustituimos \(P_n\) por \(P\) y obtenemos, primero, una estimación de \(\sigma_P\): \[\hat{\sigma}=\hat{\sigma}_{P_n} = \bigg\{\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2\bigg\}^{1/2}\]

de donde se sigue la estimación del error estándar: \[\hat{se}(\bar{x})=\hat{\sigma}_{P_n}/\sqrt{n}=\bigg\{\frac{1}{n^2}\sum_{i=1}^n(x_i-\bar{x})^2\bigg\}^{1/2}\]

Notemos que usamos el principio del plug-in en dos ocasiones, primero para estimar la esperanza \(\mu_P\) mediante \(\mu_{P_n}\) y luego para estimar el error estándar \(se_P(\bar{x})\). En el caso de la media \(\hat{\theta}=\bar{x}\) la aplicación del principio del plug-in para el cálculo de errores estándar es inmediata; sin embargo, hay estadísticas para las cuáles no es fácil aplicar este método y es ahí cuando recurrimos al bootstrap.

Antes de pasar al bootstrap podemos preguntarnos: ¿porqué tanto énfasis en el error estándar? El error estándar es la manera más común para describir la precisión de una estadística. En términos generales, esperamos que \(\bar{x}\) este a una distancia de \(\mu_P\) menor a un error estándar el 68% del tiempo, y a menos de 2 errores estándar el 95% del tiempo. Estos porcentajes están basados el teorema central del límite que nos dice que bajo ciertas condiciones (bastante generales) de \(P\) la distribución de \(\bar{x}\) se aproximará a una distribución normal: \[\bar{x} \overset{\cdot}{\sim} N(\mu_P,\sigma_P^2/n)\]

Este html puede ser interactivo, para ver la aplicación de shiny debes correr el Rmd de manera local (Run Document). Si no deseas correrlo o no tienes el Rmd puedes ver la aplicación en (shinyapps.io)[https://tereom.shinyapps.io/15-TLC/].

El estimador bootstrap del error estándar

Supongamos que tenemos una muestra aleatoria \(\textbf{x}=(x_1,x_2,...,x_n)\) proveniente de una distribución de probabilidad desconocida \(P_n\) y deseamos estimar un parámetro \(\theta = t(P)\) con base en la muestra. Para esto, calculamos una estimación \(\hat{\theta}=s(\textbf{x})\) (la estimación puede ser la estimación plug-in \(t(P_n)\) pero también puede ser otra). Entonces podemos usar bootstrap para calcular el error estándar de la estimación.

Definimos una muestra bootstrap como una muestra aleatoria de tamaño \(n\) que se obtiene de la distribución empírica \(P_n\) y la denotamos \[\textbf{x}^* = (x_1^*,...,x_n^*).\]

La notación de estrella indica que \(\textbf{x}^*\) no son los datos \(\textbf{x}\) sino una versión de remuestreo de \(\textbf{x}\).

Otra manera de frasearlo: Los datos bootsrtap \(x_1^*,...,x_n^*\) son una muestra aleatoria de tamaño \(n\) seleccionada con reemplazo de la población de \(n\) objetos \((x_1,...,x_n)\).

Ahora, a cada muestra bootstrap \(\textbf{x}^*\) le corresponde una replicación \(\hat{\theta}^*=s(\textbf{x}^*).\)

La estimación bootstrap de \(se_{P}(\hat{\theta}^*)\), esto es, el error estándar de un estadístico \(\hat{\theta}\) es una estimación plug-in en donde la distribución empírica \(P_n\) toma el lugar de la distribución desconocida \(P\): el estimador bootstrap de \(se_P(\hat{\theta})\) se define como: \[se_{P_n}(\hat{\theta}^*)\] en otras palabras, la estimación bootstrap de \(se_P(\hat{\theta})\) es el error estándar de \(\hat{\theta}\) para conjuntos de datos de tamaño \(n\) seleccionados de manera aleatoria de \(P_n\).

La fórmula \(se_{P_n}(\hat{\theta}^*)\) no existe para casi ninguna estimación que diferente de la media, por lo que recurrimos a la técnica computacional bootstrap: el algoritmo funciona seleccionando distintas muestras bootstrap, evaluando la replicación bootstrap correspondiente y estimando el error estándar de \(\hat{\theta}\) mediante la desviación estándar empírica de las replicaciones. El resultado es la estimación bootstrap del error estándar, que denotamos \(\hat{se}_B\), donde \(B\) es el número de muestras bootstrap usadas.

Algoritmo bootstrap para estimar errores estándar

  1. Selecciona \(B\) muestras bootsrtap independientes: \[\textbf{x}^{*1},..., \textbf{x}^{*B}\].

  2. Evalúa la replicación bootstrap correspondiente a cada muestra bootstrap: \[\hat{\theta}^{*b}b=s(\textbf{x}^{*b})\] para \(b=1,2,...,B.\)

  3. Estima el error estándar \(se_P(\hat{\theta})\) usando la desviación estándar muestral de las \(B\) replicaciones: \[\hat{se}_B = \bigg\{\frac{\sum_{b=1}^B[\hat{\theta}^{*}(b)-\hat{\theta}^*(\cdot)]^2 }{B-1}\bigg\}^{1/2}\]

donde \[\hat{\theta}^*(\cdot)=\sum_{b=1}^B \theta^{*}(b)/B \].

Notemos que la estimación bootstrap de \(se_{P}(\hat{\theta})\), el error estándar de una estadística \(\hat{\theta}\), es un estimador plug-in que usa la función de distribución empírica \(P_n\) en lugar de la distribución desconocida \(P\). En otras palabras, la estimación bootstrap de \(se_{P}(\hat{\theta})\) es el error estándar de \(\hat{\theta}\) para conjuntos de datos de tamaño \(n\) seleccionados aleatoriamente de \(P_n\).

Conforme el número de replicaciones \(B\) aumenta \[\hat{se}_B\approx se_{P_n}(\hat{\theta})\] este hecho equivale a decir que la desviación estándar empírica se acerca a la desviación estándar poblacional conforme crece el número de muestras. La población en este caso es la población de valores \(\hat{\theta}^*=s(x^*)\).

Al estimador de bootstrap ideal \(se_{P_n}(\hat{\theta})\) y su aproximación \(\hat{se}_B\) se les denota estimadores bootstrap no paramétricos ya que estan basados en \(P_n\), el estimador no paramétrico de la población \(P\).ºº

Ejemplo: Escribimos una función para calcular el error estándar de una media usando replicaciones bootstrap:

mediaBoot <- function(x){ 
  # x: variable de interés
  # n: número de replicaciones bootstrap
  n <- length(x)
  muestra_boot <- sample(x, size = n, replace = TRUE)
  mean(muestra_boot) # replicacion bootstrap de theta_gorro
}
thetas_boot <- rerun(500, mediaBoot(prim_muestra$esp3)) %>% flatten_dbl()
sd(thetas_boot)
#> [1] 4.05

y se compara con \(\hat{se}(\bar{x})\) (estimador plug-in del error estándar):

se <- function(x) sqrt(sum((x - mean(x)) ^ 2)) / length(x)
se(prim_muestra$esp3)
#> [1] 4.14

Nota: Conforme \(B\) aumenta \(\hat{se}_{B}(\bar{x})\to \{\sum_{i=1}^n(x_i - \bar{x})^2 / n \}^{1/2}\), se demuestra con la ley débil de los grandes números.

Considera el coeficiente de correlación muestral entre la calificación de \(y=\)español 3 y la de \(z=\)español 6: \(\hat{corr}(y,z)=0.9\). ¿Qué tan preciso es esta estimación?

¿Cuántas replicaciones bootstrap (B)?

La estimación bootstrap ideal es un resultado asintótico \(B=\infty\), en esta caso \(\hat{se}_B\) iguala la estimación plug-in \(se_{P_n}\). En la práctica para elegir el tamaño de \(B\) debemos considerar que buscamos las mismas propiedades para la estimación de un error esrándar que para cualquier estimación: poco sesgo y desviación estándar chica. El sesgo de la estimación bootstrap del error estándar suele ser bajo y el error estándar está

Reglas de dedo (Effron y Tibshirani):

  1. Incluso un número chico de replicaciones bootstrap, digamos \(B=25\) es informativo, y \(B=50\) con frecuencia es suficiente para dar una buena estimación de \(se_P(\hat{\theta})\).

  2. En pocos casos es necesario realizar más de \(B=200\) replicaciones cuando se busca estimar error estándar.

seMediaBoot <- function(x, B){
    thetas_boot <- rerun(B, mediaBoot(x)) %>% flatten_dbl()
    sd(thetas_boot)
}

B_muestras <- data_frame(n_sims = c(5, 25, 50, 100, 200, 400, 1000, 1500, 3000)) %>% 
    mutate(est = map_dbl(n_sims, ~seMediaBoot(x = prim_muestra$esp3, B = .)))
B_muestras
#> # A tibble: 9 x 2
#>   n_sims   est
#>    <dbl> <dbl>
#> 1      5  3.91
#> 2     25  3.40
#> 3     50  4.10
#> 4    100  4.50
#> 5    200  3.81
#> 6    400  3.94
#> 7   1000  3.91
#> 8   1500  4.15
#> 9   3000  4.19

Ejemplos

Componentes principales: calificaciones en exámenes

Los datos marks (Mardia, Kent y Bibby, 1979) contienen los puntajes de 88 estudiantes en 5 pruebas: mecánica, vectores, álgebra, análisis y estadística. Cada renglón corresponde a la calificación de un estudiante en cada prueba.

marks <- read_csv("datos/marks.csv")
#> Parsed with column specification:
#> cols(
#>   id = col_integer(),
#>   MECH = col_integer(),
#>   VECT = col_integer(),
#>   ALG = col_integer(),
#>   ANL = col_integer(),
#>   STAT = col_integer()
#> )
glimpse(marks)
#> Observations: 88
#> Variables: 6
#> $ id   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17...
#> $ MECH <int> 77, 63, 75, 55, 63, 53, 51, 59, 62, 64, 52, 55, 50, 65, 3...
#> $ VECT <int> 82, 78, 73, 72, 63, 61, 67, 70, 60, 72, 64, 67, 50, 63, 5...
#> $ ALG  <int> 67, 80, 71, 63, 65, 72, 65, 68, 58, 60, 60, 59, 64, 58, 6...
#> $ ANL  <int> 67, 70, 66, 70, 70, 64, 65, 62, 62, 62, 63, 62, 55, 56, 5...
#> $ STAT <int> 81, 81, 81, 68, 63, 73, 68, 56, 70, 45, 54, 44, 63, 37, 7...
marks <- select(marks, -id)

Entonces un análisis de componentes principales proseguiría como sigue:

pc_marks <- princomp(marks)
summary(pc_marks)
#> Importance of components:
#>                        Comp.1 Comp.2  Comp.3 Comp.4 Comp.5
#> Standard deviation     26.061 14.136 10.1276 9.1471  5.638
#> Proportion of Variance  0.619  0.182  0.0935 0.0763  0.029
#> Cumulative Proportion   0.619  0.801  0.8948 0.9710  1.000
loadings(pc_marks)
#> 
#> Loadings:
#>      Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
#> MECH -0.505  0.749  0.300  0.296       
#> VECT -0.368  0.207 -0.416 -0.783 -0.189
#> ALG  -0.346        -0.145         0.924
#> ANL  -0.451 -0.301 -0.597  0.518 -0.286
#> STAT -0.535 -0.548  0.600 -0.176 -0.151
#> 
#>                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
#> SS loadings       1.0    1.0    1.0    1.0    1.0
#> Proportion Var    0.2    0.2    0.2    0.2    0.2
#> Cumulative Var    0.2    0.4    0.6    0.8    1.0
plot(pc_marks, type = "lines")

biplot(pc_marks)

Los cálculos de un análisis de componentes principales involucran la matriz de covarianzas empírica \(G\) (estimaciones plug-in)

\[G_{jk} = \frac{1}{88}\sum_{i=1}^88(x_{ij}-\bar{x_j})(x_{ik}-\bar{x_k})\]

para \(j,k=1,2,3,4,5\), y donde \(\bar{x_j} = \sum_{i=1}^88 x_{ij} / 88\) (la media de la i-ésima columna).

G <- cov(marks) * 87 / 88
G
#>      MECH  VECT   ALG   ANL  STAT
#> MECH  302 125.8 100.4 105.1 116.1
#> VECT  126 170.9  84.2  93.6  97.9
#> ALG   100  84.2 111.6 110.8 120.5
#> ANL   105  93.6 110.8 217.9 153.8
#> STAT  116  97.9 120.5 153.8 294.4

Los pesos y las componentes principales no son mas que los eigenvalores y eigenvectores de la matriz de covarianzas \(G\), estos se calculan a través de una serie de de manipulaciones algebraicas que requieren cálculos del orden de p^3 (cuando G es una matriz de tamaño p\(\times\)p).

eigen_G <- eigen(G)
lambda <- eigen_G$values
v <- eigen_G$vectors
lambda
#> [1] 679.2 199.8 102.6  83.7  31.8
v
#>        [,1]    [,2]   [,3]     [,4]    [,5]
#> [1,] -0.505  0.7487  0.300  0.29618 -0.0794
#> [2,] -0.368  0.2074 -0.416 -0.78289 -0.1889
#> [3,] -0.346 -0.0759 -0.145 -0.00324  0.9239
#> [4,] -0.451 -0.3009 -0.597  0.51814 -0.2855
#> [5,] -0.535 -0.5478  0.600 -0.17573 -0.1512
  1. Proponemos el siguiente modelo simple para puntajes correlacionados:

\[\textbf{x}_i = Q_i \textbf{v}\]

donde \(\textbf{x}_i\) es la tupla de calificaciones del i-ésimo estudiante, \(Q_i\) es un número que representa la habilidad del estudiante y \(\textbf{v}\) es un vector fijo con 5 números que aplica a todos los estudiantes. Si este modelo simple fuera cierto, entonces únicamente el \(\hat{\lambda}_1\) sería positivo y \(\textbf{v} = \hat{v}_1\). Sea \[\hat{\theta}=\frac{\hat{\lambda}_1}{\sum_{i=1}^5\hat{\lambda}_i}\] el modelo propuesto es equivalente a \(\hat{\theta}=1\), inculso si el modelo es correcto, no esperamos que \(\hat{\theta}\) sea exactamente uno pues hay ruido en los datos.

theta_hat <- lambda[1]/sum(lambda)
theta_hat
#> [1] 0.619

El valor de \(\hat{\theta}\) mide el porcentaje de la varianza explicada en la primer componente principal, ¿qué tan preciso es \(\hat{\theta}\)? La complejidad matemática en el cálculo de \(\hat{\theta}\) es irrelevante siempre y cuando podamos calcular \(\hat{\theta}^*\) para una muestra bootstrap, en esta caso una muestra bootsrtap es una base de datos de 88$$5 \(\textbf{X}^*\), donde las filas \(\textbf{x_i}^*\) de \(\textbf{X}^*\) son una muestra aleatoria de tamaño 88 de la verdadera matriz de datos.

pc_boot <- function(){
    muestra_boot <- sample_n(marks, size = 88, replace = TRUE)
    G <- cov(muestra_boot) * 87 / 88 
    eigen_G <- eigen(G)
    theta_hat <- eigen_G$values[1] / sum(eigen_G$values)
}
B <- 1000
thetas_boot <- rerun(B, pc_boot()) %>% flatten_dbl()

Veamos un histograma de las replicaciones de \(\hat{\theta}\):

ggplot(data_frame(theta = thetas_boot)) +
    geom_histogram(aes(x = theta, y = ..density..), binwidth = 0.02, fill = "gray40") + 
    geom_vline(aes(xintercept = mean(theta)), color = "red") +
    labs(x = expression(hat(theta)^"*"), y = "")

Estas tienen un error estándar

theta_se <- sd(thetas_boot)
theta_se
#> [1] 0.0476

y media

mean(thetas_boot)
#> [1] 0.624

la media de las replicaciones es muy similar a la estimación \(\hat{\theta}\), esto indica que \(\hat{\theta}\) es cercano a insesgado. El intervalo de confianza estándar para el verdadero valor \({\theta}\) es:

\[\theta \in (\hat{\theta}- z^{(1-\alpha)}\cdot \hat{se}, \hat{\theta}-z^{(\alpha)}\cdot \hat{se})\] con un nivel de confianza de \(1-2\alpha\).

Entonces:

# 0.68
theta_hat - theta_se
#> [1] 0.571
theta_hat + theta_se
#> [1] 0.667

# 0.9
theta_hat - qnorm(0.90) * theta_se
#> [1] 0.558
theta_hat + qnorm(0.90) * theta_se
#> [1] 0.68
  1. El eigenvetor \(\hat{v}_1\) correspondiente al mayor eigenvalor se conoce como primera componente de \(G\), supongamos que deseamos resumir la calificación de los estudiantes mediante un único número, entonces la mejor combinación lineal de los puntajes es

\[y_i = \sum_{k = 1}^5 \hat{v}_{1k}x_{ik}\]

esto es, la combinación lineal que utiliza las componentes de \(\hat{v}_1\) como ponderadores. Si queremos un resumen compuesto por dos números \((y_i,z_i)\), la segunda combinación lineal debería ser:

\[z_i = \sum_{k = 1}^5 \hat{v}_{2k}x_{ik}\]

Las componentes principales \(\hat{v}_1\) y \(\hat{v}_2\) son estadísticos, usa bootstrap para dar una medición de su variabilidad.

Mejorando el código

El código que usamos para el ejemplo de componentes principales arriba es un poco lento comparemos el tiempo que tarda con una alternativa usando la función bootstrap() del paquete broom.

Usamos system.time() para comparar tiempos de ejecución, esta función calcula el tiempo en segundos que toma ejecutar una expresión (si hay un error, regresa el tiempo hasta que ocurre el error):

La función system.time supone que sabes donde buscar, es decir, que expresiones debes evaluar, una función que puede ser más útil cuando uno desconoce cuál es la función que alenta un programa es profvis() (paquete `profvis).

library(profvis)

profvis({
    pcBoot <- function(){
        muestra_boot <- sample_n(marks, size = 88, replace = TRUE)
        G <- cov(muestra_boot) * 87 / 88 
        eigen_G <- eigen(G)
        theta_hat <- eigen_G$values[1] / sum(eigen_G$values)
    }
    B <- 500
    thetas_boot_p <- rerun(B, pcBoot())
})

profvis() utiliza a su vez la función Rprof() de R base, este es un perfilador de muestreo que registra cambios en la pila de funciones, funciona tomando muestras a intervalos regulares y tabula cuánto tiempo se lleva en cada función.