Librerías y opciones.

# Librerias y opciones
library(pacman)
p_load(moderndive, # Para hacer el muestreo
       tidyverse)  # Para usar la %>% y los verbos dplyr

# Para mostrar más números decimales en las tibbles.
options(pillar.sigfig = 7)

Ejercicio 1

En clase analizamos si la media y la mediana eran estimadores insesgados del promedio de una población normal.

En este ejercicio harás lo mismo, pero para una distribución beta. En esta (p. 906) y esta otra¹ ligas puedes aprender más sobre esta distribución.

1. Usa el siguiente código para construir la distribución:

set.seed(1)
datos <- tibble(valores = 
                  rbeta(n = 1e6, 
                        shape1 = 6, 
                        shape2 = 1))

2. Grafica la distribución. ¿Es simétrica? Si la respuesta es no, ¿qué tipo de sesgo tiene?

datos %>%
  ggplot(aes(valores)) +
  geom_density()

# Sesgo izquierdo

3. Calcula la media y mediana poblacional

est_pob <-
  datos %>%
  summarise(media_pob = mean(valores),
            mediana_pob = median(valores))

est_pob

## # A tibble: 1 x 2
##   media_pob mediana_pob
##       <dbl>       <dbl>
## 1 0.8572228   0.8910398

set.seed(1)
datos %>%
  rep_sample_n(size = 5, reps = 1000) %>%
  summarise(media_muestral = mean(valores),
            mediana_muestral = median(valores)) %>%
  summarise(promedio_media_muestral =
              mean(media_muestral),
            promedio_mediana_muestral =
              mean(mediana_muestral)) %>%
  gather(key = variable,
         value = valor)

## # A tibble: 2 x 2
##   variable                      valor
##   <chr>                         <dbl>
## 1 promedio_media_muestral   0.8576992
## 2 promedio_mediana_muestral 0.8820205

set.seed(1)
datos %>%
  rep_sample_n(size = 5, reps = 1000) %>%
  summarise(media_muestral = mean(valores),
            mediana_muestral = median(valores)) %>%
  summarise(promedio_media_muestral =
              mean(media_muestral),
            promedio_mediana_muestral =
              mean(mediana_muestral)) %>%
  gather(key = variable,
         value = valor) %>%
  mutate(dif_absoluta_media = 
           ifelse(variable == "promedio_media_muestral",
                  yes = abs(est_pob$media_pob - valor),
                  no = abs(est_pob$mediana_pob - valor)))

## # A tibble: 2 x 3
##   variable                      valor dif_absoluta_media
##   <chr>                         <dbl>              <dbl>
## 1 promedio_media_muestral   0.8576992       0.0004763692
## 2 promedio_mediana_muestral 0.8820205       0.009019281

set.seed(1)
datos %>%
  rep_sample_n(size = 50, reps = 1000) %>%
  summarise(media_muestral = mean(valores),
            mediana_muestral = median(valores)) %>%
  summarise(promedio_media_muestral = mean(media_muestral),
            promedio_mediana_muestral = mean(mediana_muestral)) %>%
  gather(key = variable,
         value = valor) %>%
  mutate(dif_absoluta_media = 
           ifelse(variable == "promedio_media_muestral", 
                  yes = abs(est_pob$media_pob - valor), 
                  no = abs(est_pob$mediana_pob - valor)))

## # A tibble: 2 x 3
##   variable                      valor dif_absoluta_media
##   <chr>                         <dbl>              <dbl>
## 1 promedio_media_muestral   0.8573599       0.0001370857
## 2 promedio_mediana_muestral 0.8904710       0.0005687476

set.seed(1)
datos %>%
  rep_sample_n(size = 100, reps = 1000) %>%
  summarise(media_muestral = mean(valores),
            mediana_muestral = median(valores)) %>%
  summarise(promedio_media_muestral = mean(media_muestral),
            promedio_mediana_muestral = mean(mediana_muestral)) %>%
  gather(key = variable,
         value = valor) %>%
  mutate(dif_absoluta_media = 
           ifelse(variable == "promedio_media_muestral", 
                  yes = abs(est_pob$media_pob - valor), 
                  no = abs(est_pob$mediana_pob - valor)))

