Lab 3: Probit, Logit

Los datos que utilizaremos para nuestro ejemplo son los datos de Herriges & Kling (1999) usados y provistos por Cameron & Trivedi (2005).

Exploramos los datos:

Las observaciones son clientes
mode indica la elección del lugar donde pescó
El ingreso de la persona es income
La tasa de pesca es q
El precio es p

Empezamos como siempre cargando nuestras librerias:

rm(list=ls())
setwd("C:/DAVE2/CIDE/Doc/Econometria II/Lab 3")
options(scipen=999) 

library(tidyverse)

## -- Attaching packages ---------------------------------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.2
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0

## -- Conflicts ------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(janitor)

## 
## Attaching package: 'janitor'

## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(modelsummary)

Importamos los datos:

data_fishing <- read_csv("./fishing_data.csv", locale = locale(encoding = "latin1")) %>% clean_names()

## Parsed with column specification:
## cols(
##   client = col_double(),
##   mode = col_double(),
##   price = col_double(),
##   crate = col_double(),
##   dbeach = col_double(),
##   dpier = col_double(),
##   dprivate = col_double(),
##   dcharter = col_double(),
##   pbeach = col_double(),
##   ppier = col_double(),
##   pprivate = col_double(),
##   pcharter = col_double(),
##   qbeach = col_double(),
##   qpier = col_double(),
##   qprivate = col_double(),
##   qcharter = col_double(),
##   income = col_double()
## )

Hacemos una primera analisis explorativa del dataframe:

head(data_fishing, n=5)

## # A tibble: 5 x 17
##   client  mode price  crate dbeach dpier dprivate dcharter pbeach ppier pprivate
##    <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>  <dbl> <dbl>    <dbl>
## 1      1     4 183.  0.539       0     0        0        1  158.  158.     158. 
## 2      2     4  34.5 0.467       0     0        0        1   15.1  15.1     10.5
## 3      3     3  24.3 0.241       0     0        1        0  162.  162.      24.3
## 4      4     2  15.1 0.0789      0     1        0        0   15.1  15.1     55.9
## 5      5     3  41.5 0.108       0     0        1        0  107.  107.      41.5
## # ... with 6 more variables: pcharter <dbl>, qbeach <dbl>, qpier <dbl>,
## #   qprivate <dbl>, qcharter <dbl>, income <dbl>

Podemos ver los nombres que caracterizan nuestra base:

names(data_fishing)

##  [1] "client"   "mode"     "price"    "crate"    "dbeach"   "dpier"   
##  [7] "dprivate" "dcharter" "pbeach"   "ppier"    "pprivate" "pcharter"
## [13] "qbeach"   "qpier"    "qprivate" "qcharter" "income"

Escogemos nada màs las modalidad barco y puerto para nuestro analisis con probit y logit:

data_binary<-data_fishing %>% 
  filter(mode==2 | mode==4) %>% 
  mutate(charter=ifelse(mode==4,1,0)) %>% 
  mutate(pratio=100*log(pcharter/ppier),
         lnrelp=log(pcharter/ppier))

Para la estimacion del probit utilizamos la funcion glm:

mprobit <- glm(formula = charter ~ lnrelp,
               family = binomial(link = "probit"), 
               data = data_binary)

Obtenemos los coeficientes:

summary(mprobit)$coef

##              Estimate Std. Error   z value
## (Intercept)  1.194358 0.08814144  13.55048
## lnrelp      -1.055513 0.07542127 -13.99490
##                                                         Pr(>|z|)
## (Intercept) 0.00000000000000000000000000000000000000000787176739
## lnrelp      0.00000000000000000000000000000000000000000001674672

Poemos hacer lo mismo por el logit:

mlogit <- glm(charter ~ lnrelp,
              family = binomial(link = "logit"), 
              data = data_binary)

Como serìa un modelo de probabilidad lineal¿

Solución

mlineal <- lm(charter ~ lnrelp,
              data=data_binary)

Vamos ahora a crear una colecciòn de modelos para efectuar la comparaciòn:

models<-list('Logit'=mlogit, 'Probit'=mprobit, 'Lineal'=mlineal)

La lista cm nos permite controlar los nombres e los coeficientes:

cm <- c( '(Intercept)' = 'Constante', 'lnrelp' = 'ln relp')

Gm controlarà als medidad de ajuste

gm <- modelsummary::gof_map

Pedimos ahora que omita todo:

gm$omit <- TRUE

Y que despliegue solo el log de la funcion de verosimilitud:

gm <- gm %>% 
  mutate(omit=ifelse(raw=="logLik", F, omit))

la tabla final resultante serà:

msummary(models,
         statistic='statistic',
         coef_map = cm,
         gof_map = gm)

	Logit	Probit	Lineal
Constante	2.053	1.194	0.784
	(12.154)	(13.550)	(58.214)
