Se acercan las elecciones

Se vienen las elecciones en Villa Inferencia y las encuestadoras están dando todo de si para estimar la cantidad de votos que obtendrá cada partido.

• ¿Podrá el intendente ser reelecto? ¿Será vencido por los opositores?

• ¿Hay alguna diferencia en el ingreso promedio de opositores y no-opositores?

Segunda y Tercera encuesta

# Otra encuestadora
muestra_2 <- poblacion %>% 
  sample_n(tamano_muestra)

prop_opo_muestra_2 <- muestra_2 %>% 
  summarize(mean(opositora)) %>%
  as.numeric()

# Otra encuestadora
muestra_3 <- poblacion %>% 
  sample_n(tamano_muestra)

prop_opo_muestra_3 <- muestra_3 %>% 
  summarize(mean(opositora)) %>%
  as.numeric()

# Guardemos resultados en vector y hechemos un vistazo
(muestras_opo <- c(prop_opo_muestra, 
                   prop_opo_muestra_2, 
                   prop_opo_muestra_3))

[1] 0.495 0.493 0.501

Se repite el mismo patrón. La segunda da como ganador al actual intendente, mientras que la tercera anticipa un empáte técnico. Hechemos un vistazo al resto de las encuestas.

Distrubición de la media muestral

Supongamos que disponemos de 500 muestras independientes en total, de 1000 observaciones cada una.

# 500 encuestas...? Aproximando la distribución muestral.
for(i in 1:500){
  muestra <- poblacion %>% 
    sample_n(tamano_muestra)
  
  muestras_opo[i+10] <- muestra %>%
    summarize(mean(opositora)) %>%
    as.numeric()
}

media_dist_muestral <- mean(muestras_opo)
EE_dist_muestral <- sd(muestras_opo)

dist_muestral <- ggplot(data = tibble(x = muestras_opo), 
                        aes(x = x, y = stat(density*width))) +
  geom_histogram() + 
  geom_vline(xintercept = media_dist_muestral)+
  labs(title="Distribución Muestral",
       subtitle=paste0("Promedio de las estimaciones: ",round(media_dist_muestral,4)),
       x="Opositores (%)", 
       y = "Proporción de encuestas")

dist_muestral

Distrubición de la media muestral

• Todos los estimadores tienen su propia distribución.

• La ley de los grandes números postula que el promedio de una muestra lo suficientemente grande de promedios muestrales es igual al promedio poblacional. Es una de las razones por las cuales nos gusta tanto promedio cómo estimador.

• En las próximas secciones discutiremos los estimadores de la regresión lineal.

• Cómo son promedios, aplica la misma lógica.

¿Cual es la relación entre los ingresos y las preferencias políticas?

Si alguno de los comandos de campaña haya alguna relación entre ingreso y preferencias, podría sacar ventaja diseñando campañas electorales acordes. ¹

Comparemos el ingreso promedio de opositores \(E(income|opo=1)\) y no-opositores \(E(income|opo=0)\) en nuestras encuestas.

¿Cómo puedo hacer esa comparación de forma sitemática? Propongo un modelo lineal sencillo y lo estimo con el método de MCO.

La regresión lineal: una máquina de promedios condicionales.

El modelo y su estimación

El modelo sería:

\[ income_i=\alpha + \beta_1 opo_i + u_i \]

Y al estimarlo con MCO nos devolvaría:

• \(\hat\alpha\): \(\hat{E(income|opo=0)}\), el ingreso promedio del grupo de control, osea \(\bar{income_{gob}}\).

• \(\hat\beta_1\): \(\hat{E(income|opo=1)}-\hat{E(income|opo=0)}\), la diferencia entre ingresos promedio del grupo de tratamiento y control, osea \(\bar {income_{opo}}- \bar {income_{gob}}\).

Distribución de estimaciones

A continuación simularemos las estimaciones en 500 encuestas diferentes, y veremos la distribución de las mismas.

set.seed(123)
# Distribución muestral de la correlación entre ingresos y oposición al gobierno actual
muestras_diff_mean <- c()
muestras_B1 <- c()
for(i in 1:500){
  muestra <- poblacion %>% 
    sample_n(tamano_muestra)
  
  # Primero de forma manual (no sistematica)
  
  ## Ingreso promedio de los opositores (promedio de variable ingreso, ponderado por si es opositor o no.)
  mean_opo<-sum(muestra$opositora*muestra$ingresos)/sum(muestra$opositora)
  
  ## Ingreso promedio de pro-intendente (recuerde que el operador lógico "!" es de negación)
  mean_pro<-sum((!muestra$opositora)*muestra$ingresos)/sum(!muestra$opositora)
  
  # Diferencia
  muestras_diff_mean[i] <- mean_opo-mean_pro
  
  # Ahora con MCO
  ## regresion: Ing=B1*op + e
  modelo <- lm(ingresos ~ opositora, muestra)   

  muestras_B1[i] <- modelo$coefficients[2] %>%
    as.numeric()
}