## # A tibble: 2 x 3
##   variable                      valor dif_absoluta_media
##   <chr>                         <dbl>              <dbl>
## 1 promedio_media_muestral   0.8572043      0.00001850602
## 2 promedio_mediana_muestral 0.8902756      0.0007641925

set.seed(1)
datos %>%
  rep_sample_n(size = 5, reps = 10000) %>%
  summarise(media_muestral = mean(valores),
            mediana_muestral = median(valores)) %>%
  summarise(promedio_media_muestral = mean(media_muestral),
            promedio_mediana_muestral = mean(mediana_muestral)) %>%
  gather(key = variable,
         value = valor) %>%
  mutate(dif_absoluta_media = 
           ifelse(variable == "promedio_media_muestral", 
                  yes = abs(est_pob$media_pob - valor), 
                  no = abs(est_pob$mediana_pob - valor)))

## # A tibble: 2 x 3
##   variable                      valor dif_absoluta_media
##   <chr>                         <dbl>              <dbl>
## 1 promedio_media_muestral   0.8573634       0.0001406020
## 2 promedio_mediana_muestral 0.8802809       0.01075889

Ejercicio 2

A continuación verás el código que utilicé en clase para ilustrar qué significa el nivel de confianza de un intervalo de confianza. Por favor cópialo y ejecútalo en el script de tu tarea:

# Paso 1 - Generar los datos poblacionales
set.seed(33)
datos <-
  tibble(valores = rnorm(n = 1e6, mean = 0, sd = 1))

# Paso 2 - Definir el tamaño de cada muestra y el número de muestras
tamaño_muestra <- 1000
numero_muestras <- 100

# Paso 3 - Calcular el error estándar de la media muestral, asumiendo
# que el tamaño de la muestra es 1,000
error_est_con_datos_pob <- sd(datos$valores)/sqrt(tamaño_muestra)

# Paso 4 - Calcular la media de cada muestra y construye su respectivo
# intervalo de confianza usando el error estándar que calculamos
# en el paso 3
set.seed(4)
datos %>%
  rep_sample_n(size = tamaño_muestra,
               replace = T,
               reps = numero_muestras) %>%
  summarise(media_muestral = mean(valores),
            lim_sup = media_muestral + 1.96*error_est_con_datos_pob,
            lim_inf = media_muestral - 1.96*error_est_con_datos_pob) %>%  
  mutate(color_intervalos = ifelse(lim_sup < 0 | lim_inf > 0,
                                   yes = "salmon", 
                                   no = "steelblue")) %>%
  ggplot(aes(x = replicate, 
             y = media_muestral, 
             color = color_intervalos)) +
  geom_hline(yintercept = 0, 
             color = "grey40", 
             linetype = 2, 
             size = 0.5) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = lim_inf, 
                    ymax = lim_sup), 
                width = 0.01,
                alpha = 0.82, 
                size = 1) +
  scale_color_manual(values = 
                       c("salmon", "steelblue")) +
  labs(x = "\nNúmero de simulación") +
  theme_minimal() +
  theme(axis.title.x = 
          element_text(hjust = 1, face = "bold",
                       color = "grey30", size = 15),
        axis.title.y = element_blank(),
        axis.text = element_text(size = 15),
        legend.position = "none")

1. Por favor explica qué es lo que estoy haciendo en cada paso. Para el paso 3 basta con que expliques hasta el renglón en donde uso scale_color_manual().

Paso 1 Generar los datos poblacionales.

# Paso 1 - Generar los datos poblacionales
set.seed(33) # Fijamos semilla para que todo sea reproducible

# Generamos el tibble poblacional, que consiste en una variable
## aleatoria que distribuye de manera normal con parámetros 0
## y desviación estándar de 1
datos <-
  tibble(valores = 
           rnorm(n = 1e6, mean = 0, sd = 1))

Paso 2 Definir el tamaño de cada muestra y el número de muestras.

