¡Descarga Regresión Logística Multinomial y más Apuntes en PDF de Ciencia Política solo en Docsity!
Tema 5: Regresi´on Log´ıstica
Multinomial
Introducci´on
En este tema se considera un modelo de regresi´on log´ıstica donde la variable depen-
diente tiene m´as de dos categor´ıas. La respuesta puede o bien ser nominal o bien ordinal.
A su vez, las variables explicativas pueden ser categ´oricas o cuantitativas.
En los modelos de regresi´on multinomial se asume que los recuentos de las categor´ıas de
Y tienen una distribuci´on multinomial. Esta distribuci´on es, a su vez, una generalizaci´on
de la distribuci´on binomial.
Ver:
es.wikipedia.org/wiki/Distribuci´on multinomial
Modelos Logit para respuestas nominales
El n´umero de categor´ıas de Y se denota como J donde { π 1 ,... , πJ } son las probabili-
dades de las distintas categor´ıas tal que
j πj = 1.
Se parte de n observaciones independientes que se localizan en las J categor´ıas. La
distribuci´on de probabilidad del n´umero de observaciones de las J categor´ıas sigue una dis-
tribuci´on multinomial. Esta modeliza la probabilidad de cada una de las posibles maneras
en que n observaciones pueden repartirse entre las J categor´ıas.
Al ser la escala de medida nominal, el orden entre las categor´ıas es irrelevante.
Se toma una categor´ıa como respuesta base , por ejemplo la ´ultima categor´ıa ( J ), y se
define un modelo logit con respecto a ella:
log
πj
πJ
= αj + βj x
donde j = 1 ,... , J − 1.
El modelo tiene J − 1 ecuaciones con sus propios par´ametros, y los efectos var´ıan con
respecto a la categor´ıa que se ha tomado como base.
Cuando J = 2 , el modelo equivale a una ´unica ecuaci´on log ( π 1 /π 2 ) = logit( π 1 ) y se
obtiene el modelo de regresi´on log´ıstica est´andar.
La ecuaci´on general logit con respecto a la categor´ıa base J determina tambi´en los
logits para cualquier pareja de categor´ıas. As´ı, si consideramos dos categor´ıas cualesquiera
a y b ,
log
πa
πb
= log
πa/πJ
πb/πJ
= log
πa
πJ
− log
πb
πJ
= ( αa + βax ) − ( αb + βbx ) = ( αa − αb ) + ( βa − βb ) x
De este modo, la ecuaci´on para las categor´ıas a y b tiene tambi´en la forma α + βx donde
α = ( αa − αb ) y β = ( βa − βb ).
Ejemplo
Se trata de calcular las primas de seguros para los lugare˜nos de una localidad de Florida
(USA). Para ello se analizan los datos sobre los gustos alimentarios de 59 caimanes que
viven por la zona.
La variable respuesta es la elecci´on del tipo de comida principal que se encuentra en
el est´omago de los animales, con cinco posibles categor´ıas: peces , invertebrados , reptiles ,
p´ajaros y otros. Como variables predictoras se consideran:
i ) El lago en donde se captura al animal ( L );
ii ) Su g´enero ( G );
iii ) el tama˜no S (si es menor o igual a 2.3 metros o mayor que este valor).
library ( MASS )
library ( nnet )
comen.f = factor ( c ( " peces " ," invert " ," reptiles " ," pajaritos " ," otros " ) ,
levels = c ( " peces " ," invert " ," reptiles " ," pajaritos " ," otros " ))
tamano.f = factor ( c ( " <2 .3 " , " >2 .3 " ) , levels = c ( " >2 .3 " , " <2 .3 " ))
sexo.f = factor ( c ( " m " ," f " ) , levels = c ( " m " ," f " ))
Coefficients : ( Intercept ) tamano <2. invert -1 .034070 0. reptiles -1 .241705 -0. pajaritos -1 .727214 -0. otros -1 .241709 0.
Std. Errors : ( Intercept ) tamano <2. invert 0 .2910708 0. reptiles 0 .3148729 0. pajaritos 0 .3836949 0. otros 0 .3148735 0.
Residual Deviance : 589. AIC : 605.
La categor´ıa base que se toma es justamente la primera: peces.
Se pueden ajustar distintos modelos logit con la categor´ıa base en peces y compa-
rarlos en t´erminos de sus diferencias en deviance , es decir mediante las diferencias en la
verosimilitud: −2(log( L ( ˆ μ 1 , y )) − log( L ( ˆ μ 2 , y ))).
fitS = multinom ( comen ∼ charca * tamano * sexo , data = tabla.array ,
weights = temp )
fitS2 = multinom ( comen ∼ charca * tamano , data = tabla.array ,
weights = temp )
fit0 = multinom ( comen ∼ 1 , data = tabla.array , weights = temp )
fit1 = multinom ( comen ∼ charca , data = tabla.array , weights = temp )
fit2 = multinom ( comen ∼ tamano , data = tabla.array , weights = temp )
fit3 = multinom ( comen ∼ sexo , data = tabla.array , weights = temp )
fit4 = multinom ( comen ∼ charca + tamano , data = tabla.array ,
weights = temp )
fit5 = multinom ( comen ∼ charca + tamano + sexo , data = tabla.array ,
weights = temp )
# Se compara con el modelo que NO considera sexo
deviance ( fitS2 ) - deviance ( fitS )
[1] 35.
# Se usa un ANOVA con el modelo que NO considera sexo
anova ( fitS2 , fitS , test = ’ Chisq ’)
Likelihood ratio tests of Multinomial Models
Response : comen Model Resid. df Resid. Dev Test Df LR stat. Pr ( Chi ) 1 charca * tamano 288 523. 2 charca * tamano * sexo 256 487 .6018 1 vs 2 32 35 .39865 0.
deviance ( fit0 ) - deviance ( fitS )
[1] 116.
deviance ( fit1 ) - deviance ( fitS )
[1] 73.
deviance ( fit2 ) - deviance ( fitS )
[1] 101.
deviance ( fit3 ) - deviance ( fitS )
[1] 114.
deviance ( fit4 ) - deviance ( fitS )
[1] 52.
deviance ( fit5 ) - deviance ( fitS )
[1] 50.
# Se usa como modelo el que NO considera sexo
deviance ( fit4 ) - deviance ( fitS2 )
anova ( fit4 , fitS2 , test = ’ Chisq ’)
[1] 17.
Likelihood ratio tests of Multinomial Models
Response : comen Model Resid. df Resid. Dev Test Df LR stat. Pr ( Chi ) 1 charca + tamano 300 540. 2 charca * tamano 288 523 .0005 1 vs 2 12 17 .07983 0.
Se trabaja con el modelo que utiliza como covariables la charca y el tama˜no. Se obtienen
los coeficientes de los odds en los que se compara cada tipo de comida con la categor´ıa
base peces.
Tambi´en se pueden estimar las probabilidades de cada una de las categor´ıas para
distintos valores de las variables explicativas.
predictions = predict ( fit4 , type = " probs " ,
newdata = expand.grid ( tamano = tamano.f , charca = charca.f ))
cbind ( expand.grid ( tamano = tamano.f , charca = charca.f ) , predictions )
tamano charca peces invert reptiles pajaritos otros 1 <2 .3 hancock 0 .5352844 0 .09311221 0 .04745855 0 .070402770 0. 2 >2 .3 hancock 0 .5701841 0 .02307664 0 .07182898 0 .140896661 0. 3 <2 .3 oklawaha 0 .2581899 0 .60188004 0 .07723295 0 .008820525 0. 4 >2 .3 oklawaha 0 .4584248 0 .24864187 0 .19484367 0 .029424141 0. 5 <2 .3 trafford 0 .1843017 0 .51682297 0 .08877041 0 .035897986 0. 6 >2 .3 trafford 0 .2957471 0 .19296047 0 .20240167 0 .108228505 0. 7 <2 .3 george 0 .4521217 0 .41284675 0 .01156715 0 .029664777 0. 8 >2 .3 george 0 .6574619 0 .13968167 0 .02389991 0 .081046951 0.
Modelos Acumulados para datos ordinales
Cuando las respuestas de la variable categ´orica son ordinales se pueden utilizar
modelos logit acumulados.
La probabilidad acumulada de una variable Y es la probabilidad de que Y sea menor
o igual que un determinado valor j.
As´ı, para una categor´ıa dada j se define la probabilidad acumulada como
P ( Y ≤ j ) = π 1 + · · · + πj
para j = 1 ,... , J.
Las probabilidades acumuladas reflejan el orden entre las categor´ıas:
P ( Y ≤ 1) ≤ P ( Y ≤ 2) ≤ · · · ≤ P ( Y ≤ J ) = 1
Los logits de las probabilidades acumuladas son
logit [ P ( Y ≤ j )] = log
[
P ( Y ≤ j )
1 − P ( Y ≤ j )
]
log
[
π 1 + · · · + πj
πj +1 + · · · + πJ
]
= log
∑^ j
i =
πi
∑ J
i = j +
πi
para j = 1 ,... , J − 1.
El modelo logit de odds proporcionales, cuando se consideran d variables explicativas
x = ( x 1 ,... , xd ) se puede expresar como
logit [ P ( Y ≤ j | x )] = αj +
∑^ d
i =
βixi
para j = 1 ,... , J − 1.
Cada variable independiente tiene un solo coeficiente que no depende del valor j.
La dependencia con respecto al valor j aparece solo en el coeficiente o constante αj.
Teniendo en cuenta que si j ′^ < j entonces
logit [ P ( Y ≤ j ′| x )] ≤ logit [ P ( Y ≤ j | x )]
se tiene que verificar que αj ≤ αj ′ , es decir, los valores de αj aumentan junto con los
valores de j.
Ejemplo: Enfermedad mental
En el libro Categorical Data Analysis (2002) de Agresti (pag. 279) se muestran los
datos de un estudio sobre una enfermedad mental que se trata de relacionar con dos
variables explicativas. La enfermedad mental se resume en una variable categ´orica con los
siguientes niveles: buen estado , s´ıntomas leves , s´ıntomas moderados y enfermedad. Como
variables predictoras tenemos:
x 1 ≡ mide el n´umero de sucesos impactantes en la vida de la persona en los ´ultimos tres
a˜nos (divorcios, fallecimientos, etc.).
x 2 ≡ Estatus socio-econ´omico con niveles 1 ( alto ) y 0 ( bajo ).
La enfermedad mental, como variable respuesta , es un factor que presenta ordenaci´on
entre sus categor´ıas.
dat = read.table (
" http : / / www.biostat.jhsph.edu /∼bcaffo / aglm / files /
mentalestadoData.dat " )
names ( dat ) = c ( " estado " , " estatus " , " sucesos " )
dat $ estado = as.ordered ( dat $ estado )
levels ( dat $ estado ) = c ( " well " , " mild " , " mod " , " imp " )
dat $ estatus = as.factor ( dat $ estatus )
levels ( dat $ estatus ) = c ( " low " , " high " )
library ( MASS )
fit = polr ( estado ∼ estatus + sucesos , data = dat )
summary ( fit )
Call : polr ( formula = estado ∼ estatus + sucesos , data = dat )
Coefficients : Value Std. Error t value estatushigh -1 .1112 0 .6109 -1. sucesos 0 .3189 0 .1210 2.
Intercepts : Value Std. Error t value well | mild -0 .2819 0 .6423 -0. mild | mod 1 .2128 0 .6607 1. mod | imp 2 .2094 0 .7210 3.
Residual Deviance : 99. AIC : 109.
Se calculan ahora intervalos de confianza para los coeficientes originales y su expresi´on
en t´erminos de la exponencial.
# Intervalos de confianza
( ci = confint ( fit ))
estatushigh -2 .34715412 0. sucesos 0 .09203235 0.
# Exponencial de los parametros e IC
exp ( cbind ( OR = coef ( fit ) , ci ))
OR 2 .5 % 97 .5 %
estatushigh 0 .3291535 0 .09564096 1. sucesos 1 .3755605 1 .09640030 1.
# Seleccionas los de estatus 1 y mas de 5 sucesos
que = subset ( dat , estatus == " high " & sucesos >5)
# Prediccion para el grupo anterior
cbind ( que , predict ( fit , data.frame ( que ) , type = " probs " ))
estado estatus sucesos well mild mod imp 2 well high 9 0 .1150236 0 .2518325 0 .2439856 0. 10 well high 7 0 .1973878 0 .3255946 0 .2251307 0. 17 mild high 8 0 .1516700 0 .2918553 0 .2399334 0. 21 mild high 9 0 .1150236 0 .2518325 0 .2439856 0. 29 mod high 6 0 .2527800 0 .3485130 0 .2020681 0. 32 imp high 8 0 .1516700 0 .2918553 0 .2399334 0. 34 imp high 7 0 .1973878 0 .3255946 0 .2251307 0. 38 imp high 8 0 .1516700 0 .2918553 0 .2399334 0.
Otra opci´on es usar la funci´on vglm. Esta funci´on est´a en la librer´ıa VGAM.
library ( VGAM )
fit = vglm ( estado ∼ estatus + sucesos ,
family = cumulative ( parallel = TRUE ) , data = dat )
summary ( fit )
O en SAS est´andar:
OPTIONS nodate ls =65 formchar = ’ | - - - -|+| - - -+=| -/\ < > * ’; x ’ cd " c :\ dondeSea " ’;
DATA mental ;
INFILE ’ DatosMental. txt ’; INPUT estado estatus sucesos ;
PROC genmod ;
model estado = sucesos estatus / dist = multinomial link = clogit lrci type3 ; RUN ;
Ejemplo
Se tiene una muestra de 735 personas a los que se pregunta por sus preferencias en
cuanto a tres marcas ( brands ) de algunos productos. Se considera adem´as el g´enero y la
edad de las personas de la encuesta.
Los datos (en formato SAS) se pueden descargar en
http://www.ats.ucla.edu/stat/sas/dae/mlogit.sas7bdat
La variable dependiente es brand. La variable female se codifica como 0 para hombres
y 1 para mujeres.
library ( sas7bdat )
lacosa = read.sas7bdat (
" http : / / www.ats.ucla.edu / stat / sas / dae / mlogit.sas7bdat " )
colnames ( lacosa ) = c ( " brand " ," female " ," age " )
attach ( lacosa )
# Se considera un analisis descriptivo
table ( brand )
table ( female )
female 0 1 269 466
summary ( age )
Min. 1 st Qu. Median Mean 3 rd Qu. Max. 24 .0 32 .0 32 .0 32 .9 34 .0 38.
xtabs (∼ female + brand )
brand female 1 2 3 0 92 99 78 1 115 208 143
En los resultados anteriores se obtienen los coeficientes y sus p-valores. Cada variable
independiente y el intercept aparecen dos veces, una para cada una de las categor´ıas
denominadas alt2 y alt3. Corresponden a dos ecuaciones:
log
P (brand = 2)
P (brand = 1)
= b 10 + b 11 · female + b 12 · age
log
P (brand = 3)
P (brand = 1)
= b 20 + b 21 · female + b 22 · age
donde los bij son los coeficientes de regresi´on de la salida.
Por ejemplo, por cada aumento en una unidad de la variable age, el logaritmo del
ratio de las dos probabilidades, P (brand = 2) /P (brand = 1), se incrementa en 0.368, y el
logaritmo del ratio de las dos probabilidades, P (brand = 3) /P (brand = 1), se incrementa
en 0.686. Por tanto, en general, cuanto mayor sea una persona tendr´a m´as preferencia por
brand igual a 2 ´o a 3, que por brand igual a 1.
A continuaci´on, se muestran los resultados del modelo en t´erminos de las probabilida-
des. Por ejemplo, se muestra un rango de distintas edades y se calculan las probabilidades
de escoger cada categor´ıa de brand para mujeres y hombres. Se generan los valores pre-
dichos en la escala logit usando los coeficientes del modelo. En brand = 1, el valor se fija
en 0.
Las columnas etiquetadas como pred.1, pred.2, y pred.3, contienen las probabilida-
des predichas de que brand sea igual a 1, 2 y 3 respectivamente.
newdata = data.frame ( cbind ( age = rep (24:38 , 2) ,
female = c ( rep (0 , 15) , rep (1 , 15))))
logit1 = rep (0 , 30)
logit2 = mlogit.model $ coefficients [[1]] +
mlogit.model $ coefficients [[3]] * newdata $ female +
mlogit.model $ coefficients [[5]] * newdata $ age
logit3 = mlogit.model $ coefficients [[2]] +
mlogit.model $ coefficients [[4]] * newdata $ female +
mlogit.model $ coefficients [[6]] * newdata $ age
logits = cbind ( logit1 , logit2 , logit3 )
p.unscaled = exp ( logits )
p = cbind ( newdata , ( p.unscaled / rowSums ( p.unscaled )))
colnames ( p ) = c ( " age " , " female " , " pred.1 " , " pred.2 " , " pred.3 " )
print ( p )
Se pueden presentar gr´aficamente las predicciones. Por ejemplo, se dibujan las proba-
bilidades predichas de que una persona seleccione brand = 1 en funci´on de la variable age, cuando es hombre (female = 0, en azul ), y cuando es mujer (female = 1, en rojo ).
plot ( p $ pred.1 [ p $ female ==0]∼p $ age [ p $ female ==0] , type = " l " , col = " blue " ,
lwd =1 , ylab = " Probabilidad predicha para Brand = 1 " , xlab = " Edad " ) lines ( p $ pred.1 [ p $ female ==1]∼p $ age [ p $ female ==1] , col = " red " , lwd =1)
- age female pred.1 pred.2 pred.
- 1 24 0 0 .94795544 0 .05023196 0.
- 2 25 0 0 .92560550 0 .07088034 0.
- 3 26 0 0 .89429225 0 .09896621 0.
- 4 27 0 0 .85114170 0 .13611844 0.
- 5 28 0 0 .79312711 0 .18330129 0.
- 6 29 0 0 .71787603 0 .23976170 0.
- 7 30 0 0 .62506819 0 .30169310 0.
- 8 31 0 0 .51809475 0 .36137225 0.
- 9 32 0 0 .40487183 0 .40810394 0.
- 10 33 0 0 .29639567 0 .43175036 0.
- 11 34 0 0 .20299473 0 .42732003 0.
- 12 35 0 0 .13058003 0 .39723991 0.
- 13 36 0 0 .07951592 0 .34957297 0.
- 14 37 0 0 .04627661 0 .29400375 0.
- 15 38 0 0 .02598254 0 .23855066 0.
- 16 24 1 0 .91531696 0 .08189411 0.
- 17 25 1 0 .88078797 0 .11388332 0.
- 18 26 1 0 .83412301 0 .15585707 0.
- 19 27 1 0 .77287121 0 .20869457 0.
- 20 28 1 0 .69561310 0 .27144346 0.
- 21 29 1 0 .60315223 0 .34013099 0.
- 22 30 1 0 .49958719 0 .40713475 0.
- 23 31 1 0 .39239931 0 .46212850 0.
- 24 32 1 0 .29086432 0 .49503122 0.
- 25 33 1 0 .20320727 0 .49979166 0.
- 26 34 1 0 .13411370 0 .47668396 0.
- 27 35 1 0 .08404322 0 .43168573 0.
- 28 36 1 0 .05034227 0 .37368477 0.
- 29 37 1 0 .02903251 0 .31143314 0.
- 30 38 1 0 .01623160 0 .25162236 0.
Se usa ahora SAS Unversity para tratar el mismo conjunto de datos.
OPTIONS nodate ls =75;
/ * O D S l i s t i n g f i l e = ’/ f o l d e r s / m y f o l d e r s / s a l e. lst ’ ; * /
DATA market ; INFILE ’/ folders / myfolders / market. txt ’;
INPUT brand female age ;
PROC freq data = market ; tables brand ;
PROC sort data = market ; by brand ;
PROC means data = market ; by brand ; var age female ;
PROC logistic data = market ; class brand ( ref = "1") ; model brand = female age / link = glogit ; RUN ;
/ * O D S l i s t i n g c l o s e ; * /
O en SAS est´andar
OPTIONS nodate ls =65 formchar = ’ | - - - -|+| - - -+=| -/\ < > * ’; x ’ cd " c :\ deSAS " ’;
DATA market ; INFILE ’ market. txt ’; INPUT brand female age ;
... Igual que el anterior programa ...