ln relp	-1.823	-1.056	-0.243
	(-12.607)	(-13.995)	(-23.284)
Log.Lik.	-206.827	-204.411	-195.167

Podemos ahora hacer predicciones con todos los modelos:

plogit <- mlogit %>% predict(data=data_binary, type = "response")

pprobit <- mprobit %>% predict(data=data_binary, type = "response")

plineal <- mlineal %>% predict(data=data_binary, type = "response")

Ponemos ahora los resultados en un dataframe:

data_for_plot <- data.frame(lnrelp=data_binary$lnrelp,
                            plogit,
                            pprobit,
                            plineal)

Y podemos ajustar los datos para que esten en un formato adecuado

data_for_plot <- pivot_longer(data_for_plot,
                              cols= c("plogit","pprobit","plineal"),
                              names_to="Modelo",
                              values_to = "prob")

Graficamos ahora las probabilidades predictas y ajustadas:

data_for_plot %>% 
  ggplot()+
  geom_line(aes(x=lnrelp,y=prob, color=Modelo))

Ahora estamos interesados en estudiar una variable dependiente binaria/dicotomica/dummy

Vamos cargando mas librerias:

pacman::p_load(ggplot2,tidyverse,knitr,wooldridge,stargazer,
               xtable,stringr,tibble,sjPlot,sjmisc,jtools,interactions)

Analizamos la base de datos:

data("mroz")
names(mroz)

##  [1] "inlf"     "hours"    "kidslt6"  "kidsge6"  "age"      "educ"    
##  [7] "wage"     "repwage"  "hushrs"   "husage"   "huseduc"  "huswage" 
## [13] "faminc"   "mtr"      "motheduc" "fatheduc" "unem"     "city"    
## [19] "exper"    "nwifeinc" "lwage"    "expersq"

Describe la participacion en la fuerza laboral de las mujeres casadas

La variable inlf es igual a 1 si la mujer estaba trabajando en 1975

Obtenemos un primer modelo lineal:

mpl <- lm(inlf ~ nwifeinc + educ + exper + expersq
          + age + kidslt6 + kidsge6, data = mroz)
