Introducción

Modelo de Riesgos Proporcionales de Cox

Dentro del análisis de supervivencia, es frecuente querer saber, dadas qué condiciones, el individuo alcanzará máss rápida o lentamente la falla; que puede ser cualquier evento de interés. Depende del estudio que sea realizado, se buscará adelantar o ralentizar el tiempo de fallo. Durante el curso se han visto como estimar la función de supervivencia como Kaplan-Meier, y hemos visto como comparar funciones de supervivencia para determinar si siguenn un comportamiento similar o no, Log-Rank en el caso no paramétrico y F de Cox para el caso paramétrico. Los métodos mencionados anteriormente son útiles pero únicamente en el caso univariado, ignorando el impacto con otras variables. En adición Kaplan-Meier y Log-Rank son útiles cuando la variable predictora es categórica y no funcionan tan fácilmente cuando es cuantitativa (número natural,entero o real).

El Modelo de Riesgos Proporcionales de Cox nos permite establecer una relación de proporcionalidad directa o inversa entre una variable respuesta y un conjunto de variables explicativas que pueden ser tanto cuantitativas como categóricas.

Objetivo

El propósito de ajustar este modelo es evaluar el impacto de un conjunto de variables sobre la supervivencia en cuestión. Así la función de riesgo explicada sigue la siguiente expresión.

\[h(t)=h_0(t)e^{(\theta_{1}x_{1}+\theta_{2}x_{2}+...+\theta_{n}x_{n})}\] Donde: \(t\) es el tiempo. \(h(t)\) es la función de riesgo explicada por el modelo. \(h_0(t)\) se le conoce como riesgo base. \(\vec{\theta}\) es el vector de coeficientes. \(\vec{x}\) es el vector de covariables.

La idea de este modelo es a partir de un modelo análogo a una regresión lineal múltiple (o simple en el caso univariado), establecer como cambia el riesgo original (base) en función de las covariables.

Planteamiento del problema

Para este trabajo se nos dio una base de datos que contiene distintas características de los trabajadores de una compañía, con el objetivo de determinar que mejoras puede hacer esta compañía para evitar la salida temprana de sus trabajadores. Realizaremos un análisis exploratorio de la base.

## Loading required package: dplyr
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: ggplot2
## Loading required package: survival
## Loading required package: survminer
## Loading required package: ggpubr
## Loading required package: magrittr
## Loading required package: flexsurv
## Loading required package: DataExplorer

Veamos el porcentaje de valores faltantes No hay muchos datos faltantes, una buena señal. Este es el gráfico de barras para las variables categóricas evidentemente. 

## Warning: Computation failed in `stat_bin()`:
## `binwidth` must be positive

Estos son los histogramas de las variables contínuas.

Funciones de supervivencia

Ahora porcedemos a un análisis visual de las variables que podriamos usar en el modelo. Para esto graficaremos distintas supervivencias por variable, donde el tiempo será el tiempo en años de servicio y la falla si el empleado ha experimentado algún detrimento en sus capacidades que le impidan trabajar (Attrition).

Pero podemos notar que en una mayoria de variables, las supervivencias se cruzan por lo que el supuesto gráfico de riesgos proporcionales no se cumple. Además formalmente el p-value de la hipotésis nula de que ambas supervivencias son iguales es muy alto. Por lo que no podemmos hacer mucho al respecto dado que los datos mostrados no cumplen con los supuestos del modelo, por lo que dicho ajuste no será nada bueno. Por ejemlo ajustemos el modelo de Cox con las covariables de puesto de trabajo e inmersión laboral. Obtenemos el siguiente resumen

