Modelos predictivos en la justicia - Parte 1

Datos y feature engineering

Introducción

En esta serie de documentos vamos a presentar un ejemplo de implementación de modelos predictivos aplicados el servicio de justicia con el propósito de predecir la cantidad de sentencias que producirá una jurisdicción completa (es decir, una unidad territorial del Poder Judicial) en un determinado período. Entendemos que este ejemplo -con ciertos ajustes- puede extenderse a otros indicadores relevantes para la actividad judicial, por caso predecir cuántas causas se iniciará en un período y lugar determinado, brindando así información de mucho valor para la organización del servicio.

Dado que nos interesa compartir la implementación y divulgar el uso de estas herramientas (de allí el open source) nos hemos ocupado fundamentalmente en la explicación y documentación más que en en el desarrollo de modelos aptos para producción. No obstante lo cual, para ésto último no descartamos alguna publicación próxima.

Trabajos vinculados

Elaborar pronósticos o proyecciones sobre variables importantes de un negocio es algo común en el mundo privado, aunque no abundan los casos documentados en el ámbito público judicial. Sin perjuicio de ello, contamos con abundante material que aborda el tema ‘proyecciones’ (forecast) en distintos dominios a partir de modelos estadísticos y machine learning (Alsharef, A., Aggarwal, K., Sonia et al.). Además de lo anterior, una gran cantidad de recursos se encuentra disponible en formato abierto y de libre acceso en github.

Código y ambiente de trabajo en R

Los scripts completos de esta implementación y los documentos RMarkdown que sirven a esta presentación puede accederse en mi repositorio de github castillosebastian. En el README de ese repositorio se puede encontrar un detalle de la configuración del ambiente de trabajo que he utilizado en el proyecto y recomendaciones para su instalación.

Como veremos más adelante, algunas partes de las operaciones de desarrollo de modelos (vinculadas a un pipeline de MLOps) requieren cierto hardware particular (sobretodo RAM y CPU). En nuestro caso hemos trabajado en un Centos 7, con 12 vCPU y 32G de RAM, y en ocasiones esos recursos fueron insuficientes, por lo que algo más de disponibilidad será de ayuda.

Respecto del pipeline de trabajo nos hemos apoyado exclusivamente en R, y particularmente en tidymodels, gracias al aprecio que tenemos por tidyverse. Este framework nos brinda toda una batería de herramientas que facilitan el trabajo y la experimentación, aspecto fundamental en el desarrollo de modelos de machine learning. Para culquiera intersado en él dirigirse aquí: https://www.tidymodels.org/.

El repositorio está armado de tal forma que los script puedan correrse de dos formas:

• a través se sus documentos RMarkdown, o
• mediante lso script .R en /source.

Librerías y algunas funciones de ayuda

A continuación presentamos las librerías que emplearemos en este documento. Cada nueva publicación de la serie incluirá un detalle de sus recursos.

pacman::p_load(tidyverse,
               timetk,
               tsibble,
               tsibbledata,
               fastDummies,
               skimr, 
               tidymodels, 
               kableExtra)

options(scipen = 999)

myskim <- skim_with(numeric = sfl(max, min), append = TRUE)

Explorando los datos

Para este ejemplo de implementación hemos generado un dataset de sentencias mensuales dictadas por materia, información de acceso público suministrada por el Superior Tribunal de Justicia de Entre Ríos -accesible aquí. La información comprende la producción de sentencias en primera instancia a lo largo del período que se extiende entre los años 2018 y septiembre 2022.

Las materias antes referidas se vinculan a órganos jurisdiccionales, así, por ejemplo, existen 38 juzgados con materia civil-comercial en la provincia, por lo que las sentencias de dicha materia refieren a la suma de la producción de los 38 juzgados. Para consultar extensamente la estructura y estadística pública judicial remitimos al link includio más arriba. Dicha información se registra digitalmente y es procesada por el Área de Planificación, Gestión y Estadística.

Los datos temporalmente organizados constituyen una serie temporal y su análisis es un área de estudio en sí mismo. Para una buena aproximación remitimos a “Forecasting: Principles and Practice” (3rd ed), de Rob J Hyndman and George Athanasopoulos, accesible aquí.

Leemos los datos:

produccion_tbl <- readRDS(paste0(RAW_DATA_DIR, "/produccion.rds")) %>% 
  mutate(materia = as.factor(materia)) %>% 
  tk_tbl()

Estrectura de los datos

Revisaremos la distribución de los datos agrupados por materia. La función skim nos brinda un excelente recurso para el análisis exploratorio:

myskim(produccion_tbl %>% group_by(materia)) %>%  yank("numeric") %>% select(-n_missing) %>% kable() %>% kable_styling()

skim_variable	materia	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist	max	min
sentencias_dictadas	civil-comercial	1	460.67308	141.98878	198	370.00	457.5	554.25	891	▃▇▅▂▁	891	198
sentencias_dictadas	familia	1	737.73077	197.23174	306	555.25	769.0	878.50	1087	▃▃▆▇▂	1087	306
sentencias_dictadas	laboral	1	343.86538	101.76348	43	283.25	345.5	423.00	512	▁▂▆▇▆	512	43
sentencias_dictadas	paz_2°-3°	1	328.01923	567.41497	2	28.00	51.0	197.25	2161	▇▁▁▁▁	2161	2
sentencias_dictadas	penal	1	406.82692	149.57216	33	350.50	430.0	487.25	716	▂▂▇▆▂	716	33
sentencias_dictadas	quiebra-ejecuciones	1	429.26923	232.13072	105	254.00	394.0	529.75	1190	▇▇▃▁▁	1190	105
sentencias_dictadas	contenc-admn	1	26.10417	14.28768	5	14.00	25.0	34.00	63	▇▆▆▂▂	63	5
sentencias_dictadas	paz_1°	1	887.12195	466.93400	293	588.00	742.0	1094.00	2351	▇▆▃▁▁	2351	293