stargazer(mpl,type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                inlf            
## -----------------------------------------------
## nwifeinc                     -0.003**          
##                               (0.001)          
##                                                
## educ                         0.038***          
##                               (0.007)          
##                                                
## exper                        0.039***          
##                               (0.006)          
##                                                
## expersq                      -0.001***         
##                              (0.0002)          
##                                                
## age                          -0.016***         
##                               (0.002)          
##                                                
## kidslt6                      -0.262***         
##                               (0.034)          
##                                                
## kidsge6                        0.013           
##                               (0.013)          
##                                                
## Constant                     0.586***          
##                               (0.154)          
##                                                
## -----------------------------------------------
## Observations                    753            
## R2                             0.264           
## Adjusted R2                    0.257           
## Residual Std. Error      0.427 (df = 745)      
## F Statistic           38.218*** (df = 7; 745)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Hagamos otro ejemplo

library(AER)

## Loading required package: car

## Loading required package: carData

## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

## Loading required package: lmtest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

## Loading required package: survival

data(HMDA)

Usaremos datos relacionados con las aplicaciones a hipotecas en Boston en 1990

names(HMDA)

##  [1] "deny"      "pirat"     "hirat"     "lvrat"     "chist"     "mhist"    
##  [7] "phist"     "unemp"     "selfemp"   "insurance" "condomin"  "afam"     
## [13] "single"    "hschool"

Nos interesa la variable deny que nos indica si les rechazaron la aplicacion al a hipoteca

La variable pirat es el tamano que el pago anticipado mensual del prestamo relativo al ingreso del aplicante

Tenemos que cambiar deny a numerica:

as.numeric(HMDA$deny)

##    [1] 1 1 1 1 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##   [38] 1 1 1 1 1 2 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##   [75] 1 1 1 1 2 1 1 1 2 2 1 1 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1
##  [112] 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [149] 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1
##  [186] 1 1 1 1 1 2 2 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1
##  [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 2 2 1 1 1 1 1
##  [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 2
##  [371] 2 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1
##  [408] 1 1 1 2 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [445] 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1
##  [482] 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1
##  [519] 1 1 1 1 1 2 1 1 2 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 2
##  [556] 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 1 2 1 1 1 1 1
##  [593] 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1
##  [630] 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [667] 1 1 1 1 1 1 2 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2
##  [704] 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [741] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [778] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [815] 2 2 2 2 2 2 2 1 2 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1
##  [852] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [889] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1
##  [963] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1000] 1 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 1 1 1 2 2 2 1 2 2 2 1 1 2 1 1 1 1 1 1
## [1037] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 2 1 1 1 2 1 1 1 1
## [1074] 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1
## [1111] 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1148] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1
## [1185] 1 1 1 2 1 1 1 1 1 1 1 1 1 2 2 2 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1222] 1 1 1 1 1 1 2 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1259] 1 1 1 1 1 2 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 2 1 2 2 1 1 1 1 1 1 1 1 1 1 1
## [1296] 2 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1333] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1370] 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1
## [1407] 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1444] 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1481] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1
## [1518] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1
## [1555] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1
## [1592] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1629] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1666] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1703] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1740] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [1777] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2
## [1814] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 2 1
## [1851] 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 2 2 1 2 1 2 2
## [1888] 1 1 1 2 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 2 1 1 1 2 1 2 1 1 1 2 1 1 1 1 2 1 1
## [1925] 2 2 1 2 2 1 1 2 1 1 2 2 2 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 2 1 1 2 1
## [1962] 1 1 1 2 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1
## [1999] 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [2036] 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1
## [2073] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1
## [2110] 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 2 1 2 2 1 1 2 1 1 1 1
## [2147] 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1
## [2184] 1 2 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1
## [2221] 1 1 2 2 1 2 1 1 2 2 2 2 1 2 1 1 1 2 2 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1
## [2258] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1
## [2295] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [2332] 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1 1
## [2369] 1 1 1 1 1 1 1 1 1 1 2 2

Lo que pasa se convierte no -> deny = 1 yes -> deny = 2

Entonces tenemos que quitar un uno para obtener una dummie:

HMDA$deny <- as.numeric(HMDA$deny)-1

Aplicamos de nuevo un modelo lineal:

mpl2 <- lm(deny ~ pirat, data = HMDA)
mpl2

## 
## Call:
## lm(formula = deny ~ pirat, data = HMDA)
## 
## Coefficients:
## (Intercept)        pirat  
##    -0.07991      0.60353

Graficamos los que acabamos de encontrar:

plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Negativa Hipoteca vs. ratio Pago/Ingreso",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.8)

Agreguemos algunos elementos graficos:

plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Negativa Hipoteca vs. ratio Pago/Ingreso",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.8)
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Hipoteca denegada")
text(2.5, -0.1, cex= 0.8, "Hipoteca aprobada")

Agregamos la linea de regresion

plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Negativa Hipoteca vs. ratio Pago/Ingreso",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.8)
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Hipoteca denegada")
text(2.5, -0.1, cex= 0.8, "Hipoteca aprobada")
abline(mpl2, 
       lwd = 1.8, 
       col = "steelblue")

¿como obtendrian los errores estandar robustos? coeftest(x, vcov. = NULL, df = NULL, type) vcovHC = Heteroskedasticity-Consistent Covariance Matrix Estimation

coeftest(mpl2, vcov. = vcovHC, type = "HC1")

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value       Pr(>|t|)    
## (Intercept) -0.079910   0.031967 -2.4998        0.01249 *  
## pirat        0.603535   0.098483  6.1283 0.000000001036 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vamor ahora a renombrar la variable afroamericano

colnames(HMDA)[colnames(HMDA) == "afam"] <- "black"

Agregamos la variable a la regresion

mpl3 <- lm(deny ~ pirat + black, data = HMDA)
coeftest(mpl3, vcov. = vcovHC)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value          Pr(>|t|)    
## (Intercept) -0.090514   0.033430 -2.7076          0.006826 ** 
## pirat        0.559195   0.103671  5.3939 0.000000075750142 ***
## blackyes     0.177428   0.025055  7.0815 0.000000000001871 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Porcentaje correctamente predicho

y_hat <- predict(mpl3, HMDA, type="response") 
summarise(HMDA, min(y_hat), max(y_hat), mean(y_hat))

##   min(y_hat) max(y_hat) mean(y_hat)
## 1 -0.0905136    1.58707   0.1197479

Vamos a ver ahora a analizar los resultados obtenidos:

HMDA <- mutate(HMDA, y_pred_ad = case_when(y_hat >= 0.50 ~ "1",
                                   y_hat < 0.50 ~ "0" ))
HMDA$y_pred_ad <- as.numeric(HMDA$y_pred_ad)
summarise(HMDA, min(y_pred_ad), max(y_pred_ad), mean(y_pred_ad))

##   min(y_pred_ad) max(y_pred_ad) mean(y_pred_ad)
## 1              0              1     0.003361345

HMDA <- mutate(HMDA, corr_pred = if_else(y_pred_ad == deny, 1, 0))

num <- nrow(HMDA) 

summarise(HMDA, correct = sum(corr_pred)/num )

##     correct
## 1 0.8819328

Ahora calculamos con un probit

probit1 <- glm(deny ~ pirat, 
                  family = binomial(link = "probit"), 
                  data = HMDA)

Vamos a calcular errores robustos:

coeftest(probit1, vcov. = vcovHC, type = "HC1")

## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value              Pr(>|z|)    
## (Intercept) -2.19415    0.18901 -11.6087 < 0.00000000000000022 ***
## pirat        2.96787    0.53698   5.5269         0.00000003259 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hagamos otra grafica

plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Modelo Probit Prob. Hipoteca Denegada vs. P/I Ratio",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.85)


abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Hipoteca Denegada")
text(2.5, -0.1, cex= 0.8, "Hipoteca Aprobada")

x <- seq(0, 3, 0.01)
y <- predict(probit1, list(pirat = x), type = "response")

lines(x, y, lwd = 1.5, col = "steelblue")

¿Que diferencias hay con la grafica anterior?

Solución

¿Como vemos el efecto de pasar de un pirat de 0.3 a 0.4?r

Solución

predictions <- predict(probit1, 
                       newdata = data.frame("pirat" = c(0.3, 0.4)),
                       type = "response")
predictions

##          1          2 
## 0.09615344 0.15696777

Calculamos la diferencia en las predicciones

diff(predictions)

##          2 
## 0.06081433

El efecto es de un aumento de la probabilidad de 6.2%

probit2 <- glm(deny ~ pirat + black, 
                   family = binomial(link = "probit"), 
                   data = HMDA)

coeftest(probit2, vcov. = vcovHC, type = "HC1")

## 
## z test of coefficients:
## 
##              Estimate Std. Error  z value              Pr(>|z|)    
## (Intercept) -2.258787   0.176608 -12.7898 < 0.00000000000000022 ***
## pirat        2.741779   0.497673   5.5092         0.00000003605 ***
## blackyes     0.708155   0.083091   8.5227 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

¿Cual seria el efecto estimado para dos personas con un pirat = 0.3 pero con diferencia ser afroamericano?

Solución

Predicciones P/I ratio = 0.3

  predictions <- predict(probit2, 
                       newdata = data.frame("black" = c("no", "yes"), 
                                            "pirat" = c(0.3, 0.3)),
                       type = "response")

Diferencia en Probabilidades

diff(predictions)

##         2 
## 0.1578117

La diferencia es de aproximadamente 15.8%

Logit

Estimamos el modelo logit

logit1 <- glm(deny ~ pirat, 
                 family = binomial(link = "logit"), 
                 data = HMDA)

coeftest(logit1, vcov. = vcovHC, type = "HC1")

## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value              Pr(>|z|)    
## (Intercept) -4.02843    0.35898 -11.2218 < 0.00000000000000022 ***
## pirat        5.88450    1.00015   5.8836        0.000000004014 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Y hacemos otra grafica

plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Modelos Probit y Logit sobre la hipoteca",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.9)
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Hipoteca Denegada")
text(2.5, -0.1, cex= 0.8, "Hipoteca Aprobada")

Agregamos los modelos probit y logit

plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Modelos Probit y Logit sobre la hipoteca",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.9)
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Hipoteca Denegada")
text(2.5, -0.1, cex= 0.8, "Hipoteca Aprobada")

x<-seq(0, 3, 0.01)
y_probit <- predict(probit1, list(pirat = x), type = "response")
y_logit <- predict(logit1, list(pirat = x), type = "response")

lines(x, y_probit, lwd = 1.5, col = "steelblue")
lines(x, y_logit, lwd = 1.5, col = "black", lty = 2)

Anadimos leyendas

plot(x = HMDA$pirat, 
     y = HMDA$deny,
     main = "Modelos Probit y Logit sobre la hipoteca",
     xlab = "P/I ratio",
     ylab = "Deny",
     pch = 20,
     ylim = c(-0.4, 1.4),
     cex.main = 0.9)
abline(h = 1, lty = 2, col = "darkred")
abline(h = 0, lty = 2, col = "darkred")
text(2.5, 0.9, cex = 0.8, "Hipoteca Denegada")
text(2.5, -0.1, cex= 0.8, "Hipoteca Aprobada")

x<-seq(0, 3, 0.01)
y_probit <- predict(probit1, list(pirat = x), type = "response")
y_logit <- predict(logit1, list(pirat = x), type = "response")

lines(x, y_probit, lwd = 1.5, col = "steelblue")
lines(x, y_logit, lwd = 1.5, col = "black", lty = 2)
legend("topleft",
       horiz = TRUE,
       legend = c("Probit", "Logit"),
       col = c("steelblue", "black"), 
       lty = c(1, 2))

Calculamos un segundo modelo logit

logit2 <- glm(deny ~ pirat + black, 
                  family = binomial(link = "logit"), 
                  data = HMDA)

coeftest(logit2, vcov. = vcovHC, type = "HC1")

## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value              Pr(>|z|)    
## (Intercept) -4.12556    0.34597 -11.9245 < 0.00000000000000022 ***
## pirat        5.37036    0.96376   5.5723         0.00000002514 ***
## blackyes     1.27278    0.14616   8.7081 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

¿Cual es la diferencia por raza para obtener la hipotecadado un valor de pirat = 0.3? pr

Solución

predictions <- predict(logit2, 
                       newdata = data.frame("black" = c("no", "yes"), 
                                            "pirat" = c(0.3, 0.3)),
                       type = "response")
predictions

##          1          2 
## 0.07485143 0.22414592

Prob aplicante blanco sea denegado 7.5% Prob afroamericano sea denegado 22.4%

diff(predictions)

##         2 
## 0.1492945

#Fuente: https://www.econometrics-with-r.org/11-2-palr.html

Efectos Marginales

Hasta ahora hemos comparado el efecto para distintos individuos

Una pregunta que nos falta responder es el efecto marginal de x_j

library("margins")

Efecto Marginal Promedio

m <- margins(probit2)
summary(m)

##    factor    AME     SE      z      p  lower  upper
##  blackyes 0.1696 0.0242 6.9955 0.0000 0.1221 0.2171
##     pirat 0.5014 0.0690 7.2657 0.0000 0.3662 0.6367

plot(m)

En STATA margins tiene las opciones * at: calcula los efectos marginales en un valor especifico (potencialmente representativo) * atmeans: calcula los efectos marginales en la media de los datos * over: calcula los efectos marginales en subsets de los datos por ejemplo, en hombres y en mujeres

¿Como lo hacemos en R?

library("mfx")

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:wooldridge':
## 
##     cement

## The following object is masked from 'package:dplyr':
## 
##     select

## Loading required package: betareg

probitmfx(formula, data, atmean = TRUE, robust = FALSE, clustervar1 = NULL, clustervar2 = NULL, start = NULL, control = list())

con atmean=TRUE tienen el efecto parcial en el promedio si le ponen atmean=FALSE tienen el efecto parcial promedio

probitmfx(deny ~ pirat + black, data = HMDA, atmean = TRUE, robust = FALSE, clustervar1 = NULL,
          clustervar2 = NULL, start = NULL, control = list())

## Call:
## probitmfx(formula = deny ~ pirat + black, data = HMDA, atmean = TRUE, 
##     robust = FALSE, clustervar1 = NULL, clustervar2 = NULL, start = NULL, 
##     control = list())
## 
## Marginal Effects:
##             dF/dx Std. Err.      z              P>|z|    
## pirat    0.500219  0.069456 7.2019 0.0000000000005936 ***
## blackyes 0.171688  0.024695 6.9523 0.0000000000035949 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## dF/dx is for discrete change for the following variables:
## 
## [1] "blackyes"

probitmfx(deny ~ pirat + black, data = HMDA, atmean = FALSE, robust = FALSE, clustervar1 = NULL,
          clustervar2 = NULL, start = NULL, control = list())

## Call:
## probitmfx(formula = deny ~ pirat + black, data = HMDA, atmean = FALSE, 
##     robust = FALSE, clustervar1 = NULL, clustervar2 = NULL, start = NULL, 
##     control = list())
## 
## Marginal Effects:
##             dF/dx Std. Err.      z              P>|z|    
## pirat    0.501420  0.069012 7.2657 0.0000000000003711 ***
## blackyes 0.169592  0.024243 6.9955 0.0000000000026440 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## dF/dx is for discrete change for the following variables:
## 
## [1] "blackyes"

De igual forma se pueden calcular los efectos promedio para el logit.

logitmfx(formula, data, atmean = TRUE, robust = FALSE, clustervar1 = NULL, clustervar2 = NULL, start = NULL, control = list())

logitmfx(deny ~ pirat + black, data = HMDA, atmean = TRUE, robust = FALSE, clustervar1 = NULL, 
         clustervar2 = NULL, start = NULL, control = list())

## Call:
## logitmfx(formula = deny ~ pirat + black, data = HMDA, atmean = TRUE, 
##     robust = FALSE, clustervar1 = NULL, clustervar2 = NULL, start = NULL, 
##     control = list())
## 
## Marginal Effects:
##             dF/dx Std. Err.      z               P>|z|    
## pirat    0.494854  0.066287 7.4653 0.00000000000008312 ***
## blackyes 0.167080  0.024410 6.8447 0.00000000000766454 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## dF/dx is for discrete change for the following variables:
## 
## [1] "blackyes"

Odds Ratio

logitor(deny ~ pirat + black, data = HMDA, robust = FALSE, clustervar1 = NULL, 
        clustervar2 = NULL, start = NULL, control = list())

## Call:
## logitor(formula = deny ~ pirat + black, data = HMDA, robust = FALSE, 
##     clustervar1 = NULL, clustervar2 = NULL, start = NULL, control = list())
## 
## Odds Ratio:
##          OddsRatio Std. Err.      z                 P>|z|    
## pirat    214.94067 156.54450 7.3737     0.000000000000166 ***
## blackyes   3.57077   0.52204 8.7059 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ejemplo 2

La regresión probit, también llamada modelo probit, se utiliza para modelar variables de resultado dicotómicas o binarias. En el modelo probit, la distribución normal estándar inversa de la probabilidad se modela como una combinación lineal de los predictores.

Ejemplo 1

Supongamos que estamos interesados en los factores que influyen en si un candidato político gana una elección. La variable de resultado (respuesta) es binaria (0/1); ganar o perder. Las variables predictoras de interés son la cantidad de dinero gastada en la campaña, la cantidad de tiempo dedicado a la campaña de manera negativa y si el candidato es un titular.

Ejemplo 2

Un investigador está interesado en cómo las variables, como GRE (calificaciones del examen de registro de graduados), GPA (promedio de calificaciones) y el prestigio de la institución de pregrado, afectan la admisión a la escuela de posgrado. La variable de respuesta, admitir / no admitir, es una variable binaria.

Hacemos un analisis de este ultimo ejemplo, descargamos los datos y analizamolos:

mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")

## convert rank to a factor (categorical variable)
mydata$rank <- factor(mydata$rank)

## view first few rows
head(mydata)

##   admit gre  gpa rank
## 1     0 380 3.61    3
## 2     1 660 3.67    3
## 3     1 800 4.00    1
## 4     1 640 3.19    4
## 5     0 520 2.93    4
## 6     1 760 3.00    2

Tenemos admit como variable de respuesta binaria. gre y gpa son variables continuas mmientras *rank" es una categorica.

summary(mydata)

##      admit             gre             gpa        rank   
##  Min.   :0.0000   Min.   :220.0   Min.   :2.260   1: 61  
##  1st Qu.:0.0000   1st Qu.:520.0   1st Qu.:3.130   2:151  
##  Median :0.0000   Median :580.0   Median :3.395   3:121  
##  Mean   :0.3175   Mean   :587.7   Mean   :3.390   4: 67  
##  3rd Qu.:1.0000   3rd Qu.:660.0   3rd Qu.:3.670          
##  Max.   :1.0000   Max.   :800.0   Max.   :4.000

xtabs(~rank + admit, data = mydata)

##     admit
## rank  0  1
##    1 28 33
##    2 97 54
##    3 93 28
##    4 55 12

A continuación se muestra una lista de algunos métodos de análisis que puede haber encontrado. Algunos de los métodos enumerados son bastante razonables, mientras que otros han caído en desgracia o tienen limitaciones.

Regresión probit, el tema del ejemplo
Regresión logística. Un modelo logit producirá resultados de regresión probit similares. La elección de probit frente a logit depende en gran medida de las preferencias individuales.
Regresión OLS. Cuando se usa con una variable de respuesta binaria, este modelo se conoce como modelo de probabilidad lineal y se puede usar como una forma de describir probabilidades condicionales. Sin embargo, los errores (es decir, residuales) del modelo de probabilidad lineal violan los supuestos de homocedasticidad y normalidad de errores de la regresión OLS, lo que da como resultado errores estándar no válidos y pruebas de hipótesis. Para una discusión más detallada de estos y otros problemas con el modelo de probabilidad lineal, vea Long (1997, p. 38-40).
Análisis de función discriminante de dos grupos. Un método multivariado para variables de resultado dicotómicas.
T2 de Hotelling. El resultado 0/1 se convierte en la variable de agrupación y los predictores anteriores se convierten en variables de resultado. Esto producirá una prueba general de significancia, pero no dará coeficientes individuales para cada variable, y no está claro hasta qué punto cada “predictor” se ajusta al impacto de los otros “predictores”.

AHora utilizamos el probit:

myprobit <- glm(admit ~ gre + gpa + rank, family = binomial(link = "probit"), 
    data = mydata)

## model summary
summary(myprobit)

## 
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = binomial(link = "probit"), 
##     data = mydata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6163  -0.8710  -0.6389   1.1560   2.1035  
## 
## Coefficients:
##              Estimate Std. Error z value  Pr(>|z|)    
## (Intercept) -2.386836   0.673946  -3.542  0.000398 ***
## gre          0.001376   0.000650   2.116  0.034329 *  
## gpa          0.477730   0.197197   2.423  0.015410 *  
## rank2       -0.415399   0.194977  -2.131  0.033130 *  
## rank3       -0.812138   0.208358  -3.898 0.0000971 ***
## rank4       -0.935899   0.245272  -3.816  0.000136 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.98  on 399  degrees of freedom
## Residual deviance: 458.41  on 394  degrees of freedom
## AIC: 470.41
## 
## Number of Fisher Scoring iterations: 4

Podemos utilizar la funciòn confint para obtener los intervalos e confianza:

confint(myprobit)

## Waiting for profiling to be done...

##                     2.5 %       97.5 %
## (Intercept) -3.7201050682 -1.076327713
## gre          0.0001104101  0.002655157
## gpa          0.0960654793  0.862610221
## rank2       -0.7992113929 -0.032995019
## rank3       -1.2230955861 -0.405008112
## rank4       -1.4234218227 -0.459538829

Podemos testar un efecto general del rank utilizando el test de Wald

library(aod)

## Warning: package 'aod' was built under R version 4.0.5

## 
## Attaching package: 'aod'

## The following object is masked from 'package:survival':
## 
##     rats

wald.test(b = coef(myprobit), Sigma = vcov(myprobit), Terms = 4:6)

## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 21.4, df = 3, P(> X2) = 0.000089

El estadistico asociado al test chi-cuado es e Que significa¿

Podemos testar diferentes hipotesis sobre la diferencia entra coeficiente estimados. Probamos si el coeficiente de rank=2 es igual al del rank=3

l <- cbind(0, 0, 0, 1, -1, 0)
wald.test(b = coef(myprobit), Sigma = vcov(myprobit), L = l)

## Wald test:
## ----------
## 
## Chi-squared test:
## X2 = 5.6, df = 1, P(> X2) = 0.018

Podemos tambièn predecir las probabiliaes para obtener una mayor comprension

newdata <- data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100), 
    4 * 4), gpa = rep(c(2.5, 3, 3.5, 4), each = 100 * 4), rank = factor(rep(rep(1:4, 
    each = 100), 4)))

