ccs-ayudantia2-tutorial.qmd

---
title: "Datos observacionales y estrategias de identificación"
subtitle: "Ayudantía 2- Curso métodos en CCS"
#author:
#  name: Melanie Oyarzún
#  affiliation: Universidad del Desarrollo
#  email: moyarzunw@udd.cl
#date: "`r format(Sys.time(), '%d %B %Y')`"

format: 
    html:
        theme: cosmo
        toc: true
        toc-depth: 3
        embed-resources: true
        self-contained-math: true
        
keep-md: true
editor: source

execute:
  eval: true
  warning: false

number-sections: false
crossref:
  chapters: true
---

# I. Bienvenida al tutorial

Hola, este tutorial está escrito como la actividad de aplicación de la ayudantía 2 del curso

La idea es ir desarrollando el tutorial en paralelo a la sesión y que puedan tomar apuntes adicionales directamente en el notebook. Está escrito como .qmd, un notebook Quarto y puede ser editado en rstudio, vscode o cualquier otro IDE o editor de texto. El lenguaje utilizado para el procesamiento y análisis es R.

Se recomienda crear un fork al repositorio y trabajar en su propia rama.

**Como funciona**

-   El tutorial se estructura como un proyecto de investigación ficticio completo, que empezaremos cargando datos brutos hasta obtener e interpretar coficientes de regresión.

-   Lo métodos a cubir son:

    1. Regresión lineal y método de diferencias (en datos experimentales)

    2. Logit y probit

    3. Variables Instrumentales

    4. Regresión Discontinua

    5. Diferencias en Diferencias

    6. Efectos fijos, aleatorios y modelos jerárquicos

    7. Matching