## Calculamos promedios del estimador y el error estandar observado entre una muestra y otra.
media_dist_muestral_diff <- mean(muestras_diff_mean)
EE_dist_muestral_diff <- sd(muestras_diff_mean)

Distribución de estimaciones

#|output-location: column
## Manual
dist_muestral_diff <- ggplot(
  data = tibble(x = muestras_B1), 
  aes(x = x, y = stat(density*width))) +
  geom_histogram() + 
  geom_vline(xintercept = media_dist_muestral_diff)+
  labs(title="Distribución Muestral",
       subtitle = paste0("Promedio de las estimaciones: ",round(media_dist_muestral_diff,2) ),
       x = "Estimaciones de E(Income|anti)-E(Income|pro)", 
       y = "Proporción de encuestas")

dist_muestral_diff

#|output-location: column
## MCO
media_dist_muestral_B1<- mean(muestras_B1)
EE_dist_muestral_B1 <- sd(muestras_B1)

dist_muestral_B1 <- ggplot(
  data = tibble(x = muestras_B1), 
  aes(x = x, y = stat(density*width))) +
  geom_histogram() + 
  geom_vline(xintercept = media_dist_muestral_B1)+
  labs(title="Distribución Muestral",
      subtitle = paste0("Promedio de las estimaciones: ",round(media_dist_muestral_B1,2) ),
       x="Estimaciones de B1", 
       y = "Proporción de encuestas")

dist_muestral_B1

Ventajas de usar la regresión lineal

Las distrbuciones anteriores son identicas.

Lección: Podemos usar la regresión para ver las diferencias promedio entre distintos grupos de individuos.

¿Por qué es mejor usar la regresión?

1. Te permite testear hipótesis acerca de la estimación.

2. Te lo permite hacer de manera mucho más cómoda, pero esto será mas evidente luego (i.e. cuando la var X sea continua, o cuando tengamos que controlar por una tercera variable).

Situación: Tu jefe no te cree nada

Eres un analista de campaña.

Tienes una sóla encuesta a tu disposición para responder la pregunta anterior.

Calculas promedios y encuentras que los opositores son, en promedio, mas ricos por una diferencia:

# Ingreso promediod de opositores
mean_opo<-sum(muestra$opositora*muestra$ingresos)/sum(muestra$opositora)
  
# Ingreso promediod de pro
mean_pro<-sum((!muestra$opositora)*muestra$ingresos)/sum(!muestra$opositora)
  
# Diferencia
mean_opo-mean_pro

[1] 19061.91

Tu supervisor dice que tus resultados son bullshit y que el candidato es transversal (igual de popular en todos los niveles de ingreso). ¿Con qué confianza puedes defender tus resultados?

Usamos las estimaciones de MCO

# regresion: Ing=B1*op + e
modelo <- lm(ingresos ~ opositora, muestra) 

broom::tidy(modelo)

# A tibble: 2 × 5
  term          estimate std.error statistic  p.value
  <chr>            <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    998514.     1348.     741.  0       
2 opositoraTRUE   19062.     1891.      10.1 8.18e-23

Test de hipótesis con ggplot2

Obtenemos una diferencia de ingresos promedio (\(\hat\beta_1\)) alrededor de $16.000, con un error estándar de $1.880.

Recuerda, tu supervisor dice que la verdadera relación es cero. Es decir \(\beta_1=0\).

#|output-location: slide
modelo_tidy<-broom::tidy(modelo)             # Almacenamos coeficientes y SEs en un tibble

std_opositora<-modelo_tidy %>%               # Extraemos SE de B1
  filter(term=="opositoraTRUE") %>% 
  select(std.error) %>% 
  as.numeric()

B1_opositora<-modelo_tidy %>%                # Extraemos el estimado de B1
  filter(term=="opositoraTRUE") %>% 
  select(estimate) %>% 
  as.numeric()

supervisor_ho<-tibble(                       # Tu supervisor dice que B1 poblacional es cero.
  x= 1:5000,
  y=rnorm(5000,mean = 0,sd =std_opositora)   # Simulamos distribucion teorica de B1 estimado alrededor de esa hipotesis
)

