library(ggplot2)
library(foreign)
library(kableExtra)
library(scales)
library(openxlsx)
library(reshape2)

PRELIMINARES

EVALUACIÓN

2 Exámenes parciales 30%
- 1o. cuando acabemos estimación puntual (2024-10-02).
- 2o. cuando acabemos pruebas de hipótesis (2024-11-27).
Tareas: 5%
Examen final: 35%
Formulario: permitido. El formulario oficial lo pueden descargar de la página del Departamento.

CRITERIOS Y POLÍTICAS

El canal primario de comunicación entre nosotros será la Bandeja de Entrada de Canvas. Por favor, estén atentos también a la sección de Anuncios.
La asistencia no es obligatoria, pero les recomiendo enfáticamente que asistan a clase. No hay un sustituto claro para el aprendizaje activo. Aprovéchenlo.
Como cualquier materia, se espera que ustedes hagan una gran parte del trabajo por su cuenta (fuera de clase). Una parte de este trabajo es previo a la clase, revisando el material propuesto para cada tema. Otra parte, también muy importante, es el trabajo posterior a la clase (repaso y tareas).
- Una buena regla “de dedo” es realizar, al menos, el mismo número de ejercicios que los realizados en clase y en las tareas de manera independiente.
- Hagan y entreguen las tareas. La experiencia señala claramente que quienes hacen el trabajo son los que obtienen resultados sobresalientes.
Sobre las calificaciones:
- Calificaciones parciales (exámenes y tareas) no se redondean. La calificación final sí.
- El redondeo es estricto, no hay negociaciones.
- Acorde al Art. 27 del Reglamento de Alumnos, es condición necesaria aprobar el examen final.
- La calificación del examen final no se redondea. Calificaciones menores a 6 (estricto) son reprobatorias.
En la calificación de los exámenes se evaluará (en la medida de lo posible): 1) planteamiento, 2) desarrollo y 3) resultado. El resultado correcto puede resultar en una valoración positiva; un resultado incorrecto, pero con buen planteamiento o desarrollo, puede aportar puntos; resultados incorrectos, sin planteamiento o desarrollo claros no permiten valoración alguna.
El estándar para la presentación de resultados será a 4 decimales.
Esta es una materia coordinada por el Departamento de Estadística. Algunas decisiones podrán estar subordinadas al criterio que el Departamento imponga.

BIBLIOGRAFÍA BÁSICA

Aguirre et al. (2006)
Wackerly, Mendenhall, and Scheaffer (2008)
DeGroot (1988)

INTRODUCCIÓN

ESTADÍSTICA

Lecturas recomendadas

Wackerly, Mendenhall, and Scheaffer (2008), caps. 1.1, 2.12.
Aguirre et al. (2006), caps. 11.3, 11.4.

Ejemplo: Conocer el ingreso de las personas de una determinada región es importante para una gran variedad de situaciones. Imagina que se desea conocer el ingreso de los habitantes de la Ciudad de México.

¿Qué harías para conocer este dato?

A partir del ejemplo de arriba muy probablemente surgieron un conjunto de preguntas. Usamos la estadística para dar respuesta a todas ellas. ¿Qué es la estadística entonces? En primer lugar, debemos reconocer que es un término ambiguo ya que puede referirse, por un lado, a una disciplina y, por otro, a una cosa.

Como cosa, solemos usar el término para referirnos a un número, un dato, generalmente resumen, que nos permite conocer algún aspecto de un objeto de estudio o interés.

Como disciplina, la estadística es una rama de las matemáticas (código 62, según el Mathematics Subject Classification System) y su campo de estudio reside, fundamentalmente, en el análisis de datos. Más específicamente, la estadística es la rama de las matemáticas que estudia los métodos para el acopio, análisis y presentación de datos (y, algunos agregan, para la toma de decisiones en presencia de incertidumbre). En general, el desarrollo de los métodos estadísticos ha sido motivado principalmente por la necesidad de encontrar patrones en los datos (Maindonald and Braun 2010).

Las herramientas que propociona la estadística como disciplina son fundamentales, entonces, para cualquier profesionista que se encarge de conseguir, analizar, interpretar o presentar información para la toma de decisiones.

Ejemplo: Considerando lo arriba planteado, ¿puedes pensar en actividades en las que la estadística no juegue algún rol?

Si bien como disciplina surge históricamente en el ámbito de la toma de decisiones gubernamentales, de estado (de ahí su nombre) en el S. XVI., la realidad es que el uso de información para la toma de decisiones es una actividad intrínsecamente humana y podemos encontrar ejemplos o casos de uso incluso en el funcionamiento del cerebro humano. No obstante, el tratamiento formal estadístico moderno comenzaron a tomar forma en los 1920s y 1930s (Maindonald and Braun 2010).

Ahora bien, como disciplina es frecuente que se divida, en atención al uso que se hará de la información:

descripción (estadística descriptiva)
inferencia¹ (inferencia estadística)

Sin embargo, antes de comenzar a revisar las técnicas estadísticas, vamos a hablar de algunos conceptos importantes para su aplicación.

POBLACIONES, MUESTRAS, PARÁMETROS Y ESTADÍSTICOS

Cuando realizamos análisis de datos generalmente nuestro objetivo es conocer alguna característica que nos interesa de una población. Es posible que la manera en la que usamos coloquialmente la palabra población no nos haga sentido en este contexto, pero en estadística una población se define como el conjunto de los objetos (individuos o eventos) de interés para un estudio.

Ejemplo: Los siguientes son diferentes ejemplos de posibles poblaciones de interés:

Las personas que habitaron en los Estados Unidos Mexicanos en los meses de enero a diciembre del 2023;
Las células de la médula ósea de una persona en particular;
Los billetes de $50 en circulación en un momento determinado;
Los dígitos de la expansión decimal de $\pi$;
Las personas mayores de 18 años que habitaron en los Estados Unidos Mexicanos en los meses de enero a diciembre del 2023;
Las personas mayores de 18 años del sexo femenino que habitaron en los Estados Unidos Mexicanos en los meses de enero a diciembre del 2023.

¿Puedes plantear preguntas de interés para cada una de estas poblaciones?

¿Qué información crees que se necesite acopiar para dar respuesta a tu pregunta?

Ejemplo: ¿Puedes intentar definir la población de nuestro ejemplo inicial?

Observa, de los ejemplos, lo siguiente:

La población puede ser un conjunto muy amplio o muy reducido, el tamaño no es lo que define a una población.
La definición de la población es fundamental, la manera en la que se defina a la población condiciona la información que se recabe y, por lo tanto, la respuesta que se obtenga.
No necesariamente podemos obtener información de toda (a lo que llamaríamos un censo) la población: porque la población no es finita o porque diferentes restricciones (e.g. tiempo, dinero, tipo de población) no permiten recabar información de toda la población.

Para los casos en los que no es posible o deseable recabar información de toda la población es necesario recurrir a muestras. Una muestra no es otra cosa más que un subconjunto la población.

En este sentido, podemos pensar que la utilidad de la muestra consiste en que podemos observar en la muestra la característica de interés y de alguna manera utilizar esta observación para inferir el valor de dicha característica en la población. Utilizamos el nombre de parámetro para referirnos a ese valor desconocido de la característica de interés de la población, y nos referimos como estadística al valor que dicha característica toma en la muestra observada.

Ejemplo: La Encuesta Nacional de Ingresos y Gastos de los Hogares (ENIGH) tiene como objetivo proporcionar un panorama estadístico del comportamiento de los ingresos y gastos de los hogares en cuanto a su monto, procedencia y distribución. Adicionalmente ofrece información sobre las características ocupacionales y sociodemográficas de los integrantes del hogar, así como las características de la infraestructura de la vivienda y el equipamiento del hogar.

¿Cuál es la población objetivo?¿Cuáles sus parámetros?

¿Te imaginas algunas de las estadísticas que se construyen con esta encuesta?

¿Qué uso crees que se haga de esta información?

Ahora bien, una muestra puede ser obtenida o determinada de muchas maneras, pero en general vamos a distinguir entre muestras deterministas y muestras aleatorias. Una muestra determinista es determinada utilizando un criterio discrecional, mientras que una muestra aleatoria es determinada utilizando un criterio probabilístico.

En primera instancia, y para alguien que no haya reflexionado mucho al respecto, podría quizás resultar atractivo el proponer un levantamiento de muestras no aleatorias. Es decir, alguien que pudiera tener cierto conocimiento del dominio de interés podría sugerirse como mejor calificado para seleccionar a los individuos de una población que mejor representen a las características de la población, sobre todo mejor que un procedimiento aleatorio.

Sin embargo, poco tardaremos en darnos cuenta de que esto en realidad no es así. Las muestras no aleatorias tienen precisamente el problema de que no hay manera de garantizar que los individuos seleccionados en la muestra sean individuos que representen razonablemente bien a la generalidad de la población en la característica que nos interesa estudiar.

El muestreo aleatorio busca controlar no a qué individuos son incluidos en la muestra, sino la probabilidad de que un individuo perteneciente a la población sea incluido en la muestra. El hecho de conocer la probabilidad de inclusión permite entonces, paradójicamente, garantizar que la muestra recogida es representativa de la población, al mismo tiempo que nos permite desarrollar las técnicas necesarias para poder medir la incertidumbre asociada al hecho de que la muestra no es, en efecto, la población. De hecho, hay quien asegura que el objetivo de la estadística es fundamentalmente desarrollar las herramientas matemáticas para realizar inferencias sobre las características de la población a partir de la información obtenida de una muestra aleatoria.

Ahora, el levantamiento de una muestra aleatoria puede ser realizado de diferentes maneras (siguiendo diferentes procedimientos). Desde luego, la manera en la que se levanta la muestra tiene implicaciones pues la elección del método de muestreo determina la cantidad de información contenida en la muestra, así como la probabilidad de observar algún resultado en específico.

El diseño de experimentos es, por ello, de suma importancia, aún cuando la atención se centra, con frecuencia, en el proceso de inferencia.

El método más sencillo de levantamiento de muestras es conocido como muestreo aleatorio simple (MAS).

MUESTREO ALEATORIO SIMPLE

Definición: Sean $N$ y $n$ el número de elementos de una población y de una muestra de la población, respectivamente. Se le llama muestreo aleatorio simple al procedimiento de selección de muestras en el que cada una de las ${N \choose n}$ posibles muestras que se pueden seleccionar de la población tienen la misma probabilidad de ser seleccionadas.

Lo anterior, desde luego, implica que cada individuo de la población tiene la misma probabilidad de ser seleccionado para la muestra (sin reemplazo), de ahí que se pueda garantizar la representatividad de la muestra en igualdad de circunstancias de los individuos de la población cuando el tamaño de la muestra es el apropiado. El problema con el MAS es que cuando la muestra es pequeña, por pura probabilidad, podríamos estar seleccionando individuos que no representen satisfactoriamente a la generalidad de la población. Por ello, es de nuestro interés asegurar que la muestra levantada sea una muestra representativa. Cuando el muestreo es MAS, esta muestra representativa es posible que requiera ser de un tamaño considerable. Otras técnicas de muestreo han sido creadas, entonces, para lograr una mayor representatividad de la muestra, con un menor tamaño.

También puede ser el caso que interés del estudio se centre en diferentes parámetros, o bien en parámetros para sub-poblaciones de interés. En estos casos, nuevamente, es muy probable que el MAS no nos proporcione la información requerida de manera apropiada. Nuevamente, para ello es necesario requerir a métodos más complejos de muestreo.

El MAS es, por lo tanto, apropiado principalmente cuando la población bajo estudio es relativamente homogénea o bien, no contamos con información previa que nos permita asumir que es necesario llevar a cabo un esquema de muestreo más complejo.

Ejemplo: Supongamos que se desea conocer el número de eritrocitos en la sangre de una persona. Dada la naturaleza del problema, claramente no sería apropiado extraer toda la sangre de la persona para realizar la medición de la cantidad de eritrocitos (asumiendo que es requerido que la persona sobreviva al procedimiento). Sin embargo, dada la homogeneidad de la sangre en todo el cuerpo, es posible realizar el procedimiento con una muestra significativamente pequeña de sangre. Dado que no controlamos exactamente qué mililitros de sangre extraer, el método de extracción constituye una muestra aleatoria simple.

Ejemplo: Para nuestro problema de los ingresos, ¿crees que sea válido o conveniente realizar un MAS para conocer el ingreso de los habitantes de la Ciudad de México?

En la práctica, es posible que no siempre contemos con un mecanismo natural de muestreo aleatorio (como en el ejemplo de la sangre). En estos casos, tenemos el problema de determinar en forma aleatoria a los miembros de la población que formarán parte de la muestra. Para ello es necesario contar con un marco muestral.

Marco muestral: es el registro (lista, algoritmo o procedimiento) que identifica a todos los miembros de la población susceptibles de ser seleccionados en una muestra.

Toda vez que contamos con el marco muestral, el procedimiento puede volverse relativamente sencillo ya que solamente necesitamos generar valores aleatorios de un índice que recorra todo el marco muestral. Esto, hoy en día, es muy sencillo.

Ejemplo: Obtén aleatoriamente el identificador de 50 individuos de entre una población de 10,000 individuos (sin reemplazo).
R: No contamos con un marco muestral propiamente, sin embargo sabemos que se trata de 10,000 individuos que podemos numerar consecutivamente del 1 al 10,000, entonces queremos seleccionar aleatoriamente 50 de estos números sin reemplazo.

sample(x = 1:10000, size = 50, replace = FALSE)

##  [1] 6519  158 9508  309 7155 8660 8763 9016 9156 3556 9075 9251 8354 3444   69
## [16]  398 6446 3772 5852 1588 2239 9004  483 4860 3272 5434 5780 5338 1058 8325
## [31] 4803 3870 4699 1732 1526 3164  762 9647 5589 7561 7666  397 7697 4772 2263
## [46] 3413 9888 1946 4243 9576

ESTADÍSTICA DESCRIPTIVA

Lectura recomendada

Aguirre et al. (2006), cap. 1.

TIPOS DE VARIABLES Y ESCALAS DE MEDICIÓN

Para poder “explorar” los datos de interés, primero tenemos que conocer con qué tipo de datos contamos. Algunas técnicas / herramientas serán apropiadas para un tipo de datos pero no para otros.

Primero, consideremos la naturaleza de los datos, pudiendo tratarse de datos provenientes de documentos textuales (datos no estructurados) o bien de datos con algún tipo de formato estructurado (muy frecuentemente en formato tabular).

Por brevedad, centraremos nuestro interés en datos estructurados y asumiremos, en general, que adoptan (o pueden ser llevados a) una estructura tabular.

Dentro de los datos estructurados, entonces, contamos con una estructura de datos cuyos atributos (columnas) representan variables para el análisis. Estas variables o atributos pueden ser de diferentes tipos:

Cualitativos o categóricos: un dato categórico es aquel que consiste en “etiquetas” que usamos para denotar pertenencia ciertos grupos, clases o categorías. Estas etiquetas pueden o no tener un significado literal, pueden o no tener una relación entre sí y, desde luego, pueden o no estar representadas por un número. Idealmente, las categorías deben ser exhaustivas y mutuamente excluyentes (aunque, en la práctica, con frecuencia estos requisitos suelen no presentarse).
- Nominales: son datos categóricos cuyas etiquetas asignan “nombres” a los grupos a los que se encuentran asociadas. Pueden o no ser numéricas, pero no es posible establecer una relación entre las etiquetas más allá de la pertenencia o no al grupo. Por ejemplo, $\{Rojo, Negro\}$ puede ser un conjunto de valores categóricos nominales, al igual que $\{Grupo 1, Grupo2\}$ o bien $\{1,2\}$, cuando el 1 y el 2 los usamos para representar al grupo y no necesariamente una relación numérica en los datos.
- Ordinales: a diferencia de los datos categóricos nominales, los datos categóricos ordinales si nos permiten establecer un “orden” en los datos. Por ejemplo, si nuestra variable adopta los valores $\{alto, bajo\}$, podemos en general asumir que $bajo < alto$. Observa, sin embargo, que por lo general este tipo de valores no me permiten evaluar al distancia entre los valores, es decir, no me permite determinar cuál es exactamente la diferencia que existe entre $alto$ y $bajo$, simplemente que existe una relación de orden. En este mismo sentido, podríamos habernos encontrado con valores del tipo $\{1,10\}$, donde $1 = bajo; 10 = alto$. A pesar de haber usado valores numéricos para representar a las categorías, de todos modos no sería válido afirmar que existe una distancia de 9 entre el valor $alto$ y el valor $bajo$, ya que se trata de una variable cualitativa ordinal.
Cuantitativos: son variables cuya naturaleza es numérica y representan una medición de un atributo de una cosa.
- Discretos: son variables que únicamente pueden tomar valores enteros. Generalmente registran el resultado de un conteo. Por ejemplo: el número de incidentes observados durante un intervalo de tiempo determinado.
- Continuos: son variables que pueden tomar valores reales. Generalmente registran el resultado de una medición. Por ejemplo: la altura de las personas pertencientes a una población de interés. Con frecuencia damos tratamiento de variables continuas a datos que estrictamente quizá no lo son, pero que por cuestiones prácticas es conveniente darles dicho tratamiento (por ejemplo, el precio de algunos productos en la economía real).

Para las variables cuantitativas, además, podemos tener diferentes escalas de medición:

De intervalos: es una escala cuantitativa en la que, cuando comparamos dos mediciones, lo relevante es la diferencia entre esas dos mediciones (son escalas aditivas). Por este motivo, las escalas de intervalos pueden incluir entre sus valores al 0 y valores negativos, sin embargo el significado del cero es arbitrario. Dado que el cero es un valor arbitrario, las comparaciones proporcionales entre mediciones no son generalmente apropiadas o significativas.

La temperatura en grados Celsius (o Farenheit) es un ejemplo de escala de medición de intervalos. Es una medida de la cantidad de calor presente en un ambiente, pero es posible hablar de grados Celsius negativos de temperatura porque el cero es un punto arbitrario de referencia (el punto de congelación del agua).

De razones: en las escalas de razones, el 0 tiene un significado específico: la total ausencia de lo que queremos medir. Por este motivo, las escalas de razones no permiten la presencia de valores negativos y, también por este motivo, es posible caracterizarlas como escalas multiplicativas, por lo que las comparaciones en términos de proporciones aquí sí hacen sentido.

¿Qué ejemplos de escalas de medición de razones encuentras?

Ejemplo: Abajo se muestra la serie de tiempo de los precios del NAFTRAC a partir del 1 de enero del 2020. ¿Qué tipo de dato es?¿En qué escala se encuentra?

temp <- read.xlsx(xlsxFile = 'datos/NAFTRACISHRS.xlsx', detectDates = TRUE)

ggplot(data = temp) +
  aes(x = Date, y = Close) +
  geom_line() +
  theme_classic()

ggplot(data = temp) +
  aes(x = Close) +
  geom_density() +
  theme_classic()

ggplot(data = temp) +
  aes(x = log(Close)) +
  geom_density() +
  theme_classic()

ggplot(data = temp) +
  aes(x = Date, y = Yield) +
  geom_line() +
  geom_line(mapping = aes(y = Yield_2), colour = 'red') +
  theme_classic()

ggplot(data = temp) +
  aes(x = Yield) +
  geom_density() +
  theme_classic()

Cuando estudiamos el tipo de datos con los que contamos, también es importante considerar la manera en la que se nos presentan:

Univariados: observamos una única variable de interés.
Multivariados: observamos para cada registro (renglón) de nuestra tabla múltiples características (atributos) de manera simultánea.
- Datos composicionales: los datos forman parte de un todo.
- Series de tiempo: las observaciones están referidas (y pueden ser ordenadas) a una referencia de tiempo (segundos, minutos, horas, días, semanas, meses, años, …). Datos que no capturamos u observamos como series de tiempo son conocidos como datos de corte transversal (podría pensarse que, por lo tanto, serían datos acopiados en un único momento en el tiempo, pero puede ser también que simplemente la variable relativa a la temporalidad no se consideró relevante).
- Datos en/de panel: son datos que combinan características temporales y transversales. Típicamente se recogen estos datos para dar seguimiento a una población determinada (p.e., un conjunto bien definido de pacientes de un hospital sobre los que quiere conocerse su evolución en el tiempo), aunque pueden existir conjuntos de datos en panel para los cuales la definición de la población puede ser más compleja.

AED UNIVARIADO

Queremos describir el conjunto de datos con los que estamos trabajando, y comenzaremos describiendo las variables una a una.

La manera más apropiada de describir una variable en particular dependerá, desde luego, del tipo de datos que tenemos, pero en general buscamos caracterizar a la variable de interés describiendo, de alguna manera, la variabilidad observada en, o bien, mediante estadísticas, es decir, mediante un nuevo dato, “sintético”, creado a partir de nuestros datos, que describe en forma resumida algún rasgo relevante de la variable.

En términos generales, algunos de los elementos en los que con mayor frecuencia vamos a centrar nuestra atención son:

la forma de la distribución (¿Al rededor de qué valores se concentran más nuestros datos?¿En dónde menos?);
observaciones atípicas (¿Qué significa que algo sea atípico?¿Por qué es importante?);
huecos en los datos (¿Es esperado que haya huecos?¿Qué significa que haya huecos? ¿Qué implica para nuestro análisis?).

¿Por qué hacemos esto? Pues porque nos interesa conocer la forma en la que se “distribuyen” nuestros datos. Esto puede tener muchas implicaciones muy diversas:

Puede determinar los pasos a seguir en mi proyecto (p.e., identificar variables que necesitan ser transformadas)
Puede determinar el tipo de herramientas/modelos a utilizar (p.e., identificando si es razonable esperar que algunos de los supuestos se cumplan).
Permite identificar posibles relaciones entre las variables que determine o impacte en algunos de los puntos ya señalados.
Permite identificar posibles problemas en los datos.

DISTRIBUCIONES DE FRECUENCIA

Una primera forma básica de describir el comportamiento de una variable consiste simplemente en contar las veces en las que ocurre cada uno de los valores que puede tomar dicha variable. A esto le llamamos la distribución de frecuencia de la variable.

Variables cualitativas: La manera de representar la distribución de frecuencia de una variable cualitativa es utilizando una tabla de frecuencias. En una tabla de frecuencias se registra, para cada uno de los posibles valores de una variable cualitativa (su dominio) el número de veces que cada uno de esos valores se observa en los datos ($f_i$). Usualmente, se registra en la tabla de frecuencias también la frecuencia relativa ($p_i$) de cada posible valor de la variable, esto es, la frecuencia del valor correspondiente divida entre el número total de observaciones.
Ejemplo: En la ENIGH 2022 se encuestó a 105,525 viviendas, representando 106,893 hogares. La pregunta 1 de la encuesta recoge la clase de vivienda.

Clase	Descripción
1	Casa independiente
2	Departamento en edificio
3	Vivienda en vecindad
4	Vivienda en cuarto de azotea
5	Local no construido para habitación

¿Qué tipo de variable es esta? Construye su tabla de frecuencias.

viviendas.enigh <- read.dbf(file = 'datos/viviendas.dbf')

temp <- aggregate(x = list(Frecuencia = viviendas.enigh$tipo_viv), by = list(Clase = viviendas.enigh$tipo_viv), FUN = length)

temp$frec.rel <- temp$Frecuencia / sum(temp$Frecuencia)

kable(x = temp |> transform(Frecuencia = comma(Frecuencia, accuracy = 1), frec.rel = percent(frec.rel, accuracy = 0.01)), col.names = c('Clase','Frecuencia','Frec. Rel.'), align = c('l','r','r'), format = 'html') |> kable_classic()

Clase	Frecuencia	Frec. Rel.
&	107	0.12%
1	84,259	94.86%
2	2,141	2.41%
3	1,442	1.62%
4	103	0.12%
5	771	0.87%

Desde luego, a las tablas de frecuencias para datos categóricos es posible representarlas gráficamente mediante gráficas de barras.

Grafica los resultados anteriores en una gráfica de barras.

ggplot(data = viviendas.enigh) +
  aes(x = tipo_viv) +
  geom_bar() +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Número de viviendas por clase de vivienda', x = 'Clase', y = 'Número de viviendas') +
  scale_y_continuous(labels = comma)

ggplot(data = viviendas.enigh) +
  aes(x = tipo_viv) +
  geom_bar(aes(y = (after_stat(count))/sum(after_stat(count)))) +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Porcentaje de viviendas por clase de vivienda', x = 'Clase', y = 'Número de viviendas') +
  scale_y_continuous(labels = percent)

Variables cuantitativas: Para la descripción de este tipo de variables mediante una distribución de frecuencias es necesario considerar, en primer lugar, el tipo de variable cuantitativa:

Variables cuantitativas discretas: este caso es relativamente sencillo ya que podemos aplicar exactamente el mismo procedimiento que utilizamos para las variables cualitativas. Esto es, para este tipo de variables simplemente podemos registrar en forma tabular tanto la frecuencia como la frecuencia relativa observadas para cada valor posible de la variable.
Variables cuantitativas continuas: si bien para las variables continuas sería posible aplicar nuevamente el mismo procedimiento, la realidad es que muy probablemente registraríamos una sola observación para cada valor, lo cual no sería de mucha utilidad si nuestro objetivo es describir el comportamiento de la variable. Por ello, es necesario agrupar a las observaciones en rangos de valores. Surge entonces, necesariamente, la pregunta, ¿cuántos rangos tenemos que considerar para agrupar las observaciones?

Si bien no es posible dar una respuesta directa a esta pregunta, existen una multiplicidad de propuestas para poder determinar el número “ideal” de rangos a utilizar. Una de estas propuestas es la regla de Freedman-Diaconis que sugiere utilizar intervalos de longitud igual a $h = \frac{2 \times IQR}{n^{1/3}}$ donde $IQR$ es el rango intercuartil y $n$ es el número de observaciones.

A su vez, $IQR = Q3 - Q1$ y los valores $Q3$ y $Q1$ corresponden al tercer y primer cuartil de los datos, respectivamente (es decir, los valores que son mayores al 75% y 25% de las observaciones, respectivamente).

Otra alternativa, por ejemplo, es la regla de Sturges, que sugiere usar $1 + \log_2 n$ intervalos.

