🔷 Visión General del Script

Este primer script constituye la puerta de entrada del pipeline completo. Su propósito es cargar, limpiar y preparar los datos provenientes de la encuesta aplicada a consumidores de Amazon en Estados Unidos, transformando las respuestas crudas en un conjunto de variables estructuradas listas para el modelo POGIT.

El script no estima ningún parámetro del modelo. Su rol es exclusivamente de preparación y exploración: construir la base de covariables demográficas y de comportamiento que posteriormente alimentarán los vectores \(x_{1i}\) y \(x_{2i}\) del modelo.

🗺️ Conexión con el Paper

Bloque de Código	Sección del Paper
Carga del survey	Sección 2 — Data
Limpieza de nombres	Sección 3 — Construction of Predictor Variables
Race / Life dummies	Sección 3 — Construction of Predictor Variables
Merge final	Sección 3 — Construction of Predictor Variables
Gráficos descriptivos	EDA implícito (no entra al modelo)

📐 Lógica del Pipeline

El flujo de transformación de este script puede resumirse así:

CSV crudo (survey.csv)
        ↓
  clean_names()  →  nombres estandarizados
        ↓
  race dummies   →  multi-respuesta → columnas binarias
        ↓
  life dummies   →  multi-respuesta → columnas binarias
        ↓
  left_join()    →  survey_enriched (dataset final)
        ↓
  bar_plot()     →  visualización exploratoria
        ↓
  PDF export     →  documentación visual

Este dataset final (survey_enriched) es el output clave del script y el punto de partida del script 2.

📦 Bloque 1 — Carga de Paquetes

Explicación Teórica

Antes de cualquier procesamiento, R necesita tener disponibles todas las herramientas del pipeline. Este bloque carga 30 paquetes organizados en grupos funcionales.

La función load_pkg() es una mejora sobre el simple library(): en lugar de fallar silenciosamente o lanzar mensajes confusos, verifica que cada paquete se haya cargado correctamente y lanza un error explícito indicando exactamente cuál(es) fallaron. Esto es especialmente útil en entornos colaborativos o cuando el script se ejecuta en una máquina diferente.

Paquetes por función en el modelo POGIT

Grupo	Paquetes	Rol en el pipeline
Manipulación	`tidyverse`, `dplyr`, `tidyr`, `stringr`	Transformación de datos en todos los scripts
Fechas	`lubridate`	Ventanas temporales T1/T2 (script 2)
Limpieza	`janitor`, `skimr`	Estandarización de nombres y resúmenes
Encoding	`fastDummies`	One-hot encoding de variables categóricas
Visualización	`ggplot2`, `patchwork`, `scales`, `viridis`	Gráficos descriptivos y de validación
Modelado	`tidymodels`, `caret`, `pscl`	Infraestructura del modelo POGIT
Optimización	`optimParallel`, `parallel`, `furrr`	Estimación MLE en paralelo (script 3)
Derivadas	`numDeriv`	Cálculo de errores estándar del modelo
Exportación	`openxlsx`, `readxl`	Lectura/escritura de matrices Excel
Métricas	`Metrics`	MAE y otras métricas de evaluación

Nota teórica: La necesidad de optimParallel y numDeriv refleja que la estimación del modelo POGIT se hace por Máxima Verosimilitud (MLE), lo cual requiere optimización numérica iterativa y cálculo de derivadas para obtener los errores estándar de \(\hat{\beta}_1\) y \(\hat{\beta}_2\).

packages <- unique(c(
  "tidyverse", "lubridate", "janitor", "skimr", "snakecase",
  "stringr", "dplyr", "gridExtra", "scales", "viridis",
  "patchwork", "fastDummies", "readxl", "openxlsx", "tidymodels",
  "Matrix", "caret", "pscl", "xgboost", "Metrics",
  "numDeriv", "optimParallel", "parallel", "furrr", "progressr",
  "tibble", "readr", "forcats", "ggplot2", "tidyr"
))

load_pkg <- function(pkgs) {
  ok <- vapply(pkgs, function(p) {
    suppressPackageStartupMessages(require(p, character.only = TRUE))
  }, logical(1))
  if (any(!ok)) {
    stop("Packages not installed/loaded: ", paste(pkgs[!ok], collapse = ", "))
  }
  invisible(ok)
}
load_pkg(packages)

✅ Output esperado

Si todos los paquetes están instalados, este bloque no produce output visible. Si alguno falla, el error indica exactamente cuál instalar con install.packages("nombre").

🎨 Bloque 2 — Configuración Visual y Paleta de Colores

Explicación Teórica

Este bloque cumple dos funciones de soporte:

1. Configuración del locale: Se establece es_ES.UTF-8 para que los separadores decimales, nombres de meses y caracteres especiales (tildes, ñ) sean consistentes con el contexto hispanohablante del proyecto. Si el sistema operativo no soporta ese locale, el script continúa con una advertencia (no con un error), garantizando portabilidad entre máquinas.

2. Paleta de colores institucional (pal): Se define un vector nombrado de 10 colores que se usa en todos los gráficos del proyecto. Definirla como objeto global antes de las funciones es crítico porque bar_plot() la recibe como argumento por defecto (fill_pal = pal). Si se definiera después, R lanzaría un error al intentar usar la función.

Nota: Esta paleta no tiene ningún rol en la estimación del modelo. Es puramente de presentación, pero su consistencia facilita la interpretación visual de los resultados en la Sección 2 del paper (descripción de la muestra).

