Sesión 5
Profesor:Dr. José Manuel MAGALLANES REYES, Ph.D
Profesor Principal del Departamento de Ciencias Sociales, Sección de Ciencia Política y Gobierno.
Oficina 223 - Edificio CISEPA / ECONOMIA / CCSS
Telefono: (51) 1 - 6262000 anexo 4302
Correo Electrónico: jmagallanes@pucp.edu.pe
Modelando la Duración
1 Introducción
Estos datos nos dan información sobre ex convictos y el tiempo que demoran en volver a la cárcel (Allison 2010):
knitr::include_url("https://docs.google.com/spreadsheets/d/e/2PACX-1vQSlGaMI8Q8qlXI0Bp3m7BQcEh8ZLzaP7RymVtRYkg3ah1sZVlCi6-HmeKCic1RjfuH3gL_wrbMms88/pubhtml")Veámos cómo los ha traído R:
link="https://docs.google.com/spreadsheets/d/e/2PACX-1vQSlGaMI8Q8qlXI0Bp3m7BQcEh8ZLzaP7RymVtRYkg3ah1sZVlCi6-HmeKCic1RjfuH3gL_wrbMms88/pub?gid=1573532387&single=true&output=csv"
carcel=read.csv(link, stringsAsFactors = T)
str(carcel)## 'data.frame': 432 obs. of 10 variables:
## $ semanasLibre : int 20 17 25 52 52 52 23 52 52 52 ...
## $ fueArrestado : int 1 1 1 0 0 0 1 0 0 0 ...
## $ tuvoApoyoDinero : int 0 0 0 1 0 0 0 1 0 0 ...
## $ edad : int 27 18 19 23 19 24 25 21 22 20 ...
## $ esNegro : int 1 1 0 1 0 1 1 1 1 1 ...
## $ expLaboralPrevia : int 0 0 1 1 1 1 1 1 0 1 ...
## $ casado : int 0 0 0 1 0 0 1 0 0 0 ...
## $ libertadCondicional: int 1 1 1 1 1 0 1 1 0 0 ...
## $ vecesEnCarcel : int 3 8 13 1 3 2 0 4 6 0 ...
## $ nivelEduca : int 2 3 2 4 2 3 3 2 2 4 ...A partir de aquí, hagamos ajustes a los tipos de datos según la metadata:
carcel[,c(2,3,5,6,7,8)]=lapply(carcel[,c(2,3,5,6,7,8)], as.factor)
carcel$nivelEduca=as.ordered(carcel$nivelEduca)
str(carcel)## 'data.frame': 432 obs. of 10 variables:
## $ semanasLibre : int 20 17 25 52 52 52 23 52 52 52 ...
## $ fueArrestado : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 2 1 1 1 ...
## $ tuvoApoyoDinero : Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 2 1 1 ...
## $ edad : int 27 18 19 23 19 24 25 21 22 20 ...
## $ esNegro : Factor w/ 2 levels "0","1": 2 2 1 2 1 2 2 2 2 2 ...
## $ expLaboralPrevia : Factor w/ 2 levels "0","1": 1 1 2 2 2 2 2 2 1 2 ...
## $ casado : Factor w/ 2 levels "0","1": 1 1 1 2 1 1 2 1 1 1 ...
## $ libertadCondicional: Factor w/ 2 levels "0","1": 2 2 2 2 2 1 2 2 1 1 ...
## $ vecesEnCarcel : int 3 8 13 1 3 2 0 4 6 0 ...
## $ nivelEduca : Ord.factor w/ 5 levels "1"<"2"<"3"<"4"<..: 2 3 2 4 2 3 3 2 2 4 ...El propósito de este trabajo es entender qué factores afectan la reincidencia de los presos dejados en libertad. Revisando la metadata, tenemos dos variables que pueden ser usadas como dependiente. Exploremos ambas:
- semanasLibre
summary(carcel$semanasLibre)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 50.00 52.00 45.85 52.00 52.00- fueArrestado
table(carcel$fueArrestado)##
## 0 1
## 318 114Con la data disponible podemos proponer estas hipotesis:
H1: El tiempo que permanece en libertad un exreo hasta que vuelve a la carcel está afectado si tuvo financiamiento,por su nivel educativo, y por sus encarcelamientos previos.
y…
H2: El que un ex reo vuelva a la carcel está afectado si tuvo financiamiento,por su nivel educativo, y por sus encarcelamientos previos.
La 1.1 nos muestra el resultado de modelar la H1 con una regresión Gaussiana (MAGALLANES 2022a), y la H2 con una regresión logística binaria (MAGALLANES 2022b):
#
h1=formula(semanasLibre~tuvoApoyoDinero+nivelEduca+vecesEnCarcel)
h2=formula(fueArrestado~tuvoApoyoDinero+nivelEduca+vecesEnCarcel)
#
rGauss=lm(h1,data=carcel)
rLogit=glm(h2,data=carcel, family = binomial)
#
models=list('Tiempo en Libertad (Gauss)'=rGauss, "Ser Arrestado (Logit)"=rLogit)
#
library(modelsummary)
modelsummary(models,
title = "Regresiones Gauss y Logit",
stars = TRUE,
output = "kableExtra")| Tiempo en Libertad (Gauss) | Ser Arrestado (Logit) | |
|---|---|---|
| (Intercept) | 48.474*** | −1.597*** |
| (1.297) | (0.295) | |
| tuvoApoyoDinero1 | 2.159+ | −0.473* |
| (1.206) | (0.227) | |
| nivelEduca.L | 1.061 | −0.618 |
| (2.965) | (0.769) | |
| nivelEduca.Q | 4.032 | −0.717 |
| (2.566) | (0.659) | |
| nivelEduca.C | −0.555 | 0.507 |
| (1.992) | (0.493) | |
| nivelEduca^4 | 0.551 | −0.049 |
| (1.423) | (0.321) | |
| vecesEnCarcel | −0.737*** | 0.097** |
| (0.212) | (0.037) | |
| Num.Obs. | 432 | 432 |
| R2 | 0.046 | |
| R2 Adj. | 0.033 | |
| AIC | 3413.9 | 488.7 |
| BIC | 3446.5 | 517.2 |
| Log.Lik. | −1698.952 | −237.357 |
| F | 3.431 | 3.601 |
| RMSE | 12.35 | 0.43 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
De la Tabla 1.1 vemos que para la Regresión Gaussiana:
- El tener apoyo financiero aumenta el tiempo en libertad con una significancia de 0.1.
- El nivel educativo no tendría efecto en el tiempo en libertad (probabilidad de efecto cero es mayor al 0.1).
- El tener encarcelamientos previos disminuye el tiempo en libertad con una significancia de 0.001.
Para el caso de los resultados con la Logit tenemos:
- El tener apoyo financiero disminuye el Log Odds Ratio de ser arrestado con una significancia de 0.05.
- El nivel educativo no tendría efecto en el Log Odds Ratio de ser arrestado (probabilidad de efecto cero es mayor al 0.1).
- El tener encarcelamientos previos aumenta el Log Odds Ratio de ser arrestado con una significancia de 0.01.
Nota que en la Tabla 1.1 ambos casos la variable ordinal viene con unas letras L-Q, etc. Esto informa si el efecto de la ordinal es lineal, cuadrático, cúbico, etc.
¿Qué problema no estamos advirtiendo?
Que la duración de la libertad está condicionada a ser arrestado, ambas son un todo. Una situación dura hasta que algo sucede y la situación cambia.
Que el hecho de ser arrestado es un evento, NO una característica, y, lo que es más, el NO ser arrestado es algo que puede variar en el tiempo, sólo que la investigación acabó y el liberado aun seguía libre.
En situaciones que combinan duración y observación de eventos, debemos usar el EHA o Análisis de Eventos Históricos (Grace-Martin 2010). Esta técnica funciona de tal manera que puede lidiar con el hecho de no darse el evento, en este caso, no ser arrestado: esto representa un caso censurado.
2 Analizando Eventos Históricos
El primer paso para usar EHA, es indicarle a R que trate a la data de esa manera. Creemos una nueva columna survival en nuestra tabla actual.
library(survival)
# note que necesito el factor como numérico
carcel$survival=with(carcel,Surv(time = semanasLibre,event = as.numeric(fueArrestado)))
# que es:
library(magrittr) # needed for pipe %>%
carcel%>%
rmarkdown::paged_table()Como ves, la columna creada tiene valores con un +, lo que indica que están censurados.
2.1 Análisis Kaplan-Meier (KM)
KM es el procedimiento descriptivo básico que se utiliza para ver la dinámica de sobrevivencia (Goel, Khanna, and Kishore 2010). Así, veamos en la Figura 2.1 nos muestra el comportamiento genérico de permanecer libre.
library(ggplot2)
library(ggfortify)
#aqui el generico
KM.generico = survfit(survival ~ 1, data = carcel)
###graficando
ejeX='SEMANAS\n curva cae cuando alguien es arrestado'
ejeY='Probabilidad \n(PERMANECER LIBRE)'
titulo="Curva de Sobrevivencia: permanecer libre"
autoplot(KM.generico,xlab=ejeX,ylab=ejeY, main = titulo,conf.int = F)
Figure 2.1: Análisis Kaplan-Meier (simple)
La Figura 2.1 nos da una idea de cómo se comporta esta población; por ejemplo, la gráfica nos dice que si pasan 40 semanas, la probabilidad de seguir libre está cerca al 80%.
El análisis KM es más interesante para ver una comparación:
KM_H1=formula(survival ~ tuvoApoyoDinero)
KM.fondos = survfit(KM_H1, data = carcel)
###
ejeX='SEMANAS\n curva cae cuando alguien es arrestado'
ejeY="Prob ('seguir libre')"
titulo="Curva de Sobrevivencia: ¿Beneficia el apoyo financiero?"
autoplot(KM.fondos,xlab=ejeX,ylab=ejeY,
main = titulo,conf.int = F) +
labs(colour = "Apoyo Financiero?") +
scale_color_discrete(labels = c("No", "Sí"))
Figure 2.2: Análisis Kaplan-Meier (grupos)
La posición de las curvas nos hace pensar que le cuesta más a los que no tuvieron fondos permanecer libres. Para una mayor confianza en tal hipótesis, podemos hacer la prueba de Mantel-Cox (LogRanK), obteniendo este nivel de significancia:
LogRank=survdiff(KM_H1, data = carcel)
# ver p-valor
LogRank$pvalue## [1] 0.05011612Siendo la H0 de KM: no hay diferencias entre grupos, con el p-valor obtenido la diferencia no es significativa al 0.05 (pero sí al 0.1). La Figura 2.3 aclara por qué.
autoplot(KM.fondos,xlab=ejeX,ylab=ejeY,
main = titulo,conf.int = T)+
labs(colour = "Apoyo Financiero?") +
scale_color_discrete(labels = c("No", "Sí"))
Figure 2.3: Curva de Sobrevivencia (grupos con intervalo de confianza)
De nuevo, como sólo hay dos variables, es difícil saber qué más interviene. De ahí que necesitamos un modelo que permita análisis multivariado.
2.2 Regresión de Cox
Como toda regresión, esta técnica permite utilizar regresores o predictores o covariados, pero no modela la duración sino el riesgo de que el evento suceda (ser re arrestado) (Abd ElHafeez et al. 2021). Veamos la Tabla 2.1.
COX_H1= formula(survival~tuvoApoyoDinero+nivelEduca+vecesEnCarcel)
#regression
rcox1 <- coxph(COX_H1,data=carcel)
modelcox=list('Riesgo - Re arrestado'=rcox1,'Riesgo- Re arrestado (exponenciado)'=rcox1)
#f <- function(x) format(x, digits = 4, scientific = FALSE)
library(modelsummary)
modelsummary(modelcox,
#fmt=f,
exponentiate = c(F,T),
statistic = 'conf.int',
title = "Regresión Cox",
stars = TRUE,
output = "kableExtra")| Riesgo - Re arrestado | Riesgo- Re arrestado (exponenciado) | |
|---|---|---|
| tuvoApoyoDinero1 | −0.415* | 0.660* |
| [−0.789, −0.042] | [0.454, 0.959] | |
| nivelEduca.L | −0.558 | 0.573 |
| [−1.976, 0.861] | [0.139, 2.366] | |
| nivelEduca.Q | −0.694 | 0.500 |
| [−1.907, 0.519] | [0.149, 1.680] | |
| nivelEduca.C | 0.396 | 1.486 |
| [−0.506, 1.298] | [0.603, 3.662] | |
| nivelEduca^4 | −0.008 | 0.992 |
| [−0.582, 0.566] | [0.559, 1.762] | |
| vecesEnCarcel | 0.087** | 1.091** |
| [0.032, 0.142] | [1.032, 1.152] | |
| Num.Obs. | 432 | 432 |
| AIC | 1338.1 | 1338.1 |
| BIC | 1362.6 | 1362.6 |
| RMSE | 0.51 | 0.51 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
La Tabla 2.1 muestra los coeficientes originales y exponenciados (las razones de riesgo o hazard ratios) . Veamos cada uno:
- Se puede inferir que dar financiamiento disminuye el riesgo de ser re arrestado con una significancia del 0.05. La columna a la derecha (los HRs) nos da el promedio de ese efecto en el riesgo como la distancia absoluta al uno (1); en este caso:
(apoyoDinero=abs(1-exp(coef(rcox1)[1])))## tuvoApoyoDinero1
## 0.3399735De ahí que usando el hazard ratio, podemos sostener que el riesgo de volver a la cárcel para alguien con apoyo financiero es en promedio 34% menor a la de los que no lo tienen.
Se puede inferir que es muy poco probable que el nivel educativo tenga efecto en el riesgo de ser re arrestado.
Se puede inferir que haber estado previamente en la cárcel aumenta el riesgo de ser re arrestado con una significancia del 0.01. El promedio de ese riesgo es la distancia absoluta al uno (1); en este caso:
(carcelantes=abs(1-exp(coef(rcox1)[6])))## vecesEnCarcel
## 0.0907188De ahí que, podemos sostener que el riesgo de volver a la cárcel para alguien con apoyo financiero es en promedio 9.07% mayor cada vez que aumenta en uno los años de encarcelamientos previos1.
Los hazard ratios de la Tabla 2.1 (columna derecha) podemos verlos gráficamente en la Figura 2.4.
library(survminer)
ggforest(rcox1, data = carcel,main = "¿Quiénes tienen mayor riesgo de volver a ser encarcelados?")
Figure 2.4: Los Hazard Ratios
2.3 Comparación
Podemos plantear un modelo sin educación. Veamos la Tabla 2.2.
COX_H1= formula(survival~tuvoApoyoDinero+nivelEduca+vecesEnCarcel)
COX_H2= formula(survival~tuvoApoyoDinero+vecesEnCarcel)
#regression
rcox2 <- coxph(COX_H2,data=carcel)
modelcox=list('Riesgo de Re arrestado (I)'=rcox2,'Riesgo de Re arrestado (II)'=rcox1)
#f <- function(x) format(x, digits = 4, scientific = FALSE)
library(modelsummary)
modelsummary(modelcox,
#fmt=f,
exponentiate = T,
statistic = 'conf.int',
title = "Regresión Cox (sólo Hazard Ratios)",
stars = TRUE,
output = "kableExtra")| Riesgo de Re arrestado (I) | Riesgo de Re arrestado (II) | |
|---|---|---|
| tuvoApoyoDinero1 | 0.671* | 0.660* |
| [0.462, 0.974] | [0.454, 0.959] | |
| vecesEnCarcel | 1.110*** | 1.091** |
| [1.053, 1.169] | [1.032, 1.152] | |
| nivelEduca.L | 0.573 | |
| [0.139, 2.366] | ||
| nivelEduca.Q | 0.500 | |
| [0.149, 1.680] | ||
| nivelEduca.C | 1.486 | |
| [0.603, 3.662] | ||
| nivelEduca^4 | 0.992 | |
| [0.559, 1.762] | ||
| Num.Obs. | 432 | 432 |
| AIC | 1338.3 | 1338.1 |
| BIC | 1346.4 | 1362.6 |
| RMSE | 0.51 | 0.51 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Vemos que las medidas globales de ajuste no son tan diferentes, por lo que es mejor evaluar la significancia de esas diferencias. Veamos la Tabla 2.3.
anova(rcox2,rcox1)%>%
knitr::kable(caption = "Tabla anova para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)| loglik | Chisq | Df | Pr(>|Chi|) |
|---|---|---|---|
| -667.1444 | NA | NA | NA |
| -663.0728 | 8.143309 | 4 | 0.0864674 |
Así, añadir nivel educativo no es significativo al 0.05, pero sí lo es al 0.1. Esto daría espacio para sostener la inclusión de nivel educativo en el modelo.