EHA_cox.utf8.md

Análisis de Eventos Históricos (EHA)

Hablamos de EHA cuando queremos examinar y modelar la duración, es decir, el tiempo que demora algo en ocurrir.

link='https://raw.githubusercontent.com/PoliticayGobiernoPUCP/EstadisticaPoliticaGobiernoII/master/sesiones/data/reincidencia.txt'
eha=read.table(link,header = T, stringsAsFactors = T)
names(eha)

##  [1] "week"   "arrest" "fin"    "age"    "race"   "wexp"   "mar"   
##  [8] "paro"   "prio"   "educ"   "emp1"   "emp2"   "emp3"   "emp4"  
## [15] "emp5"   "emp6"   "emp7"   "emp8"   "emp9"   "emp10"  "emp11" 
## [22] "emp12"  "emp13"  "emp14"  "emp15"  "emp16"  "emp17"  "emp18" 
## [29] "emp19"  "emp20"  "emp21"  "emp22"  "emp23"  "emp24"  "emp25" 
## [36] "emp26"  "emp27"  "emp28"  "emp29"  "emp30"  "emp31"  "emp32" 
## [43] "emp33"  "emp34"  "emp35"  "emp36"  "emp37"  "emp38"  "emp39" 
## [50] "emp40"  "emp41"  "emp42"  "emp43"  "emp44"  "emp45"  "emp46" 
## [57] "emp47"  "emp48"  "emp49"  "emp50"  "emp51"  "emp52"

Los datos mostrados se refieren a casos de reincidencia criminal; y contienen lo siguiente:

week: Semana del primer arresto luego de haber salido de la prisión. También representa el momento cuando acabo de hacerle seguimiento a este caso.
arrest: Categoría lógica, este es el evento. Vale 1 si al final del periodo de observación fue arrestado en algún momento; 0 si no fue arrestado hasta que acabó el estudio.
fin: Categoría lógica que indica si el individuo recibió ayuda financiera al salir de prisión (la ayuda fue dada por el investigador aleatoriamente).
age: en años al momento de dejar la cárcel.
race: Categoría lógica que indica si el individuo es negro o no.
wexp: Categoría lógica, que indica si el individuo tenia trabajo a tiempo completo antes de ser encarcelado.
mar: Categoría lógica, que indica si el individuo estaba casado al momento de salir de la cárcel.
paro: Categoría lógica que indica si el individuo salió en libertad bajo palabra.
prio: Número de encarlamientos previos.
educ: Nivel educativo: 2 (grado 6 o menos), 3 (6 al 9), 4 (10 y 11), 5 (12), o 6 ( post escuela).
emp1 – emp52: El estudio duro 52 semanas. Cada una de estas columnas indica si el individuo obtuvo trabajo en alguna de esas semanas.

Subsetiemos la data y veamos algunas filas:

eha=eha[,c(1:10)]
head(eha,10)

##    week arrest fin age  race wexp         mar paro prio educ
## 1    20      1  no  27 black   no not married  yes    3    3
## 2    17      1  no  18 black   no not married  yes    8    4
## 3    25      1  no  19 other  yes not married  yes   13    3
## 4    52      0 yes  23 black  yes     married  yes    1    5
## 5    52      0  no  19 other  yes not married  yes    3    3
## 6    52      0  no  24 black  yes not married   no    2    4
## 7    23      1  no  25 black  yes     married  yes    0    4
## 8    52      0 yes  21 black  yes not married  yes    4    3
## 9    52      0  no  22 black   no not married   no    6    3
## 10   52      0  no  20 black  yes not married   no    0    5

Las filas 1,2,3,y 7 tiene a un individuo que fue arrestado, las otras filas tiene individuos que no fueron arrestados nuevamente, mientras el estudio se llevó a cabo.

Veamos las estructura de la data:

str(eha)

## 'data.frame':    432 obs. of  10 variables:
##  $ week  : int  20 17 25 52 52 52 23 52 52 52 ...
##  $ arrest: int  1 1 1 0 0 0 1 0 0 0 ...
##  $ fin   : Factor w/ 2 levels "no","yes": 1 1 1 2 1 1 1 2 1 1 ...
##  $ age   : int  27 18 19 23 19 24 25 21 22 20 ...
##  $ race  : Factor w/ 2 levels "black","other": 1 1 2 1 2 1 1 1 1 1 ...
##  $ wexp  : Factor w/ 2 levels "no","yes": 1 1 2 2 2 2 2 2 1 2 ...
##  $ mar   : Factor w/ 2 levels "married","not married": 2 2 2 1 2 2 1 2 2 2 ...
##  $ paro  : Factor w/ 2 levels "no","yes": 2 2 2 2 2 1 2 2 1 1 ...
##  $ prio  : int  3 8 13 1 3 2 0 4 6 0 ...
##  $ educ  : int  3 4 3 5 3 4 4 3 3 5 ...

Reformateando:

eha$educ=as.ordered(eha$educ)

eha$fin=as.numeric(eha$fin)

