library(haven)
data <- read_dta("Data1_R.dta")
View(data)clase 25 mayo
Instalar paquetes si es necesario
install.packages(“MASS”) # Para modelo probit install.packages(“car”) # Para pruebas estadísticas install.packages(“margins”) # Para calcular efectos marginales install.packages(“pROC”) # Para curva ROC install.packages(“ggplot2”) # Para visualización install.packages(“caret”) # Para matriz de confusión install.packages(“ggplot”) # Para visualización
Cargar librerías
library(MASS) library(car) library(margins) library(pROC) library(ggplot2) library(ggplot) library(caret) library(haven)
Leer datos
Ver las primeras filas de la base de datos
head(data)# A tibble: 6 × 50
area empleo region edad t_hijos nac_vivo_murieron mortinato_2
<dbl+lbl> <dbl+lbl> <dbl+l> <dbl> <dbl> <dbl+lbl> <dbl+lbl>
1 1 [Urbano] 1 [Trabajó al … 1 [Sie… 19 1 0 [No] 0 [No]
2 1 [Urbano] 0 [No trabajó] 1 [Sie… 23 1 0 [No] 0 [No]
3 1 [Urbano] 1 [Trabajó al … 1 [Sie… 38 5 0 [No] 0 [No]
4 1 [Urbano] 0 [No trabajó] 1 [Sie… 18 1 0 [No] 0 [No]
5 1 [Urbano] 0 [No trabajó] 1 [Sie… 21 1 0 [No] 0 [No]
6 1 [Urbano] 1 [Trabajó al … 1 [Sie… 22 1 0 [No] 0 [No]
# ℹ 43 more variables: depresion_pp <dbl+lbl>, intensidad_dpp <dbl+lbl>,
# etnia <dbl+lbl>, f2_s2_216_1 <dbl+lbl>, f2_s2_216_2 <dbl>,
# f2_s2_218_1_a <dbl+lbl>, tiempo_dpp <dbl+lbl>, f2_s5_504a_1 <dbl+lbl>,
# f2_s5_504b_1 <dbl+lbl>, f2_s5_504c_1 <dbl+lbl>, f2_s5_504d_1 <dbl+lbl>,
# f2_s5_504e_1 <dbl+lbl>, f2_s5_504f_1 <dbl+lbl>, f2_s5_504g_1 <dbl+lbl>,
# f2_s5_504h_1 <dbl+lbl>, f2_s5_504i_1 <dbl+lbl>, f2_s5_504j_1 <dbl+lbl>,
# f2_s5_504k_1 <dbl+lbl>, est_civil <dbl+lbl>, q_usted <dbl+lbl>, …Revisar estructura de los datos
str(data)Ejemplo 1: Modelos con variable dependiente dicotómica
Modelos logit y probit
Ajustar el modelo logit
modelo_logit <- glm(depresion_pp ~ lingrl + anios_esc + edad + t_hijos + etnia + area,
data = data, family = binomial(link = "logit"))Ajustar el modelo probit
modelo_probit <- glm(depresion_pp ~ lingrl + anios_esc + edad + t_hijos + etnia + area,
data = data, family = binomial(link = "probit"))Resumen modelo Logit
summary(modelo_logit)
Call:
glm(formula = depresion_pp ~ lingrl + anios_esc + edad + t_hijos +
etnia + area, family = binomial(link = "logit"), data = data)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.3377859 0.1015521 -23.021 < 2e-16 ***
lingrl 0.0006157 0.0071763 0.086 0.9316
anios_esc -0.0078052 0.0049109 -1.589 0.1120
edad 0.0333503 0.0032243 10.344 < 2e-16 ***
t_hijos 0.0391392 0.0189765 2.063 0.0392 *
etnia 0.3502255 0.0605997 5.779 7.5e-09 ***
area 0.1089295 0.0425378 2.561 0.0104 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 17346 on 16450 degrees of freedom
Residual deviance: 17105 on 16444 degrees of freedom
AIC: 17119
Number of Fisher Scoring iterations: 4Analisis:
- El ingreso y los años de escolaridad no son no son estadisticamente significativos, es decir no ayudan a explicar la probabilidad de que las mujeres ecuatorianas sufran de depresion postparto
- Edad si es estadisticamnete sifnificativa, lo que representa que las mujeres con mayor edad tienen probabilidad de sufrir de depresion postparto.
- El numero de hijos presenta significancia estadistica ya que ha medida que aumenta el numero de hijos incrementa la probabilidad de que las mujeres sufran depresion postparto
- La etnia es una variable estadisticamente significativa para el modelo, lo cual permite determinar que las mujeres indigenas tienen mayor probabilidad de padecer depresion postparto.
- Las mujeres del area rural tienen mayor probabilidad de sufrir depresion postparto en comparacion a las mujeres del area urbana.
Resumen modelo probit
summary(modelo_probit)
Call:
glm(formula = depresion_pp ~ lingrl + anios_esc + edad + t_hijos +
etnia + area, family = binomial(link = "probit"), data = data)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.401e+00 5.852e-02 -23.942 < 2e-16 ***
lingrl 3.942e-05 4.170e-03 0.009 0.99246
anios_esc -4.481e-03 2.861e-03 -1.566 0.11733
edad 1.958e-02 1.890e-03 10.363 < 2e-16 ***
t_hijos 2.334e-02 1.123e-02 2.078 0.03774 *
etnia 2.078e-01 3.585e-02 5.796 6.8e-09 ***
area 6.431e-02 2.452e-02 2.623 0.00872 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 17346 on 16450 degrees of freedom
Residual deviance: 17103 on 16444 degrees of freedom
AIC: 17117
Number of Fisher Scoring iterations: 4Comparar AIC y BIC de ambos modelos
Comparar criterio de información de Akaike (AIC) de los modelos
aic_logit <- AIC(modelo_logit)
aic_probit <- AIC(modelo_probit)###Comparar criterio de Información Bayesiano (BIC) de los modelos
bic_logit <- BIC(modelo_logit)
bic_probit <- BIC(modelo_probit)Mostrar resultados
cat("AIC Logit:", aic_logit, " | AIC Probit:", aic_probit, "\n")AIC Logit: 17119.38 | AIC Probit: 17117.13 El mejor es el Probit ya que es mas bajo
cat("BIC Logit:", bic_logit, " | BIC Probit:", bic_probit, "\n")BIC Logit: 17173.34 | BIC Probit: 17171.09 Conclusion El modelo con menor AIC/BIC es el preferido, ya que tiene mejor ajuste, es decir el modelo probit
Calcular efectos marginales para Logit
library(margins)
marg_logit <- margins(modelo_logit)
summary(marg_logit) factor AME SE z p lower upper
anios_esc -0.0013 0.0008 -1.5897 0.1119 -0.0029 0.0003
area 0.0184 0.0072 2.5619 0.0104 0.0043 0.0325
edad 0.0056 0.0005 10.4239 0.0000 0.0046 0.0067
etnia 0.0592 0.0102 5.7944 0.0000 0.0392 0.0793
lingrl 0.0001 0.0012 0.0858 0.9316 -0.0023 0.0025
t_hijos 0.0066 0.0032 2.0632 0.0391 0.0003 0.0129Analisis
- Los años de escolaridad y el ingreso no son estadisticamente significativos, es decir no ayudan a explicar la probabilida de que las mujeres Ecuatorianas padezcan de depresion postparto.
- El area al ser una variable estadisticamente significativa nos explica que las mujeres del area rural tienen 1.84% mas probabilidad de sufrir de depresion postparto.
- Edad un año adicional en las mujeres Ecuatorianas en promedio aumenta la probabilidad en 0.56% de sufrir de depresion postparto.
- La etnia una mujer indigenas tiene 5.92% de probabilidad de sufrir de depresion postparto - A medida que aumenta el numero de hijos aumenta la probabilidad en 0.66% de que padezcan de depresion postparto,
Análisis (multiplicar *100)
Calcular efectos marginales para Probit
marg_probit <- margins(modelo_probit)
summary(marg_probit) factor AME SE z p lower upper
anios_esc -0.0013 0.0008 -1.5664 0.1173 -0.0029 0.0003
area 0.0188 0.0072 2.6238 0.0087 0.0047 0.0328
edad 0.0057 0.0005 10.4392 0.0000 0.0046 0.0068
etnia 0.0606 0.0104 5.8096 0.0000 0.0402 0.0811
lingrl 0.0000 0.0012 0.0095 0.9925 -0.0024 0.0024
t_hijos 0.0068 0.0033 2.0783 0.0377 0.0004 0.0132Analisis
- Los años de escolaridad y el ingreso no son estadisticamente significativos, es decir no ayudan a explicar la probabilida de que las mujeres Ecuatorianas padezcan de depresion postparto.
- El area al ser una variable estadisticamente significativa nos explica que las mujeres del area rural tienen 1.88% mas probabilidad de sufrir de depresion postparto.
- Edad un año adicional en las mujeres Ecuatorianas en promedio aumenta la probabilidad en 0.57% de sufrir de depresion postparto. - Etnia, una mujer indigena tiene 6.06% de probabilidad de sufrir de depresion postparto
- A medida que aumenta el numero de hijos aumenta la probabilidad en 0.68% de que padezcan de depresion postparto.
Para el modelo logit podemos calcular la matriz de confusión
- Obtener predicciones del modelo logit
pred_logit <- ifelse(predict(modelo_logit, type = "response") > 0.5, 1, 0)Crear matriz de confusión
conf_matrix <- table(Predicho = pred_logit, Real = data$depresion_pp)Mostrar matriz de confusión
print(conf_matrix) Real
Predicho 0 1
0 12828 3623Calcular exactitud
exactitud <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat("Exactitud del modelo logit:", exactitud, "\n")Exactitud del modelo logit: 0.7797702 Las variables consideradas para explicar la probabilidad de que las mujeres con depresion postparto ocurra con El modelo logit es de 77.98%
Calcular la curva ROC para el modelo logit (ayuda a ver con cuanta exactitud el modelo está prediciendo los resultados)
library(pROC)Type 'citation("pROC")' for a citation.
Adjuntando el paquete: 'pROC'The following objects are masked from 'package:stats':
cov, smooth, varroc_logit <- roc(data$depresion_pp, predict(modelo_logit, type = "response"))Setting levels: control = 0, case = 1Setting direction: controls < casesGraficar la curva ROC
library(ggplot2)
ggplot() +
geom_line(aes(x = roc_logit$specificities, y = roc_logit$sensitivities), color = "blue") +
geom_abline(linetype = "dashed", color = "red") +
labs(title = "Curva ROC - Modelo Logit",
x = "1 - Especificidad",
y = "Sensibilidad") +
theme_minimal()
Mostrar el área bajo la curva (AUC)
auc_logit <- auc(roc_logit)
cat("Área bajo la curva (AUC) - Modelo Logit:", auc_logit, "\n")Área bajo la curva (AUC) - Modelo Logit: 0.5836005 El modelo esta especificado en un 58% (faltan de aumentar variables de control)
Curva ROC: Muestra el rendimiento del modelo en diferentes umbrales de clasificación.
medico <- ifelse(data$f2_s5_504c_1 == 3, 1, 0)table(data$f2_s5_504c_1)
1 2 3
405 2008 14038 data$medico <- ifelse(data$f2_s5_504c_1 == 3, 1, 0)MODELO LOGIT CON MAS VARIABLES DE CONTROL
modelo_logit <- glm(depresion_pp ~ lingrl + anios_esc + edad + t_hijos + etnia + area,
data = data, family = binomial(link = "logit"))