#Introducción

La siguiente práctica tiene el objetivo de introducir a los estudiantes a los modelos de regresión logística. Para consideraciones de este análisis vamos a utilizar la base de datos Titanic.

El RMS Titanic fue un transatlántico británico, el mayor barco del mundo al finalizar su construcción, que se hundió en la madrugada del 15 de abril de 1912 durante su viaje inaugural desde Southampton a Nueva York. En el hundimiento del Titanic murieron 619 personas de las 1046 que iban a bordo, lo que convierte a esta tragedia en uno de los mayores naufragios de la historia ocurridos en tiempo de paz.

La base de datos la pueden descargar acá

Nombre	Descripción	Valores	Categorías
pclass	Clase del pasajero/a	1	Clase Alta
		2	Clase Intermedia
		3	Clase Baja
survived	Estatus sobrevivencia	0	No sobrevive
		1	Sobrevive
name	Nombre
sex	Sexo pasajero/a	0	Hombre
		1	Mujer
age	Edad pasajero/a		rango: 0.2-80.0
sibsp	Número de hermanos /pareja a bordo		rango: 0-8
parch	Número de padres / niños a bordo		rango: 0-9

#Preambulo

En primer lugar cargamos los paquetes.

pacman::p_load(tidyverse,data.table,sjmisc,sjlabelled,
               dplyr,texreg,stargazer)

En segundo lugar cargamos la base de datos.

#recordar definir en base al directorio de trabajo o working directory
load("titanic.Rdata")
names(tt)

## [1] "pclass"   "survived" "sex"      "age"      "sibsp"    "parch"

1 Concepto de razón o ratio

En primer lugar revisamos los estadísticos descriptivos de las variables de interés.

table(tt$survived)

No sobrevive Sobrevive 619 427

table(tt$sex)

Hombre Mujer 658 388

#Porcentajes de sobrevientes
prop.table(table(tt$survived))

No sobrevive Sobrevive 0.5917782 0.4082218

El resultado muestra que 427 personas sobrevivieron al hundimiento del Titanic. En terminos de porcentajes esto expresa el 40.8% de los datos(427/(427+619)=0.4082218) en contraste con el 59.1% de los pasajeros que no sobrevivieron (619/(427+619)=0.5917782).

Esto se puede expresar de la misma forma en términos de razón.

427/619

## [1] 0.6898223

Existen 0.69 personas que sobreviven por cada persona que no sobrevive.

O alternativamente, por cada 100 no sobrevivientes existen 69 sobrevivientes.

##Tambien se puede hacer la interpretación al revés
619/427

## [1] 1.449649

#1.119109 existen 1.12 personas que no tienen disposición de ir a votar por cada una que esta dispuesta.

Existen 1.45 personas que no sobreviven por cada persona que sobrevive.

o alternativamente, por cada 100 sobrevivientes, existen 145 no sobrevivientes al hundimiento del Titanic.

La razón o ratio es el coenciente entre dos cantidades y señala cuantas veces una cantidad es mayor o menor respecto a la otra

2 Concepto de Odds

addmargins(table(tt$survived,tt$sex))

##               
##                Hombre Mujer  Sum
##   No sobrevive    523    96  619
##   Sobrevive       135   292  427
##   Sum             658   388 1046

round(prop.table(table(tt$survived,tt$sex),2),2) #perfil columna se le indica a prop.table con el segundo argumento "2"

##               
##                Hombre Mujer
##   No sobrevive   0.79  0.25
##   Sobrevive      0.21  0.75

Según esto el 21% de los hombres sobrevivieron al hundimiento y el 75% de las mujeres sobrevivieron.

Tambien se puede expresar en términos de razón:

la razón de los hombres que tienen de sobrevivir se puede expresar como:

\[Odds_{hombres}=\frac{p}{(p-1)}=\frac{0.21}{0.79}=0.27\] y en R

0.21/0.79

## [1] 0.2658228

Hay 0.27 hombres que sobreviven al titanic por cada uno que no sobrevive

o en otros términos

Hay 27 hombres que sobreviven al titanica por cada 100 hombres que no sobreviven

y para mujeres

\[Odds_{mujeres}=\frac{p}{(p-1)}=\frac{0.75}{0.25}=3\]

0.75/0.25

## [1] 3

Hay 3 mujeres sobrevivientes por cada mujere que no sobrevive al hundimiento del titanic

El análisis que acabamos de realizar es un analisis de Odds

\[0dd =\frac{p}{q}=\frac{p}{1-p}\]

el término odd refiere a la razón que se establece entre la ocurriencia de un evento (o su probabilidad) respecto al suceso de su no ocurrencia.

3 Concepto de Odds Ratio

La comparación de los Odds de dos grupos es conocido como Odds Ratio (OR). Esto equivale a