head(newdata)

##        gre gpa rank
## 1 200.0000 2.5    1
## 2 206.0606 2.5    1
## 3 212.1212 2.5    1
## 4 218.1818 2.5    1
## 5 224.2424 2.5    1
## 6 230.3030 2.5    1

Ahora podmos calculare los valores/probabiliades asi como sus se.

newdata[, c("p", "se")] <- predict(myprobit, newdata, type = "response", se.fit = TRUE)[-3]

ggplot(newdata, aes(x = gre, y = p, colour = rank)) + geom_line() + facet_wrap(~gpa)

Tenemos diferentes medidas de resumen y precision.

summary(myprobit)

## 
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = binomial(link = "probit"), 
##     data = mydata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6163  -0.8710  -0.6389   1.1560   2.1035  
## 
## Coefficients:
##              Estimate Std. Error z value  Pr(>|z|)    
## (Intercept) -2.386836   0.673946  -3.542  0.000398 ***
## gre          0.001376   0.000650   2.116  0.034329 *  
## gpa          0.477730   0.197197   2.423  0.015410 *  
## rank2       -0.415399   0.194977  -2.131  0.033130 *  
## rank3       -0.812138   0.208358  -3.898 0.0000971 ***
## rank4       -0.935899   0.245272  -3.816  0.000136 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.98  on 399  degrees of freedom
## Residual deviance: 458.41  on 394  degrees of freedom
## AIC: 470.41
## 
## Number of Fisher Scoring iterations: 4

