library(haven)
data <- read_dta("Data1_R.dta")
View(data)DEBER Tarea AA 3 U1
Cambiar directorio de trabajo
setwd(“C:/Users/pajv2/OneDrive/Escritorio/DEBER Tarea AA 3 U1/DEBER Tarea AA 3 U1_files”)
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)
library(margins)
Setear directorio
setwd(“C:/Users/pajv2/OneDrive/Escritorio/DEBER Tarea AA 3 U1”)
Leer los datos (ajustar la ruta del archivo)
print(names(data)) [1] "area" "empleo" "region"
[4] "edad" "t_hijos" "nac_vivo_murieron"
[7] "mortinato_2" "depresion_pp" "intensidad_dpp"
[10] "etnia" "f2_s2_216_1" "f2_s2_216_2"
[13] "f2_s2_218_1_a" "tiempo_dpp" "f2_s5_504a_1"
[16] "f2_s5_504b_1" "f2_s5_504c_1" "f2_s5_504d_1"
[19] "f2_s5_504e_1" "f2_s5_504f_1" "f2_s5_504g_1"
[22] "f2_s5_504h_1" "f2_s5_504i_1" "f2_s5_504j_1"
[25] "f2_s5_504k_1" "est_civil" "q_usted"
[28] "q_pareja" "f2_s4b_406_" "lugar"
[31] "form_parto" "c_posparto" "años_esc"
[34] "ingrl" "planf" "lingrl"
[37] "lugarx1" "lugarx2" "regionx1"
[40] "regionx2" "regionx3" "regionx4"
[43] "intensidad_dppx1" "intensidad_dppx2" "intensidad_dppx3"
[46] "tiempo_dppx1" "tiempo_dppx2" "tiempo_dppx3"
[49] "tiempo_dppx4" "tiempo_dppx5" names(data)[names(data) == "años_esc"] <- "anios_esc"print(names(data)) [1] "area" "empleo" "region"
[4] "edad" "t_hijos" "nac_vivo_murieron"
[7] "mortinato_2" "depresion_pp" "intensidad_dpp"
[10] "etnia" "f2_s2_216_1" "f2_s2_216_2"
[13] "f2_s2_218_1_a" "tiempo_dpp" "f2_s5_504a_1"
[16] "f2_s5_504b_1" "f2_s5_504c_1" "f2_s5_504d_1"
[19] "f2_s5_504e_1" "f2_s5_504f_1" "f2_s5_504g_1"
[22] "f2_s5_504h_1" "f2_s5_504i_1" "f2_s5_504j_1"
[25] "f2_s5_504k_1" "est_civil" "q_usted"
[28] "q_pareja" "f2_s4b_406_" "lugar"
[31] "form_parto" "c_posparto" "anios_esc"
[34] "ingrl" "planf" "lingrl"
[37] "lugarx1" "lugarx2" "regionx1"
[40] "regionx2" "regionx3" "regionx4"
[43] "intensidad_dppx1" "intensidad_dppx2" "intensidad_dppx3"
[46] "tiempo_dppx1" "tiempo_dppx2" "tiempo_dppx3"
[49] "tiempo_dppx4" "tiempo_dppx5" 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: 4“Análisis”
El ingreso y los años de escolaridad no son estadisticamente significativos, es decir no ayudan a explicar la probabilidad de que las mujeres ecuatorianas sufran de depresión.
La edad es estadisticamente significativa lo que se traduce en que, las mujeres con mayor edad tienen probabilidad de sufrir de depresión posparto.
El número de hijos presenta significancia estadística, es decir, que a medidad que aumenta el número de hijos incrementa la probabilidad de que las mujeres sufran de depresión posparto.
La etnia, es una variable estadísticamente significativa en el modelo lo que significa que, las mujeres indígenas tienen mayor probabilidad de padecer de depresión posparto.
Las mujeres del área rural tienen mayor probabilidad de sufrir de depresión posparto en comparación a las mujeresd el área 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 el 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 cat("BIC Logit:", bic_logit, " | BIC Probit:", bic_probit, "\n")BIC Logit: 17173.34 | BIC Probit: 17171.09 Conclusión: 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.0129###Calcular efectos marginales para Probit
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.0129“Análisis”
Los años de escolaridad y el ingreso no son estadisticamente significativos, es decir no ayudan a explicar la probabilidad de que las mujeres ecuatorianas padescan de depresión posparto.
El area al ser una variable estadisticamente significativa nos explica que las mujeres del area ruraltiene 1.84 más probabilidad de sufrir de depresión posparto.
Edad, un año adicional en las mujeres ecuatorianas en promedio aumenta la probabilidad en 0.56 % de sufrir de depresión posparto
Una mujer indígena tiene 5.92% de probabilidad de sufrir de depresión posparto en 0.66 de que padecan de depresión posparto.
###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.0132“Análisis”
Los años de escolaridad y el ingreso no son estadisticamente significativos, es decir no ayudan a explicar la probabilidad de que las mujeres ecuatorianas padescan de depresión posparto.
El area al ser una variable estadisticamente significativa nos explica que las mujeres del area ruraltiene 1.88más probabilidad de sufrir de depresión posparto.
Edad, un año adicional en las mujeres ecuatorianas en promedio aumenta la probabilidad en 0.57% de sufrir de depresión posparto
Una mujer indígena tiene 0.06%de probabilidad de sufrir de depresión posparto en 0.66 de que padecan de depresión posparto.
###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 probabilidada de que las mujeres con depresión posparto ocurra de acuerdo al modelo logic es del 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 < cases###Graficar 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 ##urva ROC: Muestra el rendimiento del modelo en diferentes umbrales de clasificación. AUC > 0.7 indica un buen modelo predictivo.
(Aumentar más variables de control)
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)clinica <- ifelse(data$f2_s5_504c_2 == 2, 1, 0)Warning: Unknown or uninitialised column: `f2_s5_504c_2`.##Modelo Logit con mas variables de control.
data$casado <- ifelse(as.numeric(data$est_civil) == 1, 1, 0)modelo_logit <- glm(depresion_pp ~ lingrl + anios_esc + edad + t_hijos + etnia + area + medico + casado,
data = data, family = binomial("logit"))###Modelo Logic 2
library(pROC)
roc_logit <- roc(data$depresion_pp, predict(modelo_logit, type = "response"))Setting levels: control = 0, case = 1Setting direction: controls < casesexactitud <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat("Exactitud del modelo logit:", exactitud, "\n")Exactitud del modelo logit: 0.7797702 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.5859612 ###Modelo Logic 3
medico <- ifelse(data$f2_s5_504c_1 == 3, 1, 0)data$estado_civil <- ifelse(data$est_civil == 1, 0, 1)modelo_logit <- glm(depresion_pp ~ lingrl + anios_esc + edad + t_hijos + etnia + area + medico + casado,
data = data, family = binomial(link = "logit"))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()
auc_logit <- auc(roc_logit)
cat("Área bajo la curva (AUC) - Modelo Logit:", auc_logit, "\n")Área bajo la curva (AUC) - Modelo Logit: 0.5859612