\[\theta=\frac{odds_{1}}{odds_{2}}=\frac{\pi_{1}/(1-\pi_{1})}{\pi_{2}/(1-\pi_{2})}\] Propiedades:

Cuando X e Y son independientes \(\theta=1\) ya que \(odds_{1}=odds_{2}\)
El rango de posibles valores es: \(0<\theta<\infty\)

-Cuando los valores van de 0 a 1, \(\theta\) indica que \(odds_{1}<odds_{2}\)

-Cuando los valores van de 1 a \(\infty\), \(\theta\) indica que \(odds_{1}>odds_{2}\)
Es una medida de magnitud de asociación simétrica: un \(\theta=4\) es una asociació positiva proporcional a la asociación negativa \(\theta=1/4=0.25\)

A través de los ratios se puede responder a la pregunta ¿Cuánto más probable es que las mujeres sobrevivan al titanic en relación a los hombres?**

prop.table(table(tt$survived,tt$sex),2)

##               
##                   Hombre     Mujer
##   No sobrevive 0.7948328 0.2474227
##   Sobrevive    0.2051672 0.7525773

odd_H<-0.2051672/0.7948328
odd_M<-0.7525773/0.2474227

odd_M/odd_H

## [1] 11.78364

La probabilidad de encontrar una mujer que sobreviva al titanic sobre una que no lo hace es de 11.78 veces respecto al caso de los varones.

o alternativamente:

La probabilidad de sobrevivir al titanic entre las mujeres es [(11.78-1)*100]=1078% más altas que entre los hombres

Este término se le llama Odd ratio (o razón de ODD, o razón de momios) y puede interpretarse como una razón de probabilidades. En nuestro caso la probabilidad de encontrar entre a una mujer sobreviviente es mucho más alta que la probabilidad de encontrar a un hombre superviviente.

3.1 Odds y Logit: algunas consideraciones

Los Odds pueden variar de 0 a \(+\infty\). Teniendo en cuenta de que la probabilidad (\(p\) o \(\pi\) según sea la nomenclatura) varia de 0 a 1, por lo tanto cuando p está muy cerca a 1:

\[odd=\frac{p}{1-p}=\frac{1}{0}=+\infty\]

Y en el caso inverso

\[odd=\frac{p}{1-p}=\frac{0}{1}=0\]

A partir del Odd es posible definir el logit, es decir el logaritmo natural del Odd:

En regresión logística se emplea el logaritmo natural o neperiano cuya base es el número *e*, un número irracional cuyo valor aproximado es 2,718. El logaritmo natural de un número x, notado como ln(x), es el exponente al que debe ser elevado el número *e* para obtener x.

El logit puede tomar cualquier valor real entre \(-\infty\) y \(+\infty\), y permite una lectura simétrica de la relación entre proporciones.

Vamos a tomar el caso de los pasajeros del Titcanis. El 63% de los pasajeros corresponden a hombres y el 37% corresponden a mujer.

prop.table(table(tt$sex))

## 
##    Hombre     Mujer 
## 0.6290631 0.3709369

(odd_h<-0.6290631/0.3709369)

## [1] 1.695876

(odd_m<-0.3709369/0.6290631)

## [1] 0.5896656

Se definomos p como la proporción de hombres, el odd será de 1.695876(0.6290631/0.3709369), mientras que si definimos p como la proporción de mujeres, el odd será de 0.5896656(0.3709369/0.6290631)

Cuando p es menor que q, el Odd se moverá entre 0 y 1. mientras que cuando \(p>q\) el odd se moverá entre 1 y \(+\infty\).

Sin embargo, mediante a la transformación logarítmica del odd los valores son comparables sin importar cuál sea la categória tomada para el cálculo de la proporición.

log(odd_h)

## [1] 0.5281996

log(odd_m)

## [1] -0.5281996

En los siguientes gráficos se puede observar la comparación en el comportamiento de las odds y de logit en relación a las distintos valores que puede asumir la probabilidad

#se cargan paquetes más especificos
library(hrbrthemes)
library(plotly)
library(babynames)
library(viridis)

#ajustamos una variable de probabilidades y sus respectivas transformaciones a odds y logit
p<-seq(0,1,by=0.01)
odd<- p/(1-p)
logit<- log(p/(1-p))
dat<- as.data.frame(cbind(p,odd,logit))


dat %>%
  ggplot( aes(x=p, y=odd)) +
    geom_line(color="#69b3a2") +
    ggtitle("odd segun valores de p") +
    ylab("odd") +
    theme_ipsum()

dat %>%
  ggplot( aes(x=p, y=logit)) +
    geom_line(color="#69b3a2") +
    ggtitle("logit segun valores de p") +
    ylab("logit") +
    theme_ipsum()