p1<-ggplot() +
  geom_histogram(data=supervisor_ho,aes(x = y))+
  geom_vline(xintercept = c(0),
             color="red",
             linetype="dashed")+
  labs(title = paste0("Distribución de estimadores si B1=0 y std.error(B1)=",round(std_opositora)))

p2<-p1+
  geom_vline(xintercept = c(0+std_opositora*1,
                            0-std_opositora*1),
             color="darkslategray1")+
  labs(subtitle="+/- 1 sd cubren 68% del area")

p3<-p1+
  geom_vline(xintercept = c(0+std_opositora*2,
                            0-std_opositora*2),
             color="darkslategray4")+
  labs(subtitle="+/- 2 sd cubren 95.4% del area")

p4<-p1+
  geom_vline(xintercept = c(0+std_opositora*3,
                            0-std_opositora*3),
             color="darkslategray")+
   labs(subtitle="+/- 3 sd cubren 99.7% del area")

p5 <-p4+
  geom_vline(xintercept = B1_opositora,
             color="red",linetype="solid")+
   labs(subtitle="Es improbable que nuestra estimación de B1 haya salido de esta distribución")

Test de hipótesis con ggplot2

p1

p2

p3

p4

p5

Test de hipótesis

Le podés decir a tu jefe: “la probabilidad de que yo haya obtenido estos resultados si tu hipótesis es verdad es menor al 0.01%”.

Eso es lo mismo que checkeamos con los p-valores en el output de regresión.

modelo_tidy %>% 
  filter(term=="opositoraTRUE") %>% 
  select(p.value) %>% 
  as.numeric()

[1] 8.183399e-23

O con el output del comando t.test

t.test(ingresos ~ opositora, muestra) 

    Welch Two Sample t-test

data:  ingresos by opositora
t = -10.076, df = 996.24, p-value < 2.2e-16
alternative hypothesis: true difference in means between group FALSE and group TRUE is not equal to 0
95 percent confidence interval:
 -22774.25 -15349.57
sample estimates:
mean in group FALSE  mean in group TRUE 
           998514.4           1017576.3 

Asi luce el output de la regresión que acabamos de fijar

modelo

Call:
lm(formula = ingresos ~ opositora, data = muestra)

Coefficients:
  (Intercept)  opositoraTRUE  
       998514          19062  

Función summary

Vemos una impresión de los coeficientes que estimamos, pero nada más. Podemos usar la función summary() para ver una impresión mas detallada.

summary(modelo)

Call:
lm(formula = ingresos ~ opositora, data = muestra)

Residuals:
   Min     1Q Median     3Q    Max 
-90184 -19630   -675  20270 100730 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)     998514       1348  740.66   <2e-16 ***
opositoraTRUE    19062       1892   10.08   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 29900 on 998 degrees of freedom
Multiple R-squared:  0.09237,   Adjusted R-squared:  0.09146 
F-statistic: 101.6 on 1 and 998 DF,  p-value: < 2.2e-16

Librería broom

Aveces necesitamos extraer los datos para calculos y presentaciones. La librería broom fácilita el proceso.

broom::tidy and broom::glance

# Tidy retorna los coeficientes, se's y T-statisitcs en un tibble
broom::tidy(modelo)

# A tibble: 2 × 5
  term          estimate std.error statistic  p.value
  <chr>            <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    998514.     1348.     741.  0       
2 opositoraTRUE   19062.     1891.      10.1 8.18e-23

# mientras que glance reporta los estadisticos del modelo en un tibble
broom::glance(modelo)

# A tibble: 1 × 12
  r.squared adj.r.squared  sigma statistic  p.value    df  logLik    AIC    BIC
      <dbl>         <dbl>  <dbl>     <dbl>    <dbl> <dbl>   <dbl>  <dbl>  <dbl>
1    0.0924        0.0915 29903.      102. 8.18e-23     1 -11724. 23453. 23468.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Librería broom

Aveces necesitamos extraer los datos para calculos y presentaciones. La librería broom fácilita el proceso.

broom::augment

# y augment genera una copia del data.frame con los fitted values (predicciones) y los residuos con intervalos de confianza. Tambien sirve para producir predicciones en nueva data.
broom::augment_columns(x = modelo,data = muestra)