El dataset tiene tres variables:

• mes: es la marca temporal, registrada como año-mes-dia,
  
• materia: contiene la materia de las sentencias dictadas. Como digimos estas materias ad-hoc se corresponden con la estructura orgánica del servicio de justicia (ej: civil-comercial corresponde a juzgados civiles y comerciales, paz_1 corresponde a juzgados de paz de primera categoría, y así para las demás), y
  
• sentencias_dictadas: contiene la cantidad de sentencias.

Sentencias mensuales dictadas en Primera Instancia en el STJER

Veamos a continuación las series temporales por materia de las sentencias dictadas en el período.

produccion_tbl %>% 
  group_by(materia) %>% 
  plot_time_series(mes, sentencias_dictadas,
                   .facet_ncol = 2,
                   .facet_scales = "free",
                   .interactive = T)

En el gráfico podrían verse tendencias positivas pronunciadas (eg.familia), tendencias casi estacionarias (civil-comercial) y tenencias negativas (paz_2_3). También podríamos agrupar las materias por similitud en su evolución, y así podríamos tratar de agrupar las siguientes materias: grupo1) familia, laboral y penal, grupo2) paz_1 y paz_2-3, y grupo3) contenc-adm, civil-comercial y quiebra-ejecuciones. Se ve también estacionalidad asociada a los periodos de actividad y recesos judiciales (que se extienden durante todo el mes de enero y mitad de julio). Estamos viendo 5 años dado que no se disponen de datos previos, lo cual es un límite importante a considerar en el análisis de la serie temporal. Vemos también la importante caida que tuvo el indicador en abril-mayo del año 2020 producto de la Pandemia COVID-19.

Serie con ajuste estacional

El presente gráfico muestra las series temporales con ajuste estacional donde buscamos resaltar la varición interanual seperándola de las fluctuaciones mensuales.

produccion_tbl %>% 
  group_by(materia) %>% 
  tk_stl_diagnostics(.date_var = mes, .value = sentencias_dictadas) %>% 
  plot_time_series(mes, seasadj,
                   .facet_ncol = 2,
                   .facet_scales = "free",
                   .interactive = T)

Sentencias dictadas por trimestre

Veamos ahora las sentencias totales agrupadas por trimestre.

produccion_tbl %>% 
  group_by(mes) %>% summarise(sentencias_dictadas = sum(sentencias_dictadas)) %>% 
  mutate(
    quarter = str_c("Trimestre ", as.character(quarter(mes)))
    ) %>%
    plot_time_series(
      .date_var = mes,
      .value = sentencias_dictadas,
      .facet_vars = quarter,
      .facet_ncol = 4, 
      .color_var = quarter, 
      .facet_scales = "fixed",
      .interactive = FALSE,
      .legend_show = FALSE,
      .title = "Seasonal Plot"
      )

Sentencias mensuales por materia

En el siguiente gráfico de densidad vemos las sentencias según su materia. Vemos que hay materias con una distribución aproximadamente normal (ej. quiebras y paz) y otras con distribuciones claramente asimétricas (penal y laboral).

produccion_tbl %>% 
  ggplot()+
  geom_density(aes(x = sentencias_dictadas, fill = materia), alpha = 0.5)+
  scale_x_log10()+
  theme(legend.position = "note") +
  labs(fill = "Circunscripcion", x = "Presentaciones Diarias") + 
  facet_wrap( ~ materia, scales = "free") +
  theme(strip.text.x = element_text(size = 16))

Autocorrelación

Así como la correlación mide el alcance de una relación lineal entre dos variables, la autocorrelación mide la relación lineal entre los valores rezagados (lags) de una serie de tiempo. Por ese motivo evaluar la autocorrelación de las series temporales es un buen primer paso para desarrollar nuestra intuición acerca de las posibilidades de modelarlas.

Para exponerlo de manera simple el gráfico ACF (Autocorrelation Function) explica la correlación que existe entre el valor presente de nuestra variable de interés y su valor pasado (t-1, t-2…t-n). El gráfico ACF muestra sobre el eje y la correlación y en el eje x los lags. El gráfico PACF (Partial autocorrelation function) explica la correlación parcial entre la serie y los lags en sí.

En el caso que estudiamos, calculamos los lags de 24 meses. Vemos que aparece en la correlación la estacionalidad anual en muchas materias (ej. civil, familia, laboral). En otras, ese patrón no aparece y la correlación decrece gradualmente (ej. paz_2_3 y cont-adm). Esto claramente es un desafío importante para nuestros modeles, que deberán aprender a diferenciar el comportamiento de las distitas materias.

produccion_tbl %>%  
  group_by(materia) %>% 
  plot_acf_diagnostics(mes, sentencias_dictadas,
                       .lags = 36,
                       .facet_scales = "free",
                       .interactive = T)

Modelo base: regresión lineal múltiple

A continuación generamos un modelo de regresión con los datos originales que nos servirá de base para evaluar las trasnformaciones que realizaremos buscando mejorar el poder predictivo de nuestros datos.

model = lm(sentencias_dictadas ~., data = produccion_tbl )
jtools::summ(model)

Observations	401
Dependent variable	sentencias_dictadas
Type	OLS linear regression

F(8,392)	36.04
R²	0.42
Adj. R²	0.41

	Est.	S.E.	t val.	p
(Intercept)	2354.88	540.63	4.36	0.00
mes	-0.10	0.03	-3.51	0.00
materiafamilia	277.06	55.67	4.98	0.00
materialaboral	-116.81	55.67	-2.10	0.04
materiapaz_2°-3°	-132.65	55.67	-2.38	0.02
materiapenal	-53.85	55.67	-0.97	0.33
materiaquiebra-ejecuciones	-31.40	55.67	-0.56	0.57
materiacontenc-admn	-427.81	56.85	-7.53	0.00
materiapaz_1°	445.27	59.53	7.48	0.00
Standard errors: OLS