# Paso 2 - Definir el tamaño de cada muestra y el número de muestras
tamaño_muestra <- 1000 # Observaciones por muestra aleatoria
numero_muestras <- 100 # Numero de muestras

Paso 3 Calcular el error estándar de la media muestral, asumiendo que el tamaño de la muestra es 1,000.

# Paso 3 - Calcular el error estándar de la media muestral, asumiendo
# que el tamaño de la muestra es 1,000

# Calculamos el error estándar de la muestra con la formula normal.
error_est_con_datos_pob <- sd(datos$valores)/sqrt(tamaño_muestra)

Paso 4 Calcular la media de cada muestra y construye su respectivo intervalo de confianza usando el error estándar que calculamos en el paso 3.

# Paso 4 - Calcular la media de cada muestra y construye su respectivo # intervalo de confianza usando el error estándar que calculamos
# en el paso 3

set.seed(4) # Fijamos semilla de reproducibilidad

# Generamos las muestras:
datos %>%
  rep_sample_n(size = tamaño_muestra,
               replace = T,
               reps = numero_muestras) %>%
# Calculamos la estimación muestral y los límites del I.C.
  summarise(media_muestral = mean(valores),
            lim_sup = media_muestral + 1.96*error_est_con_datos_pob,
            lim_inf = media_muestral - 1.96*error_est_con_datos_pob) %>%  
# Obtenemos las observaciones que caen fuera del I.C., y
##  les asignamos un color característico
  mutate(color_intervalos = ifelse(lim_sup < 0 | lim_inf > 0,
                                   "salmon", "steelblue")) %>%
# Generamos la gráfica; declaramos los canales  
  ggplot(aes(x = replicate, y = media_muestral, color = color_intervalos)) +
# Generamos una línea recta que pasa por el 0 - EL parámetro
##  poblacional  
  geom_hline(yintercept = 0, color = "grey40", linetype = 2, size = 0.5) +
# Generamos una gráfica de puntos que marquen la media muestral  
  geom_point(size = 2) +
# Generamos un errorbar; una gráfica que genera segmentos
## de recta finitos, que marquen cada uno de nuestros I.C.  
  geom_errorbar(aes(ymin = lim_inf, ymax = lim_sup), width = 0.01,
                alpha = 0.82, size = 1) +

Paso 5. Añadimos elementos estéticos a la gráfica.

# Colores de las errorbars
  scale_color_manual(values = c("salmon", "steelblue")) +
# Etiquetas del eje x  
  labs(x = "\nNúmero de simulación") +
# Tema mínimo (fondo blanco y sin vida)
  theme_minimal() +
# Personalización del tema:  
  theme(axis.title.x = element_text(hjust = 1, face = "bold",
                                    color = "grey30", size = 15),
        axis.title.y = element_blank(),
        axis.text = element_text(size = 15),
        legend.position = "none")

2. Asume que el viejo Segasi no sabe nada de estadística. Por favor explícale de forma comprensible qué SÍ es y que NO es el nivel de confianza de un intervalo de confianza.

Que no es un Intervalo de Confianza:

NO ES la probabilidad de que un parámetro poblacional se encuentre dentro de cierto rango.

Que es un Intervalo de Confianza.

Son las veces que, repitiendo n veces (donde n tiende a infinito) un proceso de muestreo, los intervalos de confianza resultantes contienen al parámetro poblacional al menos un x % de veces, siendo x el nivel de confianza.

3. Por favor explica por qué si en el código anterior usas set.seed(4) al comienzo del paso 3, hay más de cinco intervalos de confianza rojos. ¿R se equivocó?

A ver 👀.

# FIJAMOS LA SEMILLA QUE DICE SEGASI
set.seed(4)

# Paso 3 - Calcular el error estándar de la media muestral, asumiendo
# que el tamaño de la muestra es 1,000
error_est_con_datos_pob <- sd(datos$valores)/sqrt(tamaño_muestra)