# A tibble: 1,000 × 10
    index opositora ingresos  .fitted .se.fit  .resid    .hat .sigma    .cooksd
    <int> <lgl>        <dbl>    <dbl>   <dbl>   <dbl>   <dbl>  <dbl>      <dbl>
 1 658850 FALSE      943408.  998514.   1348. -55107. 0.00203 29867. 0.00347   
 2  33552 TRUE      1054661. 1017576.   1327.  37085. 0.00197 29895. 0.00152   
 3 320338 TRUE       957329. 1017576.   1327. -60247. 0.00197 29857. 0.00401   
 4 251319 TRUE       980136. 1017576.   1327. -37440. 0.00197 29895. 0.00155   
 5 513425 FALSE      996388.  998514.   1348.  -2126. 0.00203 29918. 0.00000516
 6 930865 FALSE      997515.  998514.   1348.   -999. 0.00203 29918. 0.00000114
 7  83651 TRUE       945152. 1017576.   1327. -72424. 0.00197 29830. 0.00580   
 8 203834 TRUE      1038175. 1017576.   1327.  20599. 0.00197 29911. 0.000469  
 9 277737 TRUE       978295. 1017576.   1327. -39282. 0.00197 29892. 0.00171   
10 222122 TRUE      1015827. 1017576.   1327.  -1749. 0.00197 29918. 0.00000338
# ℹ 990 more rows
# ℹ 1 more variable: .std.resid <dbl>

Inspección sanitaria en restaurantes

• Cada país dispone de agencias o departamentos que resguardan la salubridad de los alimentos vendidos al público.

• En el caso de Argentina, este lugar es ocupado por la Administración Nacional de Medicamentos, Alimentos y Tecnología Médica (ANMAT).

• Suponga que es un analista de datos en la ANMAT y le interesa saber si los restaurantes con mayor cantidad de sedes tienen puntajes peores o mejores al resto.

• Cuenta con una base de datos que muestra el puntaje que esta le asignó a una muestra de restaurantes inspeccionados.

res <- causaldata::restaurant_inspections

head(res)

# A tibble: 6 × 5
  business_name            inspection_score  Year NumberofLocations Weekend
  <chr>                               <dbl> <dbl>             <dbl> <lgl>  
1 MCGINLEYS PUB                          94  2017                 9 FALSE  
2 VILLAGE INN #1                         86  2015                66 FALSE  
3 RONNIE SUSHI 2                         80  2016                79 FALSE  
4 FRED MEYER - RETAIL FISH               96  2003                86 FALSE  
5 PHO GRILL                              83  2017                53 FALSE  
6 TACO KING #2                           95  2008                89 FALSE  

Estimemos un modelo lineal con MCO

¿Cual es la diferencia entre el promedio de restaurantes con un local y el de aquellos con dos locales?

¿Cual es la diferencia entre el promedio de aquellos con 100 y aquellos con 101?

¿Si me dan datos nuevos con un restaurante de 250 ubicaciones, que score puedo esperar que tenga?

Podemos usar la siguiente regresión lineal sobre esta base de datos de restaurantes (r) para generalizar este resultado:

\[ Score_r=\beta_0 +\beta_1 NumberofLocations_r + \epsilon_r \]

\(\beta_1\) arroja la diferencia entre el score de inspección promedio de un restaurante y el de cualquier restaurante con un local menos que este.

Estimemos un modelo lineal con MCO

# Perform the first, one-predictor regression
# use the lm() function, with ~ telling us what 
# the dependent variable varies over
m1 <- lm(inspection_score ~ NumberofLocations, data = res)

summary(m1)

Call:
lm(formula = inspection_score ~ NumberofLocations, data = res)

Residuals:
     Min       1Q   Median       3Q      Max 
-27.1673  -3.5449   0.9835   5.4362  17.3253 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       94.8656964  0.0462975 2049.05   <2e-16 ***
NumberofLocations -0.0188715  0.0004356  -43.32   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.052 on 27176 degrees of freedom
Multiple R-squared:  0.0646,    Adjusted R-squared:  0.06456 
F-statistic:  1877 on 1 and 27176 DF,  p-value: < 2.2e-16

Estimemos un modelo lineal con MCO

Sabemos que el equipo de inspección tuvo una rotación alta a lo largo del proceso de recolección de datos y que los procesos de inspección no estaban 100% estandarizados.

¿Cómo podemos hacer si queremos comparar evaluaciones que ocurrieron en el mismo año? Añadimos la variable año como control.

# Now add year as a control
# Just use + to add more terms to the regression
m2 <- lm(inspection_score ~ NumberofLocations + Year, data = res)

summary(m2)

Call:
lm(formula = inspection_score ~ NumberofLocations + Year, data = res)

Residuals:
     Min       1Q   Median       3Q      Max 