Una vez que hemos determinado la longitud de los rangos a utilizar obtener el número de rangos necesarios para la construcción de la distribución de frecuencias se sigue directamente de considerar el rango de los datos de nuestra variable. Habrá que tener solamente cuidado en considerar los límites inferior y superior apropiados para establecer el rango ya que en ocasiones estos límites estarán determinados por los valores observados y en ocasiones por los valores posibles.

Ahora bien, es importante recordar que esta es una de entre una multiplicidad de propuestas y que, sin importar la regla que se use, se recomienda no utilizar menos de cinco y no más de veinte rangos para la construcción de las tablas de frecuencias.

Por último, es importante observar que el número de rangos que escojamos puede cambiar significativamente la manera en la que se presentan los resultados. No hay recetas mágicas, es necesario probar con diferentes opciones que logren comunicar apropiadamente el contenido de nuestra información.

¿Los rangos son útiles únicamente para las variables continuas?

La pregunta 23 de la ENIGH recaba información sobre el valor de la renta en las viviendas. Construye la tabla de frecuencias para la renta registrada en la encuesta.

R: Primero determinaremos el número de rangos a utilizar. Contamos con 88,823 observaciones ($n$).

n <- nrow(viviendas.enigh)

Para obtener datos completos, vamos a juntar en un único vector los datos contenidos en los reactivos 23.1 y 23.2:

renta <- viviendas.enigh$renta
renta[is.na(renta)] <- viviendas.enigh$estim_pago[is.na(viviendas.enigh$renta)]

Obtenemos ahora el rango intercuartil:

ric <- IQR(renta)

ric

## [1] 1900

Calculamos la longitud de los intervalos:

h <- (2*ric)/(n^(1/3))
h

## [1] 85.16774

Observa que esto nos dará: 2353.0036438 intervalos. No parece una muy buena opción. Utilizaremos entonces la regla de Sturges.

Necesitamos entonces:

k <- 1 + log(n,2)
k

## [1] 17.43865

intervalos.

Creamos la variable que captura el rango:

rango <- rep(x = NA, times = n)

for (i in 1:ceiling(k)){
  rango[which(renta >= min(renta) + (i-1)*(max(renta) - min(renta))/ceiling(k) & renta < min(renta) + i*(max(renta) - min(renta))/ceiling(k))] <- i
}

y tabulamos:

rangos.rentas <- aggregate(x = list(Frecuencia = rango), by = list(Rango = rango), FUN = length)

rangos.rentas$frec.rel <- rangos.rentas$Frecuencia/sum(rangos.rentas$Frecuencia)

kable(x = rangos.rentas |> transform(Frecuencia = comma(Frecuencia, accuracy = 1), frec.rel = percent(frec.rel, accuracy = 0.01)), col.names = c('Rango','Frecuencia','Frec. Rel.'), align = c('l','r','r'), format = 'html') |> kable_classic()

Rango	Frecuencia	Frec. Rel.
1	87,378	98.37%
2	1,208	1.36%
3	177	0.20%
4	38	0.04%
5	14	0.02%
7	3	0.00%
9	2	0.00%
11	1	0.00%
17	1	0.00%

A la representación gráfica, mediante gráficas de barras, de la distribución de frecuencias de una variable cuantitativa normalmente se le conoce como un histograma. Los histogramas son usados entonces para describir de manera aproximada la densidad.

Por otro lado, es posible también graficar la frecuencia acumulada de los valores de la variable de interés, en cuyo caso estaríamos describiendo de manera aproximada la función de distribución de la variable. Este tipo de gráficas se conocen como ojivas.

Obtén el histograma de las rentas, su densidad y su ojiva.

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = Renta) +
  geom_histogram(bins = k, boundary = 1) +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Histograma del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Frecuencia') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

## Warning: Computation failed in `stat_bin()`.
## Caused by error in `bin_breaks_bins()`:
## ! `bins` must be a whole number, not the number 17.44.

hist(renta)

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = Renta) +
  geom_histogram(bins = 5, boundary = 1) +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Histograma del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Frecuencia', caption = 'Cinco intervalos.') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = Renta) +
  geom_histogram(bins = 100, boundary = 1) +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Histograma del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Frecuencia', caption = 'Cien intervalos.') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

ggplot(data = data.frame(Renta = renta[which(renta <= 20000)])) + 
  aes(x = Renta) +
  geom_histogram(bins = 1 + log(length(renta[which(renta <= 20000)]),2), boundary = 1) +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Histograma del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Frecuencia', caption = 'Rentas menores a $20,000.') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

## Warning: Computation failed in `stat_bin()`.
## Caused by error in `bin_breaks_bins()`:
## ! `bins` must be a whole number, not the number 17.43.

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = Renta) +
  geom_density() +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Densidad empírica del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Densidad') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = log(Renta)) +
  geom_density() +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Densidad empírica del logaritmo del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Densidad') +
  scale_y_continuous(labels = comma)

Para graficar la ojiva, primero calcularemos la función de distribución empírica:

ecdf.renta <- ecdf(x = renta)

Ahora creamos un data frame con los datos de las rentas y su probabilidad acumulada:

temp <- data.frame(Renta = renta, P = ecdf.renta(renta))

Y la graficamos:

ggplot(data = temp) +
  aes(x = Renta, y = P) +
  geom_point() +
  geom_line() +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Distribución empírica del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Probabilidad acumulada') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

plot(ecdf(renta))

Ahora, si ya contamos con una descripción gráfica de la distribución de nuestros datos, típicamente buscamos fijarnos en ciertos rasgos de dicha distribución que nos ayuden a describir su forma. Normalmente, los rasgos de interés son:

Simetría: es decir, si la forma de la densidad de nuestros datos es simétrica o asimétrica. En el caso de distribuciones asimétricas, típicamente las describiremos en función de su sesgo. En estos casos podemos tener, entonces, densidades sesgadas a la izquierda o a la derecha.

temp <- data.frame(Valores = rnorm(n = 10000))

ggplot(data = temp) +
  aes(x = Valores) +
  geom_density() +
  theme_classic() +
  labs(y = NULL, caption = 'Función de densidad simétrica.')

temp <- data.frame(Valores = rgamma(n = 10000, shape = 2))

ggplot(data = temp) +
  aes(x = Valores) +
  geom_density() +
  theme_classic() +
  labs(y = NULL, caption = 'Función de densidad sesgada a la derecha.')

temp <- data.frame(Valores = rbeta(n = 10000, shape = 10, shape2 = 0.75))

ggplot(data = temp) +
  aes(x = Valores) +
  geom_density() +
  theme_classic() +
  labs(y = NULL, caption = 'Función de densidad sesgada a la izquierda.')

Curtosis: la curtosis se refiere a qué tan concentrados se encuentran los valores de la distribución al rededor de un valor central. La curtosis normalmente se evalúa en relación a una distribución de referencia, la distribución normal. En este sentido, decimos que una distribución que muestra una concentración mayor que la de la distribución normal es una distribución leptocúrtica; si su concentración es menor, platicúrtica; si su concentración es igual, mesocúrtica.

temp <- data.frame(Valores = rnorm(n = 10000))
temp2 <- data.frame(Valores = rnorm(n = 10000, sd = 2))
temp3 <- data.frame(Valores = rnorm(n = 10000, sd = 1/2))
ggplot(data = temp) +
  aes(x = Valores) +
  geom_density() +
  geom_density(data = temp2, colour = 'red') +
  geom_density(data = temp3, colour = 'blue') +
  theme_classic() +
  labs(y = NULL, caption = 'Negro: función de densidad mesocúrtica\n Rojo: función de densidad platicúrtica\n Azul: función de densidad leptocúrtica')

Moda: la moda de la densidad se define como el valor con la mayor frecuencia/densidad en los datos. No obstante, en ocasiones podemos observar distribuciones multimodales (i.e., con más de una moda). En estos casos, con frecuencia, se trata de modas locales y se aprecian como funciones de densidad con dos (o más) jorobas.

temp <- data.frame(Valores = c(rnorm(n = 5000), rnorm(n = 5000, mean = 5)))

ggplot(data = temp) +
  aes(x = Valores) +
  geom_density() +
  theme_classic() +
  labs(y = NULL, caption = 'Función de densidad bi-modal.')

MEDIDAS DESCRIPTIVAS

Si bien el análisis gráfico de la distribución de los datos es con frecuencia útil y eficiente en la comunicación de los rasgos generales de los datos, el análisis gráfico puede que no transmita con mucha precisión las características de nuestros datos. Para ello, entonces, necesitamos de mediciones.

Para describir a nuestros datos mediante medidas, entonces, vamos a centrarnos en las siguientes características de los datos:

Posición
Dispersión
Sesgo
Curtosis

MEDIDAS DE POSICIÓN

Nos referimos aquí como medidas de posición a medidas que nos indican la ubicación de puntos específicos dentro de la distribución de nuestros datos. Así, las medidas de posición más frecuentemente utilizadas son:

Mínimo / máximo: son, respectivamente, el menor y el mayor de los valores observados en nuestros datos. Cuando analizamos los datos observados, es muy importante tomar en consideración los valores extremos teóricos de nuestros datos ya que la comparación de los valores observados contra los valores teóricos nos pueden ayudar a identificar problemas en la calidad de nuestros datos (e.g., valores no permitidos) o bien problemas en el acopio de los datos (e.g., valores muy alejados del mínimo teórico).
Moda: como ya se mencionó la moda corresponde al valor observado con mayor frecuencia en nuestros datos. En datos continuos es necesario considerar la densidad de la distribución, no la frecuencia observada. Para el análisis de la moda, como ya se mencionó, en ocasiones es importante preguntarse si existen modas locales (para lo cual el análisis gráfico puede ser una herramienta auxiliar muy útil).
Media: La media no es otra cosa que el promedio aritmético de las observaciones (ya sea de la población o de la muestra). De las medidas de posición, la media es quizás la más conocida y utilizada. Esto refleja, por ejemplo, en el hecho de que incluso típicamente se reserva la letra $\mu$ para representar a la media poblacional; y el símbolo $\bar{X}$ para la media muestral (suponiendo que $X$ representa a las observaciones de la muestra).

\[ \mu = \frac{\sum\limits_{i=1}^N X_i}{N};\\ \bar{X} = \frac{\sum\limits_{i=1}^n X_i}{n} \]

donde $N$ es el tamaño de la población y $n$ el tamaño de la muestra.

Es importante señalar que la media se puede calcular aún cuando no tenemos las observaciones individuales sino agregadas mediante una tabla de frecuencias. En este caso, si cada $X_i$ representa el valor observado en la variable de interés y $f_i$ la frecuencia con la que fue observado dicho valor, la media será calculada como:

\[ \mu = \frac{\sum\limits_{i=1}^k f_i \times X_i}{N};\\ \bar{X} = \frac{\sum\limits_{i=1}^k f_i \times X_i}{n} \]

donde $N = \sum\limits_{i=1}^k f_i$ o $n = \sum\limits_{i=1}^k f_i$ según se trate de la población o de una muestra, y $k$ es el número de diferentes valores de $X$ observados.

Percentiles: un percentil, el p-ésimo percentil, es aquél valor de los datos tal que el p% de las observaciones son menores a dicho valor. Para calcular cualquier percentil, entonces, es necesario, en primer lugar, ordenar los datos. Una vez ordenados los datos, calculamos el índice $i = n \times \frac{p}{100}$. $i$ puede o no ser entero. Si $i$ es entero, el valor del percentil corresponderá al promedio de los valores en las posiciones $i$ e $i+1$; si no lo es, el cuantil será el valor de la observación en la posición correspondiente al entero siguiente a $i$.
- Mediana: corresponde al percentil 50. Es el valor que, una vez que ordenamos a los datos de menor a mayor, se encuentra justo a la mitad de las observaciones. Así, exactamente el 50% de los datos serán menores o iguales a la mediana; el otro 50% mayores o iguales.

Ejemplo: Calcula el máximo, mínimo, moda, media, primer, segundo y tercer cuartil de las rentas registradas en la ENIGH.

Llamamos cuartiles a los cuantiles que dividen a las observaciones en cuartos.

R

Máximo:

max(renta)

## [1] 200450

Mínimo:

min(renta)

## [1] 50

Moda: veamos una tabla con la frecuencia observada para los 3 valores más frecuentes en el vector de las rentas:

temp <- aggregate(x = list(Frecuencia = renta), by = list(Renta = renta), FUN = length)

kable(x = temp[order(temp$Frecuencia, decreasing = TRUE), ][1:3, ] |> transform(Renta = dollar(Renta), Frecuencia = comma(Frecuencia, accuracy = 1)), row.names = FALSE, align = c('r','r')) |> kable_classic()

Renta	Frecuencia
$2,000	13,009
$1,500	12,268
$1,000	11,292

La moda se encuentra, entonces, en $2,000.

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = Renta) +
  geom_density() +
  geom_vline(xintercept = 2000, colour = 'red') +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Densidad empírica del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Densidad') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

Media: primero veamos que tenemos 88,823 observaciones, entonces, la media aritmética de las rentas registradas en la muestra será:

sum(renta)/length(renta)

## [1] 2631.156

También podemos obtenerla de la tabla de frecuencias que construimos para el cálculo de la moda:

sum(temp$Renta * temp$Frecuencia)/sum(temp$Frecuencia)

## [1] 2631.156

nota que:

sum(temp$Frecuencia)

## [1] 88823

Finalmente, podríamos habérsela pedido directamente a R:

mean(renta)

## [1] 2631.156

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = Renta) +
  geom_density() +
  geom_vline(xintercept = 2000, colour = 'red') +
  geom_vline(xintercept = mean(renta), colour = 'green') +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Densidad empírica del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Densidad') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = log(Renta)) +
  geom_density() +
  geom_vline(xintercept = log(2000), colour = 'red') +
  geom_vline(xintercept = log(mean(renta)), colour = 'green') +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Densidad empírica del valor del logaritmo la renta de la vivienda', x = 'Valor de la renta (log)', y = 'Densidad') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

Cuantiles: calculemos la posición de los cuantiles solicitados:

i.q1 <- length(renta) * 25/100
i.q2 <- length(renta) * 50/100
i.q3 <- length(renta) * 75/100

$i(Q1)$: 22,205.75
$i(Q2)$: 44,411.50
$i(Q3)$: 66,617.25

Dado que las posiciones corresponden a números fraccionales podemos recabar directamente el valor de los cuantiles correspondientes:

Q.1 <- renta[order(renta)][ceiling(i.q1)]
Q.2 <- renta[order(renta)][ceiling(i.q2)]
Q.3 <- renta[order(renta)][ceiling(i.q3)]

$Q1$: $1,100
$Q2$: $2,000
$Q3$: $3,000

En R, podemos obtener rápidamente estos valores, y otros, mediante la función quantile:

quantile(renta)

##     0%    25%    50%    75%   100% 
##     50   1100   2000   3000 200450

Observa, finalmente, que el valor del 2do. cuartil, el percentil 50, corresponde con el valor de la mediana:

median(renta)

## [1] 2000

Entonces, en el 2022, el 50% de los hogares contestaron pagar una renta menor o igual a $2,000.

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = Renta) +
  geom_density() +
  geom_vline(xintercept = 2000, colour = 'red') +
  geom_vline(xintercept = mean(renta), colour = 'green') +
  geom_vline(xintercept = median(renta), colour = 'blue') +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Densidad empírica del valor de la renta de la vivienda', x = 'Valor de la renta', y = 'Densidad') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

ggplot(data = data.frame(Renta = renta)) + 
  aes(x = log(Renta)) +
  geom_density() +
  geom_vline(xintercept = log(2000), colour = 'red') +
  geom_vline(xintercept = log(mean(renta)), colour = 'green') +
  geom_vline(xintercept = log(median(renta)), colour = 'blue') +
  theme_classic() +
  labs(title = 'ENIGH 2022', subtitle = 'Densidad empírica del valor del logaritmo la renta de la vivienda', x = 'Valor de la renta (log)', y = 'Densidad') +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = comma)

La moda, la media y la mediana son conocidas como medidas de centralidad o de tendencia central porque caracterizan la parte central de la distribución de los datos. Adicionalmente, la moda y la mediana son consideradas medidas robustas de centralidad, ya que no se ven afectadas por la presencia de valores extremos.

Ejemplo (cont.): Imagina que nos dicen que hubo un error de captura en las rentas de la ENIGH y que la renta máxima registrada, 200450, en realidad era una renta de $20,450. ¿Cómo afecta esto a los cálculos realizados anteriormente?

renta.2 <- renta
renta.2[which(renta.2 == max(renta.2))] <- 20450

max(renta.2)

## [1] 180000

min(renta.2)

## [1] 50

names(table(renta.2))[which(table(renta.2) == max(table(renta.2)))]

## [1] "2000"

mean(renta.2)

## [1] 2629.129

quantile(renta.2)

##     0%    25%    50%    75%   100% 
##     50   1100   2000   3000 180000

¿Qué crees que tiene que pasar para que se modifique la mediana o la moda?

MEDIDAS DE DISPERSIÓN

Con lo que hemos visto hasta ahora podemos hacernos una idea de dónde se encuentra ubicados algunos puntos relevantes de la distribución de nuestros datos. Otra característica que generalmente es de interés es la dispersión, es decir, qué tan concentrados o dispersos se encuentran nuestros datos respecto de ciertos puntos de referencia.

Amplitud o rango: nos dice la distancia que existe entre el máximo y el mínimo de los datos registrados.

\[ R = Max(X) - Min(X) \]

Rango intercuartílico: ya lo utilizamos cuando calculamos el valor de Freedman-Diaconis, consiste en la distancia entre el valor correspondiente al primer cuartil y el valor del tercer cuartil. Por lo tanto, en este rango se acumula el 50% de las observaciones. Dado que los cuantiles son estadísticas robustas, el rango intercuarílico también lo es.

\[ RIC = Q3(X) - Q1(X) \]

Varianza: es quizás la medida de dispersión más comúnmente utilizada. Mide la desviación o error cuadrático promedio de las observaciones alrededor de la media. Comúnmente se utiliza el símbolo $\sigma^2$ para representarla cuando se habla de la varianza poblacional:

\[ \sigma^2 = \frac{\sum\limits_{i=1}^N (X_i - \bar{X})^2}{N} \]

sin embargo, cuando nos referimos a la varianza muestral, típicamente empleamos la notación $s^2$ para representarla, y la fórmula es ligeramente diferente:

\[ s^2 = \frac{\sum\limits_{i=1}^n (X_i - \bar{X})^2}{n - 1} \]

Un problema de la varianza es que está expresada en unidades de medición al cuadrado, por lo tanto, es muy frecuente que se utilice como referencia la raíz cuadrada (positiva) de la varianza, conocida como desviación estándar.

\[ \sigma = \sqrt{\sigma^2};\\ s = \sqrt{s^2}. \]

La varianza, además, dado que es una medida de dispersión alrededor de la media, y como la media no es una medida robusta, tampoco es una medida robusta.

Coeficiente de variación: cuando es necesario comparar la dispersión o variabilidad de dos variables diferentes, tenemos el problema de que las variables muy probablemente estén expresadas en unidades diferentes. Una manera de resolver este problema es expresando la dispersión de manera relativa. El coeficiente de variación, entonces, mide la magnitud de la variabilidad en relación a la magnitud de la media de la variable utilizando la desviación estándar como medida de dispersión:

\[ CV = \frac{\sigma}{\mu}. \]

El coeficiente de variación no tiene, entonces, unidades de medición, por lo que permite comparar entre variables.

Ejemplo: Calcula las medidas de dispersión para los datos de la renta registrados en la ENIGH.

Rango:

max(renta) - min(renta)

## [1] 200400

Rango intercuartílico:

quantile(x = renta, probs = c(0.25, 0.75))

##  25%  75% 
## 1100 3000

IQR(x = renta)

## [1] 1900

Varianza:

sum((renta - mean(renta))^2)/(length(renta))

## [1] 9018951

sum((renta - mean(renta))^2)/(length(renta) - 1)

## [1] 9019052

var(x = renta)

## [1] 9019052

sqrt(var(x = renta))

## [1] 3003.174

sd(x = renta)

## [1] 3003.174

Coeficiente de variación:

sd(renta) / mean(renta)

## [1] 1.14139

DIAGRAMAS DE CAJA Y BRAZOS

Una manera muy práctica de presentar visualmente varias de las medidas presentadas hasta este punto es el conocido como diagrama de caja y brazos. Este tipo de diagramas son muy útiles porque visualmente nos permiten identificar y juzgar varias de las características de la distribución de nuestros datos:

Centralidad: la centralidad estará representada por la mediana (aunque hay versiones de estos gráficos que también incluyen a la media)
Sesgo: mediante la representación de la acumulación de los datos en cuartiles podemos observar si la distribución parece simétrica o sesgada
Dispersión: la dispersión de los datos estará representada mediante el rango intercuartílico
Valores atípicos: se establece un criterio para la determinación de valores atípicos que serán representados como puntos aislados en el gráfico

Estos elementos es más fácil explicarlos sobre un ejemplo:

Ejemplo: Construye el diagrama de caja y brazos de la renta declarada por los encuestados en la ENIGH en la delegación Álvaro Obregón de la CDMX.
R

summary(renta[which(viviendas.enigh$ubica_geo=='09010')])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     700    3000    4000    6655    7000   40000

m <- min(renta[which(viviendas.enigh$ubica_geo=='09010')])
M <- max(renta[which(viviendas.enigh$ubica_geo=='09010')])

Q.1 <- quantile(renta[which(viviendas.enigh$ubica_geo=='09010')], probs = 0.25)
Q.2 <- quantile(renta[which(viviendas.enigh$ubica_geo=='09010')], probs = 0.50)
Q.3 <- quantile(renta[which(viviendas.enigh$ubica_geo=='09010')], probs = 0.75)

ric <- IQR(renta[which(viviendas.enigh$ubica_geo=='09010')])

uw <- renta[which(viviendas.enigh$ubica_geo=='09010')]
uw <- max(uw[which(uw <= Q.3 + 1.5*ric)])

