Pipeline de Modelos Estadísticos (POGIT y benchmarks) — Adidas | Script 4: Estimación y comparación de modelos
Tesis de Maestría
22 de junio de 2026
Sobre este documento Script 4: estimación del modelo POGIT (modelo principal) y tres benchmarks estadísticos (Poisson, Binomial Negativa, ZIP), siguiendo el plan de menor riesgo acordado. Referencias:
- POGIT: Gutiérrez, I., Ramírez, S., & Jofré, L. (2026). Efficient EM Estimation for the Pogit Model via Polya-Gamma Augmentation. Entropy, 28(2), 207. https://doi.org/10.3390/e28020207 — verificar con tu directora si esta es la referencia exacta que ella entregó, o si hay una versión distinta/anterior.
- Poisson / Binomial Negativa: Cameron, A. C., & Trivedi, P. K. (2013). Regression Analysis of Count Data (2nd ed.). Cambridge University Press.
- ZIP: Lambert, D. (1992). Zero-Inflated Poisson Regression, With an Application to Defects in Manufacturing. Technometrics, 34(1), 1–14.
hurdle()/zeroinfl()en R: Zeileis, A., Kleiber, C., & Jackman, S. (2008). Regression Models for Count Data in R. Journal of Statistical Software, 27(8), 1–25.Pre-requisito: este script asume que el Script 2 ya corrió completo y que
data/matrices/matrices_data_adidas.xlsxtiene las 8 hojas (ids_train,ids_test,V_train,V_test,W_train,W_test,Y_train,Y_test).
1 Bloque 1 — Carga de paquetes
packages <- c(
"readxl", "dplyr", "tibble", "janitor",
"optimParallel", "numDeriv", "parallel",
"MASS", "pscl", "AER", "car", "Metrics",
"ggplot2", "scales", "sysfonts", "showtext"
)
invisible(lapply(packages, require, character.only = TRUE))## Cargando paquete requerido: readxl## Cargando paquete requerido: dplyr##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Cargando paquete requerido: tibble## Cargando paquete requerido: janitor##
## Adjuntando el paquete: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test## Cargando paquete requerido: optimParallel## Cargando paquete requerido: parallel## Cargando paquete requerido: numDeriv## Cargando paquete requerido: MASS##
## Adjuntando el paquete: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Cargando paquete requerido: pscl## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.## Cargando paquete requerido: AER## Cargando paquete requerido: car## Cargando paquete requerido: carData##
## Adjuntando el paquete: 'car'## The following object is masked from 'package:dplyr':
##
## recode## Cargando paquete requerido: lmtest## Cargando paquete requerido: zoo##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Cargando paquete requerido: sandwich## Cargando paquete requerido: survival## Cargando paquete requerido: Metrics## Cargando paquete requerido: ggplot2## Cargando paquete requerido: scales## Cargando paquete requerido: sysfonts## Cargando paquete requerido: showtext## Cargando paquete requerido: showtextdbOutput esperado: Ningún output visible si todos los paquetes están instalados. Si falta alguno (probablemente
AERocar, que no se usaron en los Scripts 1-2), instálalo coninstall.packages("AER")/install.packages("car").
2 Bloque 1B — Helpers de
exportación gráfica (Times New Roman vía showtext)
2.1 ¿Qué hace este bloque y por qué es necesario?
⚠️ NUEVO: este script solo tiene un gráfico (Bloque 3 — distribución de
Y_train), pero por consistencia visual con el resto del proyecto (Scripts 1, 2 y 3), se le aplica el mismo estándar: fuente Times New Roman exacta víashowtext, tema conbase_size = 16, y exportación individual a PDF tamaño carta.
Igual que en los Scripts 2 y 3, si el archivo de fuente no existe en
la ruta esperada, el bloque cae de forma segura a "serif" y
emite una advertencia, en vez de detener la ejecución.
📌 Supuesto a confirmar: no se especificó una carpeta de exportación para este script, así que se usa por defecto
Graficas_Adidas_Estimacion_Modelos, siguiendo el mismo patrón de nombres queGraficas_Adidas_Analisis_Descriptivo(Script 3). Si prefieres otra ruta, dímelo y la ajusto.
ruta_times_ttf <- "C:/Windows/Fonts/times.ttf"
FUENTE_EXPORT <- "Times New Roman"
if (file.exists(ruta_times_ttf)) {
sysfonts::font_add(family = FUENTE_EXPORT, regular = ruta_times_ttf)
showtext::showtext_auto()
showtext::showtext_opts(dpi = 300)
} else {
warning(
"No se encontr\u00f3 el archivo de fuente en ", ruta_times_ttf,
"; las gr\u00e1ficas exportadas usar\u00e1n 'serif' como respaldo en vez de ",
"'Times New Roman' exacta."
)
FUENTE_EXPORT <- "serif"
}
# [NUEVO] Carpeta de salida para las gr\u00e1ficas de este script (Estimaci\u00f3n de modelos).
# Supuesto: nombre an\u00e1logo a la carpeta del Script 3 descriptivo; ajustar si se prefiere otra ruta.
dir_export_modelos <- "C:/Users/User/Documents/Maestria Ciencia de Datos - PUJ Cali/Proyecto_Maestria_SOW_Addidas/cod_grupo2 2/cod_grupo2/Graficas_Adidas_Estimacion_Modelos"
if (!dir.exists(dir_export_modelos)) dir.create(dir_export_modelos, recursive = TRUE)
CARTA_W <- 8.5
CARTA_H <- 11
tamano_fuente <- function(n_elementos) {
if (n_elementos > 6) 16 else 17
}
tema_export_final <- function(base_size = 17) {
ggplot2::theme_minimal(base_size = base_size) +
ggplot2::theme(
text = ggplot2::element_text(family = FUENTE_EXPORT),
panel.grid = ggplot2::element_blank(),
panel.background = ggplot2::element_blank(),
plot.background = ggplot2::element_blank(),
axis.line = ggplot2::element_line(color = "grey60"),
axis.ticks = ggplot2::element_line(color = "grey60"),
axis.text = ggplot2::element_text(color = "grey20"),
axis.title = ggplot2::element_text(color = "grey20"),
plot.title = ggplot2::element_text(color = "grey20", lineheight = 1.05),
plot.margin = ggplot2::margin(t = 18, r = 18, b = 18, l = 18)
)
}
guardar_pdf_carta_gg <- function(plot, nombre_archivo, n_elementos = 1,
orientacion = c("vertical", "horizontal")) {
orientacion <- match.arg(orientacion)
bs <- tamano_fuente(n_elementos)
plot_final <- plot +
tema_export_final(base_size = bs) +
ggplot2::theme(plot.margin = ggplot2::margin(t = 20, r = 22, b = 20, l = 20))
w <- if (orientacion == "vertical") CARTA_W else CARTA_H
h <- if (orientacion == "vertical") CARTA_H else CARTA_W
ruta <- file.path(dir_export_modelos, nombre_archivo)
ggplot2::ggsave(
filename = ruta,
plot = plot_final,
width = w, height = h, units = "in",
dpi = 300,
device = grDevices::cairo_pdf,
limitsize = FALSE
)
message("Guardado: ", ruta)
invisible(ruta)
}
# Tema para visualizaci\u00f3n EN PANTALLA (serif, base_size=16 \u2014 mismo est\u00e1ndar
# que tema_descriptivo en el Script 3), distinto del tema de exportaci\u00f3n final.
tema_descriptivo <- ggplot2::theme_minimal(base_size = 16) +
ggplot2::theme(text = ggplot2::element_text(family = "serif"))Output esperado: Ningún output visible salvo, posiblemente, el
warningde fuente no encontrada.
3 Bloque 2 — Carga de matrices y limpieza de nombres de columna
3.1 ¿Qué hace este bloque y por qué es necesario?
Lee las 8 hojas del Excel del Script 2. Inmediatamente después,
aplica janitor::clean_names() a V y W (train y test) — esto
es importante porque fastDummies::dummy_cols() (Script 2,
Bloque 14) genera nombres de columna como
q_demos_income_$50,000 - $74,999 (con espacios, comas y
signos de dólar), que no son nombres válidos de R. Sin
limpiarlos, cualquier fórmula que construyamos más adelante
(Y ~ ., o las fórmulas de dos partes de
zeroinfl()) fallaría o requeriría escribir cada nombre
entre comillas invertidas. Limpiándolos una sola vez aquí, todo el resto
del script queda libre de ese problema.
clean_names() se aplica por separado a
train y test, pero como es una función determinística del texto del
nombre (no de los datos), produce exactamente los mismos nombres limpios
en ambos — no hay riesgo de que las columnas queden desalineadas entre
train y test.
# [AJUSTADO] Ruta absoluta — mismo criterio que el Script 3 (descriptivo): knitr fija
# el directorio de trabajo a la carpeta del .Rmd, no a la ra\u00edz del proyecto, lo cual
# rompe rutas relativas si el archivo no est\u00e1 guardado exactamente en esa carpeta.
ruta_excel <- "C:/Users/User/Documents/Maestria Ciencia de Datos - PUJ Cali/Proyecto_Maestria_SOW_Addidas/cod_grupo2 2/cod_grupo2/data/matrices/matrices_data_adidas.xlsx"
ids_train <- readxl::read_excel(ruta_excel, sheet = "ids_train")$response_id
ids_test <- readxl::read_excel(ruta_excel, sheet = "ids_test")$response_id
V_train <- readxl::read_excel(ruta_excel, sheet = "V_train") |> as.data.frame() |> janitor::clean_names()
V_test <- readxl::read_excel(ruta_excel, sheet = "V_test") |> as.data.frame() |> janitor::clean_names()
W_train <- readxl::read_excel(ruta_excel, sheet = "W_train") |> as.data.frame() |> janitor::clean_names()
W_test <- readxl::read_excel(ruta_excel, sheet = "W_test") |> as.data.frame() |> janitor::clean_names()
Y_train <- readxl::read_excel(ruta_excel, sheet = "Y_train")[[1]]
Y_test <- readxl::read_excel(ruta_excel, sheet = "Y_test")[[1]]
# Verificación de alineación
stopifnot(length(ids_train) == nrow(V_train), length(ids_train) == nrow(W_train), length(ids_train) == length(Y_train))
stopifnot(length(ids_test) == nrow(V_test), length(ids_test) == nrow(W_test), length(ids_test) == length(Y_test))
stopifnot(identical(names(V_train), names(V_test)))
stopifnot(identical(names(W_train), names(W_test)))
message("V_train: ", nrow(V_train), " x ", ncol(V_train), " | W_train: ", nrow(W_train), " x ", ncol(W_train))## V_train: 472 x 27 | W_train: 472 x 6## Nombres de V (primeros 5): intercept, race_black_or_african_american, race_white_or_caucasian, race_asian, race_otherOutput esperado: Confirmación de dimensiones, y que los nombres de V_train y V_test sean idénticos (mismo para W). Si algún
stopifnot()falla, el Script 2 no terminó de correr correctamente — revisa el Bloque 28 de ese script.
4 Bloque 3 — Diagnóstico de la distribución de Y (insumo para la Sección 1 de tu documento)
4.1 ¿Qué hace este bloque y por qué es necesario?
Calcula los estadísticos descriptivos de Y_train que
necesitas para justificar la elección de cada benchmark: % de ceros
(para ZIP), y la relación media/varianza (para Binomial Negativa).
⚠️ AJUSTADO: el histograma ahora usa
tema_descriptivo(serif, tamaño 16, igual que en los Scripts 1-3) en vez detheme_minimal(base_size = 14), y se exporta de forma individual a PDF (tamaño carta, Times New Roman) enGraficas_Adidas_Estimacion_Modelos.
media_y <- mean(Y_train)
var_y <- var(Y_train)
pct_cero <- mean(Y_train == 0) * 100
resumen_y <- tibble::tibble(
estadistico = c("Media", "Varianza", "Raz\u00f3n Var/Media", "% de ceros", "M\u00e1ximo"),
valor = c(media_y, var_y, var_y / media_y, pct_cero, max(Y_train))
)
print(resumen_y)## # A tibble: 5 × 2
## estadistico valor
## <chr> <dbl>
## 1 Media 0.265
## 2 Varianza 0.569
## 3 Razón Var/Media 2.15
## 4 % de ceros 85.0
## 5 Máximo 4p_dist_y <- ggplot2::ggplot(data.frame(Y = Y_train), ggplot2::aes(x = Y)) +
ggplot2::geom_histogram(binwidth = 1, fill = "grey75", color = "grey40") +
ggplot2::labs(
title = "Distribuci\u00f3n de Y_train (gasto Adidas discretizado, T2)",
x = "Y", y = "Frecuencia"
) +
tema_descriptivo
print(p_dist_y)
# [NUEVO] Exportación individual a PDF, tamaño carta, Times New Roman
guardar_pdf_carta_gg(p_dist_y, "03_distribucion_y_train.pdf", n_elementos = 1)## Guardado: C:/Users/User/Documents/Maestria Ciencia de Datos - PUJ Cali/Proyecto_Maestria_SOW_Addidas/cod_grupo2 2/cod_grupo2/Graficas_Adidas_Estimacion_Modelos/03_distribucion_y_train.pdfOutput esperado: Una tabla con 5 estadísticos y un histograma con pico pronunciado en 0, en pantalla y exportado a
Graficas_Adidas_Estimacion_Modelos. Si “Razón Var/Media” es mucho mayor a 1, confirma sobredispersión (justifica Binomial Negativa). Si “% de ceros” está cerca de 77%, confirma la necesidad de un modelo de inflación de ceros (justifica ZIP).
5 Bloque 4 — Preparación
de train_df/test_df para los benchmarks
5.1 ¿Qué hace este bloque y por qué es necesario?
⚠️ Corrección importante:
V_trainyW_trainya tienen cada una su propia columnaintercept(agregada en el Bloque 28 del Script 2). Si las combinamos sin quitarla,glm()agregaría un tercer intercepto automáticamente (vía la fórmulaY ~ .), generando columnas perfectamente colineales. Por eso se eliminan ambosinterceptmanuales aquí — los benchmarks usan el intercepto queglm()/zeroinfl()ya agregan por su cuenta.
V_train_sb <- dplyr::select(V_train, -intercept)
W_train_sb <- dplyr::select(W_train, -intercept)
V_test_sb <- dplyr::select(V_test, -intercept)
W_test_sb <- dplyr::select(W_test, -intercept)
train_df <- cbind(V_train_sb, W_train_sb)
train_df$Y <- Y_train
test_df <- cbind(V_test_sb, W_test_sb)
test_df$Y <- Y_test
stopifnot(identical(names(train_df), names(test_df)))
message("train_df: ", nrow(train_df), " x ", ncol(train_df), " (incluye columna Y)")## train_df: 472 x 32 (incluye columna Y)Output esperado: Confirmación de dimensiones.
train_dfytest_dfdeben tener exactamente las mismas columnas (mismo orden, mismos nombres).
6 Bloque 5 — Modelo Poisson (Benchmark 1)
6.1 ¿Qué hace este bloque y por qué es necesario?
Ajusta un Poisson de una sola ecuación sobre Y_train,
usando V+W combinadas como predictoras. Es el benchmark más simple: ¿qué
tan bien predice un modelo de conteo estándar, sin la separación teórica
SoW/SioW del POGIT?
pois_model <- glm(Y ~ ., data = train_df, family = poisson(link = "log"))
cat("AIC Poisson:", AIC(pois_model), "| BIC:", BIC(pois_model), "\n")## AIC Poisson: 638.6684 | BIC: 771.6918## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2984960 0.7090144 -1.8314100 0.06703937
## race_black_or_african_american -0.6616411 0.5647182 -1.1716304 0.24134546
## race_white_or_caucasian -0.3888962 0.5680603 -0.6846036 0.49359407
## race_asian -0.5854937 0.5917984 -0.9893466 0.32249359
## race_other -0.2281024 0.7642918 -0.2984494 0.76536019Output esperado: Un modelo ajustado sin errores de rango deficiente. Si ves coeficientes como
NAen el resumen completo (summary(pois_model)), hay colinealidad perfecta residual — revisa el Bloque 4.
7 Bloque 6 — Diagnósticos post-Poisson: VIF y test de sobredispersión
7.1 ¿Qué hace este bloque y por qué es necesario?
VIF (Variance Inflation Factor): complementa la correlación de Spearman que ya calculaste en el Script 2 (Bloque 26). Un VIF mayor a 5-10 para una variable indica que su información es altamente redundante con las demás.
Test de sobredispersión
(AER::dispersiontest): a diferencia del chequeo crudo
var(Y)/mean(Y) del Bloque 3 (que mezcla heterogeneidad
explicada y no explicada), este test evalúa sobredispersión
condicional — después de controlar por las covariables
del modelo Poisson ya ajustado. Es la versión rigurosa del mismo
diagnóstico, y la que deberías citar en la tesis.
vif_vals <- car::vif(pois_model)
vif_altos <- sort(vif_vals[vif_vals > 5], decreasing = TRUE)
if (length(vif_altos) > 0) {
message("Variables con VIF > 5 (multicolinealidad relevante):")
print(vif_altos)
} else {
message("Ninguna variable tiene VIF > 5.")
}## Variables con VIF > 5 (multicolinealidad relevante):## dollar_avg12_top1_ratio race_white_or_caucasian
## 6.146511 5.199223##
## Overdispersion test
##
## data: pois_model
## z = 2.9625, p-value = 0.001526
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
## alpha
## 0.7056799Output esperado: Una lista (posiblemente vacía) de variables con VIF alto, y el resultado del test de dispersión. Si el p-valor del test es < 0.05 y el coeficiente de dispersión estimado es > 1, hay sobredispersión estadísticamente significativa — esa es tu justificación formal para la Binomial Negativa.
8 Bloque 7 — Modelo Binomial Negativa (Benchmark 2)
nb_model <- MASS::glm.nb(Y ~ ., data = train_df)
cat("AIC NegBin:", AIC(nb_model), "| BIC:", BIC(nb_model), "\n")## AIC NegBin: 588.9152 | BIC: 726.0955## Parámetro de dispersión (theta): 0.3554297Output esperado: AIC/BIC, idealmente más bajos que los del Poisson (Bloque 5) si efectivamente había sobredispersión — eso confirmaría que la Binomial Negativa es una mejor especificación que el Poisson simple.
9 Bloque 8 — Modelo ZIP (Benchmark 3)
9.1 ¿Qué hace este bloque y por qué es necesario?
⚠️ Corrección importante: la fórmula de dos partes de
zeroinfl()separa explícitamente el componente de conteo (antes del|, usando las variables de W — comportamiento/intensidad) del componente de inflación de ceros (después del|, usando las variables de V — demografía/asignación). Esto preserva exactamente los mismos roles teóricos que V y W tienen en el POGIT, en vez de meter todas las variables en ambas ecuaciones.
form_zip <- as.formula(paste(
"Y ~", paste(names(W_train_sb), collapse = " + "),
"|", paste(names(V_train_sb), collapse = " + ")
))
zip_model <- pscl::zeroinfl(form_zip, data = train_df, dist = "poisson")
cat("AIC ZIP:", AIC(zip_model), "| BIC:", BIC(zip_model), "\n")## AIC ZIP: 587.8074 | BIC: 724.9877##
## Call:
## pscl::zeroinfl(formula = form_zip, data = train_df, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -0.8161 -0.4080 -0.2971 -0.1894 8.0112
##
## Count model coefficients (poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.20908 0.17396 -1.202 0.22942
## recency_days_adidas -0.48380 0.18536 -2.610 0.00905 **
## frequency_adidas 0.33332 0.14920 2.234 0.02548 *
## monetary_adidas -0.02621 0.13063 -0.201 0.84101
## dollar_avg_12m -0.16531 0.17109 -0.966 0.33394
## dollar_avg12_top1_ratio -0.04461 0.20937 -0.213 0.83128
##
## Zero-inflation model coefficients (binomial with logit link):
## Estimate Std. Error
## (Intercept) 0.04818 1.22875
## race_black_or_african_american 1.29185 0.93441
## race_white_or_caucasian 1.11247 0.94557
## race_asian 1.55103 1.04841
## race_other 0.33078 1.32705
## race_american_indian_native_american_or_alaska_native 0.69457 1.29816
## race_native_hawaiian_or_other_pacific_islander 13.97908 2361.03322
## life_lost_a_job 0.54520 0.64777
## life_moved_place_of_residence 0.93589 0.44283
## life_divorce -1.56362 1.13671
## life_had_a_child -0.76669 0.92937
## life_became_pregnant -0.03031 1.01372
## q_demos_age_25_34_years 0.37434 0.64184
## q_demos_age_35_44_years -0.07152 0.68301
## q_demos_age_45_54_years 0.12294 0.74396
## q_demos_age_55_64_years 0.26548 1.00837
## q_demos_age_65_and_older -0.90590 1.82251
## q_demos_income_50_000_74_999 -0.39102 0.56956
## q_demos_income_75_000_99_999 -0.10477 0.62443
## q_demos_income_100_000_149_999 -0.58537 0.54953
## q_demos_income_150_000_or_more -1.03766 0.59311
## q_demos_income_less_than_25_000 0.95054 0.91012
## q_demos_income_prefer_not_to_say 14.24986 2064.17327
## q_amazon_use_howmany_2 -0.20859 0.37772
## q_amazon_use_howmany_3 -1.00212 0.64379
## q_amazon_use_howmany_4 -0.37400 0.75422
## n_states -0.05271 0.17245
## z value Pr(>|z|)
## (Intercept) 0.039 0.9687
## race_black_or_african_american 1.383 0.1668
## race_white_or_caucasian 1.177 0.2394
## race_asian 1.479 0.1390
## race_other 0.249 0.8032
## race_american_indian_native_american_or_alaska_native 0.535 0.5926
## race_native_hawaiian_or_other_pacific_islander 0.006 0.9953
## life_lost_a_job 0.842 0.4000
## life_moved_place_of_residence 2.113 0.0346 *
## life_divorce -1.376 0.1690
## life_had_a_child -0.825 0.4094
## life_became_pregnant -0.030 0.9762
## q_demos_age_25_34_years 0.583 0.5597
## q_demos_age_35_44_years -0.105 0.9166
## q_demos_age_45_54_years 0.165 0.8687
## q_demos_age_55_64_years 0.263 0.7923
## q_demos_age_65_and_older -0.497 0.6191
## q_demos_income_50_000_74_999 -0.687 0.4924
## q_demos_income_75_000_99_999 -0.168 0.8668
## q_demos_income_100_000_149_999 -1.065 0.2868
## q_demos_income_150_000_or_more -1.750 0.0802 .
## q_demos_income_less_than_25_000 1.044 0.2963
## q_demos_income_prefer_not_to_say 0.007 0.9945
## q_amazon_use_howmany_2 -0.552 0.5808
## q_amazon_use_howmany_3 -1.557 0.1196
## q_amazon_use_howmany_4 -0.496 0.6200
## n_states -0.306 0.7599
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 31
## Log-likelihood: -260.9 on 33 DfOutput esperado: Dos bloques de coeficientes en el
summary(): “Count model coefficients” (las variables de W) y “Zero-inflation model coefficients” (las variables de V). Sizeroinfl()no converge, prueba agregandoEM = TRUEa la llamada, o reduce el número de variables candidatas (podrías aplicar el mismo criterio de Spearman/NZV del Script 2 a W antes de esto).
10 Bloque 9 — Log-verosimilitud del modelo POGIT
10.1 ¿Qué hace este bloque y por qué es necesario?
Implementa la verosimilitud marginal \((y_i
\mid \theta_i,\lambda_i) \sim
\text{Poisson}(E_i\,\theta_i\,\lambda_i)\) de Gutiérrez, Ramírez
y Jofré (2026), con \(\theta_i =
\text{sigmoid}(V_i\beta_1)\) y \(\lambda_i = \exp(W_i\beta_2)\). A
diferencia de los benchmarks, el POGIT sí necesita la
columna intercept manual de V y W (porque la verosimilitud
no la agrega por su cuenta), así que aquí usamos V_train/
W_train completas, no las versiones _sb del
Bloque 4.
neg_loglik_pogit <- function(par, V, W, Y, E = 1) {
kV <- ncol(V); kW <- ncol(W)
beta1 <- par[1:kV]
beta2 <- par[(kV + 1):(kV + kW)]
eta1 <- as.vector(V %*% beta1)
eta2 <- as.vector(W %*% beta2)
theta <- plogis(eta1)
lambda <- exp(eta2)
mu <- E * theta * lambda
mu <- pmax(mu, 1e-10)
ll <- sum(Y * log(mu) - mu)
-ll
}
message("Funci\u00f3n neg_loglik_pogit() definida.")## Función neg_loglik_pogit() definida.11 Bloque 10 — Estimación POGIT por máxima verosimilitud (TRAIN)
V_train_mat <- as.matrix(V_train)
W_train_mat <- as.matrix(W_train)
n_par_V <- ncol(V_train_mat)
n_par_W <- ncol(W_train_mat)
par_init <- rep(0, n_par_V + n_par_W)
cl <- parallel::makeCluster(max(1, parallel::detectCores() - 1))
fit_pogit <- optimParallel::optimParallel(
par = par_init,
fn = neg_loglik_pogit,
V = V_train_mat, W = W_train_mat, Y = Y_train,
method = "BFGS",
parallel = list(cl = cl, forward = FALSE, loginfo = TRUE)
)
parallel::stopCluster(cl)
beta1_hat <- fit_pogit$par[1:n_par_V]
beta2_hat <- fit_pogit$par[(n_par_V + 1):(n_par_V + n_par_W)]
names(beta1_hat) <- colnames(V_train_mat)
names(beta2_hat) <- colnames(W_train_mat)
message("Convergencia (0 = exitosa): ", fit_pogit$convergence)## Convergencia (0 = exitosa): 1## Log-verosimilitud en el óptimo: -237.39Output esperado:
convergence = 0. Si no converge, pruebamethod = "L-BFGS-B"o revisa colinealidad en V concar::vif()análogo al Bloque 6.
12 Bloque 11 — Errores estándar del POGIT (Hessiano)
hess <- numDeriv::hessian(
func = neg_loglik_pogit, x = fit_pogit$par,
V = V_train_mat, W = W_train_mat, Y = Y_train
)
vcov_mat <- solve(hess)
se_all <- sqrt(diag(vcov_mat))
se_beta1 <- se_all[1:n_par_V]
se_beta2 <- se_all[(n_par_V + 1):(n_par_V + n_par_W)]
tabla_beta1 <- tibble::tibble(
variable = names(beta1_hat), beta = beta1_hat, se = se_beta1,
z = beta / se, p_valor = 2 * (1 - pnorm(abs(z)))
)
tabla_beta2 <- tibble::tibble(
variable = names(beta2_hat), beta = beta2_hat, se = se_beta2,
z = beta / se, p_valor = 2 * (1 - pnorm(abs(z)))
)
cat("=== Coeficientes theta (V) ===\n"); print(tabla_beta1)## === Coeficientes theta (V) ===## # A tibble: 27 × 5
## variable beta se z p_valor
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 intercept -1.05 7.10 -0.148 0.882
## 2 race_black_or_african_american -1.02 2.16 -0.473 0.636
## 3 race_white_or_caucasian -0.668 1.47 -0.456 0.648
## 4 race_asian -1.01 2.38 -0.423 0.673
## 5 race_other -0.507 1.58 -0.320 0.749
## 6 race_american_indian_native_american_or_alaska_n… -0.776 1.41 -0.550 0.582
## 7 race_native_hawaiian_or_other_pacific_islander -2.94 4.95 -0.593 0.553
## 8 life_lost_a_job 0.265 0.518 0.512 0.608
## 9 life_moved_place_of_residence -0.430 0.370 -1.16 0.244
## 10 life_divorce 0.847 0.921 0.920 0.358
## # ℹ 17 more rows##
## === Coeficientes lambda (W) ===## # A tibble: 6 × 5
## variable beta se z p_valor
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 intercept 0.243 4.64 0.0525 0.958
## 2 recency_days_adidas -0.347 0.189 -1.83 0.0667
## 3 frequency_adidas 0.430 0.125 3.44 0.000574
## 4 monetary_adidas 0.0553 0.121 0.459 0.646
## 5 dollar_avg_12m -0.0910 0.190 -0.479 0.632
## 6 dollar_avg12_top1_ratio 0.0285 0.186 0.153 0.878Output esperado: Dos tablas con coeficiente, error estándar, z y p-valor.
13 Bloque 12 — Predicciones de los 4 modelos en TEST y tabla comparativa
13.1 ¿Qué hace este bloque y por qué es necesario?
Genera predicciones en TEST para los cuatro modelos y construye la
tabla comparativa final (MAE, RMSE, correlación de Spearman, AIC, BIC).
Para el POGIT, AIC/BIC se calculan a mano porque no es un objeto
glm/zeroinfl con método AIC()
nativo.
pred_pois <- predict(pois_model, newdata = test_df, type = "response")
pred_nb <- predict(nb_model, newdata = test_df, type = "response")
pred_zip <- predict(zip_model, newdata = test_df, type = "response")
V_test_mat <- as.matrix(V_test)
W_test_mat <- as.matrix(W_test)
theta_hat_test <- plogis(as.vector(V_test_mat %*% beta1_hat))
lambda_hat_test <- exp(as.vector(W_test_mat %*% beta2_hat))
pred_pogit <- theta_hat_test * lambda_hat_test # E_i = 1
# AIC/BIC manuales para POGIT
k_pogit <- length(fit_pogit$par)
n_train <- length(Y_train)
loglik_pogit <- -fit_pogit$value
aic_pogit <- 2 * k_pogit - 2 * loglik_pogit
bic_pogit <- k_pogit * log(n_train) - 2 * loglik_pogit
tabla_comparativa <- tibble::tibble(
modelo = c("Poisson", "NegBin", "ZIP", "POGIT"),
MAE = c(Metrics::mae(Y_test, pred_pois), Metrics::mae(Y_test, pred_nb),
Metrics::mae(Y_test, pred_zip), Metrics::mae(Y_test, pred_pogit)),
RMSE = c(Metrics::rmse(Y_test, pred_pois), Metrics::rmse(Y_test, pred_nb),
Metrics::rmse(Y_test, pred_zip), Metrics::rmse(Y_test, pred_pogit)),
Spearman = c(cor(Y_test, pred_pois, method = "spearman"),
cor(Y_test, pred_nb, method = "spearman"),
cor(Y_test, pred_zip, method = "spearman"),
cor(Y_test, pred_pogit, method = "spearman")),
AIC = c(AIC(pois_model), AIC(nb_model), AIC(zip_model), aic_pogit),
BIC = c(BIC(pois_model), BIC(nb_model), BIC(zip_model), bic_pogit)
)
print(tabla_comparativa)## # A tibble: 4 × 6
## modelo MAE RMSE Spearman AIC BIC
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Poisson 0.389 0.714 0.127 639. 772.
## 2 NegBin 0.373 0.680 0.151 589. 726.
## 3 ZIP 0.366 0.644 0.152 588. 725.
## 4 POGIT 0.385 0.700 0.117 541. 678.Output esperado: La tabla final con las 4 filas (Poisson, NegBin, ZIP, POGIT) y las 5 métricas. Esta tabla es directamente la “Sección 3 — Resultados” de tu documento.
14 Bloque 13 — POGIT: SoW/SioW predicho vs. empírico (validación final)
14.1 ¿Qué hace este bloque y por qué es necesario?
Compara \(\hat\theta_i\) (SoW
predicho) contra el SoW empírico exportado en el Script 2
(empirical_sow_adidas.xlsx) — el benchmark real de negocio,
no solo la métrica de ajuste de conteos del Bloque 12.
# [AJUSTADO] Ruta absoluta — mismo criterio que el Bloque 2 y que el Script 3 (descriptivo)
ruta_sow_empirico <- "C:/Users/User/Documents/Maestria Ciencia de Datos - PUJ Cali/Proyecto_Maestria_SOW_Addidas/cod_grupo2 2/cod_grupo2/data/matrices/empirical_sow_adidas.xlsx"
sow_emp <- readxl::read_excel(ruta_sow_empirico, sheet = "sow")[[1]]
ids_sow <- readxl::read_excel(ruta_sow_empirico, sheet = "id")[[1]]
idx_match <- match(ids_test, ids_sow)
stopifnot(!any(is.na(idx_match)))
sow_emp_alineado <- sow_emp[idx_match]
mae_theta <- Metrics::mae(sow_emp_alineado, theta_hat_test)
cor_theta <- cor(sow_emp_alineado, theta_hat_test, method = "spearman")
message("MAE (theta_hat vs SoW emp\u00edrico, TEST): ", round(mae_theta, 4))## MAE (theta_hat vs SoW empírico, TEST): 0.1668## Correlación de Spearman (theta_hat vs SoW empírico): -0.0073Output esperado: Un MAE entre 0 y 1, y una correlación de Spearman — idealmente alta y positiva, indicando que el modelo ordena correctamente a los clientes de mayor a menor SoW real, aunque el valor exacto tenga error.