Regresión con variable dependiente cualitativa

REGRESIÓN CON VARIABLE DEPENDIENTE CUALITATIVA, POR JOSÉ MANUEL ROJO ABUÍN I. INTRODUCCIÓN En muchas ocasiones estaremos interesados en predecir los valores de una variable dicotómica binaria, es decir, una variable que sólo puede tomar dos valores, los valores son complementarios y dichos valores no son comparables, como sucede en regresión lineal.

Ejemplos de variable dependiente dicotómica pueden ser: sano o enfermo, paga o no paga, …, etc.

El modelo de regresión logística se utiliza cuando estamos interesados en pronosticar la probabilidad de que ocurra o no un suceso determinado. Por ejemplo, a la vista de un conjunto de pruebas médicas, que una persona tenga una determinada enfermedad, o bien que un cliente devuelva un crédito bancario.

A diferencia del análisis discriminante que requiere la normalidad multivariante de los datos, el análisis de regresión logística sólo precisa del principio de monotonía, es decir, si el suceso A es que una determinada persona padezca de artrosis y X representa la edad, deberá de ocurrir:

xi ? xj ? P(A/xi) ? P(A/xj)

A diferencia del análisis discriminante, podremos estudiar el impacto que tiene cada una de las variables explicativas en la probabilidad de que ocurra el suceso en estudio.

El análisis de regresión logística es una herramienta muy flexible en cuanto a la naturaleza de las variables explicativas, pues éstas pueden ser de escala y categóricas. II. PLANTEAMIENTO DEL PROBLEMA Supongamos que tenemos la variable de estudio codificada de la siguiente manera:

?0 Noocurreelsuceso ?1 Si ocurreelsuceso

Además, vamos a considerar que sólo tenemos una variable explicativa X ; en estas condiciones podríamos considerar un modelo de regresión lineal con el propósito de ver qué dificultades nos van a surgir:

yi ? pi ? a?b*xi ?ui

Si estimamos este modelo y representamos gráficamente la recta de regresión:

Podemos observar que la línea de regresión no está acotada en el intervalo [0,1] y, por lo tanto, ya no va a representar una probabilidad.

Además, consideraciones de índole matemática nos llevan a la conclusión de que los residuos no van a ser homocedásticos y, por tanto, la técnica de estimación por mínimos cuadrados dejará de ser un método óptimo de estimación.

Una forma que tenemos de garantizar que los valores pronosticados estén en el intervalo [0,1] es considerar la siguiente transformación:

p(a/x) ? F(x*b)

Donde F es una función de distribución.

Nota Una función de distribución es una función real de variable real: F :R ? R

De forma que verifica:

Está acotada en el intervalo [0,1] 0? F(x) ?1 ? x

Es monótona no decreciente: x1 ? x2 ? F(x1) ? F(x2)

Y, además, está definida en todo R, tomando los siguientes valores: F(??) ? 0 F(??) ?1

En general, la grafica de una función de distribución es: Si utilizamos la función de distribución logística, el análisis se denomina Regresión Logística, y si utilizamos la función de distribución normal se denomina Regresión Probit.

III. EL MODELO DE REGRESIÓN LOGÍSTICA El modelo de regresión logística parte de la hipótesis de que los datos siguen el siguiente modelo: p 1? p ln( ) ?b0 ?b 1*x1 ?b2*x2 ?…?bk *xk ?u ? x*b?u Con el fin de simplificar la notación, definimos Z:

z ?b0 ?b 1*x1?b2*x2 ?…?bk *xk

