Carga de librerías

library(tidyverse)
library(glue)
library(furrr)
plan(multiprocess)
set.seed(42)

Escenarios posibles

Generamos las combinaciones de distribución del error y tamaño muestral \(n\) que se piden. Son 45 escenarios posibles y cada uno será repetido 1000 veces, lo que da una tabla de 45.000 filas:

datos <- cross_df(list(
  n_obs = c(10, 25, 100, 250, 500, 1000, 1500, 2000, 3000),
  tipo_de_error = c("exponencial", "lognormal", "uniforme", "chi_cuadrado", "t_student"),
  run_id = 1:1000
))

nrow(datos)

## [1] 45000

Para cada uno de esos 45 x 1.000 escenarios, generamos las observaciones correspondientes:

generar_observaciones <- function(n_obs, tipo_de_error) {
  observaciones <- tibble(
    i = 1:n_obs,
    x1 = runif(n_obs, -5, 5),
    x2 = runif(n_obs, -5, 5),
    x3 = runif(n_obs, -5, 5),
    x4 = runif(n_obs, -5, 5)
  )
  
  if (tipo_de_error == "exponencial") {
    errores <- rexp(n_obs, 0.5) - 2
  } else if (tipo_de_error == "lognormal") {
    errores <- rlnorm(n_obs, 0, 1) - exp(-0.5)
  } else if (tipo_de_error == "uniforme") {
    errores <- runif(n_obs, -3, 3)
  } else if (tipo_de_error == "chi_cuadrado") {
    errores <- rchisq(n_obs, df = 3) - 3
  } else if (tipo_de_error == "t_student") {
    errores <- rt(n_obs, df = 3)
  } else {
    stop(glue("No sé generar el tipo de error {tipo_de_error}"))
  }
  
  observaciones <- mutate(observaciones,
                          err = errores,
                          y = 4 + 2 * x1 - 3 * x2 + 0.5 * x3 + err)
  
  return (observaciones)
}

datos <- mutate(datos, observaciones = future_map2(n_obs, tipo_de_error,
                                                   generar_observaciones))

Ajuste de cada escenario

Ajustamos el modelo para cada grupo de observaciones:

ajustar_modelo <- function(observaciones) {
  return (lm(y ~ x1 + x2 + x3 + x4, data = observaciones))
}
datos <- mutate(datos, modelo = future_map(observaciones, ajustar_modelo))

Intervalos de confianza para cada grupo de observaciones

Construimos en cada caso los IC de nivel 0.90 y de nivel 0.90 asintótico para β1 y β4, a partir del elemento “modelo” que produce lm:

calcular_IC <- function(modelo, j, tipo = "exacto", alpha = 0.1) {
  j <- j + 1 # La matriz empieza en beta 0, el beta_1 está en el índice 2, etc.
  X <- model.matrix(modelo)
  n_obs <- nrow(X)
  p <- modelo$rank
  beta_j_estimado <- modelo$coefficients[[j]]
  desvio_del_beta_j <- sqrt(vcov(modelo)[[j, j]])
  if (tipo == "exacto") {
    cuantil <- qt(1 - alpha/2, n_obs - p)
  } else if (tipo == "asintotico") {
    cuantil <- qnorm(1 - alpha/2)
  } else {
    stop(glue("No conozco el tipo de IC '{tipo}'"))
  }
  delta <- cuantil * desvio_del_beta_j
  return(c(beta_j_estimado - delta, beta_j_estimado + delta))
}

datos <-
  datos %>%
  mutate(
    IC_1 = map(modelo, ~ calcular_IC(.x, j = 1)),
    lower_1 = map_dbl(IC_1, ~ .x[[1]]),
    upper_1 = map_dbl(IC_1, ~ .x[[2]]),
    
    IC_4 = map(modelo, ~ calcular_IC(.x, j = 4)),
    lower_4 = map_dbl(IC_4, ~ .x[[1]]),
    upper_4 = map_dbl(IC_4, ~ .x[[2]]),
    
    IC_1_asnt = map(modelo, ~ calcular_IC(.x, j = 1, tipo = "asintotico")),
    lower_1_asnt = map_dbl(IC_1_asnt, ~ .x[[1]]),
    upper_1_asnt = map_dbl(IC_1_asnt, ~ .x[[2]]),
    
    IC_4_asnt = map(modelo, ~ calcular_IC(.x, j = 4, tipo = "asintotico")),
    lower_4_asnt = map_dbl(IC_4_asnt, ~ .x[[1]]),
    upper_4_asnt = map_dbl(IC_4_asnt, ~ .x[[2]])
  ) %>%
  select(-IC_1, -IC_4, -IC_1_asnt, -IC_4_asnt)