m0<-lm(as.integer(survived)~sex,data=tt)
stargazer::stargazer(m0, type="text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                        as.integer(survived)    
## -----------------------------------------------
## sexMujer                     0.547***          
##                               (0.027)          
##                                                
## Constant                     1.205***          
##                               (0.016)          
##                                                
## -----------------------------------------------
## Observations                   1,046           
## R2                             0.289           
## Adjusted R2                    0.289           
## Residual Std. Error      0.415 (df = 1044)     
## F Statistic          425.272*** (df = 1; 1044) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

modelo_titanic<-glm(survived~sex,data=tt, family = "binomial")


summary(modelo_titanic)

## 
## Call:
## glm(formula = survived ~ sex, family = "binomial", data = tt)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6713  -0.6777  -0.6777   0.7540   1.7798  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.35431    0.09654  -14.03   <2e-16 ***
## sexMujer     2.46671    0.15219   16.21   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1414.6  on 1045  degrees of freedom
## Residual deviance: 1102.0  on 1044  degrees of freedom
## AIC: 1106
## 
## Number of Fisher Scoring iterations: 4

stargazer::stargazer(modelo_titanic, type="text")

## 
## =============================================
##                       Dependent variable:    
##                   ---------------------------
##                            survived          
## ---------------------------------------------
## sexMujer                   2.467***          
##                             (0.152)          
##                                              
## Constant                   -1.354***         
##                             (0.097)          
##                                              
## ---------------------------------------------
## Observations                 1,046           
## Log Likelihood             -551.004          
## Akaike Inf. Crit.          1,106.008         
## =============================================
## Note:             *p<0.1; **p<0.05; ***p<0.01

stargazer::stargazer(modelo_titanic,apply.coef = exp,
          apply.se   = exp,
          covariate.labels = "Mujer",
          dep.var.caption = "Variable dependiente",
          dep.var.labels = "Sobrevivío",
          type = "text")

## 
## =============================================
##                      Variable dependiente    
##                   ---------------------------
##                           Sobrevivío         
## ---------------------------------------------
## Mujer                      11.784***         
##                             (1.164)          
##                                              
## Constant                     0.258           
##                             (1.101)          
##                                              
## ---------------------------------------------
## Observations                 1,046           
## Log Likelihood             -551.004          
## Akaike Inf. Crit.          1,106.008         
## =============================================
## Note:             *p<0.1; **p<0.05; ***p<0.01

4 Ajuste

4.1 Devianza

Devianza =-2*log likelihood: Se utiliza como una medida de los residuos generados por el modelo, comparando con el modelo nulo (sin predictores). En general si disminuye, el modelo es mejor

modelo_titanic$null.deviance # devianza modelo sin predictores

## [1] 1414.62

modelo_titanic$deviance # devianza modelo con predictores

## [1] 1102.008

4.2 Test de razón de verosimilitud

Compara las verosimilitudes del modelo con otro con menos predictores

null_titanic=glm(survived~1, data=tt, family=binomial)

anova(modelo_titanic, null_titanic, test ="Chisq")

## Analysis of Deviance Table
## 
## Model 1: survived ~ sex
## Model 2: survived ~ 1
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1044     1102.0                          
## 2      1045     1414.6 -1  -312.61 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La diferencia entre los modelos es estadísticamente significativa con una probabilidad \(p<0.01\). Por lo tanto el modelo con predictores (sexo) ofrece un mejor ajuste a los datos que un modelo sin predictores.

4.3 McFadden (pseudo) R2

Se define como: \(1−[LL(LM)/LL(L0)]\), donde

LL es el log likelihood del modelo
LM es el modelo posterior (con predictores)
L0 es el modelo nulo

logLik(modelo_titanic); logLik(null_titanic)

## 'log Lik.' -551.0042 (df=2)

## 'log Lik.' -707.3102 (df=1)

1-(-551/-707)

## [1] 0.2206506

También se puede obtener con la función PseudoR2 de la librería DescTools, junto a otras versiones de pseudo R2s, como “Nagelkerke”, “CoxSnell” y “Effron”.

4.4 Akaike (AIC)

AIC - Akaike information criteria, evalua la calidad del modelo a través de la comparación con otros modelos penalizando por la inclusión de predictores (análogo al R2 ajustado):

\[AIC=-2(log-likelihood)+2K\]

Donde K= número de parámetros del modelo (regresores + intercepto)

logLik(modelo_titanic)

## 'log Lik.' -551.0042 (df=2)

2*551

## [1] 1102

\[AIC=-2(-551)+2(2)=1102+4=1106\]

En sí mismo no es interpretable, solo como criterio comparativo: menor AIC es mejor fit.

Practica 6: Introduccióna Logit

Estadistica Multivariada - Sociología FACSO Universidad de Chile

Juan Carlos Castillo & Alejandro Plaza