#summary(model)

Iniciamos la ingeniería de variables (feature engineering)

En primer lugar, controlemos la regularidad de la serie. Para ello utilizaremos la función tk_summary_diagnostics(), que arroja información de la estructura de la serie por materia.

produccion_tbl %>%
  group_by(materia) %>%
  tk_summary_diagnostics(.date_var = mes) %>% 
  select(1:6, 10) %>% 
  kable() %>% kable_styling()

materia	n.obs	start	end	units	scale	diff.median
civil-comercial	52	2018-02-01	2022-09-01	days	month	2678400
familia	52	2018-02-01	2022-09-01	days	month	2678400
laboral	52	2018-02-01	2022-09-01	days	month	2678400
paz_2°-3°	52	2018-02-01	2022-09-01	days	month	2678400
penal	52	2018-02-01	2022-09-01	days	month	2678400
quiebra-ejecuciones	52	2018-02-01	2022-09-01	days	month	2678400
contenc-admn	48	2018-06-01	2022-09-01	days	month	2678400
paz_1°	41	2019-02-01	2022-09-01	days	month	2678400

Vemos que la serie tiene, en efecto, estructura mensual (scale = month) y que hay materias que tiene menos observaciones: contencioso-administrativo y paz 1º (esto podría afectar las predicciones). Dada la poca cantidad de datos disponibles, utilizaremos el dataset completo, aunque podemos asumir que el período 2020 contendrá valores extremos debido al efecto de la pandemia sobre la actividad social en general y judicial en particular.

En machine learning la creación de nuevas variables es una estrategia elemental para explotar la información que contiene el dataset, aumentando sus señales en relación a una variable de interés (en general el target del modelo). Las posibilidades aquí son inagotables: desde transformaciones lineales hasta creación de nuevas variables basadas en combinaciones aleatorias. En todos los casos el conocimiento del dominio puede ser determinante. Una buena referencia para este tema puede consultarse aquí.

Como vemos en el siguiente gráfico, ésta actividad se encuentra en las primeras etapas del desarrollo de un pipeline para predicciones de series temporals (gráfico publicado en Meisenbacher et all). Esta estructura nos aporta una visión general del proyecto y representa una buena síntesis del desafío que enfrentamos.

Entonces, a continuación presentaremos una breve reseña de las transformaciones que implementaremos en los datos:

• Transformación logarítmica: La transformación logarítmica es un tipo de transformación matemática que puede ser útil para reducir la varianza en los datos. Esto puede ser pertinente cuando los datos evidencian una variación que aumenta o disminuye con la serie, y tiene como ventaja que los logaritmos son interpretables como cambios relativos (o porcentuales) en la escala original (ver pto. 3.2 Transformations and adjustments, acceso público aquí).
  
• Estandarización: Dado que tenemos tantas series temporales como materias a predecir estandarizaremos las cantidades para centrar las series (mean=0, sd=1). Esto supone que las observaciones tienen una distribución normal, cosa que como vimos no sucede, y es probable que los resultados no sean precisos. Sin embargo, en esta etapa de la experimentación procuramos mantener todas la hipótesis abiertas, y luego veremos qué funciona y qué no.
  
• Variables de calendario: Estas variables surgen como operaciones sobre la marca temporal y es muy fácil calcularlas Hay funciones que se encargan de producirlas por defecto (e.g. timetk::tk_augment_timeseries_signature()) y luego pueden desecharse aquellas que no son de utilidad o no aplican al caso bajo estudio. En nuestro caso como estamos tratando con datos mensuales, luego de aplicar la función debemos suprimir todas las columnas relacionadas con horas, días y semanas.
  