Chequeamos para cada escenario si el β1 y el β4 estimados están contenidos en sus IC y en sus IC asintóticos:

BETA_1_VERDADERO <- 2
BETA_4_VERDADERO <- 0

datos <- mutate(datos,
                beta_1_en_IC = (BETA_1_VERDADERO >= lower_1 & BETA_1_VERDADERO <= upper_1),
                beta_4_en_IC = (BETA_4_VERDADERO >= lower_4 & BETA_4_VERDADERO <= upper_4),
                beta_1_en_IC_asnt = (BETA_1_VERDADERO >= lower_1_asnt & BETA_1_VERDADERO <= upper_1_asnt),
                beta_4_en_IC_asnt = (BETA_4_VERDADERO >= lower_4_asnt & BETA_4_VERDADERO <= upper_4_asnt))

Extraemos los betas estimados que nos interesan a columnas separadas:

datos <- mutate(datos,
                betas_estimados = map(modelo, ~ .x$coefficients),
                beta_1_est = map_dbl(betas_estimados, ~ .x[[2]]),
                beta_2_est = map_dbl(betas_estimados, ~ .x[[3]]),
                beta_4_est = map_dbl(betas_estimados, ~ .x[[5]]))

Un caso solo de cada escenario

Para un solo caso por separado, en cada escenario, vemos si:

• Los IC nivel 0.90 contienen a los parámetros β1 y β4
• Los IC nivel asintótico 0.90 contienen a los parámetros β1 y β4

Nos quedamos con una sola iteración de cada escenario:

resumen_1_caso <-
  datos %>%
  filter(run_id == 1) %>%
  select(matches("n_obs|tipo_de_error|lower|upper|beta_(1|4)")) %>%
  gather(key = "name", value = "beta_value", -n_obs, -tipo_de_error) %>%
  mutate(tipo_de_intervalo = ifelse(str_detect(name, "asnt"), "asintótico", "no asintótico"),
         tipo_de_error = map_chr(tipo_de_error, ~ glue("error {.x}")))

Graficamos los IC en relación a los parámetros verdaderos:

resumen_1_caso %>%
  filter(!str_detect(name, "beta_._en_IC")) %>%
  ggplot() +
  facet_grid(tipo_de_error ~ .) +
  aes(y = as.factor(n_obs), x = beta_value) +
  geom_point(
     data = filter(resumen_1_caso, str_detect(name, "lower|upper")),
     mapping = aes(shape = str_detect(name, "lower"),
                   color = tipo_de_intervalo,
                   size = tipo_de_intervalo),
     size = 4,
  ) +
  geom_point(
     data = filter(resumen_1_caso, str_detect(name, "beta_._est")),
     shape = 4
  ) +
  geom_vline(xintercept = BETA_1_VERDADERO, size = .3, alpha = .5) +
  geom_vline(xintercept = BETA_4_VERDADERO, size = .3, alpha = .5) +
  scale_shape_manual(values = c(93, 91)) +
  guides(shape = F) +
  labs(
    title = "Intervalos de Confianza para β1 y β4",
    subtitle = glue("β1 verdadero = {BETA_1_VERDADERO}; β4 verdadero = {BETA_4_VERDADERO}"),
    x = "β",
    y = "n",
    color = "Tipo de IC:"
  ) +
  theme(legend.position = "top") +
  xlim(-1, 3)

## Warning: Removed 7 rows containing missing values (geom_point).

library(gt)

resumen_1_caso %>%
  filter(str_detect(name, "beta_._en_IC")) %>%
  rename(
    Contenido = beta_value,
    n = n_obs,
    `Tipo de IC` = tipo_de_intervalo,
    `Parámetro` = name
  ) %>%
  mutate(
    Contenido = ifelse(Contenido, "✓", "☹"),
    `Parámetro` = ifelse(str_detect(`Parámetro`, "beta_1"), "β1", "β4")
  ) %>%
  spread(n, Contenido) %>%
  group_by(tipo_de_error) %>%
  gt %>%
  tab_spanner("n", columns = matches("[1-9]+")) %>%
  tab_source_note(
    source_note = "✓: Parámetro contenido en el intervalo | ☹: Parámetro no contenido en el intervalo"
  )

1.000 iteraciones de cada escenario

Para cada tipo de escenario, repetimos 1000 veces y analizamos:

• Qué proporción de los IC nivel 0.90 contienen a los parámetros \(β_{1}\) y \(β_{4}\)
• Qué proporción de los IC nivel asintótico 0.90 contienen a los parámetros \(β_{1}\) y \(β_{4}\)

analisis_IC <-
  datos %>%
  # filter(n_obs %in% c(10, 25, 100)) %>%
  group_by(n_obs, tipo_de_error) %>%
  summarise(
    `β1` = mean(beta_1_en_IC),
    `β4` = mean(beta_4_en_IC),
    `β1 asnt` = mean(beta_1_en_IC_asnt),
    `β4 asnt` = mean(beta_4_en_IC_asnt)
  ) %>%
  ungroup %>%
  rename(
    Modelo = tipo_de_error,
    n = n_obs
  ) %>%
  mutate(Modelo = map_chr(Modelo, ~ glue("Error {.x}")))

analisis_IC %>%
  group_by(n) %>%
  group_by(Modelo) %>%
  gt(
    rowname_col = "tipo_de_error",
    groupname_col = ""
  ) %>%
  tab_header(title = "Proporción de IC que contienen a los parámetros según n y tipo de error") %>%
  tab_spanner(label = "IC 90%", columns = vars(`β1`, `β4`)) %>%
  tab_spanner(label = "IC 90% asintótico", columns = vars(`β1 asnt`, `β4 asnt`)) %>%
  data_color(
    columns = matches("β"),
    colors = scales::col_bin(palette = c("Tomato", "Black"), bins = c(0, 0.89, 1)),
    apply_to = "text"
  ) %>%
  cols_align("center") %>%
  tab_footnote(
    footnote = "En rojo los escenarios donde menos del 89% de los intervalos contuvo al parámetro.",
    locations = cells_column_labels(columns = matches("β"))
  ) %>%
  tab_options(footnote.glyph = c("*"))

datos %>%
  select(n_obs, tipo_de_error, beta_2_est) %>%
  mutate(
    tipo_de_error = map_chr(tipo_de_error, ~glue("Error\n{.x}")),
    n_obs_factor = map_chr(n_obs, ~ glue("n = {.x}"))
  ) %>%
  ggplot() +
  facet_grid(fct_reorder(n_obs_factor, n_obs) ~ tipo_de_error, scales = "free") +
  aes(sample = beta_2_est) +
  stat_qq(size = .5) +
  stat_qq_line(color = "Red", alpha = .7) +
  theme_minimal() +
  theme(
    axis.text = element_blank(), 
    axis.ticks = element_blank()
  ) +
  labs(
    title = "Distribución del β2 estimado comparada con una distribución normal",
    subtitle = "1.000 observaciones en cada grupo",
    x = "", y = ""
  )

library(bayestestR)

BETA_2_VERDADERO <- -3

beta_2_comparacion <-
  datos %>%
  select(n_obs, tipo_de_error, beta_2_est) %>%
  group_by(n_obs, tipo_de_error) %>%
  mutate(
    desvio_beta_2 = sd(beta_2_est),
    # beta_2_normal = rnorm(n(), mean = BETA_2_VERDADERO, sd = desvio_beta_2)
    beta_2_normal = distribution_normal(n(), mean = BETA_2_VERDADERO, sd = desvio_beta_2)
  ) %>%
  ungroup() %>%
  select(-desvio_beta_2) %>%
  gather("origen", "beta_value", -n_obs, -tipo_de_error) %>%
  mutate(
    escenario = map2_chr(n_obs, tipo_de_error, ~ glue("Error {.y}\nn = {.x}")),
    escenario = fct_reorder(escenario, n_obs),
    origen = ifelse(origen == "beta_2_est", "β2 estimado", "β2 de distribución normal")
  )

beta_2_comparacion %>%
  ggplot() +
  facet_wrap(escenario ~ ., scales = "free", ncol = 5) +
  aes(x = beta_value, color = origen, linetype = origen, size = origen) +
  # geom_density() +
  geom_line(stat = "density") +
  geom_vline(xintercept = BETA_2_VERDADERO, size = .1) +
  scale_color_manual(values = c("Tomato", "Black")) +
  scale_size_manual(values = c(.8, .55)) +
  scale_linetype_manual(values = c("solid", "solid")) +
  theme_minimal() +
  labs(
    title = "Distribución del β2 estimado comparada con una distribución normal",
    subtitle = "1.000 observaciones en cada grupo",
    y = "", x = "", color = "", linetype = "", size = ""
  ) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    legend.position = "top"
  )