ggplot(
  data = data.frame(Renta = renta[which(viviendas.enigh$ubica_geo=='09010')])
  ) +
  aes(x = Renta) +
  geom_boxplot() +
  theme_classic() +
  scale_x_continuous(labels = dollar, breaks = c(m, Q.1, Q.2, Q.3, uw, M)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(labels = NULL, breaks = NULL) + 
  labs(title = 'Diagrama de caja y brazos', subtitle = 'ENIGH: rentas declaradas en viviendas de Álvaro Óbregon', y = NULL)

ggplot(
  data = data.frame(Renta = renta[which(viviendas.enigh$ubica_geo=='09010')])
  ) +
  aes(x = Renta, y = 0) +
  geom_violin() +
  theme_classic() +
  scale_x_continuous(labels = dollar, breaks = c(m, Q.1, Q.2, Q.3, uw, M)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(labels = NULL, breaks = NULL) + 
  labs(title = 'Diagrama de violín', subtitle = 'ENIGH: rentas declaradas en viviendas de Álvaro Óbregon', y = NULL)

AED BI-VARIADO

Con frecuencia queremos explorar el comportamiento de una variable en presencia o interacción con otras variables de nuestros datos. Desde luego, mientras más variables queramos explorar de manera simultánea, más complejo será el análisis. Nosotros nos enfocaremos en el análisis de dos variables.

Tenemos que considerar, entonces, los tipos de datos que podemos tener. Esto nos da lugar a tres posibles diferentes alternativas².

CATEGÓRICA X CATEGÓRICA

Para analizar dos variables categóricas simultáneamente utilizaremos como herramienta principal las tablas cruzadas o de contingencia. Este tipo de tablas no son otra cosa más que tablas de frecuencia en las que los casos cuya frecuencia se cuentan son los casos en los que se observan las diferentes combinaciones de los valores de las dos variables.

Ejemplo: Para los datos de la ENIGH 2022, construye la tabla de contingencia del tipo de vivienda y el tipo de tenencia de la vivienda (registrado en la pregunta 23).

Clave	Tenencia
1	Rentada
2	Prestada
3	Propia, pagando
4	Propia, pagada
5	En litigio
6	Otra

tipo_tenencia_renta <- data.frame(Tipo = viviendas.enigh$tipo_viv, Tenencia = viviendas.enigh$tenencia, Renta = renta)

table(tipo_tenencia_renta$Tipo, tipo_tenencia_renta$Tenencia)

##    
##         1     2     3     4     5     6
##   &    52    12     4    37     2     0
##   1  9294  9967  6701 55515  2068   714
##   2   993   170   307   632    25    14
##   3   928   257    12   205    26    14
##   4    40    17     3    40     1     2
##   5   130    83    87   435    26    10

100*table(tipo_tenencia_renta$Tipo, tipo_tenencia_renta$Tenencia)/nrow(viviendas.enigh)

##    
##                1            2            3            4            5
##   &  0.058543395  0.013510014  0.004503338  0.041655877  0.002251669
##   1 10.463506074 11.221192709  7.544217151 62.500703647  2.328225797
##   2  1.117953683  0.191391869  0.345631199  0.711527420  0.028145863
##   3  1.044774439  0.289339473  0.013510014  0.230796078  0.029271698
##   4  0.045033381  0.019139187  0.003377504  0.045033381  0.001125835
##   5  0.146358488  0.093444266  0.097947604  0.489738018  0.029271698
##    
##                6
##   &  0.000000000
##   1  0.803845851
##   2  0.015761683
##   3  0.015761683
##   4  0.002251669
##   5  0.011258345

temp <- 
  aggregate(
    x = list(Frecuencia = tipo_tenencia_renta$Tipo), 
    by = 
      list(
        Tipo = tipo_tenencia_renta$Tipo, 
        Tenencia = tipo_tenencia_renta$Tenencia
        ), 
    FUN = length
    )

temp <- dcast(data = temp, formula = Tipo ~ Tenencia, fill = 0)

## Using Frecuencia as value column: use value.var to override.

Total <- sum

rownames(temp) <- temp[,1]

temp <- addmargins(A = as.matrix(temp[,-1]), FUN = Total)

## Margins computed over dimensions
## in the following order:
## 1: 
## 2:

kable(x = temp, format = 'html', format.args = list(big.mark = ',')) |> 
  kable_classic()

	1	2	3	4	5	6	Total
&	52	12	4	37	2	0	107
1	9,294	9,967	6,701	55,515	2,068	714	84,259
2	993	170	307	632	25	14	2,141
3	928	257	12	205	26	14	1,442
4	40	17	3	40	1	2	103
5	130	83	87	435	26	10	771
Total	11,437	10,506	7,114	56,864	2,148	754	88,823

Desde luego, podemos presentar esta misma información en gráficas de barras. Cuando el número de valores de la variable categórica es pequeño, es posible presentar las barras hombro con hombro, aunque esto no siempre será práctico o viable.

Ejemplo: Construye una gráfica de barras en la que se muestre la frecuencia registrada en la ENIGH 2022 del tipo de vivienda “Casa independiente” y “Departamento en edificio” por tipo de tenencia de la vivienda.

ggplot(data = tipo_tenencia_renta[which(tipo_tenencia_renta$Tipo == '1'|tipo_tenencia_renta$Tipo == '2'),]) +
  aes(x = Tenencia, group = Tipo, fill = factor(Tipo)) +
  geom_bar(position = position_dodge()) +
  labs(y = 'Frecuencia', fill = 'Tipo', title = 'Frecuencia del tipo de tenencia de vivienda por tipo de vivienda', subtitle = 'Barras hombro con hombro') +
  theme_classic() +
  scale_y_continuous(labels = comma)

ggplot(data = tipo_tenencia_renta[which(tipo_tenencia_renta$Tipo == '1'|tipo_tenencia_renta$Tipo == '2'),]) +
  aes(x = Tenencia, group = Tipo, fill = factor(Tipo)) +
  geom_bar() +
  labs(y = 'Frecuencia', fill = 'Tipo', title = 'Frecuencia del tipo de tenencia de vivienda por tipo de vivienda', subtitle = 'Barras superpuestas') +
  theme_classic() +
  scale_y_continuous(labels = comma)

CATEGÓRICA X CUANTITATIVA

Cuando tenemos una combinación de una variable categórica con una variable cuantitativa y queremos hacer el análisis de su interacción estaríamos buscando dar respuesta a la pregunta de cómo se comporta una variable condicional en la otra.

Si primero nos preguntamos cómo se comporta la variable categórica dado que observamos ciertos valores de la variable cuantitativa, estaríamos pensando en realizar una tabla de frecuencias de la variable categórica para los valores observados de la variable cuantitativa. Si los valores observados de la variable cuantitativa son muchos, entonces es necesario establecer rangos y, por lo tanto, caemos en el caso de dos variables categóricas.

Si pensamos en dar respuesta a la pregunta de cómo se comporta la variable cuantitativa dado que observamos ciertos valores de la variable categórica, entonces podemos pensar en que la variable categórica segmenta a nuestros datos para la variable cuantitativa. En estos casos, entonces, nos interesa comparar los atributos de posición, dispersión, sesgo y curtosis de la variable cuantitativa después de segmentarla usando la variable categórica.

Previamente hablamos de las gráficas de caja y brazos. Este tipo de gráficos son particularmente útiles para comparar la distribución o comportamiento de una variable cuantitativa segmentada utilizando una variable categórica.

Ejemplo: Compara la distribución de las rentas declaradas por los encuestados en la ENIGH 2022 de las alcaldías Álvaro Obregón y Benito Juárez de la CDMX.
R: Para construir estos diagramas, es necesario primero identificar que en la ENIGH las alcaldías solicitadas se identifican con las claves ‘09014’ y ‘09010’ (Benito Juárez y Álvaro Obregón, respectivamente), en el campo ‘ubica_geo’ de la tabla con las características de las viviendas.

Podemos entonces ya construir nuestros diagramas:

temp <- data.frame(
      Renta = 
        renta[
          which(
            viviendas.enigh$ubica_geo =='09014' | 
              viviendas.enigh$ubica_geo =='09010'
            )
          ]
      , Alcaldía = 
        viviendas.enigh[
          which(
            viviendas.enigh$ubica_geo =='09014' | 
              viviendas.enigh$ubica_geo =='09010'
            ),
          'ubica_geo'
        ]
    )

temp$Alcaldía <- factor(temp$Alcaldía)

by(data = temp$Renta, INDICES = temp$Alcaldía, FUN = summary)

## temp$Alcaldía: 09010
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     700    3000    4000    6655    7000   40000 
## ------------------------------------------------------------ 
## temp$Alcaldía: 09014
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    3500    8000   12000   12707   15000   30000

ggplot(data = temp) +
  aes(x = Renta, y = Alcaldía) +
  geom_boxplot() +
  theme_classic() +
  scale_x_continuous(labels = dollar) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    title = 'Diagrama de caja y brazos'
    , subtitle = 'ENIGH: rentas declaradas en viviendas de Álvaro Óbregon y Benito Juárez'
    , caption = '09010: AO; 09014: BJ'
    )

ggplot(data = temp) +
  aes(x = Renta, group = Alcaldía, colour = factor(Alcaldía)) +
  geom_density() +
  theme_classic() +
  scale_x_continuous(labels = dollar) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    title = 'Diagrama de caja y brazos'
    , subtitle = 'ENIGH: rentas declaradas en viviendas de Álvaro Óbregon y Benito Juárez'
    , caption = '09010: AO; 09014: BJ'
    , colour = 'Alcaldía'
    , y = 'Densidad'
    )

CUANTITATIVA X CUANTITATIVA

Cuando contamos con dos variables cuantitativas, un primer análisis de interés es, con frecuencia, simplemente el comparar las características del comportamiento de las dos variables. Desde luego, podemos hacer esto numéricamente (calculando sus medidas de posición y dispersión), pero la manera más inmediata de lograrlo es comparando gráficamente su distribución.

Ejemplo: Para la ENIGH 2022 surge la hipótesis de que las personas, en general, a mayor ingreso, mayor renta (o el equivalente de renta) pagan (mayor valor de su vivienda).

En la base de datos del concentrado de hogares se incluye el dato del ingreso corriente trimestral de las personas que componen al hogar. Se desea explorar la relación entre este dato y el de la renta construido anteriormente. ¿Cómo se compara la distribución de estas dos variables?

concentradohogar.enigh <- read.dbf(file = 'datos/concentradohogar.dbf')

ingreso.vivienda <- 
  aggregate(
    x = list(Ingreso = concentradohogar.enigh$ing_cor)
    , by = list(Vivienda = concentradohogar.enigh$folioviv)
    , FUN = sum
    )

ingreso.renta.vivienda <- 
  data.frame(Vivienda = viviendas.enigh$folioviv, Renta = renta) |> 
  merge(y = ingreso.vivienda, all.x = TRUE)

summary(ingreso.renta.vivienda)

##        Vivienda         Renta           Ingreso       
##  0100005002:    1   Min.   :    50   Min.   :      0  
##  0100005003:    1   1st Qu.:  1100   1st Qu.:  28798  
##  0100005004:    1   Median :  2000   Median :  46831  
##  0100012002:    1   Mean   :  2631   Mean   :  62375  
##  0100012004:    1   3rd Qu.:  3000   3rd Qu.:  75565  
##  0100012006:    1   Max.   :200450   Max.   :7153770  
##  (Other)   :88817

ingreso.renta.vivienda <- melt(data = ingreso.renta.vivienda, id.vars = c('Vivienda'), measure.vars = c('Renta', 'Ingreso'), value.name = 'Valor')

ggplot(data = ingreso.renta.vivienda) +
  aes(x = Valor) +
  facet_wrap(facets = 'variable', ncol = 1) +
  geom_boxplot() +
  theme_classic() +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = NULL)

ggplot(data = ingreso.renta.vivienda) +
  aes(x = Valor) +
  facet_wrap(facets = 'variable', scales = 'free', ncol = 1) +
  geom_boxplot() +
  theme_classic() +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = NULL)

ggplot(data = ingreso.renta.vivienda) +
  aes(x = log(Valor)) +
  facet_wrap(facets = 'variable', ncol = 1) +
  geom_boxplot() +
  theme_classic() +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = NULL)

## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

ggplot(data = ingreso.renta.vivienda) +
  aes(x = log(Valor)) +
  facet_wrap(facets = 'variable', scales = 'free', ncol = 1) +
  geom_boxplot() +
  theme_classic() +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = NULL)

## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

ggplot(data = ingreso.renta.vivienda) +
  aes(x = Valor, group = variable) +
  geom_density() +
  theme_classic() +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = NULL)

ggplot(data = ingreso.renta.vivienda) +
  aes(x = Valor) +
  facet_wrap(facets = 'variable', ncol = 1, scales = 'free_y') +
  geom_density() +
  theme_classic() +
  scale_x_continuous(labels = dollar)

ggplot(data = ingreso.renta.vivienda) +
  aes(x = Valor) +
  facet_wrap(facets = 'variable', scales = 'free', ncol = 1) +
  geom_density() +
  theme_classic() +
  scale_x_continuous(labels = dollar)

ggplot(data = ingreso.renta.vivienda) +
  aes(x = log(Valor), group = variable) +
  geom_density() +
  theme_classic() +
  scale_x_continuous(labels = dollar)

## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_density()`).

ggplot(data = ingreso.renta.vivienda) +
  aes(x = log(Valor)) +
  facet_wrap(facets = 'variable', ncol = 1) +
  geom_density() +
  theme_classic() +
  scale_x_continuous(labels = dollar)

## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_density()`).

ggplot(data = ingreso.renta.vivienda) +
  aes(x = log(Valor)) +
  facet_wrap(facets = 'variable', scales = 'free', ncol = 1) +
  geom_density() +
  theme_classic() +
  scale_x_continuous(labels = dollar)

## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_density()`).

Es importante observar que, si el objetivo es comparar una distribución con otra, para que dicha comparación sea válida, generalmente se buscará que el eje correspondiente a los valores de las variables, se encuentre expresado en las mismas unidades y cubra el mismo intervalo, de otra manera la comparación puede ser engañosa. Sin embargo, una vez analizados los datos en términos comparables, separar los rangos de valores, o incluso transformar las variables, nos puede ayudar a apreciar mejor algunas de las características y diferencias entre las distribuciones.

En el caso en el que deseamos explorar la relación existente entre dos variables cuantitativas, de igual manera, nos interesa explorar el impacto que tiene un nivel observado de una de las variables sobre el nivel observado en la otra. En otras palabras, si existe algún tipo de relación entre las variables.

Visualmente, la manera más sencilla de explorar esto es mediante un diagrama de dispersión. Para construir un diagrama de dispersión, simplemente representamos en un gráfico las parejas ordenadas de observaciones de cada una de las variables. Desde luego, esto implica que las variables cuantitativas representan dos diferentes mediciones de atributos de los mismos individuos de la población (de otra manera no haría ningún sentido ordenarlas en parejas).

ggplot(data = data.frame(x = rnorm(n = 100), y = rnorm(n = 100))) +
  aes(x = x, y = y) +
  geom_point()

La idea es, entonces, utilizar el diagrama de dispersión para visualmente intentar identificar patrones relevantes (tendencias, asociaciones, agrupaciones, observaciones atípicas, etc.)

Ahora bien, una de las relaciones que con mayor frecuencia nos interesa explorar es la asociación lineal entre las variables (esto es, nos interesa verificar si la relación entre las variables puede ser aproximada mediante una linea recta).

La intensidad de la presencia de esta relación lineal entre dos variables se puede medir mediante el coeficiente de correlación:

\[ \rho_{XY} = \frac{\frac{\sum(x_i - \mu_X)(y_i - \mu_Y)}{n-1}}{\sigma_X \sigma_Y} = \frac{Cov(X,Y)}{\sigma_X \sigma_Y} \]

El coeficiente de correlación tomará valores entre $-1$ y $1$. Los extremos corresponderán, respectivamente, a una correlación lineal perfecta negativa (es decir, que cuando una variable aumenta su valor, la otra lo disminuye) o positiva (cuando ambas variables se mueven en el mismo sentido).

Aunque no existen valores estándar para decidir cuándo dos variables están altamente correlacionadas (y, de alguna manera, los valores de referencia dependerán del contexto)

Ejemplo: Siguiendo con el ejemplo de la ENIGH para las variables de ingreso y renta, si se desea explorar esta relación para posteriormente ajustar un modelo que permita explicar el ingreso de las personas a partir de la renta que pagan por su vivienda.

En la base de datos del concentrado de hogares se incluye el dato del ingreso corriente trimestral de las personas que componen al hogar. Calcula el ingreso corriente trimestral registrado por vivienda y explora la relación entre este dato y el de la renta construido anteriormente. ¿Parece existir alguna relación como la que se sugiere en la investigación?

concentradohogar.enigh <- read.dbf(file = 'datos/concentradohogar.dbf')

ingreso.vivienda <- aggregate(x = list(Ingreso = concentradohogar.enigh$ing_cor), by = list(Vivienda = concentradohogar.enigh$folioviv), FUN = sum)

summary(ingreso.vivienda)

##        Vivienda        Ingreso       
##  0100005002:    1   Min.   :      0  
##  0100005003:    1   1st Qu.:  28798  
##  0100005004:    1   Median :  46831  
##  0100012002:    1   Mean   :  62375  
##  0100012004:    1   3rd Qu.:  75565  
##  0100012006:    1   Max.   :7153770  
##  (Other)   :88817

ggplot(data = ingreso.vivienda) +
  aes(x = Ingreso) +
  geom_boxplot() +
  theme_classic() +
  scale_x_continuous(labels = dollar)

ingreso.renta.vivienda <- data.frame(Vivienda = viviendas.enigh$folioviv, Renta = renta) |> merge(y = ingreso.vivienda, all.x = TRUE)

ggplot(data = ingreso.renta.vivienda) +
  aes(x = Renta, y = Ingreso) +
  geom_point() +
  theme_classic() +
  scale_x_continuous(labels = dollar) +
  scale_y_continuous(labels = dollar)

var(x = ingreso.renta.vivienda$Renta)

## [1] 9019052

var(x = ingreso.renta.vivienda$Ingreso)

## [1] 6235389304

cov(x = ingreso.renta.vivienda$Renta, y = ingreso.renta.vivienda$Ingreso)

## [1] 94807474

cor(x = ingreso.renta.vivienda$Renta, y = ingreso.renta.vivienda$Ingreso)

## [1] 0.3997885

Los datos no parecen sugerir una relación útil para los fines propuestos.

AED - CONCLUSIONES

Preguntas generales que nos podemos plantear a partir de este análisis:

¿Cuál es la población representada en los datos? ¿Cómo impacta esta definición de la población en las características de los datos?
- Los datos, ¿representan muestras o son resultado de censos?
- Si se trata de muestras, ¿cómo se construyeron las muestras?
¿Se observan características “inesperadas” en los datos? (en otras palabras, ¿se observar rasgos o datos atípicos o problemas de “calidad”?)
¿Es posible formular alguna hipótesis sobre el comportamiento de las variables?
¿Qué tipo de inferencia estadística es necesario/posible realizar? ¿Qué conjunto de técnicas se pueden contemplar?
¿Es necesario conseguir más datos?

Al final, todas la medidas y gráficos que vimos son simplificaciones (generalmente) de estructuras de datos complejas. Simplificaciones que se enfocan en alguna característica en particular. No son, por lo tanto, perfectas sino simplemente indicadores de elementos a explorar a mayor profundidad o verificar mediante técnicas de inferencia estadística.

ESTADÍSTICA INFERENCIAL

Lectura recomendada

Wackerly, Mendenhall, and Scheaffer (2008), cap. 1.1.

Hasta aquí, hemos visto algunas de las técnicas desarrolladas para describir o explorar un fenómeno u objeto de estudio a través de datos recolectados respecto de su comportamiento. El otro gran uso que se hace de la estadística es para realizar inferencia sobre las características del objeto de nuestro estudio. Es decir, usamos a la estadística con el objetivo de deducir o sacar conclusiones sobre nuestro objeto de estudio y, en particular, hacemos uso de las técnicas estadísticas cuando buscamos hacer inferencia en presencia de incertidumbre. En este sentido, las técnicas estadísticas son deductivas mientras que la teoría de probabilidad es inductiva.

Una de las circunstancias más frecuentes en las que esta incertidumbre se hace presente es cuando se desea sacar conclusiones sobre una población (sobre sus características) a partir de información recogida de una muestra de esa población. Típicamente, además, nos interesa no solamente expresar el resultados de dichas conclusiones sino que, dado que nos encontramos en presencia de incertidumbre, es deseable también que podamos pronunciarnos sobre qué tan confiables (qué tan buenas) son nuestras conclusiones.

DISTRIBUCIONES MUESTRALES Y TEOREMA CENTRAL DEL LÍMITE

Lecturas recomendadas

Wackerly, Mendenhall, and Scheaffer (2008), cap. 1.4, 7.
Aguirre et al. (2006), cap. 5.
DeGroot (1988), cap. 7.

Para poder realizar inferencias estadísticas, la herramienta principal es la probabilidad. La estadística hace uso de la probabilidad para que, partiendo de supuestos probabilísticos sobre la muestra, podamos sacar conclusiones sobre las características de la población. Por lo tanto, en adelante vamos a trabajar asumiendo, en general, que contamos con una muestra $X = \{X_1, X_2, \dots, X_n\}$ tomada aleatoriamente de una población de interés.

Nuestro objetivo, en general, será entonces construir funciones de esta muestra aleatoria que podamos usar como estadísticas para poder hacer inferencia sobre los parámetros poblacionales.

Si consideramos que el valor de la variable $X$ es, entonces, aleatorio, las estadísticas construidas serán por lo tanto funciones de una variable aleatoria y, en consecuencia, serán a su vez variables aleatorias. Tendrán asociada, por lo tanto, una distribución aleatoria, a la cual llamaremos la distribución muestral de la estadística.

La distribución muestral de una estadística dependerá entonces totalmente de la distribución de la variable aleatoria $X$, distribución que típicamente se nos dará como un supuesto de trabajo pero de la cual desconocemos sus parámetros, mismos que deseamos estimar (o sobre los que deseamos hacer algún otro tipo de inferencia) a través de la información contenida en la muestra recogida.

En ocasiones, lo que se busca es determinar ciertas características de la muestra (p.e., su tamaño) que permitan cumplir con determinadas restricciones impuestas al experimento (p.e., nos puede interesar que el parámetro poblacional sea estimado dentro de un determinado margen de error). En estos casos utilizamos la distribución muestral del estimador para determinar los valores de interés antes del levantamiento de la muestra.

En términos generales, para obtener la distribución muestral del estadístico podemos hacerlo analíticamente (esto es, calcular la función de distribución de la función de la variable aleatoria) o bien, podemos recurrir a la simulación.

Ejemplo: Supongamos que se tiene una moneda y se jugarán volados con esa moneda. Se hacen tres tiros y se calcula su promedio (águila = 1; sol = 0). Sea la estadística $\bar{X} = \frac{X_1 + X_2 + X_3}{3}$. ¿Cuál es la distribución muestral de $\bar{X}$?
R: Cada volado es una v.a. $Bernoulli(0.5)$. La función de distribución del promedio de los volados es la misma que la de la suma de los volados, y sabemos que la suma de vv.aa. Bernoulli corresponde a una función de distribución $Binomial(n,p)$. Por lo tanto

\[ \bar{X} \sim Binomial(3, 0.5). \]

Otra manera de resolverlo es observando que la variable puede tomar los siguientes valores:

$\sum X_i$	$\frac{\sum X_i}{3}$	p	P
0	0	$\frac{1}{8}$	$\frac{1}{8}$
1	1/3	$\frac{3}{8}$	$\frac{4}{8}$
2	2/3	$\frac{3}{8}$	$\frac{7}{8}$
3	1	$\frac{1}{8}$	$\frac{8}{8}$

Finalmente, podríamos simular:

temp <- data.frame(X.1 = rbinom(n = 10, size = 1, prob = 0.5), X.2 = rbinom(n = 10, size = 1, prob = 0.5), X.3 = rbinom(n = 10, size = 1, prob = 0.5)) |> transform(X.barra = (X.1 + X.2 + X.3)/3)

ggplot(data = temp) +
  aes(x = X.barra, y = after_stat(prop)) +
  geom_bar() +
  theme_classic() +
  labs(y = 'Frecuencia relativa') +
  scale_x_continuous(breaks = c(0,1/3,2/3,1))

temp <- data.frame(X.1 = rbinom(n = 100, size = 1, prob = 0.5), X.2 = rbinom(n = 100, size = 1, prob = 0.5), X.3 = rbinom(n = 100, size = 1, prob = 0.5)) |> transform(X.barra = (X.1 + X.2 + X.3)/3)

ggplot(data = temp) +
  aes(x = X.barra, y = after_stat(prop)) +
  geom_bar() +
  theme_classic() +
  labs(y = 'Frecuencia relativa') +
  scale_x_continuous(breaks = c(0,1/3,2/3,1))

temp <- data.frame(X.1 = rbinom(n = 1000, size = 1, prob = 0.5), X.2 = rbinom(n = 1000, size = 1, prob = 0.5), X.3 = rbinom(n = 1000, size = 1, prob = 0.5)) |> transform(X.barra = (X.1 + X.2 + X.3)/3)

ggplot(data = temp) +
  aes(x = X.barra, y = after_stat(prop)) +
  geom_bar() +
  theme_classic() +
  labs(y = 'Frecuencia relativa') +
  scale_x_continuous(breaks = c(0,1/3,2/3,1))

En particular, dada su relevancia en una gran diversidad de problemas de inferencia, nos vamos a enfocar en distribuciones muestrales de estadísticos que provienen de una muestra de observaciones independientes provenientes de una variable aleatoria normal. A continuación entonces estudiaremos algunas de las distribuciones muestrales que se originan a partir de la normal.

DISTRIBUCIONES MUESTRALES RELACIONADAS CON UNA POBLACIÓN NORMALMENTE DISTRIBUIDA

Si nuestras observaciones provienen de una muestra aleatoria de una variable normalmente distribuida con media $\mu$ y varianza $\sigma^2$, entonces, recordemos que:

$f(x|\mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma}} e^{-\frac{(x - \mu)}{2\sigma^2}}$
Si $X \sim N(\mu, \sigma^2)$ entonces $a + bX \sim N(\mu + a, a^2 \sigma^2)$.
Dado que se trata de observaciones independientes, tenemos que:

\[ f(x_1, x_2, \dots, x_n|\mu, \sigma^2) = \prod\limits_{i = 1}^n \frac{1}{\sqrt{2 \pi \sigma}} e^{-\frac{(x_i - \mu)}{2\sigma^2}} \]

La suma de vv.aa. normales independientes e idénticamente distribuidas es, a su vez, una v.a. normal con media $n \mu$ y varianza $n \sigma^2$.

Ejemplo: Calcula la distribución muestral de la media muestral de una variable aleatoria distribuida $N(\mu, \sigma^2)$.
R

$\chi^2$

Ejemplo (cont.): Para los datos del ingreso declarado por las personas encuestadas en la ENIGH 2022, asume que el logaritmo del ingreso se distribuye normal. ¿Cuál es el cuantil 97.5 de la varianza muestral del logaritmo del ingreso?
R: Para dar respuesta a esta pregunta necesitamos obtener la distribución muestral de la varianza muestral. Nos dicen que nuestra variable de interés se distribuye normal por lo que observamos que

\[ s^2 = \frac{\sum\limits_{i = 1}^n (X_i - \bar{X})^2}{n - 1} \]

donde $X_i$ es el logaritmo del ingreso reportado por la i-ésima persona y $\bar{X}$ es la media muestral correspondiente. Y sabemos que $X_i$ es una v.a. normal.

Para poder proporcionar el cuantil 97.5 de $s^2$ necesitamos, entonces, su distribución muestral. Nos preguntan entonces por la distribución de la suma del cuadrado de vv.aa. normales.

Si, como lo planteamos inicialmente tenemos $n$ observaciones independientes de una v.a. $Y \sim N(\mu, \sigma^2)$ y definimos $Z = \left(\frac{Y - \mu}{\sigma}\right)^2$ entonces decimos que una nueva v.a. $W = \sum\limits_{i=1}^n Z_i$ es una variable aleatoria $\chi_{n}^2$.

La función de densidad de $W$ es igual a