loc_ok <- base::Sys.setlocale("LC_ALL", "es_ES.UTF-8")
if (identical(loc_ok, "")) {
  warning(
    "Could not set locale 'es_ES.UTF-8'; ",
    "separators and month names will use the default locale."
  )
}

pal <- c(
  azul_1   = "#4E79A7",
  naranja  = "#F28E2B",
  rojo     = "#E15759",
  turquesa = "#76B7B2",
  verde    = "#59A14F",
  amarillo = "#EDC948",
  morado   = "#B07AA1",
  rosa     = "#FF9DA7",
  cafe     = "#9C755F",
  gris     = "#BAB0AC"
)

# Verificar que la paleta se creó correctamente
print(pal)

##    azul_1   naranja      rojo  turquesa     verde  amarillo    morado      rosa 
## "#4E79A7" "#F28E2B" "#E15759" "#76B7B2" "#59A14F" "#EDC948" "#B07AA1" "#FF9DA7" 
##      cafe      gris 
## "#9C755F" "#BAB0AC"

✅ Output esperado

Se imprime el vector pal con los 10 colores nombrados y sus códigos hexadecimales.

📊 Bloque 3 — Función `bar_plot()` y Funciones Auxiliares

Explicación Teórica

La función `bar_plot()`

Esta es la función central de visualización del script. Genera gráficos de barras horizontales estandarizados para cualquier variable categórica del survey, con tres características importantes:

Robustez: Filtra NA automáticamente y recicla la paleta de colores si hay más categorías que colores disponibles, evitando errores en variables con muchos niveles (ej: estado de residencia con 52 opciones).
Informatividad: Muestra simultáneamente frecuencias absolutas (longitud de la barra) y porcentajes (etiqueta de texto), que es la representación estándar en análisis descriptivo de encuestas.
Consistencia visual: Usa fct_reorder() para ordenar las barras de mayor a menor, facilitando la comparación entre categorías.

La función `mode_safe()`

Calcula la moda de un vector de strings de forma segura: ignora NA y strings vacíos, y retorna NA_character_ si el vector está completamente vacío. Se usa en el script 2 para identificar el estado de envío más frecuente de cada cliente.

La función `reorder_drop_id()`

Filtra y reordena un dataframe para alinear sus filas con un vector de referencia de IDs. Es esencial para garantizar que las matrices V, W e Y estén en el mismo orden de clientes antes de pasarlas al optimizador del modelo.

Conexión con el paper: Estas funciones implementan el principio de la Sección 2 de que las covariables del encuestado se usan para construir \(x_{1i}\) y \(x_{2i}\). La correcta alineación por ID garantiza que cada fila de la matriz corresponda exactamente al cliente correcto.

# Función principal de visualización
bar_plot <- function(df, var, title, x_lab,
                     fill_pal = pal, text_size = 3) {
  
  counts <- df %>%
    dplyr::filter(!is.na({{ var }})) %>%
    dplyr::count({{ var }}, name = "n") %>%
    dplyr::mutate(
      pct = n / sum(n) * 100,
      var_chr = as.character({{ var }})
    ) %>%
    dplyr::arrange(dplyr::desc(n)) %>%
    dplyr::mutate(var_chr = forcats::fct_reorder(var_chr, n))
  
  fill_vals <- rep(unname(fill_pal), length.out = nrow(counts))
  names(fill_vals) <- levels(counts$var_chr)
  
  ggplot2::ggplot(counts, ggplot2::aes(x = var_chr, y = n, fill = var_chr)) +
    ggplot2::geom_col() +
    ggplot2::geom_text(
      ggplot2::aes(label = sprintf("%.1f%%", pct), y = n),
      hjust = -0.1, size = text_size, colour = "black"
    ) +
    ggplot2::coord_flip() +
    ggplot2::scale_y_continuous(expand = ggplot2::expansion(mult = c(0, 0.15))) +
    ggplot2::scale_fill_manual(values = fill_vals, guide = "none") +
    ggplot2::labs(title = title, x = x_lab, y = "Number of participants") +
    ggplot2::theme_minimal()
}