# Paso 4 - Calcular la media de cada muestra y construye su respectivo
# intervalo de confianza usando el error estándar que calculamos
# en el paso 3
set.seed(4)
datos %>%
  rep_sample_n(size = tamaño_muestra,
               replace = T,
               reps = numero_muestras) %>%
  summarise(media_muestral = mean(valores),
            lim_sup = media_muestral + 1.96*error_est_con_datos_pob,
            lim_inf = media_muestral - 1.96*error_est_con_datos_pob) %>%  
  mutate(color_intervalos = ifelse(lim_sup < 0 | lim_inf > 0,
                                   "salmon", "steelblue")) %>%
  ggplot(aes(x = replicate, y = media_muestral, color = color_intervalos)) +
  geom_hline(yintercept = 0, color = "grey40", linetype = 2, size = 0.5) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = lim_inf, ymax = lim_sup), width = 0.01,
                alpha = 0.82, size = 1) +
  scale_color_manual(values = c("salmon", "steelblue")) +
  labs(x = "\nNúmero de simulación") +
  theme_minimal() +
  theme(axis.title.x = element_text(hjust = 1, face = "bold",
                                    color = "grey30", size = 15),
        axis.title.y = element_blank(),
        axis.text = element_text(size = 15),
        legend.position = "none")

No, R no se equivocó. Lo que pasa es que, por probabilidad, es posible que de los primeros 100 Intervalos de Confianza simulados, haya más de 5 en los que no cae el parámetro poblacional. Sin embargo, si replicamos este modelo un número muy grande de veces, la poporción de intervalos que acaparen el parámetro poblacional efectivamente va a tender al 95 %.

4. Partiendo del código que tienes arriba, ajusta el código para (i) aumentar el número de muestras a 10,000; (ii) calcular la proporción de intervalos de confianza generados con estas muestras que incluyen al parámetro poblacional. NO TIENES QUE GRAFICAR LOS 100K INTERVALOS DE CONFIANZA.

# Paso 1 - Generar los datos poblacionales
set.seed(33)
datos <-
  tibble(valores = rnorm(n = 1e6, mean = 0, sd = 1))

# Paso 2 - Definir el tamaño de cada muestra y el número de muestras
tamaño_muestra <- 1000
numero_muestras <- 10000

# Paso 3 - Calcular el error estándar de la media muestral, asumiendo
# que el tamaño de la muestra es 1,000
error_est_con_datos_pob <- sd(datos$valores)/sqrt(tamaño_muestra)

# Paso 4 - Calcular la media de cada muestra y construye su respectivo
# intervalo de confianza usando el error estándar que calculamos
# en el paso 3
set.seed(4)
datos %>%
  rep_sample_n(size = tamaño_muestra,
               replace = T,
               reps = numero_muestras) %>%
  summarise(media_muestral = mean(valores),
            lim_sup = media_muestral + 1.96*error_est_con_datos_pob,
            lim_inf = media_muestral - 1.96*error_est_con_datos_pob) %>%  
  summarise(intervalo_incluye = sum(ifelse(lim_sup < 0 | lim_inf > 0, 0, 1))) %>%
  mutate(prop_incluye = intervalo_incluye/numero_muestras)

## # A tibble: 1 x 2
##   intervalo_incluye prop_incluye
##               <dbl>        <dbl>
## 1              9488       0.9488

Como podemos ver acá, ya ahora si la proporción de intervalos se aproxima al 95 %.

Ejercicio 3

En el paso 3 del ejercicio anterior usé la desviación estándar de la variable poblacional (misma que en el paso 1 definimos como sd = 1) para calcular el error estándar de la distribución muestral de la media.

En la vida real, es muy poco probable que conozcamos la desviación estándar de la variable poblacional, así que en este ejercicio asumiremos que no la conocemos.

Por favor calcula el error estándar de la distribución muestral de la media usando la estimación de la desviación estándar de la variable poblacional con los datos de 100 muestras. Guardarás el resultado de este cálculo en una variable llamada error_est_datos_muestra.

OJO: Dado que debes estimar la desviación estándar de la variable poblacional con los datos de cada muestra por separado, al final del ejercicio deberás tener un tibble con 100 renglones (uno por muestra) y tres columnas: replicate, media_muestral y error_est_datos_muestra.