\[ f(w|n) = \left\{ \begin{array}{} \frac{w^{\frac{n}{2}-1}e^{-\frac{w}{2}}}{2^{\frac{n}{2}}\Gamma(\frac{n}{2})} & w > 0\\ 0 & \text{e.o.c} \end{array}\right.. \]

Donde $\Gamma(x) = \int\limits_{0}^{\infty} t^{x-1}e^tdt$.

Ejemplo (cont.): ¿Podemos dar ya respuesta a la pregunta inicial?
R: Observemos primero que

\[ X_i - \bar{X} = \left(1 - \frac{1}{n}\right)X_i - \frac{1}{n} \sum\limits_{j \neq i} X_j. \]

Entonces cada $X_i - \bar{X}$ es una combinación lineal de vv.aa. independientes distribuidas $N(\mu, \sigma^2)$, y por lo tanto serán, a su vez, vv.aa. pero con

\[ E(X_i - \bar{X}) = E\left[\left(1 - \frac{1}{n}\right)X_i - \frac{1}{n} \sum\limits_{j \neq i} X_j\right]\\ = E\left[\left(1 - \frac{1}{n}\right)X_i\right] - E\left[\frac{1}{n} \sum\limits_{j \neq i}X_j\right]\\ = \left(1 - \frac{1}{n}\right)E\left[X_i\right] - \frac{1}{n} E\left[ \sum\limits_{j \neq i}X_j\right]\\ = \left(1 - \frac{1}{n}\right)\mu - \frac{n-1}{n} \mu\\ = 0. \]

\[ Var(X_i - \bar{X}) = Var\left[\left(1 - \frac{1}{n}\right)X_i - \frac{1}{n} \sum\limits_{j \neq i} X_j\right]\\ = Var\left[\left(1 - \frac{1}{n}\right)X_i\right] + Var\left[\frac{1}{n} \sum\limits_{j \neq i}X_j\right]\\ = \left(1 - \frac{1}{n}\right)^2 Var\left[X_i\right] - \frac{1}{n^2} E\left[ \sum\limits_{j \neq i}X_j\right]\\ = \left(\frac{n-1}{n}\right)^2\sigma^2 - \frac{n-1}{n^2} \sigma^2\\ = \left(\frac{n^2-2n + 1 + n - 1}{n^2}\right)\sigma^2\\ = \left(\frac{n^2- n}{n^2}\right)\sigma^2\\ = \frac{n - 1}{n}\sigma^2. \]

Tenemos entonces la suma de vv.aa.i.dd. $N(0, \frac{n - 1}{n}\sigma^2)$ pero no sabemos si son independientes. De hecho, no lo son. Podemos ver el caso en el que $n = 2$. En este caso es fácil ver que $s^2 = \frac{(X_1 - X_2)^2}{2}$ por lo que podemos ver que la $s^2$ puede ser expresada como el cuadrado de una v.a. normal.

En general, es posible probar que $s^2$ representa la suma del cuadrado de $n-1$ vv.aa. independientes normales y que

\[ \frac{(n-1)s^2}{\sigma^2} \sim \chi_{n-1}^2. \]

Podemos entonces usar este resultado para contestar lo solicitado. Para los datos de la ENIGH 2022 sabemos que $n = 88823$ entonces buscamos $c$ tal que:

\[ P\left[\frac{(n-1)s^2}{\sigma^2} \leq c\right] = 0.975\\ P\left[s^2 \leq \frac{c \times \sigma^2}{n-1}\right] = 0.975. \]

Observemos que $c$ corresponde al cuantil de la $\chi_{88822}^2$ por lo que conocemos todos los datos excepto $\sigma^2$. No podemos dar respuesta, en realidad, a lo que nos piden, ya que $\sigma^2$ es desconocida. Sin embargo, para fines de ejemplo, supongamos que se nos indica que los modelos teóricos sugieren que la desviación estándar del logaritmo del ingreso es igual a 0.75, entonces ahora sí podemos determinar el cuantil 97.5 de $s^2$ ya que:

\[ P\left[s^2 \leq \frac{c \times \sigma^2}{n-1}\right] = 0.975\\ P\left[s^2 \leq \frac{89649.98 \times 0.75^2}{88822}\right] = 0.975\\ P\left[s^2 \leq 0.5677\right] = 0.975\\ \]

Podemos ver que el valor de la varianza muestral del logaritmo de los ingresos obtenidos en la muestra se encuentra, efectivamente, por debajo del valor de dicho cuantil.

$t$

Regresamos ahora a la distribución muestral de $\bar{X}$ como caso de uso. Dijimos ya que $\bar{X} \sim N(\mu, \frac{\sigma^2}{n})$ por lo que $\frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}} \sim N(0, 1)$. Si quisiéramos hacer inferencia sobre el parámetro $\mu$ a partir de la muestra, tenemos el problema de que dependemos del valor (desconocido) de $\sigma$. Sin embargo, podemos sustituir el valor de $\sigma$ por el de su estimador $s$ (la desviación estándar muestral). Desde luego, esta sustitución tiene el costo correspondiente a agregar incertidumbre (precisamente porque $s$ proviene de la muestra, no de la población), por lo que ya no podemos asegurar que $\frac{\bar{X} - \mu}{\frac{s}{\sqrt{n}}} \sim N(0, 1)$.

¿Cuál es entonces la distribución de este nuevo cociente?

Si definimos a la variable aleatoria $Y = \frac{Z}{\sqrt{\frac{W}{v}}}$ donde $Z \sim N(0,1)$ y $W \sim \chi_{v}^2$, entonces resulta que $Y \sim t_{v}$. Su función de densidad será igual a:

\[ f(y) = \frac{\Gamma((v+1)/2)}{\sqrt{2 \pi v} \ \Gamma(v/2)} \left( 1 + x^2/v\right)^{-(v+1)/2} \]

Entonces, podemos ahora ver que $\frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}} \sim N(0,1)$ y sabemos que $\frac{(n-1)s^2}{\sigma^2} \sim \chi_{n-1}^2$, por lo que

\[ \frac{\frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}}}{\sqrt{\frac{\frac{(n-1)s^2}{\sigma^2}}{(n-1)}}} \sim t_{n-1}\\ \frac{\sqrt{n}\frac{\bar{X} - \mu}{\sigma}}{\frac{s}{\sigma}} \sim t_{n-1}\\ \frac{\sqrt{n}(\bar{X} - \mu)}{s} \sim t_{n-1}\\ \frac{(\bar{X} - \mu)}{\frac{s}{\sqrt{n}}} \sim t_{n-1}. \]

Podemos ya, entonces, hacer análisis sobre $\bar{X}$ y $\mu$ sin necesidad de conocer el valor de $\sigma$.

Ejemplo: Para el logaritmo del ingreso reportado en la ENIGH 2022, obtén un intervalo dentro del que se pueda afirmar que se encuentra la diferencia de la media muestral y la media poblacional del logaritmo del ingreso con probabilidad de 0.75. Nuevamente asume que el logaritmo del ingreso sigue una distribución normal.
R: Nuestra variable de interés es el logaritmo del ingreso: $X = log(Ingreso) \sim N(\mu, \sigma^2)$. Sabemos entonces que $\frac{\bar{X} - \mu}{s/\sqrt{n}} \sim t_{n-1}$. Entonces:

log.ingreso <- 
  log(ingreso.renta.vivienda$Ingreso[which(ingreso.renta.vivienda$Ingreso>0)])

n <- length(log.ingreso)
n

## [1] 88814

x.barra <- mean(log.ingreso)
x.barra

## [1] 10.74978

s <- sd(x = log.ingreso)
s

## [1] 0.7454953

por lo que buscamos

\[ P\left[\left|\bar{X} - \mu\right| \leq c\right] = 0.75\\ P\left[\left|\frac{\bar{X} - \mu}{s/\sqrt{n}}\right| \leq \frac{c}{s/\sqrt{n}} \right] = 0.75. \]

Entonces resulta que

\[ \frac{c}{s/\sqrt{n}} = {_{0.875}}t_{n-1}\\ c = {_{0.875}}t_{n-1} \times s/\sqrt{n}\\ \]

c <- qt(p = 0.875, df = n-1) * s / sqrt(n)
c

## [1] 0.002877642

Entonces la diferencia entre la media poblacional y la media muestral del logaritmo del ingreso se encontrará en el intervalo $(-0.0028776, 0.0028776)$ con probabilidad de 0.75.

$F$

La última distribución muestral que estudiaremos que emana de la distribución normal nos permite dar respuesta a preguntas como la que se plantea en el siguiente ejemplo.

Ejemplo: Sabemos que la ENIGH 2022 cubre a todas las entidades federativas, por lo tanto, alguien nos consulta porque desea saber si consideramos que la variabilidad (poblacional) del logaritmo del ingreso es la misma entre los habitantes de la CDMX y los habitantes del estado de Oaxaca con base en las respuestas registradas en la encuesta.

¿Cómo podemos dar respuesta a esta pregunta?

Para dar respuesta a esta pregunta haremos uso de la distribución $F$.

Definición: Sean $W_1$ y $W_2$ dos vv.aa.ii. $\chi^2$ con $v_1$ y $v_2$ grados de libertad, respectivamente. Entonces decimos que la variable $Y = \frac{\frac{W_1}{v_1}}{\frac{W_2}{v_2}}$ sigue una distribución $F_{v_1,v_2}$.

Ejemplo: ¿Cómo podemos entonces usar la distribución $F$ para dar respuesta a la pregunta de nuestro cliente?
R: Lo primero a notar es que tenemos dos sub-poblaciones: habitantes de CDMX y habitantes de Oaxaca. Lo que nos piden entonces es validar si la varianza poblacional de ambas sub-poblaciones parece ser la misma ($\sigma_{CDMX}^2 = \sigma_{OAX}^2)$.

Sabemos, de secciones anteriores, que:

\[ \frac{(n_{CDMX}-1)s_{CDMX}^2}{\sigma_{CDMX}^2} \sim \chi_{CDMX}^2\\ \frac{(n_{OAX}-1)s_{OAX}^2}{\sigma_{OAX}^2} \sim \chi_{OAX}^2 \]

por lo que, si $\sigma_{CDMX}^2 = \sigma_{OAX}^2$ se debería cumplir que

\[ \frac{\frac{\frac{(n_{CDMX}-1)s_{CDMX}^2}{\sigma_{CDMX}^2}}{n_{CDMX}-1}}{\frac{\frac{(n_{OAX}-1)s_{OAX}^2}{\sigma_{OAX}^2}}{n_{OAX}-1}} \sim F_{n_{CDMX}-1,n_{OAX}-1}\\ \frac{\frac{s_{CDMX}^2}{\sigma_{CDMX}^2}}{\frac{s_{OAX}^2}{\sigma_{OAX}^2}} \sim F_{n_{CDMX}-1,n_{OAX}-1}\\ \frac{s_{CDMX}^2}{s_{OAX}^2} \sim F_{n_{CDMX}-1,n_{OAX}-1}. \]

vivienda.estado <- viviendas.enigh |> subset(select = c(folioviv, ubica_geo)) |> transform(estado = substr(x = ubica_geo, start = 1, stop = 2))

ingreso.renta.vivienda.2 <- merge(x = ingreso.renta.vivienda, y = vivienda.estado, by.x = 'Vivienda', by.y = 'folioviv', all.x = TRUE) |> subset(subset = Ingreso > 0 & (estado == '09' | estado == '20')) |> transform(log.ingreso = log(Ingreso))

Para los datos de la ENIGH 2022, obtenemos que:

n <- by(data = ingreso.renta.vivienda.2$log.ingreso, INDICES = ingreso.renta.vivienda.2$estado, FUN = length)
n

## ingreso.renta.vivienda.2$estado: 09
## [1] 2546
## ------------------------------------------------------------ 
## ingreso.renta.vivienda.2$estado: 20
## [1] 2577

s.2 <- by(data = ingreso.renta.vivienda.2$log.ingreso, INDICES = ingreso.renta.vivienda.2$estado, FUN = sd)
s.2

## ingreso.renta.vivienda.2$estado: 09
## [1] 0.7068676
## ------------------------------------------------------------ 
## ingreso.renta.vivienda.2$estado: 20
## [1] 0.7733871

De lo que se desprende que, si las varianzas poblacionales son iguales el valor $\frac{s_{CDMX}^2}{s_{OAX}^2} = \frac{0.7068676}{0.7733871} = 0.9139893$ sería una observación obtenida de una distribución $F_{2545,2577}$. ¿Cómo podemos utilizar esta información para dar respuesta a nuestro cliente?

Una manera de hacerlo es preguntarnos qué cuántil de la distribución teórica representa el valor observado. La idea de hacer esto es que, en la medida en que el valor observado represente un valor extremo de la distribución teórica, podríamos pensar que las muestras sugerirían que el valor que observamos no proviene realmente de esa distribución (dada su baja probabilidad de observarlo).

temp <- data.frame(F = rf(n = 10000, df1 = n[1], df2 = n[2]))

ggplot(data = temp) +
  aes(x = F) + 
  geom_density() +
  theme_classic() +
  labs(title = 'Función de densidad F')

pf(q = s.2[1]/s.2[2], df1 = n[1], df2 = n[2])

##         09 
## 0.01145403

ggplot(data = temp) +
  aes(x = F) + 
  geom_density() +
  geom_vline(xintercept = s.2[1]/s.2[2], colour = 'red', linetype = 'dashed') +
  theme_classic() +
  labs(title = 'Función de densidad F')

Juzgar qué tan extremo es este valor es, desde luego, subjetivo, pero estamos obteniendo como resultado que, si las varianzas poblacionales fueran iguales, deberíamos observar cocientes menores o iguales al obtenido con probabilidad de un 1%, por lo que parece haber cierta evidencia en contra de la idea de que las varianzas son iguales.

TEOREMA CENTRAL DEL LÍMITE

Ejemplo: Hasta ahora hemos estado trabajando con el logaritmo de los ingresos porque vimos en el AED que su comportamiento se asemejaba al de una distribución normal. Nuestro cliente nos dice que la escala logarítmica no es muy clara, que prefiere trabajar con la variable original.

ggplot(data = ingreso.renta.vivienda) +
  aes(x = Ingreso) +
  geom_density() +
  theme_classic() +
  labs(y = 'Densidad') +
  scale_x_continuous(labels = dollar)

Si de los datos de la encuesta tomáramos sub-muestras, ¿cómo se distribuiría la media del ingreso de las muestras?

temp <- list()

temp$temp.1 <- data.frame(Muestra <- 1:1000, Media = NA)
temp$temp.2 <- data.frame(Muestra <- 1:1000, Media = NA)
temp$temp.3 <- data.frame(Muestra <- 1:1000, Media = NA)
temp$temp.4 <- data.frame(Muestra <- 1:1000, Media = NA)

for (i in 1:4){
  for (j in 1:nrow(temp[[i]])){
    muestra <- sample(x = ingreso.renta.vivienda$Ingreso, size = 10**i)
    temp[[i]][j, 'Media'] <- mean(muestra)
  }
}

ggplot(data = temp$temp.1) +
  aes(x = Media) +
  geom_density() +
  geom_density(data = temp$temp.2, colour = 'blue') +
  geom_density(data = temp$temp.3, colour = 'purple') +
  geom_density(data = temp$temp.4, colour = 'red', fill = 'red', alpha = 0.1) +
  theme_classic() +
  scale_x_continuous(labels = dollar)

Como podemos ver, conforme aumentamos el número de muestras para las cuales calculamos las medias, la distribución de las medias comienza a tomar esa familiar forma de campana.

El fenómeno que observamos en los datos del ejemplo anterior es lo que se conoce como el Teorema Central del Límite:

Teorema: Sean $\{X_1, \dots, X_n\}$ $n$ variables aleatorias e idénticamente distribuidas con media $\mu$ y varianza $\sigma^2 < \infty$. Entonces

\[ \lim\limits_{n \rightarrow \infty} P\left[\frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}} \leq u\right] = \int\limits_{-\infty}^{u} \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} dt \]

Demostración: Primero, observemos que $U_n = \frac{\frac{\sum\limits_{i} X_i}{n} - \mu}{\frac{\sigma}{\sqrt{n}}} = \frac{\frac{\sum\limits_{i} (X_i - \mu)}{n}}{\frac{\sigma}{\sqrt{n}}} = \frac{1}{\sqrt{n}}\sum\limits_{i} \frac{(X_i - \mu)}{\sigma} = \frac{1}{\sqrt{n}} \sum\limits_{i} Z_i$.

Observa entonces que la función generadora de momentos de $U_n$ es igual a

\[ M_{U_n}(t) = E[e^{tU_n}]\\ = E\left[e^{t\frac{1}{\sqrt{n}} \sum\limits_{i} Z_i}\right]\\ = E\left[e^{\frac{t}{\sqrt{n}} \sum\limits_{i} Z_i}\right]\\ = M_{\sum\limits_{i} Z_i}\left(\frac{t}{\sqrt{n}}\right)\\ = \left[M_{Z_i}\left(\frac{t}{\sqrt{n}}\right)\right]^n. \]

Si hacemos ahora la expansión de la serie de Taylor de $M_{Z_i}\left(\frac{t}{\sqrt{n}}\right)$ obtenemos que

\[ M_{Z_i}\left(\frac{t}{\sqrt{n}}\right) = M_{Z_i}\left(0\right) + M_{Z_i}'\left(0\right)\left(\frac{t}{\sqrt{n}}\right) + M_{Z_i}''\left(\xi\right)\frac{\left(\frac{t}{\sqrt{n}}\right)^2}{2}, 0 < \xi < \frac{t}{\sqrt{n}}\\ M_{Z_i}\left(0\right) = 1\\ M_{Z_i}'\left(0\right) = E(Z_i) = 0\\ M_{Z_i}\left(\frac{t}{\sqrt{n}}\right) = 1 + M_{Z_i}''\left(\xi\right)\frac{\left(\frac{t}{\sqrt{n}}\right)^2}{2}, 0 < \xi < \frac{t}{\sqrt{n}}\\ \]

pero $\lim\limits_{n \rightarrow \infty} \xi = 0$ y $M_{Z_i}''\left(0\right) = Var(Z_i) - E^2(Z_i) = Var(Z_i) = 1$ por lo que

\[ M_{Z_i}\left(\frac{t}{\sqrt{n}}\right) = 1 + \frac{1}{2} \frac{t^2}{n} + \epsilon_n\\ [M_{Z_i}\left(\frac{t}{\sqrt{n}}\right)]^n = [1 + \frac{1}{2} \frac{t^2}{n} + \epsilon_n]^n\\ \lim\limits_{n \rightarrow n}[M_{Z_i}\left(\frac{t}{\sqrt{n}}\right)]^n = e^{\frac{t^2}{n}}. \]

La última expresión puede verse que corresponde a la f.g.m de una dist. normal estándar.

En otras palabras, el TCL nos dice que $\bar{X}$ se distribuye asintóticamente normal con media $\mu$ y varianza $\frac{\sigma^2}{n}$.

Ejemplo (cont.): Previo al levantamiento de la ENIGH, los expertos consideraban que el ingreso por vivienda debería tener una media poblacional de $60,000 y una desviación estándar de $75,000. ¿La información obtenida en la ENIGH es consistente con estas percepciones?
R: Dado que no conocemos realmente la distribución del ingreso, no podemos obtener la distribución de la media. Sin embargo, por el TCL sabemos que $\frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}} \sim N(0,1)$ entonces si para una muestra de tamaño 88823 se obtuvo una media de 6.2375382^{4} podemos calcular

\[ P\left[\left| \frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}}\right| > \frac{62375.38 - 60000}{75000/\sqrt{88823}}\right] = 2 \times \Phi(-9.44) \]

pnorm(q = -9.4439186)*2

## [1] 3.59093e-21

Los resultados de pruebas de rendimiento académico en las escuelas de un estado tienen media de 60 y varianza de 64. Se toma una muestra aleatoria de tamaño 100 es una escuela cuya media (muestral) fue de 58. ¿Crees que exista evidencia suficiente para asegurar que esta escuela tiene un rendimiento menor al de la población?

R: Utilizando el TLC sabemos que, para $n$ grande, $\frac{\bar{X} - \mu}{\frac{\sigma}{\sqrt{n}}} \sim N(0, 1)$ aproximadamente. Si asumimos que $n = 100$ es suficientemente grande, entonces podemos calcular

\[ P\left[\bar{X} \leq 58\right] = P\left[\bar{X} - \mu \leq -2\right]\\ = P\left[\frac{\bar{X} - \mu}{\frac{8}{10}} \leq \frac{-2}{\frac{8}{10}}\right]\\ = P\left[Z \leq -2.5 \right] \]

donde $Z \sim N(0,1)$. Entonces la probabilidad es igual a 0.0062097. Dado que la probabilidad de observar promedios menores o iguales a 58 es muy baja podemos decir que los datos sugieren que, efectivamente, se trata de una escuela de menor rendimiento.

APROXIMACIÓN NORMAL A LA DISTRIBUCIÓN BINOMIAL

El TCL aplica incluso para variables aleatorias discretas. En estos casos, podemos usar el TCL para aproximar la distribución de la media de la v.a. cuando el cómputo exacto puede ser complicado. Uno de los casos más útiles es el de la aproximación de la distribución binomial.

Recordemos, primero, que la distribución binomial la podemos entender como la suma de observaciones de una v.a. Bernoulli. Por lo tanto, si $Y \sim Binom(n,p)$, entonces $Y = \sum\limits_{i = 1}^n X_i$ donde $X_i \sim Bernoulli(p)$. Por lo tanto, del TCL sabemos que para $n$ suficientemente grande $\bar{X} = \frac{Y}{n} \sim N\left(p, \frac{p(1-p)}{n}\right)$ o, equivalentemente, $Y \sim N\left(np, np(1-p)\right)$.

Ejemplo: Un auditor saca una muestra de 100 vouchers para determinar qué porcentaje del total se encuentran mal documentados. ¿Cuál es la probabilidad de que más del 30% de los vouchers de la muestra hayan sido incorrectamente documentados si se sabe que la proporción de vouchers mal documentados es del 20%? Si tú fueras el auditor y observaras más del 30% de documentación incorrecta, ¿qué concluirías del supuesto de la empresa de que el 20% de vouchers son mal documentados? ¿Por qué?
R: Podemos considerar al evento de que un voucher esté mal documentado como una v.a. $Bernoulli(p = 0.2)$, asumiendo que el supuesto de la empresa es el correcto. El número de vouchers mal documentados, entonces, lo podemos modelar como una v.a. $Binomial(n = 100, p = 0.2)$, que podemos aproximar como una v.a. $N(20, 16)$.

Por lo tanto, buscamos la probabilidad:

\[ P\left[\sum X_i > 30\right] = P\left[\frac{\sum X_i - 20}{4} > \frac{30 - 20}{4}\right]\\ = P\left[\frac{\sum X_i - 20}{4} > \frac{30 - 20}{4}\right]\\ = 1 - \Phi(2.5)\\ = 0.0062 \]

Por lo tanto, si la probabilidad de un voucher estuviera mal documentado es de 0.2, el haber observado 30 vouchers mal documentados en una muestra de 100 vouchers es un evento muy poco probable.

Como auditor, entonces, la muestra da pie a cuestionar el supuesto de planteado por la empresa respecto de la proporción de vouchers mal documentados.

El uso del TCL es, desde luego, solamente una aproximación. Se considera que la aproximación será mejor en la medida de que $p \pm 3\sqrt{pq/n}$ queda en intervalo $(0,1)$, o bien, $n > 9 \left(\frac{\max(p,q)}{\min(p,q)}\right)$.

Sin embargo, en ocasiones puede ser conveniente hacer un ajuste al cálculo anterior considerando que estamos aproximando una distribución discreta mediante una distribución continua. Si consideramos que para aproximar la distribución discreta con el área bajo la curva de la distribución continua, es posible compensar parcialmente el error de estimación calculando la probabilidad para un intervalo ligeramente mayor o ligeramente menor (según corresponda) al especificado:

\[ P[X > a] \approx P[Y > \frac{a - 0.5 - \mu}{\sigma}]\\ P[X \leq a] \approx P[Y \leq \frac{a + 0.5 - \mu}{\sigma}] \]

donde $X \sim Binom(n,p)$ y $Y \sim N(0,1)$. En general es difícil determinar cuándo es apropiado utilizar la corrección por continuidad, pero típicamente la corrección funciona mejor cuando $p$ se aleja de 1/2. Ante la duda, se recomienda siempre utilizar la corrección por continuidad.

Ejemplo (cont.): Resuelve el ejercicio anterior pero aplicando la corrección por continuidad.
R: Buscamos ahora la probabilidad:

\[ P\left[\sum X_i > 29.5 \right] = P\left[\frac{\sum X_i - 20}{4} > \frac{29.5 - 20}{4}\right]\\ = P\left[\frac{\sum X_i - 20}{4} > \frac{29.5 - 20}{4}\right]\\ = 1 - \Phi(2.5)\\ = 0.00877 \]

Los tiempos de procesamiento de órdenes en el mostrador de una farmacia se distribuyen exponencialmente con una media de 10 minutos. Si 100 clientes visitan la farmacia en un periodo de 2 días, ¿cuál es la probabilidad de que al menos la mitad de los clientes necesiten esperar más de 10 minutos?

R: La probabilidad de que un cliente tenga que esperar más de 10 minutos, entonces es igual a $p = P[X > 10] = 1 - (1 - e^{-10/10}) = 0.3679$. Entonces, el evento de un cliente que espere más de 10 minutos se puede considerar una v.a. $Bernoulli(0.3679)$. Por lo tanto, la cantidad de clientes de la muestra que esperan más de 10 minutos la podemos modelar como una v.a. $Binomial(n = 100, p = 0.3679)$, entonces buscamos la probabilidad de que:

\[ P\left[\sum X_i \geq 50 \right] = P\left[\frac{\sum X_i - 36.79}{\sqrt{100*p*(1-p)}} \geq \frac{50 - 36.79}{\sqrt{100*p*(1-p)}}\right]\\ = 1 - \Phi\left[\frac{50 - 36.79}{\sqrt{100*p*(1-p)}}\right]\\ = 0.0031\\ \]

Sin embargo, considerando el ajuste a la binomial para compensar por el error del área bajo la curva:

\[ P\left[\sum X_i \geq 49.5 \right] = P\left[\frac{\sum X_i - 36.79}{\sqrt{100*p*(1-p)}} \geq \frac{49.5 - 36.79}{\sqrt{100*p*(1-p)}}\right]\\ = 1 - \Phi\left[\frac{49.5 - 36.79}{\sqrt{100*p*(1-p)}}\right]\\ = 0.0042. \]

El gerente de un supermercado quiere obtener información sobre la proporción de clientes a los que les molesta una nueva política sobre pagos. ¿A cuántos clientes necesita encuestar si quiere que la estimación de la proporción se encuentre a lo más a 0.15 de la proporción real de clientes con probabilidad de 0.98?

R: Tenemos entonces que se desea estimar la proporción $p$ que representan la proporción de clientes molestos por la nueva política (casos de éxito) en una muestra de tamaño $n$ (desconocido). Sabemos entonces que, por el TCL, si $X_i$ es la variable que captura cuando un cliente está molesto por la nueva política:

\[ P\left[\left|\frac{1}{n}\sum X_i - p \right| \leq 0.15 \right] = 0.98\\ P\left[-0.15 \leq \frac{1}{n}\sum X_i - p \leq 0.15 \right] = 0.98\\ P\left[\frac{-0.15}{\sqrt{\frac{p(1-p)}{n}}} \leq \frac{\frac{1}{n}\sum X_i - p}{\sqrt{\frac{p(1-p)}{n}}} \leq \frac{0.15}{\sqrt{\frac{p(1-p)}{n}}} \right] = 0.98\\ \]

por lo que

\[ 2.326348 = \frac{0.15}{\sqrt{\frac{p(1-p)}{n}}}\\ 2.326348\sqrt{\frac{p(1-p)}{n}} = 0.15\\ (2.326348)^2\frac{p(1-p)}{n} = (0.15)^2\\ \frac{(2.326348)^2p(1-p)}{(0.15)^2} = n \]

para la expresión anterior, $n$ alcanza su máximo, entonces, cuando $p = 0.5$. Por lo tanto, $n = 61$.

ESTIMACIÓN PUNTUAL

Lecturas recomendadas

