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