O calcular la diferencia en desviaciòn para los dos modelos. Para encontrar la diferencia en la desviación para los dos modelos (es decir, el estadístico de prueba), podemos calcular el cambio en la desviación y probarlo usando una prueba de chi cuadrado: el cambio en la desviación distribuido como chi cuadrado sobre el cambio en grados de libertad.

## change in deviance
with(myprobit, null.deviance - deviance)

## [1] 41.56335

## change in degrees of freedom
with(myprobit, df.null - df.residual)

## [1] 5

## chi square test p-value
with(myprobit, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE))

## [1] 0.00000007218932

Podemos calcular la log verosimilitud de la siguiente manera:

logLik(myprobit)

## 'log Lik.' -229.2066 (df=6)

Algunas consideraciones:

Celdas vacías o celdas pequeñas: debe verificar las celdas vacías o pequeñas haciendo una tabla de referencias cruzadas entre los predictores categóricos y la variable de resultado. Si una celda tiene muy pocos casos (una celda pequeña), el modelo puede volverse inestable o puede que no funcione en absoluto.
Separación o cuasi-separación (también llamada predicción perfecta), una condición en la que el resultado no varía en algunos niveles de las variables independientes.
Tamaño de la muestra: los modelos probit y logit requieren más casos que la regresión OLS porque utilizan técnicas de estimación de máxima verosimilitud. A veces es posible estimar modelos para resultados binarios en conjuntos de datos con solo una pequeña cantidad de casos utilizando regresión logística exacta. También es importante tener en cuenta que cuando el resultado es poco común, incluso si el conjunto de datos general es grande, puede ser difícil estimar un modelo probit.
Pseudo-R-cuadrado: existen muchas medidas diferentes de pseudo-R-cuadrado. Todos intentan proporcionar información similar a la proporcionada por R-cuadrado en la regresión MCO; sin embargo, ninguno de ellos puede interpretarse exactamente como se interpreta el R cuadrado en la regresión MCO.
Diagnósticos: los diagnósticos de la regresión probit son diferentes de los de la regresión OLS. Los diagnósticos de los modelos probit son similares a los de los modelos logit.