Wackerly, Mendenhall, and Scheaffer (2008), cap. 8.1, 8.2, 8.3, 8.4
Aguirre et al. (2006), cap. 6

Dijimos en secciones anteriores que la inferencia estadística se refiere al uso de datos (en nuestro caso, obtenidos a partir de una muestra aleatoria) para sacar conclusiones respecto de características (parámetros) de una población de interés. Una tipo específico de conclusiones tiene que ver con el valor específico, puntual, que adopta ese parámetro de la población que nos interesa. Esto, de manera informal hasta aquí, llamamos un estimador puntual.

DEFINICIÓN DE ESTIMADOR

Definición: Llamamos estimador a una función o regla (fórmula) usada para calcular un valor estimado de un parámetro con base en mediciones obtenidas mediante una muestra.

Ya hemos tratado en este curso con funciones de los datos obtenidos de muestras que podríamos considerar candidatos a estimadores de características poblacionales. ¿Cuáles puedes mencionar?

Notación: Solemos usar la letra griega $\theta$ para referirnos a un parámetro poblacional de manera genérica. Para referirnos a un estimador de este parámetro poblacional se utiliza el símbolo $\hat\theta$.

Nota que, para un determinado parámetro, es posible diseñar o definir diversos posibles estimadores. Por lo tanto, se vuelve necesario contar con criterios para establecer cuándo un estimador es mejor que otro.

SESGO, ERROR DE ESTIMACIÓN, ERROR CUADRÁTICO MEDIO Y CONSISTENCIA

Una primera propiedad que nos interesará que los estimadores cumplan parte de la idea de que quisiéramos que, si tomamos varias muestras y obtenemos el estimador del parámetro, los diferentes valores que calculemos del estimador en promedio sean iguales al parámetro que se busca estimar (es decir, dado que los datos provienen de una muestra, desde luego habrá cierto error, cierta variabilidad, pero esperamos que los valores del estimador se encuentren al rededor del valor del parámetro).

Más formalmente, se busca que:

\[ E[\hat\theta] = \theta. \]

Cuando esto se cumple decimos entonces que $\hat\theta$ es un estimador insesgado de $\theta$.

Podemos entonces proponer una medida para el sesgo de un estimador:

\[ B(\hat\theta) = E[\hat\theta] - \theta. \]

Ahora bien, que un estimador sea insesgado no es la única característica que nos interesa ya que, si contamos con dos estimadores insesgados, seguramente preferiremos a aquél que genere valores del estimador que sean más cercanos al valor del parámetro. En otras palabras, preferiremos estimadores que presenten una menor dispersión respecto del valor real del parámetro que buscan estimar.

Nota que podríamos proponer como segundo criterio, entonces, escoger a aquellos estimadores cuya varianza, $Var[\hat\theta]$ sea menor. Sin embargo, recordemos que la varianza del estimador estaría definida como la variación al rededor de su media, no del parámetro a estimar (coincidirían únicamente para estimadores insesgados). Por lo tanto, es necesario proponer otro criterio. Ese criterio es el del error cuadrático medio:

\[ MSE(\hat\theta) = E[(\hat\theta - \theta)^2]. \]

Observa que:

\[ MSE(\hat\theta) = E[\hat\theta^2 - 2\hat\theta\theta + \theta^2]\\ = E[\hat\theta^2] - 2E[\hat\theta\theta] + E[\theta^2]\\ = Var[\hat\theta] + E^2[\hat\theta] - 2 \theta E[\hat\theta] + \theta^2\\ = Var[\hat\theta] + (E[\hat\theta] - \theta)^2\\ = Var[\hat\theta] + B^2(\hat\theta). \]

Ejemplo: Supongamos que $X_1, X_2, X_3$ representan valores de una muestra aleatoria obtenida de una distribución uniforme $U(a,b)$ cuya función de densidad es:

\[ f(x) = \frac{1}{b-a} \ \ \forall \ \ x \in (a,b) \]

Encuentra la media y la varianza de la distribución.

R

\[ E[x] = \int\limits_{a}^b xf(x)dx\\ = \frac{1}{b-a} \int\limits_{a}^b xdx\\ = \left.\frac{x^2}{2(b-a)}\right|_{a}^b \\ = \frac{b^2 - a^2}{2(b-a)}\\ = \frac{(b - a)(b + a)}{2(b-a)}\\ = \frac{b + a}{2}. \]

\[ Var(x) = \int\limits_{a}^b (x - \mu)^2f(x)dx\\ = \frac{1}{b-a}\int\limits_{a}^b (x - \mu)^2dx\\ = \left.\frac{(x - \mu)^3}{3(b-a)}\right|_{a}^b\\ = \frac{(b - \mu)^3 - (a - \mu)^3}{3(b-a)}\\ = \frac{(\frac{b - a}{2})^3 - (\frac{a - b}{2})^3}{3(b-a)}\\ = \frac{\frac{(b - a)^3 - (a - b)^3}{8}}{3(b-a)}\\ = \frac{(b - a)^3 - (a - b)^3}{24(b-a)}\\ = \frac{2(b - a)^3}{24(b-a)}\\ = \frac{(b - a)^2}{12}\\ \]

Considera los estimadores: $\hat\mu_1 = \frac{1}{3}X_1 + \frac{1}{3}X_2 +\frac{1}{3}X_3$, $\hat\mu_2 = \frac{1}{6}X_1 + \frac{1}{6}X_2 +\frac{2}{3}X_3$, $\hat\mu_3 = \frac{1}{3}X_1 + \frac{1}{6}X_2 +\frac{1}{6}X_3$. Calcula la media de estos estimadores, indica si son insesgados, calcula su varianza e indica cuál es el mejor estimador de la media de la distribución.

R

\[ E[\hat\mu_1] = E\left[\frac{1}{3}X_1 + \frac{1}{3}X_2 +\frac{1}{3}X_3\right]\\ = \frac{1}{3}E\left[X_1 + X_2 + X_3\right]\\ = \frac{1}{3}(E\left[X_1\right] + E\left[X_2\right] + E\left[X_3\right])\\ = \frac{b + a}{2}. \]

\[ E[\hat\mu_2] = E\left[\frac{1}{6}X_1 + \frac{1}{6}X_2 +\frac{2}{3}X_3\right]\\ = \frac{1}{6}E\left[X_1\right] + \frac{1}{6}E\left[X_2\right] + \frac{2}{3}E\left[X_3\right]\\ = \frac{b + a}{2}. \]

\[ E[\hat\mu_3] = E\left[\frac{1}{3}X_1 + \frac{1}{6}X_2 +\frac{1}{6}X_3\right]\\ = \frac{1}{3}E\left[X_1\right] + \frac{1}{6}E\left[X_2\right] + \frac{1}{6}E\left[X_3\right]\\ = \frac{2}{3}\frac{b + a}{2}\\ = \frac{b + a}{3}\\ \]

Por lo tanto, únicamente $\hat\mu_1$ y $\hat\mu_2$ son insesgados.

\[ Var[\hat\mu_1] = Var\left[\frac{1}{3}X_1 + \frac{1}{3}X_2 +\frac{1}{3}X_3\right]\\ = \frac{1}{9}Var\left[X_1 + X_2 + X_3\right]\\ = \frac{1}{3} \frac{(b - a)^2}{12}\\ = \frac{(b - a)^2}{36}. \]

\[ Var[\hat\mu_2] = Var\left[\frac{1}{6}X_1 + \frac{1}{6}X_2 +\frac{2}{3}X_3\right]\\ = \frac{1}{36}Var\left[X_1\right] + \frac{1}{36}Var\left[X_2\right] + \frac{4}{9}Var\left[X_3\right]\\ = \left(\frac{1}{36} + \frac{1}{36} + \frac{4}{9}\right) \frac{(b - a)^2}{12}\\ = \frac{18}{36} \frac{(b - a)^2}{12}\\ = \frac{(b - a)^2}{4}\\ \]

\[ Var[\hat\mu_3] = Var\left[\frac{1}{3}X_1 + \frac{1}{6}X_2 +\frac{1}{6}X_3\right]\\ = \frac{1}{9}E\left[X_1\right] + \frac{1}{36}E\left[X_2\right] + \frac{1}{36}E\left[X_3\right]\\ = \left(\frac{1}{9} + \frac{1}{36} + \frac{1}{36}\right) \frac{(b - a)^2}{12}\\ = \frac{6}{36} \frac{(b - a)^2}{12}\\ = \frac{(b - a)^2}{72}\\ \]

Claramente, el estimador 3 es el de menor varianza, sin embargo, es sesgado, ¿cómo escogemos entonces?

Para dar respuesta la pregunta anterior, recordemos que el MSE incluye tanto a la varianza como al sesgo por lo que deberíamos comparar los MSE de los estimadores 1 y 3. Para ello, utilizaremos el cociente:

\[ \frac{MSE\{\hat\mu_3\}}{MSE\{\hat\mu_1\}} = \frac{\frac{Var(X)}{6} + \frac{(b+a)^2}{36}}{\frac{Var(X)}{3}}\\ = \frac{\frac{Var(X)}{2} + \frac{(b+a)^2}{12}}{Var(X)}\\ = \frac{6Var(X) + (b+a)^2}{12Var(\mu)}\\ = \frac{1}{2} + \frac{(b+a)^2}{12Var(\mu)}\\ = \frac{1}{2} + \frac{(b+a)^2}{12\frac{(b - a)^2}{12}}\\ = \frac{1}{2} + \frac{(b+a)^2}{(b - a)^2}\\ = \frac{1}{2} + \left(\frac{b+a}{b - a}\right)^2\\ \]

\[ \frac{1}{2} + \left(\frac{b+a}{b - a}\right)^2 = 1\\ \left(\frac{b+a}{b - a}\right)^2 = \frac{1}{2}\\ \frac{b+a}{b - a} = \frac{1}{\sqrt{2}}\\ \sqrt{2}(b+a) = b - a\\ a(1 + \sqrt{2}) = b(1 - \sqrt{2})\\ \frac{a}{b} = \frac{1 - \sqrt{2}}{1 + \sqrt{2}}. \]

Entonces cuando $a/b$ es mayor que $\frac{1 - \sqrt{2}}{1 + \sqrt{2}} = -0.1716$, resulta que el estimador $\hat\mu_1$ es mejor que $\hat\mu_3$.

Al cociente que empleamos en el ejemplo anterior para comparar los MSE de los estimadores se le conoce como la eficiencia relativa.

Otra manera de determinar qué tan bueno es un estimador es fijándonos en el error de estimación:

Definición: Se llama error de estimación a $\epsilon = |\hat\theta - \theta|$

Desde luego, no podemos conocer de antemano el error de estimación pero, dado que se trata de una variable aleatoria, podemos calcular $P[\epsilon > b]$:

\[ P[\epsilon \leq b] = P[|\hat\theta - \theta| \leq b]\\ = 1 - P[|\hat\theta - \theta| > b]. \]

Si $\theta$ es un estimador insesgado, entonces podemos utilizar la desigualdad de Tchebysheff:

\[ P[\epsilon \leq b] = P[|\hat\theta - \theta| \leq b]\\ 1 - P[|\hat\theta - \theta| > a\sigma_{\hat\theta}] \geq 1 - \frac{1}{a^2}. \]

Entonces, para estimadores insesgados, podemos utilizar esta propiedad para calcular la probabilidad de que el estimador se encuentre dentro de cierta distancia respecto del valor real del parámetro.

Por último, otra propiedad que quisiéramos observar en nuestros estimadores puntuales es que, si realizamos una serie de mediciones gradualmente más grandes (incrementando $n$), quisiéramos observar no solamente que nuestro estimador proporciona de manera repetida mediciones similares, sino que adicionalmente quisiéramos que estas mediciones se acerquen cada vez más al valor real del parámetro estimado. Llamaremos a esta propiedad consistencia y podemos expresarla de manera más formal como sigue:

Definición: Decimos que un estimador $\hat\theta$ de $\theta$ es consistente si para cualquier valor positivo $\epsilon$ $\lim\limits_{n \rightarrow \infty} P[|\hat\theta - \theta| \leq \epsilon] = 1$.

Se puede probar que si $\hat\theta$ es un estimador insesgado de $\theta$, entonces $\hat\theta$ será consistente si $\lim\limits_{n \rightarrow \infty} Var[\hat\theta] = 0$. Algunas propiedades muy útiles de los estimadores consistentes son que si $\hat\theta$ y $\hat\beta$ son estimadores consistentes de $\theta$ y $\beta$, respectivamente, entonces:

$\hat\beta + \hat\theta$ es un estimador consistente de $\beta + \theta$;
$\hat\beta \times \hat\theta$ es un estimador consistente de $\beta \times \theta$;
$\hat\beta / \hat\theta$ es un estimador consistente de $\beta / \theta$, para $\hat\theta \neq 0$;
$g(\hat\theta)$ es un estimador consistente de $g(\theta)$, para $g(\cdot)$ cualquier función real continua en $\theta$.

Ejemplo

Sean $X_1 , \dots , X_n$ y $Y_1, \dots , Y_m$ muestras aleatorias independientes provenientes de poblaciones con media $\mu_1$ y $\mu_2$ y varianzas $\sigma_1^2$ y $\sigma_2^2$, respectivamente. Demuestre que $\bar{X} - \bar{Y}$ es un estimador consistente de $\mu_1 - \mu_2$.

R: Sabemos que $Var[\bar{X}] = \frac{1}{n} \sigma_1^2$ y $Var[\bar{Y}] = \frac{1}{m} \sigma_2^2$ por lo que $\bar{X}$ y $\bar{Y}$ convergen en probabilidad a $\mu_1$ y $\mu_2$, respectivamente. Por lo tanto $\bar{X} - \bar{Y}$ converge en probabilidad a $\mu_1 - \mu_2$.

Podemos aplicar la idea de consistencia mediante lo que se conoce como la Ley de los grandes números:

Teorema (LGN): Sea $X_1, X_2, \dots$ una secuencia de variables aleatorias independientes con media $\mu$ y varianza $\sigma^2 < \infty$, entonces para $\epsilon > 0$

\[ \lim\limits_{n \rightarrow \infty}P\left[\left|\frac{1}{n}(X_1 + X_2 + \dots + X_n) - \mu\right| > \epsilon\right] = 0 \]

Esto, coloquialmente, quiere decir que la media muestral se aproximará (en probabilidad) más al valor real de la media poblacional a medida que incrementamos el tamaño de la muestra.

Ejemplo: Volvamos al caso del ingreso por vivienda capturado en la ENIGH. Sabemos que la ENIGH es una muestra, pero vamos a usarla (sin conocer la media poblacional) para ilustrar la LGN. Supongamos que se desea sacar una sub-muestra de la ENIGH para estimar el ingreso capturado en ella. Veamos cómo se comporta la media del ingreso conforme aumentamos el tamaño de la muestra:

temp <- data.frame(n = 1:nrow(ingreso.vivienda))

for (i in 1:nrow(temp)){
  temp$media[i] <- mean(x = sample(x = ingreso.vivienda$Ingreso, size = i))
}

ggplot(data = temp) +
  aes(x = n, y = media) +
  geom_line() +
  geom_hline(yintercept = mean(ingreso.vivienda$Ingreso), colour = 'red') +
  theme_classic()

TEOREMA DE SLUTSKY

Teorema: Sean $X_n$ y $Y_n$ dos secuencias de variables aleatorias. Si $X_n$ converge en distribución a $X$ y $Y_n$ converge en probabilidad a una constante $c$, entonces

$X_n$ + $Y_n$ converge en distribución a $X + c$
$X_n \times Y_n$ converge en distribución a $cX$
$X_n / Y_n$ converge en distribución a $X/c$

Ejemplo: Sea $Y_1, Y_2, \dots, Y_n$ una muestra aleatoria tal que $E[Y_i] = \mu$, $E[Y_i^2] = \mu'$ y $E[Y_i^4] = \mu''$ son finitas. Demuestra que

\[ s^2 = \frac{1}{n-1} \sum\limits_{i=1}^n (Y_i - \bar{Y})^2 \]

es un estimador consistente de $\sigma^2 = Var[Y_i]$.

R: Primero, observemos que:

\[ s^2 = \frac{1}{n-1} \sum\limits_{i=1}^n (Y_i - \bar{Y})^2\\ = \frac{1}{n-1} \sum\limits_{i=1}^n (Y_i^2 - 2Y_i\bar{Y} + \bar{Y}^2)\\ = \frac{1}{n-1} \left[\sum\limits_{i=1}^n Y_i^2 - 2 \bar{Y} \sum\limits_{i=1}^n Y_i + \sum\limits_{i=1}^n \bar{Y}^2 \right]\\ = \frac{1}{n-1} \left[\sum\limits_{i=1}^n Y_i^2 - 2 n \bar{Y}^2 + n \bar{Y}^2 \right]\\ = \frac{1}{n-1} \left[\sum\limits_{i=1}^n Y_i^2 - n \bar{Y}^2\right]\\ = \frac{n}{n-1} \left[\frac{1}{n}\sum\limits_{i=1}^n Y_i^2 - \bar{Y}^2\right]. \]

Por el TCL sabemos que $\hat\beta = \frac{1}{n}\sum\limits_{i=1}^n Y_i^2$ converge en probabilidad a $\mu'$ y sabemos que $\hat\theta = \bar{Y}$ converge en probabilidad a $\mu$. Entonces, $g(\bar{Y}) = \bar{Y}^2$ converge en probabilidad a $\mu^2$, y $\hat\beta - g(\hat\theta)$ converge, entonces, a $\mu' - \mu^2 = \sigma^2$.

Finalmente, veamos que $\lim\limits_{n \rightarrow \infty} \frac{n}{n-1} = 1$. Entonces $s^2$ es un estimador consistente de $\sigma^2$.

Podemos utilizar el teorema de Slutsky para probar también el siguiente teorema:

Teorema: Supóngase que $U_n$ converge en distribución a una distribución Normal estándar y que $W_n$ converge en probabilidad a 1. Entonces $\frac{U_n}{W_n}$ converge en distribución a una distribución normal estándar.

Esto nos permite entonces ahora dar argumentos para, por ejemplo, afirmar que $\frac{\bar{X} - \mu}{\frac{s}{\sqrt{n}}}$ se distribuye aproximadamente normal estándar, sin importar la distribución que genera a los datos de la muestra con la que se calculan $\bar{X}$ y $s$.

Ejemplo

$Y \sim Binomial(n,p)$, entonces $\hat{p}_n = \frac{Y}{n}$ es un estimador insesgado de $p$. Demuestra que $\frac{\hat{p}_n - p}{\sqrt{\frac{\hat{p}_n \hat{q}_n}{n}}}$ converge en distribución a una distribución normal estándar.

R: Sabemos que $Y = \sum X_i$ donde $X_i \sim Bernoulli(p)$ por lo que $\frac{1}{n} \sum X_i = \bar{X}$ y sabemos que $E[\bar{X}] = p$ y $Var[\bar{X}] = np(1-p)$.

Por el teorema central del límite sabemos entonces que $\frac{\hat{p}_n - p}{\sqrt{\frac{p(1-p)}{n}}}$ converge en probabilidad a una $N(0,1)$.

Por otra parte, por las propiedades de convergencia en probabilidad sabemos que $\hat{p}(1 - \hat{p})$ converge en probabilidad a $p(1-p)$ por lo que $\frac{\frac{\hat{p}_n \hat{q}_n}{n}}{\frac{pq}{n}}$ converge en probabilidad a 1.

Entonces, por el teorema arriba enunciado $\frac{\hat{p}_n - p}{\sqrt{\frac{\hat{p}_n \hat{q}_n}{n}}}$ converge en distribución a una distribución normal estándar.

MÉTODO DE MÁXIMA VEROSIMILITUD

Lecturas recomendadas

Wackerly, Mendenhall, and Scheaffer (2008), cap. 9.7.

Hemos definido ya a un estimador y hemos estado analizando algunas características deseables de un estimador puntual sin embargo, ¿cómo fue que escogimos al estimador? ¿Cómo lo construimos?

A falta de un método, es posible que sea un proceso muy errático el proceso de construcción de un estimador puntual. Por lo tanto, estudiaremos ahora uno de los métodos más comúnmente empleados para construir estimadores puntuales: el método de máxima verosimilitud.

Definición (verosimilitud): Sea $\underline{X}_n = \{X_1, X_2, \dots, X_n\}$ una muestra aleatoria cuya distribución muestral depende de un vector de parámetros $\underline{\theta}$. Entonces, definimos a la verosimilitud de la muestra como la función de distribución muestral evaluada en los valores de las observaciones de la muestra y la denotamos como $L(\underline{X}_n|\underline{\theta})$, o bien, $L(\underline{\theta})$.

Observa que, dado que $\underline{\theta}$ por lo general se considera desconocido, $L$ es una función de $\underline{\theta}$.

El método de máxima verosimilitud, por lo tanto, consiste en proponer como estimador puntual del parámetro a aquel valor que maximiza la función de verosimilitud de la muestra.

¿Cómo podemos encontrar aquél valor en el cual una función toma su valor máximo?

Ejemplo: La función de probabilidad de la distribución geométrica está dada por:

\[ p(y|p) = p(1-p)^{y-1}, y = 1, 2, 3, \dots \]

Si se toma una muestra aleatoria de tamaño $n$ de una población que sigue a esta distribución, obtén el estimador de máxima verosimilitud de $p$.

R: Si $p(y|p) = p(1-p)^{y-1}, y = 1, 2, 3, \dots$, entonces:

\[ \begin{aligned} L(p|Y) &= \prod_{i = 1}^n p(1-p)^{y_i - 1}\\ &= p^n (1 - p)^{\sum y_i - n}\\ l(p|y) &= n \log(p) + \left(\sum y_i - n\right)\log(1-p)\\ \frac{d l(p|y)}{dp} &= \frac{n}{p} - \frac{\sum y_i - n}{1-p}\\ 0 &= \frac{n}{p} - \frac{\sum y_i - n}{1-p}\\ \frac{\sum y_i - n}{1-p} &= \frac{n}{p}\\ p\left(\sum y_i - n\right) &= n(1-p)\\ p\left(\sum y_i - n\right) + pn&= n\\ p\left(\sum y_i - n + n\right)&= n\\ p\sum y_i&= n\\ \hat{p} &= \frac{n}{\sum y_i}\\ \end{aligned} \]

INVARIANZA

Los EMV tienen una propiedad particularmente útil, conocida como invarianza. Si $f(\theta)$ es una función de $\theta$ y si $\hat\theta$ es el EMV de $\theta$, entonces la propiedad de invarianza de los EMVs permite afirmar que $f(\hat\theta)$ es el EMV de $f(\theta)$, es decir, $\hat{f(\theta)} = f(\hat\theta)$.

Ejemplo: ¿Cuál es el EMV de la varianza de la distribución geométrica del ejemplo anterior?
R: Sabemos que la varianza de una v.a. con distribución exponencial es igual a $\frac{1-p}{p^2}$ y, por la propiedad de invarianza de los EMVs, entonces el EMV de la varianza es igual a:

\[ \begin{aligned} \frac{1 - \hat{p}}{\hat{p}^2} &= \frac{(\sum y_i)^2 - n \sum y_i}{n^2} \end{aligned} \]

PROPIEDADES ASINTÓTICAS DE LOS EMV

Bajo ciertas condiciones, es posible demostrar que si $\hat\theta$ es el EMV de $\theta$, entonces $t(\hat\theta)$ es un estimador consistente de $t(\theta)$. Es decir que, para cualquier $\epsilon > 0$:

\[ \lim_{n \rightarrow \infty} P[|t(\hat\theta) - t(\theta)| \leq \epsilon] = 1. \]

Adicionalmente, para n grande

\[ Z = \frac{t(\hat\theta) - t(\theta)}{\sqrt{\frac{\left(\frac{dt}{d\theta}\right)^2}{nE\left[-\frac{d^2 l(\theta|X)}{d\theta^2}\right]}}} \sim N(0,1) \]

donde $l(\theta|X)$ representa la log-verosimilitud de $\theta$ para un valor de $X$.

También se puede probar que:

\[ Z = \frac{t(\hat\theta) - t(\theta)}{\left.\sqrt{\frac{\left(\frac{dt}{d\theta}\right)^2}{nE\left[-\frac{d^2 l(\theta|X)}{d\theta^2}\right]}}\right|_{\hat\theta}} \sim N(0,1) \]

Ejemplo: ¿Cuál es la varianza aproximada del EMV del parámetro $p$ de la distribución geométrica obtenido anteriormente?
R: De acuerdo con las propiedades asintóticas de los EMV, la varianza aproximada de $t(\hat\theta)$ es igual a $\frac{\left(\frac{dt}{d\theta}\right)^2}{nE\left[-\frac{d^2 l(\theta|Y)}{d\theta^2}\right]}$.

En este caso, determinamos que $l(p|Y) = n \log(p) + \left(\sum y_i - n\right)\log(1-p)$, por lo que

\[ \begin{aligned} \frac{d^2 l(p|Y)}{dp^2} &= \frac{d}{dp} \left(\frac{1}{p} - \frac{Y - 1}{1-p}\right)\\ &= \left(-\frac{1}{p^2} - \frac{Y - 1}{(1-p)^2}\right)\\ \end{aligned} \]

y, por lo tanto:

\[ \begin{aligned} nE\left[-\frac{d^2 l(p|Y)}{dp^2}\right] &= nE\left[-\left(-\frac{1}{p^2} - \frac{Y - 1}{(1-p)^2}\right)\right]\\ &= nE\left[\left(\frac{1}{p^2} + \frac{Y - 1}{(1-p)^2}\right)\right]\\ &= n\left\{\frac{1}{p^2} + E\left[\frac{Y - 1}{(1-p)^2}\right]\right\}\\ &= n\left\{\frac{1}{p^2} + \frac{E\left[Y\right] - 1}{(1-p)^2}\right\}\\ &= n\left\{\frac{1}{p^2} + \frac{1/p - 1}{(1-p)^2}\right\}\\ &= n \left\{\frac{1}{p^2} + \frac{1-p}{p(1-p)^2}\right\}\\ &= n \left\{\frac{1}{p^2} + \frac{1}{p(1-p)}\right\}\\ &= n \left\{\frac{1 - p + p}{p^2(1-p)}\right\}\\ &= \frac{n}{p^2(1-p)}\\ \end{aligned} \]