-27.4407  -3.5210   0.9142   5.3014  18.0722 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        2.253e+02  1.241e+01   18.16   <2e-16 ***
NumberofLocations -1.919e-02  4.358e-04  -44.03   <2e-16 ***
Year              -6.489e-02  6.173e-03  -10.51   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.039 on 27175 degrees of freedom
Multiple R-squared:  0.06838,   Adjusted R-squared:  0.06832 
F-statistic: 997.4 on 2 and 27175 DF,  p-value: < 2.2e-16

Librería modelsummary

La libreria modelsummary dispone de alternativas para generar estos resumenes de forma práctica y flexible.¹

# Default significance stars are: +/*/**/*** .1/.05/.01/.001. 
# Social science standard */**/*** .1/.05/.01 can be restored with
modelsummary::msummary(list(m1, m2),
    stars=c('*' = .1, '**' = .05, '***' = .01),
    gof_omit = "AIC|BIC|Log")

Predicciones del modelo

Podemos presentar los promedios condicionales provenientes de ambas regresiones (fitted values) en nuestro gráfico de la siguiente manera:

Modelo 1

predicciones<-broom::augment_columns(x=m1,
                                     data=res)

# A tibble: 6 × 12
  business_name inspection_score  Year NumberofLocations Weekend .fitted .se.fit
  <chr>                    <dbl> <dbl>             <dbl> <lgl>     <dbl>   <dbl>
1 MCGINLEYS PUB               94  2017                 9 FALSE      94.7  0.0440
2 VILLAGE INN …               86  2015                66 FALSE      93.6  0.0367
3 RONNIE SUSHI…               80  2016                79 FALSE      93.4  0.0372
4 FRED MEYER -…               96  2003                86 FALSE      93.2  0.0379
5 PHO GRILL                   83  2017                53 FALSE      93.9  0.0371
6 TACO KING #2                95  2008                89 FALSE      93.2  0.0382
# ℹ 5 more variables: .resid <dbl>, .hat <dbl>, .sigma <dbl>, .cooksd <dbl>,
#   .std.resid <dbl>

Predicciones del modelo

Podemos presentar los promedios condicionales provenientes de ambas regresiones (fitted values) en nuestro gráfico de la siguiente manera:

Modelo 1

p3<-p1 + 
    geom_line(data=predicciones  ,
               aes(x=NumberofLocations,y=.fitted,
                   color="B1*NumberofLocations\n(M1)"))+
    labs(title="Fitted values del modelo 1",
         subtitle=NULL)

Predicciones del modelo

Podemos presentar los promedios condicionales provenientes de ambas regresiones (fitted values) en nuestro gráfico de la siguiente manera:

Modelo 2

predicciones_05_25<-broom::augment_columns(x=m2,
                                           data=res) %>% 
  filter(Year %in% seq(2005,2020,by=5))

p4<-p1 + 
    geom_line(data = predicciones_05_25,
               aes(x=NumberofLocations,y=.fitted,group=Year,
                   color=paste("B1*NumberofLocations\n + B2*",Year)))+
    labs(title="Fitted values del modelo 2",
         subtitle="Para 2005 y 2015",
         color="Fitted values")

Note que los fitted values del modelo 2 los score promedio del año 2005 son mayores a los de 2015, pero la relación lineal es la misma.

Ejercicio

1. Encontramos que el coeficiente del Número de Sedes sobre el Score de Inspección era negativo y significativo. En sus palabras, ¿Que significa este resultado? ¿Nos permite afirmar que los restaurantes relajan sus estándares de calidad e higiene cuando se expanden?

2. Tratemos de replicar los resultados de Lalonde (1986).[Lalonde, Robert. 1986. “Evaluating the Econometric Evaluations of Training Programs with Experimental Data.” American Economic Review 76 (4): 604–20 ]

El estudio busca identificar el efecto de un programa de entrenamiento laboral sobre los ingresos de los asistentes tras dos años de su culminación, en el período 1975-1978.

La primera parte estudia el efecto sobre personas que fueron asignadas a programas de entrenamiento experimentalmente. Esos datos se encuentran en el data.frame jtrain2.

Luego se analizan datos donde las personas podian escoger si asistir o no de forma voluntaria. Estos datos se encuentran en el data.frame jtrain3.

Ejercicio

El ejercicio consite en:

• Reproduzca ambos resultados usando un modelo de regresión lineal.

• Presente ambos en la misma tabla de regresión e interprete los coeficientes. Si las hubieren, discuta las diferencias entre los coefientes de cada regresión.