Multinomial

Toca el turno de revisar Logit Multinomial

Se utiliza nnet package

Descripcion de los datos

library(foreign)
ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
View(ml)

Estadisticas descriptivas

with(ml, table(ses, prog))

##         prog
## ses      general academic vocation
##   low         16       19       12
##   middle      20       44       31
##   high         9       42        7

with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))

##                 M       SD
## general  51.33333 9.397775
## academic 56.25714 7.943343
## vocation 46.76000 9.318754

pacman::p_load(ggplot2,nnet,foreign,reshape2,xtable,stargazer,tibble,knitr)

ml$prog2<- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)

## # weights:  15 (8 variable)
## initial  value 219.722458 
## iter  10 value 179.982880
## final  value 179.981726 
## converged

#summary(test): salida para cada una de la ecuaciones
stargazer(test,type = "text")  #interpretacion

## 
## ==============================================
##                       Dependent variable:     
##                   ----------------------------
##                      general       vocation   
##                        (1)            (2)     
## ----------------------------------------------
## sesmiddle             -0.533         0.291    
##                      (0.444)        (0.476)   
##                                               
## seshigh              -1.163**       -0.983*   
##                      (0.514)        (0.596)   
##                                               
## write               -0.058***      -0.114***  
##                      (0.021)        (0.022)   
##                                               
## Constant             2.852**       5.218***   
##                      (1.166)        (1.164)   
##                                               
## ----------------------------------------------
## Akaike Inf. Crit.    375.963        375.963   
## ==============================================
## Note:              *p<0.1; **p<0.05; ***p<0.01