Por otra parte, $\frac{dt}{dp}$

\[ \frac{dt}{dp} = \frac{d}{dp} (p) = 1 \]

entonces:

\[ \sigma^2\{\hat{p}\} = \frac{p^2(1-p)}{n}. \]

ESTIMACIÓN POR INTERVALOS

Lecturas recomendadas

Wackerly, Mendenhall, and Scheaffer (2008), cap. 8.5 a 8.9.

Hasta ahora hemos comunicado nuestras estimaciones por medio de valores puntuales. Sin embargo, si consideramos la incertidumbre asociada a nuestra estimación, habrá quien considere proponer únicamente un estimador puntual como arriesgado. ¿Cómo podemos entonces comunicar nuestras estimaciones incorporando de alguna manera la incertidumbre o variabilidad asociada a ellas? Una manera de hacer es usando la estimación por intervalos y, en particular, los intervalos de confianza.

Antes de definir formalmente a los intervalos de confianza pensemos en lo que este tipo de estimadores representa. En primer lugar, observemos que la intensión es proponer un intervalo dentro del cual creemos que el verdadero valor del parámetro, digamos $\theta$, se encuentra. Sin embargo, recordemos que este intervalo lo vamos a construir con la información obtenida de una muestra aleatoria, por lo que sigue existiendo incertidumbre sobre el hecho de si, efectivamente, el intervalo propuesto incluye a $\theta$. Adicionalmente, observemos que dado que se trata de un intervalo, bajo condiciones fijas, preferiremos siempre intervalos lo más cerrados posible.

Nuestro objetivo será, por lo tanto, tratar de encontrar expresiones para construir intervalos lo más pequeños posibles que incluyan al verdadero valor de $\theta$ con la mayor probabilidad posible.

DEFINICIÓN DE INTERVALO DE CONFIANZA

Sea $\underline{X}_n$ una muestra aleatoria proveniente de una distribución poblacional con parámetro $\theta$. Decimos que el intervalo $\left(u(\underline{X}_n), v(\underline{X}_n)\right)$ es un intervalo de confianza para $\theta$ con coeficiente $\gamma$ si

\[ P\left[u(\underline{X}_n) \leq \theta \leq v(\underline{X}_n)\right] = \gamma. \]

Es importante tomar nota del nombre que usamos para este intervalo. Le llamamos intervalo de confianza y no de probabilidad, ¿por qué? Observemos, en primer lugar que el intervalo está planteado para $\theta$ que es un valor que, aunque es desconocido, está fijo y no lo consideramos una variable aleatoria. Adicionalmente, construimos el intervalo a partir de la regla de probabilidad asociada a $X$ pero, una vez que observamos la muestra, la probabilidad ya no aplica, simplemente el parámetro está dentro del intervalo o no lo está.

Entonces, ¿cómo debemos interpretar al intervalo de confianza? Lo podemos interpretar como la fracción de veces que un intervalo de esa naturaleza incluirá al verdadero valor del parámetro si levanto repetidamente muestras de esa distribución.

Al coeficiente de confianza $\gamma$ también se le conoce como el nivel de confianza. En ocasiones $\gamma$ es expresado alternativamente como $1 - \alpha$ o como $(1-\alpha)\%$.

Los límites inferior y superior, $u(\underline{X}_n)$ y $v(\underline{X}_n)$ respectivamente, pueden ser aleatorios ambos o puede uno de ellos estar fijo, también es posible encontrar intervalos de uno o de dos lados. Todo dependerá de las propiedades deseadas para el intervalo de confianza y, en cierta medida, de las características de la distribución.

¿Cómo podemos construir intervalos de confianza? Un primer método para su construcción se basa en la identificación de una cantidad pivotal que 1) sea una función de $\underline{X}_n$ y de $\theta$ y 2) cuya función de distribución no depende de $\theta$. De manera que, si $Y$ es una cantidad pivotal para $\theta$ podemos obtener los cuantiles $a$ y $b$ tales que $P[a \leq Y \leq b] = \gamma$ y, finalmente, despejar a $\theta$ de la desigualdad para encontrar los intervalos deseados.

Ejemplo: Se obtiene una muestra de tamaño $n = 1$ de una distribución $U(0,\theta)$. Encuentra un intervalo al 95% de confianza para $\theta$.
R: Sabemos que $X \sim U(0,\theta)$ por lo que:

\[ \begin{aligned} Y = \frac{X}{\theta} \sim U(0,1). \end{aligned} \]

Observa que $Y$ es una función de $X$ y de $\theta$ y que la distribución de $Y$ no depende de $\theta$. Entonces, usando a $Y$ como cantidad pivotal, buscamos un intervalo para $Y$ que acumule un 95% de probabilidad. Tenemos tres posibilidades: $(a,\infty)$, $(0,b)$o $(a,b)$.

Veamos la primera opción:

\[ P[Y \leq a] = 0.95\\ P[Y \leq 0.95] = 0.95\\ P[\frac{X}{\theta} \leq 0.95] = 0.95\\ P[\frac{X}{0.95} \leq \theta] = 0.95. \]

Entonces, podemos afirmar que $\theta$ se encuentra entre $(\frac{X}{0.95}, \infty)$ con un 95% de confianza.

Por otra parte, podemos también proponer el intervalo del tipo:

\[ P[Y > b] = 0.95\\ P[\frac{X}{\theta} > 0.05] = 0.95\\ P[\frac{X}{0.05} > \theta] = 0.95\\ \]

Por lo que $(0, \frac{X}{0.05})$ es otro intervalo al 95% de confianza para $\theta$.

¿Cómo construirías el intervalo del tercer tipo?

INTERVALOS DE CONFIANZA CON MUESTRAS GRANDES

Cuando estudiamos estimadores puntuales vimos que, para estimadores insesgados, para $n$ suficientemente grande podemos afirmar que:

\[ \frac{\hat\theta - \theta}{\hat\sigma^2\{\hat\theta\}} \sim N(0,1) \]

¿Qué condiciones debemos cumplir para que esto sea cierto?

y, en el caso de EMVs estudiamos una relación equivalente como parte de sus propiedades asintóticas.

Por lo tanto, asumiendo que $\hat\sigma^2\{\hat\theta\}$ depende únicamente de la muestra y de $\theta$, contamos con una expresión candidata para utilizar como cantidad pivotal cuando se cuenta con una muestra suficientemente grande.

Ejemplo: Dos marcas de refrigeradores dan garantía de 1 año para sus productos. En una muestra aleatoria de tamaño 50 de refrigeradores de la marca 1 se observó que 12 de ellos fallaron antes del año. Otra muestra aleatoria, independiente, de tamaño 60 de refrigeradores de la marca 2 mostró que también 12 de ellos fallaron antes del año. Estima la verdadera diferencia en las probabilidades de falla antes del año, $p_1$ y $p_2$ mediante un intervalo del 0.98 de confianza.
R: Sabemos que $\hat{p}_i = \frac{\sum\limits_j X_{ij}}{n_i} = \bar{X}_i$ y sabemos que

\[ \frac{(\hat{p}_1 - \hat{p}_2) - (p_1 - p_2)}{\sqrt{\frac{\hat{p}_1(1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2(1 - \hat{p}_2)}{n_2}}} \sim N(0,1) \]

para $n = n_1 + n_2$ suficientemente grande.

¿Cómo construimos entonces el intervalo al 98% de confianza? Seleccionamos un intervalo simétrico al rededor de la media de la distribución por lo que buscamos el intervalo en el que, al rededor de la media, la $N(0,1)$ acumula el 98% de probabilidad, esto es, buscamos el intervalo acotado entre los cuantiles 0.99 y 0.01. Dado que la distribución normal es simétrica, estos valores son los mismos salvo por el signo por lo que el intervalo de confianza estará dado por:

\[ (\hat{p}_1 - \hat{p}_2) \pm _{0.99}z \times \sqrt{\frac{\hat{p}_1(1 - \hat{p}_1)}{n_1} + \frac{\hat{p}_2(1 - \hat{p}_2)}{n_2}}\\ (12/50 - 12/60) \pm 2.33 \times \sqrt{(12/50 \times 38/50)/50 + (12/60 \times 48/60)/60}\\ [-0.1451, 0.2252] \]

DETERMINACIÓN DEL TAMAÑO DE MUESTRA PARA LA MEDIA O PROPORCIÓN

Un problema de especial interés para nosotros por su relevancia en la práctica es el problema de determinar el tamaño de muestra apropiado para obtener un estimador de determinadas características para la media de una población o para una proporción de individuos en una población (que, como hemos visto en ejemplos anteriores, podemos re-interpretar como una media de vv.aa. Bernoulli).

Para dar respuesta a la pregunta “¿De qué tamaño tiene que ser la muestra?”, sin caer en respuestas obvias, es necesario poner restricciones al problema. Una manera de acotar el problema es estableciendo qué tan precisa se desea que sea la estimación. Esto lo podemos establecer definiendo límites para el error de estimación.

Por lo tanto, típicamente esperaremos alguna especificación respecto del tamaño deseado para el error de estimación: $\hat\theta - \theta$. Sin embargo, recordemos que el error de estimación es desconocido, por lo que solamente podemos calcular (y en ocasiones aproximar) probabilidades asociadas al error de estimación. Así que, además de requerirse especificar un tamaño deseado para el error de estimación será necesario también especificar la probabilidad (confianza) asociada a dicha probabilidad.

De esta manera, generalmente el problema será planteado en términos de encontrar el tamaño de muestra $n$ requerido para que el error de estimación $\hat\theta - \theta$ se encuentre dentro de ciertos valores (intervalo) con una confianza de $(1-\alpha)\%$.

A partir de estas restricciones, y dependiendo de la información disponible en el planteamiento del problema, es posible obtener $n$ despejando de la desigualdad de la probabilidad utilizando 1) los cuantiles de la distribución muestral exacta o bien, 2) los cuantiles aproximados utilizando propiedades de convergencia de los estimadores.

Ejemplo: Se desea comparar la efectividad de dos métodos de entrenamiento para empleados de una planta de ensamblado. Se divide a los empleados seleccionados en dos grupos de igual tamaño, y a cada grupo se le capacita con uno de los métodos. Después de la capacitación, cada empleado realiza las tareas de ensamblado y sus tiempos son registrados. Se sabe que los tiempos de ensamblado presentan una desviación estándar de 2 minutos. Si se desea que el valor estimado de la diferencia en tiempos promedio de ensamblado no tenga un error mayor a 1 minuto con probabilidad de 0.95, ¿cuántos trabajadores deben incluirse en cada grupo?
R: Se nos pide entonces que $P[|(\hat\mu_1 - \hat\mu_2) - (\mu_1 - \mu_2)| \leq 1] = 0.95$. Por propiedades de convergencia sabemos que

\[ \begin{aligned} \frac{(\bar{X}_1 - \bar{X}_2) - (\mu_1 - \mu_2)}{\sigma\sqrt{2/n}} &\approx N(0,1)\\ P\left[_{0.025}z \leq \frac{(\bar{X}_1 - \bar{X}_2) - (\mu_1 - \mu_2)}{\sigma\sqrt{2/n}} \leq _{0.975}z\right] &\approx 0.95\\ _{0.025}z \sigma\sqrt{2/n} \leq (\bar{X}_1 - \bar{X}_2) - (\mu_1 - \mu_2) \leq _{0.975}z \sigma\sqrt{2/n}\\ 1 = _{0.975}z \sigma\sqrt{2/n}\\ \sqrt{n} = sqrt{2} _{0.975}z \sigma\\ n = 2 ( _{0.975}z \sigma)^2\\ n = 30.73. \end{aligned} \]

Entonces $n$, el tamaño de muestra de cada grupo, debe ser igual a 31.

INTERVALO DE CONFIANZA PARA LA MEDIA Y LA VARIANZA DE UNA POBLACIÓN CON DISTRIBUCIÓN NORMAL

Si se sabe (o se supone) que los datos de la muestra provienen de una población que se distribuye $N(\mu. \sigma^2)$ y deseamos construir intervalos de confianza para la media y la varianza, entonces ya contamos con las herramientas necesarias para construir estos intervalos de manera exacta y no con base en aproximaciones asintóticas, ya que ya hemos trabajado las distribuciones muestrales necesarias.

Para la media $\mu$, con $\sigma^2$ conocida, sabemos que $\bar{X} \sim N\left(\mu, \frac{\sigma^2}{n}\right)$ (como es el caso del último ejemplo de la sección anterior).
Para la media $\mu$ con $sigma^2$ desconocida, sabemos que $\frac{\bar{X} - \mu}{\frac{s}{\sqrt{n}}} \sim t_{n-1}$.
Para la diferencia de medias entre dos muestras aleatorias, sabemos que $\frac{(\bar{X}_1 - \bar{X}_2)-(\mu_!-\mu_2)}{\frac{s_p}{\sqrt{1/n_1 + 1/n_2}}} \sim t_{n_1 + n_2 - 2}$.
Para $\sigma^2$, sabemos que $\frac{(n-1)s^2}{\sigma^2} \sim \chi_{n-1}^2$.
Para el cociente $\frac{\sigma_1^2}{\sigma_2^2}$, sabemos que $\frac{\frac{s_1^2}{\sigma_1^2}}{\frac{s_2^2}{\sigma_2^2}} \sim F_{n_1 - 1, n_2 - 1}$.

Si bien podemos decir que estos son los intervalos de confianza exactos, recuerda que, por propiedades asintóticas, para valores grandes de $n$ la distribución $t$ puede ser aproximada razonablemente bien mediante la normal estándar.

Ejemplo: Se desea verificar mediante experimentos la variabilidad en la medición del volumen de un dispositivo de audio. Se levanta una muestra aleatoria de tres mediciones con los siguientes resultados: 4.1, 5.2 y 10.2. Si se supone que el volumen se distribuye normal con media $\mu$ y varianza $\sigma^2$, obtén un intervalo de confianza con un coeficiente de 0.9 para $\sigma^2$.
R: Dado que la población se distribuye $N(\mu, \sigma^2)$, sabemos entonces que

\[ \begin{aligned} \frac{(n-1)s^2}{\sigma^2} &\sim \chi_{n_1}^2\\ P\left[_{0.05} \chi_{2}^2 \leq \frac{2s^2}{\sigma^2} \leq _{0.95} \chi_{2}^2\right] &= 0.9\\ P\left[\sigma^2 _{0.05} \chi_{2}^2 \leq 2s^2 \leq \sigma^2_{0.95} \chi_{2}^2\right] &= 0.9\\ \sigma^2_{0.05} \chi_{2}^2 &\leq 2s^2\\ \sigma^2 &\leq \frac{2s^2}{_{0.05}\chi_{2}^2}\\ 2s^2 &\leq \sigma^2_{0.95} \chi_{2}^2\\ \frac{2s^2}{_{0.95} \chi_{2}^2} &\leq \sigma^2. \end{aligned} \]

Si:

sigma.2.hat <- (sd(c(4.1,5.2,10.2)))**2
q.05.csq <- qchisq(p = 0.05, df = 2)
q.95.csq <- qchisq(p = 0.95, df = 2)
a <- 2*sigma.2.hat / q.95.csq
b <- 2*sigma.2.hat / q.05.csq

$s^2 = 10.57$
$_{0.05} \chi_{2}^2 = 0.1025866$
$_{0.95} \chi_{2}^2 = 5.9914645$

Entonces, el intervalo correspondiente es $(3.5283527, 206.0698211)$.

Ejemplo: Dos muestras provenientes de una distribución normal:

Muestra 1: tamaño 5

muestra.1 <- rnorm(n = 5, mean = 15, sd = 2)
s.2.1 <- (sd(muestra.1))**2

Muestra 2: tamaño 17

muestra.2 <- rnorm(n = 17, mean = 80, sd = 3)
s.2.2 <- (sd(muestra.2))**2

Intervalo de confianza al 99% para el cociente de varianzas:

Sabemos que $\frac{\frac{s_1^2}{\sigma_1^2}}{\frac{s_2^2}{\sigma_2^2}} \sim F_{4,16}$ entonces:

b <- s.2.2*qf(p = 0.995, df1 = 4, df2 = 16)/(s.2.1)
b

## [1] 27.56726

a <- s.2.2*qf(p = 0.005, df1 = 4, df2 = 16)/(s.2.1)
a

## [1] 0.2400319

\[ b = \frac{s_2^2 \times _{0.995} F_{4,16}}{ s_1^2 }\\ a = \frac{s_2^2 \times _{0.005} F_{4,16}}{ s_1^2 }\\ \]

Entonces $\frac{\sigma_2^2}{\sigma_1^2}$ se encuentra en el intervalo $(a,b)$ con un 99% de confianza.

Nota que el verdadero cociente es 2.25.

¿Qué pasaría si aumentamos significativamente los tamaños de muestra?

Muestra 1: tamaño 101

muestra.1 <- rnorm(n = 101, mean = 15, sd = 2)
s.2.1 <- (sd(muestra.1))**2

Muestra 2: tamaño 182

muestra.2 <- rnorm(n = 182, mean = 80, sd = 3)
s.2.2 <- (sd(muestra.2))**2

Intervalo de confianza al 99% para el cociente de varianzas:

b <- s.2.2*qf(p = 0.995, df1 = 100, df2 = 181)/(s.2.1)
b

## [1] 4.250288

a <- s.2.2*qf(p = 0.005, df1 = 100, df2 = 181)/(s.2.1)
a

## [1] 1.704458

\[ b = \frac{s_2^2 \times _{0.995} F_{100,181}}{s_1^2 }\\ a = \frac{s_2^2 \times _{0.005} F_{100,181}}{s_1^2 }\\ \]

PRUEBAS DE HIPÓTESIS PARAMÉTRICAS

Lecturas recomendadas

Wackerly, Mendenhall, and Scheaffer (2008), cap. 10.

INTRODUCCIÓN

Hasta ahora hemos realizado inferencia estadística mediante estimación, ya sea puntual o por intervalos. Es decir, hemos realizado cálculos para determinar (un conjunto de) valores que consideramos como valores razonablemente factibles de un determinado parámetro poblacional desconocido. Otra forma de hacer inferencia estadística es mediante pruebas de hipótesis.

Una hipótesis (en general) es una suposición cuya validez generalmente se desea probar, típicamente como parte de un proceso de investigación. En nuestro contexto estaremos hablando, entonces, de hipótesis de carácter estadístico. Así pues, se trata de proposiciones sobre un parámetro poblacional de interés para las cuales nos interesa determinar si los datos obtenidos (mediante un procedimiento de muestreo) aportan evidencia que nos permita aceptarlas como válidas o bien, rechazarlas por falta de evidencia empírica a su favor.

Nota que las pruebas de hipótesis cobran particular relevancia en el contexto de la investigación científica ya que, precisamente, un investigador con frecuencia planteará su trabajo a partir de hipótesis de investigación que desea probar.

HIPÓTESIS NULA Y ALTERNATIVA

Ahora, si bien coloquialmente hablamos de una hipótesis, en el procedimiento estadístico de prueba o contraste de hipótesis, en realidad consideraremos dos hipótesis, mutuamente excluyentes y complementarias, a las que llamaremos hipótesis nula (generalmente denotada como $H_0$ o $H_n$) e hipótesis alternativa ($H_1$ o $H_a$). Normalmente se dice que las hipótesis típicamente serán planteadas de manera que la hipótesis en la que se centra nuestro interés teórico corresponda a la hipótesis alternativa, aunque esto no necesariamente es así.

Como veremos una vez que expliquemos el procedimiento, quizá lo más relevante sea plantear las hipótesis de manera que sepamos cuál es la distribución muestral asociada a la hipótesis nula.

Ejemplo: Un inversionista ha determinado que invertirá en un determinado fondo únicamente si los rendimientos mensuales promedio del fondo son superiores a 1.5%. Cuenta con información histórica del desempeño del fondo por lo que se plantea como hipótesis de investigación la pregunta: ¿se puede afirmar que el rendimiento mensual promedio del fondo es mayor a 1.5%?

Si deseamos realizar una prueba de hipótesis estadística, entonces la hipótesis que nos interesa probar es $\mu > 1.5$, donde $\mu$ representa al rendimiento promedio del fondo. Es decir, plantearía la prueba de hipótesis como:

\[ \begin{aligned} H_0&: \mu = 1.5\\ H_1&: \mu > 1.5. \end{aligned} \]

¿Por qué expresamos la hipótesis nula de esta manera?

Así, con la información disponible, si los datos ofrecen evidencia de que, efectivamente, los rendimientos mensuales son en promedio mayores a 1.5%, rechazaríamos la hipótesis nula. Por motivos que estudiaremos más adelante, si sucede el caso contrario, no decimos que aceptamos la hipótesis nula, simplemente que no la rechazamos.

ESTADÍSTICA DE PRUEBA Y REGIÓN DE RECHAZO

Pero, entonces, ¿cómo podemos saber si los datos ofrecen evidencia a favor o en contra de una hipótesis? En general, podemos pensar que las hipótesis estarán planteadas respecto de un parámetro poblacional, por lo que el contraste de las hipótesis, esto es, la manera en la que usaremos la información de la muestra para juzgar el apoyo o no a la hipótesis, lo haremos mediante el uso de alguna estadística, a la que en este contexto llamaremos la estadística de prueba.

Ahora bien, recordemos que el valor de la estadística de prueba, al provenir de una muestra aleatoria, se ve afectada por la aleatoridad, por la incertidumbre. Por lo tanto, es necesario establecer un criterio bajo el cual estaremos dispuestos a aceptar que la estadística respalda o no a la hipótesis. Este criterio se expresa como un rango de valores dentro de los cuales, si observamos el valor de la muestra, se considerará que no existe evidencia para apoyar a la hipótesis nula, motivo por el cual se le conoce como la región de rechazo.

Esto no parece haber resuelto en gran medida nuestro problema porque, ahora, ¿cómo decidimos cuál es la región de rechazo más apropiada? Para ello es necesario que definamos los errores a los que estamos expuestos cuando realizamos pruebas de hipótesis.

ERROR TIPO I Y ERROR TIPO II

Ahora, recordemos nuevamente que estamos trabajando con información proveniente de muestras aleatorias, por lo que sabemos que existe incertidumbre asociada al resultado de nuestra inferencia. En este caso, podemos considerar cuatro posibles resultados de nuestro proceso de prueba de hipótesis:

Que se rechace $H_0$ cuando $H_0$ es falsa;
Que no se rechace $H_0$ cuando $H_0$ es verdadera;
Que se rechace $H_0$ cuando $H_0$ es verdadera (Error Tipo I);
Que no se rechace $H_0$ cuando $H_0$ es falsa (Error Tipo II).

Las pruebas de hipótesis se dice que son como los juicios legales, la hipótesis nula es inocente hasta que se demuestre lo contrario. La carga de la prueba se encuentra del lado de la hipótesis alterna

Desde luego, nos interesa mucho no equivocarnos, por lo que medir la probabilidad de equivocarnos es un aspecto particularmente relevante de las pruebas de hipótesis. A la probabilidad de cometer un ETI se le denota comúnmente como $\alpha$ y se le conoce como el nivel de la prueba; a la probabilidad de cometer un ETII se le denota como $\beta$.

Ejemplo: Un candidato a un puesto de elección popular asegura que conseguirá más del 50% de los votos en las próximas elecciones. El equipo del candidato competidor no cree que esta afirmación sea cierta. Por lo tanto, se realiza una encuesta a 15 votantes (aleatoriamente seleccionados) y se les pregunta si votarán por el candidato. Se desea usar una prueba de hipótesis para validar la opinión del candidato sobre las preferencias de los electores. Plantea las hipótesis nula y alternativa.
R: En primer lugar, podemos interpretar la intención de voto del electorado, como una proporción de votantes, o como una probabilidad, $p$. Dado que la hipótesis que nos interesa probar consiste en la hipótesis contraria a lo que piensa el candidato, entonces:

\[ \begin{aligned} H_0 &: p = 0.5\\ H_1 &: p < 0.5. \end{aligned} \]

Si un consultor de estrategia política sugirió que la hipótesis nula se rechace si el número de personas que favorecen al candidato es menor o igual a 2, calcula el nivel de la prueba.

R: Si se nos especifica que la hipótesis nula se rechaza cuando el número de personas que favorecen al candidato sea menor o igual que 2, y se nos pide calcular el nivel de la prueba, esto es, la probabilidad de rechazar la hipótesis nula siendo que la hipótesis nula es verdadera, entonces buscamos:

\[ \alpha = P[X \leq 2 | p = 0.5, n = 15]\\ X \sim Binomial(n, p) \]

pbinom(q = 2, size = 15, prob = 0.5)

## [1] 0.003692627

Supón ahora que el candidato recibirá realmente el 30% de los votos. ¿Cuál es la probabilidad $\beta$ de que la muestre sugiera, erróneamente, que $H_0$ es verdadera?

R: Ahora buscamos determinar

\[ \begin{aligned} \beta = P[X > 2 | p = 0.3, n = 15]\\ \beta = 1 - P[X \leq 2 | p = 0.3, n = 15]\\ X \sim Binomial(n, p) \end{aligned} \]

1 - pbinom(q = 2, size = 15, prob = 0.3)

## [1] 0.8731723

Finalmente, definiremos la potencia de una prueba. La potencia es definida como la probabilidad de rechazar $H_0$ cuando $H_1$ es verdadera. Es decir, $1 - \beta$.

Volvemos ahora a nuestra pregunta anterior, ¿cómo podemos usar estas definiciones de errores para definir la región de rechazo más apropiada? Idealmente, quisiéramos realizar la prueba con los menores valores de $\alpha$ y $\beta$ posibles sin embargo esto no es posible. En la práctica, por lo tanto, se opta por 1) fijar el valor de $\alpha$ en un valor aceptablemente pequeño, 2) el valor de $\beta$ se sabe que (en la mayoría de los casos) decrecerá cuando $n$ se incrementa.

