assignment_5

## Load data
#| echo: false
#| message: false
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.1     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(haven)
library(gt)
library(modelsummary)
library(fixest)
library(plm)

Attaching package: 'plm'

The following objects are masked from 'package:dplyr':

    between, lag, lead

library(kableExtra)

Attaching package: 'kableExtra'

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

    group_rows

library(sandwich)
library(lmtest)

Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

library(tidymodels)

── Attaching packages ────────────────────────────────────── tidymodels 1.0.0 ──
✔ broom        1.0.4     ✔ rsample      1.1.1
✔ dials        1.2.0     ✔ tune         1.1.1
✔ infer        1.0.4     ✔ workflows    1.1.3
✔ modeldata    1.1.0     ✔ workflowsets 1.0.1
✔ parsnip      1.1.0     ✔ yardstick    1.1.0
✔ recipes      1.0.5     
── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ plm::between()           masks dplyr::between()
✖ scales::discard()        masks purrr::discard()
✖ dplyr::filter()          masks stats::filter()
✖ recipes::fixed()         masks stringr::fixed()
✖ kableExtra::group_rows() masks dplyr::group_rows()
✖ plm::lag()               masks dplyr::lag(), stats::lag()
✖ plm::lead()              masks dplyr::lead()
✖ yardstick::spec()        masks readr::spec()
✖ recipes::step()          masks stats::step()
• Use tidymodels_prefer() to resolve common conflicts.

library(ivreg)

ajr_df <- haven::read_dta(here::here("data", "AJR.dta"))
mdve_df <- haven::read_dta(here::here("data", "mdve.dta"))

Part 1

Lgdp = Log GDP per capita in 1995 Prot = protection against risk of expropriation index as a proxy for institutions (e.g., property rights and checks against government power) Logmort = Log settler mortality

1. Regress the log GDP per capita on the average protection against expropriation risk.

\(log(GDP_i) = \alpha + protection_i + \epsilon+i\)

mod_1 <- lm(lgdp ~ prot, data = ajr_df)
modelsummary(mod_1, stars = TRUE)

	(1)
(Intercept)	4.660***
	(0.409)
prot	0.522***
	(0.061)
Num.Obs.	64
R2	0.540
R2 Adj.	0.533
AIC	142.3
BIC	148.8
Log.Lik.	−68.168
F	72.816
RMSE	0.70
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

