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 siguen 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.

Para esto intentaremos ver que variables provocan una diferencia en el tiempo de servicio de los empleados, esto de manera descriptiva. Despúes veremos si esto provoca diferencias lo suficientemente significativas mediante el modelo de riesgos proporcionales de Cox, si es que se cumplen las hipótesis. De no ser posible esto veremos que resultados descriptivos pueden ofrecer alguna estrategia de retención.

Realizamos un análisis exploratorio de la base.

Aquí se muestra un breve resumen del tipo de variables qeu conforman la base de datos que se nis proporcionó. La cual esta compuesta en su mayoría por datos continuos. Observamos que el porcentaje de observaciones con datos faltantes es relativamente muy bajo, del 4.6% para ser precisos. Lo cual es bueno para comenzar a trabajar. Veamos el porcentaje de valores faltantes por variable.

No hay muchos datos faltantes,salvo en las variables de edad(Age), falla disciplinaria (Disciplinary Failure) y Tiempo de ausencia en horas (Absenteeism time in hours) que fue el del procentaje más alto pero de igual manera, sigue siendo relativamente bajo.Una buena señal.

Este es el gráfico de barras para las variables categóricas evidentemente. Podemos notar lo siguiente:

1.Hay pocos empleados que han presentado un desgaste significativo (Attrition), lo cual es un indicador a simple vista de que las condiciones laborales en general son buenas en esta empresa.

2.La empresa se esfuerza en la inclusión laboral de sus empleados. Como prueba de ello se ve que la mayor parte viaja raramente.

3.Predomina el género masculino en los empleados mas no de manera abrupta, hay 20% más empleados que empleadas aproximadamente.

4.Los puestos con mayor número de empleados científicos investigadores, ejecutivos de ventas y técnicos laboratoristas. Los puestos con menor número de empleados son los representantes de ventas, directores de investigación y recursos humanos (nuestro departamento “favorito”).

5.La mayoría de los empleados están casados.

6.La mayoría de los empleados no trabaja horas extra. Lo que puede contribuir mucho a un ambiente laboral sano.

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

Estos son los histogramas de algunas variables. Observamos lo siguiente:

  1. La mayoria de los empleados claifica el ambiente laboral de la empresa de una buena manera.

  2. Como habíamos supuesto anteriormente, la perpeción de la inclusión laboral es en su mayoria buena.

  3. La satisfacción laboral es vista como buena o muy buena en su mayoría.

  4. El ingreso mensual de los trabajadores, empíricamente sigue una distribución sesgada hacia la derecha, es decir, la mayor parte de los empleados tiene un ingreso bajo y pocos con un salario alto.

  5. De igual manera, la distancia de la empresa a los hogares de los trabajadores es empíricamente sesgada a la derecha.

  6. La mayoría de los empleados ha obtenido una calificación suficiente en su última evaluacion.

  7. La media de edad de los trabajadores es de 36 años. La empresa tiene empleados con desde 27 años de edad hasta 58.

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).

Esta es es la “curva” de supervivencia de los empleados en general, obtenida a través de la estimación de Kaplan-Meier que como podemos ver, gracias a los datos censurados no podemos obtener una vista completa del comportamiento. Observamos que la gráfica tiene una pendiente negativa no tan brusca,, es decir la probabilidad de salir de la empresa decae a una tasa muy baja con el paso del tiempo. Si tomamos como ejemplo un empleado de 27 años de edad (el mínimo de todas las edades) tiene probabilidad de poco menos del 50% permanecer en la empresa, aún a sus 52 años de edad. Esto refleja que en general la empresa SÍ hace un esfuerzo por mantener a sus empleados. Antes de continuar ubiquemos la mediana (la linea punteada), que se encuentra a los 24 años de servicio aproximadamente.

Ahora veamos como se comporta el tiempo de permanencia de acuerdo a algunas variables. Decidimos omitir a los indicadores de encuesta externa ya que no disponemos de información que los explique.

Ahora queremos observar lso individuos en riesgo mientras pasa el tiempo.

Vemos que la función de riesgo de salida tiene un comportamiento exponencial, que empieza a acentuarse a partir de aproximadamente de los 23 años de tiempo de servicio, donde el número de individuos en riesgo cae casi a la mitad. Obtuvimos la función de riesgo acumulada vía estimación de Nelson Aalen que a grandes rasgos es ir sumando el cociente de los empleados que salieron entre aquellos en riesgo por cada año de servicio. Se optó por este método ya que no necesita un supuesto de distribución paramétrica, cosa que el método de máxima verosimilitud sí requiere.