Observemos que una región de rechazo “amplia” hará, por definición, que se rechace la hipótesis nula con mayor frecuencia, incrementando $\alpha$ y reduciendo $\beta$. Por otro lado, una región de rechazo “pequeña” incrementará $\beta$, reduciendo $\alpha$. De ahí la idea de fijar uno de los valores. Dado que construimos las hipótesis de manera que la hipótesis de mayor relevancia (para nosotros) sea la hipótesis alterna, el valor que vamos a fijar como un valor aceptablemente pequeño (queremos aceptarla contundentemente). Con el valor de $\alpha$ fijo podemos determinar la región de rechazo correspondiente. Para la mayoría de las pruebas $\beta$ decrecerá, consecuentemente, con el incremento en tamaños de muestra.

Antes de estudiar los procedimientos que nos van a permitir determinar las regiones de rechazo para las pruebas, necesitamos todavía definir algo más: vamos a distinguir entre hipótesis simples e hipótesis compuestas.

Decimos que una hipótesis es simple cuando dicha hipótesis determina de manera unívoca la distribución poblacional de la que se extrae la muestra con la que se realizará el contraste.

$H_0 : \theta = \theta^*$ es una hipótesis simple.

Vamos a encontrar, dentro de las hipótesis compuestas, también:

Hipótesis de una cola: $H_1 : \theta \leq \theta^*$, $H_1 : \theta \geq \theta^*$.
Hipótesis de dos colas: $H_1 : \theta \neq \theta^*$.

Entonces, ¿cómo podemos construir intervalos de confianza razonablemente buenos? Una manera de hacerlo es utilizando el lema de Neyman-Person.

LEMA DE NEYMAN-PERSON

Supongamos que se cuenta con dos hipótesis simples (para fines de construcción, esto es, ya que claramente este tipo de planteamiento de hipótesis no se hará normalmente en la práctica). Entonces el lema de Neyman-Person se enuncia de la siguiente manera:

Lema de Neyman-Person: Supongamos que se desea contrastar las hipótesis $H_0: \theta = \theta_0$ y $H_1 : \theta = \theta_1$ a partir de una muestra aleatoria $\underline{X}_n$ a partir de una distribución poblacional $F$ con parámetro $\theta$. Entonces, para un determinado nivel de prueba $\alpha$, la prueba de máxima potencia (mínima $\beta$) en $\theta_1$ tiene una región de rechazo ($RR$) determinada por

\[ \frac{L(\theta_0)}{L(\theta_1)} < k \]

donde $k$ se obtiene de manera que la prueba tenga el nivel especificado $\alpha$.

Veamos con un ejemplo cómo podemos usar el lema de NP para determinar la región de rechazo que correspondería a la prueba más potente.

Ejemplo: Encuentra la prueba de máxima potencia para un nivel de prueba $\alpha = 0.05$ para las hipótesis $H_0: \theta = 2$ y $H_1: \theta = 1$ si se cuenta con una muestra de tamaño 1 proveniente de una población cuya función de densidad es:

\[ f(x|\theta) = \theta y^{\theta-1}, 0 < y < 1. \]

Calcula la potencia de la prueba.

R: Como se nos indica que es una muestra de tamaño 1 podemos obtener

\[ \begin{aligned} L(y|\theta = 2) &= 2y\\ L(y|\theta = 1) &= 1. \end{aligned} \]

Por lo tanto, por el lema de NP sabemos que la prueba más potente cumplirá con

\[ \begin{aligned} \frac{L(y|\theta = 2)}{L(y|\theta = 1)} &< k\\ \frac{2y}{1} &< k\\ 2y &< k\\ y &< k/2. \end{aligned} \]

Dado que $\alpha = 0.05$ entonces

\[ \begin{aligned} P[y < k/2] &= 0.05\\ \int\limits_{0}^{k/2} 2y dy &= 0.05\\ y^2|_0^{k/2} &= 0.05\\ k^2/4 &= 0.05\\ k^2 &= 0.2\\ k &= 0.4472 \end{aligned} \]

Entonces la prueba de máxima potencia para un nivel de prueba $\alpha = 0.05$ es la que tiene una RR correspondiente a $y < 0.2236$.

Ahora, para calcular la potencia de la prueba, calculamos entonces:

\[ \begin{aligned} P[y < 0.2236 | \theta = 1] &= \int\limits_{0}^{0.2236} dy\\ &= 0.2236. \end{aligned} \]

Como puedes observar del ejemplo anterior, el hecho de haber encontrado la prueba más potente, no necesariamente significa que la probabilidad de equivocarnos sea pequeña. En este caso, la probabilidad de cometer un ETII sigue siendo significativamente alta, a pesar de que la probabilidad de cometer un ETI está fijada en un nivel bajo y que la prueba es la más potente.

Adicionalmente, observa las siguientes particularidades de las pruebas de hipótesis:

La forma de la región de rechazo depende de $H_0$ y de $H_1$.
La forma de la región de rechazo depende también del valor de $\alpha$.
Para distribuciones discretas es posible que no se encuentre una prueba determinada para un valor específico de $\alpha$, en cuyo caso se utilizará el valor de la probabilidad disponible más cercana menor a $\alpha$.

Ahora bien, en la realidad no es de mucha utilidad que ambas hipótesis sean simples. ¿Cómo hacemos entonces para encontrar la región de rechazo si el lema de Neyman-Pearson no aplica a hipótesis compuestas? Pensemos, de inicio, en que tenemos una hipótesis nula simple y una hipótesis alterna compuesta de una cola, entonces (en teoría) podríamos aplicar el lema de NP a todos los valores del parámetro de interés en la hipótesis alterna. En muchas ocasiones, resulta que la región de rechazo depende únicamente del valor del parámetro establecido en la hipótesis nula. Cuando la región de rechazo obtenida en mediante el lema de NP maximiza la potencia de la prueba para cualquier valor del parámetro en la hipótesis alterna se dice que se cuenta con una prueba uniformemente más potente.

Ilustraremos estos conceptos aplicándolos a un caso en particular: las pruebas de hipótesis para la media de una distribución normal con varianza conocida.

PRUEBA DE HIPÓTESIS PARA LA MEDIA DE UNA DISTRIBUCIÓN NORMAL CON VARIANZA CONOCIDA

Ejemplo: Sea $\underline{X}_n$ una muestra aleatoria de una distribución normal $N(\mu, \sigma^2)$ con varianza conocida. Encuentra la prueba uniformemente más potente al nivel $\alpha$ para las hipótesis $H_0: \mu = \mu_0$ y $H_1: \mu > \mu_0$.
R: Siguiendo lo expuesto en la última parte de la sección anterior, primero veamos que:

\[ \begin{aligned} \frac{L(\mu_0)}{L(\mu_1)} &= \frac{\prod \left(\frac{1}{\sigma\sqrt{2\pi}}\right)e^{-\frac{(x_i - \mu_0)^2}{2\sigma^2}}}{\prod \left(\frac{1}{\sigma\sqrt{2\pi}}\right)e^{-\frac{(x_i - \mu_1)^2}{2\sigma^2}}}\\ &= \prod e^{-\frac{(x_i - \mu_0)^2 - (x_i - \mu_1)^2}{2\sigma^2}}\\ \prod e^{-\frac{(x_i - \mu_0)^2 - (x_i - \mu_1)^2}{2\sigma^2}} &< k\\ \sum{-\frac{(x_i - \mu_0)^2 - (x_i - \mu_1)^2}{2\sigma^2}} &< \ln k\\ \sum{- (x_i - \mu_0)^2 + (x_i - \mu_1)^2} &< 2\sigma^2\ln k\\ \sum{- x_i^2 + 2 x_i \mu_0 - \mu_0^2 + x_i^2 - 2x_i\mu_1 + \mu_1^2} &< 2\sigma^2\ln k\\ \sum{2 x_i \mu_0 - \mu_0^2 - 2x_i\mu_1 + \mu_1^2} &< 2\sigma^2\ln k\\ 2 n \mu_0 \bar{x} - n \mu_0^2 - 2n \mu_1 \bar{x} + n \mu_1^2 &< 2\sigma^2\ln k\\ 2 n (\mu_0 - \mu_1)\bar{x} - n (\mu_0^2 - \mu_1^2) &< 2\sigma^2\ln k\\ 2 n (\mu_0 - \mu_1)\bar{x} &< 2\sigma^2\ln k + n (\mu_0^2 - \mu_1^2) \\ \bar{x} &< \frac{2\sigma^2\ln k + n (\mu_0^2 - \mu_1^2)}{2 n (\mu_0 - \mu_1)}\\ \bar{x} &> - \frac{2\sigma^2\ln k + n (\mu_0^2 - \mu_1^2)}{2 n (\mu_1 - \mu_0)}\\ \end{aligned} \]

Tenemos entonces ya la forma de la RR. Por el momento llamemos a $-\frac{2\sigma^2\ln k + n (\mu_1^2 - \mu_0^2)}{2 n (\mu_0 - \mu_1)}$ $k'$ entonces, para una $\alpha$ determinada sabemos que:

\[ \begin{aligned} P[ \bar{x} > k' | \mu = \mu_0 ] &= \alpha \\ P[ \frac{\bar{x} - \mu_0}{\sigma/\sqrt{n}} > \frac{k' - \mu_0}{\sigma/\sqrt{n}} | \mu = \mu_0 ] &= \alpha \\ \frac{k' - \mu_0}{\sigma/\sqrt{n}} &= {_{1-\alpha}z}\\ k' &= \mu_0 + {_{1-\alpha}z} \times \sigma/\sqrt{n}. \end{aligned} \]

Es decir, que la prueba más potente (para las hipótesis simples) rechaza a la hipótesis nula siempre que $\bar{x} > \mu_0 + {_{1-\alpha}z} \times \sigma/\sqrt{n}$. Observemos ahora que el valor crítico no depende de $\mu_1$ por lo que será la prueba más potente para cualquier valor de $\mu_1$. Encontramos, entonces, la prueba uniformemente más potente.

Ejemplo: El funcionario encargado de las ventas en una empresa asegura que sus vendedores no están logrando más de 15 contactos de venta al mes (se cree que el funcionario está en un error). Para verificar esta hipótesis, se entrevistó a 36 vendedores aleatoriamente y se registró el número de contactos de venta realizados durante el mes. La media muestral de los contactos registrados fue 17 y se sabe que la varianza poblacional es igual a 9. ¿Respalda la evidencia a lo asegurado por el funcionario? Plantea la prueba de hipótesis correspondiente a un nivel del 5%.
R: La prueba de hipótesis a realizar es $H_0: \mu = 15$ vs. $H_1: \mu > 15$. Con base en los resultados obtenidos en clase, entonces, la prueba uniformemente más potente rechaza la hipótesis nula si:

\[ \begin{aligned} \bar{X} &> \mu_0 + {_{1-\alpha}z} \times \sigma / \sqrt{n}\\ \bar{X} &= 17\\ \mu_0 &= 15\\ \alpha &= 0.05\\ {_{0.95}z} &= 1.6449\\ \sigma &= 3\\ n &= 36. \end{aligned} \]

Dado que 17 es mayor que 15.8225 rechazamos la hipótesis nula, desmintiendo lo propuesto por el funcionario de la compañía.

PRUEBAS DE RAZÓN DE VEROSIMILITUDES (PRUEBA DE WILKS)

Las pruebas diseñadas utilizando el lema de NP tienen la ventaja de ser pruebas uniformemente más potentes pero tienen las limitantes de 1) aplican únicamente a hipótesis simples (aunque extendimos la idea a una hipótesis compuesta de una cola, en general no es posible encontrar pruebas uniformemente más potentes para hipótesis compuestas de dos colas) y 2) aplican únicamente a pruebas de hipótesis en las que la hipótesis nula determina completamente la distribución poblacional.

De entre los casos especiales que hemos señalado que nos interesan, tenemos casos en los que tenemos otros parámetros desconocidos, por lo que no nos es posible aplicar el lema de NP.

En estos casos, entonces, es necesario utilizar algún otro método.

Pruebas de razón de verosimilitudes: Sea $\Theta$ el vector de parámetros de una distribución poblacional y $\Omega_0$ y $\Omega_1$ los valores correspondientes a las hipótesis nula y alternativa, respectivamente ($\Omega = \Omega_0 \cup \Omega_1$ es el espacio paramétrico). Si $L(\Theta|\underline{X}_n)$ es la verosimilitud de $\Theta$, entonces es posible construir la región de rechazo para la prueba a partir del cociente

\[ \lambda = \frac{\max\limits_{\Theta \in \Omega_0} L(\Theta|\underline{X}_n)}{\max\limits_{\Theta \in \Omega} L(\Theta|\underline{X}_n)} \leq k \]

$\lambda < 1$.

En otras palabras, podemos construir la región de rechazo siguiendo una idea muy parecida a lo que hicimos cuando construimos la región de rechazo usando el lema de NP pero ahora utilizando los valores de los parámetros que maximizan la verosimilitud en cada una de las hipótesis. Al igual que en el caso anterior, $k$ se escogerá de manera que se obtenga el valor especificado (fijo) de $\alpha$. Valores de $\lambda$ cercanos a cero sugieren que los datos proporcionan evidencia a favor de la hipótesis alterna.

Desafortunadamente, la distribución muestral de la estadística que resulte del cociente de verosimilitudes no siempre va a ser conocida, lo que dificulta la determinación de la constante $k$. En estos casos, es posible utilizar el siguiente resultado asintótico bajo ciertas condiciones.

Teorema: Sea $\underline{X}_n$ una muestra aleatoria proveniente de una población con función de distribución $F_X(x|\Theta)$. Si la función de verosimilitud $L(\Theta|\underline{X}_n)$ es diferenciable respecto de $\Theta$ y la región para la que $L(\Theta|\underline{X}_n)$ no depende de parámetros desconocidos entonces, para $n$ grande, $-2\ln{\lambda} \overset{d}{\rightarrow} \chi_{r_0 - r}^2$ donde

\[ \begin{aligned} \lambda &= \frac{\max\limits_{\Theta \in \Omega_0} L(\Theta|\underline{X}_n)}{\max\limits_{\Theta \in \Omega} L(\Theta|\underline{X}_n)}\\ \end{aligned} \]

$r_0$ es el número de parámetros libres especificados en la hipótesis nula, $r$ es el número de parámetros libres especificados por la unión de las dos hipótesis y $H_0: \Theta \in \Omega_0; H_1: \Theta \in \Omega_1$ como en el caso anterior.

PRUEBA DE HIPÓTESIS PARA LA MEDIA DE UNA DISTRIBUCIÓN NORMAL CON VARIANZA DESCONOCIDA

Modificaremos entonces ahora el contexto para la prueba de la media de la distribución normal a una situación un poco más frecuente en la que se conoce que la población sigue una distribución normal, de la cual se desconocen sus parámetros, y se desea realizar construir una prueba de hipótesis uniformemente más potente para la media a partir de una muestra aleatoria.

La prueba de hipótesis entonces se construye plantea como:

\[ \begin{aligned} H_0: \mu = \mu_0\\ H_1: \mu > \mu_1 \end{aligned} \]

pero, recordemos, en este caso $\sigma^2$ es desconocida. El procedimiento de la prueba de Wilks, entonces, nos indica que, primero, debemos encontrar los valores de los parámetros que maximizan la verosimilitud.

La función de verosimilitud en este caso es

\[ L(\mu, \sigma^2|\underline{X}_n) = \prod \frac{1}{\sigma \sqrt{2\pi}}e^{-\frac{(x_i - \mu)^2}{2\sigma^2}} \]

Si recordamos el tema de EMV, necesitamos obtener $\frac{\partial L}{\partial \mu}$ y $\frac{\partial L}{\partial \sigma^2}$, aunque sabemos que, en este caso (distribución normal), es más apropiado trabajar con la log-verosimilitud, $\frac{\partial l}{\partial \mu}$ y $\frac{\partial l}{\partial \sigma^2}$.

\[ \begin{aligned} \frac{\partial l}{\partial \mu} &= \frac{1}{\sigma^2} \sum (x_i - \mu)\\ \frac{\partial l}{\partial \sigma^2} &= -\frac{n}{2\sigma^2} + \frac{1}{2\sigma^4}\sum (x_i - \mu)^2. \end{aligned} \]

Para la hipótesis nula, entonces, tenemos que $\mu = \mu_0$ por lo que

\[ \begin{aligned} \sigma_0^2 &= \frac{1}{n} \sum (X_i - \hat\mu_0)^2. \end{aligned} \]

Entonces, el cociente de verosimilitudes será:

\[ \begin{aligned} \lambda &= \frac{\max\limits_{\Theta \in \Omega_0} L(\Theta|\underline{X}_n)}{\max\limits_{\Theta \in \Omega} L(\Theta|\underline{X}_n)}\\ &= \left(\frac{\hat\sigma^2}{\hat\sigma_0^2}\right)^{n/2}\\ &= \left\{ \begin{array}{} \left[\frac{\sum (x_i - \bar{X})^2}{\sum (x_i - \mu_0)^2}\right]^{n/2} & \bar{X} > \mu_0 \\ 1 & \bar{X} \leq \mu_0 \end{array} \right.\\ \end{aligned} \]

Recordemos que $\lambda < 1$ por lo que

\[ \begin{aligned} \lambda &< k\\ \left[\frac{\sum (x_i - \bar{X})^2}{\sum (x_i - \mu_0)^2}\right]^{n/2} &< k\\ \frac{\sum (x_i - \bar{X})^2}{\sum (x_i - \mu_0)^2} &< k^{2/n}\\ \frac{\sum (x_i - \bar{X})^2}{\sum (x_i - \bar{X})^2 + n(\bar{X} - \mu_0)^2} &< k^{2/n}\\ \frac{\sum (x_i - \bar{X})^2 + n(\bar{X} - \mu_0)^2}{\sum (x_i - \bar{X})^2} &> k^{-2/n}\\ \frac{n(\bar{X} - \mu_0)^2}{\sum (x_i - \bar{X})^2} &> k^{-2/n} - 1\\ \frac{n(\bar{X} - \mu_0)^2}{\frac{1}{n-1}\sum (x_i - \bar{X})^2} &> (n-1)(k^{-2/n} - 1)\\ \frac{\sqrt{n}(\bar{X} - \mu_0)}{s} &> \sqrt{(n-1)(k^{-2/n} - 1)}\\ \frac{\sqrt{n}(\bar{X} - \mu_0)}{s} &> k'\\ \end{aligned} \]

Como sabemos que $\frac{\sqrt{n}(\bar{X} - \mu_0)}{s} \sim t_{n-1}$, entonces podemos escoger $k'$ de manera que corresponda al cuantil apropiado de la distribución $t$ para el valor especificado de $\alpha$. En otras palabras, rechazaremos la hipótesis nula cuando

\[ \frac{\sqrt{n}(\bar{X} - \mu_0)}{s} > {_{1 - \alpha}t_{n-1}} \]

Ejemplo: Un fabricante de pólvora ha desarrollado un nuevo producto que fue probado en 8 diferentes cartuchos. Se obtuvieron las siguientes velocidades (en pies por segundo):

datos <- c(3005, 2925, 2935, 2965, 2995, 3005, 2937, 2905)
datos

## [1] 3005 2925 2935 2965 2995 3005 2937 2905

Encuentra un intervalo al 95% de confianza para la velocidad media ($\mu$) para este tipo de cartuchos si se asume que las velocidades se distribuyen normales.

R: Sabemos entonces que

\[ \frac{\bar{X} - \mu}{s / \sqrt{n}} \sim t_{n-1} \]

por lo que el intervalo de confianza para $\mu$ estará dado por

\[ \bar{X} \pm {_{0.975}t_{7}} \times s/\sqrt{8} \]

x.barra <- mean(datos)
x.barra

## [1] 2959

s <- sd(datos)
s

## [1] 39.08964

x.barra + qt(p = 0.975, df = 7)*s/sqrt(8)

## [1] 2991.68

x.barra - qt(p = 0.975, df = 7)*s/sqrt(8)

## [1] 2926.32

El fabricante argumenta que la nueva pólvora produce una velocidad promedio de no menos de 3,000 pies por segundo. ¿Proporciona la muestra suficiente evidencia para contradecir lo sugerido por el fabricante? Utiliza un nivel de contraste del 0.025.

R: Se nos pide entonces, ahora, realizar una prueba de hipótesis en la que

\[ \begin{aligned} H_0 &: \mu = 3000\\ H_1 &: \mu < 3000 \end{aligned} \]

Ya que nos interesa poder verificar si la información contradice al fabricante. La prueba, entonces, es una prueba de la media de una población normal con varianza desconocida.

Si utilizamos el contraste desarrollado en los apuntes sabemos entonces que si $\bar{X} < \mu_0 + {_{\alpha}t_{7}} \times \frac{s}{\sqrt{8}}$ rechazamos la hipótesis nula. En este caso, $\bar{X} = 2,959, \mu_0 = 3,000, \alpha = 0.025, {_{\alpha}t_{7}} = -2.3646, s = 39.1$, por lo que podemos ver que $2,959 < 2967.312$. Por lo tanto, efectivamente, la evidencia sugiere que la afirmación del fabricante es falsa.

PRUEBA DE HIPÓTESIS PARA LA VARIANZA DE UNA DISTRIBUCIÓN NORMAL

Ahora nos interesa contrastar hipótesis del tipo:

\[ \begin{aligned} H_0 &: \sigma^2 = \sigma_0^2\\ H_1 &: \sigma^2 > \sigma_0^2\\ \end{aligned} \]

cuando los datos provienen de una población que se distribuye normal con media y varianza desconocidas. Nuevamente, entonces, tenemos una hipótesis nula que no determina completamente a la distribución poblacional. Por lo tanto, tenemos que utilizar la prueba de Wilks.

Siguiendo exactamente el mismo procedimiento que la sección anterior llegamos a que

\[ \begin{aligned} \lambda &= \left(\frac{\hat\sigma^2}{\sigma_0^2}\right)^{n/2}\\ \left(\frac{\hat\sigma^2}{\sigma_0^2}\right)^{n/2} &< k\\ \frac{\hat\sigma^2}{\sigma_0^2} &< k^{2/n}\\ \frac{\sum (x_i - \bar{X})^2}{\sigma_0^2} &< nk^{2/n}\\ \frac{(n-1)s^2}{\sigma_0^2} &< nk^{2/n}\\ \frac{(n-1)s^2}{\sigma_0^2} &< k' \end{aligned} \]

Sabemos entonces que, si $H_0$ es verdadera, entonces el lado izquierdo de la desigualdad sigue una distribución $\chi_{n-1}^2$ por lo que podemos encontrar $k'$ como el cuantil $(1-\alpha)\%$ que determina la región de rechazo de la prueba.

Ejemplo: Una compañía produce partes de maquinaria que se supone que no deben tener una varianza en su diámetro de más de 0.0002 $in^2$. Una muestra aleatoria de 10 partes arroja una varianza muestral de 0.0003. Asume que el diámetro de las piezas se distribuye normal.

Prueba, con un nivel de 0.05, $H_0: \sigma^2 = 0.0002$ vs $H_1: \sigma^2 > 0.0002$.

R: Como sabemos que el diámetro de las piezas sigue una distribución normal, entonces sabemos que

\[ c = \frac{(n-1)s^2}{\sigma^2} \sim \chi_{n-1}^2 \]

por lo que, para probar las hipótesis señaladas a un nivel de 0.05, podemos utilizar como criterio de contraste $c > {_{0.95}}\chi_{9}^2$. Entonces, dado que ${_{0.95}}\chi_{9}^2 = 16.92$ y $c = 9 \times 0.0003 / 0.0002 = 13.5$ no existe evidencia que permita rechazar la hipótesis nula.

PRUEBA DE HIPÓTESIS PARA LA DIFERENCIA DE MEDIAS DE POBLACIONES NORMALMENTE DISTRIBUIDAS

Ahora queremos probar hipótesis del tipo

\[ \begin{aligned} H_0 &: \mu_1 = \mu_2\\ H_1 &: \mu_1 > \mu_2 \end{aligned} \]

cuando los datos provienen de muestras independientes obtenidas de dos poblaciones distribuidas normales con media y varianza desconocida ($n_1 \neq n_2, \sigma_1^2 = \sigma_2^2 = \sigma^2$). Tenemos entonces nuevamente un caso en el que la distribución poblacional bajo la hipótesis nula no se encuentra completamente determinada.

Si seguimos un procedimiento similar a los anteriores pero ahora considerando que

\[ L(\mu_1, \mu_2, \sigma^2 | \underline{X}_1, \underline{X}_2) = \frac{1}{(\sigma^2 2 \pi)^{(n_1 + n_2)/2}} e^{-\frac{\sum (x_{i1} - \mu_1)^2 + \sum (x_{j2} - \mu_2)^2}{2\sigma^2}} \]

entonces obtendremos que la región de rechazo se puede construir mediante la desigualdad

\[ \begin{aligned} \frac{\bar{X}_1 - \bar{X}_2}{s_p\sqrt{1/n_1 + 1/n_2}} &< k;\\ \\ s_p &= \sqrt{\frac{(n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1+n_2 - 2}} \end{aligned} \]

y recordemos que la estadística de prueba se distribuye $t_{n_1 + n_2 - 2}$ por lo que podemos escoger la $k'$ apropiada para el nivel $\alpha$ especificado.

Ejemplo: Un nuevo método de capacitación para empleados de una línea de ensamblaje está siendo evaluado en una fábrica. Para ello, se formaron dos grupos de 9 empleados cada uno que fueron capacitados durante tres semanas. Un grupo fue capacitado con el método usual y el otro con el nuevo método. Se recogieron los tiempos de ensamblado de los empleados en cada grupo después de haber recibido su capacitación. Asumen que los tiempos de ensamblaje se distribuyen normales con varianza igual a 22 minutos cuadrados.

datos <- data.frame(Usual = c(32,37,35,28,41,44,35,31,34), Nuevo = c(35,31,29,25,34,40,27,32,31))

kable(x = datos, format = 'simple')

Usual	Nuevo
32	35
37	31
35	29
28	25
41	34
44	40
35	27
31	32
34	31

Realiza una prueba de hipótesis de dos colas para un nivel de prueba de 0.05. Plantea las hipótesis, la estadística de prueba y el criterio de rechazo. ¿Cuál es tu conclusión?

R: Las hipótesis a contrastar son, entonces:

\[ \begin{aligned} H_0 &: \mu_{U} = \mu{N}\\ H_1 &: \mu_{U} \neq \mu{N}. \end{aligned} \]

Sabemos que la estadística

\[ T = \frac{\mu_U - \mu_N}{\sqrt{\frac{8 \times s_U^2 + 8 \times s_N^2}{16}\sqrt{\frac{2}{9}}}} \sim t_{16} \]

cuando $H_0$ es verdadera. Por lo tanto, nuestro criterio de rechazo consistirá en rechazar $H_0$ cuando $|T| > {_{0.975}t_{16}} = 2.1199$. Dado que

$\bar{X}_U = 35.22$;
$\bar{X}_N = 31.56$;
$s_U^2 = 24.44$;
$s_N^2 = 20.03$;
$T = 1.6483$.

Entonces, dado que $|T| = 1.6483 < 2.1199$ concluimos que no existe evidencia suficiente para rechazar la hipótesis nula. Es decir, no podemos rechazar que las medias sean iguales. No queda claro, por lo tanto, que el nuevo método de capacitación logre una diferencia significativa en los tiempos de ensamblaje.

PRUEBA DE HIPÓTESIS PARA LA DIFERENCIA DE VARIANZAS DE POBLACIONES NORMALMENTE DISTRIBUIDAS

Otro caso general que nos interesa contrastar es cuando contamos con muestras provenientes de dos diferentes poblaciones (independientes) y queremos investigar si las varianzas poblaciones son iguales o no. Esto es, nos interesa contrastar las hipótesis:

\[ \begin{aligned} H_0&: \sigma^2_1 = \sigma_2^2\\ H_1&: \sigma^2_1 > \sigma_2^2\\ \end{aligned} \]

En este caso, dado que ambas poblaciones provienen de distribuciones normales, sabemos que la estadística

\[F = \frac{\frac{s_1^2}{\sigma_1^2}}{\frac{s_2^2}{\sigma_2^2}}\]

sigue una distribución $F_{n_1 - 1, n_2 - 1}$. Por lo tanto, si la hipótesis nula es verdadera, se cumple que

\[F = \frac{s_1^2}{s_2^2} \sim F_{n_1 - 1, n_2 - 1}.\]

Entonces, podemos establecer la región de rechazo como $F > k$ donde nuevamente $k$ será determinado utilizando el valor del nivel de prueba especificado.

Ejemplo (cont.): Una compañía produce partes de maquinaria que se supone que no deben tener una varianza en su diámetro de más de 0.0002 $in^2$. Una muestra aleatoria de 10 partes arroja una varianza muestral de 0.0003. Asume que el diámetro de las piezas se distribuye normal.

Se desea ahora comparar la varianza del diámetro de las piezas de la compañía con la varianza del diámetro de las piezas de un competidor. Si se tomó una muestra de 20 piezas del competidor, con una varianza muestral de 0.0001 $in^2$, ¿podemos afirmar que el competidor genera piezas con menor varianza en su diámetro? Realiza la prueba correspondiente a un nivel de 0.05.

R: Las hipótesis a contrastar ahora son, entonces, $H_0: \sigma_1^2 = \sigma_2^2$ vs $H_1: \sigma_1^2 > \sigma_2^2$. Para probar esta hipótesis utilizaremos la estadística $F = \frac{s_1^2 / \sigma_1^2}{s_2^2 / \sigma_2^2}$. Sabemos que, si $H_0$ es verdadera, $F$ se distribuye $F_{9,19}$ por lo que el criterio de rechazo, para un nivel de 0.05 será $F > {_{0.95}}F_{9,19}$.

Si $H_0$ es verdadera, $F = 0.0003/0.0001 = 3$. Dado que ${_{0.95}}F_{9,19} = 2.4227$, los datos proporcionan evidencia suficiente para afirmar que el competidor produce piezas de menor varianza en su diámetro.

PRUEBA DE HIPÓTESIS PARA UNA PROPORCIÓN (MUESTRAS PEQUEÑAS)

Volvamos ahora al caso en el que nos interesa hacer inferencia sobre la proporción de casos en una población. Típicamente, en estos casos, nos interesará contrastar las hipótesis:

\[ \begin{aligned} H_0 &: p = p_0\\ H_0 &: p < p_1\\ \end{aligned} \]

para un determinado nivel de prueba $\alpha$.

Si contamos con una muestra aleatoria $\underline{X}_n$ donde cada $X_i$ de la muestra representa un ensayo Bernoulli con valor 1 si la observación coincide con el caso que nos interesa y 0 en otro caso, entonces, podemos considerar a la estadística

\[ Y = \sum\limits_{i=1}^n X_i \]

como nuestra estadística de prueba ya que sabemos que, si $H_0$ es cierta, $Y \sim Bin(n,p_0)$. Así, nuestra región de rechazo será

\[ RR = \{Y \leq C \} \]

donde $C$ es el cuantil $\alpha$ de la distribución $Bin(n,p_0)$. Desde luego, para las hipótesis alternas $p>p_0$ y $p \neq p_0$, la región de rechazo deberá adecuarse a los cuantiles apropiados de la distribución.

Ejemplo: Una máquina en una fábrica debe ser reparada si produce más del 10% de productos defectuosos en un día. Una muestra aleatoria de 100 artículos del día contiene 15 defectuosos y el supervisor dice que la máquina debe ser reparada. ¿La muestra respalda esta decisión? Usa una prueba a un nivel de 0.01.
R: Nos interesa probar la hipótesis

\[ \begin{aligned} H_0 &: p = 0.1\\ H_0 &: p > 0.1\\ \end{aligned} \]

Sabemos que, bajo la hipótesis nula $Y = \sum\limits_{i=1}^{100} X_i \sim Bin(100, 0.1)$ por lo que, para un nivel de prueba de 0.01, la región de rechazo es

\[ Y \geq 18 \]

ya que 18 es el cuantil 99% de una $Bin(100, 0.1)$. Dado que $Y = 15$, podemos concluir que no, dada la muestra obtenida, los datos no respaldan la decisión del supervisor.

PRUEBA DE HIPÓTESIS PARA UNA PROPORCIÓN (MUESTRAS GRANDES)

Desde luego, cuando $n$ es grande, puede en ocasiones ser complicado utilizar la distribución binomial para realizar la prueba solicitada (usando tablas). Por lo tanto, es conveniente recordar que para las propociones sabemos que

\[ Y = \frac{\hat{p} - p_0}{\sqrt{\frac{p_0(1-p_0)}{n}}} \overset{d}{\rightarrow} N(0,1) \]

por lo que podemos usar el resultado asintótico de nuestra estadística de prueba para sustituir el uso de los cuantiles de la distribución binomial por los cuantiles de la distribución normal estándar en nuestra región de rechazo.

Ejemplo: ¿Cambia el resultado si utilizas el hecho de que la muestra es una muestra grande?
R: Dado que el número de observaciones es suficientemente grande, podemos entonces utilizar el hecho de que

\[ Y = \frac{\hat{p} - p_0}{\sqrt{\frac{p_0(1-p_0)}{n}}} \overset{d}{\rightarrow} N(0,1) \]

por lo que la región de rechazo corresponderá a

\[ Y \geq {_{0.99}z} \]

Dado que $Y = \frac{0.15 - 0.1}{\sqrt{\frac{0.1(1-0.1)}{100}}} = 1.667$ y ${_{0.99}z} = 2.33$ entonces, concluimos nuevamente que los datos no aportan evidencia a favor de la opinión del supervisor.

PRUEBA DE HIPÓTESIS PARA MUESTRAS GRANDES

En todos los casos anteriores hemos sido capaces de construir la prueba de hipótesis a partir de la distribución muestral de una estadística de prueba. ¿Qué sucede si no conocemos la distribución muestral pero contamos con una muestra suficientemente grande?

En estos casos, si para el parámetro $\theta$ conocemos una estadística (un estimador) $\hat\theta$ para el cual conocemos su distribución asintótica, podemos entonces hacer uso de esta distribución para construir la estadística de prueba y, por lo tanto, la región de rechazo.

La media poblacional (distribución desconocida);
La proporción;
La diferencia de medias poblacionales (distribución desconocida);
La diferencia de proporciones.

Todos estos, son casos que en el pasado vimos que, estandarizados, asintóticamente se distribuyen normales, por lo que podemos hacer uso de esta propiedad para realizar la prueba de hipótesis correspondiente.

EN RESUMEN: PASOS PARA REALIZAR UN CONTRASTE DE HIPÓTESIS

En términos generales podríamos decir que las pruebas anteriormente estudiadas se pueden sintetizar en los siguientes pasos:

Plantear las hipótesis nula y alternativa correspondientes.
Establecer el nivel de la prueba, es decir, la probabilidad de ETI que estamos dispuestos a aceptar.
Determinar (y calcular el valor de) una estadística de prueba para la cual es conocida en forma determinada su distribución muestral si la hipótesis nula es cierta.
Establecer la región de rechazo mediante la distribución muestral de la estadística de prueba y el nivel de la prueba pre-especificado.
Verificar si la estadística de prueba se encuentra en la región de rechazo.

VALIDACIÓN DE LA HIPÓTESIS DE NORMALIDAD

Lecturas sugeridas

Aguirre et al. (2006), cap. 9.1.

En esta sección vamos a estudiar cómo podemos contrastar una hipótesis muy particular. El supuesto de normalidad en los datos, como hemos visto a lo largo de este curso, puede ser muy frecuente. Entonces, con frecuencia también nos interesará preguntarnos si, una vez observados los datos, podemos afirmar que dichos datos respaldan la hipótesis de normalidad planteada como un supuesto.

En términos generales, podríamos decir que las pruebas siguen un planteamiento general similar, es decir, todas las pruebas buscan comparar características de los datos con características de una distribución normal y, con base en ello, determinar si los datos parecen ajustarse a esas características (a esa descripción) de la distribución normal.

ANÁLISIS GRÁFICO

Veremos primero, no tanto una prueba formalmente, sino una herramienta de exploración que nos puede ayudar en nuestro contraste de hipótesis. Para ello, haremos uso de los datos de un ejemplo.

Ejemplo: Se cuenta con 1,239 observaciones del precio de cierre del NAFTRAC y se desea comprobar si el precio y el rendimiento se comportan como una variable aleatoria normal.

datos <- read.xlsx(xlsxFile = './datos/Datos históricos NAFTRACISHRS.xlsx')
datos$Fecha <- as.Date(datos$Fecha, format = '%d.%m.%Y')

ggplot(data = datos) +
  aes(x = Fecha, y = Cierre) +
  geom_line() +
  theme_classic()

ggplot(data = datos[complete.cases(datos),]) +
  aes(x = Fecha, y = Variacion) +
  geom_line() +
  theme_classic()

Gráficamente, ¿qué podemos hacer? Una primera podría ser trazar su histograma:

ggplot(data = datos) +
  aes(x = Cierre) +
  geom_histogram(bins = 30) +
  theme_classic()

ggplot(data = datos[complete.cases(datos),]) +
  aes(x = Variacion) +
  geom_histogram(bins = 30) +
  theme_classic()

Pero, ¿cómo podemos decir con algo de certeza si esto se parece a una normal o no? Podemos tener algunas directrices pero en muchos casos no va a ser muy claro.

La herramienta más usada se basa en la frecuencia acumulada de las observaciones, en los cuantiles de la distribución. La idea básicamente consiste en comparar, no la densidad (que es lo que aproxima el histograma) sino la distribución.

Cierre <- aggregate(x = list(Frecuencia = datos$Cierre), by = list(Precio = datos$Cierre), FUN = length)
Cierre$FrecuenciaRel <- Cierre$Frecuencia / sum(Cierre$Frecuencia)

ecdf.cierre <- ecdf(x = datos$Cierre)

Cierre$FrecuenciaRelAc <- ecdf.cierre(Cierre$Precio)

kable(x = head(Cierre), format = 'simple')

Precio	Frecuencia	FrecuenciaRel	FrecuenciaRelAc
33.00	1	0.0008071	0.0008071
33.12	1	0.0008071	0.0016142
33.63	1	0.0008071	0.0024213
33.76	1	0.0008071	0.0032284
33.78	1	0.0008071	0.0040355
33.85	1	0.0008071	0.0048426

Rendimiento <- aggregate(x = list(Frecuencia = datos$Variacion[complete.cases(datos)]), by = list(Tasa = datos$Variacion[complete.cases(datos)]), FUN = length)

Rendimiento$FrecuenciaRel <- Rendimiento$Frecuencia / sum(Rendimiento$Frecuencia)

ecdf.rendimiento <- ecdf(x = datos$Variacion[complete.cases(datos)])

Rendimiento$FrecuenciaRelAc <- ecdf.rendimiento(Rendimiento$Tasa)

kable(x = head(Rendimiento), format = 'simple')

Tasa	Frecuencia	FrecuenciaRel	FrecuenciaRelAc
-0.0630892	1	0.0008078	0.0008078
-0.0606005	1	0.0008078	0.0016155
-0.0540465	1	0.0008078	0.0024233
-0.0434221	1	0.0008078	0.0032310
-0.0387534	1	0.0008078	0.0040388
-0.0379184	1	0.0008078	0.0048465

Y ahora, la comparamos contra el valor del cuantil correspondiente para una distribución normal. Claro, tendremos un problema porque no conocemos los parámetros de la distribución normal, por lo que generalmente se comparan contra los datos estandarizados:

ggplot(data = datos) +
  aes(sample = Cierre) +
  geom_qq() +
  geom_qq_line() +
  theme_classic()

datos$CierreStd <- scale(x = datos$Cierre)

ggplot(data = datos) +
  aes(sample = CierreStd) +
  geom_qq() +
  geom_qq_line() +
  theme_classic()

ggplot(data = datos[complete.cases(datos),]) +
  aes(sample = Variacion) +
  geom_qq() +
  geom_qq_line() +
  theme_classic()

datos$VariacionStd <- scale(x = datos$Variacion)

ggplot(data = datos[complete.cases(datos),]) +
  aes(sample = VariacionStd) +
  geom_qq() +
  geom_qq_line() +
  theme_classic()

Como podemos ver en el ejemplo, ambos casos parecen de alguna manera no ajustarse a la distribución normal. ¿Cómo podemos caracterizar estas desviaciones?

A manera de ayuda utilizaremos los siguientes gráficos:

PRUEBA DE KOLMOGOROV-SMIRNOFF

La prueba de KS está basada en la idea (similar al QQ plot) de comparar los cuantiles de dos muestras de datos para verificar si existe evidencia de que provienen de la misma distribución (en nuestro caso, queremos entonces comparar contra una distribución normal).

La estadística de prueba propuesta por KS es:

\[ D = \sup\limits_x |F_n(x) - F(x)| \]

donde

\[ F_n(x) = \frac{\text{número de observaciones menores o iguales a x}}{n} \]

y $F(x)$ es el valor de la función de distribución de comparación (teórica o empírica) evaluada en $x$. Bajo la hipótesis nula de igualdad de distribuciones, $D$ converge a 0.

El problema de la prueba de KS para los residuales de regresión es que está formulada para el caso en el que se conocen los parámetros de la distribución teórica de comparación. Dado que se están haciendo estimaciones de estos parámetros, la prueba estrictamente hablando no aplica.

En 1967, Hubert W, Lilliefors propuso una modificación a la prueba KS para solventar este problema (usando, nuevamente, simulación de Monte Carlo).

Ejemplo (cont.): Continuando con nuestro ejemplo, si queremos ahora realizar una prueba más formal y no basarnos solamente en el análisis gráfico de los datos, entonces podríamos utilizar la prueba de KS:

ks.test(x = datos$CierreStd, y = 'pnorm')

## Warning in ks.test.default(x = datos$CierreStd, y = "pnorm"): ties should not
## be present for the one-sample Kolmogorov-Smirnov test

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  datos$CierreStd
## D = 0.16485, p-value < 2.2e-16
## alternative hypothesis: two-sided

ks.test(x = datos$VariacionStd, y = 'pnorm')

## Warning in ks.test.default(x = datos$VariacionStd, y = "pnorm"): ties should
## not be present for the one-sample Kolmogorov-Smirnov test

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  datos$VariacionStd
## D = 0.043959, p-value = 0.01671
## alternative hypothesis: two-sided

PRUEBA DE JARQUE-BERA

Alternativamente, podemos utilizar la prueba de Jarque-Bera, que es una prueba cuya idea se basa en comparar otras características de la distribución normal: su coeficiente de asimetría y su coeficiente de kurtosis. Así, se plantea la hipótesis nula de que el coeficiente de asimetría de la distribución poblacional es igual a 0 y su excedente de kurtosis igual a 3 (los valores de referencia para la distribución normal). Bajo la hipótesis nula, entonces, la estadística de prueba:

\[ JB = \frac{n}{6}\left[A^2 + \frac{1}{4}(K - 3)^2\right] \overset{d}\rightarrow \chi_2^2 \]

donde

$A = \frac{\frac{1}{n}\sum (x_i - \bar{x})^3}{\left[\frac{1}{n}\sum (x - \bar{x})^2\right]^{3/2}}$;
$K = \frac{\frac{1}{n}\sum (x_i - \bar{x})^4}{\left[\frac{1}{n}\sum (x - \bar{x})^2\right]^{2}}$.

Recuerda: la prueba de JB es una prueba asintótica por lo que, en principio, requiere de muestras grandes.

Ejemplo (cont.): Si aplicamos ahora la prueba de JB a nuestros datos, para un nivel de prueba de 0.01:

Para el precio de cierre:

n <- length(datos$Cierre)
A <- (sum((datos$Cierre - mean(datos$Cierre))^3)/n)/(sd(datos$Cierre)^(3))
K <- (sum((datos$Cierre - mean(datos$Cierre))^4)/n)/(sd(datos$Cierre)^(4))
JB <- n*(A^2 + ((K - 3)^2)/4)/6

JB

## [1] 206.9423

qchisq(p = 0.99, df = 2)

## [1] 9.21034

pchisq(q = JB, df = 2)

## [1] 1

Para el rendimiento:

n <- length(datos$Variacion[complete.cases(datos)])
A <- (sum((datos$Variacion[complete.cases(datos)] - mean(datos$Variacion[complete.cases(datos)]))^3)/n)/(sd(datos$Variacion[complete.cases(datos)])^3)
K <- (sum((datos$Variacion[complete.cases(datos)] - mean(datos$Variacion[complete.cases(datos)]))^4)/n)/(sd(datos$Variacion[complete.cases(datos)])^(4))
JB <- n*(A^2 + ((K - 3)^2)/4)/6

JB

## [1] 421.0742

qchisq(p = 0.99, df = 2)

## [1] 9.21034

pchisq(q = JB, df = 2)

## [1] 1

¿Qué notas de peculiar en las pruebas de hipótesis de normalidad si las comparas con las otras pruebas de hipótesis que estuvimos haciendo en este curso?

COMENTARIOS FINALES

Si bien es comprensible que nuestro interés sea absorbido por los temas “de moda”, la realidad es que las funciones y preocupaciones tradicionales del analista de datos son esencialmente las mismas: la cantidad de datos no sustituye a la calidad, la diversidad de fuentes de información no sustituye a un buen esquema de muestreo, el uso de técnicas sofisticadas no garantiza el cumplimiento de supuestos teóricos.

Maindonald and Braun (2010): “Lo mejor que un buen análisis puede hacer es resaltar la información contenida en los datos. Ninguna tecnología estadística o computacional puede sustituir al buen diseño de acopio de datos, al entendimiento del contexto en el que los datos serán utilizados o la habilidad en el uso de las metodologías de análisis estadístico.”

ANEXOS

DESIGUALDAD DE TCHEBYSHEFF

Sea $Y$ una variable aleatoria con media $\mu$ y varianza $\sigma^2 < \infty$. Entonces, para toda $\alpha > 0$

\[ P[|Y - \mu| \geq \alpha] \leq \frac{\sigma^2}{\alpha^2} \]

Demostración: Sea $X = (Y - \mu)^2$, entonces, usando la desigualdad de Markov:

\[ P[|Y - \mu| \geq \alpha] = P[X \geq \alpha^2]\\ \leq \frac{E[X^2]}{\alpha^2}\\ = \frac{\sigma^2}{\alpha^2} \]

REFERENCIAS

Aguirre, Víctor, Alejandro Alegría Hernández, Begoña Artaloitia, Beatriz Balmaseda, Juan José Fernández Durán, Graciela Garza Ruiz Esparza, Víctor Manuel Guerrero, et al. 2006. Fundamentos de Probabilidad y Estadística. Jit Press.

DeGroot, Morris H. 1988. Probabilidad y Estadística. 2nd ed. Addison-Wesley Iberoamericana.

Instituto Nacional de Estadística y Geografía. 2022a. “Encuesta Nacional de Ingresos y Gastos de Los Hogares (ENIGH) 2022: Cuestionario de Hogares y Vivienda.” Instituto Nacional de Estadística y Geografía.

———. 2022b. “Encuesta Nacional de Ingresos y Gastos de Los Hogares (ENIGH) 2022: Nueva Serie, Diseño Conceptual.” Instituto Nacional de Estadística y Geografía.

Maindonald, John, and W. John Braun. 2010. Data Analysis and Graphics Using r. 3rd ed. Cambridge University Press.

Wackerly, Dennis D, William Mendenhall, and Richard L Scheaffer. 2008. Mathematical Statistics with Applications. Thomson Brooks/Cole Belmont, CA.

Inferir se refiere al acto de sacar conclusiones o deducciones a partir de otra cosa (la información o los datos, en nuestro caso).↩︎
Se considera únicamente los tipos de variables categóricas y cuantitativas, no los subtipos, ya que no dan lugar a casos especiales de análisis.↩︎

\(\sum X_i\)	\(\frac{\sum X_i}{3}\)	p	P
0	0	\(\frac{1}{8}\)	\(\frac{1}{8}\)
1	1/3	\(\frac{3}{8}\)	\(\frac{4}{8}\)
2	2/3	\(\frac{3}{8}\)	\(\frac{7}{8}\)
3	1	\(\frac{1}{8}\)	\(\frac{8}{8}\)

INFERENCIA ESTADÍSTICA

Apuntes de estudio

Enrique Cuervo Guzmán

2024-12-04

PRELIMINARES

EVALUACIÓN

CRITERIOS Y POLÍTICAS

BIBLIOGRAFÍA BÁSICA

INTRODUCCIÓN

ESTADÍSTICA

POBLACIONES, MUESTRAS, PARÁMETROS Y ESTADÍSTICOS

MUESTREO ALEATORIO SIMPLE

ESTADÍSTICA DESCRIPTIVA

TIPOS DE VARIABLES Y ESCALAS DE MEDICIÓN

AED UNIVARIADO

DISTRIBUCIONES DE FRECUENCIA

MEDIDAS DESCRIPTIVAS

MEDIDAS DE POSICIÓN

MEDIDAS DE DISPERSIÓN

DIAGRAMAS DE CAJA Y BRAZOS

AED BI-VARIADO

CATEGÓRICA X CATEGÓRICA

CATEGÓRICA X CUANTITATIVA

CUANTITATIVA X CUANTITATIVA

AED - CONCLUSIONES

ESTADÍSTICA INFERENCIAL

DISTRIBUCIONES MUESTRALES Y TEOREMA CENTRAL DEL LÍMITE

DISTRIBUCIONES MUESTRALES RELACIONADAS CON UNA POBLACIÓN NORMALMENTE DISTRIBUIDA

\(\chi^2\)

\(t\)

\(F\)

TEOREMA CENTRAL DEL LÍMITE

APROXIMACIÓN NORMAL A LA DISTRIBUCIÓN BINOMIAL

ESTIMACIÓN PUNTUAL

DEFINICIÓN DE ESTIMADOR

SESGO, ERROR DE ESTIMACIÓN, ERROR CUADRÁTICO MEDIO Y CONSISTENCIA

TEOREMA DE SLUTSKY

MÉTODO DE MÁXIMA VEROSIMILITUD

INVARIANZA

PROPIEDADES ASINTÓTICAS DE LOS EMV

ESTIMACIÓN POR INTERVALOS

DEFINICIÓN DE INTERVALO DE CONFIANZA

INTERVALOS DE CONFIANZA CON MUESTRAS GRANDES

DETERMINACIÓN DEL TAMAÑO DE MUESTRA PARA LA MEDIA O PROPORCIÓN

INTERVALO DE CONFIANZA PARA LA MEDIA Y LA VARIANZA DE UNA POBLACIÓN CON DISTRIBUCIÓN NORMAL

PRUEBAS DE HIPÓTESIS PARAMÉTRICAS

INTRODUCCIÓN

HIPÓTESIS NULA Y ALTERNATIVA

ESTADÍSTICA DE PRUEBA Y REGIÓN DE RECHAZO

ERROR TIPO I Y ERROR TIPO II

LEMA DE NEYMAN-PERSON

PRUEBA DE HIPÓTESIS PARA LA MEDIA DE UNA DISTRIBUCIÓN NORMAL CON VARIANZA CONOCIDA

PRUEBAS DE RAZÓN DE VEROSIMILITUDES (PRUEBA DE WILKS)

PRUEBA DE HIPÓTESIS PARA LA MEDIA DE UNA DISTRIBUCIÓN NORMAL CON VARIANZA DESCONOCIDA

PRUEBA DE HIPÓTESIS PARA LA VARIANZA DE UNA DISTRIBUCIÓN NORMAL

PRUEBA DE HIPÓTESIS PARA LA DIFERENCIA DE MEDIAS DE POBLACIONES NORMALMENTE DISTRIBUIDAS

PRUEBA DE HIPÓTESIS PARA LA DIFERENCIA DE VARIANZAS DE POBLACIONES NORMALMENTE DISTRIBUIDAS

PRUEBA DE HIPÓTESIS PARA UNA PROPORCIÓN (MUESTRAS PEQUEÑAS)

PRUEBA DE HIPÓTESIS PARA UNA PROPORCIÓN (MUESTRAS GRANDES)

PRUEBA DE HIPÓTESIS PARA MUESTRAS GRANDES

EN RESUMEN: PASOS PARA REALIZAR UN CONTRASTE DE HIPÓTESIS

VALIDACIÓN DE LA HIPÓTESIS DE NORMALIDAD

ANÁLISIS GRÁFICO

PRUEBA DE KOLMOGOROV-SMIRNOFF

PRUEBA DE JARQUE-BERA

COMENTARIOS FINALES

ANEXOS

DESIGUALDAD DE TCHEBYSHEFF

REFERENCIAS