Grafiquemos la duración:

hist(eha$week)

Uno usa valores discreto para medir duración, pero la distribución de la duración se diferencia de los modelos de conteos (Poisson, etc.). Nótese que el histograma hace pensar que la mayoría de los reos vuelve a la cárcel al final, pero no es cierto, pues esto significa que muchos siguen libres cuando acaba el estudio. Si modelásemos sólo a los que vuelven, perderíamos mucha filas:

table(eha$arrest)

## 
##   0   1 
## 318 114

¿Por qué hacer el esfuerzo de estudiar tantos individuos, y perder tantos en el análisis?

Por eso, la técnica EHA, también conocida en otras profesiones como análisis de supervivencia, usa todos esos datos, denominando a los que no han experimentado el evento de interés (aquí: ser arrestado) como valores censurados.

El primer paso para usar EHA, es indicarle a R que trate a la data de esa manera:

Creación del objeto survival:

library(survival)
eha$survival=with(eha,Surv(week, arrest))
# que es:
eha$survival

##   [1] 20  17  25  52+ 52+ 52+ 23  52+ 52+ 52+ 52+ 52+ 37  52+ 25  46  28 
##  [18] 52+ 52+ 52+ 52+ 52+ 24  52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+
##  [35] 52+ 50  52+ 52+ 52+ 52+ 52+ 52+ 10  52+ 52+ 52+ 52+ 20  52+ 52+ 52+
##  [52] 52+ 52+ 50  52+ 52+ 52+ 52+ 52+ 52+  6  52+ 52+ 52+ 52  52+ 52+ 52+
##  [69] 49  52+ 52+ 52+ 52+ 52+ 52+ 52+ 43  52+ 52+  5  27  52+ 52+ 52+ 22 
##  [86] 52+ 52+ 18  52+ 52+ 52+ 52+ 52+ 52+ 52+ 24  52+ 52+ 52+ 52+  2  26 
## [103] 52+ 49  52+ 21  48  52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+  8 
## [120] 52+ 52+ 49  52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+  8  52+ 52+ 13  52+ 52+
## [137] 52+ 52+  8  52+ 33  52+ 52+ 11  52+ 52+ 52+ 52+ 52+ 37  52+ 52+ 52+
## [154] 44  52+ 52  52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+
## [171] 52+ 52+  9  17  52+ 52+ 52+ 52+ 52+ 52+ 16  52+ 52+  3  52+ 52+ 52+
## [188] 52+ 52+ 52+ 52+ 52+ 52+ 52+ 45  52+ 52+ 52+ 52+ 52+ 52+ 28  52+ 16 
## [205] 15  52+ 52+ 52+ 52+ 14  52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+
## [222] 52+ 52+ 52+ 52+  7  52+ 52+ 43  46  40  52+ 14  52+ 52+  8  52+ 52+
## [239] 52+ 52+ 52+ 25  52+ 52+ 17  37  52+ 52+ 52+ 32  52+ 52+ 52+ 52+ 52+
## [256] 52+ 52+ 52+ 12  52+ 18  52+ 52+ 14  52+ 52+ 52+ 52+ 38  52+ 24  20 
## [273] 32  52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 31  20  40  52+ 52+
## [290] 52+ 52+ 42  52+ 26  52+ 52+ 52+ 52+ 52+ 47  52+ 52+ 40  52+ 52+ 52+
## [307] 52+ 21  52+ 52+ 52+ 52+ 52+  1  43  24  11  52+ 52+ 52+ 52+ 33  52+
## [324] 46  36  52+ 52+ 18  52+ 52+ 50  52+ 34  52+ 35  52+ 52+ 39   9  52+
## [341] 52+ 52+ 34  52+ 52+ 52+ 44  52+ 52+ 35  30  39  52+ 52+ 52+ 52+ 19 
## [358] 52+ 43  52+ 48  37  20  52+ 52+ 36  52+ 52+ 52+ 52+ 52+ 52+ 52+ 30 
## [375] 52+ 52+ 52+ 52+ 52+ 52+ 52+ 52+ 42  52+ 52+ 52+ 52+ 26  40  52+ 52+
## [392] 52+ 35  52+ 46  52+ 49  52+ 52+ 49  52+ 52+ 52+ 52+ 52+ 35  52+ 52+
## [409] 52+ 52+ 27  52+ 52+ 52+ 52  45   4  52  36  52+ 52+  8  15  52+ 19 
## [426] 52+ 12  52+ 52+ 52+ 52+ 52+

Como ves, la columna creada tiene valores con un +, lo que indica que están censurados.

Análisis Kaplan-Meier (KM)

KM es el procedimiento descriptivo básico que se utiliza para ver la dinámica de sobrevivencia. Veamos el comportamiento genérico de permanecer libre:

library(ggplot2)

## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang

library(ggfortify)
KM.generico = survfit(survival ~ 1, data = eha)