Ahora nos interesa ver la curva de permanencia por variables. Omitimos aquellas que corresponden a indicadores de encuesta externa. Pues no disponemos de información de estos para su interpretación.

En cuanto a estas gráficas observamos lo siguiente:

  1. No hay una diferencia en el tiempo de permanencia en las siguientes variables: Satisfacción laboral, Inclusión laboral, Ambiente laboral, Horas extra, Ingreso mensual,Última evaluación,Fallas disciplinarias y número de pasos en la compañia. Lo anterior se constata gráficamente u vía el p-value arrojado, que corresponde al de la prueba Log Rank, cuya hipótesis nula es suponer un comportamiento similar entre dos o más supervivencias. Las variables anteriores obtuvieron los p-values más altos.

  2. Hay diferencias notables en las supervivencias de las variables: Puesto de trabajo y Disposición de mudarse. La primera se constata via gráficamente y la segunda mediante el p-value reportado que es el más bajo.

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.En consecuencia no podemos hacer mucho al respecto dado que los datos proporcionados no cumplen con los supuestos del modelo, por lo que dicho ajuste no será nada bueno.

Por lo que intentamos la siguiente metodología (esto con el fin de comprobar que no será un buen ajuste pero se obtuvieron resultados interesantes):

  1. Tomar en cuenta a las covariables con los p-values más pequeños, fijamos como límite máximo \(p=0.15\).

  2. Con la tabla de frecuencias de individuos en riesgo , revisar por cada 5 años la tendencia por estrato, y si en algún punto el orden cambia, descartar a la variable.

Después de verificar estas condiciones, de un total de 19 variables analizadas únicamente 5 los cumplieron .Probamos con distintos modelos de distintas variables y no obtuvimos resultados diferentes, en total 8.

Para la elección del mejor modelo, calculamos los coeficientes de información de Akaike. Este técnica estadística en resumidas cuentas esta basada en estimar la verosimilitud de un modelo para predecir valores futuros. Su expresión matemática es la siguiente:

\[AIC= 2k-2log(L)\] Donde \(k\) es el número de parámetros a estimar, que en nuestro caso corresponde a las entradas del vector \(\hat{\theta}=(\theta_1,\theta_2,...,\theta_n)\), y \(L\) es el máximo de la función de verosimilitud.

Existe también el criterio de información Bayesiano. Cuya técnica se basa en medir el sacrificio de ajuste por complejidad del modelo.Su expresión matemática es la siguiente: \[BIC=klog(n)-2log(L)\] Donde \(k\) es el número de parámetros a estimar, \(n\) el tamaño de muestra, y \(L\) es el máximo de la función de verosimilitud.

Confiamos más en el coeficiente Bayesiano. Pues en este caso la penalización por un sobreajuste es mayor, 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.Es importante recordar que tanto el criterio de información de Akaike como el Bayesiano ofrecenn una medida de calidad relativa, no absoluta. Estos fueron los coeficientes del modelo:

## [1] 22184.58
## [1] 22231.62

Veamos un resumen del modelo.

## Call:
## coxph(formula = Surv(data2.length_of_service, Attrition) ~ JobRole + 
##     Will_Relocate, data = data_p6)
## 
##   n= 8000, number of events= 1376 
## 
##                                   coef exp(coef) se(coef)      z Pr(>|z|)
## JobRoleHuman Resources        -0.22451   0.79891  0.15615 -1.438   0.1505
## JobRoleLaboratory Technician  -0.25472   0.77513  0.10786 -2.361   0.0182
## JobRoleManager                -0.34827   0.70591  0.14021 -2.484   0.0130
## JobRoleManufacturing Director -0.04120   0.95964  0.11885 -0.347   0.7288
## JobRoleResearch Director      -0.19256   0.82484  0.14285 -1.348   0.1777
## JobRoleResearch Scientist     -0.22958   0.79487  0.10444 -2.198   0.0279
## JobRoleSales Executive        -0.23373   0.79158  0.10532 -2.219   0.0265
## JobRoleSales Representative   -0.19736   0.82089  0.14173 -1.393   0.1637
## Will_Relocate                 -0.11643   0.89009  0.05424 -2.147   0.0318
##                                
## JobRoleHuman Resources         
## JobRoleLaboratory Technician  *
## JobRoleManager                *
## JobRoleManufacturing Director  
## JobRoleResearch Director       
## JobRoleResearch Scientist     *
## JobRoleSales Executive        *
## JobRoleSales Representative    
## Will_Relocate                 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                               exp(coef) exp(-coef) lower .95 upper .95
## JobRoleHuman Resources           0.7989      1.252    0.5883    1.0850
## JobRoleLaboratory Technician     0.7751      1.290    0.6274    0.9576
## JobRoleManager                   0.7059      1.417    0.5363    0.9292
## JobRoleManufacturing Director    0.9596      1.042    0.7602    1.2114
## JobRoleResearch Director         0.8248      1.212    0.6234    1.0914
## JobRoleResearch Scientist        0.7949      1.258    0.6477    0.9754
## JobRoleSales Executive           0.7916      1.263    0.6439    0.9731
## JobRoleSales Representative      0.8209      1.218    0.6218    1.0837
## Will_Relocate                    0.8901      1.123    0.8003    0.9899
## 
## Concordance= 0.529  (se = 0.009 )
## Rsquare= 0.002   (max possible= 0.938 )
## Likelihood ratio test= 16.59  on 9 df,   p=0.06
## Wald test            = 16.96  on 9 df,   p=0.05
## Score (logrank) test = 17.02  on 9 df,   p=0.05

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.