# Paso 1 - Generar los datos poblacionales
set.seed(33)
datos <-
  tibble(valores = rnorm(n = 1e6, mean = 0, sd = 1))

# Paso 2 - Definir el tamaño de cada muestra y el número de muestras
tamaño_muestra <- 1000
numero_muestras <- 100

# Paso 3 - Calcular la media de cada muestra y el error estándar de la media muestral, estimando la desviación estándar de la variable poblacional.

set.seed(4)
datos %>%
  rep_sample_n(size = tamaño_muestra,
               replace = T,
               reps = numero_muestras) %>%
  summarise(media_muestral = mean(valores),
            error_est_datos_muestra = sd(valores)/sqrt(tamaño_muestra))

## # A tibble: 100 x 3
##    replicate media_muestral error_est_datos_muestra
##        <int>          <dbl>                   <dbl>
##  1         1    0.01835448               0.03142291
##  2         2    0.03685524               0.03065791
##  3         3    0.003059765              0.03218418
##  4         4    0.03721368               0.03173951
##  5         5    0.05931238               0.03180650
##  6         6    0.01661044               0.03241960
##  7         7   -0.02390758               0.03190266
##  8         8    0.02488777               0.03060030
##  9         9    0.01402477               0.03225584
## 10        10   -0.04237659               0.03127836
## # … with 90 more rows

1. ¿Los valores en la columna error_est_datos_muestra todos iguales? ¿Esto está bien o es señal de que hay un problema?

Los errores estándar de los datos de la muestra son todos distintos, porque es un error estándar por cada una de las estimaciones de la media de la muestra varía para cada muestra diferente.

2. ¿Los valores en la columna error_est_datos_muestra son exactamente iguales al valor que obtuvimos en el paso 3 del ejercicio anterior? ¿Esto está bien o es señal de que hay un problema?

No son iguales (aunque si muy similares), ya que el error estándar del ejercicio 3 se calculó con el valor de la desviación estándar de todos los datos de la población, mientras que el error estándar calculado en este ejercicio se realizó con la estimación muestral de esta desviación estándar.

Ejercicio 4

Repasemos la diferencia entre error muestral y error estándar. Para ello usaremos los datos de una distribución gamma. En esta (p. 910) y esta liga encontrarás más detalles sobre esta distribución.

1. Usa el siguiente código para construir la distribución:

datos <- tibble(valores = rgamma(n = 1e5, shape = 2, scale = 1))

2. Toma 100 muestras aleatorias de tamaño 120, calcula la media de cada muestra y el o los error(es) muestral(es)

datos %>%
  rep_sample_n(size = 120, reps = 100) %>%
# Calculamos la media de cada muestra
  summarise(media_muestral = mean(valores)) %>%
# Calculamos el error muestral de cada muestra con el Par. Pob.
  mutate(error_muestral = media_muestral - mean(datos$valores))

## # A tibble: 100 x 3
##    replicate media_muestral error_muestral
##        <int>          <dbl>          <dbl>
##  1         1       1.995074    0.004581609
##  2         2       1.912890   -0.07760184 
##  3         3       2.076706    0.08621392 
##  4         4       1.865318   -0.1251737  
##  5         5       1.732433   -0.2580587  
##  6         6       1.922717   -0.06777534 
##  7         7       2.043008    0.05251626 
##  8         8       2.016195    0.02570304 
##  9         9       2.123461    0.1329694  
## 10        10       2.241213    0.2507209  
## # … with 90 more rows

3. Toma 100 muestras aleatorias de tamaño 120, calcula la media de cada muestra y el o los error(es) estándar.

datos %>%
  rep_sample_n(size = 120, reps = 100) %>%
  summarise(media_muestral = mean(valores),
            error_estandar_media_muestral = sd(valores)/sqrt(120)) %>%
  summarise(error_estandar = sd(media_muestral))

## # A tibble: 1 x 1
##   error_estandar
##            <dbl>
## 1      0.1122577

