Tema 4: Regresión Logı́stica

Introducción
Es uno de los procedimientos estadı́sticos más utilizados en la práctica usando datos
categóricos.
Se considera una respuesta binaria Y junto con una variable explicativa X de manera
que
π(x) = P (Y = 1|X = x) = 1 − P (Y = 0|X = x)

El modelo de regresión logı́stica asume que

exp [α + βx]
π(x) =
1 + exp [α + βx]
o de manera equivalente
 
π(x)
log = logit (π(x)) = α + βx
1 − π(x)

La interpretación del parámetro β se hace mediante la función eβ . Este valor es un

cociente de los odds de X = x + 1 dividido entre los odds de X = x.
El parámetro β determina la tasa de incremento o decremento de la curva en forma
de S para π(x).
En la gráfica siguiente se ve la curva en forma de S para π(x).

1
 El cambio en π(x) depende del valor de cada x. La tangente en un punto particular
de x describe la tasa de cambio en ese punto. La tangente tiene una pendiente igual a
βπ(x)(1 − π(x)). Por ejemplo, la tangente en el punto x tal que π(x) = 0,50 tiene una
pendiente igual a β · 0,50 · 0,50 = 0,25 · β; sin embargo cuando π(x) = 0,90 ó π(x) = 0,10
la pendiente es 0,09β.
La pendiente se aproxima a 0 cuando la probabilidad se aproxima a 1 o a 0.
La pendiente más pronunciada se obtiene cuando π(x) = 0,50; el correspondiente
valor de x es
α
x=− ,
β
ya que  
0,5
log = 0 = α + βx
1 − 0,5
y se denomina nivel mediano efectivo (EL50 ). Representa el valor al que los resultados
tienen una probabilidad igual a 0.50.
El parámetro α no suele tener un interés especial.

2
Ejemplo

En los cangrejos cacerola hembra: ¿el hecho de que tengan algún amante depende de
la anchura de su caparazón? ¿El tamaño importa?
La variable respuesta que interesa es si las hembras de los cangrejos cacerola tienen,
al menos, un amante. Como predictor solamente se utiliza la anchura del caparazón. Se
construye la variable a partir de la variable original que recoge el número de amantes.

tabla <- [Link] ( " http : / / [Link] / stat557 / data / [Link] " ,
header =T , sep = " \ t " )
dimnames ( tabla )[[2]] = c ( " color " ," spine " ," width " ," satell " ," weight " )
names ( tabla )

# Se crea una variable de tipo 0 -1

# Es 1 si el numero de amantes es mayor que 0
( tabla $ [Link] = ifelse ( tabla $ satell > 0 , 1 , 0))

# Se aplica un modelo de regresion logit

[Link] = glm ( tabla $ [Link] ∼ tabla $ width ,
family = binomial , data = tabla )

summary ( [Link] )

# Se discretiza la variable width

tabla $ [Link] = cut ( tabla $ width , breaks = c (0 , seq (23 .25 , 29 .25 ) , Inf ))

# Se calcula la proporcion de cangrejos que tienen al menos

# un amante para cada grupo de anchura de caparazon
prop = aggregate ( tabla $ [Link] , by = list ( W = tabla $ [Link] ) , mean ) $ x

# Se calcula la anchura media del caparazon en cada grupo

# de anchura de caparazon
plot.x = aggregate ( tabla $ width , by = list ( W = tabla $ [Link] ) , mean ) $ x

X11 ()
plot ( y = prop , x = plot.x , pch =16 , type = " p " ,
xlab = expression ( paste ( " Anchura , " , italic ( x ) , " ( cm ) " )) ,
ylab = expression ( paste ( " Proporcion de amantes que tiene " ,
{∼pi } , " ( x ) " )) , ylim = c (0 ,1) , xlim = c (20 ,34))

ind = order ( tabla $ width )

lines ( x = tabla $ width [ ind ] ,
y = predict ( [Link] , type = " response " )[ ind ] , type = " l " , lty =3)

# Intervalos de confianza para los parametros.

library ( MASS )
confint ( [Link] )

# Se estima las probabilidades de que la respuesta sea igual a 1

[Link] = predict ( [Link] , type = " response " , se = T )

# Se representan los intervalos de confianza

ind = order ( tabla $ width )

3
X11 ()
par ( family = " sans " )
plot ( tabla $ width [ ind ] , [Link] $ fit [ ind ] ,
axes =F , type = " l " , xlim = c (20 ,33) ,
ylab = " Probabilidad de tener al menos un amante " ,
xlab = expression ( paste ( " Anchura " , italic ( x ) , " ( cm ) " )) )
axis (2 , at = seq (0 ,1 ,0 .2 ))
axis (1 , at = seq (20 ,32 ,2))
lines ( tabla $ width [ ind ] , [Link] $ fit [ ind ] -
1 .96 * [Link] $ se [ ind ] , lty =3 , col = " blue " , lwd =2)
lines ( tabla $ width [ ind ] , [Link] $ fit [ ind ] +
1 .96 * [Link] $ se [ ind ] , lty =3 , col = " blue " , lwd =2)

summary ( [Link] )

4
Call :
glm ( formula = [Link] ∼ tabla $ width , family = binomial , data = tabla )

Deviance Residuals :
Min 1Q Median 3Q Max
-2 .0281 -1 .0458 0 .5480 0 .9066 1 .6942

Coefficients :
Estimate Std. Error z value Pr ( >| z |)
( Intercept ) -12 .3508 2 .6287 -4 .698 2 .62e -06 * * *
tabla $ width 0 .4972 0 .1017 4 .887 1 .02e -06 * * *
---
Signif. codes : 0 ’* * * ’ 0 .001 ’* * ’ 0 .01 ’* ’ 0 .05 ’. ’ 0 .1 ’ ’ 1

( Dispersion parameter for binomial family taken to be 1)

Null deviance : 225 .76 on 172 degrees of freedom

Residual deviance : 194 .45 on 171 degrees of freedom
AIC : 198 .45

Number of Fisher Scoring iterations : 4

confint ( [Link] )

2 .5 % 97 .5 %
( Intercept ) -17 .8100090 -7 .4572470
tabla $ width 0 .3083806 0 .7090167

5
 En este caso, el odds de la respuesta para x es
π(x)
x
= eα eβ
1 − π(x)

de modo que β interviene en los odds multiplicando por eβ cada incremento de 1 unidad
de x. Ası́, los odds al nivel x + 1 son iguales que los odds al nivel x multiplicados por eβ .
Cuando β = 0, eβ = 1, de modo que los odds no cambian cuando x cambia.
En el ejemplo se tiene que
 