###
ejeX='SEMANAS HASTA SER RE-ARRESTADO\n curva cae cuando alguien es arrestado'
ejeY='Prob(PERMANECER LIBRE)'
titulo="Curva de Sobrevivencia (permanecer libre)"
autoplot(KM.generico,xlab=ejeX,ylab=ejeY, main = titulo,conf.int = F)

La gráfica anterior nos da una idea de como se comporta esta población; pero sería interesante ver una comparación:

KM.fondos = survfit(survival ~ fin, data = eha)

###
ejeX='SEMANAS HASTA SER RE-ARRESTADO\n curva cae cuando alguien es arrestado'
ejeY='Prob(PERMANECER LIBRE)'
titulo="Curva de Sobrevivencia (permanecer libre)"
autoplot(KM.fondos,xlab=ejeX,ylab=ejeY, main = titulo,conf.int = F)

La posición de las curvas nos hace pensar que le cuesta más a los que no tuvieron fondos permanecer libres. Para evitar afirmar tal cosa, podemos hacer la prueba de Mantel-Cox (LogRanK):

survdiff(survival ~ fin, data = eha)

## Call:
## survdiff(formula = survival ~ fin, data = eha)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## fin=1 216       66     55.6      1.96      3.84
## fin=2 216       48     58.4      1.86      3.84
## 
##  Chisq= 3.8  on 1 degrees of freedom, p= 0.05

Siendo H0: no hay diferencias entre grupos, con el p-valor obtenido dudamos afirmar que hay diferencias. Este gráfico aclara por qué:

autoplot(KM.fondos,xlab=ejeX,ylab=ejeY, main = titulo,conf.int = T)

Regresión Cox

Como toda regresión, esta técnica permite utilizar regresores, pero no modela la duración sino el riego de que el evento suceda (ser re arrestado):

modelo.REARRESTAR <- coxph(survival ~ fin + race + wexp + mar + paro + prio,
                           data=eha)
summary(modelo.REARRESTAR)

## Call:
## coxph(formula = survival ~ fin + race + wexp + mar + paro + prio, 
##     data = eha)
## 
##   n= 432, number of events= 114 
## 
##                    coef exp(coef) se(coef)      z Pr(>|z|)   
## fin            -0.42052   0.65671  0.19101 -2.202  0.02770 * 
## raceother      -0.27595   0.75885  0.30742 -0.898  0.36938   
## wexpyes        -0.34808   0.70604  0.20331 -1.712  0.08688 . 
## marnot married  0.54873   1.73106  0.37953  1.446  0.14823   
## paroyes        -0.01977   0.98042  0.19406 -0.102  0.91886   
## prio            0.08880   1.09286  0.02832  3.136  0.00171 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## fin               0.6567     1.5228    0.4516    0.9549
## raceother         0.7589     1.3178    0.4154    1.3862
## wexpyes           0.7060     1.4163    0.4740    1.0517
## marnot married    1.7311     0.5777    0.8227    3.6422
## paroyes           0.9804     1.0200    0.6702    1.4342
## prio              1.0929     0.9150    1.0339    1.1552
## 
## Concordance= 0.634  (se = 0.027 )
## Likelihood ratio test= 25.28  on 6 df,   p=3e-04
## Wald test            = 26.98  on 6 df,   p=1e-04
## Score (logrank) test = 27.8  on 6 df,   p=1e-04

Podemos ver que dar financiamiento es significativo, al igual que condenas previas. Una tiene una relación inversa con el riesgo de ser re arrestado (reducen el riesgo), pero la otra sí es directa. Podríamos leer el segundo bloque con los coeficientes exponenciados para dar una interpretación. Esto es diferente si se trata de una variable categórica o una numérica. Así, el riesgo de ser arrestado de alguien con fondos es 33% (1-0.656) menor de alguien sin tener fondos .

ci <- confint(modelo.REARRESTAR)
exp(cbind(coef(modelo.REARRESTAR), ci))["fin",]

##               2.5 %    97.5 % 
## 0.6567057 0.4516261 0.9549102

Por otro lado, cada condena previa aumenta el riesgo de ser re arrestado en 9.29% (1.0929-1).

ci <- confint(modelo.REARRESTAR)
exp(cbind(coef(modelo.REARRESTAR), ci))["prio",]

##             2.5 %   97.5 % 
## 1.092860 1.033860 1.155227

Curso: Estadística Para el Análisis Político II
Semestre 2019-I

Dr. Jose Manuel MAGALLANES, PhD

Departamento de Ciencias Sociales, Sección Ciencia Politica y Gobierno

Análisis de Eventos Históricos (EHA)

Análisis Kaplan-Meier (KM)

Regresión Cox

Curso: Estadística Para el Análisis Político II Semestre 2019-I

Dr. Jose Manuel MAGALLANES, PhD

Departamento de Ciencias Sociales, Sección Ciencia Politica y Gobierno

Análisis de Eventos Históricos (EHA)

Análisis Kaplan-Meier (KM)

Regresión Cox

Curso: Estadística Para el Análisis Político II
Semestre 2019-I