Por lo tanto, el modelo se puede representar como: ) ? z ?u p 1? p ln( Donde p es la probabilidad de que ocurra el suceso de estudio.

Operando algebráicamente sobre el modelo: ) ? z p 1? p ln( ? ez p 1? p p ?(1? p)*ez p ?ez ? p*ez p? p*ez ?ez p(1?ez) ?ez p ? ez 1?ez Como la función de distribución logística es: ex 1?ex F(x) ? Por tanto, podemos reescribir el modelo de forma mucho más compacta: ? F(z) ? F(x*b) ez 1?ez p ? De donde se deduce que el modelo de regresión logística es, en principio, un modelo de regresión no lineal, pero es lineal en escala logarítmica atendiendo a su definición original:

pi logP(Y) ??yi ?Log( )??log(1? pi) ) ? z p 1? p ln( ln(p)?ln(1? p) ? z

ln(p)?ln(1? p) ?b0 ?b 1*x1 ?b2 *x2 ?…?bk *xk

Es decir, la diferencia de la probabilidad de que ocurra un suceso respecto de que no ocurra es lineal pero en escala logarítmica. Por tanto, el significado de los coeficientes, aunque guardando una cierta relación con el modelo de regresión lineal, va a ser algo más complejo de interpretar.

Recordemos las dos formas más importantes de expresar el modelo de regresión logística:

ln(p)?ln(1? p) ?b0 ?b 1*x1 ?b2 *x2 ?…?bk *xk ? eb0 *eb1*X1 *eb2*X 2…ebk*X k p 1? p La primera expresión se llama logit y a la segunda Odds ratio o cociente de probabilidades. IV. ESTIMACIÓN DE LOS PARÁMETROS. i Brevemente, vamos a ver en esquema el problema que ofrece, en el caso de regresión logística, la estimación de los parámetros.

Sea una muestra de n elementos, donde se ha observado la variable respuesta Y (que sólo puede tomar dos valores: cero y uno) y la variable X .

La función de probabilidad de una observación cualquiera es:

P(Y ?1/x) ? p P(Y ? 0/x) ?1? p Por tanto:

P(Y /x) ? py *(1? p)1?y

Por tanto la función de probabilidades de la muestra es: P(y1,y2,…,yn) ??pYi *(1? pi)1?yi i

Esta expresión recibe el nombre de verosimilitud de la muestra (likelihood).

Tomando logaritmos: n n

i i 1? pi Expresando pi en función de los parámetros que deseamos estimar:

n

i

Resulta obvio que aunque derivemos y establezcamos la condición de máximo, no vamos a poder despejar los coeficientes B .

La solución que vamos a obtener es:

?1

? ?B*?B ? ?B ?

Esta solución establece cómo encontrar una solución ( Ba ) a partir de un punto próximo cualquiera, denominado B0 . Por lo tanto, deberemos de hacer una estimación inicial del valor de los verdaderos parámetros y mediante un procedimiento recursivo encontrar el verdadero valor de los mismos. Para encontrar los verdaderos valores se suele utilizar el algoritmo de Newton-Raphson.

Gráficamente: V. EJEMPLO 1 Vamos a ir introduciendo los elementos de esta técnica a través de un sencillo ejemplo.

El tratamiento y pronóstico del cáncer depende de cuánto se haya extendido la enfermedad.

Unas de las zonas propensas a ser afectadas por la enfermedad son los ganglios linfáticos.

Si los ganglios linfáticos están afectados el tratamiento pierde efectividad.

Para ciertos tipos de cáncer es preciso realizar una intervención quirúrgica para determinar si la enfermedad se ha extendido al sistema linfático, y así determinar qué tratamiento se deberá de aplicar.

Si en función a una serie de pruebas médicas no invasivas se pudiera determinar si los ganglios linfáticos están afectados o no se ahorraría tiempo y molestias a los pacientes.

Los datos que vamos a analizar pertenecen a una muestra aleatoria de 53 pacientes masculinos con cáncer de próstata. A cada paciente se le han medido las siguientes variables o características: ? ? ? ? ? ? Xray: Resultado de la prueba de rayos X Grado: Grado de agresividad del tumor. Estado: Cómo está de extendida la enfermedad. Nodos: Indicador de si los ganglios linfáticos están afectados o no por la enfermedad. Edad: edad del paciente. Acido: Prueba de laboratorio del nivel de ácido phosphatase. A continuación mostramos los estadísticos descriptivos de las variables involucradas en el análisis. Es de particular importancia asegurarse que las variables del tipo ausencia/presencia estén codificadas como cero y uno.

xray Prueba de rayos X grado Grado de agresividad

estado Estado de la enfermedad

nodos Estado de los ganglios linfaticos

Statistics a. Multiple modes exist. The smallest value is shown

V.1. Coeficientes estimados del modelo logístico

W(bj) ? ? ? b ˆj ? En la columna segunda se muestran los coeficientes estimados B. Para poder dichos que tener interpretar coeficientes hay en cuenta como están codificadas las variables, pues las dos primeras: edad y acido son contínuas y el resto están codificadas como 0 o 1, para indicar ausencia o presencia de una determinada característica.

En la tercera columna es muestra la desviación típica del estimador.

La cuarta columna muestra el estadístico de Wald; el estadístico de Wald es: ? ? 2 ? ??(bj)?

y dicho estadístico se distribuye de acuerdo con una 2 ?1 ; por tanto, todos los coeficientes que tengan un W(bj) > 4 serán significativos.

La sexta columna (sig) es el p-value del coeficiente.

La séptima columna es el exponencial del coeficiente. El interés del exponencial de los coeficientes es el estudio del impacto de las variables cualitativas.

Ejemplo.

En este ejemplo hemos codificado la prueba de rayos x de forma dicotómica (0,1). Por tanto: ? ez ? k *e2.045*xray p 1? p Si la prueba de rayos X es negativa, la variable vale 0, y si es positiva la variable vale 1; por tanto, si la prueba de rayos x es positiva el cociente de probabilidades aumenta:

e2.045 ? 7.73

Pues: ? ez ? k *e2.045*xray p 1? p Resultado negativo de la prueba: ? ez ? k *e2.045*0 ? k *1 p 1? p Resultado positivo de la prueba: Con las variables anteriores vamos a intentar construir un modelo de regresión logística para tratar de pronosticar en qué pacientes se encuentran los ganglios linfáticos (nodos) afectados por la enfermedad. Coeficientes del modelo de regresión logística. Variables in the Equation a. Variable(s) entered on step 1: edad, acid, xray, grado, estado.

? ez ? k *e2.045*1 ? k *7.73 p 1? p Luego, si la prueba de rayos x es positiva, la probabilidad de tener el sistema linfatico afectado queda multiplicada por 7.73.

V.2. Estimando probabilidades Con los coeficientes estimados ya es posible predecir la probabilidad de que una persona tenga los ganglios linfáticos afectados por el cáncer simplemente construyendo la función de probabilidad: 1 1?e?z P(nodo ?1/x) ? Donde:

Z ? 0.62?1.56*estado ?0.76*grado?2*xray ?0.024*acido ?0.07*edad

A la vista de esta ecuación, podemos estimar la probabilidad de que un hombre con determinadas características tenga el sistema linfático afectado.

Por ejemplo, la probabilidad de que un hombre de 60 años, con un nivel de ácido de 50 y negativo en el resto de las pruebas es de:

Z ? 0.62?1.56*estado ?0.76*grado?2*xray ?0.024*acido ?0.07*edad

Z ? 0.62?1.56*0? 0.76*0? 2*0? 0.024*50?0.07*60

Z ? ?2.38 ? 0.085 1 ?(?2.38) 1?e P(nodo ?1/x) ? En cambio la misma persona dando positivo en todas las pruebas va a tener una probabilidad estimada de:

Z ? 0.62?1.56*1?0.76*1? 2*1?0.024*50?0.07*60

Z ?1.94 ? 0.87 1 1?e?(1.94) P(nodo ?1/x) ? V.3. Interpretando los coeficientes Si bien en regresión lineal la interpretación de los coeficientes de regresión es simple e intuitiva: Bk es el incremento producido en la variable dependiente por un incremento unitario en la Xk . variable En la regresión logística no va a ser tan intuitiva, al depender tanto del valor de Xk donde se produzca el incremento como del valor del resto de las variables, pues la pendiente de la curva de regresión va a ir variando.

Para ayudar a interpretar los coeficientes de regresión logística definimos el Odds Ratio como el cociente de probabilidades entre que ocurra un suceso respecto de que no ocurra: ? OddRatio ? P 1? P P(Y ?1) P(Y ? 0) Teniendo en cuenta que el modelo de regresión logística puede ser escrito como:

ln(p)?ln(1? p) ?b0 ?b 1*x1 ?b2 *x2 ?…?bk *xk ln( ) ? b0 ?b 1*x1 ?b2 *x2 ?…?bk *xk p 1? p Los coeficientes B indican el incremento de la probabilidad de que ocurra el suceso, es decir, la probabilidad de que el sistema linfático esté afectado respecto de que no esté afectado pero en escala logarítmica.

Si el coeficiente p-esimo vale cero, indica que la variable p-esima no afecta a la ocurrencia del suceso.

Si el coeficiente p-esimo es negativo indica que a media que dicha variable va aumentando va a ir disminuyendo el logaritmo del cociente de probabilidades y al revés si es positivo.

Si tomamos exponenciales: ? eB0 *eB1*x1 *eB2*x2 *…*eBk *xk p 1? p Por tanto el coeficiente Bp e va a significar por cuánto se multiplica el Odds Ratio. Otra forma de verlo algo más intuitiva es considerar la derivada de la función de regresión respecto de la p-esima variable.

Tenemos que la probabilidad de ocurrencia del evento es una función de X y B:

P(Y ?1/ X) ? ?(X,B) ? ?(X *B)

Si derivamos respecto de la p-esima variable: ? ??(X *B)bp ? ??(b0 ?b 1×1 ?…bkxk)bp ?? ?xp El problema es que la derivada va a depender de qué valor tomamos para las k variables, es decir, en qué punto vamos a evaluar la curva.

Podemos evaluarla en el punto medio:

??(b0 ?b 1×1 ?…bkxk)bp

O bien podemos considerar la media de los incrementos: compute edadB= compute acidB= compute xrayB= compute gradoB= compute estadoB= 1/(1+exp(-z))**2*exp(-z)*(-0.069). 1/(1+exp(-z))**2*exp(-z)*(0.024). 1/(1+exp(-z))**2*exp(-z)*(2.045). 1/(1+exp(-z))**2*exp(-z)*(0.761). 1/(1+exp(-z))**2*exp(-z)*(1.56). execute. mean edadB to estadoB.

Incremento en una persona media: El significado de la primera expresión es el incremento en la probabilidad de ocurrencia del suceso en una persona media por un incremento unitario en la p-esima variable. La segunda expresión indica cuál es la media de los incrementos de la probabilidad de ocurrencia del suceso por un incremento unitario en la p-esima variable.

Las dos últimas formas no están implementadas en la aplicación SPSS y deberemos realizarlas a mano. En el caso que nos ocupa:

Media de los incrementos:

Report

Código en SPSS

compute z= -0.069*edad+0.024*acid+2.045*xray+0.761*grado+1.564*estado+0.062. execute. Report ? i ??(b0 ?b1xi,1 ?…bkxi,k)bp n

Código en SPSS

compute z= -0.069*59.4+0.024*69.42+2.045*0.28+0.761*0.38+1.564*0.51+0.062. compute edadB= compute acidB= compute xrayB= compute gradoB= compute estadoB= 1/(1+exp(-z))**2*exp(-z)*(-0.069). 1/(1+exp(-z))**2*exp(-z)*(0.024). 1/(1+exp(-z))**2*exp(-z)*(2.045). 1/(1+exp(-z))**2*exp(-z)*(0.761). 1/(1+exp(-z))**2*exp(-z)*(1.56). execute. mean edadB to estadoB. VI. EJEMPLO COMPLETO Seguimos con el ejemplo anterior, pero mostrando tanto estadísticos de bondad de ajuste como los de contraste de regresión. En primer lugar se muestran los esquemas de codificación de las variables, tanto la variable respuesta como las variables categóricas: A la vista del esquema de codificación de la variable respuesta, el modelo va a tratar de predecir la probabilidad de que una persona tenga el sistema linfático afectado. En el resto de las variables categóricas vemos que el esquema de codificación interno coincide con el externo.

VI.1. Historial de las iteracciones

Iteration History a,b,c,d,e a. Method: Forward Stepwise (Conditional) b. Constant is included in the model. c. Initial -2 Log Likelihood: 70,252 d. Estimation terminated at iteration number 4 because parameter estimates changed by less than ,001. e. Estimation terminated at iteration number 5 because parameter estimates changed by less than ,001. Realiza dos pasos, por lo tanto se han introducido dos variables en el modelo. En cada paso va aumentando la verosimilitud del modelo, lo cual implica que disminuye la siguiente expresión: -2 log(verosimilitud) (-2LL). En los dos pasos el algoritmo termina correctamente porque se alcanza el criterio de parada, es decir, el cambio entre los coeficientes estimados en a ultima iteración es inferior a 0.001.

VI.2. Contraste de regresión El contraste de regresión en estos modelos no se realiza sobre la descomposición de la suma de cuadrados como en regresión lineal sino sobre el incremento de la verosimilitud, mas exactamente sobre la disminución de -2LL.

Hipótesis

H0 b 1 ?b2 ?…?bk ?0 H1 ? bp ? 0

Construcción del contraste

C2LL? ?2LL(b0)?(?2LL(b0,b 1,b2,…,bk) j ? 2 , donde J es con sólo la constante es de La verosimilitud 70.252. La diferencia de verosimilitudes se distribuye de acuerdo con una distribución la diferencia del número de parámetros en el modelo. Iteration History a,b,c 1 2 3 Iteration Step 0 -2 Log likelihood 70,253 70,252 70,252 Coefficients Constant -,491 -,501 -,501 parameter estimates changed by less than ,001. a. Constant is included in the model. b. Initial -2 Log Likelihood: 70,252 c. Estimation terminated at iteration number 3 because

P(?1 La verosimilitud (-2LL) con una sola variable es de 59.001.

La verosimilitud (-2LL) con dos variables es de 53.353. La primera variable que entra produce una disminución en -2LL de:

70.252-59.001= 11.251. 2

Por lo tanto rechazamos la hipótesis nula y aceptamos que la primera variable es significativa. sobre la primera produce La segunda variable una reducción de -2LL de:

59.001-53.353=5.647 2 ?5.647) = 0.01747 Por lo tanto la introducción de la segunda variable sigue siendo significativa.

VI.3. Medidas de bondad del ajuste En este tipo de modelos no se emplea el R2 para mostrar la bondad del ajuste, sino que se 2 calcula el incremento de la verosimilitud, aunque reciben el nombre de R no van a tener el L(b0) L(b0,b 1,…bk) R2 ?1? significado geométrico que tienen en regresión lineal por lo tanto deberían de llamarse pseudos R2 .

Model Summary

Cox y Snell: a. Estimation terminated at iteration number 4 because parameter estimates changed by less than ,001. b. Estimation terminated at iteration number 5 because parameter estimates changed by less than ,001.

70.25 53 RMax RMax ?1??L(b0)?N 2 ? L(b0) ?N ? L(b0,b 1,…bk)?

En este ejemplo: 53.353 2 ) ? 0.273 R2 ?1?( Este coeficiente está acotado:

0 ? R2 ?1

Es decir, no puede alcanzar el valor 1.

R cuadrado de Nagelkerke Nagelkerke prefiere definir R2 como: R2 2 R2 ? Donde 2 2 El test de Hosmer y Lemeshow es un constaste de distribución.

La hipótesis nula es que no hay diferencias entre los valores observados y los valores pronosticados (probabilidades); la alternativa es que sí las hay. Por tanto, el rechazo de este test indica que el modelo no está bien ajustado.

En este caso la significatividad de este es de 0.798, no rechazamos la hipótesis nula y por tanto no rechazamos que el modelo tiene falta de ajuste. Para así poder alcanzar el valor 1.

Aunque estos coeficientes tratan de medir la variabilidad explicada, en general, van a ser mucho más bajos que en regresión lineal y deberán de ser complementados con otras medidas de bondad de ajuste.

Test de Hosmer y Lemeshow

Contingency Table for Hosmer and Lemeshow Test 29 4 18 11 3 1 29,000 4,000 18,593 10,407 2,407 1,593 9 11 3 6 2 9 9,000 11,000 2,407 6,593 2,593 8,407 38 15 21 17 5 10 1 2 1 2 3 4 Step 1 Step 2 Observed Expected nodos Estado de los ganglios linfaticos = 0 No afectados Observed Expected nodos Estado de los ganglios linfaticos = 1 Afectados Total

De las 9 + 11 personas que sí tienen los ganglios afectados, 11 han sido pronosticados como afectados, un porcentaje de aciertos del El porcentaje global de aciertos es del

Tabla de clasificación Si bien los coeficientes de bondad de ajuste no son del todo fiables, la tabla de clasificación es normalmente el criterio que debemos de seguir para indicar la bondad de ajuste del modelo.

En esta tabla se muestran los casos bien clasificados en la diagonal principal, y los casos mal clasificados en la segunda diagonal.

Classification Table a. The cut value is ,500

De las 29 + 4 personas que no tienen los ganglios afectados, 29 han sido pronosticados como sanos, es decir, un porcentaje de aciertos del 29 33

11 20 ?87%

? 55% ? 75.5% 29?11 29?11?4?9