SOC285: Ejercicios regresión logística - Violencia Contra la Mujer
David Sulmont
1/6/2021
Introducción
Para estos ejercicios trabajaremos con la base de datos sobre violencia física y sexual intraconyugal “endes_viol.rda”, elaborada a partir de los datos de la Encuesta Nacional de Demografía y Salud Familiar – ENDES, 2004-2007, realizada por el Instituto Nacional de Estadística e Informática.
Se recomienda revisar el texto de Matos y Sulmont 2009.
Carga de paquetes y datos
library(ggplot2)
library(stargazer)
library(DescTools)
load("endes_viol.rda")Trabajaremos con una base de datos que excluya todos los casos que tienen valores perdidos en las variables que utilizaremos en este ejercicio
endes.viol.a <- na.omit(endes_viol)Incidencia de la Violencia Física
prop.table(table(endes.viol.a$vio_gen))##
## 0 1
## 0.6093315 0.3906685Modelos logísticos
Primer Modelo: Violencia física según edad de la primera relación sexual de la mujer
m1 <- glm(vio_gen ~ edadsex1, family = "binomial", data = endes.viol.a)
summary(m1)##
## Call:
## glm(formula = vio_gen ~ edadsex1, family = "binomial", data = endes.viol.a)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2039 -1.0203 -0.9065 1.3104 2.1285
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.875433 0.119739 7.311 2.65e-13 ***
## edadsex1 -0.073928 0.006641 -11.133 < 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: 11529 on 8615 degrees of freedom
## Residual deviance: 11397 on 8614 degrees of freedom
## AIC: 11401
##
## Number of Fisher Scoring iterations: 4PseudoR2(m1, "Nagel")## Nagelkerke
## 0.02065614¿Cuál sería la probabilidad estimada que una mujer que tuvo su primera relación sexual a los 18 años sea víctima de violencia física por parte de su pareja?
Usamos la fórmula para estimar esa probabilidad:
\[\pi = \frac{e^{b_{0}+b_{1}X}}{1+e^{b_{0}+b_{1}X}}\] Usando los parámetros del modelo:
\[\pi = \frac{e^{(0.875-(0.074*18))}}{1+e^{(0.875-(0.074*18))}}\] En el R, la sintaxis para esa fórmula es:
exp(0.875-(0.074*18))/(1+(exp(0.875-(0.074*18))))## [1] 0.3876977Usando el modelo podemos estimar la probabilidad de que una mujer sea víctima de violencia por parte de su pareja para todas las edades en que las mujeres tuvieron su primera relación sexual:
table(endes.viol.a$edadsex1)##
## 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## 1 8 70 278 685 1133 1216 1139 1183 654 649 368 344 208 170 143
## 26 27 28 29 30 31 32 33 34 35 36 37 38 40 41
## 90 75 55 36 42 13 16 15 10 5 2 5 1 1 1Creamos una base de datos que tiene el rango de edades de 10 a 41 años
new.data1.a <- data.frame(edadsex1 = c(10:41))Estimamos la probabilidad de nuestra variable dependiente para cada una de esas edades:
new.data1.a$pred.y <- predict(m1, new.data1.a, type = "response")
new.data1.a## edadsex1 pred.y
## 1 10 0.5339847
## 2 11 0.5155501
## 3 12 0.4970730
## 4 13 0.4786040
## 5 14 0.4601932
## 6 15 0.4418904
## 7 16 0.4237442
## 8 17 0.4058013
## 9 18 0.3881065
## 10 19 0.3707019
## 11 20 0.3536267
## 12 21 0.3369169
## 13 22 0.3206051
## 14 23 0.3047201
## 15 24 0.2892871
## 16 25 0.2743272
## 17 26 0.2598581
## 18 27 0.2458936
## 19 28 0.2324439
## 20 29 0.2195156
## 21 30 0.2071124
## 22 31 0.1952347
## 23 32 0.1838802
## 24 33 0.1730441
## 25 34 0.1627192
## 26 35 0.1528965
## 27 36 0.1435650
## 28 37 0.1347125
## 29 38 0.1263254
## 30 39 0.1183890
## 31 40 0.1108879
## 32 41 0.1038061Podemos graficar el resultado
ggplot(new.data1.a, aes(x=edadsex1, y=pred.y)) + geom_line() +
ylab("Probabilidad estimada") + xlab("Edad de 1ra relación sexual") +
ggtitle("Probabilidad de que una mujer sea víctima de violencia física \nsegún edad de su primera relación sexual") +
theme_bw()
¿Qué otras variables pueden influir?
m2 <- glm(vio_gen ~ edadsex1 + qingres, family = "binomial",
data = endes.viol.a)
m3 <- glm(vio_gen ~ edadsex1 + ledu_m, family = "binomial",
data = endes.viol.a)
m4 <- glm(vio_gen ~ edadsex1 + muj_trab, family = "binomial",
data = endes.viol.a)
stargazer(m1, m2, m3, m4, type = "text")##
## =============================================================
## Dependent variable:
## -------------------------------------------
## vio_gen
## (1) (2) (3) (4)
## -------------------------------------------------------------
## edadsex1 -0.074*** -0.075*** -0.074*** -0.074***
## (0.007) (0.007) (0.007) (0.007)
##
## qingres 0.013
## (0.018)
##
## ledu_m 0.003
## (0.015)
##
## muj_trab 0.256***
## (0.048)
##
## Constant 0.875*** 0.864*** 0.876*** 0.710***
## (0.120) (0.121) (0.120) (0.124)
##
## -------------------------------------------------------------
## Observations 8,616 8,616 8,616 8,616
## Log Likelihood -5,698.355 -5,698.119 -5,698.336 -5,684.154
## Akaike Inf. Crit. 11,400.710 11,402.240 11,402.670 11,374.310
## =============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Veamos este modelo:
m.x <- glm(vio_gen ~ dif_edu + lindig_m + edad_m + edadsex1 + muj_trab +
control, family = "binomial", data = endes.viol.a)
summary(m.x)##
## Call:
## glm(formula = vio_gen ~ dif_edu + lindig_m + edad_m + edadsex1 +
## muj_trab + control, family = "binomial", data = endes.viol.a)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3585 -0.8808 -0.6584 1.0795 2.0915
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.021948 0.153618 -6.653 2.88e-11 ***
## dif_edu -0.045131 0.020036 -2.253 0.02429 *
## lindig_m -0.118256 0.069250 -1.708 0.08770 .
## edad_m 0.029909 0.003113 9.608 < 2e-16 ***
## edadsex1 -0.075332 0.007371 -10.220 < 2e-16 ***
## muj_trab 0.171981 0.053418 3.220 0.00128 **
## control 0.539564 0.017233 31.310 < 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: 11529 on 8615 degrees of freedom
## Residual deviance: 10068 on 8609 degrees of freedom
## AIC: 10082
##
## Number of Fisher Scoring iterations: 4PseudoR2(m.x, "Nagel")## Nagelkerke
## 0.2114478Vamos a estimar las probabilidades de ser víctima de violencia, según los diferentes niveles de la variable de “control machista”, manteniendo constantes las otras variables. Para ello usaremos la media o la moda del resto de variables independientes:
summary(endes.viol.a)## vio_gen vio_sex fregion
## Min. :0.0000 Min. :0.00000 Norte :1660
## 1st Qu.:0.0000 1st Qu.:0.00000 Sierra Sur :1138
## Median :0.0000 Median :0.00000 Sur :1189
## Mean :0.3907 Mean :0.08798 Sierra Centro : 902
## 3rd Qu.:1.0000 3rd Qu.:0.00000 Sierra centro-sur: 762
## Max. :1.0000 Max. :1.00000 Selva :1934
## Costa central :1031
## fresid1 ledu_m lindig_m edad_m
## Capital, ciudad: 649 Min. :0.000 Min. :0.0000 Min. :15.00
## Pequeña ciudad :2723 1st Qu.:1.000 1st Qu.:0.0000 1st Qu.:27.00
## Pueblo :1220 Median :3.000 Median :0.0000 Median :33.00
## Rural :4024 Mean :2.807 Mean :0.1508 Mean :33.42
## 3rd Qu.:4.000 3rd Qu.:0.0000 3rd Qu.:40.00
## Max. :5.000 Max. :1.0000 Max. :49.00
##
## edad_1c edadsex1 qingres ledu_c
## Min. :10.00 Min. :10.00 Min. :1.000 Min. :0.000
## 1st Qu.:17.00 1st Qu.:15.00 1st Qu.:2.000 1st Qu.:2.000
## Median :19.00 Median :17.00 Median :3.000 Median :4.000
## Mean :19.93 Mean :17.98 Mean :2.872 Mean :3.191
## 3rd Qu.:22.00 3rd Qu.:20.00 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :45.00 Max. :41.00 Max. :5.000 Max. :5.000
##
## muj_trab dif_age dif_edu m_gana
## Min. :0.0000 Min. :-25.00 Min. :-4.0000 Min. :0.00000
## 1st Qu.:0.0000 1st Qu.: 1.00 1st Qu.: 0.0000 1st Qu.:0.00000
## Median :1.0000 Median : 3.00 Median : 0.0000 Median :0.00000
## Mean :0.6798 Mean : 3.96 Mean : 0.3843 Mean :0.06801
## 3rd Qu.:1.0000 3rd Qu.: 7.00 3rd Qu.: 1.0000 3rd Qu.:0.00000
## Max. :1.0000 Max. : 41.00 Max. : 5.0000 Max. :1.00000
##
## trago1 trago2 control just_vio
## Min. :0.0000 Min. :0.0000 Min. :0.000 Min. :0.00000
## 1st Qu.:1.0000 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.:0.00000
## Median :1.0000 Median :1.0000 Median :1.000 Median :0.00000
## Mean :0.7637 Mean :0.6815 Mean :1.506 Mean :0.06012
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:2.000 3rd Qu.:0.00000
## Max. :1.0000 Max. :1.0000 Max. :6.000 Max. :1.00000
##
## vio_fam vio_mam year_en
## Min. :0.00000 Min. :0.0000 Min. :2005
## 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:2005
## Median :0.00000 Median :0.0000 Median :2006
## Mean :0.05989 Mean :0.4836 Mean :2006
## 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:2007
## Max. :1.00000 Max. :1.0000 Max. :2007
## new.data.x1 <- data.frame(dif_edu = .38, lindig_m= 0, edad_m = 30.42,
edadsex1 = 17.98, muj_trab = 1,
control = c(0:6))
new.data.x1$pred.y <- predict(m.x, new.data.x1, type = "response")
new.data.x1## dif_edu lindig_m edad_m edadsex1 muj_trab control pred.y
## 1 0.38 0 30.42 17.98 1 0 0.2121947
## 2 0.38 0 30.42 17.98 1 1 0.3160070
## 3 0.38 0 30.42 17.98 1 2 0.4421059
## 4 0.38 0 30.42 17.98 1 3 0.5761391
## 5 0.38 0 30.42 17.98 1 4 0.6998340
## 6 0.38 0 30.42 17.98 1 5 0.7999643
## 7 0.38 0 30.42 17.98 1 6 0.8727655Y graficamos nuestras predicciones:
ggplot(new.data.x1, aes(x=control, y = pred.y)) + geom_line() +
ylab("Probabilidad estimada") + xlab("Escala de control machista") +
ggtitle("Probabilidad de que una mujer sea víctima de violencia física \nsegún escala de control machista") +
theme_bw()