π
b(x)
log π (x)) = −12,35 + 0,497x
= logit (b
1−π b(x)
De este modo, los odds estimados para cada cm de incremento de anchura en el caparazón
implica multiplicar por eβ = exp(0,497) = 1,64. Es decir, hay un incremento del 64 %
b

de la probabilidad de tener 1 amante por cada incremento de 1 cm de anchura.

Por ejemplo, con valor medio de anchura igual a x = 26,3 se obtiene

6
 exp (α + βx)
π
b(x) = =
1 + exp (α + βx)
exp (−12,35 + 0,497x)
=
1 + exp (−12,35 + 0,497x)
exp (−12,35 + 0,497 · 26,3)
= 0,6729
1 + exp (−12,35 + 0,497 · 26,3)
Ası́ los odds son
0,6729
= 2,057
1 − 0,6729
Si aumentamos una unidad x, es decir se toma x = 27,3, entonces
exp (−12,35 + 0,497 · 27,3)
π
b(x) = = 0,772
1 + exp (−12,35 + 0,497 · 27,3)
de modo que los odds son
0,772
= 3,386
1 − 0,772
de manera que  
3,386
log = log(1,646) = 0,497 = β
2,057

7
Modelos Logit con predictores categóricos
Supongamos que se tiene una variable independiente categórica (factor) X con I
categorı́as. Representamos entonces los datos en una tabla I × 2 donde en la primera
columna están los yi (i = 1, . . . , I) el número de éxitos en las ni pruebas correspondientes
a cada categorı́a y en la segunda el número de fracasos.
Se tiene, entonces, que cada yi se distribuye como una binomial Bin(πi , ni ).
El modelo logit con un factor serı́a
 
πi
logit[P (Y = 1)] = log = α + βi ,
1 − πi
de modo que un valor mayor de βi implica un mayor valor de la probabilidad πi .
La parte derecha de la expresión es semejante a un modelo ANOVA unifactorial. En
este caso hay tantos parámetros βi como categorı́as, pero uno es redundante. Es decir,
con I categorı́as basta fijar I − 1 parámetros.
Una opción habitual que se considera es fijar el valor de uno de los parámetros en 0,
por ejemplo, βI = 0.
Si los valores originales no satisfacen la restricción anterior, entonces se pueden re-
codificar para que esto se cumpla. Por ejemplo, se fijan βei = βi − βI y α e = α + βI que
satisfacen βeI = 0. De este modo,

α − βI ) + (βei + βI ) = α
logit(πi ) = α + βi = (e e + βei

Si el factor X no tiene efectos, entonces

β1 = · · · = βI = 0

o equivalentemente
π1 = · · · = πI

Una formulación alternativa con variables dummy

En la tabla original I × 2, se puede plantear una formulación alternativa si conside-

ramos que xi = 1 para las observaciones de la fila i−ésima y xi = 0 para las restantes
filas i = 1, . . . , I − 1.
De este modo, el modelo se puede expresar como

logit(πi ) = α + β1 x1 + · · · + βI−1 xI−1

Esta formulación es equivalente a la que se usa con la restricción βI = 0 para una

categorı́a arbitraria.

8
 Otra posible opción que se puede considerar es tomar
X
βi = 0
i

En el caso de I = 2, se tiene β1 = −β2 que equivale a definir una variable dummy

que vale x = 1 en la categorı́a 1 y x = −1 en la categorı́a 2.

Ejemplo

Se considera la siguiente tabla en la que se plantea la relación entre el consumo de

alcohol, con la proporción de malformaciones en el feto durante el embarazo.

Consumo de alcohol Presente Ausente

0 48 17066
<1 38 14464
1−2 5 788
3−5 1 126
≥6 1 37

Para evitar la redundancia entre los parámetros, se utiliza el comando options.

La opción de que la suma de los parámetros sea nula lo indicamos del siguiente modo:

options ( contrasts = c ( " [Link] " , " [Link] " ))

También podemos indicar que el parámetro asociado a la primera categorı́a sea nulo:

options ( contrasts = c ( " [Link] " , " [Link] " ))

Vamos a estudiar la posible asociación entre el consumo de alcohol y las malforma-

ciones en los fetos.

Alcohol = factor ( c ( " 0 " , " <1 " , " 1 -2 " , " 3 -5 " , " >=6 " ) ,
levels = c ( " 0 " ," <1 " , " 1 -2 " , " 3 -5 " , " >=6 " ))

malformaciones = c (48 , 38 , 5 , 1 , 1)

n = c (17066 , 14464 , 788 , 126 , 37) + malformaciones

Primero, indicamos que β1 ha de ser nula y luego se define un modelo logit.

options ( contrasts = c ( " [Link] " , " [Link] " ))

# Ajustamos un modelo logit

( tabla.logit1 = glm ( malformaciones / n ∼ Alcohol ,
family = binomial , weights = n ))

Se obtiene

9
 Call : glm ( formula = malfor macione s / n ∼ Alcohol , family = binomial ,
weights = n )

Coefficients :
( Intercept ) Alcohol <1 Alcohol1 -2 Alcohol3 -5 Alcohol >=6
-5 .87364 -0 .06819 0 .81358 1 .03736 2 .26272

Degrees of Freedom : 4 Total ( i.e. Null ) ; 0 Residual

Null Deviance : 6 .202
Residual Deviance : 5 .995e -15 AIC : 28 .63

Ahora, de manera alternativa, se fija que la última categorı́a sea nula (βI = 0)

revAlcohol = factor ( c ( " 0 " , " <1 " , " 1 -2 " , " 3 -5 " , " >=6 " ) ,
levels = rev ( c ( " 0 " ," <1 " , " 1 -2 " , " 3 -5 " , " >=6 " )))

( tabla.logit2 = glm ( malformaciones / n ∼ revAlcohol ,

family = binomial , weights = n ))

Se obtiene

Call : glm ( formula = malfor macione s / n ∼ revAlcohol , family = binomial ,

weights = n )

Coefficients :
( Intercept ) revAlcohol3 -5 revAlcohol1 -2 revAlcohol <1 revAlcohol0
-3 .611 -1 .225 -1 .449 -2 .331 -2 .263

Degrees of Freedom : 4 Total ( i.e. Null ) ; 0 Residual

Null Deviance : 6 .202
Residual Deviance : 6 .084e -14 AIC : 28 .63

Aunque los valores de los parámetros son diferentes, las predicciones de las probabi-
lidades son las mismas con ambos modelos; es decir, no influye el hecho de fijar en 0 un
efecto u otro.

cbind ( logit = predict ( tabla.logit1 ) ,

[Link] = predict ( tabla.logit1 , type = " response " ))

logit [Link]
1 -5 .873642 0 .002804721
2 -5 .941832 0 .002620328
3 -5 .060060 0 .006305170
4 -4 .836282 0 .007874016
5 -3 .610918 0 .026315789

cbind ( logit = predict ( tabla.logit2 ) ,

[Link] = predict ( tabla.logit2 , type = " response " ))

10
 logit [Link]
1 -5 .873642 0 .002804721
2 -5 .941832 0 .002620328
3 -5 .060060 0 .006305170
4 -4 .836282 0 .007874016
5 -3 .610918 0 .026315789

Podemos observar cómo la proporción de malformaciones se incrementa conforme

aumenta el consumo de alcohol.
Pero, en el análisis previo no se ha tenido en cuenta la naturaleza ordinal de la
covariable consumo de alcohol. Ası́, alternativamente, se puede considerar un modelo
alternativo escrito como
logit(πi ) = α + βxi

de modo que la hipótesis de independencia equivale a β = 0.

Fijamos como valores de xi :
x1 = 0; x2 = 0,5; x3 = 1,5; x4 = 4; x5 = 7
Donde el último valor es arbitrario.

puntuaciones = c (0 , 0 .5 , 1 .5 , 4 , 7)
[Link] = glm ( malformaciones / n ∼ puntuaciones ,
family = binomial , weights = n )
summary ( [Link] )

Se obtiene

Call :
glm ( formula = malf ormacio nes / n ∼ puntuaciones , family = binomial ,
weights = n )

Deviance Residuals :
1 2 3 4 5
0 .5921 -0 .8801 0 .8865 -0 .1449 0 .1291

Coefficients :
Estimate Std. Error z value Pr ( >| z |)
( Intercept ) -5 .9605 0 .1154 -51 .637 <2e -16 * * *
puntuaciones 0 .3166 0 .1254 2 .523 0 .0116 *
---
Signif. codes : 0 ’* * * ’ 0 .001 ’* * ’ 0 .01 ’* ’ 0 .05 ’. ’ 0 .1 ’ ’ 1

( Dispersion parameter for binomial family taken to be 1)

Null deviance : 6 .2020 on 4 degrees of freedom

Residual deviance : 1 .9487 on 3 degrees of freedom
AIC : 24 .576

Number of Fisher Scoring iterations : 4

Se pueden obtener los logits y las proporciones ajustadas:

11
 cbind ( logit = predict ( [Link] ) ,
[Link] = predict ( [Link] , type = " response " ))

logit [Link]
1 -5 .960461 0 .002572090
2 -5 .802181 0 .003011861
3 -5 .485620 0 .004128844
4 -4 .694219 0 .009065077
5 -3 .744538 0 .023100302

También se pueden obtener los intervalos de confianza para los parámetros:

library ( MASS )
confint ( [Link] )

2 .5 % 97 .5 %
( Intercept ) -6 .19302366 -5 .7396968
puntuaciones 0 .01868149 0 .5234947

Regresión logı́stica con múltiples predictores categóricos

Tanto en la regresión logı́stica, como en la regresión ordinaria, habitualmente se
trabaja con varias variables predictoras o explicativas, es decir, en general se estudian
varios factores al mismo tiempo.
Supongamos que una variable de respuesta binaria Y tiene dos variables predictoras,
X y Z. Los datos se pueden presentar, entonces, en forma de una tabla de contingencia
de orden 2 × 2 × 2.
Supongamos que x y z toman valores 0 ó 1 para representar las dos categorı́as de
cada variable. El modelo, entonces, se escribe como

logit [P (Y = 1)] = α + β1 x + β2 z

y se puede visualizar en forma de una tabla como

x z Logit
0 0 α
1 0 α + β1
0 1 α + β2
1 1 α + β1 + β2

12
 Este modelo supone que no hay interacción, es decir, el efecto de un factor es igual
a lo largo de cada categorı́a del otro factor.
Para un valor fijo de la categorı́a z de Z el efecto del cambio de x = 0 a x = 1 es

[α + β1 (1) + β2 z] − [α + β1 (0) + β2 z] = β1

Ası́, para Z fijada, el odds para el éxito en x = 1 es igual a exp(β1 ) veces al odds del
éxito para x = 0.
Existe independencia condicional entre X e Y , controlando o fijando Z, si β1 = 0.
En este caso, la razón de odds común es igual a 1. Se puede aplicar, entonces el modelo
más simple a la tabla de tres variables:

logit [P (Y = 1)] = α + β2 z.

Ejemplo: SIDA y AZT

Se considera un estudio antiguo (1991) que trataba el efecto de un retroviral (AZT ) en

el desarrollo de los sı́ntomas del SIDA. Se obtuvo una muestra de 338 personas afectadas
de SIDA en donde como respuesta consideramos si desarrollan sı́ntomas de SIDA y como
predictores la raza y si se les administra AZT de modo inmediato o cuando los linfocitos
tipo T muestran debilidad inmune.

Sı́ntomas
Raza Uso AZT Si No
Blanca Si 14 93
No 32 81
Negra Si 11 52
No 12 43

Denominamos como X ≡ AZT, Z ≡ Raza e Y ≡ SIDA.

En el modelo, consideramos que y = 1 si se desarrollan los sı́ntomas del SIDA,
y = 0 si no se desarrollan. Del mismo modo, x = 1 son aquellos que tomaron AZT
inmediatamente, x = 0 para los que no, y z = 1 son personas de raza blanca y z = 0
son personas de color. Ası́, el modelo se puede expresar como
AZT
logit [P (Y = 1)] = α + βSi + βblanc

En este caso, α es la razón de odds correspondiente a desarrollar sı́ntomas de SIDA

para personas de color a los que no se les administró inmediatamente AZT.
AZT
Por otro lado, βSi es el incremento en la razón de odds para los que usan inmedia-
tamente AZT y βblanc es el incremento de los odds para las personas blancas.

13
 [Link] = [Link] ( AZT = factor ( c ( " Si " , " No " ) ,
levels = c ( " No " , " Si " )) ,
Raza = factor ( c ( " Blanca " , " Negra " ) ,
levels = c ( " Negra " , " Blanca " )))

tabla = [Link] ( [Link] , Si = c (14 , 32 , 11 , 12) ,

No = c (93 , 81 , 52 , 43))

options ( contrasts = c ( " [Link] " , " [Link] " ))

fit1 = glm ( cbind ( Si , No ) ∼ AZT + Raza , family = binomial , data = tabla )

summary ( fit1 )

Call :
glm ( formula = cbind ( Si , No ) ∼ AZT + Raza , family = binomial ,
data = tabla )

Deviance Residuals :
1 2 3 4
-0 .5547 0 .4253 0 .7035 -0 .6326

Coefficients :
Estimate Std. Error z value Pr ( >| z |)
( Intercept ) -1 .07357 0 .26294 -4 .083 4 .45e -05 * * *
AZTSi -0 .71946 0 .27898 -2 .579 0 .00991 * *
RazaBlanca 0 .05548 0 .28861 0 .192 0 .84755
---
Signif. codes : 0 ’* * * ’ 0 .001 ’* * ’ 0 .01 ’* ’ 0 .05 ’. ’ 0 .1 ’ ’ 1

( Dispersion parameter for binomial family taken to be 1)

Null deviance : 8 .3499 on 3 degrees of freedom

Residual deviance : 1 .3835 on 1 degrees of freedom
AIC : 24 .86

Number of Fisher Scoring iterations : 4

Por tanto, la razón de odds para los que usan AZT inmediatamente después de
desarrollar sı́ntomas de SIDA serı́a exp(−0,71946) = 0,49. Un intervalo de confianza se
puede obtener con

confint ( fit1 )

2 .5 % 97 .5 %
( Intercept ) -1 .6088054 -0 .5734959
AZTSi -1 .2773237 -0 .1798769
RazaBlanca -0 .5022982 0 .6334104

Veamos si la raza y la aparición de sı́ntomas son condicionalmente independientes

dado o fijado el tratamiento por AZT.

fitSimple = update ( object = fit1 , formula = ∼. - Raza )

anova ( fit1 , fitSimple , test = " Chisq " )

14
Analysis of Deviance Table

Model 1: cbind ( Si , No ) ∼ AZT + Raza

Model 2: cbind ( Si , No ) ∼ AZT
Resid. Df Resid. Dev Df Deviance Pr ( > Chi )
1 1 1 .3835
2 2 1 .4206 -1 -0 .037084 0 .8473

No es significativa la raza en el modelo.

Opción con SAS:

OPTIONS nodate ls =65 formchar = ’ | - - - -|+| - - -+=| -/\ < >* ’;

x ’ cd " c :\ DeSAS " ’;

DATA sida ;
input raza $ azt $ si no @@ ;
casos = si + no ;
DATALINES ;
blanca y 14 93 blanca n 32 81
negra y 11 52 negra n 12 43
;
PROC genmod order = data ;
class raza azt ;
model si / casos = raza azt / dist = bin link = logit obstats type3 ;
RUN ;

15
 The GENMOD Procedure

Class Level Information

Class Levels Values

raza 2 blanca negra
azt 2 y n

Criteria For Assessing Goodness Of Fit

Criterion DF Value Value / DF

Deviance 1 1 .3835 1 .3835

Scaled Deviance 1 1 .3835 1 .3835
Pearson Chi - Square 1 1 .3910 1 .3910
Scaled Pearson X2 1 1 .3910 1 .3910
Log Likelihood -167 .5756

Algorithm converged.

Analysis Of Parameter Estimates

Standard Wald 95 %
Parameter DF Estimate Error Confidence Limits

Intercept 1 -1 .0736 0 .2629 -1 .5889 -0 .5582

Analysis Of Parameter Estimates

Chi -
Parameter Square Pr > ChiSq

Intercept 16 .67 < .0001

Analysis Of Parameter Estimates

Standard Wald 95 %
Parameter DF Estimate Error Confidence Limits

raza blanca 1 0 .0555 0 .2886 -0 .5102 0 .6212

raza negra 0 0 .0000 0 .0000 0 .0000 0 .0000
azt y 1 -0 .7195 0 .2790 -1 .2662 -0 .1727
azt n 0 0 .0000 0 .0000 0 .0000 0 .0000
Scale 0 1 .0000 0 .0000 1 .0000 1 .0000

16
 Analysis Of Parameter Estimates

Chi -
Parameter Square Pr > ChiSq

raza blanca 0 .04 0 .8476

raza negra . .
azt y 6 .65 0 .0099
azt n . .
Scale

NOTE : The scale parameter was held fixed.

LR Statistics For Type 3 Analysis

Chi -
Source DF Square Pr > ChiSq

raza 1 0 .04 0 .8473

azt 1 6 .87 0 .0088

17
Ejemplo: Cangrejos cacerola

Se considera el ejemplo de los cangrejos cacerola y se aplica el siguiente modelo:

logit [π] = α + β1 c1 + β2 c2 + β3 c3 + β4 x

donde

– π = P (Y = 1)

– x ≡ anchura (cm).

– c1 = 1 para color medio-claro y 0 en otro caso.

– c2 = 1 para color medio y 0 en otro caso.

– c3 = 1 para color medio-oscuro y 0 en otro caso.

Se trata de ver la relación entre el color del caparazón y la probabilidad de tener

algún amante.

tabla <- [Link] ( " http : / / [Link] / stat557 / data / [Link] " ,
header =T , sep = " \ t " )
dimnames ( tabla )[[2]] = c ( " color " ," spine " ," width " ," satell " ," weight " )
names ( tabla )

options ( contrasts = c ( " [Link] " , " [Link] " ))

# Definimos como factor la variable que nos indica el color.

tabla $ [Link] <- factor ( tabla $ color , levels = c ( " 5 " ," 4 " ," 3 " ," 2 " ) ,
labels = c ( " oscuro " , " med - oscuro " , " med " , " med - claro " ))

# Se define una variable binaria : tiene o no tiene amantes

tabla $ [Link] = ifelse ( tabla $ satell > 0 , 1 , 0)

[Link] <- glm ( [Link] ∼ [Link] + width ,

family = binomial , data = tabla )
summary ( [Link] , cor = F )

library ( MASS )
confint ( [Link] )

res1 = predict ( [Link] , type = " response " ,

newdata = [Link] ( width = seq (18 , 34 , 1) , [Link] = " med - claro " ))

res2 = predict ( [Link] , type = " response " ,

newdata = [Link] ( width = seq (18 , 34 , 1) , [Link] = " med " ))

res3 = predict ( [Link] , type = " response " ,

newdata = [Link] ( width = seq (18 , 34 , 1) , [Link] = " med - oscuro " ))

res4 = predict ( [Link] , type = " response " ,

18
newdata = [Link] ( width = seq (18 , 34 , 1) , [Link] = " oscuro " ))

plot ( seq (18 , 34 , 1) , res1 , type = " l " , bty = " L " ,
ylab = " Probabilidad Predicha " , axes =F ,
xlab = expression ( paste ( " Ancho , " , italic ( x ) , " ( cm ) " )) , col = " yellow " )

axis (2 , at = seq (0 , 1 , 0 .2 ))
axis (1 , at = seq (18 , 34 , 2))

lines ( seq (18 , 34 , 1) , res2 , col = " orange " )

lines ( seq (18 , 34 , 1) , res3 , col = " red " )

lines ( seq (18 , 34 , 1) , res4 , col = " darkred " )

arrows ( x0 =29 , res1 [25 - 17] , x1 =25 , y1 = res1 [25 - 17] , length =0 .09 )

text ( x =29 .1 , y = res1 [25 - 17] , " Color 1 " , adj = c (0 , 0) , col = " yellow " )

arrows ( x0 =23 , res2 [26 - 17] , x1 =26 , y1 = res2 [26 - 17] , length =0 .09 )

text ( x =21 .1 , y = res2 [26 - 17] , " Color 2 " , adj = c (0 , 0) , col = " orange " )

arrows ( x0 =28 .9 , res3 [24 - 17] , x1 = 24 , y1 = res3 [24 - 17] ,

length =0 .09 )

text ( x =29 , y = res3 [24 - 17] , " Color 3 " , adj = c (0 , 0) , col = " red " )

arrows ( x0 =25 .9 , res4 [23 - 17] , x1 =23 , y1 = res4 [23 - 17] ,

length =0 .09 )

text ( x =26 , y = res4 [23 - 17] , " Color 4 " , adj = c (0 , 0) ,

col = " darkred " )

19
Call :
glm ( formula = [Link] ∼ [Link] + width , family = binomial ,
data = tabla )

Deviance Residuals :
Min 1Q Median 3Q Max
-2 .1124 -0 .9848 0 .5243 0 .8513 2 .1413

Coefficients :
Estimate Std. Error z value Pr ( >| z |)
( Intercept ) -12 .7151 2 .7617 -4 .604 4 .14e -06 ***
[Link] - oscuro 1 .1061 0 .5921 1 .868 0 .0617 .
[Link] 1 .4023 0 .5484 2 .557 0 .0106 *
[Link] - claro 1 .3299 0 .8525 1 .560 0 .1188
width 0 .4680 0 .1055 4 .434 9 .26e -06 ***
---
Signif. codes : 0 ’* * * ’ 0 .001 ’* * ’ 0 .01 ’* ’ 0 .05 ’. ’ 0 .1 ’ ’ 1

( Dispersion parameter for binomial family taken to be 1)

Null deviance : 225 .76 on 172 degrees of freedom

Residual deviance : 187 .46 on 168 degrees of freedom
AIC : 197 .46

Number of Fisher Scoring iterations : 4

20
 2 .5 % 97 .5 %
( Intercept ) -18 .45674069 -7 .5788795
[Link] - oscuro -0 .02792233 2 .3138635
[Link] 0 .35269965 2 .5260703
[Link] - claro -0 .27377584 3 .1356611
width 0 .27128167 0 .6870436

Se puede considerar un modelo con interacción color y anchura del caparazón.

cr ab .f it .l og is t. ia = update ( object = [Link] ,

formula = ∼. + width : [Link] )
anova ( [Link] , [Link] , test = " Chisq " )

Analysis of Deviance Table

Model 1: [Link] ∼ [Link] + width

Model 2: [Link] ∼ [Link] + width + [Link] : width
Resid. Df Resid. Dev Df Deviance Pr ( > Chi )
1 168 187 .46
2 165 183 .08 3 4 .3764 0 .2236

Ası́, la interacción no es significativa.

21
Consideramos ahora el modelo, solo con la variable anchura, usando SAS:

OPTIONS nodate ls =65 formchar = ’ | - - - -|+| - - -+=| -/\ < >* ’;

x ’ cd " c :\ deSAS " ’;

DATA crabs ;
input anchura casos amante ;
DATALINES ;
22.69 14 5
23.84 14 4
24.77 28 17
25.84 39 21
26.79 22 15
27.74 24 20
28.67 18 15
30.41 14 14
;

PROC genmod ;
model amante / casos = anchura / dist = bin link = logit waldci
lrci alpha =.01;

PROC logistic ;
model amante / casos = anchura / influence ;

OUTPUT out = predict p = pi _ hat lower = LCL upper = UCL ;

PROC print data = predict ;

RUN ;

22
 The GENMOD Procedure

Criteria For Assessing Goodness Of Fit

Criterion DF Value Value / DF

Deviance 6 5 .9696 0 .9949

Scaled Deviance 6 5 .9696 0 .9949
Pearson Chi - Square 6 5 .0292 0 .8382
Scaled Pearson X2 6 5 .0292 0 .8382
Log Likelihood -98 .8470

Algorithm converged.

Analysis Of Parameter Estimates

Standard Wald 99 % Confidence

Parameter DF Estimate Error Limits

Intercept 1 -11 .5140 2 .5489 -18 .0796 -4 .9485

anchura 1 0 .4647 0 .0985 0 .2108 0 .7185
Scale 0 1 .0000 0 .0000 1 .0000 1 .0000

Likelihood Ratio 99 % Chi -

Parameter Confidence Limits Square Pr > ChiSq

Intercept -18 .5807 -5 .3345 20 .41 < .0001

anchura 0 .2267 0 .7388 22 .23 < .0001
Scale 1 .0000 1 .0000

NOTE : The scale parameter was held fixed.

Response Profile

Ordered Binary Total

Value Outcome Frequency

1 Event 111
2 Nonevent 62

Model Fit Statistics

Intercept
Intercept and
Criterion Only Covariates

AIC 227 .759 201 .694

SC 230 .912 208 .001
-2 Log L 225 .759 197 .694

23
 Testing Global Null Hypothesis : BETA =0

Test Chi - Square DF Pr > ChiSq

Likelihood Ratio 28 .0644 1 < .0001

Score 25 .6828 1 < .0001
Wald 22 .2312 1 < .0001

Analysis of Maximum Likelihood Estimates

Standard Wald
Parameter DF Estimate Error Chi - Square Pr > ChiSq

Intercept 1 -11 .5128 2 .5488 20 .4031 < .0001

anchura 1 0 .4646 0 .0985 22 .2312 < .0001

Odds Ratio Estimates

Point 95 % Wald
Effect Estimate Confidence Limits

anchura 1 .591 1 .312 1 .930

Association of Predicted Probabilities and Observed Responses

Percent Concordant 66 .3 Somers ’ D 0 .454

Percent Discordant 20 .9 Gamma 0 .520
Percent Tied 12 .8 Tau - a 0 .210
Pairs 6882 c 0 .727

Obs anchura casos amante pi _ hat LCL UCL

1 22 .69 14 5 0 .27481 0 .15971 0 .43036

2 23 .84 14 4 0 .39268 0 .28007 0 .51800
3 24 .77 28 17 0 .49901 0 .40218 0 .59592
4 25 .84 39 21 0 .62086 0 .53879 0 .69655
5 26 .79 22 15 0 .71801 0 .63346 0 .78953
6 27 .74 24 20 0 .79835 0 .70525 0 .86756
7 28 .67 18 15 0 .85913 0 .76133 0 .92102
8 30 .41 14 14 0 .93192 0 .84095 0 .97256

24
Ejemplo admsiones en postgrado de una universidad con SAS

Un investigador/a está interesado en analizar cómo las variables gre (Graduate Re-
cord Exam), gpa (promedio de calificaciones) y el prestigio de la universidad, afectan
en la admisión en unos estudios de posgrado. La variable en estudio admitir/no admitir,
es binaria.
Este conjunto de datos tiene una variable respuesta binaria llamada admit, que es
igual a 1 si el individuo fue admitido en los estudios de posgrado, y 0 en caso contrario.
Hay tres variables predictoras: gre, gpa, y rank. Se consideran gre y gpa como conti-
nuas. La variable rank toma los valores de 1 a 4. Las instituciones con un rango igual
a 1 tienen el prestigio más alto, mientras que aquellos con un rango de 4 tienen el más
bajo.
Los datos se encuentran en

[Link]

Se consideran primero algunos estadı́sticos descriptivas.

OPTIONS nodate ls =65 formchar = ’ | - - - -|+| - - -+=| -/\ < >* ’;

x ’ cd " c :\ deSAS " ’;

PROC means data =" binary ";

var gre gpa ;

PROC freq data =" binary ";

TABLES rank admit admit * rank ;
RUN ;

25
 The MEANS Procedure

Variable N Mean Std Dev Minimum

---------------------------------------------------------------
GRE 400 587 .7000000 115 .5165364 220 .0000000
GPA 400 3 .3899000 0 .3805668 2 .2600000
---------------------------------------------------------------

Variable Maximum
------------------------
GRE 800 .0000000
GPA 4 .0000000
------------------------

The FREQ Procedure

Cumulative Cumulative
RANK Frequency Percent Frequency Percent
---------------------------------------------------------
1 61 15 .25 61 15 .25
2 151 37 .75 212 53 .00
3 121 30 .25 333 83 .25
4 67 16 .75 400 100 .00

Cumulative Cumulative
ADMIT Frequency Percent Frequency Percent
----------------------------------------------------------
0 273 68 .25 273 68 .25
1 127 31 .75 400 100 .00

Table of ADMIT by RANK

ADMIT RANK

Frequency |
Percent |
Row Pct |
Col Pct | 1| 2| 3| 4| Total
- - - - - - - - -+ - - - - - - - -+ - - - - - - - -+ - - - - - - - -+ - - - - - - - -+
0 | 28 | 97 | 93 | 55 | 273
| 7 .00 | 24 .25 | 23 .25 | 13 .75 | 68 .25
| 10 .26 | 35 .53 | 34 .07 | 20 .15 |
| 45 .90 | 64 .24 | 76 .86 | 82 .09 |
- - - - - - - - -+ - - - - - - - -+ - - - - - - - -+ - - - - - - - -+ - - - - - - - -+
1 | 33 | 54 | 28 | 12 | 127
| 8 .25 | 13 .50 | 7 .00 | 3 .00 | 31 .75
| 25 .98 | 42 .52 | 22 .05 | 9 .45 |
| 54 .10 | 35 .76 | 23 .14 | 17 .91 |
- - - - - - - - -+ - - - - - - - -+ - - - - - - - -+ - - - - - - - -+ - - - - - - - -+
Total 61 151 121 67 400
15 .25 37 .75 30 .25 16 .75 100 .00

26
 Para modelizar los 1’s en lugar de los 0’s, se usa la opción descending. Esto se hace
ası́, porque de manera predeterminada, el PROC logistic modeliza los 0’s en lugar de
los 1’s. En este caso, eso significarı́a modelizar la predicción de la probabilidad de no
entrar en el curso (admit=0) frente a entrar (admit=1).
Ambos modelos son equivalentes pero tiene conceptualmente más sentido modelizar
la probabilidad de entrar el curso que de no hacerlo. La opción param=ref indica la
codificación mediante variables dummy para los niveles de rank.

OPTIONS nodate ls =70;

proc logistic data =" binary " descending ;
class rank / param = ref ;
model admit = gre gpa rank ;
RUN ;

The LOGISTIC Procedure

Response Profile

Ordered Total
Value ADMIT Frequency

1 1 127
2 0 273

Probability modeled is ADMIT =1 .

Class Level Information

Class Value Design Variables

RANK 1 1 0 0
2 0 1 0
3 0 0 1
4 0 0 0

Model Fit Statistics

Intercept
Intercept and
Criterion Only Covariates

AIC 501 .977 470 .517

SC 505 .968 494 .466
-2 Log L 499 .977 458 .517

27
 Testing Global Null Hypothesis : BETA =0

Test Chi - Square DF Pr > ChiSq

Likelihood Ratio 41 .4590 5 < .0001

Score 40 .1603 5 < .0001
Wald 36 .1390 5 < .0001

Type 3 Analysis of Effects

Wald
Effect DF Chi - Square Pr > ChiSq

GRE 1 4 .2842 0 .0385

GPA 1 5 .8714 0 .0154
RANK 3 20 .8949 0 .0001

En la salida muestra el ajuste del modelo. El valor de -2 Log L = 499.977 se puede

usar para hacer comparaciones en modelos anidados.
La razón de verosimilitudes chi-cuadrado igual a 41.4590 con un p-valor igual a 0.0001
indica que el modelo completo tiene un ajuste mejor que el modelo con una constante.
En la parte del análisis de efectos de Type 3, se muestran los contrastes de hipótesis
para cada una de las variables de manera individual. Los p-valores de los estadı́sticos
muestran que cada una de las variables aumenta el ajuste del modelo.

Analysis of Maximum Likelihood Estimates

Standard Wald
Parameter DF Estimate Error Chi - Square Pr > ChiSq

Intercept 1 -5 .5414 1 .1381 23 .7081 < .0001

GRE 1 0 .00226 0 .00109 4 .2842 0 .0385
GPA 1 0 .8040 0 .3318 5 .8714 0 .0154
RANK 1 1 1 .5514 0 .4178 13 .7870 0 .0002
RANK 2 1 0 .8760 0 .3667 5 .7056 0 .0169
RANK 3 1 0 .2112 0 .3929 0 .2891 0 .5908

Se muestran los coeficientes, los errores estándar, los estadı́sticos y sus p-valores. Los
coeficientes para gre y gpa son significativos, ası́ como los términos correspondientes a
rank=1 y rank=2 (con respecto a la categorı́a base rank=4).
Los coeficientes de la regresión logı́stica muestran el cambio en el logaritmo de los
odds cuando se incrementa en una unidad la variable predictora correspondiente.

28
 Por cada cambio en una unidad de gre, el logaritmo del odds de la admisión frente
a la no admisión se incrementa en 0.002.
Por cada cambio en una unidad de gpa, el logaritmo del odds de la admisión frente
a la no admisión se incrementa en 0.804.
Los coeficientes para las categorı́as de rank tienen una interpretación algo diferente.
Por ejemplo, el hecho de venir de una universidad con rango igual a 1, frente a otra con
rango igual a 4, incrementa el logaritmo del odds de la admisión por 1.55.

Odds Ratio Estimates

Point 95 % Wald
Effect Estimate Confidence Limits

GRE 1 .002 1 .000 1 .004

GPA 2 .235 1 .166 4 .282
RANK 1 vs 4 4 .718 2 .080 10 .701
RANK 2 vs 4 2 .401 1 .170 4 .927
RANK 3 vs 4 1 .235 0 .572 2 .668

Association of Predicted Probabilities and Observed Responses

Percent Concordant 69 .1 Somers ’ D 0 .386

Percent Discordant 30 .6 Gamma 0 .387
Percent Tied 0 .3 Tau - a 0 .168
Pairs 34671 c 0 .693

Se obtienen las razones de odds, (la exponencial de los coeficientes de la regresión

logı́stica) que se pueden interpretar como el cambio multiplicativo en los odds por cada
cambio de una unidad en una variable explicativa. Por ejemplo, para el incremento de
una unidad en gpa, los odds de ser admitidos en el curso (frente a no ser admitidos) se
incrementa en un factor de 2.24.
Se muestra un test para el efecto total de la variable rank, ası́ como los coeficientes
que describen la diferencia entre el grupo de referencia (rank=4) y cada uno de los tres
otros grupos.
Se contrasta también las diferencias entre los otros niveles de rank. Por ejemplo, se
puede contrastar la diferencia entre los coeficientes para rank=2 y rank=3, es decir,
comparar los odds de la admisión para estudiantes que proceden de una universidad de
rango 2 frente a una de rango 3.

29
 proc logistic data =" binary " descending ;
class rank / param = ref ;
model admit = gre gpa rank ;
contrast ’ rank 2 vs 3 ’ rank 0 1 -1 / estimate = parm ;
run ;

Wald
Contrast DF Chi - Square Pr > ChiSq

rank 2 vs 3 1 5 .5052 0 .0190

Contrast Rows Estimation and Testing Results

Standard
Contrast Type Row Estimate Error Alpha

rank 2 vs 3 PARM 1 0 .6648 0 .2833 0 .05

Observando el p-valor, el coeficiente para rank=2 es significativamente diferente del

coeficiente para rank=3.
Se obtiene también la estimación de de la diferencia que es 0.6648, lo que indica que
el hecho de asistir a una universidad con rango igual a 2 frente a una universidad con
rango igual a 3, incrementa el logaritmo de los odds de la admisión por 0.6648, lo que
corresponderı́a a una razón de odds igual a 1.944102.

Se pueden usar también las probabilidades predichas para interpretar mejor el mo-
delo. Se estiman las probabilidades predichas de admisión cuando gre cambia de 200 a
800 (en incrementos de 100). Cuando se estiman las probabilidades predichas se man-
tiene gpa constante en 3.39 (que es su media), y rank en 2 (poniendo rank 0 1 0).
El término intercept seguido de un 1 indica que se incluye también en el modelo la
constante.

30
proc logistic data =" binary " descending ;
class rank / param = ref ;
model admit = gre gpa rank ;
contrast ’ gre =200 ’ intercept 1 gre 200 gpa 3.3899
rank 0 1 0 / estimate = prob ;
contrast ’ gre =300 ’ intercept 1 gre 300 gpa 3.3899
rank 0 1 0 / estimate = prob ;
contrast ’ gre =400 ’ intercept 1 gre 400 gpa 3.3899
rank 0 1 0 / estimate = prob ;
contrast ’ gre =500 ’ intercept 1 gre 500 gpa 3.3899
rank 0 1 0 / estimate = prob ;
contrast ’ gre =600 ’ intercept 1 gre 600 gpa 3.3899
rank 0 1 0 / estimate = prob ;
contrast ’ gre =700 ’ intercept 1 gre 700 gpa 3.3899
rank 0 1 0 / estimate = prob ;
contrast ’ gre =800 ’ intercept 1 gre 800 gpa 3.3899
rank 0 1 0 / estimate = prob ;
run ;

Contrast Test Results

Wald
Contrast DF Chi - Square Pr > ChiSq

gre =200 1 9 .7752 0 .0018

gre =300 1 11 .2483 0 .0008
gre =400 1 13 .3231 0 .0003
gre =500 1 15 .0984 0 .0001
gre =600 1 11 .2291 0 .0008
gre =700 1 3 .0769 0 .0794
gre =800 1 0 .2175 0 .6409

Contrast Rows Estimation and Testing Results

Standard
Contrast Type Row Estimate Error Alpha Confidence Limits

gre =200 PROB 1 0 .1844 0 .0715 0 .05 0 .0817 0 .3648

gre =300 PROB 1 0 .2209 0 .0647 0 .05 0 .1195 0 .3719
gre =400 PROB 1 0 .2623 0 .0548 0 .05 0 .1695 0 .3825
gre =500 PROB 1 0 .3084 0 .0443 0 .05 0 .2288 0 .4013
gre =600 PROB 1 0 .3587 0 .0399 0 .05 0 .2847 0 .4400
gre =700 PROB 1 0 .4122 0 .0490 0 .05 0 .3206 0 .5104
gre =800 PROB 1 0 .4680 0 .0685 0 .05 0 .3391 0 .6013

31
 Contrast Rows Estimation and Testing Results

Wald
Contrast Type Row Chi - Square Pr > ChiSq

gre =200 PROB 1 9 .7752 0 .0018

gre =300 PROB 1 11 .2483 0 .0008
gre =400 PROB 1 13 .3231 0 .0003
gre =500 PROB 1 15 .0984 0 .0001
gre =600 PROB 1 11 .2291 0 .0008
gre =700 PROB 1 3 .0769 0 .0794
gre =800 PROB 1 0 .2175 0 .6409

En las probabilidades predichas se observa que la estimación de la probabilidad de

ser admitido es solo de 0.18 si la puntuación de alguien en gre es 200, pero se incrementa
hasta 0.47 si la puntuación de gre es 800, manteniendo gpa en su media (3.39) y rank
en 2.

32
Estrategias en la selección de modelos
Para un conjunto de datos con una respuesta binaria y múltiples variables predic-
toras, una pregunta fundamental es ¿cómo seleccionar un modelo de regresión logı́stica
adecuado?
Se presentan los mismos problemas que en la regresión ordinaria. Ası́, el proceso de
selección de un modelo adecuado se hace difı́cil cuando el número de variables explica-
tivas aumenta, porque aumentan a su vez los posibles efectos e interacciones.
Hay dos objetivos contrapuestos: el modelo debe ser lo suficientemente complejo
como para adaptarse bien a los datos, pero los modelos más simples son más fáciles de
interpretar. Es decir, es mejor encontrar una versión suavizada en lugar de un sobreajuste
de los datos.

¿Cuántas variables predictoras se deben introducir en el modelo? Los datos

se dice que son no balanceados en Y si y = 1 ó y = 0 aparecen relativamente pocas
veces. Esto limitarı́a el número de variables predictoras cuyos efectos se pueden estimar
de manera precisa.
Como una regla heurı́stica, se considera que deberı́a haber al menos 10 observaciones
de 1 ó 0 por cada variable predictora. Por ejemplo, si y = 1 sólo 30 veces en n = 1000
observaciones, el modelo no deberı́a tener más de tres variables predictoras, aunque el
tamaño total de la muestra fuera grande.
En general, modelos con muchas variables predictoras pueden presentar colinealida-
des lo cual hace necesario eliminar las variables que sean redundantes.
Como ocurre en regresión ordinaria, se pueden usar diferentes algoritmos para incluir
o eliminar variables en un modelo de regresión logı́stica. Ası́, se pueden usar métodos
de:

– Selección hacia adelante (forward ).

– Selección hacia atrás (backward ).

– Selección stepwise.

– Información de Akaike:

AIC = −2×(máximo de la log-verosimilitud) – número de parámetros.

33
Ejemplo: Cangrejos cacerola

En los datos de los cangrejos cacerola, se tienen cuatro variables predictoras: color (4
categorı́as), tipo de espina (3 categorı́as), peso, y anchura del caparazón. Se plantea un
modelo de regresión logı́stica para predecir si una hembra tiene amantes o no. Es decir,
y = 1 si tiene al menos 1 amante e y = 0 en caso contrario.
Se denominan (c1 , c2 , c3 ) las variables indicadoras de los tres primeros colores de
más oscuro a más claro (respecto al total de 4 colores) y (s1 , s2 ) denota a las variables
indicadoras de las dos primeras condiciones de la espina (de las 3 que hay).
Se considera el modelo

logit(P (Y = 1)) = α + β1 P eso + β2 Ancho + β3 c1 + β4 c2 + β5 c3 + β6 s1 + β7 s2

tabla <- [Link] ( " http : / / [Link] / stat557 / data / [Link] " ,
header =T , sep = " \ t " )
dimnames ( tabla )[[2]] = c ( " color " ," spine " ," width " ," satell " ," weight " )
names ( tabla )

# Se define como factor la variable que nos indica el color.

tabla $ [Link] <- factor ( tabla $ color , levels = c ( " 5 " ," 4 " ," 3 " ," 2 " ) ,
labels = c ( " oscuro " , " med - oscuro " , " med " , " med - claro " ))

# Se define una variable binaria : tiene o no tiene amantes

tabla $ [Link] = ifelse ( tabla $ satell > 0 , 1 , 0)

options ( contrasts = c ( " [Link] " , " [Link] " ))

tabla $ [Link] <- factor ( tabla $ color , levels = c ( " 5 " , " 4 " , " 3 " , " 2 " ))
tabla $ [Link] <- factor ( tabla $ spine , levels = c ( " 3 " , " 2 " , " 1 " ))

# Se ajusta el modelo con todos los predictores

# El peso se divide por 1000 para evitar overflows

# Nota : I () indica que se toma la operacion aritmetica

# y se evita el sentido de la formula simbolica

c a n g r e . f i t . l o g i s t . t o d o <- glm ( [Link] ∼ [Link] + [Link] + width

+ I ( weight / 1000) , family = binomial , data = tabla )

summary ( [Link] , cor = T )

# Existe alta correlacion entre peso y anchura del caparazon.

# Por ello no se considera el peso en lo que sigue.

c a n g r e . f i t . l o g i s t . o p t i m o <- glm ( [Link] ∼ [Link] * [Link] * width ,

family = binomial , data = tabla )

res = step ( [Link] , list ( lower = ∼ 1 ,

upper = formula ( c a n g r e . f i t . l o g i s t . o p t i m o )) , scale =1 , trace =F ,
direction = " backward " )
res $ anova
res

34
 # Se repite el analisis con el comando stepAIC de la libreria MASS
library ( MASS )
stepAIC ( c a n g r e . f i t . l o g i s t . o p t i m o )

Call :
glm ( formula = [Link] ∼ [Link] + [Link] + width + I ( weight / 1000) ,
family = binomial , data = tabla )

Deviance Residuals :
Min 1Q Median 3Q Max
-2 .1977 -0 .9424 0 .4849 0 .8491 2 .1198

Coefficients :
Estimate Std. Error z value Pr ( >| z |)
( Intercept ) -9 .2734 3 .8378 -2 .416 0 .01568 *
C.fac4 1 .1198 0 .5933 1 .887 0 .05910 .
C.fac3 1 .5058 0 .5667 2 .657 0 .00788 * *
C.fac2 1 .6087 0 .9355 1 .720 0 .08552 .
S.fac2 -0 .4963 0 .6292 -0 .789 0 .43024
S.fac1 -0 .4003 0 .5027 -0 .796 0 .42588
width 0 .2631 0 .1953 1 .347 0 .17788
I ( weight / 1000) 0 .8258 0 .7038 1 .173 0 .24069
---
Signif. codes : 0 ’* * * ’ 0 .001 ’* * ’ 0 .01 ’* ’ 0 .05 ’. ’ 0 .1 ’ ’ 1

( Dispersion parameter for binomial family taken to be 1)

Null deviance : 225 .76 on 172 degrees of freedom

Residual deviance : 185 .20 on 165 degrees of freedom
AIC : 201 .2

Number of Fisher Scoring iterations : 4

Correlation of Coefficients :
( Intercept ) C.fac4 C.fac3 C.fac2 S.fac2 S.fac1 width
C.fac4 -0 .13
C.fac3 -0 .07 0 .72
C.fac2 0 .00 0 .45 0 .55
S.fac2 -0 .22 -0 .07 -0 .17 -0 .25
S.fac1 -0 .01 -0 .03 -0 .21 -0 .37 0 .24
width -0 .96 0 .02 -0 .03 -0 .03 0 .19 0 .02
I ( weight / 1000) 0 .67 -0 .01 0 .00 -0 .04 -0 .09 -0 .04 -0 .83

Se observa que en la matriz de correlaciones asintótica existe una alta correlación

entre el peso y la anchura del caparazón.

35
 Los resultados que se obtienen con la fución step y con la función stepAIC de la
librerı́a MASS son:

Step Df Deviance Resid. Df Resid. Dev AIC

1 NA NA 152 170 .4404 212 .4404
2 - [Link] : [Link] : width 3 3 .23334337 155 173 .6738 209 .6738
3 - [Link] : [Link] 6 7 .88506583 161 181 .5588 205 .5588
4 - [Link] : width 2 0 .07819029 163 181 .6370 201 .6370
5 - [Link] 2 1 .44361451 165 183 .0806 199 .0806
6 - [Link] : width 3 4 .37640503 168 187 .4570 197 .4570

---------------------------------------------------------------------

Step : AIC =197 .46

[Link] ∼ [Link] + width

Df Deviance AIC
< none > 187 .46 197 .46
- [Link] 3 194 .45 198 .45
- width 1 212 .06 220 .06

Call : glm ( formula = [Link] ∼ [Link] + width , family = binomial , data = tabla )

Coefficients :
( Intercept ) C.fac4 C.fac3 C.fac2 width
-12 .715 1 .106 1 .402 1 .330 0 .468

Degrees of Freedom : 172 Total ( i.e. Null ) ; 168 Residual

Null Deviance : 225 .8
Residual Deviance : 187 .5 AIC : 197 .5

Por tanto, en el modelo final solo se incluyen las variables anchura y colores del
caparazón.

36
Curvas ROC: resumen de la potencia predictiva
Para comprobar la efectividad de un modelo en la clasificación de observaciones,
se puede construir una tabla de clasificación donde se cruza el verdadero valor de la
observación (1 ó 0), con la predicción de la misma según el modelo que se considera.
La predicción se suele hacer con respecto a un valor de referencia arbitrario π0 :

ybi = 1 si π
bi > π0

mientras que
bi ≤ π0
ybi = 0 si π

El valor habitual que se suele tomar como valor de corte es π0 = 0,5

Capacidad Predictiva de los modelos Se resume la capacidad predictiva de un

modelo de regresión logı́stica mediante el concepto de sensibilidad:

P (b
y = 1|y = 1)

y mediante el concepto de especificidad:

P (b
y = 0|y = 0)

Es decir, la predicción de éxito cuando es cierto se denomina sensibilidad y la predicción

de un fracaso cuando es, a su vez cierto, se denomina especificidad.
Todo esto es muy sensible a las frecuencias relativas de y = 1 e y = 0.

Curvas ROC

Una curva de tipo receiver operating characteristic (ROC) es un gráfico en el

que se representa la sensibilidad en función de (1 – especificidad).
Si vamos modificando los valores del valor de corte π0 y representamos la sensibilidad
(en ordenadas) frente a (1 – especificidad) (en abscisas) tenemos la curva ROC. Es una
curva cóncava que conecta los puntos (0, 0) y (1, 1). Cuanto mayor sea el área bajo la
curva mejores serán las predicciones.
La curva ROC ofrece un mejor resumen de la capacidad predictiva que una tabla de
clasificación, porque presenta la potencia predictiva para todos los posibles valores de
referencia π0 .
Cuando π0 está cerca de 0 casi todas las predicciones serán yb = 1, con lo cual
la sensibilidad estará próxima a 1 y la especificidad estará cerca de 0. Ası́, el punto
(1 – especificidad, sensibilidad) tendrá coordenadas (1, 1)

37
 Cuando π0 está cerca de 1 casi todas las predicciones serán yb = 0, con lo cual
la sensibilidad estará próxima a 0 y la especificidad estará cerca de 1. Ası́, el punto
(1 – especificidad, sensibilidad) tendrá coordenadas (0, 0). Para una especificidad dada
(fijando un valor en el eje de abscisas), la mayor potencia predictiva corresponde a la
sensibilidad más alta (mayor valor en el eje de ordenadas), de modo que cuanto mayor
sea el área bajo la curva ROC mayor será la potencia de predicción.

En R existen numerosas librerı́as que consieran curvas ROC. Por ejemplo,

library ( verification )
# Se consideran valores 0 -1 frente a probabilidades pi _ i
suceso = c (0 ,0 ,0 ,1 ,1 ,1 ,0 ,1 ,1 ,0 ,0 ,0 ,0 ,1 ,1)
p = c (0 .928 ,0 .576 ,0 .008 ,0 .944 ,0 .832 ,0 .816 ,0 .136 ,
0 .584 ,0 .032 ,0 .016 ,0 .28 ,0 .024 ,0 ,0 .984 ,0 .952 )

A = [Link] ( suceso , p )
[Link] ( A $ suceso , A $ p ) # Dibujo la ROC
[Link] ( A $ suceso , A $ p ) # Area bajo ROC ( Ho : A =0 .5 )

Se obtiene

38
 A
[1] 0 .875
[Link]
[1] 15
[Link]
[1] 7
[Link]
[1] 8
[Link]
[1] 0 .006993007

Se pueden dibujar varias curvas ROC de distintos modelos para comparar su potencia
de predicción. Esta potencia se mide en términos del área debajo de la curva. Un área
igual a 0.5 representa al peor modelo y un área igual a 1 al mejor.

39
40
Ejemplo de curva ROC con SAS

En el ejemplo, se dibuja una curva ROC, se estiman los parámetros y se presentan

las medidas de bondad de ajuste en un modelo de regresión logı́stica.
Los datos consisten en tres variables: n (número de personas en una muestra), en-
fermo (número de enfermos de la muestra) y edad. Se considera un modelo de regresión
logı́stica para estudiar el efecto de la edad sobre la probabilidad de contraer una enfer-
medad.

OPTIONS nodate ls =65 formchar = ’ | - - - -|+| - - -+=| -/\ < >* ’;

x ’ cd " c :\ deSAS " ’;

DATA cosas ;
INPUT enfermo n edad ;
DATALINES ;
0 14 25
0 20 35
0 19 45
7 18 55
6 12 65
17 17 75
;
PROC logistic data = cosas ;
model enfermo / n = edad / scale = none
clparm = wald clodds = pl rsquare outroc = roc1 ;
units edad =10;
RUN ;

La opción scale=none se pone para calcular los estadı́sticos de ajuste sin corregir
para overdispersion, La opción rsquare se usa para calcular los estadı́sticos de tipo R2
generalizados. La opción clparm=wald muestra los intervalos para los parámetros de la
regresión logı́stica. La opción units muestra los estimadores de la razón de odds para
el cambio de 10 años en la variable edad. La opción clodds=pl se pone para obtener los
intervalos de confianza para las razones de odds. La opción outroc=roc1 muestra una
curva ROC.
El área que hay bajo la curva ROC aparece en el estadı́stico c en la tabla deno-
minada Association of Predicted Probabilities and Observed Responses (Aso-
ciación de probabilidades predichas y respuestas observadas). En este ejemplo, el área
bajo la curva ROC es igual a 0.953.

41
symbol1 i = join v = none c = lightblue ;
goptions device = png gsfname = graphout ;
filename graphout " ROC _ SAS . png ";
PROC gplot data = roc1 ;
title ’ Curva ROC ’;
plot _ sensit _ * _ 1 mspec _ =1 / vaxis =0 to 1 by .1
cframe = lightgrey ;
RUN ;

Analysis of Maximum Likelihood Estimates

Standard Wald
Parameter DF Estimate Error Chi - Square Pr > ChiSq

Intercept 1 -12 .5016 2 .5555 23 .9317 < .0001

edad 1 0 .2066 0 .0428 23 .3475 < .0001

Odds Ratio Estimates

Point 95 % Wald
Effect Estimate Confidence Limits

edad 1 .229 1 .131 1 .337

Association of Predicted Probabilities and Observed Responses

Percent Concordant 92 .6 Somers ’ D 0 .906

Percent Discordant 2 .0 Gamma 0 .958
Percent Tied 5 .4 Tau - a 0 .384
Pairs 2100 c 0 .953

Wald Confidence Interval for Parameters

Parameter Estimate 95 % Confidence Limits

Intercept -12 .5016 -17 .5104 -7 .4929

edad 0 .2066 0 .1228 0 .2904

Profile Likelihood Confidence Interval for Adjusted Odds Ratios

Effect Unit Estimate 95 % Confidence Limits

edad 10 .0000 7 .892 3 .881 21 .406

42
43