Predicción del ingreso mensual según trayectoria y cargo: evaluación de diferencias ajustadas por género
- 1 Resultado principal y decisión
- 2 Problema, usuarios y alcance
- 3 Fuente, auditoría y preparación
- 4 Exploración
- 5 Tres regresiones simples
- 6 Modelos múltiples y AIC
- 7 Diagnósticos y transformaciones
- 8 Validación y elección predictiva
- 9 Diferencias ajustadas por género
- 10 Escenarios e intervalos
- 11 Análisis crítico y recomendación
- 12 Referencias y reproducibilidad
1 Resultado principal y decisión
La trayectoria y el cargo permiten mejorar considerablemente la predicción del ingreso dentro de esta base. El procedimiento seleccionado es A1, una regresión obtenida mediante AIC que conserva experiencia, antigüedad y cargo. Su MAE de validación es $1.296.849, frente a $3.378.131 de la referencia que predice la mediana del entrenamiento. La reducción es de 61,6%.
La comparación por género ofrece una conclusión más incierta. En el modelo múltiple con controles definidos previamente, la diferencia estimada de mujeres frente a hombres es de $-170.605 mensuales y su intervalo de confianza robusto incluye cero. La magnitud y la significancia cambian al modificar la especificación o retirar observaciones influyentes. Estos resultados no establecen igualdad salarial ni discriminación.
La aplicación defendible es académica. La base corresponde a una adaptación de registros sintéticos de IBM y el ingreso original fue multiplicado por 1.000. El modelo permite ilustrar cómo revisar referencias de perfiles comparables, pero no sustenta fijar remuneraciones reales ni diagnosticar una empresa colombiana. La utilidad práctica también queda limitada por el error, la información omitida y la amplitud de los intervalos individuales.
2 Problema, usuarios y alcance
El usuario hipotético es un equipo de compensación de Recursos Humanos. La decisión sería evaluar si una referencia salarial propuesta merece una revisión más detallada y estimar órdenes de magnitud para ejercicios presupuestales. Interesan tanto una predicción como su incertidumbre y el costo de equivocarse.
La pregunta principal es qué tan confiable resulta predecir
monthly_income a partir de trayectoria y cargo. La pregunta
complementaria es cómo cambia la diferencia observada por género al
considerar características laborales.
La regresión lineal ofrece una referencia interpretable para aproximar la media del ingreso condicionada a esas características. Su conveniencia depende de validación y diagnósticos. La relación entre experiencia e ingreso es transversal: no representa el aumento que recibiría una misma persona al acumular un año adicional de trabajo.
3 Fuente, auditoría y preparación
3.1 Procedencia y escala
El archivo de entrada es
dataset_03_employee_income_and_job_characteristics.xlsx,
hoja employee_income, con 1470 filas y 24 variables
originales. La hoja Data Definitions contiene las
definiciones educativas y de escalas. No hay fecha de referencia,
ubicación de levantamiento ni diseño muestral para identificar una
población laboral colombiana.
La comparación con el archivo publicado por IBM verifica coincidencia fila a fila de 15 variables numéricas, la traducción de género y el ingreso original multiplicado por 1.000. Posit Academy describe el conjunto IBM como datos sintetizados. Esta identificación combina verificación computacional con documentación del conjunto de referencia.
Se presentan ingresos en pesos, como solicita la guía. El factor 1.000 es una adaptación de escala; no documenta una conversión cambiaria. No se incorporan al modelo variables adicionales del archivo externo, utilizado exclusivamente para verificar procedencia y escala.
tab(data.frame(check=c("Comparaciones de procedencia", "Coincidencias verificadas"),
result=c(nrow(analysis$provenance), sum(analysis$provenance$matches))))| Verificación | Resultado |
|---|---|
| Comparaciones de procedencia | 18 |
| Coincidencias verificadas | 18 |
3.2 Calidad y tratamiento
quality_labels <- c("Filas originales", "Columnas originales", "Celdas faltantes", "Filas idénticas",
"Ingresos no positivos", "Experiencia negativa", "Antigüedad mayor que experiencia",
"Antigüedad en cargo mayor que en empresa", "Experiencia incompatible con inicio desde 14 años",
"Ingresos señalados por la regla de Tukey")
tab(data.frame(check=quality_labels, result=analysis$audit$value))| Verificación | Resultado |
|---|---|
| Filas originales | 1470 |
| Columnas originales | 24 |
| Celdas faltantes | 0 |
| Filas idénticas | 0 |
| Ingresos no positivos | 0 |
| Experiencia negativa | 0 |
| Antigüedad mayor que experiencia | 0 |
| Antigüedad en cargo mayor que en empresa | 0 |
| Experiencia incompatible con inicio desde 14 años | 0 |
| Ingresos señalados por la regla de Tukey | 114 |
Se conservaron las 1.470 filas. No se imputaron valores, corrigieron ingresos ni aplicaron recortes o winsorización. Los 114 ingresos señalados por la regla de 1,5 veces el rango intercuartílico no acreditan errores: la parte alta incluye cargos directivos. Las observaciones influyentes se analizan como sensibilidad, sin retirarlas de la muestra principal. Un control de coherencia de edad y experiencia tampoco acredita la historia real de registros sintéticos.
Se preservaron las columnas y valores del archivo original. En el
análisis se creó record_id, que identifica la fila de
origen y no constituye un identificador de empleado.
education_category es un factor no ordenado derivado de los
códigos educativos. No se asume que pasar entre niveles representa
incrementos educativos o salariales iguales.
Las referencias son hombre (M), sin horas extra
(No), pregrado (bachelor, código 3) e
investigación científica (Investigador_Cientifico).
below_college y college conservan el sentido
del diccionario, sin equipararlos a niveles específicos del sistema
colombiano.
3.3 Muestras y validación
Se reservaron aleatoriamente 294 filas para calibrar intervalos y se utilizaron 1176 para desarrollo. La semilla es 2026. Dentro de desarrollo se crearon 10 particiones equilibradas por género, repetidas cinco veces. Cada modelo genera 5880 predicciones, cinco por cada registro de desarrollo.
La selección por AIC y el factor de retransformación de los logaritmos se estiman dentro de cada entrenamiento. Todos los candidatos se evalúan en los mismos registros y particiones. La muestra de calibración no interviene en la selección predictiva.
Los modelos se reajustan también sobre las 1.470 filas para interpretar coeficientes y realizar diagnósticos. Estos reajustes se distinguen del ajuste de desarrollo que predice los escenarios con intervalos calibrados. La exploración previa en la conversación impide presentar el ejercicio como evaluación externa o completamente prospectiva.
4 Exploración
z <- analysis$descriptive
tab(data.frame(variable=c("Ingreso mensual", "Experiencia (años)", "Antigüedad (años)", "Edad (años)"),
mean=num(z$mean), median=num(z$median), standard_deviation=num(z$sd), minimum=num(z$minimum),
Q1=num(z$q1), Q3=num(z$q3), maximum=num(z$maximum)))| Variable | Media | Mediana | DE | Mínimo | Q1 | Q3 | Máximo |
|---|---|---|---|---|---|---|---|
| Ingreso mensual | 6.502.931,29 | 4.919.000,00 | 4.707.956,78 | 1.009.000,00 | 2.911.000,00 | 8.379.000,00 | 19.999.000,00 |
| Experiencia (años) | 11,28 | 10,00 | 7,78 | 0,00 | 6,00 | 15,00 | 40,00 |
| Antigüedad (años) | 7,01 | 5,00 | 6,13 | 0,00 | 3,00 | 9,00 | 40,00 |
| Edad (años) | 36,92 | 36,00 | 9,14 | 18,00 | 30,00 | 43,00 | 60,00 |
La media del ingreso supera su mediana y la distribución tiene una cola hacia valores altos. Es pertinente comparar con una referencia de mediana al utilizar MAE: una predicción constante basada en la media puede ser menos favorable bajo ese criterio.
Distribución del ingreso en la escala académica de la base.
Asociaciones transversales entre ingreso y trayectoria.
tab(data.frame(Predictor=c("Experiencia", "Antigüedad", "Edad"),
Pearson=num(analysis$correlations$pearson, 3), Spearman=num(analysis$correlations$spearman, 3)))| Predictor | Pearson | Spearman |
|---|---|---|
| Experiencia | 0,773 | 0,710 |
| Antigüedad | 0,514 | 0,464 |
| Edad | 0,498 | 0,472 |
Experiencia tiene una asociación positiva más fuerte que antigüedad. Las diferencias entre las tendencias rectas y locales y la dispersión alrededor de ellas indican que una pendiente única no agota la estructura salarial. Edad, experiencia y antigüedad comparten información; se mantiene el alcance de trayectoria acordado en lugar de incorporar todas las variables sin justificación.
Ingreso por cargo, educación y horas extra.
Ingreso por cargo, educación y horas extra.
Los cargos distinguen distribuciones salariales, especialmente los directivos. Las comparaciones educativas agregadas no separan el papel del cargo o de la experiencia. Registrar horas extra tampoco equivale a conocer su cantidad ni el componente salarial que remuneran.
z <- analysis$group_summaries[analysis$group_summaries$group_variable == "gender", ]
tab(data.frame(gender=z$group, records=z$n, mean=money(z$mean_income), median=money(z$median_income)))| Género | Registros | Media | Mediana |
|---|---|---|---|
| M | 882 | $6.380.508 | $4.837.500 |
| F | 588 | $6.686.566 | $5.081.500 |
Ingreso y composición de los cargos por género.
Ingreso y composición de los cargos por género.
La diferencia inicial de medias favorece a las mujeres. Los grupos no tienen idéntica composición de cargo y trayectoria, de modo que esa diferencia no responde cuánto difieren perfiles comparables.
5 Tres regresiones simples
S1 utiliza experiencia, S2 antigüedad y S3 género binario. Coeficientes y R² corresponden al reajuste descriptivo completo; MAE corresponde a validación en desarrollo.
tab(coefficient_display(c("S1", "S2", "S3")))| Modelo | Término | Coeficiente | IC 95% HC3 | p HC3 |
|---|---|---|---|---|
| S1 | Constante | $1.227.935 | [$1.013.484; $1.442.387] | <0,001 |
| S1 | Experiencia (un año) | $467.658 | [$443.433; $491.884] | <0,001 |
| S2 | Constante | $3.733.273 | [$3.415.209; $4.051.338] | <0,001 |
| S2 | Antigüedad en la empresa (un año) | $395.205 | [$357.605; $432.804] | <0,001 |
| S3 | Constante | $6.380.508 | [$6.068.916; $6.692.100] | <0,001 |
| S3 | Mujer frente a hombre | \(306.058 |[\)-185.490; $797.607] | 0,222 |
z <- analysis$fit_summary[analysis$fit_summary$model %in% c("S1", "S2", "S3"), ]
tab(data.frame(model=z$model, r_squared=num(z$r_squared, 4), hc3_joint_p=p_text(z$hc3_joint_p),
validation_mae=money(vapply(z$model, metric, numeric(1))), check.names=FALSE))| Modelo | R² | p conjunta HC3 | MAE validación |
|---|---|---|---|
| S1 | 0,5974 | <0,001 | $2.212.231 |
| S2 | 0,2645 | <0,001 | $2.955.735 |
| S3 | 0,0010 | 0,222 | $3.663.512 |
En S1, un año adicional de experiencia se asocia con $467.658 de ingreso mensual adicional entre personas de esta base. En S2, la asociación por año de antigüedad es $395.205. Ambas son estadísticamente distinguibles de cero, pero no representan efectos de capacitación o permanencia sobre salarios.
En S3, la constante es el promedio masculino y el coeficiente de mujer es la diferencia de medias, $306.058, con p = 0,222. Su R² es pequeño y su MAE supera la referencia de mediana. Significancia y precisión predictiva son criterios distintos.
S1 es el mejor simple por menor error, coherencia e interpretación. Las constantes de S1 y S2 describen el ajuste cuando el predictor es cero; no son salarios universales de entrada.
6 Modelos múltiples y AIC
6.1 Especificaciones y coeficientes
M1 incluye experiencia, antigüedad, educación, cargo y horas extra. Experiencia representa trayectoria general; antigüedad, vinculación con la empresa. Educación y cargo distinguen perfil y función. Horas extra es una característica disponible, limitada por su codificación binaria. M2 añade género para evaluar su contribución y la diferencia ajustada sin delegar esa pregunta a la selección automática.
z <- analysis$fit_summary[analysis$fit_summary$model %in% c("M1", "M2", "A1", "A2"), ]
tab(data.frame(model=z$model, parameters=z$parameters, adjusted_r_squared=num(z$adjusted_r_squared, 4),
aic=num(z$aic), hc3_joint_p=p_text(z$hc3_joint_p), check.names=FALSE))| Modelo | Parámetros | R² ajustado | AIC | p conjunta HC3 |
|---|---|---|---|---|
| M1 | 16 | 0,8764 | 46.288,62 | <0,001 |
| M2 | 17 | 0,8766 | 46.286,88 | <0,001 |
| A1 | 11 | 0,8762 | 46.286,52 | <0,001 |
| A2 | 12 | 0,8764 | 46.284,90 | <0,001 |
tab(coefficient_display(c("M1", "M2")))| Modelo | Término | Coeficiente | IC 95% HC3 | p HC3 |
|---|---|---|---|---|
| M1 | Constante | $1.733.021 | [$1.527.654; $1.938.387] | <0,001 |
| M1 | Experiencia (un año) | $189.740 | [$168.319; $211.162] | <0,001 |
| M1 | Antigüedad en la empresa (un año) | \(19.945 |[\)-2.872; $42.762] | 0,087 | |
| M1 | Educación: inferior a college | \(56.021 |[\)-197.776; \(309.818] |0,665 | |M1 |Educación: college |\)-213.428 | [$-446.182; \(19.327] |0,072 | |M1 |Educación: maestría |\)-214.361 | [$-438.256; \(9.534] |0,061 | |M1 |Educación: doctorado |\)-222.200 |
| M1 | Cargo: Dirección de investigación | $10.105.447 | [$9.530.888; $10.680.007] | <0,001 |
| M1 | Cargo: Dirección de manufactura | $3.065.515 | [$2.713.564; $3.417.467] | <0,001 |
| M1 | Cargo: Ejecutivo de ventas | $3.021.122 | [$2.773.884; $3.268.360] | <0,001 |
| M1 | Cargo: Gerencia | $10.584.359 | [$10.060.249; $11.108.469] | <0,001 |
| M1 | Cargo: Recursos humanos | $919.856 | [$352.579; $1.487.133] | 0,001 |
| M1 | Cargo: Representante de salud | $3.037.349 | [$2.668.768; \(3.405.930] |<0,001 | |M1 |Cargo: Representante de ventas |\)-13.254 | [$-266.208; $239.700] |
| M1 | Cargo: Técnico de laboratorio | \(20.246 |[\)-172.515; $213.007] | 0,837 | |
| M1 | Horas extra: sí frente a no | \(95.090 |[\)-94.623; $284.802] | 0,326 | |
| M2 | Constante | $1.797.325 | [$1.581.264; $2.013.386] | <0,001 |
| M2 | Experiencia (un año) | $189.988 | [$168.491; $211.485] | <0,001 |
| M2 | Antigüedad en la empresa (un año) | \(19.998 |[\)-2.836; $42.832] | 0,086 | |
| M2 | Educación: inferior a college | \(45.874 |[\)-208.553; \(300.302] |0,724 | |M2 |Educación: college |\)-213.080 | [$-445.347; \(19.187] |0,072 | |M2 |Educación: maestría |\)-219.454 | [$-443.562; \(4.654] |0,055 | |M2 |Educación: doctorado |\)-215.294 |
| M2 | Cargo: Dirección de investigación | $10.105.129 | [$9.532.910; $10.677.348] | <0,001 |
| M2 | Cargo: Dirección de manufactura | $3.082.366 | [$2.729.862; $3.434.869] | <0,001 |
| M2 | Cargo: Ejecutivo de ventas | $3.022.639 | [$2.775.637; $3.269.641] | <0,001 |
| M2 | Cargo: Gerencia | $10.591.735 | [$10.068.515; $11.114.955] | <0,001 |
| M2 | Cargo: Recursos humanos | $905.304 | [$333.440; $1.477.167] | 0,002 |
| M2 | Cargo: Representante de salud | $3.036.033 | [$2.667.173; $3.404.894] | <0,001 |
| M2 | Cargo: Representante de ventas | \(391 |[\)-256.797; $257.580] | 0,998 | |
| M2 | Cargo: Técnico de laboratorio | \(10.241 |[\)-182.521; $203.003] | 0,917 | |
| M2 | Horas extra: sí frente a no | \(102.730 |[\)-87.020; \(292.480] |0,288 | |M2 |Mujer frente a hombre |\)-170.605 | [$-347.821; $6.610] | 0,059 |
Los coeficientes de trayectoria son diferencias condicionadas al resto de predictores. Cada cargo se compara con investigación científica; cada nivel educativo, con pregrado. La constante combina las referencias con cero años de experiencia y antigüedad: su utilidad es algebraica y no constituye recomendación salarial.
En M2, la asociación por año de experiencia es $189.988 y por año de antigüedad es $19.998. Las diferencias entre cargos apoyan su inclusión, aunque no separan jerarquía ni responsabilidades dentro de ellos. Las tablas permiten evaluar cada contraste; un resultado no significativo no demuestra igualdad entre categorías.
El R² ajustado cercano a 0,876 describe ajuste interno. La significancia conjunta indica asociaciones en el conjunto de predictores; no garantiza significancia individual ni utilidad para una decisión salarial.
6.2 Recorridos de selección
Se aplica stats::step en ambas direcciones, desde M1 y
M2, con términos contenidos en cada especificación. Cada factor se
conserva o retira como término conjunto. La
documentación de R exige comparar sobre las mismas
observaciones.
tab(analysis$aic_paths)| Modelo | step_number | Step | Df | Deviance | Resid..Df | Resid..Dev | AIC |
|---|---|---|---|---|---|---|---|
| A1 | 0 | NA | NA | 1454 | 3.983224e+15 | 42114.94 | |
| A1 | 1 | - overtime | 1 | 2.673089e+12 | 1455 | 3.985897e+15 | 42113.93 |
| A1 | 2 | - education_category | 4 | 1.879026e+13 | 1459 | 4.004687e+15 | 42112.84 |
| A2 | 0 | NA | NA | 1453 | 3.973111e+15 | 42113.20 | |
| A2 | 1 | - education_category | 4 | 1.895381e+13 | 1457 | 3.992065e+15 | 42112.20 |
| A2 | 2 | - overtime | 1 | 2.775170e+12 | 1458 | 3.994840e+15 | 42111.22 |
tab(data.frame(model=c("A1", "A2"), formula=vapply(analysis$fits[c("A1", "A2")],
function(fit) paste(deparse(formula(fit)), collapse=" "), character(1))))| Modelo | Fórmula |
|---|---|
| A1 | monthly_income ~ years_of_experience + company_tenure + job_role |
| A2 | monthly_income ~ years_of_experience + company_tenure + job_role + gender |
En el reajuste completo, A1 retira educación y horas extra. A2 retira
esos términos y conserva género. Se reduce complejidad con cambios
pequeños de ajuste. El recorrido muestra el AIC interno de
step, distinto por una constante del AIC de la tabla
anterior; las diferencias se comparan dentro de cada tabla y escala.
La selección se repite dentro de cada entrenamiento. La fórmula elegida en toda la muestra no se reutiliza retrospectivamente en los pliegues. Los valores p posteriores a selección son exploratorios; la inferencia principal de género conserva M2.
7 Diagnósticos y transformaciones
z <- analysis$diagnostics
tab(data.frame(model=z$model, breusch_pagan_p=p_text(z$breusch_pagan_p),
shapiro_p=p_text(z$shapiro_p), robust_reset_p=p_text(z$robust_reset_p),
cooks_count=z$cooks_above_4_n, leverage_count=z$leverage_above_2_p_n,
studentized_count=z$studentized_above_3, check.names=FALSE))| Modelo | p Breusch-Pagan | p Shapiro-Wilk | p RESET HC3 | Cook > 4/n | Apalancamiento > 2p/n | Residuo studentizado > 3 |
|---|---|---|---|---|---|---|
| S1 | <0,001 | <0,001 | <0,001 | 139 | 112 | 7 |
| M1 | <0,001 | <0,001 | <0,001 | 113 | 94 | 4 |
| M2 | <0,001 | <0,001 | <0,001 | 113 | 83 | 5 |
| L2 | <0,001 | 0,299 | <0,001 | 72 | 83 | 5 |
| Q2 | <0,001 | <0,001 | <0,001 | 110 | 86 | 6 |
| A1 | <0,001 | <0,001 | <0,001 | 121 | 139 | 4 |
Residuos frente al ajuste y cuantiles normales.
Residuos frente al ajuste y cuantiles normales.
Dispersión e influencia: los umbrales identifican registros para revisión.
Dispersión e influencia: los umbrales identifican registros para revisión.
- Linealidad y especificación: RESET robusta detecta estructura no recogida. Q2 añade experiencia al cuadrado, pero no resuelve el problema ni mejora MAE. La recta interpretable mantiene sesgos locales.
- Varianza constante: Breusch-Pagan rechaza homocedasticidad en los modelos revisados. Se utiliza covarianza HC3 para coeficientes. Esto no corrige una media mal especificada ni valida intervalos clásicos individuales.
- Normalidad: pruebas y gráficos muestran desviaciones en la escala original. L2 mejora la distribución residual, pero conserva problemas de varianza y especificación. Los gráficos y magnitudes complementan los valores p.
- Independencia: no hay estructura temporal, identificadores de equipos ni diseño de levantamiento. Las filas no forman una serie temporal; Durbin-Watson sobre ese orden no demostraría independencia. Se asume intercambiabilidad para la validación y calibración, con el alcance limitado de datos sintéticos.
- Influencia: M2 señala 113 filas por Cook > 4/n. El umbral no acredita errores. Se conserva la muestra principal y se calcula una sensibilidad de género retirándolas conjuntamente.
tab(data.frame(term=analysis$vif$term, gvif=num(analysis$vif$gvif, 3), degrees=analysis$vif$df,
adjusted_gvif=num(analysis$vif$adjusted_gvif, 3), check.names=FALSE))| Término | GVIF | Grados | GVIF ajustado |
|---|---|---|---|
| years_of_experience | 2,423 | 1 | 1,557 |
| company_tenure | 1,678 | 1 | 1,295 |
| education_category | 1,059 | 4 | 1,007 |
| job_role | 1,865 | 8 | 1,040 |
| overtime | 1,011 | 1 | 1,005 |
| gender | 1,015 | 1 | 1,008 |
El GVIF ajustado es GVIF^(1/(2*df)); no se compara
directamente con umbrales del VIF ordinario. Los valores no indican
inflación extrema en M2, aunque experiencia y antigüedad comparten
información.
L1 y L2 usan logaritmo natural del ingreso. exp(beta)-1
expresa una diferencia proporcional en la escala geométrica
condicionada. En Q2 la pendiente de experiencia es
beta_experience + 2*beta_squared*experience. Las
predicciones logarítmicas regresan a pesos mediante el factor de Duan
estimado en entrenamiento. No se comparan directamente sus R² o AIC con
los de ingreso sin transformar.
Se mantiene la escala original porque los logaritmos aumentan MAE en pesos y la alternativa cuadrática no mejora el desempeño. Persisten limitaciones de especificación; la elección no supone que todos los supuestos se satisfacen.
8 Validación y elección predictiva
z <- analysis$cv_summary
tab(data.frame(model=z$model, mae_cop=money(z$mae), repeat_sd=money(z$repeat_sd),
rmse_cop=money(z$rmse), mean_parameters=num(z$mean_parameters),
improvement_median=pct(z$improvement_over_median_pct), check.names=FALSE))| Modelo | MAE COP | DE entre repeticiones | RMSE COP | Parámetros medios | Mejora frente a mediana |
|---|---|---|---|---|---|
| M2 | $1.296.756 | $2.788 | $1.655.279 | 17,00 | 61,6% |
| A1 | $1.296.849 | $1.875 | $1.651.939 | 11,00 | 61,6% |
| M1 | $1.297.376 | $2.656 | $1.655.753 | 16,00 | 61,6% |
| Q2 | $1.297.933 | $2.486 | $1.657.291 | 18,00 | 61,6% |
| A2 | $1.298.374 | $1.419 | $1.653.692 | 11,74 | 61,6% |
| L1 | $1.370.117 | $2.535 | $1.930.689 | 16,00 | 59,4% |
| L2 | $1.371.164 | $2.682 | $1.934.757 | 17,00 | 59,4% |
| S1 | $2.212.231 | $374 | $2.975.180 | 2,00 | 34,5% |
| S2 | $2.955.735 | $1.078 | $4.026.330 | 2,00 | 12,5% |
| median | $3.378.131 | $1.644 | $5.007.471 | 1,00 | 0,0% |
| S3 | $3.663.512 | $1.141 | $4.759.520 | 2,00 | -8,4% |
MAE en registros y particiones comunes de desarrollo.
El menor MAE es de M2. A1 difiere del mínimo en apenas $92. Se utiliza una regla explícita de parsimonia: entre candidatos a menos de 0,5% del menor MAE, elegir el de menor número medio de parámetros en entrenamiento y desempatar por MAE. Esta tolerancia pragmática no es una prueba de equivalencia ni una política institucional. Se selecciona A1.
La desviación entre repeticiones describe estabilidad frente a particiones; no es un error estándar independiente, porque se reutilizan registros. La similitud de los mejores modelos tampoco acredita igualdad fuera de esta base.
Añadir género a M1 reduce MAE en $620. Ese aporte pequeño no resuelve si existe una diferencia ajustada relevante; esa pregunta requiere coeficientes e incertidumbre.
En las 294 filas reservadas, el procedimiento elegido tiene MAE de $1.373.696, frente a $3.241.976 de la mediana de desarrollo. No se utiliza ese resultado para volver a elegir candidato y no representa validación en otra empresa.
9 Diferencias ajustadas por género
z <- g[g$model != "L2", ]
tab(data.frame(specification=z$model, female_minus_male=money(z$estimate),
hc3_ci=paste0("[", money(z$ci_lower), "; ", money(z$ci_upper), "]"),
hc3_p=p_text(z$p_value), check.names=FALSE))| Especificación | Mujer menos hombre | IC 95% HC3 | p HC3 |
|---|---|---|---|
| S3 | \(306.058 |[\)-185.490; \(797.607] |0,222 | |M2 |\)-170.605 | [$-347.821; \(6.610] |0,059 | |Q2 |\)-169.455 | [$-346.785; \(7.875] |0,061 | |career_adjustment |\)-44.449 |
Diferencias estimadas y su incertidumbre según los controles.
La diferencia inicial es positiva, se aproxima a cero al ajustar por trayectoria y educación y se vuelve negativa al añadir cargo y horas extra. El cambio muestra que las especificaciones responden preguntas distintas y dependen de la composición de los grupos. No demuestra que cargo explique causalmente una desigualdad ni examina acceso a los cargos.
En M2, el intervalo es [$-347.821; $6.610], con p = 0,059. Incluye cero y diferencias negativas potencialmente relevantes. No rechazar al 5% no demuestra igualdad. En L2, la diferencia geométrica es -1,0%, con intervalo [-3,9%; 1,9%].
La sensibilidad que excluye las filas señaladas por Cook obtiene una diferencia significativa. Es una advertencia de fragilidad: no se acreditaron errores en esas filas y retirarlas modifica la muestra. No se elige una especificación por producir significancia.
Faltan nivel jerárquico, jornada contratada, cantidad de horas extra y responsabilidades dentro del cargo. La Office for National Statistics advierte que el componente no explicado no mide automáticamente discriminación. Aquí se añade la limitación del origen sintético.
10 Escenarios e intervalos
10.1 Método de calibración
Se utiliza conformal dividido. El modelo se ajusta en desarrollo y
allí se estima una escala residual por cargo. En calibración se calcula
abs(actual_income - prediction) / role_scale. El
estadístico de orden 281 entre 294 errores normalizados determina el
intervalo prediction ± quantile * role_scale. El límite
inferior se trunca en cero por el soporte positivo del ingreso.
La amplitud cambia por cargo sin exigir normalidad u homocedasticidad. El 95% nominal es cobertura marginal bajo intercambiabilidad, no garantía para cada género, cargo o persona. El método se fundamenta en Lei et al. (2018). La proporción cubierta en calibración no se presenta como evaluación independiente de cobertura futura.
Los coeficientes descriptivos se reajustaron con todos los registros. Los escenarios conservan el ajuste de desarrollo, sin sustituirlo por el ajuste completo después de calibrar. Los objetos guardados identifican ambos usos.
10.2 Perfiles y resultados
Los perfiles reproducen combinaciones de predictores de desarrollo, cercanas a la trayectoria central dentro de dos cargos, elegidas sin consultar sus ingresos. Son escenarios hipotéticos comparables, no recomendaciones para empleados identificables.
z <- analysis$scenarios
tab(data.frame(profile=c("Laboratorio", "Gerencia"), job_role=unname(role_labels[as.character(z$job_role)]),
experience=z$years_of_experience, company_tenure=z$company_tenure,
education=unname(education_labels[as.character(z$education_category)]),
overtime=as.character(z$overtime), gender=as.character(z$gender), check.names=FALSE))| Perfil | Cargo | Experiencia | Antigüedad | Educación | Horas extra | Género |
|---|---|---|---|---|---|---|
| Laboratorio | Técnico de laboratorio | 7 | 3 | Pregrado | No | M |
| Gerencia | Gerencia | 23 | 13 | Pregrado | No | F |
tab(data.frame(profile=c("Laboratorio", "Gerencia"), prediction=money(z$prediction),
prediction_lower=money(z$lower), prediction_upper=money(z$upper), check.names=FALSE))| Perfil | Predicción | Límite inferior 95% | Límite superior 95% |
|---|---|---|---|
| Laboratorio | $3.108.597 | $670.002 | $5.547.192 |
| Gerencia | $16.894.114 | $12.858.252 | $20.929.976 |
Predicciones e intervalos conformales individuales.
En laboratorio, el intervalo abarca referencias salariales muy diferentes. En gerencia, el ingreso y la incertidumbre absoluta son mayores. Para un presupuesto ilustrativo puede ser útil comparar magnitudes; comprometer un salario requiere información adicional y una base real validada. Las diferencias entre estos perfiles reúnen cambios de cargo y trayectoria y no estiman un efecto de género.
11 Análisis crítico y recomendación
Usuario y decisión. Un equipo hipotético de compensación podría utilizar el ejercicio para aprender a revisar referencias salariales y comunicar incertidumbre. El uso operativo requiere validación con registros reales de la organización.
Error y aceptabilidad. El MAE de $1.296.849 es la desviación absoluta promedio en validación interna. Como comparación de magnitudes, equivale a 26,4% de la mediana de ingreso; esa razón no es MAPE ni un error porcentual individual. No hay tolerancia institucional aportada. Si se exigiera, como ejemplo, precisión individual de ±10%, no habría evidencia para prometerla. Los intervalos desaconsejan tratar el punto como salario exacto.
Consecuencias de equivocarse. Sobreestimar puede elevar compromisos presupuestales o distorsionar comparaciones internas. Subestimar puede producir ofertas insuficientes o señalar como excesivo un ingreso justificado por responsabilidades omitidas. MAE no valora esos costos diferentes.
Limitaciones y usos inadecuados. El origen sintético, la escala adaptada, la especificación incompleta y las variables omitidas impiden recomendarlo para salarios reales, evaluar discriminación o extrapolar al mercado colombiano. Tampoco cubre cargos ausentes ni trayectorias fuera de soporte. La validación aleatoria no evalúa cambios económicos futuros.
Recomendación e información faltante. Se recomienda como ejercicio reproducible y discusión metodológica. Para uso salarial se requieren registros reales y autorizados, período y moneda verificables, nivel de responsabilidad, jornada, componentes de remuneración y una tolerancia acordada. Después, validar desempeño y cobertura en datos externos o posteriores. La pregunta de género necesita un diseño acorde con la interpretación buscada y no se resuelve optimizando MAE.
12 Referencias y reproducibilidad
- Guía del Caso 02 - Regresión lineal simple y múltiple. Material del curso, tres páginas. Referencia local del alcance y los entregables.
- IBM. Archivo público de referencia. Verificación de correspondencia y escala, sin incorporar predictores externos.
- Posit Academy. IBM HR Analytics. Documentación de la naturaleza sintética.
- R Core Team. Selección por AIC y predicción lineal.
- Office for National Statistics (2018). Understanding the gender pay gap in the UK. Referencia metodológica; no se trasladan sus resultados británicos a esta base.
- Lei, J., G’Sell, M., Rinaldo, A., Tibshirani, R. J., y Wasserman, L. (2018). Distribution-Free Predictive Inference for Regression. DOI: 10.1080/01621459.2017.1307116.
Se reproduce abriendo el proyecto y ejecutando
source("render_report.R", encoding="UTF-8"). Requiere
readxl, knitr, rmarkdown,
ggplot2, lmtest, car y
jsonlite. La referencia IBM se conserva localmente para
evitar descargas en cada ejecución. La huella MD5 del archivo de entrada
es 52633bafd9a06fb09c9aa5441b7dff2c.
Las tablas, asignaciones de muestra y pliegues, predicciones de
validación, modelos y figuras quedan en analysis_outputs/.
analysis_manifest.json registra semilla, tamaños, selección
y versiones. Los datos originales permanecen intactos.
Código completo de preparación, estimación, validación y figuras
# Reproducible analysis functions. Monetary calculations retain source COP units.
income_formulas <- list(
S1 = monthly_income ~ years_of_experience,
S2 = monthly_income ~ company_tenure,
S3 = monthly_income ~ gender,
M1 = monthly_income ~ years_of_experience + company_tenure + education_category + job_role + overtime,
M2 = monthly_income ~ years_of_experience + company_tenure + education_category + job_role + overtime + gender,
L1 = log(monthly_income) ~ years_of_experience + company_tenure + education_category + job_role + overtime,
L2 = log(monthly_income) ~ years_of_experience + company_tenure + education_category + job_role + overtime + gender,
Q2 = monthly_income ~ years_of_experience + I(years_of_experience^2) + company_tenure + education_category + job_role + overtime + gender
)
fit_candidate <- function(id, data) {
parent <- if (id == "A1") "M1" else if (id == "A2") "M2" else id
model_formula <- income_formulas[[parent]]
environment(model_formula) <- environment()
fit <- stats::lm(model_formula, data = data, na.action = na.fail)
fit$call$formula <- model_formula
if (id %in% c("A1", "A2")) {
fit <- stats::step(fit, scope = list(lower = monthly_income ~ 1,
upper = income_formulas[[parent]]),
direction = "both", trace = 0)
}
stopifnot(fit$rank == ncol(model.matrix(fit)))
fit
}
predict_candidate <- function(fit, newdata) {
prediction <- as.numeric(predict(fit, newdata = newdata))
if (grepl("^log\\(", deparse(formula(fit)[[2]]))) {
prediction <- exp(prediction) * mean(exp(residuals(fit)))
}
stopifnot(all(is.finite(prediction)))
prediction
}
coefficient_table <- function(fit, id) {
covariance <- car::hccm(fit, type = "hc3")
estimates <- coef(fit)
standard_errors <- sqrt(diag(covariance))
statistic <- estimates / standard_errors
critical <- qt(0.975, df.residual(fit))
data.frame(model = id, term = names(estimates), estimate = unname(estimates),
hc3_se = unname(standard_errors), statistic = unname(statistic),
p_value = 2 * pt(-abs(statistic), df.residual(fit)),
ci_lower = unname(estimates - critical * standard_errors),
ci_upper = unname(estimates + critical * standard_errors), row.names = NULL)
}
model_summary <- function(fit, id) {
s <- summary(fit)
f <- s$fstatistic
covariance <- car::hccm(fit, type = "hc3")[-1, -1, drop = FALSE]
b <- coef(fit)[-1]
robust_f <- as.numeric(crossprod(b, solve(covariance, b)) / length(b))
data.frame(model = id, n = nobs(fit), parameters = length(coef(fit)),
response_scale = if (id %in% c("L1", "L2")) "log_income" else "income_cop",
r_squared = s$r.squared, adjusted_r_squared = s$adj.r.squared,
aic = AIC(fit), classical_f_p = pf(f[1], f[2], f[3], lower.tail = FALSE),
hc3_joint_p = pf(robust_f, length(b), df.residual(fit), lower.tail = FALSE),
formula = paste(deparse(formula(fit)), collapse = " "), row.names = NULL)
}
diagnose_model <- function(fit, id) {
bp <- lmtest::bptest(fit)
reset <- lmtest::resettest(fit, power = 2:3, type = "fitted", vcov = function(x) car::hccm(x, type = "hc3"))
sw <- shapiro.test(residuals(fit))
n <- nobs(fit)
p <- length(coef(fit))
data.frame(model = id, breusch_pagan_p = bp$p.value,
shapiro_p = sw$p.value, robust_reset_p = reset$p.value,
max_cooks_distance = max(cooks.distance(fit)),
cooks_above_4_n = sum(cooks.distance(fit) > 4/n),
leverage_above_2_p_n = sum(hatvalues(fit) > 2*p/n),
studentized_above_3 = sum(abs(rstudent(fit)) > 3), row.names = NULL)
}
save_csv <- function(data, name, directory) {
utils::write.csv(data, file.path(directory, paste0(name, ".csv")), row.names = FALSE,
fileEncoding = "UTF-8", na = "")
}
run_income_analysis <- function() {
options(contrasts = c("contr.treatment", "contr.poly"))
output <- "analysis_outputs"
dir.create(output, recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(output, "figures"), showWarnings = FALSE)
source_file <- "source_data/dataset_03_employee_income_and_job_characteristics.xlsx"
source_hash <- unname(tools::md5sum(source_file))
source_data <- as.data.frame(readxl::read_excel(source_file, sheet = "employee_income"))
definitions <- as.data.frame(readxl::read_excel(source_file, sheet = "Data Definitions", col_names = FALSE))
required <- c("monthly_income", "years_of_experience", "company_tenure", "education_level", "job_role", "overtime", "gender")
stopifnot(all(required %in% names(source_data)), nrow(source_data) == 1470L)
data <- source_data
data$record_id <- seq_len(nrow(data))
data$education_category <- factor(data$education_level, levels = c(3, 1, 2, 4, 5),
labels = c("bachelor", "below_college", "college", "master", "doctorate"))
data$gender <- factor(data$gender, levels = c("M", "F"))
data$overtime <- factor(trimws(data$overtime), levels = c("No", "Si"))
data$job_role <- relevel(factor(data$job_role), ref = "Investigador_Cientifico")
stopifnot(!anyNA(data[, c(required, "education_category")]), all(data$monthly_income > 0))
audit <- data.frame(check = c("source_rows", "source_columns", "missing_cells", "exact_repeated_rows",
"nonpositive_income", "negative_experience", "tenure_above_experience",
"role_tenure_above_company_tenure", "experience_above_age_minus_14"),
value = c(nrow(data), ncol(source_data), sum(is.na(source_data)),
sum(duplicated(source_data)), sum(data$monthly_income <= 0),
sum(data$years_of_experience < 0),
sum(data$company_tenure > data$years_of_experience),
sum(data$role_tenure > data$company_tenure),
sum(data$years_of_experience > data$age - 14)))
q <- quantile(data$monthly_income, c(.25, .75))
outlier_ids <- data$record_id[data$monthly_income < q[1] - 1.5*diff(q) | data$monthly_income > q[2] + 1.5*diff(q)]
# Flags describe observations; no records are removed or winsorized.
audit <- rbind(audit, data.frame(check = "income_tukey_flags_retained", value = length(outlier_ids)))
source_columns <- data.frame(column = names(source_data),
data_type = vapply(source_data, function(x) class(x)[1], character(1)),
missing = colSums(is.na(source_data)),
distinct = vapply(source_data, function(x) length(unique(x)), integer(1)))
numeric_names <- c("monthly_income", "years_of_experience", "company_tenure", "age")
descriptive <- do.call(rbind, lapply(numeric_names, function(name) {
x <- data[[name]]
data.frame(variable = name, n = length(x), mean = mean(x), median = median(x), sd = sd(x),
minimum = min(x), q1 = unname(quantile(x, .25)), q3 = unname(quantile(x, .75)), maximum = max(x))
}))
grouped_summary <- function(variable) do.call(rbind, lapply(split(data, data[[variable]]), function(z) {
data.frame(group_variable = variable, group = as.character(z[[variable]][1]), n = nrow(z),
mean_income = mean(z$monthly_income), median_income = median(z$monthly_income), sd_income = sd(z$monthly_income))
}))
group_summaries <- do.call(rbind, lapply(c("gender", "education_category", "job_role", "overtime"), grouped_summary))
correlations <- do.call(rbind, lapply(c("years_of_experience", "company_tenure", "age"), function(name) {
data.frame(predictor = name, pearson = cor(data$monthly_income, data[[name]]),
spearman = cor(data$monthly_income, data[[name]], method = "spearman"))
}))
# Verify source provenance without importing reference-only predictors into models.
reference_file <- file.path(output, "reference_data", "ibm_employee_reference.csv")
if (!file.exists(reference_file)) {
stop("Restore analysis_outputs/reference_data/ibm_employee_reference.csv before running provenance verification.")
}
provenance <- data.frame(check = "reference_file_available", matches = file.exists(reference_file))
if (file.exists(reference_file)) {
reference <- read.csv(reference_file, fileEncoding = "UTF-8-BOM")
reference_pairs <- c(age = "Age", education_level = "Education", years_of_experience = "TotalWorkingYears",
company_tenure = "YearsAtCompany", distance_from_home = "DistanceFromHome",
environment_satisfaction = "EnvironmentSatisfaction", job_satisfaction = "JobSatisfaction",
previous_jobs = "NumCompaniesWorked", salary_increase_pct = "PercentSalaryHike",
performance_rating = "PerformanceRating", training_sessions = "TrainingTimesLastYear",
work_life_balance = "WorkLifeBalance", role_tenure = "YearsInCurrentRole",
years_since_last_promotion = "YearsSinceLastPromotion", years_with_current_manager = "YearsWithCurrManager")
stopifnot(nrow(reference) == nrow(data))
provenance <- rbind(provenance, data.frame(check = names(reference_pairs), matches = vapply(names(reference_pairs),
function(name) all(source_data[[name]] == reference[[reference_pairs[[name]]]]), logical(1))),
data.frame(check = c("income_is_reference_times_1000", "gender_translation"),
matches = c(all(data$monthly_income == reference$MonthlyIncome * 1000),
all(as.character(data$gender) == ifelse(reference$Gender == "Female", "F", "M")))))
stopifnot(all(provenance$matches))
}
# Random calibration reservation precedes fitting and predictive selection.
set.seed(2026)
calibration_ids <- sort(sample(seq_len(nrow(data)), size = ceiling(.20*nrow(data))))
development <- data[-calibration_ids, ]
calibration <- data[calibration_ids, ]
partitions <- data.frame(record_id = data$record_id,
partition = ifelse(data$record_id %in% calibration_ids, "calibration", "development"))
set.seed(2026)
folds <- do.call(rbind, lapply(seq_len(5), function(repetition) {
fold_number <- integer(nrow(development))
for (group in levels(development$gender)) {
positions <- which(development$gender == group)
fold_number[positions] <- sample(rep(seq_len(10), length.out = length(positions)))
}
data.frame(record_id = development$record_id, repetition = repetition, fold = fold_number)
}))
ids <- c("S1", "S2", "S3", "M1", "M2", "A1", "A2", "L1", "L2", "Q2")
cv_predictions <- list()
selection_history <- list()
index <- 0L
for (repetition in seq_len(5)) {
assignment <- folds[folds$repetition == repetition, ]
for (fold in seq_len(10)) {
test_ids <- assignment$record_id[assignment$fold == fold]
train <- development[!development$record_id %in% test_ids, ]
test <- development[development$record_id %in% test_ids, ]
for (variable in c("gender", "overtime", "education_category", "job_role")) {
stopifnot(all(unique(as.character(test[[variable]])) %in% unique(as.character(train[[variable]]))))
}
for (id in c("median", ids)) {
if (id == "median") {
prediction <- rep(median(train$monthly_income), nrow(test))
} else {
fit <- fit_candidate(id, train)
prediction <- predict_candidate(fit, test)
if (id %in% c("A1", "A2")) {
selection_history[[length(selection_history)+1L]] <- data.frame(
model = id, repetition = repetition, fold = fold,
formula = paste(deparse(formula(fit)), collapse = " "))
}
}
index <- index + 1L
cv_predictions[[index]] <- data.frame(model = id, repetition = repetition, fold = fold,
record_id = test$record_id, gender = as.character(test$gender),
actual = test$monthly_income, prediction = prediction, error = test$monthly_income - prediction,
parameters = if (id == "median") 1L else length(coef(fit)))
}
}
message("Completed cross-validation repetition ", repetition, "/5")
}
cv_predictions <- do.call(rbind, cv_predictions)
selection_history <- do.call(rbind, selection_history)
repeat_metrics <- do.call(rbind, lapply(split(cv_predictions, interaction(cv_predictions$model, cv_predictions$repetition)), function(z) {
data.frame(model = z$model[1], repetition = z$repetition[1], n = nrow(z),
mae = mean(abs(z$error)), rmse = sqrt(mean(z$error^2)), bias_actual_minus_predicted = mean(z$error))
}))
cv_summary <- do.call(rbind, lapply(split(repeat_metrics, repeat_metrics$model), function(z) {
data.frame(model = z$model[1], mae = mean(z$mae), repeat_sd = sd(z$mae),
min_repeat_mae = min(z$mae), max_repeat_mae = max(z$mae), rmse = mean(z$rmse),
unique_records = nrow(development), predictions = sum(z$n),
mean_parameters = mean(cv_predictions$parameters[cv_predictions$model == z$model[1]]))
}))
cv_summary <- cv_summary[order(cv_summary$mae), ]
baseline_mae <- cv_summary$mae[cv_summary$model == "median"]
cv_summary$improvement_over_median_pct <- 100*(1-cv_summary$mae/baseline_mae)
# A transparent parsimony rule treats MAE differences below 0.5% as negligible
# for this academic comparison, without claiming statistical equivalence.
minimum_mae <- min(cv_summary$mae[cv_summary$model != "median"])
near_best <- cv_summary[cv_summary$model != "median" & cv_summary$mae <= 1.005*minimum_mae, ]
near_best <- near_best[order(near_best$mean_parameters, near_best$mae), ]
selected_id <- near_best$model[1]
best_simple <- cv_summary$model[cv_summary$model %in% c("S1", "S2", "S3")][1]
predictive_fit <- fit_candidate(selected_id, development)
# Full-sample refits serve coefficient interpretation and diagnostic sensitivity.
fits <- setNames(lapply(ids, function(id) fit_candidate(id, data)), ids)
coefficients <- do.call(rbind, lapply(ids, function(id) coefficient_table(fits[[id]], id)))
fit_summary <- do.call(rbind, lapply(ids, function(id) model_summary(fits[[id]], id)))
diagnostics <- do.call(rbind, lapply(c("S1", "M1", "M2", "L2", "Q2", selected_id) |> unique(),
function(id) diagnose_model(fits[[id]], id)))
aic_paths <- do.call(rbind, lapply(c("A1", "A2"), function(id) {
x <- as.data.frame(fits[[id]]$anova)
names(x) <- make.names(names(x))
cbind(model = id, step_number = seq_len(nrow(x))-1L, x)
}))
vif_result <- car::vif(fits[["M2"]])
vif_table <- data.frame(term = rownames(vif_result), gvif = vif_result[,1],
df = vif_result[,2], adjusted_gvif = vif_result[,3], row.names = NULL)
influence <- data.frame(record_id = data$record_id, income = data$monthly_income,
cooks_distance = cooks.distance(fits[["M2"]]),
leverage = hatvalues(fits[["M2"]]), studentized_residual = rstudent(fits[["M2"]]))
influence$flag <- influence$cooks_distance > 4/nrow(data)
sensitivity_fit <- lm(income_formulas$M2, data = data[!influence$flag, ])
sensitivity_gender <- coefficient_table(sensitivity_fit, "M2_influence_sensitivity")
sensitivity_gender <- sensitivity_gender[sensitivity_gender$term == "genderF", ]
career_fit <- lm(monthly_income ~ gender + years_of_experience + company_tenure + education_category, data = data)
gender_results <- rbind(coefficients[coefficients$term == "genderF" & coefficients$model %in% c("S3", "M2", "L2", "Q2"), ],
coefficient_table(career_fit, "career_adjustment")["genderF" == names(coef(career_fit)), ],
sensitivity_gender)
# The calibration response is used only after the candidate and fitting rule are locked.
training_residual <- development$monthly_income - predict_candidate(predictive_fit, development)
role_scale <- tapply(training_residual, development$job_role, function(x) sqrt(mean(x^2)))
role_scale <- pmax(role_scale, .10*sqrt(mean(training_residual^2)))
names(role_scale) <- levels(development$job_role)
calibration_prediction <- predict_candidate(predictive_fit, calibration)
calibration_scores <- abs(calibration$monthly_income-calibration_prediction) / role_scale[as.character(calibration$job_role)]
conformal_rank <- ceiling((nrow(calibration)+1)*.95)
conformal_quantile <- unname(sort(calibration_scores)[conformal_rank])
calibration_results <- data.frame(record_id = calibration$record_id, gender = as.character(calibration$gender),
job_role = as.character(calibration$job_role), actual = calibration$monthly_income,
prediction = calibration_prediction, scale = unname(role_scale[as.character(calibration$job_role)]),
score = unname(calibration_scores))
calibration_results$lower <- pmax(0, calibration_results$prediction - conformal_quantile*calibration_results$scale)
calibration_results$upper <- calibration_results$prediction + conformal_quantile*calibration_results$scale
# Select representative observed predictor profiles; outcomes are not used to choose them.
scenarios <- do.call(rbind, lapply(c("Tecnico_Laboratorio", "Gerente"), function(role) {
pool <- development[as.character(development$job_role) == role, ]
mode_education <- names(which.max(table(pool$education_category)))
pool <- pool[as.character(pool$education_category) == mode_education, ]
distance <- abs(pool$years_of_experience-median(pool$years_of_experience)) +
abs(pool$company_tenure-median(pool$company_tenure))
pool[which.min(distance), ]
}))
scenarios$scenario <- c("laboratory_profile", "management_profile")
scenarios$prediction <- predict_candidate(predictive_fit, scenarios)
scenarios$role_scale <- unname(role_scale[as.character(scenarios$job_role)])
scenarios$lower <- pmax(0, scenarios$prediction-conformal_quantile*scenarios$role_scale)
scenarios$upper <- scenarios$prediction+conformal_quantile*scenarios$role_scale
scenarios <- scenarios[, c("scenario", "record_id", "years_of_experience", "company_tenure", "education_category",
"job_role", "overtime", "gender", "prediction", "lower", "upper")]
calibration_metrics <- data.frame(n = nrow(calibration), mae = mean(abs(calibration$monthly_income-calibration_prediction)),
median_reference_mae = mean(abs(calibration$monthly_income-median(development$monthly_income))),
conformal_rank = conformal_rank, conformal_quantile = conformal_quantile,
calibrated_fraction = mean(calibration_results$actual >= calibration_results$lower & calibration_results$actual <= calibration_results$upper))
group_cv <- do.call(rbind, lapply(split(cv_predictions, interaction(cv_predictions$model, cv_predictions$gender)), function(z) {
data.frame(model = z$model[1], gender = z$gender[1], mae = mean(abs(z$error)), mean_error = mean(z$error))
}))
# Verify the adjusted coefficient with an independent partial-regression calculation.
control_matrix <- model.matrix(fits$M1)
gender_indicator <- as.numeric(data$gender == "F")
partial_gender <- lm.fit(control_matrix, gender_indicator)$residuals
partial_income <- lm.fit(control_matrix, data$monthly_income)$residuals
partial_coefficient <- sum(partial_gender*partial_income)/sum(partial_gender^2)
stopifnot(abs(partial_coefficient-coef(fits$M2)[["genderF"]]) < 1e-5)
stopifnot(unname(tools::md5sum(source_file)) == source_hash,
all(is.finite(cv_predictions$prediction)),
all(table(cv_predictions$model, cv_predictions$record_id) == 5),
length(intersect(calibration$record_id, development$record_id)) == 0,
all(scenarios$lower <= scenarios$prediction), all(scenarios$prediction <= scenarios$upper))
tables <- list(data_quality = audit, source_columns = source_columns, numeric_summary = descriptive,
group_summary = group_summaries, income_correlations = correlations, provenance_checks = provenance,
sample_partitions = partitions, cv_folds = folds, cv_predictions = cv_predictions,
cv_repeat_metrics = repeat_metrics, cv_model_comparison = cv_summary, cv_aic_selections = selection_history,
model_coefficients_hc3 = coefficients, model_fit_summary = fit_summary, aic_paths = aic_paths,
diagnostics = diagnostics, multicollinearity = vif_table, influence = influence,
gender_comparison = gender_results, group_prediction_errors = group_cv,
calibration_predictions = calibration_results, calibration_summary = calibration_metrics,
prediction_scenarios = scenarios)
for (name in names(tables)) save_csv(tables[[name]], name, output)
save_csv(data, "analysis_data", output)
role_composition <- as.data.frame.matrix(table(data$job_role, data$gender))
role_composition$job_role <- rownames(role_composition)
save_csv(role_composition, "role_gender_composition", output)
model_bundle <- list(selected_id = selected_id, predictive_model = predictive_fit,
development_ids = development$record_id, calibration_ids = calibration$record_id,
role_scale = role_scale, conformal_quantile = conformal_quantile,
full_sample_models = fits, source_md5 = source_hash)
saveRDS(model_bundle, file.path(output, "income_models.rds"))
results <- list(data = data, source_data = source_data, definitions = definitions, audit = audit,
descriptive = descriptive, group_summaries = group_summaries, correlations = correlations,
provenance = provenance, development = development, calibration = calibration,
cv_summary = cv_summary, repeat_metrics = repeat_metrics, cv_predictions = cv_predictions,
fits = fits, coefficients = coefficients, fit_summary = fit_summary,
diagnostics = diagnostics, aic_paths = aic_paths, vif = vif_table, influence = influence,
gender_results = gender_results, group_cv = group_cv, scenarios = scenarios,
calibration_metrics = calibration_metrics, calibration_results = calibration_results,
selected_id = selected_id, best_simple = best_simple, source_md5 = source_hash,
role_composition = role_composition)
saveRDS(results, file.path(output, "analysis_results.rds"))
manifest <- list(source_file = source_file, source_md5 = source_hash, n_total = nrow(data),
n_development = nrow(development), n_calibration = nrow(calibration), seed = 2026,
folds = 10, repetitions = 5, candidate_models = c("median", ids), selected_id = selected_id,
best_simple = best_simple, conformal_rank = conformal_rank, r_version = R.version.string,
package_versions = setNames(lapply(c("readxl", "knitr", "rmarkdown", "ggplot2", "lmtest", "car", "jsonlite"),
function(x) { requireNamespace(x, quietly = TRUE); as.character(packageVersion(x)) }),
c("readxl", "knitr", "rmarkdown", "ggplot2", "lmtest", "car", "jsonlite")))
jsonlite::write_json(manifest, file.path(output, "analysis_manifest.json"), auto_unbox = TRUE, pretty = TRUE)
results
}
create_income_figures <- function(result) {
library(ggplot2)
directory <- "analysis_outputs/figures"
theme_set(theme_minimal(base_size = 12, base_family = "sans") +
theme(plot.title = element_text(face = "bold"), panel.grid.minor = element_blank(),
legend.position = "bottom", plot.margin = margin(12, 18, 12, 12)))
save_plot <- function(name, plot, height = 5.2, width = 8.8) {
ggsave(file.path(directory, paste0(name, ".png")), plot, width = width, height = height,
dpi = 160, bg = "white")
}
d <- result$data
save_plot("income_distribution", ggplot(d, aes(monthly_income/1e6)) +
geom_histogram(binwidth = 1, fill = "#315B7D", color = "white") +
geom_vline(xintercept = median(d$monthly_income)/1e6, linetype = "dashed", color = "#B14B35") +
labs(title = "Distribución del ingreso mensual", subtitle = "La línea señala la mediana; valores en la escala académica de la base",
x = "Ingreso mensual (millones de pesos)", y = "Empleados"))
scatter <- rbind(data.frame(income = d$monthly_income/1e6, x = d$years_of_experience, predictor = "Experiencia total"),
data.frame(income = d$monthly_income/1e6, x = d$company_tenure, predictor = "Antigüedad en la empresa"))
save_plot("income_career_relationships", ggplot(scatter, aes(x, income)) +
geom_point(alpha = .20, size = 1, color = "#315B7D") + geom_smooth(method = "lm", formula = y~x, se = FALSE, color = "#B14B35") +
geom_smooth(method = "loess", formula = y~x, se = FALSE, color = "#467D62", linetype = "dashed") +
facet_wrap(~predictor) + labs(title = "Ingreso y trayectoria laboral", subtitle = "Recta OLS en rojo; tendencia local en verde",
x = "Años", y = "Ingreso mensual (millones de pesos)"))
save_plot("income_gender", ggplot(d, aes(gender, monthly_income/1e6, fill = gender)) +
geom_boxplot(width = .5, outlier.alpha = .2) + scale_fill_manual(values = c("M"="#A8BDC9", "F"="#D9AA8A")) +
labs(title = "Distribución del ingreso por género", x = "Género registrado", y = "Ingreso mensual (millones de pesos)") + guides(fill = "none"))
role_labels <- c("Investigador_Cientifico"="Investigación científica", "Director_Investigación"="Dirección de investigación",
"Director_Manofactura"="Dirección de manufactura", "Ejecutivo_Ventas"="Ejecutivo de ventas", "Gerente"="Gerencia",
"Recursos_Humanos"="Recursos humanos", "Representante_Salud"="Representante de salud",
"Representante_Ventas"="Representante de ventas", "Tecnico_Laboratorio"="Técnico de laboratorio")
d$role_label <- unname(role_labels[as.character(d$job_role)])
save_plot("income_job_role", ggplot(d, aes(reorder(role_label, monthly_income, median), monthly_income/1e6)) +
geom_boxplot(fill = "#A8BDC9", outlier.alpha = .25) + coord_flip() +
labs(title = "El cargo distingue distribuciones salariales", x = NULL, y = "Ingreso mensual (millones de pesos)"), height = 6)
save_plot("income_education_overtime", ggplot(d, aes(education_category, monthly_income/1e6, fill = overtime)) +
geom_boxplot(outlier.alpha = .15) + scale_fill_manual(values = c("No"="#A8BDC9", "Si"="#D9AA8A")) +
scale_x_discrete(labels = c(bachelor="Pregrado", below_college="Inferior a\ncollege", college="College",
master="Maestría", doctorate="Doctorado")) +
labs(title = "Ingreso por educación y horas extra", x = "Categoría educativa", y = "Ingreso mensual (millones de pesos)", fill = "Horas extra"))
save_plot("gender_role_composition", ggplot(d, aes(role_label, fill = gender)) + geom_bar(position = "fill") +
coord_flip() + scale_y_continuous(labels = scales::label_percent(accuracy = 1)) +
scale_fill_manual(values = c("M"="#A8BDC9", "F"="#D9AA8A")) +
labs(title = "Composición por género dentro de cada cargo", x = NULL, y = "Participación", fill = "Género"), height = 6)
comparison <- result$cv_summary
save_plot("cross_validation_mae", ggplot(comparison, aes(reorder(model, mae), mae/1e6)) +
geom_point(size = 2.8, color = "#315B7D") +
geom_errorbar(aes(ymin = min_repeat_mae/1e6, ymax = max_repeat_mae/1e6), width = .15, color = "#315B7D") +
coord_flip() + labs(title = "Error de predicción con particiones comunes", subtitle = "MAE promedio; barras: mínimo y máximo entre cinco repeticiones",
x = "Modelo", y = "MAE (millones de pesos)"), height = 6)
diagnostic_data <- do.call(rbind, lapply(unique(c("M2", "L2", "Q2", result$selected_id)), function(id) {
fit <- result$fits[[id]]
data.frame(model = id, fitted = fitted(fit), residual = residuals(fit)/sigma(fit),
scale_location = sqrt(abs(rstandard(fit))), cooks_distance = cooks.distance(fit), record_id = d$record_id)
}))
diagnostic_data$fitted_display <- ifelse(diagnostic_data$model == "L2", diagnostic_data$fitted,
diagnostic_data$fitted/1e6)
save_plot("residual_diagnostics", ggplot(diagnostic_data, aes(fitted_display, residual)) +
geom_point(alpha = .2, size = .8, color = "#315B7D") + geom_hline(yintercept = 0, linetype = "dashed") +
geom_smooth(method = "loess", formula = y~x, se = FALSE, color = "#B14B35") + facet_wrap(~model, scales = "free_x") +
labs(title = "Linealidad y residuos", x = "Ajustado (millones; L2 en logaritmos)", y = "Residuo / desviación residual"), height = 6.3)
save_plot("residual_qq", ggplot(diagnostic_data, aes(sample = residual)) + stat_qq(alpha = .25, size = .8) +
stat_qq_line(color = "#B14B35") + facet_wrap(~model) +
labs(title = "Normalidad de los residuos", x = "Cuantil normal teórico", y = "Cuantil observado"), height = 6.3)
save_plot("scale_location", ggplot(diagnostic_data, aes(fitted_display, scale_location)) +
geom_point(alpha = .2, size = .8) + geom_smooth(method = "loess", formula = y~x, se = FALSE, color = "#B14B35") +
facet_wrap(~model, scales = "free_x") + labs(title = "Dispersión residual según el nivel ajustado",
x = "Ajustado (millones; L2 en logaritmos)", y = "Raíz del residuo estandarizado absoluto"), height = 6.3)
save_plot("influence_diagnostics", ggplot(diagnostic_data, aes(record_id, cooks_distance)) +
geom_col(fill = "#315B7D") + geom_hline(yintercept = 4/nrow(d), color = "#B14B35", linetype = "dashed") +
facet_wrap(~model) + labs(title = "Observaciones influyentes", subtitle = "La línea 4/n es una señal de revisión, no una regla de eliminación",
x = "Fila original", y = "Distancia de Cook"), height = 6.3)
g <- result$gender_results
g <- g[g$model %in% c("S3", "career_adjustment", "M2", "Q2"), ]
g$label <- factor(g$model, levels = c("Q2", "M2", "career_adjustment", "S3"),
labels = c("Con experiencia cuadrática", "Con cargo y horas extra", "Con trayectoria y educación", "Sin ajuste"))
save_plot("adjusted_gender_differences", ggplot(g, aes(label, estimate/1e6)) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_errorbar(aes(ymin = ci_lower/1e6, ymax = ci_upper/1e6), width = .15, color = "#315B7D") +
geom_point(size = 3, color = "#B14B35") + coord_flip() + labs(title = "Diferencia estimada: mujeres menos hombres",
subtitle = "Intervalos de confianza del 95% con covarianza HC3", x = NULL, y = "Diferencia mensual (millones de pesos)"))
s <- result$scenarios
s$label <- c("Perfil de laboratorio", "Perfil de gerencia")
save_plot("scenario_predictions", ggplot(s, aes(label, prediction/1e6)) +
geom_errorbar(aes(ymin = lower/1e6, ymax = upper/1e6), width = .18, color = "#315B7D") +
geom_point(size = 3.5, color = "#B14B35") + coord_flip() +
labs(title = "Dos perfiles: predicción e incertidumbre", subtitle = "Intervalos conformales del 95% nominal, calibrados en datos reservados",
x = NULL, y = "Ingreso mensual (millones de pesos)"))
invisible(NULL)
}# Save the complete functions above as income_analysis.R in the project folder.
source("income_analysis.R", encoding="UTF-8")
analysis <- run_income_analysis()
create_income_figures(analysis)sessionInfo()## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26300)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8 LC_CTYPE=Spanish_Colombia.utf8
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_4.0.3
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 generics_0.1.4 lattice_0.23-1 digest_0.6.39
## [5] magrittr_2.0.5 evaluate_1.0.5 grid_4.6.1 RColorBrewer_1.1-3
## [9] fastmap_1.2.0 Matrix_1.7-6 cellranger_1.1.0 jsonlite_2.0.0
## [13] Formula_1.2-6 mgcv_1.9-4 scales_1.4.0 textshaping_1.0.5
## [17] jquerylib_0.1.4 abind_1.4-8 cli_3.6.6 rlang_1.3.0
## [21] splines_4.6.1 withr_3.0.3 cachem_1.1.0 yaml_2.3.12
## [25] otel_0.2.0 tools_4.6.1 dplyr_1.2.1 vctrs_0.7.3
## [29] R6_2.6.1 stats4_4.6.1 zoo_1.9-0 lifecycle_1.0.5
## [33] car_3.1-5 ragg_1.5.2 pkgconfig_2.0.3 pillar_1.11.1
## [37] bslib_0.12.0 gtable_0.3.6 glue_1.8.1 systemfonts_1.3.2
## [41] xfun_0.61 tibble_3.3.1 lmtest_0.9-40 tidyselect_1.2.1
## [45] rstudioapi_0.19.0 knitr_1.52 farver_2.1.2 nlme_3.1-171
## [49] htmltools_0.5.9 rmarkdown_2.32 carData_3.0-6 labeling_0.4.3
## [53] compiler_4.6.1 S7_0.2.2 readxl_1.5.0.1