• Transformaciones de Fourier: Con esta transformación de los datos temporales se busca convertir la temporalidad en frecuencias, obteniendo un espectro con los distintos armónicos que están presentes en la serie (resaltando sus elementos y dinámica) (ver Series Fourier, en https://otexts.com/fpp3/useful-predictors.html). Esto claramente será de utilidad en nuestro trabajo dada la estacionalidad que vimos en la serie total. Para elegir los parámetros de transformación hicimos un grupo de experimentos, en lo que advertimos lo gran incidencia que tiene esta transformación en los modelos lineales con los que testeamos significancia. Nos quedamos con dos transformaciones (k = 3 y k = 4) que nos dieron el Adj_R2 más alto.
• Lags y Rolling Lags: Los lags permiten construir una versión desfasada de la serie temporal, retardando (o adelantando) los valores de la serie en la magnitud del desfase especificado.

materia <- unique(produccion_tbl$materia)

groups <- lapply(X = 1:length(materia), FUN = function(x){
  
  produccion_tbl %>%
    filter(materia == materia[x]) %>%
    arrange(mes) %>%
    mutate(sentencias_dictadas =  log1p(x = sentencias_dictadas)) %>%
    # estandarizacion
    mutate(sentencias_dictadas =  standardize_vec(sentencias_dictadas)) %>%
    # agregamos meses a futuro
    future_frame(mes, .length_out = "12 month", .bind_data = TRUE) %>%
    mutate(materia = materia[x]) %>%
    #tk_augment_fourier(.date_var = mes, .periods = c(2,3,4,6,12), .K = 1) %>%
    #tk_augment_fourier(.date_var = mes, .periods = c(2,3,4,6,12), .K = 2) %>%
    tk_augment_fourier(.date_var = mes, .periods = c(2,3,4,6,12), .K = 3) %>%
    tk_augment_fourier(.date_var = mes, .periods = c(2,3,4,6,12), .K = 4) %>%
    # tk_augment_fourier(.date_var = mes, .periods = 6, .K = 2) %>%
    # tk_augment_fourier(.date_var = mes, .periods = 12, .K = 1) %>%
    tk_augment_lags(.value = sentencias_dictadas, .lags = c(2,3,4,5,6,12,13)) %>%
    tk_augment_slidify(.value   = c(sentencias_dictadas_lag12, sentencias_dictadas_lag13),
                       .f       = ~ mean(.x, na.rm = TRUE), 
                       .period  = c(3, 6, 9, 12),
                       .partial = TRUE,
                       .align   = "center")
})

groups_fe_tbl <- bind_rows(groups) %>%
  rowid_to_column(var = "rowid")

El resultado de aplicar las transformaciones precedentes nos deja otro dataset (aumentado) con la siguiente conformación (vemos las últimas filas que corresponden con las fechas a futuro):

groups_fe_tbl %>% tail() %>% glimpse()

## Rows: 6
## Columns: 59
## $ rowid                             <int> 507, 508, 509, 510, 511, 512
## $ mes                               <date> 2023-04-01, 2023-05-01, 2023-06-01, …
## $ materia                           <fct> familia, familia, familia, familia,…
## $ sentencias_dictadas               <dbl> NA, NA, NA, NA, NA, NA
## $ mes_sin2_K1                       <dbl> -0.8978045, 0.8486443, -0.8486443, 0…
## $ mes_cos2_K1                       <dbl> -0.4403942, 0.5289640, -0.5289640, 0…
## $ mes_sin2_K2                       <dbl> 0.7907757, 0.8978045, 0.8978045, 0.9…
## $ mes_cos2_K2                       <dbl> -0.6121060, -0.4403942, -0.4403942, …
## $ mes_sin2_K3                       <dbl> 0.2012985, 0.1011683, -0.1011683, 0.…
## $ mes_cos2_K3                       <dbl> 0.9795299, -0.9948693, 0.9948693, -0…
## $ mes_sin3_K1                       <dbl> 0.6766273, 0.3630939, -0.9884683, 0.…
## $ mes_cos3_K1                       <dbl> 0.73632570, -0.93175257, 0.15142778,…
## $ mes_sin3_K2                       <dbl> 0.9964361, -0.6766273, -0.2993631, 0…
## $ mes_cos3_K2                       <dbl> 0.08435107, 0.73632570, -0.95413926,…
## $ mes_sin3_K3                       <dbl> 0.7907757, 0.8978045, 0.8978045, 0.9…
## $ mes_cos3_K3                       <dbl> -0.6121060, -0.4403942, -0.4403942, …
## $ mes_sin4_K1                       <dbl> -0.8486443, 0.4853020, 0.8743466, -0…
## $ mes_cos4_K1                       <dbl> 0.5289640, 0.8743466, -0.4853020, -0…
## $ mes_sin4_K2                       <dbl> -0.8978045, 0.8486443, -0.8486443, 0…
## $ mes_cos4_K2                       <dbl> -0.4403942, 0.5289640, -0.5289640, 0…
## $ mes_sin4_K3                       <dbl> -0.10116832, 0.99871651, -0.05064917…
## $ mes_cos4_K3                       <dbl> -0.99486932, 0.05064917, 0.99871651,…
## $ mes_sin6_K1                       <dbl> -0.3630939, -0.9827901, -0.6513725, …
## $ mes_cos6_K1                       <dbl> -0.9317526, -0.1847261, 0.7587581, 0…
## $ mes_sin6_K2                       <dbl> 0.6766273, 0.3630939, -0.9884683, 0.…
## $ mes_cos6_K2                       <dbl> 0.73632570, -0.93175257, 0.15142778,…
## $ mes_sin6_K3                       <dbl> -0.8978045, 0.8486443, -0.8486443, 0…
## $ mes_cos6_K3                       <dbl> -0.4403942, 0.5289640, -0.5289640, 0…
## $ mes_sin12_K1                      <dbl> 0.9827901, 0.7696512, 0.3473053, -0.…
## $ mes_cos12_K1                      <dbl> -0.1847261, -0.6384645, -0.9377521, …
## $ mes_sin12_K2                      <dbl> -0.3630939, -0.9827901, -0.6513725, …
## $ mes_cos12_K2                      <dbl> -0.9317526, -0.1847261, 0.7587581, 0…
## $ mes_sin12_K3                      <dbl> -0.8486443, 0.4853020, 0.8743466, -0…
## $ mes_cos12_K3                      <dbl> 0.5289640, 0.8743466, -0.4853020, -0…
## $ mes_sin2_K4                       <dbl> -0.9680771, -0.7907757, -0.7907757, …
## $ mes_cos2_K4                       <dbl> -0.2506525, -0.6121060, -0.6121060, …
## $ mes_sin3_K4                       <dbl> 0.1681009, -0.9964361, 0.5712682, 0.…
## $ mes_cos3_K4                       <dbl> -0.98576980, 0.08435107, 0.82076344,…
## $ mes_sin4_K4                       <dbl> 0.7907757, 0.8978045, 0.8978045, 0.9…
## $ mes_cos4_K4                       <dbl> -0.6121060, -0.4403942, -0.4403942, …
## $ mes_sin6_K4                       <dbl> 0.9964361, -0.6766273, -0.2993631, 0…
## $ mes_cos6_K4                       <dbl> 0.08435107, 0.73632570, -0.95413926,…
## $ mes_sin12_K4                      <dbl> 0.6766273, 0.3630939, -0.9884683, 0.…
## $ mes_cos12_K4                      <dbl> 0.73632570, -0.93175257, 0.15142778,…
## $ sentencias_dictadas_lag2          <dbl> NA, NA, NA, NA, NA, NA
## $ sentencias_dictadas_lag3          <dbl> NA, NA, NA, NA, NA, NA
## $ sentencias_dictadas_lag4          <dbl> NA, NA, NA, NA, NA, NA
## $ sentencias_dictadas_lag5          <dbl> NA, NA, NA, NA, NA, NA
## $ sentencias_dictadas_lag6          <dbl> NA, NA, NA, NA, NA, NA
## $ sentencias_dictadas_lag12         <dbl> 0.7374918, 0.6863639, 0.8373779, -0.…
## $ sentencias_dictadas_lag13         <dbl> 1.1549213, 0.7374918, 0.6863639, 0.8…
## $ sentencias_dictadas_lag12_roll_3  <dbl> 0.8595923, 0.7537445, 0.3136831, 0.5…
## $ sentencias_dictadas_lag13_roll_3  <dbl> 0.6397977, 0.8595923, 0.7537445, 0.3…
## $ sentencias_dictadas_lag12_roll_6  <dbl> 0.4767404, 0.7031172, 0.6501933, 0.6…
## $ sentencias_dictadas_lag13_roll_6  <dbl> 0.7479946, 0.4767404, 0.7031172, 0.6…
## $ sentencias_dictadas_lag12_roll_9  <dbl> 0.7035614, 0.6808771, 0.6353826, 0.7…
## $ sentencias_dictadas_lag13_roll_9  <dbl> 0.6311865, 0.7035614, 0.6613145, 0.6…
## $ sentencias_dictadas_lag12_roll_12 <dbl> 0.7184816, 0.7169431, 0.6808771, 0.6…
## $ sentencias_dictadas_lag13_roll_12 <dbl> 0.7064215, 0.7065919, 0.7035614, 0.6…

Prueba de significancia con modelo lineal

Hacemos una rápida exploración a partir de una regresión lineal múltiple para ver qué tanto podemos modelar los datos con una aproximación lineal, exploración que también nos dará un resultado de base para comparar los próximos modelos y una referencia de la importancia de los features generados. Advertimos que la transformación de Fourier y los Lags son significativas en este modelo (p_value < 0.05), aunque el modelo en general tenga un Adj-Rsq modesto de 0.65. Esto definitivamente indica la oportunidad de mejorar la obtención y selección de variables. Sin perjuicio de esta performance (lo que podría indicar tempranamente la necesidad de buscar mejores datos), seguiremos adelante con todos los features generados a ver si modelos de machine learning confirman o rechazan lo obtenido por el modelo lineal.

model = lm(sentencias_dictadas ~., data = groups_fe_tbl )
jtools::summ(model, confint = TRUE)

Observations	312 (200 missing obs. deleted)
Dependent variable	sentencias_dictadas
Type	OLS linear regression

F(44,267)	14.15
R²	0.70
Adj. R²	0.65

	Est.	2.5%	97.5%	t val.	p
(Intercept)	87.09	43.86	130.31	3.97	0.00
rowid	0.00	-0.00	0.00	0.00	1.00
mes	-0.00	-0.01	-0.00	-3.95	0.00
materiafamilia	0.38	0.20	0.57	4.16	0.00
materialaboral	0.30	0.09	0.52	2.76	0.01
materiapaz_2°-3°	-0.44	-0.66	-0.22	-3.98	0.00
materiapenal	0.46	0.24	0.69	4.05	0.00
materiaquiebra-ejecuciones	-0.04	-0.25	0.18	-0.34	0.73
mes_sin2_K1	-0.21	-0.35	-0.08	-3.10	0.00
mes_cos2_K1	0.06	-0.07	0.18	0.91	0.36
mes_sin2_K2	-0.74	-1.55	0.08	-1.78	0.08
mes_cos2_K2	-2.52	-3.32	-1.72	-6.19	0.00
mes_sin2_K3	0.06	-0.03	0.15	1.26	0.21
mes_cos2_K3	0.17	0.04	0.31	2.49	0.01
mes_sin3_K1	-0.32	-0.48	-0.16	-3.94	0.00
mes_cos3_K1	0.07	-0.06	0.20	1.04	0.30
mes_sin3_K2	-0.05	-0.17	0.07	-0.87	0.38
mes_cos3_K2	0.07	-0.04	0.18	1.26	0.21
mes_sin3_K3	NA	NA	NA	NA	NA
mes_cos3_K3	NA	NA	NA	NA	NA
mes_sin4_K1	0.21	0.03	0.39	2.31	0.02
mes_cos4_K1	-0.08	-0.22	0.06	-1.16	0.25
mes_sin4_K2	NA	NA	NA	NA	NA
mes_cos4_K2	NA	NA	NA	NA	NA
mes_sin4_K3	-0.13	-0.24	-0.03	-2.44	0.02
mes_cos4_K3	-0.15	-0.26	-0.03	-2.48	0.01
mes_sin6_K1	0.09	-0.05	0.23	1.32	0.19
mes_cos6_K1	-0.06	-0.18	0.06	-1.04	0.30
mes_sin6_K2	NA	NA	NA	NA	NA
mes_cos6_K2	NA	NA	NA	NA	NA
mes_sin6_K3	NA	NA	NA	NA	NA
mes_cos6_K3	NA	NA	NA	NA	NA
mes_sin12_K1	-0.12	-0.27	0.03	-1.58	0.12
mes_cos12_K1	0.70	0.50	0.90	6.93	0.00
mes_sin12_K2	NA	NA	NA	NA	NA
mes_cos12_K2	NA	NA	NA	NA	NA
mes_sin12_K3	NA	NA	NA	NA	NA
mes_cos12_K3	NA	NA	NA	NA	NA
mes_sin2_K4	-0.32	-0.56	-0.07	-2.56	0.01
mes_cos2_K4	0.28	-0.10	0.66	1.45	0.15
mes_sin3_K4	-0.02	-0.14	0.09	-0.36	0.72
mes_cos3_K4	0.16	0.04	0.27	2.70	0.01
mes_sin4_K4	NA	NA	NA	NA	NA
mes_cos4_K4	NA	NA	NA	NA	NA
mes_sin6_K4	NA	NA	NA	NA	NA
mes_cos6_K4	NA	NA	NA	NA	NA
mes_sin12_K4	NA	NA	NA	NA	NA
mes_cos12_K4	NA	NA	NA	NA	NA
sentencias_dictadas_lag2	0.17	0.04	0.29	2.67	0.01
sentencias_dictadas_lag3	-0.05	-0.19	0.09	-0.69	0.49
sentencias_dictadas_lag4	-0.15	-0.28	-0.02	-2.21	0.03
sentencias_dictadas_lag5	0.10	-0.04	0.25	1.39	0.16
sentencias_dictadas_lag6	-0.04	-0.30	0.22	-0.30	0.77
sentencias_dictadas_lag12	-0.13	-0.46	0.20	-0.77	0.44
sentencias_dictadas_lag13	-0.04	-0.36	0.27	-0.27	0.78
sentencias_dictadas_lag12_roll_3	0.15	-0.46	0.76	0.48	0.63
sentencias_dictadas_lag13_roll_3	0.06	-0.45	0.58	0.24	0.81
sentencias_dictadas_lag12_roll_6	0.10	-0.50	0.69	0.32	0.75
sentencias_dictadas_lag13_roll_6	-0.41	-1.21	0.39	-1.00	0.32
sentencias_dictadas_lag12_roll_9	-0.10	-1.43	1.24	-0.14	0.89
sentencias_dictadas_lag13_roll_9	0.58	-0.34	1.50	1.25	0.21
sentencias_dictadas_lag12_roll_12	-1.00	-2.98	0.99	-0.99	0.32
sentencias_dictadas_lag13_roll_12	0.94	-0.33	2.22	1.46	0.15
Standard errors: OLS

#summary(model)

Antes de pasar al próximo paso vamos a guardar los parametros de estandarización de las series para futuro tratamientos. Estamos trabajando con datos correspondientes a múltiples materias y buscaremos predecir 8 series de tiempo (una para cada materia) lo que nos deja con 16 parámetros de estandarización (media, desviación estándar) para las series. Guardamos estos parámetros para un futuro retratamiento de los datos.

tmp <- produccion_tbl %>%
  group_by(materia) %>% 
  arrange(mes) %>%
  mutate(sentencias_dictadas = log1p(x = sentencias_dictadas)) %>% 
  group_map(~ c(mean = mean(.x$sentencias_dictadas, na.rm = TRUE),
                sd = sd(.x$sentencias_dictadas, na.rm = TRUE))) %>% 
  bind_rows()

std_mean <- tmp$mean
std_sd <- tmp$sd
rm('tmp')

Particionamos el dataset para las próximas etapas

Dividimos el data set para entrenamiento y testeo, una estrategia fundamental en el ML para probar los modelos finales con datos no vistos durante su entrenamiento. Al mismo tiempo, generaremos los datos de los futuros meses que se intentan predecir y cuyas cantidades se desconoce. Los datos futuros tienen los valores estimados para todas las variables predictoras de nuestro dataset, y tienen nuestra variable target sin dato.

Un aspecto importante a tener en cuenta cuando particionamos es la proporción en la división de entrenamiento y testeo del dataset. Como sostienen Hyndman-Athanasopoulos: “el tamaño del conjunto de prueba suele ser aproximadamente el 20 % de la muestra total, aunque este valor depende de la duración de la muestra y de la anticipación con la que desea pronosticar. Idealmente, el conjunto de datos de prueba debería ser al menos tan grande como el horizonte de pronóstico máximo requerido”. Como nosotros buscaremos predecir 12 meses a futuro, dejaremos el dataset de test con 12 meses.

# preparo datasets futuro
data_prepared_tbl <- groups_fe_tbl %>%
  filter(!is.na(sentencias_dictadas)) %>%
  drop_na()

future_tbl <- groups_fe_tbl %>%
  filter(is.na(sentencias_dictadas))

splits <- data_prepared_tbl %>%
  time_series_split(mes, 
                    assess = "12 months", 
                    cumulative = TRUE)

splits

## <Analysis/Assess/Total>
## <224/88/312>

Según la partición propuesta vemos cómo quedan divididos los datos de entrenamiento y testeo para la materia civil-comercial:

splits %>%
  tk_time_series_cv_plan() %>%
  filter(materia == materia[1]) %>%
  plot_time_series_cv_plan(.date_var = mes, 
                           .value = sentencias_dictadas)

Del gráfico se puede advertir, a simple vista, la dificultad de lograr modelar esta serie en particular debido a sus importantes variaciones e irregularidad.

Creamos una receta para nuestros modelos

La creación de una receta (mediante la función recipe()) es un paso importante para automatizar el trabajo con nuestros modelos y nos ayuda a manipular los datos de manera segura. La receta, al igual que una receta de cocina, es una especificación sobre cómo tratar los ingredientes -en este caso los datos- antes y durante la preparación de un plato -aquí un modelo predictivo-. La diferencia estaría dada en que, mediante la receta no se altera las datos originales, que siempre permenecen accesibles. Como se detalla en la documentación de la librería tidymodels las recetas se construyen como una serie de pasos de preprocesamiento para realizar las acciones de transformación pensando en lo que requieren los modelos con los que estoy entrenando. Por ejemplo algunos modelos no toleran datos faltantes o bien datos categóricos o campos fecha, en tales casos mediante las funciones asociadas a recipe() se pueden excluir tales columnas de manera segura y práctica. Entre las posibilidades de tratamiento están:

• convertir predictores cualitativos en variables indicadoras (también conocidas como variables ficticias),
  
• transformar datos para que estén en una escala diferente (por ejemplo, tomando el logaritmo de una variable),
  
• transformando grupos enteros de predictores juntos,
• extraer características clave de variables sin procesar (por ejemplo, obtener el día de la semana de una variable de fecha),
  
• agregar características basadas en calendario: step_timeseries_signature()
• eliminar columnas: step_rm()
• realizar una codificación one-hot en variables categóricas: step_dummy()
• normalizar variables numéricas: step_normalize()

Comparados con las fórmulas que se emplean para definir los modelos en R, las recetas se pueden usar para hacer muchas de las mismas cosas, pero tienen una gama mucho más amplia de posibilidades.

recipe_spec <- recipe(sentencias_dictadas ~ ., data = training(splits)) %>%
  update_role(rowid, new_role = "indicator") %>%  
  step_other(materia) %>%
  step_timeseries_signature(mes) %>%
  step_rm(matches("(.xts$)|(.iso$)|(hour)|(minute)|(second)|(day)|(week)|(am.pm)")) %>%
  step_dummy(all_nominal(), one_hot = TRUE) %>%
  step_normalize(mes_index.num, mes_year)

El retratamiento de los datos nos deja un dataset como el que sigue (las 6 primeras filas)

# Recipe summary:
recipe_spec %>% 
  prep() %>%
  juice() %>% 
  head() %>% 
  glimpse() 

## Rows: 6
## Columns: 82
## $ rowid                             <int> 14, 78, 142, 206, 270, 334
## $ mes                               <date> 2019-04-01, 2019-04-01, 2019-04-01, …
## $ mes_sin2_K1                       <dbl> 0.6513725, 0.6513725, 0.6513725, 0.…
## $ mes_cos2_K1                       <dbl> 0.7587581, 0.7587581, 0.7587581, 0.7…
## $ mes_sin2_K2                       <dbl> 0.9884683, 0.9884683, 0.9884683, 0.9…
## $ mes_cos2_K2                       <dbl> 0.1514278, 0.1514278, 0.1514278, 0.1…
## $ mes_sin2_K3                       <dbl> 0.8486443, 0.8486443, 0.8486443, 0.8…
## $ mes_cos2_K3                       <dbl> -0.528964, -0.528964, -0.528964, -0.…
## $ mes_sin3_K1                       <dbl> 0.5432217, 0.5432217, 0.5432217, 0.5…
## $ mes_cos3_K1                       <dbl> -0.8395893, -0.8395893, -0.8395893, …
## $ mes_sin3_K2                       <dbl> -0.9121663, -0.9121663, -0.9121663, …
## $ mes_cos3_K2                       <dbl> 0.4098203, 0.4098203, 0.4098203, 0.4…
## $ mes_sin3_K3                       <dbl> 0.9884683, 0.9884683, 0.9884683, 0.9…
## $ mes_cos3_K3                       <dbl> 0.1514278, 0.1514278, 0.1514278, 0.1…
## $ mes_sin4_K1                       <dbl> 0.3473053, 0.3473053, 0.3473053, 0.3…
## $ mes_cos4_K1                       <dbl> 0.9377521, 0.9377521, 0.9377521, 0.9…
## $ mes_sin4_K2                       <dbl> 0.6513725, 0.6513725, 0.6513725, 0.6…
## $ mes_cos4_K2                       <dbl> 0.7587581, 0.7587581, 0.7587581, 0.7…
## $ mes_sin4_K3                       <dbl> 0.8743466, 0.8743466, 0.8743466, 0.8…
## $ mes_cos4_K3                       <dbl> 0.485302, 0.485302, 0.485302, 0.4853…
## $ mes_sin6_K1                       <dbl> -0.9590592, -0.9590592, -0.9590592, …
## $ mes_cos6_K1                       <dbl> -0.2832055, -0.2832055, -0.2832055, …
## $ mes_sin6_K2                       <dbl> 0.5432217, 0.5432217, 0.5432217, 0.5…
## $ mes_cos6_K2                       <dbl> -0.8395893, -0.8395893, -0.8395893, …
## $ mes_sin6_K3                       <dbl> 0.6513725, 0.6513725, 0.6513725, 0.6…
## $ mes_cos6_K3                       <dbl> 0.7587581, 0.7587581, 0.7587581, 0.7…
## $ mes_sin12_K1                      <dbl> 0.8010011, 0.8010011, 0.8010011, 0.8…
## $ mes_cos12_K1                      <dbl> -0.5986629, -0.5986629, -0.5986629, …
## $ mes_sin12_K2                      <dbl> -0.9590592, -0.9590592, -0.9590592, …
## $ mes_cos12_K2                      <dbl> -0.2832055, -0.2832055, -0.2832055, …
## $ mes_sin12_K3                      <dbl> 0.3473053, 0.3473053, 0.3473053, 0.3…
## $ mes_cos12_K3                      <dbl> 0.9377521, 0.9377521, 0.9377521, 0.9…
## $ mes_sin2_K4                       <dbl> 0.2993631, 0.2993631, 0.2993631, 0.2…
## $ mes_cos2_K4                       <dbl> -0.9541393, -0.9541393, -0.9541393, …
## $ mes_sin3_K4                       <dbl> -0.7476486, -0.7476486, -0.7476486, …
## $ mes_cos3_K4                       <dbl> -0.6640946, -0.6640946, -0.6640946, …
## $ mes_sin4_K4                       <dbl> 0.9884683, 0.9884683, 0.9884683, 0.9…
## $ mes_cos4_K4                       <dbl> 0.1514278, 0.1514278, 0.1514278, 0.1…
## $ mes_sin6_K4                       <dbl> -0.9121663, -0.9121663, -0.9121663, …
## $ mes_cos6_K4                       <dbl> 0.4098203, 0.4098203, 0.4098203, 0.4…
## $ mes_sin12_K4                      <dbl> 0.5432217, 0.5432217, 0.5432217, 0.5…
## $ mes_cos12_K4                      <dbl> -0.8395893, -0.8395893, -0.8395893, …
## $ sentencias_dictadas_lag2          <dbl> 1.31276320, 0.23255187, 0.40506319, …
## $ sentencias_dictadas_lag3          <dbl> 0.78148764, 0.26207786, 0.59135275, …
## $ sentencias_dictadas_lag4          <dbl> 0.589759717, 0.413805271, 0.70427530…
## $ sentencias_dictadas_lag5          <dbl> 1.0233139, 0.6679058, 0.9152942, 1.7…
## $ sentencias_dictadas_lag6          <dbl> 0.3035026, 0.3037943, 0.1618912, 1.4…
## $ sentencias_dictadas_lag12         <dbl> 0.4407968, -1.5932630, -0.4052926, 1…
## $ sentencias_dictadas_lag13         <dbl> 0.1464362, -1.4118520, -4.7741312, 1…
## $ sentencias_dictadas_lag12_roll_3  <dbl> 0.2754190, -1.4467564, -1.9066597, 1…
## $ sentencias_dictadas_lag13_roll_3  <dbl> 0.2936165, -1.5025575, -2.5897119, 1…
## $ sentencias_dictadas_lag12_roll_6  <dbl> 0.2166046, -1.2441192, -1.2943691, 1…
## $ sentencias_dictadas_lag13_roll_6  <dbl> 0.2808430, -1.3140124, -1.5055773, 1…
## $ sentencias_dictadas_lag12_roll_9  <dbl> -0.06722874, -1.48588734, -1.3568458…
## $ sentencias_dictadas_lag13_roll_9  <dbl> 0.2166046, -1.2441192, -1.2943691, 1…
## $ sentencias_dictadas_lag12_roll_12 <dbl> 0.03646295, -1.18046292, -0.97457860…
## $ sentencias_dictadas_lag13_roll_12 <dbl> -0.001685573, -1.392499675, -1.13693…
## $ sentencias_dictadas               <dbl> 1.430246153, 0.453616119, 0.64847990…
## $ mes_index.num                     <dbl> -1.638512, -1.638512, -1.638512, -1.…
## $ mes_year                          <dbl> -1.236077, -1.236077, -1.236077, -1.…
## $ mes_half                          <int> 1, 1, 1, 1, 1, 1
## $ mes_quarter                       <int> 2, 2, 2, 2, 2, 2
## $ mes_month                         <int> 4, 4, 4, 4, 4, 4
## $ materia_civil.comercial           <dbl> 1, 0, 0, 0, 0, 0
## $ materia_familia                   <dbl> 0, 1, 0, 0, 0, 0
## $ materia_laboral                   <dbl> 0, 0, 1, 0, 0, 0
## $ materia_paz_2..3.                 <dbl> 0, 0, 0, 1, 0, 0
## $ materia_penal                     <dbl> 0, 0, 0, 0, 1, 0
## $ materia_quiebra.ejecuciones       <dbl> 0, 0, 0, 0, 0, 1
## $ materia_other                     <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_01                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_02                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_03                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_04                  <dbl> 1, 1, 1, 1, 1, 1
## $ mes_month.lbl_05                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_06                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_07                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_08                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_09                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_10                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_11                  <dbl> 0, 0, 0, 0, 0, 0
## $ mes_month.lbl_12                  <dbl> 0, 0, 0, 0, 0, 0

Rescatamos los parámetros de la normalización min-max

Obtenemos los parámetros del dataset original para la normalización min-max de las variables mes_index.num y mes_year.

# mes_index.num and yar normalization Parameters
myskim <- skim_with(numeric = sfl(max, min), append = TRUE)

norm_param = recipe(sentencias_dictadas ~ ., data = training(splits)) %>%
  step_timeseries_signature(mes) %>% prep() %>% juice() %>%  myskim() %>% 
  dplyr::filter(skim_variable == "mes_index.num" | 
                  skim_variable == "mes_year") %>% 
  yank("numeric") %>% as_tibble()

# mes index
mes_index.num_limit_lower = norm_param$min[norm_param$skim_variable == "mes_index.num"]
mes_index.num_limit_upper = norm_param$max[norm_param$skim_variable == "mes_index.num"]
# mes_year normalization Parameters
mes_year_limit_lower = norm_param$min[norm_param$skim_variable == "mes_year"]
mes_year_limit_upper = norm_param$max[norm_param$skim_variable == "mes_year"]

Guardamos el trabajo

Finalmente metemos los datos, recetas, divisiones y parámetros en una lista y la guardaremos como un archivo RDS.

dir.create(paste0(HOME_DIR, "/exp/001/"), showWarnings = FALSE )
archivo_salida  <-  paste0(HOME_DIR,"/exp/001/feature_engineering_artifacts_list.rds")

feature_engineering_artifacts_list <- list(
  # Data
  data = list(
    data_prepared_tbl = data_prepared_tbl,
    future_tbl      = future_tbl,
    materia = materia
  ),
  
  # Recipes
  recipes = list(
    recipe_spec = recipe_spec
  ),
  
  # Splits
  splits = splits,
  
  # Inversion Parameters
  standardize = list(
    std_mean = std_mean,
    std_sd   = std_sd
  ),
  
  normalize = list(
    mes_index.num_limit_lower = mes_index.num_limit_lower, 
    mes_index.num_limit_upper = mes_index.num_limit_upper,
    mes_year_limit_lower = mes_year_limit_lower,
    mes_year_limit_upper = mes_year_limit_upper
  )  
)


feature_engineering_artifacts_list %>% 
  write_rds(archivo_salida)

Conclusion

A lo largo del documento presentamos los primeros pasos en un pipeline para la creación de modelos predictivos de indicadores judiciales, en este caso sentencias dictadas por materia.

Hemos trabajando con el framework que nos brinda tidymodels y otras librerías importantes de R que manipulan series temporales. Con ello, hemos aumentado nuestros datos y probado la significancia de las nuevas variables.

Aplicamos una serie de transformaciones que tuvieron buenos resultados en la evaluación realizada en un modelo lineal de regresión múltiple, que llevaron a mejorar el R2-Ajustado de 0.4 a 0.65. Entre esas trasnformaciones la series de Fourier y los lags fueron los que tuvieron mayor valor predictivo.

En los próximos documentos avanzaremos en el armado de workflows de procesamiento con modelos de machine learning: alphabet, xgboost, random forest y ligthgbm. Luego seguiremos con el tuneo de hiperparámetros mediante optimizaciones bayesianas. A no desesperar en el camino; una vez que se realiza no hay vuelta atrás :)