Pues ya hemos visto, podríamos decir hasta el cansancio que no es posible implementar un modelo de riesgos proporcionales pues las hipotésis de riesgos proporcionales no se cumplen. A pesar de verificarlo gráficamente, este es un primer filtro. Por lo que si este no se cumplió, un criterio más formal habría fallado también.

Sin embargo el intento de la implementeción arrojó un resultado importante. Si revisamos el resumen, al parecer si el trabajador esta en el puesto de Administrador (Manager), director de investigación o ejecutivo en ventasm tendrá cierta relevancia (significancia) en el tiempo de permanencia. De igual manera si el empleado esta dispuesto a mudarse (Will Relocate). Además graficamente si seguimos la línea punteada en las supervivencias de las variables mencionadas (la mediana), vemos que hay una diferencia importante en el caso de los puestos de trabajo, donde la mediana se alcanza antes en los puestos que no resulltaron significativos ejemplo director de manufactura (Manufacturing Director). Por lo que con esto llegamos a los siguientes resultados.

Conclusiones

  1. En cuanto a puestos de trabajo, aquellos que reportan una salida más tardía son los técnicos laboratoristas, los administradores, los cientificios de investigación y los ejecutivos de ventas. Vale la pena analizar más a fondo sus condiciones laborales (salario, prestaciones, ambiente laboral, viajes de negocios, inclusión laboral, etc.) y pueden tomarse como referencia para mejorar las codiciones de aquellos que salen más temprano como los de recursos humanos, representantes de salud y directores de investigación.

  2. Aquellos que están dispuestos a cambiar de residencia, suelen salir más tarde de la empresa que aquellos que no lo hacen. Habría que analizar ¿Qué puestos ocupan? ¿Cual es su salario? ¿Que tan bueno es su desmpeño?,¿Cómo califican su satisfacción inclusión y ambiente laboral? Solicitamos a la compañía llevar a cabo este estudio, identificando a aquellos que sí están dispuestos a mudarse vía su ID.

  3. Como vimos en un principio, la permanencia de los empleados es buena, por lo que cabe mencionar que aún considerando las conclusiones anteriores dicho tiempo no mejorará sustancialmente (reiteramos, porque ya de por si es bueno).En otras palabras,la empresa puede mejorar en los rubros que considere necesarios como por ejemplo aumentar sueldos, aumentar prestaciones laborales, mejorar el ambiente de trabajo, involucrar más a los empleados, etcétera. Más no es probable, y por ende no es garantía de que las decisiones tomadas al respecto, provoquen una mejoría sustancial en el tiempo en que los empleados se queden en la empresa. Sin embargo esto también puede ser una señal que las cosas dentro de la compañía se hacen bien, tanto así que factores que parecen influir en la carrera profesional de los empleados en la compañía como el salario, puesto de trabajo, ambiente laboral, satisfacción laboral, entre otros; no parecen tener suficiente impacto en el tiempo que laboran dentro de la empresa. En pocas palabras la empresa va por buen camino, cuida de sus empleados.

Somos conscientes, mejor que nadie que resultados gráficos no pueden ser usados como una prueba formal para una hipótesis en el comportamiento de los riesgos, mas proveen una fuerte evidencia de tales hipótesis. Hubieramos preferido dar una respuesta basada en un modelo estadístico, sin embargo aunque no ofrecimos la respuesta que buscabamos pudimos ofrecer un resultado vía análisis descriptivo en su mayoria.