# Funciones auxiliares
mode_safe <- function(x) {
  x <- x[!is.na(x) & x != ""]
  if (length(x) == 0) return(NA_character_)
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

reorder_drop_id <- function(df, ref_ids) {
  df %>%
    dplyr::filter(response_id %in% ref_ids) %>%
    dplyr::slice(match(ref_ids, response_id))
}

# Verificar que las funciones quedaron definidas
cat("✓ Funciones definidas correctamente:\n")

## ✓ Funciones definidas correctamente:

cat("  - bar_plot()        → visualización de variables categóricas\n")

##   - bar_plot()        → visualización de variables categóricas

cat("  - mode_safe()       → moda robusta para strings\n")

##   - mode_safe()       → moda robusta para strings

cat("  - reorder_drop_id() → alineación de IDs entre matrices\n")

##   - reorder_drop_id() → alineación de IDs entre matrices

✅ Output esperado

Mensaje de confirmación de que las tres funciones fueron definidas sin errores.

📂 Bloque 4 — Carga del Survey

Explicación Teórica

Este bloque carga el archivo principal de la encuesta, que es la primera fuente de datos del estudio (Sección 2 del paper).

Origen de los datos

El dataset proviene de un protocolo de crowdsourcing con consentimiento informado, aprobado por el IRB del MIT (protocolo #2205000649). Contiene respuestas de 5,027 participantes en Estados Unidos que voluntariamente compartieron sus datos demográficos, de estilo de vida y sus historiales de compras en Amazon.

¿Por qué `clean_names()`?

Los archivos de encuesta exportados de plataformas como Qualtrics o SurveyMonkey suelen tener nombres de columnas inconsistentes: - Mezcla de mayúsculas y minúsculas: Q_Demos_Age - Espacios: Q Demos Age - Caracteres especiales: Q_Demos_Age?

clean_names() del paquete janitor los convierte todos a snake_case minúsculo (q_demos_age), eliminando espacios y caracteres especiales. Esto es fundamental para que el resto del código funcione de forma predecible.

Variables clave del survey

Variable	Descripción	Rol en el modelo
`q_demos_age`	Grupo de edad	Dummy en \(x_{1i}\) (SoW)
`q_demos_race`	Raza (multi-select)	Dummy en \(x_{1i}\) (SoW)
`q_demos_income`	Rango de ingreso	Dummy en \(x_{2i}\) (SioW)
`q_amazon_use_howmany`	Personas que comparten cuenta	Dummy en \(x_{2i}\) (SioW)
`q_life_changes`	Eventos de vida 2021 (multi-select)	Dummy en \(x_{2i}\) (SioW)
`q_demos_hispanic`	Origen hispano/latino	Dummy en \(x_{1i}\) (SoW)

survey <- readr::read_csv(
  "data/survey-data/survey.csv",
  show_col_types = FALSE
) %>%
  janitor::clean_names()

# Verificación del dataset cargado
cat("============================================================\n")

## ============================================================

cat("SURVEY CARGADO EXITOSAMENTE\n")

## SURVEY CARGADO EXITOSAMENTE

cat("============================================================\n")

## ============================================================

cat("Dimensiones:", nrow(survey), "filas x", ncol(survey), "columnas\n")

## Dimensiones: 6325 filas x 41 columnas

cat("Participantes únicos:", n_distinct(survey$response_id), "\n\n")

## Participantes únicos: 6325

cat("Primeras 5 columnas del dataset:\n")

## Primeras 5 columnas del dataset:

print(head(survey[, 1:5], 3))

## # A tibble: 3 × 5
##   duration_in_seconds recorded_date     response_id q_prolific_mturk q_demos_age
##                 <dbl> <chr>             <chr>       <chr>            <chr>      
## 1                 332 9/21/2022 10:00:… R_1ou69fj4… <NA>             35 - 44 ye…
## 2                 488 9/21/2022 10:02:… R_24dboHVO… <NA>             25 - 34 ye…
## 3                 309 9/21/2022 10:10:… R_2UbJL30H… <NA>             45 - 54 ye…

cat("\nNombres de columnas (primeras 20):\n")

## 
## Nombres de columnas (primeras 20):

print(names(survey)[1:20])

##  [1] "duration_in_seconds"  "recorded_date"        "response_id"         
##  [4] "q_prolific_mturk"     "q_demos_age"          "q_demos_hispanic"    
##  [7] "q_demos_race"         "q_demos_education"    "q_demos_income"      
## [10] "q_demos_gender"       "q_sexual_orientation" "q_demos_state"       
## [13] "q_amazon_use_howmany" "q_amazon_use_hh_size" "q_amazon_use_how_oft"
## [16] "q_substance_use_1"    "q_substance_use_2"    "q_substance_use_3"   
## [19] "q_personal_1"         "q_personal_2"

✅ Output esperado

Confirmación de carga con dimensiones (≈5,027 filas)
Muestra de las primeras filas y columnas con nombres en snake_case

🔢 Bloque 5 — Procesamiento de Raza (Multi-respuesta → Dummies)

Explicación Teórica

El problema de las variables multi-respuesta

La pregunta de raza en la encuesta permite selección múltiple: un participante puede identificarse con más de una categoría racial simultáneamente. Esto genera un problema de representación: la celda contiene un string con múltiples valores separados por coma, como "White or Caucasian, Asian".

Este formato no puede entrar directamente al modelo porque los algoritmos de regresión necesitan columnas numéricas individuales.

La solución: Expand → Encode

El pipeline de transformación tiene dos etapas:

Etapa 1 — Formato largo (race_long):

\[\text{"White, Asian"} \xrightarrow{\text{separate\_rows()}} \begin{cases} \text{id=1, race="White"} \\ \text{id=1, race="Asian"} \end{cases}\]

Etapa 2 — One-hot encoding (race_dummies):

\[\xrightarrow{\text{pivot\_wider()}} \begin{array}{|c|c|c|} \hline \text{id} & \text{race\_white} & \text{race\_asian} \\ \hline 1 & 1 & 1 \\ \hline \end{array}\]

¿Por qué NO se elimina una categoría de referencia aquí?

El Apéndice 7 del paper establece una regla clara:

“Las variables de selección múltiple se codifican como un conjunto completo de indicadores binarios, ya que las categorías no son mutuamente excluyentes.”

A diferencia de variables de selección única (donde se elimina una categoría para evitar multicolinealidad), en variables multi-select las categorías son independientes entre sí: ser “White” no implica no ser “Asian”. Por eso se conservan todas las columnas binarias.

Las 6 categorías raciales

race_levels <- c(
  "White or Caucasian",
  "Black or African American",
  "Asian",
  "American Indian/Native American or Alaska Native",
  "Native Hawaiian or Other Pacific Islander",
  "Other"
)

# Paso 1: Formato largo (una fila por raza declarada)
race_long <- survey %>%
  dplyr::select(response_id, race_raw = q_demos_race) %>%
  dplyr::filter(!is.na(race_raw) & race_raw != "") %>%
  dplyr::mutate(
    race_raw = race_raw %>%
      stringr::str_replace_all("\\s*/\\s*", "/") %>%
      stringr::str_replace_all("\\s*,\\s*", ", ") %>%
      stringr::str_squish()
  ) %>%
  tidyr::separate_rows(race_raw, sep = ",\\s*") %>%
  dplyr::mutate(race_raw = stringr::str_trim(race_raw))

# Paso 2: One-hot encoding (formato wide con columnas binarias)
race_dummies <- race_long %>%
  dplyr::mutate(
    race  = factor(race_raw, levels = race_levels),
    value = 1L
  ) %>%
  dplyr::select(-race_raw) %>%
  tidyr::pivot_wider(
    names_from   = race,
    values_from  = value,
    names_prefix = "race_",
    values_fill  = 0L
  ) %>%
  janitor::clean_names()

# Garantizar que todas las columnas de raza existen (aunque nadie la haya seleccionado)
needed_race <- paste0("race_", janitor::make_clean_names(race_levels))
race_dummies[setdiff(needed_race, names(race_dummies))] <- 0L

# Verificación
cat("============================================================\n")

## ============================================================

cat("RACE DUMMIES — VERIFICACIÓN\n")

## RACE DUMMIES — VERIFICACIÓN

cat("============================================================\n")

## ============================================================

cat("Filas en race_long (una por raza declarada):", nrow(race_long), "\n")

## Filas en race_long (una por raza declarada): 6690

cat("Filas en race_dummies (una por participante):", nrow(race_dummies), "\n")

## Filas en race_dummies (una por participante): 6325

cat("Columnas generadas:\n")

## Columnas generadas:

print(names(race_dummies))

## [1] "response_id"                                          
## [2] "race_black_or_african_american"                       
## [3] "race_white_or_caucasian"                              
## [4] "race_asian"                                           
## [5] "race_other"                                           
## [6] "race_american_indian_native_american_or_alaska_native"
## [7] "race_native_hawaiian_or_other_pacific_islander"

cat("\nDistribución por raza (conteo individual):\n")

## 
## Distribución por raza (conteo individual):

print(
  race_long %>%
    dplyr::count(race_raw, sort = TRUE) %>%
    dplyr::mutate(pct = scales::percent(n / nrow(survey), accuracy = 0.1))
)

## # A tibble: 6 × 3
##   race_raw                                             n pct  
##   <chr>                                            <int> <chr>
## 1 White or Caucasian                                5126 81.0%
## 2 Asian                                              682 10.8%
## 3 Black or African American                          553 8.7% 
## 4 Other                                              176 2.8% 
## 5 American Indian/Native American or Alaska Native   130 2.1% 
## 6 Native Hawaiian or Other Pacific Islander           23 0.4%

cat("\nPrimeras filas de race_dummies:\n")

## 
## Primeras filas de race_dummies:

print(head(race_dummies, 5))

## # A tibble: 5 × 7
##   response_id       race_black_or_african_am…¹ race_white_or_caucas…² race_asian
##   <chr>                                  <int>                  <int>      <int>
## 1 R_1ou69fj4DQGsVcp                          1                      0          0
## 2 R_24dboHVOzohx1kw                          0                      1          0
## 3 R_2UbJL30HRjK1sdD                          0                      1          0
## 4 R_UPXamGKtmf4RVIZ                          0                      1          0
## 5 R_2dYk5auG9Fv5Qve                          0                      1          0
## # ℹ abbreviated names: ¹race_black_or_african_american,
## #   ²race_white_or_caucasian
## # ℹ 3 more variables: race_other <int>,
## #   race_american_indian_native_american_or_alaska_native <int>,
## #   race_native_hawaiian_or_other_pacific_islander <int>

✅ Output esperado

Dimensiones de race_long (mayor que survey por las multi-respuestas)
6 columnas binarias en race_dummies
Distribución de razas con conteos y porcentajes

🔢 Bloque 6 — Procesamiento de Eventos de Vida (Multi-respuesta → Dummies)

Explicación Teórica

Este bloque aplica exactamente el mismo pipeline que el bloque anterior, pero para la pregunta de cambios de vida en 2021. La lógica es idéntica porque ambas variables son de selección múltiple.

Relevancia teórica para el modelo

Los eventos de vida tienen un rol particular en el componente de intensidad del modelo POGIT (\(\lambda_i\), SioW). La Sección 8 del paper reporta uno de los hallazgos más interesantes:

“El evento ‘quedar embarazada’ tiene un efecto negativo significativo sobre \(\lambda_i\)”

Esto es intuitivamente consistente: un evento de vida tan significativo como un embarazo implica reorientación del gasto hacia otras categorías (bebé, salud, hogar), reduciendo el gasto en tecnología/electrónica.

Los 5 eventos de vida

Variable	Efecto esperado en \(\lambda_i\)	Justificación
`life_lost_a_job`	Negativo	Reducción de capacidad de gasto
`life_moved_place_of_residence`	Ambiguo	Puede generar compras tech (setup)
`life_divorce`	Negativo	Reorganización financiera
`life_had_a_child`	Negativo	Reorientación del gasto
`life_became_pregnant`	Negativo (significativo)	Confirmado en Sección 8

life_levels <- c(
  "Lost a job",
  "Moved place of residence",
  "Divorce",
  "Had a child",
  "Became pregnant"
)

# Paso 1: Formato largo
life_long <- survey %>%
  dplyr::select(response_id, life_raw = q_life_changes) %>%
  dplyr::filter(!is.na(life_raw) & life_raw != "") %>%
  dplyr::mutate(
    life_raw = life_raw %>%
      stringr::str_replace_all("\\s*/\\s*", "/") %>%
      stringr::str_replace_all("\\s*,\\s*", ", ") %>%
      stringr::str_squish()
  ) %>%
  tidyr::separate_rows(life_raw, sep = ",\\s*") %>%
  dplyr::mutate(life_raw = stringr::str_trim(life_raw))

# Paso 2: One-hot encoding
life_dummies <- life_long %>%
  dplyr::mutate(
    life  = factor(life_raw, levels = life_levels),
    value = 1L
  ) %>%
  dplyr::select(-life_raw) %>%
  tidyr::pivot_wider(
    names_from   = life,
    values_from  = value,
    names_prefix = "life_",
    values_fill  = 0L
  ) %>%
  janitor::clean_names()

# Garantizar que todas las columnas de eventos existen
needed_life <- paste0("life_", janitor::make_clean_names(life_levels))
life_dummies[setdiff(needed_life, names(life_dummies))] <- 0L

# Verificación
cat("============================================================\n")

## ============================================================

cat("LIFE DUMMIES — VERIFICACIÓN\n")

## LIFE DUMMIES — VERIFICACIÓN

cat("============================================================\n")

## ============================================================

cat("Filas en life_long (una por evento declarado):", nrow(life_long), "\n")

## Filas en life_long (una por evento declarado): 2493

cat("Filas en life_dummies (una por participante):", nrow(life_dummies), "\n")

## Filas en life_dummies (una por participante): 2012

cat("Columnas generadas:\n")

## Columnas generadas:

print(names(life_dummies))

## [1] "response_id"                   "life_lost_a_job"              
## [3] "life_moved_place_of_residence" "life_divorce"                 
## [5] "life_had_a_child"              "life_became_pregnant"

cat("\nDistribución por evento de vida:\n")

## 
## Distribución por evento de vida:

print(
  life_long %>%
    dplyr::count(life_raw, sort = TRUE) %>%
    dplyr::mutate(pct = scales::percent(n / nrow(survey), accuracy = 0.1))
)

## # A tibble: 5 × 3
##   life_raw                     n pct  
##   <chr>                    <int> <chr>
## 1 Moved place of residence  1330 21.0%
## 2 Lost a job                 713 11.3%
## 3 Had a child                193 3.1% 
## 4 Became pregnant            183 2.9% 
## 5 Divorce                     74 1.2%

cat("\nPrimeras filas de life_dummies:\n")

## 
## Primeras filas de life_dummies:

print(head(life_dummies, 5))

## # A tibble: 5 × 6
##   response_id       life_lost_a_job life_moved_place_of_residence life_divorce
##   <chr>                       <int>                         <int>        <int>
## 1 R_1ou69fj4DQGsVcp               1                             0            0
## 2 R_1LBvKYfg8hgoloJ               0                             1            0
## 3 R_1kUFSbSAbc8sEME               0                             0            1
## 4 R_1OK73VTU4KtCpkI               0                             0            0
## 5 R_11Y3ZtLQrhccXC7               1                             0            0
## # ℹ 2 more variables: life_had_a_child <int>, life_became_pregnant <int>

✅ Output esperado

5 columnas binarias de eventos de vida
Distribución de eventos con la frecuencia de cada uno
Nota: life_became_pregnant es la variable con efecto negativo significativo en \(\lambda_i\)

🔗 Bloque 7 — Merge Final: `survey_enriched`

Explicación Teórica

Este es el bloque más importante del script. Aquí se integran las tres fuentes de información en un único dataset a nivel de cliente:

\[\underbrace{\text{survey}}_{\text{respuestas originales}} + \underbrace{\text{race\_dummies}}_{\text{6 columnas binarias}} + \underbrace{\text{life\_dummies}}_{\text{5 columnas binarias}} = \underbrace{\text{survey\_enriched}}_{\text{matriz de covariables del encuestado}}\]

¿Por qué `distinct()` antes del join?

Antes de hacer el join, se aplica distinct(response_id, .keep_all = TRUE) a los tres dataframes. Esto garantiza que no haya response_id duplicados, lo que convertiría el join en una relación many-to-many produciendo filas duplicadas silenciosamente. El argumento relationship = "one-to-one" en el left_join() es una verificación adicional: si hubiera duplicados, R lanzaría un error en lugar de continuar con datos incorrectos.

Estructura del dataset final

El survey_enriched tiene la siguiente estructura conceptual:

Columna	Tipo	Descripción
`response_id`	ID	Identificador único del participante
`q_demos_age`	Categórica	Grupo etario (se convertirá en dummies en script 2)
`q_demos_income`	Categórica	Rango de ingreso
`q_amazon_use_howmany`	Categórica	Personas que comparten cuenta
`race_white_or_caucasian`	Binaria	Dummy de raza (ya procesada)
…	Binaria	Otras 5 dummies de raza
`life_became_pregnant`	Binaria	Dummy de evento de vida (ya procesada)
…	Binaria	Otros 4 eventos de vida

Conexión con el paper (Sección 3): Este dataset es la “customer-level matrix” descrita en el paper. Cada fila es un cliente, cada columna es una covariable. Esta es exactamente la estructura que el modelo POGIT necesita para estimar los vectores \(\beta_1\) y \(\beta_2\).

# Garantizar unicidad antes del merge
survey       <- survey       %>% dplyr::distinct(response_id, .keep_all = TRUE)
race_dummies <- race_dummies %>% dplyr::distinct(response_id, .keep_all = TRUE)
life_dummies <- life_dummies %>% dplyr::distinct(response_id, .keep_all = TRUE)

# Merge en dos pasos: survey + race + life
survey_enriched <- survey %>%
  dplyr::left_join(race_dummies, by = "response_id", relationship = "one-to-one") %>%
  dplyr::left_join(life_dummies, by = "response_id", relationship = "one-to-one")

# Verificación completa del dataset final
cat("============================================================\n")

## ============================================================

cat("SURVEY_ENRICHED — OUTPUT FINAL DEL SCRIPT 1\n")

## SURVEY_ENRICHED — OUTPUT FINAL DEL SCRIPT 1

cat("============================================================\n")

## ============================================================

cat("Dimensiones:", nrow(survey_enriched), "filas x", ncol(survey_enriched), "columnas\n\n")

## Dimensiones: 6325 filas x 52 columnas

cat("Columnas de raza (dummies):\n")

## Columnas de raza (dummies):

print(names(survey_enriched)[stringr::str_detect(names(survey_enriched), "^race_")])

## [1] "race_black_or_african_american"                       
## [2] "race_white_or_caucasian"                              
## [3] "race_asian"                                           
## [4] "race_other"                                           
## [5] "race_american_indian_native_american_or_alaska_native"
## [6] "race_native_hawaiian_or_other_pacific_islander"

cat("\nColumnas de eventos de vida (dummies):\n")

## 
## Columnas de eventos de vida (dummies):

print(names(survey_enriched)[stringr::str_detect(names(survey_enriched), "^life_")])

## [1] "life_lost_a_job"               "life_moved_place_of_residence"
## [3] "life_divorce"                  "life_had_a_child"             
## [5] "life_became_pregnant"

cat("\nVerificación de NAs en dummies de raza:\n")

## 
## Verificación de NAs en dummies de raza:

race_cols <- names(survey_enriched)[stringr::str_detect(names(survey_enriched), "^race_")]
print(colSums(is.na(survey_enriched[race_cols])))

##                        race_black_or_african_american 
##                                                     0 
##                               race_white_or_caucasian 
##                                                     0 
##                                            race_asian 
##                                                     0 
##                                            race_other 
##                                                     0 
## race_american_indian_native_american_or_alaska_native 
##                                                     0 
##        race_native_hawaiian_or_other_pacific_islander 
##                                                     0

cat("\nVerificación de NAs en dummies de vida:\n")

## 
## Verificación de NAs en dummies de vida:

life_cols <- names(survey_enriched)[stringr::str_detect(names(survey_enriched), "^life_")]
print(colSums(is.na(survey_enriched[life_cols])))

##               life_lost_a_job life_moved_place_of_residence 
##                          4313                          4313 
##                  life_divorce              life_had_a_child 
##                          4313                          4313 
##          life_became_pregnant 
##                          4313

cat("\nResumen de las columnas más relevantes:\n")

## 
## Resumen de las columnas más relevantes:

survey_enriched %>%
  dplyr::select(dplyr::all_of(c(race_cols, life_cols))) %>%
  dplyr::summarise(dplyr::across(dplyr::everything(), sum)) %>%
  tidyr::pivot_longer(dplyr::everything(),
                      names_to = "variable",
                      values_to = "n_positivos") %>%
  dplyr::mutate(pct = scales::percent(n_positivos / nrow(survey_enriched), accuracy = 0.1)) %>%
  print(n = Inf)

## # A tibble: 11 × 3
##    variable                                              n_positivos pct  
##    <chr>                                                       <int> <chr>
##  1 race_black_or_african_american                                553 8.7% 
##  2 race_white_or_caucasian                                      5126 81.0%
##  3 race_asian                                                    682 10.8%
##  4 race_other                                                    176 2.8% 
##  5 race_american_indian_native_american_or_alaska_native         130 2.1% 
##  6 race_native_hawaiian_or_other_pacific_islander                 23 0.4% 
##  7 life_lost_a_job                                                NA <NA> 
##  8 life_moved_place_of_residence                                  NA <NA> 
##  9 life_divorce                                                   NA <NA> 
## 10 life_had_a_child                                               NA <NA> 
## 11 life_became_pregnant                                           NA <NA>

✅ Output esperado

Dimensiones del dataset final
Lista de columnas binarias de raza y eventos de vida
Conteo de valores positivos por dummy (cuántos participantes pertenecen a cada grupo)
Verificación de que no hay NAs en las dummies

📊 Bloque 8 — Gráficos Descriptivos del Survey

Explicación Teórica

Este bloque genera 20 gráficos de barras horizontales, uno por variable del survey. Corresponde al análisis exploratorio de datos (EDA) de la Sección 2 del paper, donde se describe la composición de la muestra.

¿Por qué es importante el EDA antes del modelado?

Valida la representatividad: Permite verificar que la muestra tiene distribución razonable en variables clave (no todos hombres, no todos de un solo estado, etc.)
Detecta problemas: Una variable con 99% de respuestas en una sola categoría es una señal de que podría ser near-zero-variance y debería considerarse para eliminación en el preprocesamiento.
Contextualiza los resultados: Cuando el modelo reporta que “income_100_149k” tiene efecto positivo en \(\lambda_i\), el gráfico de distribución de ingresos nos dice qué porcentaje de la muestra está en ese rango.

Nota importante: Estos gráficos no entran al modelo. Son herramientas de comprensión de la muestra. Sin embargo, las distribuciones que muestran son exactamente las que determinan qué dummies tendrán suficiente varianza para ser útiles como predictores.

Variables visualizadas

Gráfico	Variable	Conexión con el modelo
Edad	`q_demos_age`	Dummy en \(x_{1i}\)
Raza	`q_demos_race` / `race_raw`	Dummy en \(x_{1i}\)
Ingresos	`q_demos_income`	Dummy en \(x_{2i}\) — significativo en SioW
Cuenta compartida	`q_amazon_use_howmany`	Dummy en \(x_{2i}\) — efecto positivo fuerte
Eventos de vida	`q_life_changes` / `life_raw`	Dummy en \(x_{2i}\) — embarazo negativo
Educación	`q_demos_education`	Dummy en \(x_{1i}\)

# Gráficos ejecutados individualmente (se muestran en el documento)
bar_plot(survey, q_prolific_mturk, "Participation in Amazon MTurk", "Response")

bar_plot(survey, q_demos_age, "Age group distribution", "Age group")

bar_plot(survey, q_demos_hispanic, "Hispanic/Latino origin", "Response",
         fill_pal = pal[c("azul_1", "naranja")])

bar_plot(survey, q_demos_race, "Race combinations (declared)", "Combination",
         text_size = 2.8)

bar_plot(race_long, race_raw, "Participants by race (individual count)", "Race")

bar_plot(survey, q_demos_education, "Education level", "Education", text_size = 2.8)

bar_plot(survey, q_demos_gender, "Declared gender identity", "Gender")

bar_plot(survey, q_sexual_orientation, "Declared sexual orientation", "Orientation")

bar_plot(survey, q_demos_state, "State distribution", "State", text_size = 2.4)

bar_plot(survey, q_amazon_use_howmany, "People sharing Amazon account", "Number of people")

bar_plot(survey, q_amazon_use_hh_size, "Household size", "People in household")

bar_plot(survey, q_amazon_use_how_oft, "Amazon order frequency", "Frequency")

bar_plot(survey, q_substance_use_1, "Household cigarette use", "Response")

bar_plot(survey, q_substance_use_2, "Household marijuana use", "Response")

bar_plot(survey, q_substance_use_3, "Household alcohol use", "Response")

bar_plot(survey, q_personal_1, "Diabetes in household", "Response")

bar_plot(survey, q_personal_2, "Wheelchair use in household", "Response")

bar_plot(survey, q_life_changes, "Life changes in 2021 (combinations)", "Combination",
         text_size = 2.8)

bar_plot(life_long, life_raw, "Life changes in 2021 (individual)", "Life change")

bar_plot(survey, q_demos_income, "Income range", "Response")

✅ Output esperado

20 gráficos de barras horizontales, uno por variable. Cada gráfico muestra barras ordenadas de mayor a menor frecuencia con porcentajes etiquetados.

💾 Bloque 9 — Exportación a PDF

Explicación Teórica

Este bloque exporta todos los gráficos a un único archivo PDF de 20 páginas. La estrategia de dos pasos (primero mostrar en pantalla, luego guardar en lista) es un patrón estándar en R:

Bloque 8 → exploración interactiva (imprime en el viewer de RStudio)
Este bloque → documentación permanente (guarda en disco)

El flujo pdf() → loop print() → dev.off() es necesario porque los objetos ggplot dentro de un loop no se imprimen automáticamente a menos que se llame explícitamente a print(). dev.off() cierra el dispositivo gráfico y escribe físicamente el archivo en disco.

# Recolectar todos los gráficos en una lista
plots_survey <- list(
  bar_plot(survey, q_prolific_mturk, "Participation in Amazon MTurk", "Response"),
  bar_plot(survey, q_demos_age, "Age group distribution", "Age group"),
  bar_plot(survey, q_demos_hispanic, "Hispanic/Latino origin", "Response",
           fill_pal = pal[c("azul_1", "naranja")]),
  bar_plot(survey, q_demos_race, "Race combinations (declared)", "Combination",
           text_size = 2.8),
  bar_plot(race_long, race_raw, "Participants by race (individual count)", "Race"),
  bar_plot(survey, q_demos_education, "Education level", "Education", text_size = 2.8),
  bar_plot(survey, q_demos_gender, "Declared gender identity", "Gender"),
  bar_plot(survey, q_sexual_orientation, "Declared sexual orientation", "Orientation"),
  bar_plot(survey, q_demos_state, "State distribution", "State", text_size = 2.4),
  bar_plot(survey, q_amazon_use_howmany, "People sharing Amazon account", "Number of people"),
  bar_plot(survey, q_amazon_use_hh_size, "Household size", "People in household"),
  bar_plot(survey, q_amazon_use_how_oft, "Amazon order frequency", "Frequency"),
  bar_plot(survey, q_substance_use_1, "Household cigarette use", "Response"),
  bar_plot(survey, q_substance_use_2, "Household marijuana use", "Response"),
  bar_plot(survey, q_substance_use_3, "Household alcohol use", "Response"),
  bar_plot(survey, q_personal_1, "Diabetes in household", "Response"),
  bar_plot(survey, q_personal_2, "Wheelchair use in household", "Response"),
  bar_plot(survey, q_life_changes, "Life changes in 2021 (combinations)", "Combination",
           text_size = 2.8),
  bar_plot(life_long, life_raw, "Life changes in 2021 (individual)", "Life change"),
  bar_plot(survey, q_demos_income, "Income range", "Response")
)

# Exportar a PDF
out_pdf <- "survey_barplots.pdf"
grDevices::pdf(out_pdf, width = 10, height = 7)
for (p in plots_survey) print(p)
grDevices::dev.off()
message("Saved: ", out_pdf)

cat("============================================================\n")
cat("EXPORTACIÓN COMPLETADA\n")
cat("============================================================\n")
cat("Archivo generado:", out_pdf, "\n")
cat("Número de gráficos:", length(plots_survey), "\n")
cat("Dimensiones por página: 10 x 7 pulgadas\n")

✅ Output esperado

Mensaje “Saved: survey_barplots.pdf”
Archivo survey_barplots.pdf con 20 páginas en el directorio de trabajo

🏁 Resumen Final del Script 1

Lo que se construyó

cat("============================================================\n")

## ============================================================

cat("RESUMEN FINAL — SCRIPT 1\n")

## RESUMEN FINAL — SCRIPT 1

cat("============================================================\n\n")

## ============================================================

cat("📦 Paquetes cargados:", length(packages), "\n\n")

## 📦 Paquetes cargados: 30

cat("📂 Dataset principal:\n")

## 📂 Dataset principal:

cat("   survey_enriched:", nrow(survey_enriched), "filas x",
    ncol(survey_enriched), "columnas\n\n")

##    survey_enriched: 6325 filas x 52 columnas

cat("🔢 Variables dummy creadas:\n")

## 🔢 Variables dummy creadas:

cat("   Race dummies:", sum(str_detect(names(survey_enriched), "^race_")), "columnas\n")

##    Race dummies: 6 columnas

cat("   Life dummies:", sum(str_detect(names(survey_enriched), "^life_")), "columnas\n\n")

##    Life dummies: 5 columnas

cat("📊 Gráficos generados:", length(plots_survey), "\n\n")

## 📊 Gráficos generados: 20

cat("💾 Archivo exportado: survey_barplots.pdf\n\n")

## 💾 Archivo exportado: survey_barplots.pdf

cat("➡️  OUTPUT CLAVE PARA SCRIPT 2:\n")

## ➡️  OUTPUT CLAVE PARA SCRIPT 2:

cat("   survey_enriched → fuente de covariables x1 y x2 del modelo POGIT\n")

##    survey_enriched → fuente de covariables x1 y x2 del modelo POGIT

cat("============================================================\n")

## ============================================================

Conexión con el Script 2

El objeto survey_enriched que se creó en este script es el punto de entrada del script 2. Allí se integrará con los datos transaccionales de Amazon (1.8M+ compras) para construir las matrices finales del modelo:

\[\underbrace{\text{survey\_enriched}}_{\text{Script 1}} + \underbrace{\text{amazon-purchases.csv}}_{\text{Script 2}} \longrightarrow \underbrace{V, W, Y}_{\text{Matrices del modelo POGIT}}\]

Documentación generada con R Markdown para reproducibilidad académica. Conecta con: Secciones 2 y 3 del paper — Data & Construction of Predictor Variables.

Code 1 — Survey Data Processing & Exploratory Analysis

POGIT Model Pipeline | Master’s Thesis Documentation

Documentation for Academic Reproducibility

2026-04-07

🔷 Visión General del Script

🗺️ Conexión con el Paper

📐 Lógica del Pipeline

📦 Bloque 1 — Carga de Paquetes

Explicación Teórica

Paquetes por función en el modelo POGIT

✅ Output esperado

🎨 Bloque 2 — Configuración Visual y Paleta de Colores

Explicación Teórica

✅ Output esperado

📊 Bloque 3 — Función bar_plot() y Funciones Auxiliares

Explicación Teórica

La función bar_plot()

La función mode_safe()

La función reorder_drop_id()

✅ Output esperado

📂 Bloque 4 — Carga del Survey

Explicación Teórica

Origen de los datos

¿Por qué clean_names()?

Variables clave del survey

✅ Output esperado

🔢 Bloque 5 — Procesamiento de Raza (Multi-respuesta → Dummies)

Explicación Teórica

El problema de las variables multi-respuesta

La solución: Expand → Encode

¿Por qué NO se elimina una categoría de referencia aquí?

Las 6 categorías raciales

✅ Output esperado

🔢 Bloque 6 — Procesamiento de Eventos de Vida (Multi-respuesta → Dummies)

Explicación Teórica

Relevancia teórica para el modelo

Los 5 eventos de vida

✅ Output esperado

🔗 Bloque 7 — Merge Final: survey_enriched

Explicación Teórica

¿Por qué distinct() antes del join?

Estructura del dataset final

✅ Output esperado

📊 Bloque 8 — Gráficos Descriptivos del Survey

Explicación Teórica

¿Por qué es importante el EDA antes del modelado?

Variables visualizadas

✅ Output esperado

💾 Bloque 9 — Exportación a PDF

Explicación Teórica

✅ Output esperado

🏁 Resumen Final del Script 1

Lo que se construyó

Conexión con el Script 2

📊 Bloque 3 — Función `bar_plot()` y Funciones Auxiliares

La función `bar_plot()`

La función `mode_safe()`

La función `reorder_drop_id()`

¿Por qué `clean_names()`?

🔗 Bloque 7 — Merge Final: `survey_enriched`

¿Por qué `distinct()` antes del join?