library("stats")
library("psych")
library("readxl")
library("MASS")
library("ISLR")
library("fRegression")
library("vcd")
library("openxlsx")

Base1 = read.xlsx("test2.xlsx")

# race = cateogria de raza
# captial gain = si tiene o no ingresos
# hours-per-week = horas de trabajo por semana

Base2 = read.xlsx("URBI.xlsx")
# y = si se dio o no lo que se esta analizando 
# z = un numero que esta atado a la y (si se les dio o no (1 y 0))

Modelo Logit

modelo_logit1 = glm(capitalgain~hoursweek+race,data=Base1,family=binomial(link="logit"))
summary(modelo_logit1)

## 
## Call:
## glm(formula = capitalgain ~ hoursweek + race, family = binomial(link = "logit"), 
##     data = Base1)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.605073   0.526591  -6.846 7.59e-12 ***
## hoursweek    0.027962   0.009424   2.967  0.00301 ** 
## race        -0.154478   0.224193  -0.689  0.49080    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 456.40  on 898  degrees of freedom
## Residual deviance: 447.46  on 896  degrees of freedom
## AIC: 453.46
## 
## Number of Fisher Scoring iterations: 5

# La cantidad de horas si tiene significanciía en si generan ingresos o no. Mientras que la raza no tiene significancia en sus ingresos.

modelo_logit2 = glm(Y~Z,data=Base2,family=binomial(link="logit"))
summary(modelo_logit2)

## 
## Call:
## glm(formula = Y ~ Z, family = binomial(link = "logit"), data = Base2)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.0703     0.5004  -2.139 0.032449 *  
## Z            -0.4428     0.1170  -3.786 0.000153 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 62.183  on 44  degrees of freedom
## Residual deviance: 30.669  on 43  degrees of freedom
## AIC: 34.669
## 
## Number of Fisher Scoring iterations: 6

Modelo Probit

modelo_probit1 = glm(capitalgain~hoursweek+race,data=Base1,family=binomial(link="probit"))
summary(modelo_probit1)

## 
## Call:
## glm(formula = capitalgain ~ hoursweek + race, family = binomial(link = "probit"), 
##     data = Base1)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.003983   0.260522  -7.692 1.45e-14 ***
## hoursweek    0.014436   0.004933   2.926  0.00343 ** 
## race        -0.074307   0.103155  -0.720  0.47131    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 456.4  on 898  degrees of freedom
## Residual deviance: 447.2  on 896  degrees of freedom
## AIC: 453.2
## 
## Number of Fisher Scoring iterations: 5

modelo_probit2 = glm(Y~Z,data=Base2,family=binomial(link="probit"))
summary(modelo_probit2)

## 
## Call:
## glm(formula = Y ~ Z, family = binomial(link = "probit"), data = Base2)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.67883    0.27671  -2.453   0.0142 *  
## Z           -0.25095    0.05836  -4.300 1.71e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 62.183  on 44  degrees of freedom
## Residual deviance: 30.835  on 43  degrees of freedom
## AIC: 34.835
## 
## Number of Fisher Scoring iterations: 6

Criterios de Información

Estimamos los criterios de información para cada modelo utilizado en la Base1

CIA_Logit1 = AIC(modelo_logit1)
CIA_Probit1 = AIC(modelo_probit1)
CIA_Logit1

## [1] 453.4582

CIA_Probit1

## [1] 453.2016

El modelo con mejor ajuste con base al CIA es el probit al tener un valor menor.

Ahora estimamos los criterios de información para cada modelo utilizado en la Base2.

CIA_Logit2 = AIC(modelo_logit2)
CIA_Probit2 = AIC(modelo_probit2)
CIA_Logit2

## [1] 34.66863

CIA_Probit2

## [1] 34.83505

Aquí el modelo con mejor ajuste con base al CIA es el logit al tener un valor menor.

El objetivo de hacer esto es predecir

¿Cómo generar predicciones?

log.odds <- predict(modelo_probit1, data.frame(hoursweek = 99, race = 2))
Prob1 = pnorm(log.odds)
Prob1

##         1 
## 0.2347176

# Equivalentemente
predict(modelo_probit1, data.frame(hoursweek = 3, race = 1), type="response")

##          1 
## 0.02092637

### lo que arroja es la probabilidad de que tengan ingresos, ya que el numero que arroja esta más cerca al 0.

log.odds2 <- predict(modelo_logit2, data.frame(Z = mean(Base2$Z)))
Prob2 = exp(log.odds2)/(1+exp(log.odds2))
Prob2

##         1 
## 0.6753464

#Equivalentemente
predict(modelo_logit2, data.frame(Z = mean(Base2$Z)), type="response")

##         1 
## 0.6753464

### el 66% de los datos son rechazados y lo restante es lo aceptado --> esta mas cerca al 1

Interpretación de Coeficientes

En este modelo los coeficientes indican cómo una unidad de esta unidad afecta al logaritmo de cociente de probabilidades. Lo más relevante en sí es el signo de esta relación para así saber si el incremento de la variable en cuestión aumenta o disminuye la probabilidad de tener el atributo que se está prediciendo.

Prob1*coef(modelo_probit1)

##  (Intercept)    hoursweek         race 
## -0.470370191  0.003388467 -0.017441175

# B0x seria hoursweek / B1 seria race / intercept =b2

Prob2*coef(modelo_logit2)

## (Intercept)           Z 
##  -0.7228405  -0.2990640

Matriz de confusiones

La Matriz de Confusiones es la herramienta principal para generar diferentes métricas que sirven para evaluar el desmepeño de nuestro modelo. Estas métricas las veremos con mayor detalle en las próximas clases.

predicciones1 <- ifelse(test = modelo_probit1$fitted.values > 0.5, yes =1, no = 0)
matriz_confusion1 <- table(predicciones1,modelo_probit1$model$capitalgain,dnn = c( "predicciones","observaciones"))
matriz_confusion1

##             observaciones
## predicciones   0   1
##            0 836  63

# si es mayor a 0.5 es si y si es menor a 0.5 es no

predicciones2 <- ifelse(test = modelo_logit2$fitted.values > 0.5, yes =1, no = 0)
matriz_confusion2 <- table(modelo_logit2$model$Y, predicciones2,dnn = c("observaciones", "predicciones"))
matriz_confusion2

##              predicciones
## observaciones  0  1
##             0 19  2
##             1  3 21

A partir de esta matriz se puede estimar una medida del ajuste

(836+0)/(899)*100

## [1] 92.99221

# (calificacion buena )/ total de calificaciones * 100

(19+21)/(21+24)*100

## [1] 88.88889