library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(haven)
library(haven)
Casen_2017 <- read_dta("C:/Users/anavc/OneDrive/Escritorio/Casen 2017.dta")
View(Casen_2017)
dataset <- Casen_2017%>%
select(e6a, edad, ch1)%>%
sample_n(size = 5000)%>%
na.omit()
names(dataset) <- c("Educacion", "Edad", "Ocupacion")
head(dataset)
## # A tibble: 6 × 3
## Educacion Edad Ocupacion
## <dbl+lbl> <dbl> <dbl+lbl>
## 1 7 [Educación Básica] 7 4 [Menores …
## 2 7 [Educación Básica] 11 4 [Menores …
## 3 4 [Prekinder/Kinder (Transición menor y Transición Mayor)] 6 4 [Menores …
## 4 11 [Educación Media Técnica Profesional] 44 2 [Patrón o…
## 5 9 [Educación Media Científico-Humanista] 15 3 [Familiar…
## 6 2 [Sala cuna] 1 4 [Menores …
## 1- COMPROBAR LOS SUPUESTO DE NORMALIDAD Y HOMOGENEIDAD
# Realizar la prueba de normalidad para cada variable
resultados_shapiro <- sapply(dataset, function(x) shapiro.test(x)$p.value)
# Verificar los resultados
resultados_shapiro
## Educacion Edad Ocupacion
## 1.074636e-80 1.019581e-30 3.203250e-59
# Realizar la prueba de homogeneidad de varianzas para todas las variables
resultados_bartlett <- bartlett.test(dataset[, c("Educacion", "Edad", "Ocupacion")], dataset$Ocupacion)
# Verificar los resultados
resultados_bartlett
##
## Bartlett test of homogeneity of variances
##
## data: dataset[, c("Educacion", "Edad", "Ocupacion")]
## Bartlett's K-squared = 28143, df = 2, p-value < 2.2e-16
#### No se cumple la normalidad, por que no se cumple una distribucion normal y tampoco se cumple la homogeneidad, ya que la varianza entre las variables no son iguales. Es por esto que vamos a utilizar los modelos de logit y probit.
# Creación del rango etario
dataset$edad_cat <- cut(dataset$Edad, breaks = 3, labels = c("joven", "adulto", "mayor"))
# Convertir variables categóricas a factor
dataset$Educacion <- factor(dataset$Educacion)
dataset$edad_cat <- factor(dataset$edad_cat)
# Crear una nueva columna en el dataframe llamada "posgrado", la cual será dicotómica
dataset$Ocupacion <- ifelse(dataset$Ocupacion == 1 , 1, 0)
dataset$Ocupacion <- factor(dataset$Ocupacion)
summary(dataset)
## Educacion Edad Ocupacion edad_cat
## 9 :1339 Min. : 0.00 0:3450 joven :2242
## 7 :1248 1st Qu.:19.00 1:1550 adulto:2055
## 15 : 388 Median :37.00 mayor : 703
## 11 : 377 Mean :37.91
## 14 : 282 3rd Qu.:55.00
## 6 : 272 Max. :96.00
## (Other):1094
# Ajustar modelo logit
modelo_logit <- glm(Ocupacion ~ edad_cat + Educacion , family = binomial(link = "logit"), data = dataset)
# Ajustar modelo probit
modelo_probit <- glm(Ocupacion ~ edad_cat + Educacion , family = binomial(link = "probit"), data = dataset)
# Resumen de los modelos
summary(modelo_logit)
##
## Call:
## glm(formula = Ocupacion ~ edad_cat + Educacion, family = binomial(link = "logit"),
## data = dataset)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0006 -0.8362 -0.5585 1.1006 2.9275
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.58196 0.38572 -9.287 < 2e-16 ***
## edad_catadulto 0.96088 0.07324 13.120 < 2e-16 ***
## edad_catmayor -0.68923 0.16758 -4.113 3.91e-05 ***
## Educacion2 0.58623 1.09489 0.535 0.592355
## Educacion3 -12.98410 288.87150 -0.045 0.964149
## Educacion4 -12.99467 217.84757 -0.060 0.952434
## Educacion5 0.10248 1.10124 0.093 0.925857
## Educacion6 1.26319 0.46822 2.698 0.006978 **
## Educacion7 1.80276 0.39211 4.598 4.27e-06 ***
## Educacion8 1.83405 0.48263 3.800 0.000145 ***
## Educacion9 2.60677 0.39008 6.683 2.35e-11 ***
## Educacion10 2.11432 0.59212 3.571 0.000356 ***
## Educacion11 2.80435 0.40019 7.008 2.42e-12 ***
## Educacion12 3.09464 0.42789 7.232 4.75e-13 ***
## Educacion13 3.55110 0.40939 8.674 < 2e-16 ***
## Educacion14 2.61266 0.40651 6.427 1.30e-10 ***
## Educacion15 3.40027 0.40008 8.499 < 2e-16 ***
## Educacion16 4.47711 0.77837 5.752 8.83e-09 ***
## Educacion17 3.61869 0.52258 6.925 4.37e-12 ***
## Educacion99 2.83689 0.67936 4.176 2.97e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 6191.0 on 4999 degrees of freedom
## Residual deviance: 5123.3 on 4980 degrees of freedom
## AIC: 5163.3
##
## Number of Fisher Scoring iterations: 15
summary(modelo_probit)
##
## Call:
## glm(formula = Ocupacion ~ edad_cat + Educacion, family = binomial(link = "probit"),
## data = dataset)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0119 -0.8558 -0.5363 1.0987 2.9982
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.96212 0.16813 -11.670 < 2e-16 ***
## edad_catadulto 0.60260 0.04399 13.698 < 2e-16 ***
## edad_catmayor -0.32249 0.09241 -3.490 0.000483 ***
## Educacion2 0.29373 0.49776 0.590 0.555126
## Educacion3 -3.43425 72.84592 -0.047 0.962398
## Educacion4 -3.46080 54.90775 -0.063 0.949743
## Educacion5 -0.05354 0.52277 -0.102 0.918427
## Educacion6 0.59722 0.21499 2.778 0.005471 **
## Educacion7 0.85418 0.17287 4.941 7.77e-07 ***
## Educacion8 0.87706 0.22824 3.843 0.000122 ***
## Educacion9 1.35475 0.17183 7.884 3.17e-15 ***
## Educacion10 1.06182 0.30331 3.501 0.000464 ***
## Educacion11 1.47719 0.18052 8.183 2.77e-16 ***
## Educacion12 1.65631 0.20353 8.138 4.03e-16 ***
## Educacion13 1.91941 0.18784 10.218 < 2e-16 ***
## Educacion14 1.35970 0.18513 7.344 2.07e-13 ***
## Educacion15 1.82955 0.18034 10.145 < 2e-16 ***
## Educacion16 2.47585 0.42673 5.802 6.56e-09 ***
## Educacion17 1.96043 0.27054 7.246 4.28e-13 ***
## Educacion99 1.48950 0.38272 3.892 9.95e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 6191.0 on 4999 degrees of freedom
## Residual deviance: 5116.2 on 4980 degrees of freedom
## AIC: 5156.2
##
## Number of Fisher Scoring iterations: 15
# Comprobación de multicolinealidad (criterio vif < 10)
library(faraway)
vif(modelo_logit)
## edad_catadulto edad_catmayor Educacion2 Educacion3 Educacion4
## 6.492136e+00 1.696758e+01 2.506860e+01 5.678367e+06 5.603400e+06
## Educacion5 Educacion6 Educacion7 Educacion8 Educacion9
## 2.415779e+01 5.638608e+01 1.439883e+02 3.323369e+01 1.491799e+02
## Educacion10 Educacion11 Educacion12 Educacion13 Educacion14
## 1.253111e+01 5.582405e+01 2.179232e+01 3.829317e+01 4.397182e+01
## Educacion15 Educacion16 Educacion17 Educacion99
## 5.728591e+01 7.855729e+00 1.056736e+01 6.443295e+00
vif(modelo_probit)
## edad_catadulto edad_catmayor Educacion2 Educacion3 Educacion4
## 2.342507e+00 5.158827e+00 5.181185e+00 3.610976e+05 3.559701e+05
## Educacion5 Educacion6 Educacion7 Educacion8 Educacion9
## 5.443829e+00 1.188826e+01 2.798725e+01 7.432372e+00 2.894794e+01
## Educacion10 Educacion11 Educacion12 Educacion13 Educacion14
## 3.287993e+00 1.135939e+01 4.930638e+00 8.061409e+00 9.120340e+00
## Educacion15 Educacion16 Educacion17 Educacion99
## 1.163910e+01 2.361122e+00 2.832234e+00 2.044869e+00
# Comprobación de la normalidad de los residuos
residuos_logit <- residuals(modelo_logit, type = "pearson")
residuos_probit <- residuals(modelo_probit, type = "pearson")
shapiro.test(residuos_logit)
##
## Shapiro-Wilk normality test
##
## data: residuos_logit
## W = 0.88031, p-value < 2.2e-16
shapiro.test(residuos_probit)
##
## Shapiro-Wilk normality test
##
## data: residuos_probit
## W = 0.87574, p-value < 2.2e-16
# Comprobación de la homocedasticidad
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
bptest(modelo_logit)
##
## studentized Breusch-Pagan test
##
## data: modelo_logit
## BP = 1122.4, df = 19, p-value < 2.2e-16
bptest(modelo_probit)
##
## studentized Breusch-Pagan test
##
## data: modelo_probit
## BP = 1122.4, df = 19, p-value < 2.2e-16
# Cargar paquete pROC
library(pROC)
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
# Evaluación de la bondad de ajuste utilizando la curva ROC
roc_logit <- roc(dataset$Ocupacion, fitted(modelo_logit))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc_probit <- roc(dataset$Ocupacion, fitted(modelo_probit))
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
# Mostrar el área bajo la curva (AUC) para cada modelo
cat("AUC del modelo logit:", auc(roc_logit), "\n")
## AUC del modelo logit: 0.7674855
cat("AUC del modelo probit:", auc(roc_probit), "\n")
## AUC del modelo probit: 0.7657415
library(ggplot2)
library(pROC)
# Calcular las curvas ROC de los modelos
roc_logit <- roc(dataset$Ocupacion, modelo_logit$fitted.values)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
roc_probit <- roc(dataset$Ocupacion, modelo_probit$fitted.values)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
# Crear un dataframe con los valores de las curvas ROC
df_roc_logit <- data.frame(fpr = roc_logit$specificities, tpr = roc_logit$sensitivities, model = "Logit")
df_roc_probit <- data.frame(fpr = roc_probit$specificities, tpr = roc_probit$sensitivities, model = "Probit")
df_roc <- rbind(df_roc_logit, df_roc_probit)
# Graficar las curvas ROC utilizando ggplot2
ggplot(df_roc, aes(x = fpr, y = tpr, color = model)) +
geom_line() +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
ggtitle("Curvas ROC de los modelos logit y probit") +
xlab("Tasa de falsos positivos") +
ylab("Tasa de verdaderos positivos")
# Comparación de modelos
logLik(modelo_logit)
## 'log Lik.' -2561.66 (df=20)
logLik(modelo_probit)
## 'log Lik.' -2558.12 (df=20)
# Comparar los modelos utilizando anova()
resultado <- anova(modelo_logit, modelo_probit, test = "Chisq")
# Imprimir el resultado
print(resultado)
## Analysis of Deviance Table
##
## Model 1: Ocupacion ~ edad_cat + Educacion
## Model 2: Ocupacion ~ edad_cat + Educacion
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 4980 5123.3
## 2 4980 5116.2 0 7.0794
##ANALISIS: Con los resultados de los obtenido podemos concluir que el ser asalariado en Chile está relacionado con el nivel de educación y el rango etario. Esto debido a que los modelos logit y probit muestran que la educación y el rango etario tienen un efecto significativo en la probabilidad de ser asalariado.