Ejercicio 5

Basta de simulaciones. Bueno, casi. En este ejercicio tendrás que seguir el caminito lleno de incertidumbre de toda persona que trata de entender las características de una población a partir de una muestra.

1. Asume que 67% de la población de México está a favor del paro nacional del 9 de marzo “Un día sin mujeres” (liga) y el 33% restante está en contra. Construye una población hipotética de 1 millón de personas con estas proporciones. Puedes representar a los que están a favor del paro con un 1 y a los que están en contra con un 0.

2. Olvida que sabes las proporciones poblacionales.

3. Planta un set.seed(1) y toma una sola muestra de tamaño 500.

4. Usando los datos de tu muestra, estima la proporción de personas que está a favor del paro.

5. Usando los datos de tu muestra, construye un intervalo de confianza de 99%.

6. Reporta los resultados mostrando la/el pro que eres.

Muestra de tamaño 500

poblacion <- tibble(preferencia = c(rep(1, 670000), rep(0, 330000)))

tamaño_muestra <-  500

# Mi propuesta :3
poblacion %>%
  rep_sample_n(size = tamaño_muestra, reps = 1) %>%
  summarise(proporcion_muestral = mean(preferencia),
            EE = sd(preferencia)/sqrt(tamaño_muestra)) %>%
  mutate(lim_inf = proporcion_muestral - 2.58*EE,
         lim_sup = proporcion_muestral + 2.58*EE,
         mult = 2.58*EE)

## # A tibble: 1 x 6
##   replicate proporcion_muestral         EE   lim_inf   lim_sup       mult
##       <int>               <dbl>      <dbl>     <dbl>     <dbl>      <dbl>
## 1         1             0.68800 0.02074060 0.6344893 0.7415107 0.05351074

Reportar como pro:

Según una encuesta realizada por los alumnos de CPRI del CIDE a inicios de marzo del 2020, un 68.3 % (mas/menos 3.7%) de la población mexicana se encuentra a favor del paro del 9 de marzo. El intervalo de confianza abarca desde 66.46 a un 73.93 %, con un nivel de confianza del 99 %.

7. Repite los pasos c), d) y e) usando una muestra de tamaño 100, y posteriormente una muestra de tamaño 1,000. ¿Qué efecto tiene el tamaño de la muestra en el tamaño del intervalo de confianza?

Muestra de tamaño 100

poblacion <- tibble(preferencia = c(rep(1, 670000), rep(0, 330000)))

tamaño_muestra <-  100

# Mi propuesta :3
poblacion %>%
  rep_sample_n(size = tamaño_muestra, reps = 1) %>%
  summarise(proporcion_muestral = mean(preferencia),
            EE = sd(preferencia)/sqrt(tamaño_muestra)) %>%
  mutate(lim_inf = proporcion_muestral - 2.58*EE,
         lim_sup = proporcion_muestral + 2.58*EE)

## # A tibble: 1 x 5
##   replicate proporcion_muestral         EE   lim_inf   lim_sup
##       <int>               <dbl>      <dbl>     <dbl>     <dbl>
## 1         1                0.52 0.05021167 0.3904539 0.6495461

Muestra de tamaño 1,000

poblacion <- tibble(preferencia = c(rep(1, 670000), rep(0, 330000)))

tamaño_muestra <-  1000

# Mi propuesta :3
poblacion %>%
  rep_sample_n(size = tamaño_muestra, reps = 1) %>%
  summarise(proporcion_muestral = mean(preferencia),
            EE = sd(preferencia)/sqrt(tamaño_muestra)) %>%
  mutate(lim_inf = proporcion_muestral - 2.58*EE,
         lim_sup = proporcion_muestral + 2.58*EE)

## # A tibble: 1 x 5
##   replicate proporcion_muestral         EE   lim_inf   lim_sup
##       <int>               <dbl>      <dbl>     <dbl>     <dbl>
## 1         1             0.69700 0.01453968 0.6594876 0.7345124

R: A medida en que el tamaño de la muestra va creciendo, el intervalo de confianza se va volviendo más angosto.