-   Los datos y notebook son descargables en el siguiente [repositorio de proyecto](https://github.com/ClasesMOW/ayudantia2ccs)  así que se recomienda correr todo en tu propio computador o en una instancia en la nube.

-   Este tutorial es una adaptación del trabajo de [Hans H. Sievertsen](https://hhsievertsen.shinyapps.io/applied_econ_with_R_dynamic), agregando la estimación con métodos de de métodos de matching y algunos cambios a la situación a analizar.

<br>

## 1 Pregunta de investigación y datos

Vamos a usar una situación ficticia, con datos simulados para poder aplicar de manera consisa los diferentes tipos de estrategias de identificación revisadas en clase. Algunas (que se supone ya manejan) las revisaremos muy rápidamente y en otras, vamos a tener mayor énfasis.

### 1.1 Pregunta de investigación

Nuestro objetivo es responder la siguiente pregunta **ficticia** de investigación:

> Asistir a cursos de verano mejora los resultados académicos?

Para responder esta pregunta, usaremos unos datos **ficticios y simulados**

### 1.2 Contexto ficticio

La pregunta de investigación se inspira en trabajos como el de [Matsudaira (2007)](https://www.sciencedirect.com/science/article/pii/S0304407607001194?casa_token=hnnF764CKPoAAAAA:5b9WhCManNDsdW4SmOHnnzNr0fZIarW8s6EsvpQW7MdUt470eNPmN2T8IFCsNc6Iajew5tEeNA) e intervenciones en estudiantes de bajo nivel socioeconómico por [Dietrichson et al ( 2017)](https://journals.sagepub.com/doi/abs/10.3102/0034654316687036).

El **escenario ficticio** es el siguiente:

-   Para un conjunto de colegios en una comuna, existe la opción de asistir a un curso de verano intensivo durante el verano entre 5 y 6to básico.
-   El curso de verano se enfoca en mejorar las habilidades académicas de preparar la prueba de admisión a la universidad vigente (PSU en ese momento)
-   El curso de verano es gratuito, pero para ser matriculados requiere que los padres se involucren en un proceso.
-   Estamos interesados en testear el impacto de la participación en el curso en los resultados académicos de los estudiantes.

### 1.3 Datos ficticios dispobibles

1.  school_data_1.csv

-   Usamos esta data para ejemplificar como cargar data guardada en formato csv.
-   Este dataset tiene información sobre cada individuo (con identificador id), la escuela a la que asiste, un indicador si participó en el curso de verano, sexo, ingreso del hogar (en logaritmo), educación de los padres, resultados en una prueba estandarizada que se realiza a nivel de la comuna tanto para el año 5 como para el año 6.

2.  school_data_2.dta

-   Usamos esta data para ejemplificar como cargar data guardada en formato STATA.
-   Este dataset tiene información de cada individuo (con identificador id).
-   Este dataset tiene la información si el individuo recibió la carta de invitación para participar del curso de verano.

3.  school_data_3.xlsx

-   Usamos este dataset para practicar como cargar datos guardados en formato Microsoft Excel.
-   Este dataset incluye datos sobre cada individuo (con identificador id)
-   Este dataset tiene información de rendimiento académico antes y después del curso de verano.

# II. Preparación de los datos

En esta sección vamos a preparar los datos para el análisis. Este proceso generalmente incluye cargarlos, inspeccionarlos, limpiar y dar la estructura deseada. También revisaremos estadísticas descriptivas que nos den una idea antes de estimar cualquier modelo.

En la vida real, este proceso suele ser bastante largo, laborioso e implica volver sobre lso pasos anteriores múltiples veces. También invlocura tomar desiciones por parte de los investigadores, por lo cual la documentación de esta fase es especialmente importante.

En nuestro caso, será bastante lineal y directo, ya que son datos ficticios simulados. Pero en la realidad, no suele ser así.

## 0. Crear el proyecto y clonar el repositorio

Este proyecto está sustentado en el repositorio de github. Clónelo y luego carguelo en su computador.

## 1 Cargar los datos

### 1.1 Intalar y cargar paquetes

Para poder cargar los datos, necesiamos los paquetes adecuados. Siempre hay multiples formas de hacer las cosas en cualquier lenguaje o programa. En este caso, usaremos la función `read_csv()` del paquete *readr*. Para poder usarlo, debemos estar seguros de que está instalado.

Si no está instalado, podemos hacerlo con la función `install.packages("[nombre paquete a instalar]")`

```{r}
#| eval: false

install.packages("readr")
```

Si el paquete ya está instalado, para poder usarlo necesitamos tenerlo cargado en nuestra librería. Para esto usamos la función `library("[nombre paquete a cargar]")`.

```{r}
library(readr)
```

Cada paquete lo debemos instalar solo una vez por computador, pero debemos cargarlo en cada sesión para poder utilizarlo.

### 1.2 Cargar datos csv

Con *readr* estamos en condiciones de usar la función `read_csv()` para cargar la primera base de datos.

Vamos a cargar *school_data_1.csv* agregando el path a los datos en paréntesis. Puede reemplazar por el path correspondiente o usar el paquete auxiliar `here`

PS. notemos que tambien estamos incluyendo comentarios en los bloques de código. Las líneas que empiezan con el símbolo `#` son ignoradas por R.

```{r}

# cargar readr package
library("readr", "here")

#prueba el directorio que te dice con here
here::here()
```

Entonces, dentro de `here("[Escribes el directorio relativo]")` actuará como el directorio relativo sin errores. Entonces, cargamos los datos y los asignamos a un data frame.

```{r}

school_data_1 <- read.csv(here::here("data_raw/school_data_1.csv"))

```

Usualmente, después de cargar un dataset es útil visualizarlo. Empleamos la función `head()` para ver sus primeras 6 observaciones.

```{r}
#| warning: false
#| eval: true

head(school_data_1)

```

Una segunda alternativa es descargar y cargar los datos directamente en R desde internet. Puede ser cualquier link directo o, si está alojado en github, tienes que asegurarte de que sea un repositorio público.

```{r}
#| eval: false

#school_data_1 <- read.csv("https://github.com/ClasesMOW/ayudantiasccs/blob/main/data_raw/school_data_1.csv")
```

### 1.3 Cargar datos STATA

En Ciencias Sociales y Economía es muy comun contar con datos para ser utilizados en el programa STATA, estos son archivos que terminan en *.dta*. Para cargarlos, vamos a usar *HAVEN* del paquete *tidyverse*. También usaremos muchas otras funciones de ese paquete, asi que es buen momento para cargalo.

```{r}
library(tidyverse)

school_data_2<- haven::read_dta(here::here("data_raw/school_data_2.dta"))
```

Usemos el comando `tail()` para ver las últimas 10 entradas.

```{r}
# print las últimas 8 filas
tail(school_data_2,n=10)
```

### 1.4 Cargar datos Microsoft Excel

Finalmente, cargaremos el tercer dataset guardado como hoja de cálculo de Excel *.xlsx*. Para esto usaremos el paquete *readxl* que viene incluido en *tidyverse*. Luego de asignarlo a un dataframe, démosle una mirada con `glimpse()` (también del tidyverse).

```{r}
school_data_3 <- readxl::read_xlsx(here::here("data_raw/school_data_3.xlsx"))

glimpse(school_data_3)
```

## 2. Unir los datasets

Tenemos 3 bases de datos con información diferente de los mismos individuos. Generalmente es buena idea tener una sola gran tabla con toda esta información, especialmente si estimaremos modelos en base a ésta.

La base de datos 1 y 2 tienen una forma similar: los individuos son filas y las variables columnas y hay una sola fila para cada individuo.

Para hacerlo, podriamos usar varias alternativas.

Una alternativa es usar la función nativa `merge( )`. En esta función, primero mencionamos los datasets a unir, luego informamos cual es la(s) columnas(s) que debe usar para unir ambos datasets con `by=....`. Por defecto, R incluye todas las filas que están en ambos datasets (basados en la variable *by*), pero podemos fijar `all=TRUE` para mantener todas las filas que están en ambos datasets o `all.x=TRUE` para mantener todas las filas coincidentes y las del primer dataset or `all.y=TRUE` para guardar todas las filas del segundo dataset.

Veamos un ejemplo con los dos primeros datasets. Luego usemos `dim()` para conocer las dimensiones del nuevo dataset unido en términos de filas y columnas.

```{r }
# Unir school_data_1 con school_data_2
school_data_merged <-merge(school_data_1,school_data_2, by="person_id")

# Revisamos las dimensiones
dim(school_data_merged)
```

Notemos que el dataset unido tiene 3491 filas y 9 columnas. Unimos todas las filas y agregamos al dataset de información del estudiante si recibió o no la carta (school_data_2)

¿Qué ocurre si las columnas tienen igual nombre? R va a renombrarlas automáticamente agregando un sufijo *.x* (a la columna del primer dataset) e *.y* ( a la columna del segundo dataset).

Entonces, en el siguiente bloque:

1.  Unimos *school_data_1* y *school_data_2* usando como variable de unión *person_id* y guardamos el dataset unido como *school_data*.
2.  Unimos *school_data_3* con *school_data* y sobre-escribimos *school_data*.

Notar que acá unimos por las columnas *person_id* y *school_id*. Esto no es realmente necesario porque cada estudiante con id única tiene un solo colegio, pero sirve de ejemplo en como usar más de una columna mediante `c()`. 3. Usamos la función `summary()` para obtener una estadística descriptiva de las variables en el dataset unido.

```{r }
# Merge school_data_1 y school_data_2 y guardamos como school_data_merged 
school_data_merged<-merge(school_data_1,school_data_2,by="person_id")

# Merge school_data_3 con school_data_merged
school_data_merged<-merge(school_data_merged,school_data_3,by=c("person_id","school_id"))

# Estadísticas descriptivas de cada variable
summary(school_data_merged)
```

Otra opción de hacer lo mismo es con join del paquete `dplyr` del `tidyverse`.

[![Diferentes tipos de join en dplyr](images/jointypes.png){fig-alt="Diferentes tipos de join en dplyr" fig-align="center"}](https://statisticsglobe.com/r-dplyr-join-inner-left-right-full-semi-anti) (foto sacada de [link](https://statisticsglobe.com/r-dplyr-join-inner-left-right-full-semi-anti), puedes ver ejemplos detallados de los diferentes tipos de join)

En este caso, podríamos hacer en un pipe `%>%` dos left joins seguidos.

```{r}
school_data_merged <- school_data_1 %>%
                        left_join(school_data_2, by="person_id") %>%
                        left_join(school_data_3, by=c("person_id", "school_id"))

```

## 3. Lipiar los datos

### 3.1 Tidyng los datos

Ahora que hemos unido las bases de datos, trataremos de que satisfazgan los principios de [Tidy Data](https://vita.had.co.nz/papers/tidy-data.pdf).

Un data frame se considera "tidy" (Según Hadley) si se cumplen las siguientes condiciones:

-   Cada columna (variable) contiene todas las medidas de la misma data/feature/caracteristica de las observaciones.
-   Cada fila contiene medidas de la misma unidad de observación.

(puedes profundizar y ver más ejemplos aplicados en <https://sscc.wisc.edu/sscc/pubs/DWE/book/6-2-tidy-data.html> )

Uno de estos, es que cada columna debe ser una variable y cafa fila una unidad de observación.

Si inspeccionamos el número de columnas:

```{r}
# Alternativamente podemos usar la función nrow() para obtener el número de filas

ncol(school_data_merged)

```

Son 17, pero tenemos 9 variables.

```{r}
head(school_data_merged)
```

Es porque tenemos puntajes de las pruebas del año 2 al 10. Este tipo de datos son de *panel*

![](images/paste-4C87CB1E.png)

Generalmente, que hagamos con este tipo de datos depende del tipo de modelos que queramos usar. Si bien el formato wide es facil de entender, generlamente para modelos y análisi preferimos que esté en formato long. Especialmente cuando modelamos incluyendo efectos fijos También es este el que adhiere a los principios tidy de mejor manera.

Para cambiar a long, usamos `pivot_longer()`

```{r}
# lCargamos el paquete tidyr package. Ya está incluido en tidyverse, pero  también se puede llamar por si solo.

library("tidyr")

# make data tidy (make long)
school_data_tidy<-school_data_merged%>%
       tidyr::pivot_longer(
                            cols = starts_with("test_year"),
                            names_to = "year",
                            names_prefix = "test_year_",
                            names_transform = list(year = as.integer),
                            values_to = "test_score",
                             )

# ncol nos da el número de columnas del nuevo dataset
ncol(school_data_tidy)
```

Ahora tenemos nuestros datos listos para que los inspeccionemos.

### 3.2 Selección de muestra

Ya que contamos con datos que siguen los principios de tidy data, lo siguiente es seleccionar la muestra apropiada. En este trabajo, los unicos problemas que podríamos enfrentar son relacionados con valores faltantes o missing. Para inspeccionarlos vamos a usar la función `skim()` del paquete `skimr`, esta función nos muestra los vaores faltantes en nuestro dataset de una manera global.

```{r}
# Cargar skimr
library("skimr")

# Usamos skim() para inspeccionar los datos
skim(school_data_tidy)
```

Con esta función podemos ver facilmente cuantas filas y comunas son, los tipos de varables y número de missing values. Además la media, desviación estándar, percentiles e incluso un histograma para cada variable.

En estos casos, para *parental_schooling* tenemos 45 missing y para *test_score* 11. Asumamos que estos valores missing son random y deseamos remover estas filas. Para esto usamos `filter()`. Esta funcion toma dos argumentos, el dataset a filtrar y la condición para que se mantenga en el dataset, en este caso que no sea na o `!is.na(partental_schoolin)`. La función `is.na()` es verdad cuando el elemento en `()` es missing y usamos `!`para mostrar que queremos lo contrario a esto condición. En otras palabras queremos que la educaión parental no esté missing.

```{r}
# Seleccionamos las columnas sin missing values
school_data_selected<-dplyr::filter(school_data_tidy,!is.na(parental_schooling),!is.na(test_score))

# Usamos skim() para revisar los datos nuevamente
skim(school_data_selected)

```

Hemos removido todos los missing.

### 3.3 Modificar los datos

Un último paso que haremos antes de hacer estadística decsriptiva es modificar los nombres de algunas columnas para que se vean bien en la tablas.

Vambos renombrar la variable summpercap a summerschoolo. Lo hacemos con `rename()` del paquete dplyr. Esta función tiene una sintaxix similar a `filter()`

```{r}
# renombremos summercamp a summerschool
analysisdata<-rename(school_data_selected, summerschool=summercamp)

# usamos head para visualizar las primeras 6 observaciones
head(analysisdata)
```

En un siguiente paso, vamos a transformar los puntajes en la pruebas a una variable que tenga media 0 y desviación estándar 1. Es mejor trabajar con variables estandarizadas, ya que no requieren conocer el conexto específico de la medida y es más facil de comunicar y comparar.

```{r}
# Estandarizamos los resultados de las pruebas
# Agrupamos analysisdata por year
analysisdata<-group_by(analysisdata,year)

# Creamos una nueva variable con mutate. Como queremos que reemplace a la anterior, usamos su mismo nombre. 
analysisdata<-mutate(analysisdata, test_score=(test_score-mean(test_score))/sd(test_score))

# mostremos la media
print(paste("Mean of test score:",mean(analysisdata$test_score)))
```

Esto tambien podemos hacerlo dentro de un pipe:

```{r}
# Estandarizamos los resultados de las pruebas

analysisdata<- analysisdata %>% 
        group_by(year) %>% # Agrupamos analysisdata por year
        mutate(test_score=(test_score-mean(test_score))/sd(test_score))

# mostremos la media
print(paste("Mean of test score:",mean(analysisdata$test_score)))
```

Podemos comprobar que efectivamente la media y desviación estándar corresponden a dichos valores:

```{r}
# ver la media de test_score
print(paste("Mean of test score:",mean(analysisdata$test_score)))

# Ver la desviación estándar de test_score
print(paste("SD of test score:",sd(analysisdata$test_score)))
```

Notar que si aplicamos `mean()` o `sd()` o cualquier otra función matemática a columnas que tienen valores missing, también dará un valor missing a menos que usemos la opción na.rm=FALSE.

Ya estamos bien, ahora pasamos a conocer mejor nuestros datos con estadística descriptiva.

## 4. Estadística descriptiva

Hasta ahora, cargamos datos en diversos formatos (csv, dta y xlsx) los unimos, re-estructuramos el dataset, removimos valores missing y generamos algunas transformaciones. El siguiente paso es empezar a conocer nuestros datos. Para esto haremos tablas de estadísticas descriptivas y también algunos graficos descriptivos.

### 4.1 Tablas de estadística descriptiva

Hasta ahora, ya conocemos dos maneras de calcular estadísticas resumen:

1.  `sumary()` de R base. Esta función en realidad funciona en muchos tipos de objetos de R y suele dar un bien resumen. Pero no en el formato de una tabla exportable a un documento latex, word o etc. que podamos presentar en nuestra investigación o resultados.

2.  `skim()` del paquete *skimr*

#### 4.1.1 Tabla de estadísticas descriptivas "lista para llevar"

Una forma rápida de obtener una tabla de estadísticas descriptivas es con un primo de `skim()` del paquete *modelsummary*.

```{r}
# cargar modelsummary
library("modelsummary")

# creamos una tabla de estádisticas descriptivas
analysisdata%>%
  filter(year==2)%>%
  select(female,starts_with("paren"),letter,summerschool,test_score)%>%
  datasummary_skim( fmt="%.2f")

```

Esta se puede exportar a varios formatos, como word o latex con el parámetro `output=["ruta donde guardar la tabla"]`. Primero hagámoslo en word:

```{r}
# load modelsummary
library("modelsummary")

# create a summary stat table in Latex format
analysisdata%>%
  filter(year==2)%>%
  select(female,starts_with("paren"),letter,summerschool,test_score)%>%
  datasummary_skim( fmt="%.2f",
                 histogram=FALSE, output="tab_summary_statistics.docx")

```

La guardó en la carpeta por default. Si quieremos que esté en nuestra carpeta de output, podemos usar el paquete here::

```{r}
# load modelsummary
library("modelsummary")

# create a summary stat table in Latex format
analysisdata%>%
  filter(year==2)%>%
  select(female,starts_with("paren"),letter,summerschool,test_score)%>%
  datasummary_skim( fmt="%.2f",
                 histogram=FALSE, output=here::here("output/tab_summary_statistics.docx") )

```

Ahora la hacemos en formato latex:

```{r}

# load modelsummary
#library("modelsummary")

# create a summary stat table in Latex format
#analysisdata%>%
#  filter(year==2)%>%
#  ungroup() %>%
#  select(female,starts_with("paren"),letter,summerschool,test_score)%>%
#  datasummary_skim( fmt="%.2f", 
#                    histogram=FALSE, output=here::here("output/tab_summary_statistics.tex"))
    
```

#### 4.1.2 Tablas customizadas

Para customizar nuestra tabla aun más, podemos usar la función `datasummary()` tambien del pquete modelsummary. Esta función perimte que definamos una *fórmula* de la estructura de la tabla.

```{r}
# creamos una tabla de estadísticas descriptivas resumen

datasummary(female+parental_schooling+
              letter+test_score~Factor(summerschool)*(Mean+SD),
            sparse_header = FALSE,
            data=filter(analysisdata,year==2))

```

En este ejemplo: - Listamos las variables a incluir separadas por un + similar to in a `female+parental_schooling+pa...` - Usamos a `~` para separar la lista de variables en la fórmula - Usamos la formula `Factor(summerschool)*(Mean+SD)` para mostror la media y desviación estándar por separado para cada grupo creado por la variable *summerschool*. - Usamos `Factor()` para indicarle a R que debería considerar *summerschool* como una variable binaria.

- También podríamos haber hecho esto al limpiar y procesar la base de datos. - Podemos tambien invertir el orden `(Mean+SD)*Factor(summerschool)`, lo que entonces daría primero la media y devsicacion estándar y luego separar por los valores de la escuela de verano.

-   Usamos la opción `sparce_header=FALSE` para especificar que queremos incluir e *summerschool* como título.

#### 4.1.3 Nombres de variables

Hasta ahora hemos utilizado los nombres de variables directamente en las tablas. Estó no es muy estético, podemos cambiarle el nombre directamente con espacios y mayúsculas en el nombre para darle un emjor aspecto. También es posible asignanrle una "label" o etiqueta cuando creamos la tabla, como lo vemos en el ejemplo:

```{r}
# load modelsummary

library("modelsummary")
# create a summary stat table

datasummary((`Female`=female)+
            (`Parental schooling (years)`=parental_schooling)+
            (`Parental income (log)`=parental_lincome)+
            (`Received reminder letter`=letter)+
            (`Test Score`=test_score)~
            (`Attended summer school`=Factor(summerschool))*
              (Mean+SD),
            sparse_header = FALSE,
            data=filter(analysisdata,year==2))

```

#### 4.1.4 Exportando nuestras tablas

Podemos exportar nuestras tablas a word o Latex podemos usar la expresión `output="[nombre del archivo y ruta]"`. En `datasummary()` se ve así:

```{r}
# load modelsummary
library("modelsummary")
# create a summary stat table

datasummary((`Female`=female)+
            (`Parental schooling (years)`=parental_schooling)+
            (`Parental income (log)`=parental_lincome)+
            (`Received reminder letter`=letter)+
            (`Test Score`=test_score)~
            (Mean+SD+P25+P50+P75),
            sparse_header = FALSE,
            data=filter(analysisdata,year==2),
            output = here::here("output/tab_descriptive_statistics.docx"))
```

### 4.2 Gráficos de estadística descriptiva

Vamos a usar principalmente la librería ggplot2 para crear nuestros gráficos.

#### 4.2.1 Scatter plot (o gráfico de dispersión)

Nuestro primer grafico es un gráfico de dispersión. En este queremos ver como dos variables se relacionan en los datos. En estos podemos inlcuir curvas que describan la relación.

En este gráfico

1.  Iniciamos un objeto de grafico `ggplot()` usando los datos *analysisdata* que ya procesamos y la vamos a filtrar solo para incluir el año 5.
2.  Especificamos que *parental_income* sea el eje x y *test_score* el eje y en `aes()`
3.  Usamos `geom_point()` para incluir los puntos que describen la dispersión.
4.  Usamos `geom_smooth()` para agregar una linea que describa la relación.
5.  Usamos `theme()` para darle formato a los elementos
6.  Usamos `labs()` para incorporar etiquetas a los ejes y al título.

```{r}
# load ggplot2
library("ggplot2")
# creamos un scatter plot entre ingreso parental y resultados academicos en el año 5
ggplot(analysisdata%>%filter(year==5),  
       aes(x=parental_lincome,y=test_score))+
       geom_point(alpha=0.1,size=0.85,color="#63a668")+
       geom_smooth(color="#145c21") +
       theme(panel.background = element_rect(fill="#ededed",color="#ededed"),
             plot.background = element_rect(fill="#ededed",color="#ededed"),
             panel.grid.major = element_line(colour="#a3a3a3",size=0.1))+
  labs(x="Log(Parental Income)",y="Test Score (Mean=0,SD=1)", title="Test scores & Parental income")
```

#### 4.2.2 Graficos de barras y boxplot

Del grafico anteriro podemos observar que los resultados de los test se correlacionan con ingreso parental. Esto no es una sorpresa. Veamos si tambien el asistir a la escuea de verano se correlaciona con estas características individuales.

Primero, creamos un scatter plot de educación de los padres y test score en el año anterior a la escuela de verano. Usemos el mismo código de arriba pero en lugar de ingreso, usamos educación parental. Segundo, creemos un gráfico de barras que muestre que los resultados previos a la asistencia a la escuela de verano. Tercero, creamos un box plot de ingreso de los padres y si asistieron o no a la escula de verano.

Nuestro codigo ahora tiene estos elementos.

-   Creamos el objeto `ggplot()` y cargamos la data y el tema. Este objeto es llamado un *rawchart*
-   Creamos 3 gráficos basados en *rawchart*, cada uno lo guardamos con un nombre.
    -   Para el de barras, usamos `geom_bar()` y fijamos `stat="summary", fun="mean"` para decirle a R que cree un grafico de barras con la media de *test_score*\
    -   usamos `labs()` para decir los ejes y titulos.
    -   Usamos `geom_boxplot()` para crear el boxplot.
-   Usamos el paquete *patchwork* para combinar varios gráficos en uno.
-   usamos *ggsave()* para guardar el gráfico combinado como un archivo *png*

```{r}
# Load patchwork 
library("patchwork")
# Create raw chart element
rawchart<-ggplot(analysisdata%>%filter(year==4),x=as.factor(fill))+
          theme_classic()

# Create bar chart of pre summer school test score and summer school 
p1<-rawchart+
       geom_smooth(aes(x=parental_schooling,y=test_score)) +
       geom_point(aes(x=parental_schooling,y=test_score),alpha=0.1)+
       labs(x="Parental schooling", y="Test Score Year 5")

# Create bar chart of pre summer school test score and summer school 
p2<-rawchart+
       geom_bar(aes(x=as.factor(summerschool),y=test_score),
                    stat="summary",fun="mean")+
       labs(y="Test Score Year 5", x="Attended Summer School")

# Create bar chart of parental schooling and summer school attendance
p3<-rawchart+
              geom_boxplot(aes(x=as.factor(summerschool),y=parental_lincome))+
       labs(y="Parental Income (log)", x="Attended Summer School")

# Combine charts
p1/(p2+p3)

# Export chart
ggsave("fig1.png")
```

Estos tres gráficos nos muestran que los puntajes en las pruebas se correlacionan con características de los padres (ambos scatters) y que aquellos que asistieron a la escuela de verano tenían mejores puntajes incluso antes de ir a la escuela de verano y que las caraterísticas de los padres se relaciona con la asistencia al éste 

Es decir, tenemos **SESGO DE SELECCIÓN**

#### 4.1.3 Histogramas y gráficos de densidad

Comaperemos la distribución de los puntajes depués de asistir a la escuela de verano por quienes fueron y los que no fueron. Creamos un histograma con una linea que muestre las densidades estimadas.

```{r}
# create a histogram and density chart
ggplot(filter(analysisdata,year==6),
       aes(x=test_score,fill=as.factor(summerschool)))+
       geom_histogram(aes(y=..density..),bins = 50,alpha=0.5,
                       position="identity",color="white")+
       geom_density(alpha=0.0,size=1,show.legend= FALSE)+
       theme_minimal()+
       labs(y="Densidad",x="Puntaje en prueba año 6 (estandarizado)",fill=" ")+
       scale_fill_brewer(palette="Set2",labels=c("No asistió","Asistío a la escuela de verano"))+
       theme(legend.position="top")
```

Claramente hay diferencias en sus resultados. Pero la estimación directa va a confundir cuanto de esto proviene del sesgo de selección y cuanto es el efecto real de ir a la escuela de verano.

#### 4.1.4 Correalograma

Muchas veces queremos saber que tán correlacionadas estan las variables en una muestra. Podríamos simplementa calular una tabla de correlaciones, pero el paquete \[PerformanceAnalytics\] tiene una función muy conveniente: `chart.Correlation()` que nos presenta un gráfico con las correlaciones de a pares, su significancia estadística, gráficos de dispersión y distribución.

```{r}

#Cargamos el paquete
library("PerformanceAnalytics")

correlations_graph <-  analysisdata %>%
            dplyr::filter(year==6) %>%
            ungroup() %>%
            dplyr::select(female, parental_schooling, parental_lincome, summerschool, letter )  #variables para el correalograma

chart.Correlation(correlations_graph, histogram=TRUE, pch=19)


```

# III. Estimación e identifiación

Tenemos nuestro problema base claro, ahora usemos algunas estrategias de identificaión para encontrar el efecto de asitir a una escuela de verano en esta muestra.

Para esto, partiremos del mejor caso: contamos con un experimento y luego iremos variando la información disponible y la estrategia que es factible emplear.

## 1. Regresión lineal ingenua

Nuestro objetivo es estimar la relación en los resultados académicos y asistira la escuela de verano. Para esto podríamos pensar en un modelo que directamente incluya a estos dos elementos.

$$ Resultado_i = \beta_0 + \beta_1 I(Summerschool=1)_i + u_i $$ 
El modelo estimado sería:

```{r}
ols <- lm ( test_score ~ summerschool, data= analysisdata )
ols
```

```{r}
summary(ols)
```

Podemos hacer una mejor tabla usando el paquete *stargazer*

```{r}
library(stargazer)

stargazer::stargazer(ols, type = "text", out="tabla1.tex")

```

Stargazer es util, pero esta siendo reemplazado por otros paquetes, particularmente [modelsummary](https://vincentarelbundock.github.io/modelsummary/index.html). Para una refrencia extra recomiendo revisar este [link](https://elbersb.de/public/pdf/web-7-regression-tables-graphs.pdf)

```{r}

# Podemos hacer un mapeo de nombres de variables


cm <- c( "test_score" = "Test Score",
         "summerschool" = "Escuela de verano" )
```

```{r}
library(modelsummary)

modelsummary::modelsummary(ols, stars= TRUE, fmt=3, 
                           estimate = "{estimate} ({std.error}){stars}", statistic=NULL,  gof_omit = "AIC|BIC|Lik",
                             coef_map = cm)

```

También permite graficar:

```{r}
modelsummary::modelplot(ols, coef_omit = "intercept")
```

No se ve muy bien para solo 1 modelo. Otra opcion es con el paquete parameters (<https://easystats.github.io/parameters/>)

```{r}
library(parameters)

m <- parameters::model_parameters(ols)
m
  
```

Este modelo, como ya lo revisamos de la estadística descriptiva, podemos pensar que sufre de sesgo de selección. Es posible que los estudiantes con padres más involucrados, con más educación y mayor ingreso

Si solo fueran el problema la eduación de los padres y su ingreso, podríamos controlar por estos. Incluyámoslos al modelo en "cascada". Una manera de hacerlo es con varios modelos en una lista.

```{r}

#modelos

naive_ols_models <- list(
   " (1)" = lm( test_score ~ summerschool , data= analysisdata), 
    "(2)" = lm ( test_score ~ summerschool + parental_lincome, data= analysisdata ),
    "(3)" = lm ( test_score ~ summerschool + parental_schooling, data= analysisdata ),
    "(4)" = lm ( test_score ~ summerschool + parental_lincome + parental_schooling, data= analysisdata )      )

# Mapeo de nombre de variables

cm <- c( "test_score" = "Test Score",
         "summerschool" = "Summer school indicator",
         "parental_lincome" = "Parental Income (log)",
         "parental_schooling" = "Parental schooling" )
```

Tabla con stargazer:

```{r}
library(stargazer)

stargazer::stargazer(naive_ols_models, type = "text", out="tabla1.tex")

```

Tabla con modelsummary:

```{r}

modelsummary::modelsummary(naive_ols_models, stars= TRUE, fmt=3, 
                           estimate = "{estimate} ({std.error}){stars}", statistic=NULL,  gof_omit = "AIC|BIC|Lik",
                             coef_map = cm)

```

Sin embargo, si creemos que lo que también correlacona es una inobservable como la preocupación / involucramiento de los padres aun preservan el problema de la autoselección, todavía tenemos presencia de sesgo por variables omitidad.

¿Qué hacer?

Vamos a revisra varias alternativas: 

- Experimental

    1. Método de diferencias (en datos experimentales)

    2. Logit y probit

 - Datos observacionales: 

    3. Variables Instrumentales

    4. Regresión Discontinua

    5. Diferencias en Diferencias

    6. Efectos fijos, aleatorios y modelos jerárquicos

    7. Matching

## 2. Método de diferencias (en datos experimentales)

Un primer caso, es que antes de asistir al curso de verano se selecciona al azar un grupo de estudiantes para enviarles una carta de invitación que recuerda sus características y detalla como participar.

¿Podría ser este un buen experimento?

-   Asignación aleatoria
-   Tratamiento: recibir la carta -\> no es directamente la escuela de verano. Nos permite ver el efecto causal de recibir la carta, pero no de asistir a la escuela de verano. Vamos a usar una técnica para enfrentar esto posteriormente.

Por ahora, supongamos que conocer el efecto de la carta en el puntaje es una pregunta lo suficientemente interesante.

$testscore_i=\beta_0+\beta_1Letter_i+u_i$.

Este tipo de modelos podemos estimarlos con el paquete *linear modelos* `lm()` o también, con `feols()` del paquete *fixtest*

### 2.1 Estimación OLS

#### 2.1.2 Usando lm()

```{r}

#modelos

ols_letter_experiment <- list(
   " (1)" = lm( test_score ~ letter , data= analysisdata)    )

# Mapeo de nombre de variables

cm <- c( "test_score" = "Test Score",
         "summerschool" = "Summer School",
         "parental_lincome" = "Parental Income (log)",
         "parental_schooling" = "Parental schooling" )
```

Tabla con stargazer:

```{r}
library(stargazer)

stargazer::stargazer(ols_letter_experiment, type = "text", out=here::here("output/tabla2.tex"))

```

Podemos ver que los que recibieron la carta, tienen 0.13 más puntaje que los que no, manteniendo lo demás constante.

Si le creemos a la asignción aleatoria, entonces este efecto podría ser interpetado como casual.

Para saber si la signación aletaoria efectivamente separó a dos grupos comparables de tratados y controles lo que debemos hacer es probar el **balance en obbservables**. Esto quiere decir que al menos en las variables observables, efectivamente ambos grupos no son estadísticamente diferentes antes del tratamiento.

### 2.1 Balance en observables

#### 2.1 Comparación de a pares

Una primera opción es realizar una prueba t en las variables observables en ambos grupos, usamos el comando `t.tes()` . Por ejemplo podriamos comparar el ingreso de los padres con `parental_income~summerschool` para probar si la media de test_score es significantemente diferente en los dos grupos especificados por la asistencia al curso de veranp

```{r}
# Realizamos la prueba T
t.test(parental_lincome~summerschool,data=analysisdata)

```

Podemos notar que no hay diferencias significativas entre los grupos.

Esto tendríamos que realizarlo en todas las variables y se suele presentar en una **tabla de balance**

#### 2.3 Tabla de balance

Una tabla de balance presenta la diferencia de medias para grupos tratados y de control. Lo que deseamos es que no haya diferencias significativas antes del tratamiento en los grupos.

Podemos hacer una tabla que se vea bien usando `datasummary_balance()` del paquete *modelsummary*. Cargamos tambien el paquete *estimatr* porque le permite a `datasummary_balance()` calcular e inluir pruebas t de a pares directamente en la tabla.

```{r}
# Load libraries
library(modelsummary)
library(estimatr)
# Filter and modify data
testdata<-filter(analysisdata,year==5)
testdata<-ungroup(testdata)
testdata<-mutate(testdata,Treated=ifelse(letter==1,"Letter","No Letter"))
testdata<-select(testdata,female,parental_schooling,parental_lincome,test_score,Treated)
testdata<-rename(testdata,`Female`=female,
          `Parental schooling (years)`=parental_schooling,
          `Parental income (log)`=parental_lincome,
          `Test Score`=test_score)
     
# Table with balancing test
datasummary_balance(~Treated,
                    data = testdata,
                    title = "Balance of pre-treatment variables",
                    notes = "Notes: This is a brilliant table!",
                    fmt= '%.5f',
                    dinm_statistic = "p.value")
```

Los reslutados de la tabla sugieren que no hay diferencias significativas entre tratados y controles ANTES del tratamiento. ¡ Muy bien!

Y como antes, podríamos directaente agregar la opción `outpu=tab_balancing.docx` para exportarla en word (o latex). (no olvidar especificar la carperta con `here::here("output/`)

#### 2.4 Aproximación gráfica

Podemos ver el efecto causal propuesto en la siguiente gráfica. En los años antes del tratamiento son iguales y en los años posteriores, las distribuciones se empiezan a separar.

```{r}
# Load ggridges
library("ggridges")

# create a ggridges  chart
ggplot(analysisdata,aes(y=as.factor(year),x=test_score,fill=as.factor(letter) ))+
        geom_density_ridges(  alpha = .7, scale=1.5,color = "white", from = -2.5, to = 2.5)+
        theme_minimal()+
        theme_ridges(grid = FALSE)+
  scale_y_discrete(expand = c(0, 0)) +
  scale_x_continuous(expand = c(0, 0)) +
        scale_fill_brewer(palette="Set1",labels=c("No letter","Letter"))+
        labs(x="Test Score",y="Year",fill=" ",
            title="Test score distribution by reminder letter status  & year")+
        theme(legend.position="top",aspect.ratio=4/3,plot.title = element_text(hjust = 0.5))
```

## 3. Efectos fijos, aleatorios y modelos jerárquicos

Muchas veces, los datos estan agrupados en diferentes categorías o dimensiones.

En nuestro caso los estudiantes están agrupados en colegios y es esperable que dentro de un colegio los resultados de dichos estudiantes estén correlacionados. Esto también ocurre al tener datos temporales o cualquier forma de agrupación lógica de la información.

Podemos ignorar esta estructura en los datos y simplemente considerarlos como un conjunto de unidades comparables, es decir tratarlo como un **pooled cross-section**. Si las diferencias entre grupos no son sistemáticas, no deberíamos observar diferencias, sin embargo esto puede estar haciendo que caigamos en problemas como la *paradoja de Simpson* o el *Mathew effect*

Para profundizar recomiendo: https://ds4ps.org/pe4ps-textbook/docs/p-040-fixed-effects.html o el capítulo del libro

![](images/FEproblems.png)

En nuestro caso, tenemos por ejemplo 30 colegios en los datos. 

Otro punto es considerar que no solo el efecto está modulado por la organziación de agrupación de la data, sino también los errores. Esto lo llamamos clusterizar.

#### 3.1.2 Estimación usando lm()

Si cada grupo está afectando los resultados, calculemos el efecto el un OLS solo del año 5, un pooled en los años hasta el 5 y en este ultimo, incluyamos un efecto fijo por colegio  y luego por año.

```{r}

# No siempre necesitamos hacer una lista de
OLS_year5_data <-filter(analysisdata,year==5)

OLS_year5 = lm( test_score ~ letter , data= OLS_year5_data )


# Para el pooled, usemos todos los años antes del 5

POLS_data <-filter(analysisdata,year<=5)

Pooled_OLS = lm (  test_score ~ letter , data= POLS_data )

# Fixed effect por colegio

FE = lm (  test_score ~ letter + factor(school_id), data= POLS_data )

# Usemos stargazer para la tabla esta vez

stargazer( OLS_year5, Pooled_OLS, FE,
           type = "text", 
           dep.var.labels = ("Test score"),
           column.labels = c("Cross-Sectional OLS (year 5)", "Pooled OLS", "School FE"),
          covariate.labels = c("Intercept", "letter"),
           keep =c("Intercept", "letter"),
           omit.stat = "all", 
           digits = 2, intercept.bottom = FALSE )

```



#### 3.1.2 Estimación usando feols()

En el caso de que deseemos hacer un modelamiento más complejo, como opr ejemplo incluir efectos fijos (lo revisaremos al final de la sesión) o calcular los errores estándar en cluster podemos usar la función `feols()` del paquete *fixest*.

El paquete *fixest* permite la inclusión de múltiples efectos fijos y clusterización de errores es, nosotros usaremos la función `feols()`

Esta función espera una fórmula de especificación de la forma:

`y~x1+x2+..|fixed effects|IV specification,cluster=..`

Por ejemplo, en ese caso queremos clusterizar los errores a nivel de colegio:

```{r}
# Load packages
library(fixest)

# Select data
regdata<-filter(analysisdata,year==6)

# Regression
m1<-feols(test_score~letter+parental_lincome+female+parental_schooling, cluster="school_id",data=regdata)

# Summary of regression
summary(m1)

```

Quisieramos compararlo con el modelo sin dicha clusterización, podriamos hacer una mejor tabla. Aprovechamos la función `modelsummary()` y:

1.  Corremos ambas regresiones y guardamos los modelos en una lista. La primera es sin ajustar los errores, la segunda con cluster, la tercera con un efecto fijo por colegio.

2.  Usamos `modelsummary()`para el output de las regresiones.

3.  Especificamos que entre paréntesis esté el error estándar y que omita la constante en la tabla.

4.  Definimos el formato de los números con `fmt=...`

```{r}
# Load flextable and modelsummary
library(modelsummary)
# Select data

regdata<-filter(analysisdata,year==6)

# Regressions
models<-list(
  m1<-lm(test_score~letter+parental_lincome+female+parental_schooling ,data=regdata),
  m2<-feols(test_score~letter+parental_lincome+female+parental_schooling, cluster="school_id",data=regdata),
  m3<-feols(test_score~letter+parental_lincome+female+parental_schooling|school_id, cluster="school_id",data=regdata)
  )

# Generate table
modelsummary(models, stars = TRUE, statistic = 'std.error',
             fmt= '%.4f',
              coef_omit= '(Intercept)')
```

Podemos obsercar que no cambian los estimadores entre el modelo 1 y 2, pero si el error. Podemos ver que en el modelo 3 también cambia el estimador.

#### 3.1.3 Elementos clave de los efectos fijos:

Los elementos centrales de efectos fijos son:

- Usar efectos fijos nos permite controlar por todas las características invariantes en la unidad de analisis.

- Asumimos que algunas de las características de la entidad que estamos estudiando afectan nuestra prediccion y necesitamos controlar por estas.  Notar que efectos fijos nos permite entonces, capturar todas las caracteristicas que no varian en el tiempo en dicha entidad reduciendo indirectamente el sesgo de variables omitidas.

- Importante saber que no podemos estudiar el efecto de estas caracter+pisticas que no varian en el tiempo usando efectos fijos. Por ejemplo, no podríamos estudiar ni el impacto del sexo, ni de la educación de los padres en este ejemplo. Como dichas variables no cambian en el tiempo, la constante que agregamos en el efecto fijo absorve dicho efecto tambien, pues son perefectaente colineares.

-  Por esto, es importante tener en cuenta que los efectos fijos estan diseñados para estudiar características  que cambian en el tiempo para una entidad.

_ finalmente. efectos fijos tambien asume que cada entidad es única, tal que el error de dicha entidad es constante y no se correlaciona con el error de otra entidad. Es decir, asume independencia entre entidades. 

### 3.2 Random vs fixed effects

CUn supuesto subyacente de la implementación de efectos fijos es que no hay correlación entre los errores de las entidades de agrucación y las constantes, es decir, que cada colegio es indepediente del otro (en nuestro caso).

Este supuesto se puede testear estadísticamente mediante un modelo aletrnativo que, si bien calcula agrupando los datos permite que los efectos sean aleatorios, es decir calculamos asumiendo que las diferencias entre las entidades afectan a la variable dependiente y que hau correlación entre el error de las entidades y sus constantes. Podemos usar la función `plm()` del paquete plm especificando model="random".

Necesitamos que el panel esté bien conformado para poder aplicar este paquete, es decir tener una observación por unidad dimensionada y variables que cambien intra y entre grupos.

```{r}
#regdata <- analysisdata %>%
#            group_by(school_id, year)%>%
#            dplyr::summarize( test_score = mean(test_score), summerschool= mean(letter), female= mean(female), suma=n() ) %>%
#            ungroup()

#head(regdata)

```



```{r}
# cargamos la librería
#library(plm)

# Para especificar un modelo correctaente, necesitamos variables entre e intra entidades. EN este caso, agruparemos la data por colegio y año.

#regdata <- analysisdata %>%
#           group_by(school_id, year) %>%
#           summary( test_score = mean(test_score), summerschool=test_score ) %>%
#           ungroup()


#Fixed effects
#fixed <- plm( test_score ~ summerschool, data=regdata, index=c("school_id", "year"), model="within" )

#Random effects
#random <- plm( test_score ~ summerschool, data=regdata, index=c("school_id", "year"), model="random" )


```

Podemos comparar ambos modelos mediante la prueba de Hausman para decidir que modelo es mejor. Generalmente se prefiere efecto fijo ya que es más eficiente, si ambos están correctos.

Algunos otros elementos a tener en cuenta:
- Los efectos aleatorios permiten icnluir elementos que son invariables en el tiempo dado que asume que el error de la entidad no está correlacionado con el predictor.
- Los efectos aleatorios están influenciados por sesgo de variable omitida, ya que no estamos incluyendo directamente las dummies.


## 3. Logit y probit {#logit-y-probit}

Un segundo elemento es revisar que ocurre con probabilidad de entrar al preuniversitario dado que se recibió la carta.

En este caso, el modelo a estimar tendría como outcome de interés a una variable dicotómica.

$$ summerschool_i = g(letter_i)+u_i $$

### 3.1 Probabilidad lineal

Una primera aproximación sería sencillamente estimar un modelo lineal:

$$ summerschool_i = \beta_0+ \beta_1 letter_i + \gamma' X_i+u_i $$ Ya sabemos como hacerlo, simplemente usamos lm() o felm(). Aprovechemos esta oportunidad para agregar una fila con la media de la variable dependiente, este es genralemnte un elemento importante a incluir cuando estamos rabajando con variables dicotómicas ya que pone en contexto el análisis. Esto lo hacemos con el parámetro `add_rows=...`

```{r}

# Estimate LPM

models<-list(
  m1<-feols(summerschool ~letter, cluster="school_id",data=regdata),
  m2<-feols(summerschool ~letter+parental_schooling+parental_lincome+female,cluster="school_id",data=regdata)
  )

# Store the mean of dependent variable in a data frame
added_stats<-tibble("Mean of Dep. ",m1=mean(regdata$summerschool),m2=mean(regdata$summerschool))

# Generate table
modelsummary(models, stars = TRUE,statistic = 'std.error',  
             fmt= '%.4f',add_rows = added_stats,
              coef_omit= '(Intercept)', output = 'flextable')
```

Podemos ver que recubir la carta aumenta la probabilidad de asistir a la escuela de verano en 44%, lo que es un incremento de casi un 100% considerando que la media de dicha variable era 0.46.

A pesar de su gran interpretabilidad, los modelos de probabilidad lineal tienen varios problemas: predicciones pueden salir del rango \[0,1\] y que son heterocedasticos por definición.

Una alternativa a este modelamiento es usar una forma funciónal que asegure que la predicciones estén en el rango \[0,1\] es con logit o probit.

### 3.2 Logit y probit

Para estimar logit usamos `glm()` para estimar modelos logit y probit. Espeificamos la función en `family=binomial(link="probit")` o `family=binomial(link="logit")`.

Estimemos con logit:

```{r}
#  Estimate a binary outcomes model using a logit
logit_results <- glm(summerschool ~letter, data = regdata, family = binomial(link="logit"))

# Print the results
summary(logit_results)
```

El problema es que estos parámetros no son directamente interpretables porque son las pendientes en la función logistica. Para esto, calculamos los efectos marginales.

### 3.3 Efectos marginales

Usamos la librería margins para obtener los efectos marginales promedio. Simplemente usamos la función `margins()` en un objeto que tenga la estimación por probit o logit.

```{r}
# Load margins package
library("margins")

# Select data
regdata<-filter(analysisdata,year==6)

##  Estimate a binary outcomes model using a logit and probit
logit_results <- glm(summerschool ~letter+parental_lincome+female+parental_schooling, data = regdata, family = binomial(link="logit"))
probit_results <- glm(summerschool ~letter+parental_lincome+female+parental_schooling, data = regdata, family = binomial(link="probit"))

# Estimate linear probability model
lpm<-feols(summerschool ~letter+parental_schooling+parental_lincome+female, cluster="school_id",data=regdata)

#  Compute marginal effects
mfx_logit<- margins(logit_results)
mfx_probit<- margins(probit_results)

# Include in table
modelsummary(list("LPM"=lpm,"Logit"=mfx_logit,"Probit"=mfx_probit))
```

Acá podemos ver una tabla con los tres modelos y compararlos.

## 4. Variables Instrumentales

Volviendo a nuestra pregunta de investigación inicial queríamos conocer elefecto de asistir a la escuela de verano en los puntajes, no de la carta.

Sin embargo, sabemos que la carta se asignó aleatoriamente y que afecta a la probabilidad de asitir a la escuela de verano. En otras palabras: **podemos usar la carta como instrumento**

Es más, podríamos no realizar ninguna estimación y recurrir a lo ya obtenido:

-   Sabemos que el estimador de la forma reducida: $$ testscore_i = \beta_0 + \beta_1 letter_i + u_i $$ es 0.2 SD.}
-   En el punto anterior calculamos la primera etapa: $$ summerschool_o = \alpha_0 + \alpha_1 letter_i + \mu_i$$ que corresponde aproximadamente a 0.4.

Entonces, el **Estimador de Wald** de variables instrumentales es: IV=0.2/0.4 =0.5

Podríamos obtener este estimador directamente ahora desde los datos. Vamos a revisar dos métodos.

### 4.2 IV usando ivreg()

Primero usemos la función del paquete AER (Applied Economics with R). La sintaxis es: - `ivreg(outcome equation| exogenous variables, data)`

```{r}
# Load the AER package
library(AER)

# Select data
regdata<-filter(analysisdata,year==6)

# Estimate IV specification with
summary(ivreg(test_score~summerschool+female+parental_lincome+parental_schooling|
                         female+parental_lincome+parental_schooling+letter,data=regdata))
```

Observamos que el estimador es cercano a 0.5.

### 4.2 IV usando feols()

Podemos usar nuevamente `feols()` e incluir variables instrumentales. Además `felm()` también permite estimar especificaciones. En este caso no incluirmos efecto fijo, así que le poneos 0.

```{r}
# Estimate IV specification with felm 
m1<-feols(test_score~parental_lincome+female+parental_schooling| # Outcome eq.
           0|                                                    # Fixed effects
           summerschool~letter                                   # First stage
           ,cluster="school_id"                                  # Cluster var
            ,data=regdata)
# Summary of results
summary(m1)
```

Los coeficientes son iguales que antes, pero los errores estándar cambiaron ya que los clusterizamos al nivel del colegio. Notemos que `feols()` nos indica que el F-stat del isntrumento excluido es 687 (cumpliendo la condición de relevancia del isntrumento).

### 4.3 Tabla de reporte IV

Terminemos la sección de IV con una tabla que muestre la correlación pira entre la asistencia a la escuela de verano y puntajes (ingenuo), la forma reducida, la primera etapa y la estimación de variable instrumental

Notemos que usaremos `coef_map= c( )` para renombrar las variables en la tabla, espeficicar los coeficientes a incluiry permitir que dos coeficientes con diferente nombre aparezcan en la misma fila (summerschool y fit_summerschool)

```{r}
# Estimate OLS
OLS<-feols(test_score~summerschool+parental_lincome+female+parental_schooling,cluster="school_id",data=regdata)
# Estimate reduced form 
RF<-feols(test_score~letter+parental_lincome+female+parental_schooling,cluster="school_id",data=regdata)
# Estimate first stage 
FS<-feols(summerschool~letter+parental_lincome+female+parental_schooling,cluster="school_id",data=regdata)
# Estimate IV specification 
IV<-feols(test_score~parental_lincome+female+parental_schooling| # Outcome eq.
           0|                                                    # Fixed effects
           summerschool~letter                                   # First stage
           ,cluster="school_id"                                  # Cluster var
            ,data=regdata)
# Combine results
IVresults<-list("OLS"=OLS,"RF"=RF,"FS"=FS,"IV"=IV)
# Coefficients
cm <- c('fit_summerschool' = 'Summer School',
        'summerschool' = 'Summer School','parental_lincome' = 'Parental Income', 
        'letter' = 'Reminder letter', "female"="Female","parental_schooling"="Parental Schooling")
# Output Table
modelsummary(IVresults, stars = TRUE,statistic = 'std.error',  
             fmt= '%.4f',coef_map=cm, output =      )
```

La primera columna nos muestra que si ignoramos el problema de selección y etsimamos la regresión directamente entre el puntaje y la asistencia a la escuela de verano (condicional en los controles) nos dá un efecto promedio de 0.78 SD. Como sería esperable, este efecto está **sobrestimado** pueso incluye el sesgo por selección, lo que podemos comprobar al comprarlo con el estimador de IV de 0.47SD.

Las estadísticas descriptivas que estudiamos, sugieren que quienes asistieron a la escuela de verano ya tenían mejores puntajes antes del la escuela de verano y que sus padres tenían mayores ingresos y años de educación. Entonces, el coeficiente en la primera columna no solo captura el efecto propio de la escuela de verano, sino todas estas caracteristicas inobservables correlacionadas positivamente con los puntajes y con la asistencia la escuela de verano.

## 6. Diferencias en Diferencias

La estrategia que tomamos con IV es bastante convincente, porque se basa en un exterimento aleatorio (RCT) y el coeficiente "ingenuo" era mayor que el de IV, como era esperable.

¿Pero que hacemos si no contamos con una carta? Imaginemos que no tenemos la carta, pero si contamos con mediciones antes y despues de la escuela de verano (como lo hacemos). Podemos usar la estrategia de **DiD** .

Estaría dada por la siguiente especificación:

$$test_score=\beta1+\beta_2summerschool_i+\beta_3after_i+\beta_4after_i\times summerschool+\gamma'X_i+u_i$$

en el cual after es una variable dictotómica que considera el periodo después de asistir la escuela de verano.

### 6.1 Tendencias paralelas

Un supuesto clave en la estrategia de diferencias en diferencias es que, en ausencia del tratamiento los grupos habrian estado evolucionando "en paralelo". Con esto argumentamos que la evolución en los controles son un buen contrafactual a la evolución de los tratados.

```{r}
# Collapse data on the year X summer school level

df<-group_by(analysisdata,year,summerschool)

# Calculate the mean test score by year
df<-summarise(df,test_score=mean(test_score))

# Create a chart
ggplot(df,aes(x=year,y=test_score,color=as.factor(summerschool)))+
  geom_line(size=1)+
  theme_classic()+
  theme(legend.position = "top")

```

Observamos que previo a la escuela de verano iban "mas o menos" en paralelo. Pero después del tratamiento, estos caen para los controles. Esto es un poco preocupante, pero fuera de eso, se ven similares.

Generalmente, no contamos con tantas mediciones para tener una evaluación empirica de las tendencias.

### 6.2 Estimación modelo de diferencias en diferencias.

Primero, en un caso más similar, usemos solo las estimaciones justo antes y una después de la escuela de verano usando la función `lm()`.

Primero, creamos la variable dummy de después usando `mutate()` y `ifelse()`.

```{r}
# Select  year 5 and year 6 observations
df<-analysisdata%>%filter(year%in%c(5,6))


# Generate an after indicator and the interaction between after and summer school
df<-mutate(df,after=ifelse(year==6,1,0),afterXsumerschool=after*summerschool)

# Estimate the DiD regression and save the results in an element called m1
m<-lm(test_score~summerschool+after+afterXsumerschool,data=df)

# Print the results
summary(m)
```

Bien, observamos que el estiamodr de DIF es de 0.75SD. Observamos que es mayor que el que obtuvimos en IV y menor que el de OLS, pero por poco.

Veamos que ocurre si ahora usamos todos los datos disponibles. Podemos usar la función `feols()` para estimar con todos los años (también podriamos usar `lm()`). Este es el enfoque tradiciona, con una variable que indica después del tratamiento y la interacción entre después y tratamiento y una dummy por tratamiento. En el segundo modelo ponemos la otra opción, que es agregar un efecto fijo temporal en vez de la dummy de después. En el tercer modelo, agregamos controles.

```{r}

# Modify data
df<-mutate(analysisdata,after=ifelse(year>5,1,0),afterXsummerschool=after*summerschool)
df<-rename(df,school=summerschool)
# Estimate DiD without controls and fixed effects
m1<-feols(test_score~school+afterXsummerschool+after,cluster="school_id",data=df)
# Estimate DiD without controls
m2<-feols(test_score~school+afterXsummerschool|school_id+year,cluster="school_id",data=df)
# Estimate DiD  with controls
m3<-feols(test_score~afterXsummerschool+school+parental_lincome+female+parental_schooling
         |school_id+year,cluster="school_id",data=df)
# Leave year 5 out to avoid dip
m4<-feols(test_score~afterXsummerschool+school+parental_lincome+female+parental_schooling
         |school_id+year,cluster="school_id",data=filter(df,year!=5))

# Coefficients
cm <- c('afterXsummerschool'  = 'Summer School X After',
        'school'        = 'Summer School',
        'after'='After')
# Output Table
modelsummary(list(m1,m2,m3,m4), stars = TRUE,statistic = 'std.error',  
             statistic_override = vcov,
             fmt= '%.4f',coef_map=cm, output = 'gt'
        )

```

Podemos ver que el efecto de tratamiento es cercano a 0.6SD, lo que se acerca más al de IV y se aleja aun más que del OLS ingenuo.

La última columna (4) excluye el año 5 donde veiamos la desviación en las tendencias. Esto nos da un estimador aun más cerca del de IV.



## 7. Matching

Otra alternativa de análisis es mediante el matching. EN este caso, buscaremos para cada individuo que fue a la escuela de verano algún compañero similar en observables que no haya asisitido. 

Hay dos formas de hacerlo, buscar alguien que sea exactamente igual (exact matching) o que sean aproximadamente igual (approximate matching). Una de las formas más usadas de encontrar estos aproximadamente iguales es quela probabilidad de asistir a la escuela de verano basado en un set de observables sea igual (propensity score matching). Vamos a ver ambos tipos.

Una revisión completa de matching está acá: [ds4ps](https://ds4ps.org/pe4ps-textbook/docs/p-080-matching#matchit-function-in-r) o en [
Mixtape cap 5] (https://mixtape.scunning.com/05-matching_and_subclassification) la mayoría basados en el trabajo de [Abadie e Imbens (2006)](https://econpapers.repec.org/article/ecmemetrp/v_3a74_3ay_3a2006_3ai_3a1_3ap_3a235-267.htm)

La estrategia de matching sigue los siguientes dos pasos:
1 -  Para cada tratado encontramos a su "clon". Para esto hay dos alternativos, usamos uno (o varios) match exactos o usamos el propensity score matching.

2 -  Usando los "clones" calculamos calculamos el ATE.

Lo primero que deberíamos hacer entonces es obtener los clones. 

### 7.1 Hayando clones manualmente

Primero lo haremos manualmente. 

#### 7.1.1 Exact matching manual

Un exact match, lo que hace es identificar uno (varios) contrafactuales que se encuentren más cercanos a cada individuo tratado en la muestra.

Lo que hacemos es definir en que variables vamos a considerar el match. En nuestro caso quisieramos hacerlo por colegio, sexo, educacion de los padres e ingreso de los padres.

Para simplificar como se hace el algortmo, primero usemos una única característica en un único colegio, educación de los padres.


```{r}
# Realizamos un filtro de la data para ilustrar el ejemplo en un solo colegio
# Dado que la educaión de los padres la tenemos solo para un año, quedémonos con los datos del año 6 para ilustrar con el outcome de interes el test_score después del tratamiento.

exact_macth_df_treated <- analysisdata %>%
                          dplyr::filter(school_id==1&year==6&summerschool==1) %>%
                          dplyr::select(person_id, summerschool, parental_schooling, test_score)
                    
exact_macth_df_controls <- analysisdata %>%
                          dplyr::filter(school_id==1&year==6&summerschool==0) %>%
                          dplyr::select(person_id, summerschool, parental_schooling, test_score)



head(exact_macth_df_treated)
head(exact_macth_df_controls)

```

Entonces, para cada individuo tratado encontremos el control en su colegio que sería más cercano en términos de educación de sus padres.

Varias cosas que debemos decidir:

1. Como definiremos "cercano" : mediremos una distancia, una diferencia, otra cosa.
2. Cuantos matchs "más cercanos" buscaremos
    2.1 Si es solo 1, cual elegimos. (ventaja clones más similares)
    2.2 Si es más de uno: Los guardamos todos o realizamos una agregación. (permite mayor muestra)
        2.2.1: ¿Cómo agregamos?: la media, la mediana... otro?
3. Si usamos un miembro de la muestra, ¿lo podemos volver a usar?

En la [documentación del paquete](https://cran.r-project.org/web/packages/MatchIt/vignettes/matching-methods.html) `MatchIT` de Noah Greifer hay un buen listado y resumen de los diferentes métodos posibles para hacer el matching.

Hagamos un match con el más cercano en el valor absoluto de la diferencia. Esto se llama *" nearest neighbor matching"*


Si hay más de un candidato, elegimos el primero.

Pro ejemplo, calculamos la distancia para cada elemento de los controles con el individuo 1 de los tratados.

```{r}

exact_macth_df_controls_t1 <-  
    exact_macth_df_controls %>%
    mutate(distance_1 = abs(11-parental_schooling) )

head(exact_macth_df_controls_t1)
```
Acá podemos elegir el que la distancia sea 0 o cualquier valor arbitrariamente pequeño. Por ejemplo, elijamos todos los que el valor es 0.

```{r}
exact_macth_df_controls_t1 <- exact_macth_df_controls_t1 %>%
                              dplyr::filter(distance_1==0)

head(exact_macth_df_controls_t1)
dim(exact_macth_df_controls_t1)
```

Podemos ver que hay 24 candidatos. ¿Cómo elejimos? Imaginemos que elegimos uno al azar y lo defifinimos como el clon para el individuo 1 de los tratados.


```{r}
set.seed(0) #para que todos tengamos el mismo resultado
slice_sample(exact_macth_df_controls_t1,n=1)

```

Entonces nuestro clon para el individuo 66 sería el 2421. Si recordamos su resultado era 3.571861 y del match, decimos que su contrafactual sería 2.372809, entonces su ATE sería:

```{r}
3.571861 - 2.372809
```

Luego pasamos al segundo individuo tratado. Para esto primero debemos saber si la persona 2421 estará disponible como match potencial. Una vez definido eso, replicamos los pasos anteriores.

Con esta idea es posible crear algorthmos adhoc a cada problema, pero nosotros aprovecharemos un paquete que ya los tiene implementados: `MatchIT`


#### 7.1.2 Propensity score matching

Del ejemplo anterior es facil imaginarse como sería una alternativa para multiples dimensiones, como una función con múltiples imputs que describa alguna distancia.

Una segunda manera de hacer el match es con el propensity score. Una forma comun es mediante una función logit. A cada individuo vamos a calcular su probabilididad de ser tratado dadas las covariables en las que deseamos hacer el match.


```{r}
# Realizamos un filtro de la data para ilustrar el ejemplo en un solo colegio
# Dado que la educaión de los padres la tenemos solo para un año, quedémonos con los datos del año 6 para ilustrar con el outcome de interes el test_score después del tratamiento.

props_macth_df<- analysisdata %>%
                          dplyr::filter(school_id==1&year==6) 
                    

head(props_macth_df)

```

Calculamos el propensity score estimando un modelo logit y luego calculamos la probabilidad estimada.

```{r}
## Estimamos el logit model

log_matching = glm( summerschool ~ female + parental_schooling + parental_lincome, data = props_macth_df, family = "binomial")

## Calculate the propensity score by estimating the probability for each child to attend a private or a public school.

props_macth_df$prop_score <- predict(log_matching, newdata = props_macth_df, type = "response")

# creamos un indice para identificar alumnos más facil que su id (es muy larga)
props_macth_df <- props_macth_df %>% dplyr::mutate(number = row_number())


head(props_macth_df %>% dplyr::select(person_id, summerschool, female, parental_schooling, parental_lincome, prop_score, number ) )

```

Podemos ver el rango en el que se mueve prop_score como referencia

```{r}
summary(props_macth_df$prop_score)

```

Podriamos también hacer una gráfica de cada estudiante y su propensity score

```{r}

palette(c("darkred", "dodgerblue1")) #azules tratados y rojos no tratados

plot(prop_score ~ number, data=subset(props_macth_df, number<20), lwd = 3, ylab = "Predicted probabilities", xlab = "Observations", xaxt = "n", col=summerschool+1 )

axis(1, at=1:20)

for (i in 1:nrow(props_macth_df)){
  segments(props_macth_df$number[i], 0, props_macth_df$number[i], props_macth_df$prop_score[i], lty = 2, col = "gray70")}

```

Entonces, **¿cómo buscamos los clones?** 

Ordenemos nuestro dataframe por propensity score, para ver cuales son similares:

```{r}
## Ordenemos el data frame por propensity score 

props_macth_df = as.data.frame(props_macth_df[order(props_macth_df$prop_score), ])

## We just look at the top 5 observations
props_macth_df_2 = round(props_macth_df[5:9,], 2)
props_macth_df_2
```
Por ejemplo en este grupo todos tienen propensity score similar y solo uno de estos fue tratado.


Para elegir su match podriamos calcular la distancia y definir bajo que valor definiremos un match.


### 7.2 Clones con paquete MatchIt

Una alternatia al trabajo manual es usar el paquete [*MatchIT* ](https://kosukeimai.github.io/MatchIt/articles/MatchIt.html) (Muy recomendable revisar su documentación)


```{r}
#cargar el paquete

library("MatchIt")
```

Partiremos con definir un proceso de matching por defecto con todas las observaciones y definiendo como covariables school_id + female + parental_schooling + parental_lincom + test_score (año 5)

Subset para matching solo con el periodo antes de la escuela de verano

```{r}
match_data <- school_data_1 %>%  #usamos este porque quisieramos tener en una fila el outcome del test en año 6 y del 5 para usar de control 
             dplyr::filter(!is.na(parental_schooling),!is.na(test_year_5),!is.na(test_year_6)) %>%
             dplyr::mutate(summerschool=summercamp)  
            

match_data
```

De estos, en total 1598 individuos son tratados.

```{r}
mean(match_data$summerschool)
sum(match_data$summerschool)
```
Recordemos que los grupos tratados y tratamiento eran bastante diferentes:


En ingreso de los padres

```{r}

ggplot(filter(match_data),
       aes(x=parental_lincome,fill=as.factor(summerschool)))+
       geom_histogram(aes(y=..density..),bins = 50,alpha=0.5,
                       position="identity",color="white")+
       geom_density(alpha=0.0,size=1,show.legend= FALSE)+
       theme_minimal()+
       labs(y="Densidad",x="Ingreso de los padres (log)",fill=" ")+
       scale_fill_brewer(palette="Set2",labels=c("No asistió","Asistío a la escuela de verano"))+
       theme(legend.position="top")

```


En la educación de los padres


```{r}

t.test(match_data$parental_schooling[match_data$summerschool==1],match_data$parental_schooling[match_data$summerschool==0], paried=T)  

```

En los resultados previos


```{r}

ggplot(filter(match_data),
       aes(x=test_year_5, fill=as.factor(summerschool)))+
       geom_histogram(aes(y=..density..),bins = 50,alpha=0.5,
                       position="identity",color="white")+
       geom_density(alpha=0.0,size=1,show.legend= FALSE)+
       theme_minimal()+
       labs(y="Densidad",x="Resultados año 5",fill=" ")+
       scale_fill_brewer(palette="Set2",labels=c("No asistió","Asistío a la escuela de verano"))+
       theme(legend.position="top")

```


Usemos ahora la función `matchit()`:esta crea el propensity score y asigna los clones segun las especificaciones que le demos.

Notar el código de la función: primero especificamos el modelo logit en base al cual se calcularán los propensity scores. Luego especificamos el método para emparejar los scores. Nosotros usaremos vecino más cercano y por defecto trabaja sin reemplazo.


```{r}

# Definamos un matching con todo por defecto.

match<- MatchIt::matchit( summerschool ~ as.factor(school_id) + female + parental_schooling + parental_lincome+ test_year_5,  method = "nearest", data= match_data )

match

```

Observemos que lo que obtenemos es un objeto de matchit usando por defecto el vecino más cercano 1:1 y sin reemplazo.

Nos ha obtenido 3196 matchs del total de 3480 filas. 

Una vez el agortmo ha calculado los propensity scores, podemos crear un nuevo dataset que contenga las observaciones con matching usando la función `match.data()` donde ahora tenemos una medida de distancia, subclass que indica los pares.


```{r}
matched_data <- match.data(match)

head(matched_data)

```

Los elementos clave que agrega este proceso son weights (pesos computados), subclass (que indica los pares de match), distance (el propensity score estimador) y match.matrix (una matriz de maching).

###7.3 Evaluando la calidad del match

Matching suele ser un proceso efectivo para eliminar las diferencias entre grupo tratado y de control para alcanzar el balance en las covariables, su desemepeño debe ser evaluado.

Si las covariables siguen desvalanceadas después de emparejar, el matching no fue exitoso y se deben probar diferentes especificaciones. 

Lo primero es revisar si hay diferencias significativas entre tratados y controles:

```{r}
t.test(matched_data$parental_schooling[matched_data$summerschool==1],matched_data$parental_schooling[matched_data$summerschool==0], paried=T)  

```

Y vemos que efectivamente hay...

**Otras opciones en `matchit()`:

La función matchit tiene varias opciones que pueden explorar mejir en su [documentación](https://r.iq.harvard.edu/docs/matchit/2.4-15/User_s_Guide_to.html). Algunas de las más importantes son:

- replace: perimte que se realice el procedimiento con y sin reemplzao. Si deseas que sea con reemplazo se espcifica replace=TRUE
- ratio: si el matching se realiza con reempzalo, se puede especificar el número de controles para caso tratado (uno a muchos). Por ejemplo, ratio=5
- caliper: Se puede identificar una distancia espcifica para el match especificando las desviaciones etsándar en la distancia que son aceptadas. Tambien se puede indicar que, si no hay matches dentro de la distancia especificada, se use el vecino más cercano usando calclosest=TRUE


**Otra opción del matching:**


```{r}

# Definamos un matching con todo por defecto.

match<- MatchIt::matchit( summerschool ~ as.factor(school_id) + female + parental_schooling + parental_lincome+ test_year_5,  method = "nearest", data= match_data, replace=TRUE, ratio=5 )

match

```




```{r}
matched_data2 <- match.data(match)

head(matched_data2)

```


Chequemos que pasó ahora:


```{r}
t.test(matched_data2$parental_schooling[matched_data2$summerschool==1],matched_data2$parental_schooling[matched_data2$summerschool==0], paried=T)  

```


Podemos ver que todavia hay diferencias significativas entre tratados y controles.

El paquete ofrece varias herramientas para evaluar el balance después del match.

Podemos inpeccionar mejor este objeto con summary. Esta nos entrega la tabla de balance para los datos y mayores datos del match:

```{r}
summary(match)
```


Tambien podemos graficar el balance:


```{r}
plot(summary(match))
```


Otro elemento a considerar es que ocurre con el tamaño de muestra después del match, si se descartan demasiado unidades puede ser que reduzcamos las estimaciones y pasen a ser imprecisas, esto porque hay un trade off entre balance y tamaño de muestra resultante.


```{r}
plot(match, type = "jitter", interactive = FALSE)
```

Tambien podemos examinar bvisualmente usando gráficos de densidad.

```{r}
plot(match, type = "density", interactive = FALSE,
     which.xs = ~female + parental_schooling + parental_lincome)
```


### 7.4 Efecto tratamiento de matching

Ahora que tenemos un dataset emparejado el dataset podemos calcular el efecto del tratamiento. 

Una opción es estimar una regression, teniendo cuidado de incluir los pesos de los pares.

```{r}

# Estimate LPM

models<-list(
regression_match_nocontrols <-lm (test_year_6 ~ summerschool, data=matched_data ),
regression_match_controls  <- lm (test_year_6 ~ summerschool + as.factor(school_id) + female + parental_schooling + parental_lincome+ test_year_5, data=matched_data, weights=weights )
  )


# Generate table
modelsummary(models, stars = TRUE,statistic = 'std.error',  
             fmt= '%.4f',add_rows = added_stats,
              coef_omit= '(Intercept)', output = 'flextable')

```

Podemos usar la librería `marginaleffects`para  computar el ATT, se recomienda emplear errores robustos cluster para la mayoría de los análisis usando el par como clustering.

```{r}
library("marginaleffects")

comp  <- comparisons(regression_match_controls,
                     variables = "summerschool",
                     vcov = ~subclass,
                     newdata = subset(matched_data, summerschool == 1),
                     wts = "weights")

summary(comp)
```

Ojo:  los coeficientes de este modelo no deberían ser reportados ni interpretados directamente sin consideraciones. Se recomienda revisar [la viñeta del paquete](https://kosukeimai.github.io/MatchIt/articles/estimating-effects.html) para mayores detalles.

## 9. Regresión Discontinua

La última técnica que usaremos el la regresión discontinua. Vamos a hacer un ejemplo muy corto, pero necesitaremos más datos.

Pensemos que para asistir al curso de verano, los estudiantes tenían que dar una prueba diagnóstico y se invitaron a los que scaaron hasta 500 puntos. Estos datos los vamos a simular.

```{r}

```





# IV. Presentación de los resultados

El ultimo elemento que vamos a preparar son un conjunto de tablas que puedan dar cuenta de nuestros resultados.

Para esto usaremos dos paquetes: stargazer & modelsummary