## Call:
## coxph(formula = Surv(data2.length_of_service, Attrition) ~ JobRole + 
##     JobInvolvement, data = data_p6)
## 
##   n= 8000, number of events= 1376 
## 
##                                   coef exp(coef) se(coef)      z Pr(>|z|)
## JobRoleHuman Resources        -0.21255   0.80852  0.15612 -1.361   0.1734
## JobRoleLaboratory Technician  -0.25174   0.77745  0.10787 -2.334   0.0196
## JobRoleManager                -0.34736   0.70655  0.14020 -2.478   0.0132
## JobRoleManufacturing Director -0.03465   0.96594  0.11886 -0.292   0.7706
## JobRoleResearch Director      -0.18948   0.82739  0.14285 -1.326   0.1847
## JobRoleResearch Scientist     -0.23043   0.79419  0.10445 -2.206   0.0274
## JobRoleSales Executive        -0.23704   0.78896  0.10530 -2.251   0.0244
## JobRoleSales Representative   -0.19011   0.82686  0.14168 -1.342   0.1796
## JobInvolvement                -0.04002   0.96077  0.03777 -1.059   0.2894
##                                
## JobRoleHuman Resources         
## JobRoleLaboratory Technician  *
## JobRoleManager                *
## JobRoleManufacturing Director  
## JobRoleResearch Director       
## JobRoleResearch Scientist     *
## JobRoleSales Executive        *
## JobRoleSales Representative    
## JobInvolvement                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                               exp(coef) exp(-coef) lower .95 upper .95
## JobRoleHuman Resources           0.8085      1.237    0.5954    1.0979
## JobRoleLaboratory Technician     0.7774      1.286    0.6293    0.9605
## JobRoleManager                   0.7065      1.415    0.5368    0.9300
## JobRoleManufacturing Director    0.9659      1.035    0.7652    1.2193
## JobRoleResearch Director         0.8274      1.209    0.6253    1.0947
## JobRoleResearch Scientist        0.7942      1.259    0.6472    0.9746
## JobRoleSales Executive           0.7890      1.267    0.6418    0.9698
## JobRoleSales Representative      0.8269      1.209    0.6264    1.0915
## JobInvolvement                   0.9608      1.041    0.8922    1.0346
## 
## Concordance= 0.527  (se = 0.009 )
## Rsquare= 0.002   (max possible= 0.938 )
## Likelihood ratio test= 13.08  on 9 df,   p=0.2
## Wald test            = 13.47  on 9 df,   p=0.1
## Score (logrank) test = 13.53  on 9 df,   p=0.1

Vemos que casi todas las covariables no son significativas. Intentamos tomar en cuenta a las covariables con los p-values más pequeños. Probamos con distintos modelos de distintas variables y no obtuvimos resultados diferentes, en total 8. Dichas combinaciones de covariables fueron: 1. Puesto de trabajo con Inmersión laboral. 2. Puesto de trabajo con Emp_Position 3. Puesto de trabajo con Viajes de Negocios. 4. Puesto de trabajo con si el empleado esta dispuesto a mudarse. 5. Puesto de trabajo con si el empleado esta dispuesto a mudarse con Emp_Sat_Remote. 6. Puesto de trabajo con si el empleado esta dispuesto a mudarse con pasos en el trabajo. 7. Puesto de trabajo con si el empleado esta dispuesto a mudarse con Promedio de presión arterial mínima. 8. Puesto de trabajo con si el empleado esta dispuesto a mudarse con Promedio de presión arterial mínima con Viajes de Negocios.

Para la elección del mejor modelo, calculamos los coeficientes de Akaike.

## [1] 22188.09
## [1] 22188.24
## [1] 22187.27
## [1] 22184.58
## [1] 22184.64
## [1] 22186.58
## [1] 22183.92
## [1] 22184.08

Ahora calculamos los coeficientes Bayesianos

## [1] 22235.13
## [1] 22235.29
## [1] 22239.54
## [1] 22231.62
## [1] 22236.91
## [1] 22238.85
## [1] 22236.19
## [1] 22246.8

Confiamos más en el coeficiente Bayesiano y el criterio es elegir al de menor valor, por lo que el modelo 4 es el indicado, es decir aquel que relaciona puesto de trabajo y si el empleado está dispuesto a mudarse.

Ahora realizaremos una gráfica de diagnóstico del modelo que se presenta a continuación.

Esta gráfica indica la distribución de los residuales, puesto que el modelo de riesgos proporcionales de Cox involucra una regresión, los residuales deben poder ser ajustados via una regresión lineal es decir linealidad de los residuos, cosa que evidentemente no se cumple pese a elegir el “mejor” modelo.

Conclusiones

Desafortunadamente desde que graficamos las supervivencias de distintas covariables nos dimos cuenta que no hay una relación proporcional entre los riesgos de cada atributo, lo vimos graficamente y a la hora de experimentar con el ajuste de distintos modelos no obtuvimos una respuesta favorable.