options(scipen = 999)
set.seed(123)
# cargamos las librerias a utilizar
pacman::p_load(tidyverse, magrittr)
## Iniciaremos trabajando con toda la población ##
poblacion <- as.data.frame(read.csv('data/covid.csv')) # Cargamos la base poblacional2 Repaso de Estadística
La base necesaria para la econometría consiste en tener un entendimiento importante de algunos fundamentos de probabilidad e inferencia estadística. Esta sección hace un repaso no exhaustivo de algunos conceptos fundamentales que emplearemos. Se recomienda, sin embargo, que aquellos que no se sientan cómodos con alguos fundamentos básicos hagan una revisión de dichos conceptos para fortalecer su entendimmento de los temas vistos posteriormente. Crash course statistics es una referencia útil y fácilmente accesible para aquellos que quieran hacer un repaso de conceptos básicos de estadística.
2.1 Probabilidad
Esta sección consiste es un repaso muy breve y general de conceptos de probabilidad que se utilizarán a lo largo del curso1. En general, cuando hablamos de probabilidad es usual distinguir entre variables discretas y continuas. En términos prácticos, las variables discretas son aquellas que toman un pequeño número de valores posibles. Por su parte, las variables continuas son aquellas que pueden tener un número infinito de valores posibles. Usualmente, se les identifica como variables que toman un valor dentro de cierto rango de valores posibles.
Las variables suelen ser descritas por una función de densidad (pdf), la cual describe la probabilidad de que una variable tome cierto valor: \(f(x_j) = P(X=x_j) = p_j\) . Dicha función de densidad suele estar ligada a la función de densidad acumulada (cdf), que describe la probabilidad de que una variable tenga un valor menor o igual a cierto número: \(F(x_j) = P(X\leq x_j)\).
Cuando involucramos en nuestro análisis más de una variable, es común hacer referencia a la función de densidad conjunta: \(f_{X,Y}(x_j,y_k)= P(X = x_j,Y = y_k)\) y a la función de densidad marginal: \(f_X(x_j) = P(X = x_j )\). En el caso en que las variables \(X\) y \(Y\) sean independientes: \(f_{X,Y}=(x_j, y_k) = f_X(x_j)f_Y(y_k)\). Asimismo, con más de una variable, será común referirnos a la densidad condicional: \(f_{Y|X}(y_k|x_j) = \frac{f_{X,Y}(x_j,y_k)}{f_X(x_j)}\) . [En adelante, para simplificar únicamente utilizaremos \(x\) en vez de \(x_j\) y \(y\) en vez de \(y_k\). El único propósito de utilizar \(x_j\) y \(y_k\) era para indicar que son valores específicos, i.e. algún número. En adelante, cuando utilicemos minúsculas estaremos haciendo referencia a valores específicos.]
Valor esperado. El valor esperado es una medida de tendencia central. Si \(X\) es una variable discreta que toma \(k\) distintos posibles valores, tendremos que: \[E[X]=\sum_{j=1}^kx_jf(x_j)\]
Asimismo, en el caso continuo: \[E[X]=\int_{-\infty}^{\infty}xf(x)\] Propiedades del valor esperado:
\(E(a) = a\), donde \(a\) es una constante
\(E(aX + bY + c) = aE(X) + bE(Y ) + c\), donde \(a\), \(b\), \(c\) son constantes
Varianza. La varianza describe, en promedio, que tan lejos suele estar una variable del valor esperado (\(\mu_X = E(X)\)). \[\sigma^2_X= Var(X) = E[(X-\mu_X)^2] = E(X^2)-\mu_X^2\] La desviación estándar es simplemente: \[\sigma_X=\sqrt{Var(X)}\]
Propiedades de la varianza:
\(Var(a) = 0\)
\(Var(aX+b)=a^2Var(X)\)
\(Var(aX\pm bY)=a^2Var(X)+b^2Var(Y)\pm ab Cov(X,Y)\)
Covarianza. La covarianza mide la relación entre dos variables. Una covarianza positiva indica que ambas variables suelen moverse en la misma dirección. Una covarianza negativa indica lo contrario \((\mu_Y = E(Y ))\) \[\sigma_{X,Y} = Cov(X, Y ) = E[(X-\mu_X)(Y-\mu_Y)]\] Propiedades de la covarianza:
\(Cov(X,Y) = 0\) si \(X\) y \(Y\) son independientes
\(Cov(aX + b, cY + d) = acCov(X, Y )\), donde \(a\), \(b\), \(c\), \(d\) son constantes
El coeficiente de correlación es simplemente: \(\rho_{X,Y} = corr(X,Y) = \frac{Cov(X,Y )}{\sigma_X\sigma_Y}\)
2.2 Notación
El análisis econométrico empírico generalmente involucra llevar a cabo una inferencia a partir de los datos a los cuales se tiene acceso (la muestra) acerca del comportamiento de cierta población. Obtener muestras representativas de la población acerca de la cual se busca aprender algo es fundamental y suele ser en un tema a desarrollar por sí solo. A lo largo del curso daremos esta pregunta por resuelta y supondremos que tenemos acceso a datos de muestras representativas de cierta población.
La población es un grupo definido de unidades de observación que pueden estar o no restringidas de distintas formas (por ejemplo: geográfica, económica o demográficamente). Una unidad de observación puede ser desde una acción, como por ejemplo, un mensaje de texto enviado vía celular o una ruta específica de un punto a otro que una persona se traslada, hasta un ente con una identidad específica que puede ser desde individuos, empresas, ciudades o países. En resumen, la población es el grupo objetivo que se pretende estudiar. En gran parte de los estudios suele establecerse que nos interesa conocer algún aspecto específico de la población, típicamente un estadístico. Por ejemplo, es común indicar que nos interesa conocer la media de una variable o la correlación entre dos variables de dicha población. Siempre que nos interese un estadístico en específico, nos referiremos a éste como el parámetro poblacional de interés.
Por su parte, la muestra es un subconjunto de la población el cual será utilizado en el análisis empírico para poder decir algo acerca de la población. La muestra será aquella de la cual recopilaremos información que tendremos concentrada en una base de datos que utilizaremos con programas estadísticos para llevar a cabo el análisis. La inferencia generalmente conlleva hacer una estimación utilizando datos de la muestra con el objetivo de obtener conclusiones acerca de la población. Por ejemplo, utilizando diversas estrategias se plantean estimadores del parámetro poblacional, mismos que emplean datos de una muestra para dar evidencia e información específica de dicho parámetro.
Es importante conocer los detalles de la estrategia utilizada para recabar la muestra, ya que esto es indicativo de si dicha muestra es representativa de la población que se pretende2. Algunas causas de sesgos muestrales comunes incluyen sesgo de respuesta (i.e. el hecho de que las personas que no contestan no son personas al azar), problemas comprendiendo las preguntas y problemas encontrando ciertos tipos de individuos en la población. Estos problemas suelen llevar a sesgos muestrales, lo cual significa que la muestra puede no ser representativa de la población objetivo. Para evaluar esto, es útil tener una base de datos de comparación (i.e. un benchmark) que te permita determinar si las características de tu muestra son similares a aquellas características de la población objetivo. Los censos nacionales son un buen ejemplo de este tipo de benchmarks.
En estadística, generalmente nos referiremos al uso de muestras aleatorias debido a que eso nos asegura representatividad. En las muestras aleatorias cada uno de sus componentes se selecciona de forma independiente y proviene de una distribución común \(\{Y_1,\dots,Y_n\}\). En este caso se dice que \(Y_i\) es una variable aleatoria independiente e idénticamente distribuida (i.i.d.). Las variables aleatorias \(\{Y_1,\dots,Y_n\}\) son variables desconocidas. Una vez que la muestra es recabada tendremos un conjunto de números \(\{y_1,\dots, y_n\}\), los cuales utilizaremos para llevar a cabo inferencia.
Ejemplo: Supón que nuestra variable de interés es la media de la edad de los alumnos cursando educación superior en México.
¿Cuál es la población?
¿Cómo recabarías una muestra para poder obtener/hacer una estimación válida?
¿Qué sucedería si yo recabo una muestra basada en los alumnos del ITAM o de la clase?
A continuación buscamos utilizar nuestra muestra para conocer algo acerca del parámetro. Para ello plantearemos un estimador. El estimador en sí es una fórmula que utiliza como insumo los elementos de la muestra y entrega como resultado un valor que pretende aproximar al parámetro poblacional objetivo. Es incorrecto pensar en que el estimador tiene valores específicos (a pesar de que el parámetro poblacional si lo tenga) dado que el estimador en sí es una estrategia que se plantea para aproximar lo mejor posible al parámetro poblacional. Al utilizar los valores específicos de la base de datos a la cual tenemos acceso y aplicar la fórmula del estimador para obtener un valor específico, estaremos hablando de un valor estimado. Para calificar si un estimador es apropiado y si es mejor o peor que otro, solemos hablar de propiedades que nos ayudan a calificarlo, como insesgadez y eficiencia.
Ejemplo: Continuando con el ejemplo de la media de edad de alumnos en educación superior en México: supongamos a partir de aquí que se recaba una muestra aleatoria \(\{Y_1,\dots,Y_n\}\)
Parámetro: \(\mu\). Es un número específico que típicamente no es conocido. En este caso, es la media de la edad de todos los alumnos cursando educación media superior en México.
Estimador: \(W=h(Y_1,\dots,Y_n)\). Es una fórmula que genera un estimador del parametro.
Valor estimado: \(w=h(y_1,\dots,y_n)\). Es un número específico que resulta de aplicar la fórmula del estimador a la muestra recabada.
Definimos a continuación las propiedades deseables de los estimadores:
El estimador (\(W\)) de un parámetro (\(\mu\)) es insesgado si su valor esperado es igual al parámetro: \(E(W)=\mu\). Si un estimador es insesgado, esto no quiere decir que el valor estimado será igual al parámetro (o incluso cercano a éste), ya que esto dependerá de la muestra que sea recabada.
Un estimador (\(W_1\)) es más eficiente relativamente a otro (\(W_2\)) si \(Var(W_1)\leq Var(W_2)\).
Ejemplo: Consideremos dos estimadores para el parámetro \(\mu\):
El valor promedio: \(W_1=\frac{1}{n}\sum_{i=1}^n Y_i=\bar{Y}\)
El promedio de edad de dos personas de mi muestra elegidas al azar: \(W_2=\frac{Y_a+Y_b}{2}\)
¿Son estos estimadores insesgados?, ¿cuál es más eficiente?, ¿qué valor estimado estará más cerca de la media poblacional?
2.3 Propiedades asintóticas de los estimadores
En el ejemplo anterior podemos deducir que el estimador \(W_1\) se vuelve más eficiente conforme el tamaño de muestra aumenta. Las propiedades asintóticas de los estimadores son aquellas que aplican cuando se tienen muestras grandes. Sin embargo, no es claro de qué tamaño necesita ser el número de observaciones (\(n\)) para que la muestra sea considerada como grande y sea correcto aplicar las propiedades asintóticas a los estimadores. Generalmente, esto depende de la distribución poblacional de la variable de interés, pero en la mayoría de los casos en los que utilizamos encuestas, aplicar propiedades asintóticas será razonable.
Ejemplo: Dar un ejemplo mostrando dos distribuciones normales, una más dispersa que la otra.
Consistencia. La consistencia se refiere al comportamiento de la distribución muestral del estimador conforme el tamaño de la muestra se incrementa. Conforme aumentamos el tamaño de la muestra, la distribución de \(W_1\) se volverá más concentrada alrededor de \(\mu\). Por lo tanto, menos probable será que un valor estimado se ubique lejos de \(\mu\).
Ejemplo: ¿Es el estimador \(W_2\) consistente?
Ley de Grandes Números (LGN). La ley de grandes números nos dice que si queremos aproximarnos a la media poblacional, podemos hacerlo en gran medida si elegimos muestras lo suficientemente grandes y utilizamos el estimador del promedio. Sin embargo, utilizando la LGN únicamente obtenemos estimadores puntuales y no tenemos información acerca de su distribución.
Teorema Central del Límite (TCL). Sea \(\{Y_1,\dots,Y_n\}\) una muestra aleatoria con media \(\mu\) y varianza \(\sigma^2\). Entonces, \[\begin{equation} Z_n=\frac{\sqrt{n}(\bar{Y}_n-\mu)}{\sigma}\xrightarrow{d}N(0,1) \text{ conforme $n\xrightarrow{}\infty$} \end{equation}\]
Intuitivamente, este resultado indica que, sin importar la distribución poblacional de \(Y\), la distribución de la variable \(Z_n\) (que es la versión estandarizada de \(\bar{Y}_n\)) se aproxima en gran medida a una distribución normal estándar (\(N(0,1)\)) conforme el tamaño de la muestra (\(n\)) aumenta3.
2.4 Pruebas de hipótesis
En la gran mayoría de las aplicaciones empíricas de econometría tendremos que llevar a cabo pruebas de hipótesis. En dichas pruebas de hipótesis es importante notar que evaluamos si el parámetro es igual a cierto valor en el caso de la hipótesis nula. Esto es particularmente adecuado debido a que el parámetro es un número específico que es (usualmente) desconocido para el econometrista. Generalmente, durante el curso asumiremos que las muestras son grandes y por tanto podemos aplicar propiedades asintóticas, Esto quiere decir que en la mayoría de los casos podremos utilizar la distribución normal y ji-cuadrada para llevar a cabo las pruebas de hipótesis.
Para repasar cómo llevar a cabo pruebas de hipótesis supongamos que estamos interesados en evaluar si la media de edad de los alumnos en educación media superior en México es igual a 20. Cabe señalar que nuestra hipótesis es acerca del valor de un parámetro poblacional y utilizaremos una muestra para evaluar dicha hipótesis. La hipótesis nula debe establecerse como igualdad debido a que en la evaluación de la hipótesis se asume que es cierta y eso genera una distribución para el estadístico que se utilizará. En nuestro ejemplo, establecemos la siguiente hipótesis nula:
\[\begin{equation*} H_0: \mu=20 \end{equation*}\]La hipótesis alternativa se establece para especificar la zona de rechazo de la hipótesis nula. Generalmente, en una prueba se busca rechazar la hipótesis nula en favor de la alternativa. En nuestro caso, la hipótesis alternativa será:
\[\begin{equation*} H_1: \mu\neq 20 \end{equation*}\]La hipótesis alternativa puede también ser establecida como \(\mu>\) 20 (o \(\mu<\) 20). Este sería el caso si lo que nos interesa es evaluar si la media poblacional es mayor (menor) a 20. Es importante recordar que como resultado de la prueba de hipótesis, la hipótesis nula puede ser rechazada o no rechazada. Sin embargo, es incorrecto decir que es aceptada. Formalmente, en el ejemplo anterior podríamos concluir ya sea que: (i) hay evidencia suficiente para rechazar que la media poblacional es igual a 20 con \(x\%\) de significancia, o que (ii) no hay evidencia suficiente para rechazar que la media poblacional es igual a 20 con x% de significancia.
Para establecer el nivel de significancia \(x\%\) (o alternativamente el nivel de confianza \([1-x]\%\)), hay que tomar en cuenta los dos tipos de errores que podemos cometer al evaluar pruebas de hipótesis:
Error tipo I: Podemos rechazar la hipótesis nula siendo que esta es verdadera
Error tipo II: Podemos no rechazar la hipótesis nula siendo esta falsa
Típicamente, el nivel de significancia se establece basado en el error tipo I, que generalmente busca reducirse en las pruebas de hipótesis. Dada nuestra notación, el nivel de significancia se define como:
\[x\%=Pr(\text{Rechazar}\phantom{t}H_0|H_0)=Pr(\text{Error tipo I})\]
El valor de \(x\%\) será una valor que tendremos que asumir para llevar a cabo la prueba de hipótesis. El valor más común es de 0.05 (o 5%) de significancia, seguido de 0.01 y 0.1 (lo cual es equivalente a 95%, 99% y 90% de nivel de confianza, respectivamente).
El error tipo II suele estar ligado al poder estadístico. Más adelante en el curso discutiremos cómo utilizar el poder estadístico para determinar el número de observaciones que se requieren para llevar a cabo un análisis estadístico experimental.
Supongamos por el momento que elegimos un nivel de significancia de 5% y que nuestra muestra de estudiantes mexicanos es aleatoria y consta de 1000 individuos. Supongamos que el promedio muestral de edad es de 21.5 y la varianza de las edades de 500.
Tenemos 3 alternativas para llevar a cabo la prueba de hipótesis:
2.4.1 Estadístico \(t\)
Para utilizar este método utilizamos la media y la desviación estándar estimada. La idea es asumir que la hipótesis nula es verdadera. Dado que asumimos esto, querremos determinar qué tan probable es que obtengamos un valor estimado \(\bar{y}=\) 21.5, dado que proviene de una distribución de la variable aleatoria \(\bar{Y}\) que es normal (por propiedades asintóticas) con media 20 y desviación estándar \(\sqrt{500/1000}\).
Tomando dichos supuestos y aplicando el TCL estandarizamos \(\bar{Y}\), con lo cual derivamos nuestro estadístico \(t\), mismo que tendrá una distribución \(N(0,1)\):
\[\begin{equation} t=\frac{(\bar{Y}-\mu_0)}{\sqrt{S^2/n}} \longrightarrow~N(0,1) \end{equation}\]Y finalmente, de dicha distribución obtenemos un valor: \[\begin{equation*} \widehat{t}=\frac{(21.5-20)}{\sqrt{500/1000}}=2.1213 \end{equation*}\]
Una vez que tenemos dicho valor utilizamos la distribución de la normal estándar para comparamos este valor con un referente que esté en el límite de ser razonable (el valor crítico). Para ello empleamos el nivel de significancia. Utilizando un 5% de significancia determinamos que valor crítico (en términos absolutos) representa el 95% del cdf de la distribución normal estándar. Dicho valor (1.96) se compara con el valor estimado para determinar qué tan probable es observarlo dada la distribución de la cual asumimos que proviene (debido a que supusimos que la hipótesis nula es cierta). Dado que el estadístico-\(t\) es mayor que el valor crítico rechazamos la hipótesis nula con un 5% de significancia a favor de la hipótesis alternativa.
2.4.2 Valor-p
El valor-\(p\) nos dice hasta qué nivel de significancia la hipótesis nula sería rechazada. Siempre que el nivel de significancia sea mayor al valor-\(p\), la hipótesis nula sería rechazada. Para determinar significancia a los niveles usuales, este valor usualmente se compara con 0.01, 0.05 y 0.1. Sin embargo, el valor-\(p\) tiene el significado en si mismo de informar que probabilidad existe de observar un valor igual o más extremo que el obtenido de la muestra. En nuestro ejemplo:
\[\begin{equation} \textit{valor-p}=2\cdot(1-F(|t|))=0.034 \end{equation}\]Es importante, tener en consideración que en el caso en que la hipótesis alternativa sea unilateral (one-sided), \(\textit{valor-p}=(1-F(|t|))\).
2.4.3 Intervalo de confianza
Si nuevamente utilizamos un nivel de significancia de 5%, necesitaremos determinar un intervalo de confianza del 95%. Dicho intervalo de confianza se genera de la siguiente forma:
\[\begin{equation} \begin{split} Pr\biggl(-1.96<\frac{\sqrt{n}(\bar{Y}-\mu)}{S}<1.96\biggl)&=0.95\\ Pr\biggl(\bar{Y}-1.96\cdot\frac{S}{\sqrt{n}}<\mu<\bar{Y}+1.96\cdot\frac{S}{\sqrt{n}}\biggl)&=0.95\\ \end{split} \end{equation}\]En el caso del intervalo de confianza es importante recordar que la incertidumbre radica en el intervalo dado que \(\bar{Y}\) es una variable aleatoria. El intervalo lo podemos interpretar como: de cada 100 muestras aleatorias que obtengamos, 95% de ellas tendrán al valor real del parámetro poblacional \(\mu\). No podemos decir que una vez que calculemos el intervalo, con 95% de probabilidad éste contendrá el valor real del parámetro. Recordemos que el parámetro es un valor específico (no aleatorio), por lo tanto, se encuentra o no en el intervalo.
En nuestro caso, el intervalo de 95% será: [20.114,22.88]. Dado que 20 no se encuentra dentro del intervalo, rechazamos la hipótesis nula con un 5% de significancia.
Para las pruebas unilaterales la alternativa al intervalo de confianza es el límite (o cota) inferior o superior. Imaginemos que en el caso de nuestro ejemplo nuestra hipótesis alternativa es:
\[\begin{equation*} H_1: \mu > 20 \end{equation*}\]En este caso nos interesaría comparar a el valor propuesto (20) con la cota inferior (\(C_I\)), dado que nuestro “intervalo” consistiría en: \([C_I,\infty)\). Para calcular la cota:
\[\begin{equation} \begin{split} Pr\biggl(\mu > \bar{Y}-1.64\cdot\frac{S}{\sqrt{n}} \biggl)&=0.95\\ Pr\biggl(\frac{\sqrt{n}(\bar{Y}-\mu)}{S} < 1.64 \biggl)&=0.95\\ \end{split} \end{equation}\]En el caso de nuestro ejemplo la cota inferior sería 20.34, y como dicho valor es mayor a 20, concluiríamos que la hipóitesis nula se rechazaría a favor de nuestra nueva hipótesis alternativa unilateral.
2.5 Bootstrap
En el caso de las pruebas de hipótesis anteriores estamos basando los resultados en el teorema central del límite. Para aplicar el TCL calculamos analíticamente la varianza de la media como:
\[\begin{equation} Var(\bar{Y})=\frac{\sigma^2}{n} \end{equation}\] donde \(\sigma^2\) es la varianza de \(Y_i\).
En el caso de TCL asumimos que una versión estandarizada de \(\bar{Y}\) (\(Z_n\)) converge en distribución a una normal estándar. Una alternativa a este procedimiento consiste en generar una distribución empírica de \(\bar{Y}\) y utilizar dicha distribución para calcular la varianza. Un problema con esta idea radica en que para generar una distribución empírica de \(\bar{Y}\) necesitamos varias observaciones de \(\bar{Y}\).
El método de bootstrap genera diversas observaciones partiendo de una muestra aleatoria \(\{Y_1,\dots,Y_n\}\) siguiendo los siguientes pasos:
Utilizando las observaciones de la muestra, elige una submuestra aleatoria de tamaño \(n\) (mismo tamaño que la muestra) con reemplazo. Esto quiere decir que habrá observaciones que se repitan más de una vez
Usando la submuestra calcula el estimador (\(\bar{Y}\) en nuestro ejemplo)
Repite los pasos anteriores \(M\) veces. Con esto tendrás \(M\) observaciones para \(\bar{Y}\): \(\{\bar{Y}_1,\dots,\bar{Y}_M\}\) y habrás generado una distribución empírica
Genera los estimadores: \[\begin{equation} \begin{split} E(\bar{Y})=&\frac{1}{M}\sum\limits_{k=1}^M \bar{Y}_k\\ Var(\bar{Y})=&\frac{1}{M}\sum\limits_{k=1}^M (\bar{Y}_k-E(\bar{Y}))^2 \end{split} \end{equation}\]
Utiliza dichos estimadores para llevar a cabo pruebas de hipótesis
Este método puede ser aplicado con la mayoría de los estimadores que veremos durante el curso. Es un método de gran utilidad siempre que sea dificil calcular una varianza para llevar a cabo pruebas de hipótesis. En particular podría utilizarse para calcular errores estándar de los coeficientes en una regresión. En dicho caso los pasos a seguir son los mismos que los descritos anteriormente. Lo que sucedería en el caso de una regresión de mínimos cuadrados ordinarios es que se llevaría a cabo una regresión con cada una de las submuestras elegidas en el primer paso. Con ello se obtendrían \(M\) posibles coeficientes para cada variable. El error estándar podría calcularse como la desviación estándar para cada uno de los coeficientes.
Existe también la posibilidad de utilizar una submuestra de tamaño menor al tamaño de la muestra original (\(n\)). En dicho caso, tendría que llevarse a cabo un ajuste al cálculo de la varianza. Supongamos que se eligen submuestras de tamaño \(L\). Todos los pasos serían los mismos que los descritos anteriormente, con la diferencia que el estimador de la varianza se calcularía como:
\[\begin{equation} Var(\bar{Y})=\frac{L}{n}\frac{1}{M}\sum\limits_{k=1}^M (\bar{Y}_k-E(\bar{Y}))^2 \end{equation}\]
Ejemplo: Bootstrap e intervalos de confianza para la media muestra en
R
A continuación se presenta un ejemplo de Bootstrap realizado en R. En el ejemplo utilizaremos datos sobre la pandemia de COVID-19 provenientes del sitio Our World in Data. La base contiene información sobre la pandemia hasta el día 27 de agosto de 2020 para 188 países. En este ejemplo nos interesará estimar la media de la tasa de contagios \(\left(TC=\frac{\text{Casos Confirmados}}{\text{Población}}\right)\).
Antes de comenzar revisemos de forma rápida las variables con las que contamos:
glimpse(poblacion)Rows: 188
Columns: 25
$ ISO_code <chr> "ABW", "AFG", "AGO", "ALB", "AND", "ARE", "…
$ continent <chr> "North America", "Asia", "Africa", "Europe"…
$ country <chr> "Aruba", "Afghanistan", "Angola", "Albania"…
$ confirmed <int> 1760, 38126, 2332, 8927, 1098, 68020, 37017…
$ deaths <int> 8, 1401, 103, 263, 53, 378, 7839, 861, 3, 5…
$ confirmed_per_million <dbl> 16484.649, 979.389, 70.954, 3102.022, 14210…
$ deaths_per_million <dbl> 74.930, 35.989, 3.134, 91.389, 685.951, 38.…
$ stringency_index <dbl> 62.04, NA, NA, 53.70, 41.67, NA, NA, NA, NA…
$ population <int> 106766, 38928341, 32866268, 2877800, 77265,…
$ population_density <dbl> 584.800, 54.422, 23.890, 104.871, 163.755, …
$ median_age <dbl> 41.2, 18.6, 16.8, 38.0, NA, 34.0, 31.9, 35.…
$ aged_65_older <dbl> 13.085, 2.581, 2.405, 13.188, NA, 1.144, 11…
$ aged_70_older <dbl> 7.452, 1.337, 1.362, 8.643, NA, 0.526, 7.44…
$ gdp_pc <dbl> 35973.781, 1803.987, 5819.495, 11803.431, N…
$ extreme_poverty <dbl> NA, NA, NA, 1.1, NA, NA, 0.6, 1.8, NA, 0.5,…
$ cardiovasc_death_rate <dbl> NA, 597.029, 276.045, 304.195, 109.135, 317…
$ diabetes_prevalence <dbl> 11.62, 9.59, 3.94, 10.08, 7.97, 17.26, 5.50…
$ female_smokers <dbl> NA, NA, NA, 7.1, 29.0, 1.2, 16.2, 1.5, NA, …
$ male_smokers <dbl> NA, NA, NA, 51.2, 37.8, 37.4, 27.7, 52.1, N…
$ handwashing_facilities <dbl> NA, 37.746, 26.664, NA, NA, NA, NA, 94.043,…
$ hospital_beds_per_thousand <dbl> NA, 0.500, NA, 2.890, NA, 1.200, 5.000, 4.2…
$ life_expectancy <dbl> 76.29, 64.83, 61.15, 78.57, 83.73, 77.97, 7…
$ people_tested <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ tests_performed <int> NA, NA, NA, NA, NA, 6755457, NA, NA, NA, 59…
$ tests_uu <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…Creamos la variable de Tasa de Contagios (TC) y visualizamos su distribución
# Creamos la variable de Tasa de Contagios (TC)
poblacion %<>% mutate(TC = confirmed / population)
poblacion %>% ggplot(aes(TC)) +
geom_histogram(aes(y = ..count.. / sum(..count..)),
fill = "palegreen4",
color = "white",
bins = 25) +
labs(x = "Tasa de Contagio",
y = "Proporción",
title = "Histograma para TC (Utilizando la población)") +
theme_minimal()
La media poblacional es igual a 0.004.
Trabajaremos ahora con una muestra de la población. Para esto crearemos una muestra aleatoria de tamaño 100 sin reemplazo a partir de la base poblacional y visualizamos la distribución muestral:
muestra <- sample_n(poblacion, size = 100, replace = F)
# Realizamos nuevamente un histograma
muestra %>% ggplot(aes(TC)) +
geom_histogram(aes(y = ..count.. / sum(..count..)),
fill = "slategray3",
color = "white",
bins = 25) +
labs(x = "Tasa de Contagio Muestral",
y = "Proporción",
title = "Histograma para TC (Utilizando la muestra)") +
theme_minimal()
La media muestral de \(TC\) es 0.004, la varianza muestral es 0.00004 y la varianza de la media muestral es 0.0000004. Realizamos ahora un Bootstrap de 1,000 submuestras del tamaño de la muestra original (\(n\) = 100):
meanBoot <- c() # Creamos un vector vacío que iremos
# llenando en el siguiente for
for (n in 1:1000){
meanBoot <- c(meanBoot, # n-1 medias de submuestras
mean(sample(muestra$TC, # calculamos la media de
size = 100, # la submuestra n
replace = T)))
}
meanBoot <- data.frame(mean = meanBoot) # Creamos un data frame con
# la media de cada una de las
# 1,000 submuestrasAhora creamos un histograma para la media muestral de \(TC\) (calculada a partir de bootstrap)
meanBoot %>% ggplot(aes(mean)) +
geom_histogram(aes(y = ..count../sum(..count..)),
fill = "lightblue",
color = "white",
bins = 20) +
geom_vline(xintercept = mean(meanBoot$mean), colour = "red4") +
geom_label(mapping = aes(x = mean(meanBoot$mean),
y = .05,
label = paste("Media pob =",
round(mean(poblacion$TC), 6)),
hjust = -.5, vjust =-.5),
colour = "red4") +
labs(x = "Medias muestrales (Bootstrap)",
y = "Proporción",
title = "Histograma de la media muestral") +
theme_minimal()
En la gráfica anterior la línea vertical muestra la media poblacional de \(TC\), dicho numero es nuestro objetivo.
Utilizando los resultados del bootstrap, podemos estimar la varianza de la media muestral de \(TC\): 0.0000004. Comparemos los valores de la varianza de la media muestral de \(TC\) que hemos obtenido:
var(muestra$TC) / length(muestra$TC)= 0.0000004var(meanBoot$mean)= 0.0000004
Ahora calcularemos intervalos de confianza al 95% para la media de TC considerando dos estrategias:
Usemos el estimador de la media muestral y el de la varianza de la media muestral:
llm_1 <- mean(muestra$TC, na.rm = T) - 1.96 * sqrt(var(muestra$TC) / length(muestra$TC))
ulm_1 <- mean(muestra$TC, na.rm = T) + 1.96 * sqrt(var(muestra$TC) / length(muestra$TC))El intervalo al 95% es igual a (llm_1, ulm_1) = (0.00314, 0.00562)
Usemos los resultados que obtuvimos en Bootstrap:
llm_2 <- mean(meanBoot$mean) - 1.96 * sqrt(var(meanBoot$mean))
ulm_2 <- mean(meanBoot$mean) + 1.96 * sqrt(var(meanBoot$mean))El intervalo al 95% es igual a (llm_2, ulm_2) = (0.00315, 0.00562)
Grafiquemos los intervalos de confianza anteriores:
lm <- c(llm_1, llm_2)
um <- c(ulm_1, ulm_2)
meanm <- c(mean(muestra$TC, na.rm = T),
mean(meanBoot$mean))
tICm <- c("Tradicional", "Bootstrap")
ICm <- data.frame(tipo = tICm,
lower = lm,
media = meanm,
upper = um)
ICm %>% ggplot(aes(x = media, y = tipo)) +
geom_point() +
geom_errorbarh(aes(xmin = lower,
xmax = upper),
size = 1,
height = 0.1) +
geom_vline(xintercept = mean(poblacion$TC),
colour = "red4") +
labs(x = "Tasa de Contagios",
y = "Método",
title = "Intervalos de Confianza al 95% para la media de TC") +
theme_minimal()
La línea vertical indica el valor de la media de \(TC\) poblacional
Ejemplo: Bootstrap e intervalos de confianza para otros estadísticos
Consideremos nuevamente la base de datos que utilizamos en el ejemplo anterior. Ahora estaremos interesados en calcular el ratio entre el percentil 90 y el percentil 10 \(\left(\varphi=\frac{P_{90}}{P_{10}}\right)\) para la tasa de mortalidad \(\left(TM=\frac{\text{Muertes}}{\text{Casos Confirmados}}\right)\). Este ratio nos indica cuántas veces mayor es la tasa de mortalidad del pais en el percentil 90 contra aquel del percentil 10.
Iniciaremos trabajando nuevamente con toda la población. Creamos la variable de Tasa de Mortalidad (TM) y visualizamos su distribución:
poblacion %<>% mutate(TM = deaths / confirmed)
poblacion %>% ggplot(aes(TM)) +
geom_histogram(aes(y = ..count.. / sum(..count..)),
fill = "palegreen4",
color = "white",
bins = 25) +
labs(x = "Tasa de Mortalidad",
y = "Proporción",
title = "Histograma para TM (Utilizando la población)") +
theme_minimal()
Calculamos ahora el ratio del percentil 90 sobre el percentil 10 (p90p/p10p):
p90p <- quantile(poblacion$TM,
probs = 0.9,
type = 2)
p10p <- quantile(poblacion$TM,
probs = 0.1,
type = 2)El cual es igual a 9.514. Trabajemos ahora con la muestra de la población creada en el ejemplo previo.
muestra %<>% mutate(TM = deaths / confirmed)
# Graficamos el histograma de TM para la muestra
muestra %>% ggplot(aes(TM)) +
geom_histogram(aes(y = ..count.. / sum(..count..)),
fill = 'slategray3',
color = "white",
bins = 25) +
labs(x = "Tasa de Mortalidad", y = "Proporción",
title = "Histograma para TM (Utilizando la muestra)") +
theme_minimal()
A partir de la muestra estimaremos el ratio:
p90m <- quantile(muestra$TM,
probs = 0.9,
type = 2)
p10m <- quantile(muestra$TM,
probs = 0.1,
type = 2)El valor estimado es igual a: 8.625.
A diferencia del ejemplo anterior, para el caso de este ratio es complicado obtener una expresión analítica para su varianza. Por tanto, para estimar dicha varianza utilizaremos Bootstrap:
rat90_10Boot <- c() # Creamos un vector vacío que iremos
# llenando en el siguiente for
for (n in 1:1000){
sampleboot <- sample(muestra$TM, # Extraemos una submuestra
size = 100, # aleatoria con reemplazo
replace = T)
rat90_10Boot<-c(rat90_10Boot, # n-1 ratios de submuestras
quantile(sampleboot, # calculamos el ratio para la
probs = 0.9, # submuestra n
type = 2) / quantile(sampleboot,
probs = 0.1,
type = 2))
}
rat90_10Boot <- data.frame(ratio = rat90_10Boot)Ahora crearemos un histograma para el ratio calculadas a partir de bootstrap:
rat90_10Boot %>% ggplot(aes(ratio)) +
geom_histogram(aes(y = ..count.. / sum(..count..)),
fill = "lightblue",
color = "white",
bins = 20) +
geom_vline(xintercept = p90p / p10p,
color = "red4") +
geom_label(mapping = aes(x = p90p / p10p,
y = .125,
label = paste("Ratio pob =",
round(p90p / p10p, 6)),
hjust = -.5, vjust = -.5),
colour = "red4") +
labs(x = "Ratio P_90/P_10 (Bootstrap)",
y = "Proporción",
title = "Histograma del ratio P_90/P_10") +
theme_minimal()
La línea vertical muestra el valor del ratio P_90/P_10 poblacional el cual es el número que tenemos como objetivo.
Al observar el histograma, podemos notar que la distribución del ratiomuestral no corresponde a la de una normal. Considerando esto último, calculemos intervalos de confianza al 95% utilizando las submuestras del Bootstrap a través de dos estrategias:
Asumamos que el ratio \(\varphi\) muestral se distribuye normal:
llr_1 <- mean(rat90_10Boot$ratio) - 1.96 * sd(rat90_10Boot$ratio)
ulr_1 <- mean(rat90_10Boot$ratio) + 1.96 * sd(rat90_10Boot$ratio)El intervalo de confianza al 95% es (llr_1, ulr_1) = (3.82824, 12.89114).
La segunda estrategia consiste en utilizar la distribución empírica que obtuvimos a través del bootstrap para el ratio \(\varphi\) muestral. Para esto calculamos el valor del percentil 2.5 y el del percentil 97.5. De esta forma obtenemos el intervalo de valores en los cuales se acumula el 95%:
llr_2 <- quantile(rat90_10Boot$ratio,
probs = .025,
type = 2)
ulr_2 <- quantile(rat90_10Boot$ratio,
probs = .975,
type = 2)El intervalo de confianza al 95% es (llr_2, ulr_2) = (5.30662, 14.26597).
Grafiquemos los intervalos de confianza anteriores para compararlos:
lr <- c(llr_1, llr_2)
ur <- c(ulr_1, ulr_2)
meanr <- c(mean(rat90_10Boot$ratio),
mean(rat90_10Boot$ratio))
tICr <- c("Normal", "Dist. Empírica")
ICr <- data.frame(tipo = tICr,
lower = lr,
media = meanr,
upper = ur)
ICr %>% ggplot(aes(x = media, y = tipo)) +
geom_point() +
geom_errorbarh(aes(xmin = lower,
xmax = upper),
size = 1,
height = 0.1) +
geom_vline(xintercept = p90p / p10p,
colour = "red4") +
labs(x = "Tasa de Letalidad",
y = "Método",
title = "Intervalos de Confianza al 95% para el ratio P_90/P_10") +
theme_minimal()
La línea vertical indica el valor del ratio \(\varphi\) poblacional.
Aquellos alumnos que sientan que necesitan un repaso más a detalle, se recomienda que revisen el capítulo 2 del Stock y Watson o el apéndice B del Wooldridge.↩︎
Para conocer mas acerca de muestreo se recomienda ver el video https://youtu.be/Rf-fIpB4D50↩︎
En clase haremos simulaciones utilizando la página http://faculty.carrollu.edu/ckuster/CT/Central\%20Limit\%20Theorem\%20Simulation.html para dar una mayor intuición acerca del tema.↩︎