ajr_df |> 
  ggplot(aes(lgdp, prot)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

2. Regress the average protection against expropriation risk on (log) settler mortality. Call this π

mod_2 <- lm(prot ~ logmort, data = ajr_df)
modelsummary(mod_2, stars = TRUE)

	(1)
(Intercept)	9.400***
	(0.612)
logmort	−0.621***
	(0.127)
Num.Obs.	64
R2	0.278
R2 Adj.	0.266
AIC	215.0
BIC	221.5
Log.Lik.	−104.504
F	23.816
RMSE	1.24
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

ajr_df |> 
  ggplot(aes(prot, logmort)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

3. Regress the log GDP per capita on settler mortality. Call this γ.

mod_3 <- lm(lgdp ~ logmort, data = ajr_df)
modelsummary(mod_3, stars = TRUE)

	(1)
(Intercept)	10.708***
	(0.375)
logmort	−0.570***
	(0.078)
Num.Obs.	64
R2	0.462
R2 Adj.	0.454
AIC	152.3
BIC	158.8
Log.Lik.	−73.160
F	53.342
RMSE	0.76
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

ajr_df |> 
  ggplot(aes(lgdp, logmort)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

pie <- broom::augment(mod_2)$.fitted
gamma <- broom::augment(mod_3)$.fitted

ratio = (gamma / pie)

4. Compute the 2SLS regression of log GDP per capita on the average expropriation risk, using settler mortality as an instrumental variable

1st stage

first_stage <- lm(logmort ~ prot, data = ajr_df)
fs_resid <- first_stage$residuals
ajr_df <- bind_cols(ajr_df, fs_resid) |> 
  rename(fs_resid = ...9)

New names:
• `` -> `...9`

2nd stage

second_stage <- lm(lgdp ~ fs_resid, data = ajr_df)
summary(second_stage)

Call:
lm(formula = lgdp ~ fs_resid, data = ajr_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.54288 -0.66451 -0.03718  0.75594  2.02504 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.0622     0.1234   65.33   <2e-16 ***
fs_resid     -0.3397     0.1175   -2.89   0.0053 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9873 on 62 degrees of freedom
Multiple R-squared:  0.1187,    Adjusted R-squared:  0.1045 
F-statistic: 8.353 on 1 and 62 DF,  p-value: 0.005299

5. Compute a 2SLS regression of log GDP per capita on the average expropriation risk, using settler mortality as an instrumental variable

mod_iv <- ivreg(lgdp ~ prot | logmort, data = ajr_df)
summary(mod_iv)

Call:
ivreg(formula = lgdp ~ prot | logmort, data = ajr_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.40124 -0.55430  0.07021  0.66539  1.64168 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.0869     0.9853   2.118   0.0382 *  
prot          0.9171     0.1502   6.106 7.35e-08 ***

Diagnostic tests:
                 df1 df2 statistic  p-value    
Weak instruments   1  62     23.82 7.76e-06 ***
Wu-Hausman         1  61     21.23 2.13e-05 ***
Sargan             0  NA        NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9223 on 62 degrees of freedom
Multiple R-Squared: 0.231,  Adjusted R-squared: 0.2186 
Wald test: 37.29 on 1 and 62 DF,  p-value: 7.354e-08 

Part 2

We will use mdve.dta, which contains the case-level information. fail = A dummy variable for re-occurrence of a domestic assault within 6 months z_coddle = A dummy variable to indicate that the random assignment is advice or separate d_coddle = A dummy variable to indicate that the delivered assignment is advice or separate

1. Estimate the effect of coddle on recidivism by the OLS estimation of the equation

mod_4 <- lm(fail ~ d_coddle, data = mdve_df)
modelsummary(mod_4, stars = TRUE)

	(1)
(Intercept)	0.147***
	(0.025)
d_coddle	0.112*
	(0.048)
Num.Obs.	314
R2	0.017
R2 Adj.	0.014
AIC	288.6
BIC	299.8
Log.Lik.	−141.290
F	5.496
RMSE	0.38
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

2. Estimate the first-stage effect of the coddle assignment on delivering coddle by estimating

mod_5 <- lm(d_coddle ~ z_coddle, data = mdve_df)
modelsummary(mod_5, stars = TRUE)

	(1)
(Intercept)	0.011
	(0.043)
z_coddle	0.386***
	(0.052)
Num.Obs.	314
R2	0.152
R2 Adj.	0.149
AIC	344.9
BIC	356.2
Log.Lik.	−169.471
F	55.748
RMSE	0.42
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

3. Estimate the intention-to-treatment (ITT) effect of coddle relative to arrest by estimating:

mod_6 <- lm(fail ~ z_coddle, data = mdve_df)
modelsummary(mod_6, stars = TRUE)

	(1)
(Intercept)	0.098*
	(0.040)
z_coddle	0.114*
	(0.047)
Num.Obs.	314
R2	0.018
R2 Adj.	0.015
AIC	288.3
BIC	299.5
Log.Lik.	−141.127
F	5.827
RMSE	0.38
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

4. Compute the Wald estimate

beta <- broom::augment(mod_5)$.fitted
gamma <- broom::augment(mod_6)$.fitted

tau = (beta / gamma)

5. Estimate and interpret the 2SLS estimate of the effect of coddle on recidivism

mod_iv_2 <- ivreg(fail ~ z_coddle | d_coddle, data = mdve_df)
summary(mod_iv_2)

Call:
ivreg(formula = fail ~ z_coddle | d_coddle, data = mdve_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.26162 -0.26162 -0.26162  0.02261  1.02261 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.02261    0.09014  -0.251   0.8022  
z_coddle     0.28423    0.12368   2.298   0.0222 *

Diagnostic tests:
                 df1 df2 statistic p-value    
Weak instruments   1 312    55.748 8.3e-13 ***
Wu-Hausman         1 311     2.339   0.127    
Sargan             0  NA        NA      NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3884 on 312 degrees of freedom
Multiple R-Squared: -0.02268,   Adjusted R-squared: -0.02596 
Wald test: 5.281 on 1 and 312 DF,  p-value: 0.02221