PEEA UNI Mod 1 Tarea 1
Felix A. Garcia Aguilar
29/4/2022
Paquete Foreing
Foreing es un paquete que permite la carga de base de datos de diferentes fuentes, como por ejemplo importar base de datos de STATA o SPSS.
Primero llamamos al paquete Foreing para poder operaralo:
library(foreign)
library(tidyverse)Luego, cambiamos el directorio de trabajo donde estaran los archivos trabajados y las bases de datos:
setwd("~/PEEA")Luego el paquete foreing, permite subir base de datos en formatos Stata y SPSS. Para este caso se importara data desde Stata (Base de Enaho):
df <- read.dta("sumaria-2020.dta",
convert.dates = TRUE, convert.factors = TRUE,
missing.type = FALSE, convert.underscore = FALSE,
warn.missing.labels = TRUE)Lo mismo aplicaria para la carga de base de datos desde un archivo spss, cambiando unicamente la palabra read.dta por read.spss. Los argumentos de convertir a fechas, lo hace siempre y cuando tengan el formato de fechas “yyyy-mm-dd”, los factores se convierten automaticamente al distinguir entre character de multiples opciones(Si se desea convertir al factor en otra clase, se coloca “numerico”, “charcater” u otro). Missing types(Si quisieramos identificar los NA como valor de 999, colocariamos 999), y connvert underscore son ajustes por cada columna de la base de datos, colocando NA si es missing value y eliminando los underscore.
Paquete Stats
stats es una funcion base que contine las estadiscas funcionales basicas del algebra matricial, lineal y calculo de regresiones, predicciones, etc. Dentro de las funciones estadisticas credas, se explicara las funciones para el modelamiento de modelos Logit y Probit.
De la funcion Logit
En este caso, se usa la base de datos de Enaho, como ejemplo determinar la probabilidad de que una persona sea pobre o no pobre. Primero seleccionamos las variables en un nuevo data frame y cargamos la libreria Tidiverse:
variable dependiente: * Pobreza (1 = Pobre, 0 = No Pobre) Variables independiente: * Miembros por hogar * Perceptores en el hogar * Ingreso neto * Ingreso percapita por hogar * Gasto monetario total
library(tidyverse)
df_pobreza <- df %>% select(pobreza, mieperho, percepho, inghog2d, gashog1d) %>% mutate(pobreza = ifelse(pobreza=="no pobre",0,1))
head(df_pobreza)## pobreza mieperho percepho inghog2d gashog1d
## 1 0 2 2 37274.43 9875.892
## 2 0 2 2 41967.57 10971.022
## 3 0 6 4 42292.54 30125.684
## 4 0 4 3 37866.21 17645.979
## 5 0 4 3 91681.96 46508.578
## 6 0 4 3 29749.26 58692.770Luego de ello realizamos un analisis descriptivo de las variables:
df_pobreza %>% mutate(pobreza = as.factor(pobreza)) %>% ggplot(aes(mieperho, color = pobreza)) + geom_boxplot()
df_pobreza %>% mutate(pobreza = as.factor(pobreza)) %>% ggplot(aes(inghog2d, color = pobreza, fill= pobreza)) + geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
df_pobreza %>% mutate(pobreza = as.factor(pobreza)) %>% ggplot(aes(inghog2d,gashog1d, color = pobreza, fill= pobreza)) + geom_point() + facet_grid(.~pobreza)
inmediatamente, pasamos a crear el modelo lineal como base de comparacion, teniendo en cuenta que “~.” implica que la variable pobreza esta en funcion de todas las variables de la base de datos actual:
fit_lm <- lm(pobreza ~. , data = df_pobreza)
summary(fit_lm)##
## Call:
## lm(formula = pobreza ~ ., data = df_pobreza)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0701 -0.2290 -0.1223 0.1639 2.8729
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.937e-02 4.607e-03 21.569 < 2e-16 ***
## mieperho 1.156e-01 1.365e-03 84.731 < 2e-16 ***
## percepho -3.720e-02 2.542e-03 -14.635 < 2e-16 ***
## inghog2d -5.453e-07 9.316e-08 -5.854 4.86e-09 ***
## gashog1d -1.132e-05 1.891e-07 -59.874 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.353 on 34485 degrees of freedom
## Multiple R-squared: 0.2832, Adjusted R-squared: 0.2831
## F-statistic: 3406 on 4 and 34485 DF, p-value: < 2.2e-16Luego realizamos las predicciones del modelo, donde la prediccion sera de “response” debido a que la variable es una variable de respuesta binaria.
pred_lm <- predict(fit_lm, type = "response")
summary(pred_lm)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.8729 0.1111 0.2097 0.2239 0.3455 1.7343Se ve que el minimo es negativo, lo cual no es aceptable y el maximo es mayor a 1 lo que implica que los valores no se ajustan adecuadamente para esta base de datos.
Por ello se aplica el modelo Probit:
fit_probit <- glm(pobreza ~. , data = df_pobreza,
family = binomial(link = "probit"))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(fit_probit)##
## Call:
## glm(formula = pobreza ~ ., family = binomial(link = "probit"),
## data = df_pobreza)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5.0633 -0.4415 -0.0731 0.0000 5.4157
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.203e+00 2.678e-02 -44.922 <2e-16 ***
## mieperho 7.751e-01 1.033e-02 75.041 <2e-16 ***
## percepho -1.981e-02 1.485e-02 -1.334 0.182
## inghog2d -2.225e-05 1.206e-06 -18.452 <2e-16 ***
## gashog1d -1.756e-04 2.818e-06 -62.319 <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: 36688 on 34489 degrees of freedom
## Residual deviance: 18228 on 34485 degrees of freedom
## AIC: 18238
##
## Number of Fisher Scoring iterations: 8Con ello realizamos las predicciones del modelo Probit:
pred_probit <- predict(fit_probit, type = "response")
summary(pred_probit)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.001228 0.073072 0.220130 0.339104 1.000000Luego similarmente, desarrollamos el modelo logit:
fit_logit <- glm(pobreza ~. , data = df_pobreza,
family = binomial(link = "logit"))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(fit_logit)##
## Call:
## glm(formula = pobreza ~ ., family = binomial(link = "logit"),
## data = df_pobreza)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0638 -0.4272 -0.1084 -0.0002 4.4697
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.012e+00 4.892e-02 -41.125 <2e-16 ***
## mieperho 1.431e+00 2.031e-02 70.466 <2e-16 ***
## percepho -5.740e-02 2.659e-02 -2.158 0.0309 *
## inghog2d -4.613e-05 2.273e-06 -20.298 <2e-16 ***
## gashog1d -3.244e-04 5.470e-06 -59.302 <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: 36688 on 34489 degrees of freedom
## Residual deviance: 18088 on 34485 degrees of freedom
## AIC: 18098
##
## Number of Fisher Scoring iterations: 8Nuevamente realizamos las predicciones del modelo, pero esta vez bajo la forma Logit.
pred_logit <- predict(fit_logit, type = "response")
summary(pred_logit)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00374 0.06902 0.22395 0.33871 1.00000Sus interpretaciones en lineas generales estan basadas en el signo, donde para el modelo Logit, mayores mimbros del hogar aumenta la probabilidad de ser pobre, mayores miembros de perceptores disminuye la probabilidad de ser pobre, mayores niveles de ingreso disminuye la probabilidad de ser pobre, y mayores niveles de gastos tambien disminuye la probabilidad de ser pobre, pero siendo menor a al impacto de los ingresos.
Paquete ROCR
ROCR ha existido durante casi 14 años y ha sido un caballo de batalla sólido como una roca para dibujar curvas ROC. Esto no solo es tranquilizadoramente transparente, sino que muestra la flexibilidad para calcular casi todas las medidas de rendimiento para un clasificador binario ingresando el parámetro apropiado.
Su uso se basa en determinar los valores de la curva ROC para los modelos probabilisticos, siendo el primer paso el de comparar los resultados predichos con los reales:
fit_logit <- glm(pobreza ~. , data = df_pobreza,
family = binomial(link = "logit"))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(fit_logit)##
## Call:
## glm(formula = pobreza ~ ., family = binomial(link = "logit"),
## data = df_pobreza)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0638 -0.4272 -0.1084 -0.0002 4.4697
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.012e+00 4.892e-02 -41.125 <2e-16 ***
## mieperho 1.431e+00 2.031e-02 70.466 <2e-16 ***
## percepho -5.740e-02 2.659e-02 -2.158 0.0309 *
## inghog2d -4.613e-05 2.273e-06 -20.298 <2e-16 ***
## gashog1d -3.244e-04 5.470e-06 -59.302 <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: 36688 on 34489 degrees of freedom
## Residual deviance: 18088 on 34485 degrees of freedom
## AIC: 18098
##
## Number of Fisher Scoring iterations: 8library(ROCR)## Warning: package 'ROCR' was built under R version 4.0.5pred_logit <- predict(fit_logit, df_pobreza, type = "response")
hist(pred_logit)
Sin emabrgo debemos conocer el umbral de probabilidad para determinar si es pobre o no. COn el paquete ROCR esto se puede calcular de manera simple y entendible. Lo primero es calcular los valores de la prediccion bajo el formato ROC con la funcion “prediction”, para luego evaluar el umbral en funcion al accuracy del modelo que a la vez arrogara los cortes o umbrales implicitos:
pred_logit <- prediction(pred_logit, df_pobreza$pobreza)
eval <- performance(pred_logit, "acc")
plot(eval)
Luego se puede identificar el maximo umbral que brinda el maximo accuracy del modelo logistico:
ind_max <-which.max(slot(eval, "y.values")[[1]])
acc_max <- slot(eval, "y.values")[[1]][ind_max]
cut_max <- slot(eval, "x.values")[[1]][ind_max]
paste("El maximo accuracy del modelo es",round(acc_max,3), "con un umbral de corte del", round(cut_max,3), "de probabilidad de ser pobre")## [1] "El maximo accuracy del modelo es 0.885 con un umbral de corte del 0.376 de probabilidad de ser pobre"Viendo si clasificamos correctamente el modelo: Donde: * “tpr” es True Positive Rate * “fpr” es False Positve Rate
roc_logit <- performance(pred_logit,"tpr", "fpr")
plot(roc_logit, colorize = TRUE, main = "ROC Curve", ylab = "Sensitivity",
xlab = "1- Specificity")
Finalmente el paquete ROCR tambien puede hallar el area bajo la curva, AUC por sus siglas en ingles del modelo logit que se planteo con la encuesta ENAHO 2020. El indicador AUC nos indica la capacidad que tiene el modelo de comportarse como un todo, mientras mas mayor sea el AUC mucho mejor, lo que asegura que el modelo tiene mas true positive que false negative.
auc_logit <- performance(pred_logit, "auc")
auc_logit <- unlist(slot(auc_logit, "y.values"))
paste("The final AUC del modelo es", round(auc_logit, 3))## [1] "The final AUC del modelo es 0.941"Paquete MLMetrics
Una colección de métricas de evaluación, incluyendo pérdida, puntuación y funciones de utilidad, que miden el rendimiento de regresión y clasificación.
Este paquete usualmente es usado para comparar medidas de rendimiento por cada modelos que se presenta con las mediciones del MSE, MAE, Accuracy, Gini, AUC, Confusion Matrix,
library(MLmetrics)##
## Attaching package: 'MLmetrics'## The following object is masked from 'package:base':
##
## Recallpred_logit <- predict(fit_logit, df_pobreza, type = "response")
MAE_logit <- MAE(pred_logit, df_pobreza$pobreza)
MAE_probit <- MAE(pred_probit, df_pobreza$pobreza)
MAE_logit## [1] 0.1636656MAE_probit## [1] 0.167004Por el lado del MAE, el modelo Logit es mejor que el probit.
MSE_logit <- MSE(pred_logit, df_pobreza$pobreza)
MSE_probit <- MSE(pred_probit, df_pobreza$pobreza)
MSE_logit## [1] 0.0815586MSE_probit## [1] 0.08219999Por el lado del MSE, el logit sigue siendo mejor que el Probit.
El paquete es usado tambien para la prediccion de modelos de Machine Learning, bajo un modelamiento partido en set de entrenamiento y sets de testeo.