El paquete multinom no arroja p-values. Los calculamos usando Wald tests (aca les llaman z-tests)

valor Z:

z <- summary(test)$coefficients/summary(test)$standard.errors
z

##          (Intercept)  sesmiddle   seshigh     write
## general     2.445214 -1.2018081 -2.261334 -2.705562
## vocation    4.484769  0.6116747 -1.649967 -5.112689

# test z a dos colas: p-valor
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p

##           (Intercept) sesmiddle    seshigh           write
## general  0.0144766100 0.2294379 0.02373856 0.0068189015731
## vocation 0.0000072993 0.5407530 0.09894976 0.0000003176045

Extraer los coeficientes del model y exponenciarlos

exp(coef(test))

##          (Intercept) sesmiddle   seshigh     write
## general     17.32582 0.5866769 0.3126026 0.9437172
## vocation   184.61262 1.3382809 0.3743123 0.8926116

Algunos valores predichos

head(pp <- fitted(test))

##    academic   general  vocation
## 1 0.1482764 0.3382454 0.5134781
## 2 0.1202017 0.1806283 0.6991700
## 3 0.4186747 0.2368082 0.3445171
## 4 0.1726885 0.3508384 0.4764731
## 5 0.1001231 0.1689374 0.7309395
## 6 0.3533566 0.2377976 0.4088458

dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
predict(test, newdata = dses, "probs")

##    academic   general  vocation
## 1 0.4396845 0.3581917 0.2021238
## 2 0.4777488 0.2283353 0.2939159
## 3 0.7009007 0.1784939 0.1206054

dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70),
                                                                                   3))

almacena el predicho probabilidades para cada valor de ses and write

pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))

calcular las probabilidades medias dentro de cada nivel de ses

by(pp.write[, 3:5], pp.write$ses, colMeans)

## pp.write$ses: high
##  academic   general  vocation 
## 0.6164315 0.1808037 0.2027648 
## ------------------------------------------------------------ 
## pp.write$ses: low
##  academic   general  vocation 
## 0.3972977 0.3278174 0.2748849 
## ------------------------------------------------------------ 
## pp.write$ses: middle
##  academic   general  vocation 
## 0.4256198 0.2010864 0.3732938

usando melt y ggplot2

lpp  <-  melt (pp.write,  id.vars  =  c ( "ses" ,  "write" ),  value.name  =  "probabilidad" ) 
head (lpp)   # ver las primeras combinaciones

##   ses write variable probabilidad
## 1 low    30 academic   0.09843588
## 2 low    31 academic   0.10716868
## 3 low    32 academic   0.11650390
## 4 low    33 academic   0.12645834
## 5 low    34 academic   0.13704576
## 6 low    35 academic   0.14827643

graficar las probabilidades predichas a trav?s de valores write para cada nivel de ses, plot por tipo de programa

ggplot(lpp, aes(x = write, y = probabilidad, colour = ses)) + geom_line(size=1)+ facet_grid(variable ~., scales = "free")

Fuente: https://stats.idre.ucla.edu/r/dae/multinomial-logistic-regression

Laboratorio 3

Federico Daverio

Lab 3: Probit, Logit

Obtenemos un primer modelo lineal:

Logit

Efectos Marginales

Odds Ratio

Ejemplo 2

Multinomial