Sesión 4
Profesor:Dr. José Manuel MAGALLANES REYES, Ph.D
Profesor Principal del Departamento de Ciencias Sociales, Sección de Ciencia Política y Gobierno.
Oficina 223 - Edificio CISEPA / ECONOMIA / CCSS
Telefono: (51) 1 - 6262000 anexo 4302
Correo Electrónico: jmagallanes@pucp.edu.pe
Modelando Variables Categóricas
1 Presentación
Michael Cowles y Carlone Davis prepararon una data para investigar quien es más proclive a ser voluntario en experimentos sociales (Cowles and Davis 1987). La data y la metadata están disponibles en (Arel-Bundock 2022). Veamos los contenidos:
rm(list = ls()) # clear memory
knitr::knit_hooks$set(inline = as.character) # inline as string
gitLink="https://vincentarelbundock.github.io/Rdatasets/csv/carData/Cowles.csv"
vol=read.csv(gitLink)[,-1] #no first column
# formatting:
vol[,c(3,4)]=lapply(vol[,c(3,4)],as.factor)
# display table
library(magrittr) # needed for pipe %>%
vol%>%
rmarkdown::paged_table()Según la información proporcionada en el repositorio, esta es la metadata:
- neuroticism: Una medida de inestabilidad emocional, los valores están en la escala de Eysenck (Eysenck and Eysenck 1975).
- extraversion: Una medida del interes por el mundo exterior de la gente y de las cosas, los valores están en la escala de Eysenck (Eysenck and Eysenck 1975).
- sex: una categoría con las modalidades: female y male.
- volunteer: Un factor dicotómico que representa ser o no (yes o no) voluntario en experimentos sociales.
Esta investigación hipotetiza que:
El ofrecerse como voluntario a experimentos sociales está relacionado con el sexo de la persona, su nivel de neuroticismo y su nivel de extraversión.
Hasta ahora hemos querido modelar hipotesis donde la variable dependiente eran valores continuos, y de conteos pero, como vemos, en esta hipótesis la variable dependiente es volunteer (ofrecerse como voluntario), la cuál es categórica (dicotómica). Esto no puede ser tratado con las técnicas anteriores (MAGALLANES 2022a, 2022b).
Veamos un resumen univariado en la Tabla 1.1.
library(summarytools)## Warning: package 'summarytools' was built under R version 4.4.1library(kableExtra)
# usemos la función "dfSummary" de paquete "summarytools"
dfSummary(vol,
plain.ascii = FALSE,
varnumbers = FALSE,
style = "grid",
graph.col=F,
na.col = FALSE) %>%
kable(caption = "Descriptivos Univariados")%>%
kable_styling(full_width = F)| Variable | Stats / Values | Freqs (% of Valid) | Valid |
|---|---|---|---|
|
neuroticism [integer] |
|Mean (sd) : 11.5 (4.9) min < med < max: 0 < 11 < 24 IQR (CV) : 7 (0 |
|
|1421 (100 |
|
extraversion [integer] |
|Mean (sd) : 12.4 (3.9) min < med < max: 2 < 13 < 23 IQR (CV) : 5 (0 |
|
|1421 (100 |
|
sex [factor] |
|1. female 2. male |
|\780 (54.9%) \641 (45.1 |
) |1421 (100. |
|
volunteer [factor] |
|1. no 2. yes |
|\824 (58.0%) \597 (42.0 |
) |1421 (100. |
2 Probabilidades y ODDS
Hagamos un breve repaso de cómo usar tablas de contigencia para analizar variables categóricas. En este caso, analicemos la relación entre sexo del entrevistado y ser voluntario. Al ser ambas categóricas, y ser la hipótesis asimétrica o direccional, preparemos una tabla de contingencia que muestre ello. Veamos la Tabla 2.1.
dep=vol$volunteer # a la fila
ind=vol$sex # a la columna
volsexTable=table(dep,ind,dnn = c('volunteer','sex'))
### suma por fila y columna
addmargins(volsexTable)%>%
kable(caption = "Tabla de Contingencia: 'Sexo' y 'Ser Voluntario'")%>%
kableExtra::kable_styling(full_width = F)| female | male | Sum | |
|---|---|---|---|
| no | 431 | 393 | 824 |
| yes | 349 | 248 | 597 |
| Sum | 780 | 641 | 1421 |
Recuérdese que se acostumbra que la variable independiente se sitúe en columnas, y la variable dependiente en filas. A partir de esta información tratemos de explorar la hipótesis:
Ser mujer está relacionado con ser voluntario en proyectos de investigación.
Saber si ser mujer está relacionado con ser voluntario implica responderse:
- Si elijo una mujer al azar, ¿qué probabilidad se tiene que ésta sea voluntaria?.
- ¿Cómo comparo el “voluntarismo” de la mujer con el “voluntarismo” del hombre?
Para ello usaremos los conceptos de probabilidad, y en particular de Odds Ratio (Szumilas 2010). Y todos los cálculos usarán los valores de la Tabla 2.1.
A. Probabilidad y Odds para caso de la mujer:
- La probabilidad (entre 0 y 1) que una mujer sea voluntaria es la división del ser voluntaria entre el total de mujeres:
ProbMujVol=volsexTable[2,1]/sum(volsexTable[,1])
MASS::fractions(ProbMujVol)## [1] 349/780Es decir:
ProbMujVol## [1] 0.4474359- Representemos el odds que sea voluntaria, la división entre ser o no ser voluntaria:
OddsMujVol=volsexTable[2,1]/volsexTable[1,1]
MASS::fractions(OddsMujVol)## [1] 349/431Es decir:
OddsMujVol## [1] 0.8097448El odds suele representarse además como la razón entre dos probabilidades: la probabilidad que SÍ ocurra un evento dividido por la probabilidad que NO ocurra ese evento:
ProbMujVol/(1-ProbMujVol)## [1] 0.8097448B. Probabilidad y Odds para caso del hombre:
- Probabilidad que un hombre sea voluntario
ProbHomVol=volsexTable[2,2]/sum(volsexTable[,2])
MASS::fractions(ProbHomVol)## [1] 248/641O…
ProbHomVol## [1] 0.3868955Igual que antes, calculamos el odds:
OddsHomVol=ProbHomVol/(1-ProbHomVol)
OddsHomVol## [1] 0.6310433C. Comparando Mujer y Hombre:
Hasta aquí ya sabemos cuál odds ha salido mayor. Eso es información clara. Para precisar esa información, podemos dividir los odds, lo que nos da el odds ratio, que en este caso nos dirá qué diferencia produce el sexo en el voluntariado:
(OR_MujHom=OddsMujVol/OddsHomVol)## [1] 1.283184Con ese valor, ya sabemos que el odds de ser mujer supera al odds del hombre por un factor de 0.28. El odds ratio (OR) puede ir de 0 a infinito. Un OR de 1 implica que no hay diferencias.
Veámoslo como fracción:
MASS::fractions(OddsMujVol/OddsHomVol)## [1] 137157/106888Si encuentras la cantidad de voluntarias del numerador, se espera la cantidad de voluntarios del denominador.
Veamos la tabla de contingencia (Tabla 2.1) de manera gráfica en la Figura 2.1.
mosaicplot( t(volsexTable),col = c("orange", "green"),main = "")
Figure 2.1: Tabla de Contingencia como Mosaico
La diferencia es clara, pero ¿será ésta estadísticamente significativa?
3 Regresión Logística Binaria
La regresión logística modela el comportamiento de la probabilidad del evento de interés, según vemos en la Ecuación (3.1).
\[\begin{equation} {log}(\frac{{Y}}{{1 – Y}}) ={log}(\frac{{p(X)}}{{1 – p(X)}}) =\alpha + \beta \cdot X + \epsilon \tag{3.1} \end{equation}\]La parte izquierda de la Ecuación (3.1) representa esa probabilidad pero en términos del odds; y la parte derecha es igual a la ecuación de una recta, pero estamos modelando al logaritmo del odds. Veamos la Tabla 3.1.
### semilla
set.seed(2019)
### first hypothesis
h1=formula(volunteer~sex)
#regression
rlog1=glm(h1, data=vol,family = binomial)
modelrl=list('Ser Voluntario (I)'=rlog1)
#f <- function(x) format(x, digits = 4, scientific = FALSE)
library(modelsummary)
modelsummary(modelrl,
title = "Regresión Logística",
stars = TRUE,
output = "kableExtra")| Ser Voluntario (I) | |
|---|---|
| (Intercept) | −0.211** |
| (0.072) | |
| sexmale | −0.249* |
| (0.108) | |
| Num.Obs. | 1421 |
| AIC | 1932.2 |
| BIC | 1942.7 |
| Log.Lik. | −964.101 |
| F | 5.286 |
| RMSE | 0.49 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
En general, se nota que el sexo (ser hombre) tiene efecto negativo y significativo sobre el logaritmo del odds de ser voluntario. Pero antes de ir a la interpretación detallada veamos otras consideraciones.
3.1 Consideraciones Previas
La regresión es sencilla de calcular, pero la presencia de variables categóricas puede traer complicaciones. Veamos cómo R está utilizando la variable dependiente:
# la representación numérica
as.numeric(as.factor(c("no","yes")))## [1] 1 2Qué pasaría si tuvieramos otras columnas:
vol$IsVolunteer=ifelse(vol$volunteer=="yes",T,F)
vol$IsVolunteer_1=as.integer(ifelse(vol$volunteer=="yes",1,0))
vol$IsVolunteer_claro=as.factor(ifelse(vol$volunteer=="yes",
"claro","nunca"))
# nos queda:
library(magrittr)
vol[,c("volunteer","IsVolunteer","IsVolunteer_1","IsVolunteer_claro")]%>%
rmarkdown::paged_table(options = )Resscribamos la hipótesis inicial con las nuevas variables:
h1b=formula(IsVolunteer~sex)
h1c=formula(IsVolunteer_1~sex)
h1d=formula(IsVolunteer_claro~sex)Veamos:
- Dicotómica como Booleana (lógica)
coef(glm(h1b, data=vol,family = binomial))## (Intercept) sexmale
## -0.2110362 -0.2493447- Dicotómica como Entero:
coef(glm(h1c, data=vol,family = binomial))## (Intercept) sexmale
## -0.2110362 -0.2493447- Dicotómica como otras palabras
coef(glm(h1d, data=vol,family = binomial))## (Intercept) sexmale
## 0.2110362 0.2493447Como se ve, R usa bien los casos Booleanos y de enteros, pero ante la presencia de otros textos en la respuesta sólo usará el orden alfabético.
Otro tema es el uso de categóricas en las variables independientes. El resultado anterior ha usado a mujer como referencia: nos informa cuánto afecta ser hombre al logaritmo del odds de ser voluntario, en comparación con ser mujer (la referencia). La Tabla 3.1 nos confirma que ser hombre disminuye (aquí el signo ayuda) el logaritmo del odds (y por ende la probabilidad) de ser voluntario. Veamos alternativamente la Tabla 3.2 con la misma regresión pero con la categoría hombre como referencia.
# cambio de referencia
vol$sex=relevel(vol$sex,ref = "male")
#regresion
rlog1=glm(h1, data=vol,family = binomial)
modelrl=list('Ser Voluntario (I)'=rlog1)
modelsummary(modelrl,
title = "Regresión Logística (con cambio de referencia en sexo)",
stars = TRUE,
output = "kableExtra")| Ser Voluntario (I) | |
|---|---|
| (Intercept) | −0.460*** |
| (0.081) | |
| sexfemale | 0.249* |
| (0.108) | |
| Num.Obs. | 1421 |
| AIC | 1932.2 |
| BIC | 1942.7 |
| Log.Lik. | −964.101 |
| F | 5.286 |
| RMSE | 0.49 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
La Tabla 3.2, nos muestra el efecto (significativo al 0.05) positivo (nuevamente el signo ayuda) de ser mujer en ser voluntario .
Recuerda que este tema de la categoría de referencia se aplica a todo tipo de regresión.
3.2 Interpretación
Los coeficientes de las Tablas 3.1 y Tabla 3.2, como modelan al logaritmo del odds, no son fáciles de interpretar. Pero, si le aplicamos exponencial obtendrás un odds ratio: veamos el caso del coeficiente 0.249 de la Tabla 3.2:
sexF=coef(rlog1)["sexfemale"]
exp(sexF)## sexfemale
## 1.283184Este valor ya lo conocíamos, pero ahora sabemos que el efecto de sexo es significativo gracias a la regresión logística, algo que no sabíamos con las tablas de contingencia. Lo que es más, sabemos que el hecho de ser mujer eleva en promedio en un factor de 0.28 al odds ratio de ser voluntario que si el entrevistado fuera hombre.
Algo importante ahora es que con la regresión logística ya podemos tener predictores numéricos, lo que escapa de la utilidad de las tablas de contingencia. Así, podemos en el modelo II añadir una numérica a los predictores:
El ofrecerse como voluntario está relacionado con el sexo de la persona, y su nivel de neuroticismo.
La Tabla 3.3 muestra los resultados para el modelo I y el II, ya con los coeficientes exponenciados:
### second hypothesis
h2=formula(volunteer~sex + neuroticism)
rlog2=glm(h2, data=vol,family = binomial)
modelsrl=list('Ser Voluntario (I)'=rlog1,
'Ser Voluntario (II)'=rlog2)
# formato creado para modelsummary
formatoNumero = function(x) format(x, digits = 4, scientific = FALSE)
modelsummary(modelsrl,
fmt=formatoNumero, # usa función que creé antes
exponentiate = T, # coeficientes sin logaritmo
statistic = 'conf.int', # mostrar ICs
title = "Regresión Logísticas (Coeficientes Exponenciados)",
stars = TRUE,
output = "kableExtra")| Ser Voluntario (I) | Ser Voluntario (II) | |
|---|---|---|
| (Intercept) | 0.631*** | 0.6268** |
| [0.5377, 0.7391] | [0.4733, 0.8281] | |
| sexfemale | 1.283* | 1.2817* |
| [1.0378, 1.5878] | [1.0329, 1.5917] | |
| neuroticism | 1.0007 | |
| [0.9789, 1.0228] | ||
| Num.Obs. | 1421 | 1421 |
| AIC | 1932.2 | 1934.2 |
| BIC | 1942.7 | 1950.0 |
| Log.Lik. | −964.101 | −964.099 |
| F | 5.286 | 2.645 |
| RMSE | 0.49 | 0.49 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
En el modelo II notamos que el efecto de ser mujer disminuye, pero no mucho. Probemos ahora el modelo III, el que completa la hipótesis de los investigadores:
El ofrecerse como voluntario está relacionado con el sexo de la persona, su nivel de neuroticismo y su nivel de extraversión.
Veamos lo en simple:
h3=formula(volunteer~sex + neuroticism + extraversion)
rlog3=glm(h3, data=vol,family = binomial)
summary(rlog3)##
## Call:
## glm(formula = h3, family = binomial, data = vol)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.351658 0.239818 -5.636 1.74e-08 ***
## sexfemale 0.235161 0.111185 2.115 0.0344 *
## neuroticism 0.006362 0.011357 0.560 0.5754
## extraversion 0.066325 0.014260 4.651 3.30e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1933.5 on 1420 degrees of freedom
## Residual deviance: 1906.1 on 1417 degrees of freedom
## AIC: 1914.1
##
## Number of Fisher Scoring iterations: 4Ahora, a todos los modelos juntos 3.4:
modelsrl=list('Ser Voluntario (I)'=rlog1,
'Ser Voluntario (II)'=rlog2,
'Ser Voluntario (III)'=rlog3)
f <- function(x) format(x, digits = 4, scientific = FALSE)
modelsummary(modelsrl,
fmt=formatoNumero,
exponentiate = T,
statistic = 'conf.int',
title = "Comparando Regresión Logísticas (Coeficientes Exponenciados)",
stars = TRUE,
gof_map = c("nobs","aic","bic","rmse","logLik"), #comparar
gof_omit = c("F"),
output = "kableExtra")| Ser Voluntario (I) | Ser Voluntario (II) | Ser Voluntario (III) | |
|---|---|---|---|
| (Intercept) | 0.631*** | 0.6268** | 0.2588*** |
| [0.5377, 0.7391] | [0.4733, 0.8281] | [0.1611, 0.4127] | |
| sexfemale | 1.283* | 1.2817* | 1.2651* |
| [1.0378, 1.5878] | [1.0329, 1.5917] | [1.0177, 1.5738] | |
| neuroticism | 1.0007 | 1.0064 | |
| [0.9789, 1.0228] | [0.9842, 1.0291] | ||
| extraversion | 1.0686*** | ||
| [1.0393, 1.0991] | |||
| Num.Obs. | 1421 | 1421 | 1421 |
| AIC | 1932.2 | 1934.2 | 1914.1 |
| BIC | 1942.7 | 1950.0 | 1935.1 |
| RMSE | 0.49 | 0.49 | 0.49 |
| Log.Lik. | −964.101 | −964.099 | −953.031 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Del modelo III tenemos que:
- Tres coeficientes son mayores que uno (podrían tener un efecto positivo, pero falta revisar otra información).
- El intervalo de confianza de la neuroticism incluye al 1.1
- El p-valor de neuroticism es mayor a 0.1 (la probabilidad que sea 1 -no tenga efecto- es mayor al 10%, lo que contribuye a aceptar que no tendría efecto), por lo que lo consideramos no significativo.
El coeficiente de la variable sexo con hombre como referencia es:
round(exp(coef(rlog3)["sexfemale"]),4)## sexfemale
## 1.2651Ese número es el factor por el que se multiplica en promedio el odds ratio de ser voluntario cuando se es mujer, o dicho de otra manera el odds ratio de ser voluntario se eleva en en promedio en ese porcentaje:
100*round(abs(1-exp(coef(rlog3)["sexfemale"])),4)## sexfemale
## 26.51De manera similar, el coeficiente de la variable extraversión es:
round(exp(coef(rlog3)["extraversion"]),4)## extraversion
## 1.0686Ese número es el factor por el que se multiplica en promedio el odds ratio de ser voluntario cada vez que la escala de extroversión aumenta en uno, o dicho de otra manera el odds ratio de ser voluntario se eleva en promedio en este porcentaje cada vez que la escala de extroversión aumenta en uno2:
100*round(abs(1-exp(coef(rlog3)["extraversion"])),4)## extraversion
## 6.86Podemos calcular los coeficientes estandarizados (IBM SPSS 2020) para saber cuál de los predictores del modelo III tiene mayor efecto (ver Tabla 3.5).
vol$sexFem_num=ifelse(vol$sex=='male',0,1) # 1 es para mujer
sdVIs=apply(vol[,c("sexFem_num","neuroticism", "extraversion")],2,sd)
DF=list(LogitSt=sdVIs*coef(rlog3)[c(2,3,4)])%>%
data.frame()
# DF tiene los coeficientes estandarizados
DF%>% kable(caption = "Coeficientes Estandarizados (ordenar vía valores absolutos)")%>%
kableExtra::kable_styling(full_width = F)| LogitSt | |
|---|---|
| sexFem_num | 0.1170580 |
| neuroticism | 0.0311676 |
| extraversion | 0.2582503 |
Usando valores absolutos notamos cuál es el orden de importancia, en este caso extraversion es la que más efecto tiene. Pero no usemos esos valores para interpretar el efecto numérico, sólo saber cuál afecta más que otro.
3.3 Decidiendo mejor modelo
R entrega algunos índices muy para comparar los modelos logísticos, el BIC, AIC, RMSE3, y Log.Lik. Los primeros tres indican que el modelo es mejor si tienen valores bajos; mientras que el último será mejor si su valor aumenta.
Sin embargo, las diferencias pueden ser mínimas entre estos índices, por lo que es apropiado usar el test de razón de verosimilitud (likelihood ratio test -LRT), cuyo resultado lo vemos en la Tabla 3.7.
library(lmtest)
lrtest(rlog1,rlog2, rlog3) %>%
kable(caption = "Tabla LRT para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 2 | -964.1009 | NA | NA | NA |
| 3 | -964.0993 | 1 | 0.0033761 | 0.9536657 |
| 4 | -953.0306 | 1 | 22.1372287 | 0.0000025 |
Como vemos en la Tabla 3.7, pasar del modelo 1 al modelo 2 es no significativo (probabilidad de similitud entre los modelos es mayor que 90%), por lo que no se rechaza la hipotesis nula de igualdad de modelos. Pero sí es significativo pasar del modelo 2 al modelo 3.
Podemos comparar los modelos I, II, y III de manera gráfica en la Figura 3.1.
library(ggplot2)
dotwhisker::dwplot(list(Logit_I=rlog1,Logit_II=rlog2,Logit_III=rlog3),
exp=T) + #exponenciados!
scale_y_discrete(labels=c("extraversion",
"neuroticism",
"sex (ref: male)")) +
scale_color_discrete(name="Modelos para:\nSer Voluntario") +
geom_vline(xintercept = 1,
colour = "grey60",
linetype = 2)
Figure 3.1: Comparación visual de modelos
En la Figura 3.1 se nota claramente por qué el neuroticismo no es significativo.
3.4 Evaluando modelo seleccionado
Cuando usamos una regresión logística queremos predecir probabilidades, en este caso, ante una combinación de los predictores X, qué probabilidad se obtiene con esta regresión para ser voluntario. Esos son los valor ajustados (fitted) que te entrega el modelo, y nos servirá para evaluar si nos hemos acercado al Y original. Calculemos tales probabilidades:
vol$predictedProbs <- predict(rlog3, vol,type = "response")
# veamos algunos de ellos
head(vol$predictedProbs)## [1] 0.4619546 0.4080078 0.4356976 0.5648592 0.4914463 0.4210156Aquí podemos usar estos valores para crear el 0 o 1:
vol$predicted=ifelse(vol$predictedProbs>0.5,"yes","no")Veamos algunas filas de estas dos columnas:
head(vol[,c("predicted","volunteer")],10)## predicted volunteer
## 1 no no
## 2 no no
## 3 no no
## 4 yes no
## 5 no no
## 6 no no
## 7 no no
## 8 no no
## 9 no no
## 10 no noLa columns predicted busca ser la misma que la columna volunteer, pero en algunos casos fallará.
Para saber qué tanto se acertó creamos una Matriz de Confusión (Suresh 2021) usando el paquete cvms 4 en la Figura 3.2.
library(cvms)
ConfMatrix=confusion_matrix(predictions = vol$predicted,
targets= vol$volunteer)
library(cvms)
plot_confusion_matrix(ConfMatrix,
class_order=c("yes","no"))
Figure 3.2: Matriz de confusión
La diagonal secundaria o antidiagonal debería tener sólo 0s si tuviésemos una predicción perfecta, y los demás valores estar en la diagonal principal. El modelo, debería servir para reducir sustantivamente ese error.
Note que en la matriz de confusión se puede calcular valores clave (Trevethan 2017):
Verdaderos Positivos (VP)
Verdaderos Negativos (VN)
Falsos Positivo (FP)
Falsos Negativos (FN)
Y a partir de ellos:
- La Sensitividad (VP/(VP+FN)): Utilidad del modelo para determinar el “voluntario” (modelo predice a un voluntario que realmente era voluntario en la data original’).
ConfMatrix$Sensitivity## [1] 0.1474037- La Especificidad (VN/(VN+FP)): Utilidad del modelo para determinar el “no voluntario” (modelo predice a un no voluntario que realmente no era voluntario en la data original’).
ConfMatrix$Specificity## [1] 0.8932039Una medida general de resumen es el índice de Kappa: Este índice vale 1 si la predicción es perfecta, y 0 si fallan todas las predicciones.
ConfMatrix$Kappa## [1] 0.044976353.5 Prediciendo valores concretos
Si nos quedamos con el modelo III, podemos utilizarlo para predicciones. Veamos:
- Ejemplo 1: Predecir probabilidad de voluntariado de una mujer con nivel de neuroticismo=13 y extraversion=8:
ToInput1 <- data.frame(sex=factor(c("female")),
neuroticism=13, extraversion=8)
predict(object = rlog3, newdata = ToInput1, type = "response")## 1
## 0.3767916- Ejemplo 2: Predecir probabilidad de voluntariado de una mujer y hombre con los mismos niveles anteriores de neuroticismo y extraversion.
ToInput2 <- data.frame(sex=factor(c("female","male")),
neuroticism=13, extraversion=8)
predict(object = rlog3, newdata = ToInput2, type = "response")## 1 2
## 0.3767916 0.3233650- Ejemplo 3: Predecir probabilidad de voluntariado de dos mujeres, con nivel de neuroticismo de 13 y 10, y extraversion de 8 y 21, respectivamente:
ToInput3 <- data.frame(sex=factor(c("female","female")),
neuroticism=c(13,10), extraversion=c(8,21))
predict(object = rlog3, newdata = ToInput3, type = "response")## 1 2
## 0.3767916 0.58417983.6 Efectos Marginales
Recuérdese qué tenemos en volunteer, predicted y predictedProbs :
head(vol[,c('volunteer','predicted','predictedProbs')])## volunteer predicted predictedProbs
## 1 no no 0.4619546
## 2 no no 0.4080078
## 3 no no 0.4356976
## 4 no yes 0.5648592
## 5 no no 0.4914463
## 6 no no 0.4210156Cuando analizábamos le regresión Gaussiana decíamos que si el covariado aumentaba en una unidad, el coeficiente del covariado representaba el cambio promedio en la dependiente. Ese es el efecto marginal del covariado.
Esa interpretación no se puede saber con ninguno de los valores calculados hasta ahora para la regresión logística. Para ello, necesitamos calcular de manera expresa esos efectos marginales (Ford 2022):
# interpracion usando marginal effects:
library(margins)
#
margins(rlog3)## Average marginal effects## glm(formula = h3, family = binomial, data = vol)## neuroticism extraversion sexfemale
## 0.00152 0.01585 0.05623Estos coeficientes sí indican cuanto aumentaría en promedio la probabilidad de ser voluntario, cuando el predictor aumenta en una unidad. Sin embargo, recordemos que sabemos de antemano que no todos estos covariados son significativos. Veamos los intervalos de confianza en la Tabla 3.9:
marginalsData=summary(margins(rlog3))
marginalsData%>% kable(caption = "Efectos Marginales Promedio (AME)- Modelo III") %>%kableExtra::kable_styling(full_width = T)| factor | AME | SE | z | p | lower | upper |
|---|---|---|---|---|---|---|
| extraversion | 0.0158493 | 0.0033068 | 4.7928770 | 0.0000016 | 0.0093680 | 0.0223306 |
| neuroticism | 0.0015202 | 0.0027128 | 0.5603963 | 0.5752092 | -0.0037968 | 0.0068373 |
| sexfemale | 0.0562289 | 0.0264868 | 2.1229032 | 0.0337620 | 0.0043157 | 0.1081420 |
Los valores de la Tabla 3.9 podemos verlos gráficamente en la Figura 3.3.
library(ggplot2)
base= ggplot(marginalsData,aes(x=factor, y=AME)) + geom_point()
base + geom_errorbar(aes(ymin=lower, ymax=upper))
Figure 3.3: Efectos Marginales Promedio (AME)
Ya sabemos que si el intervalo incluye al cero, ese covariado no será significativo.
Bibliografía
Al igual que en la Poisson, no queremos que el intervalo de confianza incluya el uno.↩︎
Nota que si el coeficiente fuese menor que UNO, el valor que se obtiene debe interpretarse como el porcentaje de reducción.↩︎
Nota que RMSE se mantiene constante (al menos con las cifras mostradas).↩︎
Para evitar mensajes de advertencia, instale además ‘ggimage’ y ‘rsvg’, luego de intala ‘cvms’.↩︎