An increase of thousand slave exports from an ethnic group is associated with an expected -510^{-4} decrease in trust in neighbors on a scale of 0–3. This effect is statistically significant (\(p < -.001\)).

summary(mod_1_0)

## 
## Call:
## lm(formula = trust_neighbors ~ exports, data = nunn)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7891 -0.7889  0.2109  1.2109  1.6761 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.789e+00  7.719e-03  231.77   <2e-16 ***
## exports     -5.441e-04  3.343e-05  -16.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.003 on 20578 degrees of freedom
##   (1242 observations deleted due to missingness)
## Multiple R-squared:  0.01271,    Adjusted R-squared:  0.01267 
## F-statistic:   265 on 1 and 20578 DF,  p-value: < 2.2e-16

In order to plot the fitted values against the scatterplot, I use broom::augment() with a newdata argument to calculate predicted values for each value of exports from the minimum to the maximum.

mod_1_0_fitted <-
  augment(mod_1_0,
          newdata = data_frame(exports = seq(min(nunn$exports, na.rm = TRUE),
                                             max(nunn$exports, na.rm = TRUE),
                                             length.out = 100)))

Now I plot the predicted values of the regression, the 95% confidence interval, and the original data. Since there are over 20,000 observations that only take a few values of trust_neighbors, a scatter plot will not work well since there will be a lot of overplotting. There are several options to deal with this (alpha, jittering), but for this example I will use hexagon binning; which is like a 2D histogram. I also use the viridis color scale.

ggplot() +
  geom_line(data = mod_1_0_fitted,
                     mapping = aes(x = exports, y = .fitted)) +
  geom_ribbon(data = mod_1_0_fitted,
              mapping = aes(x = exports, y = .fitted,
                            ymin = .fitted - 2 * .se.fit,
                            ymax = .fitted + 2 * .se.fit),
              alpha = 0.3) +
  geom_hex(data = na.omit(select(nunn, exports, trust_neighbors)),
             aes(x = exports, y = trust_neighbors)) +
  ylab("Trust in Neighbors") +
  xlab("Slave exports (thousands)") + 
  scale_fill_viridis(direction = -1, option = "inferno") +
  theme_minimal()

The first noticeable thing about the residuals is that they are arranged in lines; this is a result of the discrete values that the outcome variable takes.

ggplot(augment(mod_1_0), aes(x = .fitted, y = .std.resid)) +
  geom_point() + 
  ylab("Studentized residuals") +
  xlab(expression(hat(y))) +
  theme_minimal()

Since the outcome variable is bounded at 0 and 3, the residuals are not symmetric. However, since the fitted values only cover a range between 1 and 2, these asymmetry is not too great. Although we are using an OLS regression for a discrete and bounded outcome variable, this does not appear to be too big of an issue.

The null hypothesis of the t-test reported by summary.lm() is \(\beta_{\mathtt{exports}} = 0\). The \(p\)-value is the probability of observing a \(t\)-statistic equal to or more extreme than what was observed in the sample under repeated samples drawn from the same population, if the null hypothesis were true. No, the \(p\)-value is not the probability that the null hypothesis is correct. In Frequentist statistics, only samples have probabilities, hypotheses do not have a probability associated with them. They are either true of false, but the truth value is unknown. Even in Bayesian inference, the \(p\)-value gives the probability of the data given the hypothesis, and to calculate the probability of the hypothesis given the data we would need the prior probability of the null hypothesis.

Anything in the article which increases the plausibility of their theory, e.g. prior literature and the mechanisms of their theory, increases the prior probability of the research hypothesis, \(p(H_a)\). This is the difference between this paper, and a paper which tests ESP.

P-value hacking makes is so that the probability of the observing the data in the paper does not depend on the truth of the hypothesis which would generate these data. In other words the probability of observing a statistically significant result is 1, \(p(\text{data}) = 1\). This means that a p-hacked paper should have no affect on our beliefs about a hypothesis since there is no relationship between the data observed and the truth of the hypotheses.

An increase of one thousand exported slaves is associated with an expected decrease of -710^{-4} in the response to trust in neighbors, controlling for individual, district, and country level differences.

The regression coefficients in the version with and without controls are not too different, -0.0005 vs. -0.00067. Coefficient statibility to the inclusion of control variables is used as a heuristic for the potential importance of OVB. More formally, Nunn and Wantchekon use the method developed in Bellows and Miguel (2009) in the Section " Using Selection on Observables to Assess the Bias from Unobservables" on p. 3237. The statistic itself is simple: \(\hat{\beta}_f / (\hat{\beta}_r - \hat{\beta}_f)\) where \(\hat{\beta}_f\) is the coefficient from the full regression (including controls), and \(\hat{\beta}_r\) is the coefficient from the restricted regession (excluding controls):

```r
beta_f <- coef(mod_1_1)[["exports"]]
beta_r <- coef(mod_1_0)[["exports"]]
beta_f / (beta_r - beta_f)
```

```
## [1] -5.030797
```
The interpretation of this is that the covariation between the omitted variables and slave exports would have to be five times larger than the covariation between slave exports and the controls included in this regression. 

The \(R^2\) and number of observations match those reported in the table.

round(summary(mod_1_1)$r.squared, 2)

## [1] 0.16

length(fitted(mod_1_1))

## [1] 20027

I had meant to ask whether the standard errors also match those reported in Table 1, but I moved that to a later question. They do not. The coefficients match, but the standard errors are much smaller for the reasons discussed below.

The purpose of asking you to calculate the fitted values manually is so that you understand some of the matrix algebra in You can calculate the design matrix, \(X\), of a regression in R using model.matrix(). This will expand the formula for the regression, e.g. turning factors into dummy variables, and return a numeric matrix.

X <- model.matrix(mod_1_1)
beta <- coef(mod_1_1)

Now calculate \[ \hat{y} = X \hat\beta = \begin{bmatrix} \hat\beta_0 \cdot 1 + \hat\beta_1 x_{1,1} + \dots + \hat\beta_{K,1} \\ \hat\beta_0 \cdot 1 + \hat\beta_1 x_{1,2} + \dots + \hat\beta_{K,2} \\ \vdots \\ \hat\beta_0 \cdot 1 + \hat\beta_1 x_{1,N} + \dots + \hat\beta_{K,N} \end{bmatrix}. \]

yhat <- X %*% beta
head(yhat)

##       [,1]
## 1 1.342690
## 2 1.259782
## 3 1.322990
## 4 1.292238
## 5 1.379541
## 6 1.411596

We can check that these are the fitted values,

all.equal(fitted(mod_1_1), as.numeric(yhat), check.attributes = FALSE)

## [1] TRUE

Note, I use all.equal so that the comparison will allow the two vectors to be nearly equal (something which should always be done with floating point numbers due to rounding error). The argument check.attributes ignores that one vector has names and the other does not, because we only care that they have the same numerical values.

If we were to plot the fitted values against exports in a scatter plot we would have to consider how to handle the control variables. Without control variables we could plot the fitted values of the regresion directly against actual data. Each fitted value has different values of \(x\). One way would be to plot the fitted values at representative values of the control varaibles, e.g. means of numeric variables, and modes of factor variables. Alternatively, we could average over all the control variables. If you are only interested in how well the regresion fits at different values of the outcome variable, plotting the residuals may make more sense.

The answers of individuals in the same country, ethnicity, district, or living in close geographic proximity may be correlated even after controlling for the variables included in these models.

There wasn’t a particularly good way to do this. The idea was to get an intuition that clustered standard errors means that groups have correlated residuals. In this case, it is likely that individuals within ethnicities or districts are correlated.

ggplot(augment(mod_1_1, data = nunn) %>%
         group_by(ethnicity) %>%
         summarize(err = mean(.std.resid),
                   std.err = mean(.std.resid) / (sd(.std.resid) / sqrt(n())), obs = n()) %>%
         filter(obs > 5),
       aes(x = ethnicity, y = err, ymin = err - 2 * std.err, ymax = err + 2 * std.err)) + geom_pointrange()

In Table 1 (p. 3232), three types of standard errors are reported

1. Clustered by ethnic groups
2. Clustering within ethnic groups and within districts
3. Standard erorrs adjusted for two-dimensional spatial autocorrelation.

The standard errors reported in Table 1 are much larger than those in the OLS regression that we ran earlier. The reason is that the classical OLS errors do not account for correlation within these clusters. Effectively, by ignoring that correlation we are pretending we have many more observations that we do. See A Practitioner’s Guide to Cluster-Robust Inference or the discussion of clustered standard errors in Angrist and Pischke “Mostly Harmless Econometrics”. The easiest way to implement cluster robust errors in R is to use the plm package.

Roughly, a one percentage point increase in slave exports in an ethnic groups is associated with an expected decrease of 0.0074 (-0.74/100) in the the response for trust in neighbors of members of that ethnicity, controlling for individual, district, and country characteristics.

Many ethnicities have values of zero for slave exports, and \(\log(0)\) is undefined. Adding a small positive number is a common way to adjust variables that seem to have a logarithmic (decreasing returns) effect, but include zeros.

The rest of this question was somewhere between poorly worded and incoherent. My intent was to have you plot the functional forms in each regression.

The sampling distribution of \(\hat{\beta}_{\mathtt{ln\_export\_pop}}\) is approximately normal and centered at its OLS value of around -0.75.

ggplot(filter(results, term == "ln_export_pop"),
       aes(x = estimate)) +
  geom_density() + geom_rug()

The standard deviation of the estimated sampling distribution of the coefficient of ln_export_pop is

filter(results, term == "ln_export_pop") %>%
  magrittr::extract2("estimate") %>%
  sd()

## [1] 0.06328606

This is close to the standard standard error reported in the regression,

filter(tidy(mod), term == "ln_export_pop")

##            term   estimate  std.error statistic      p.value
## 1 ln_export_pop -0.7434568 0.06460478 -11.50777 1.551239e-30

The bootstrap is generating values assuming that the estimated model is the “true” model. But in order to calculate a p-value you need to generate the sampling distribution assuming that the null hypothesis is true.

beta_cor <- filter(results,
       term %in% c("ln_export_pop", "age", "age2", "male", "urban_dum") |
         str_detect(term, "education") |
         str_detect(term, "occupation") |
         str_detect(term, "living_conditions") |
         str_detect(term, "religion")) %>%
  select(term, estimate, .iter) %>%
  rename(iter = .iter) %>%
  spread(term, estimate, -iter) %>%
  select(-iter) %>%
  cor() %>%
  tidy() %>%
  rename(term = .rownames) %>%
  gather(term2, cor, -term)

ggplot(beta_cor, aes(x = term, y = term2, fill = cor)) +
  geom_raster() +
  scale_fill_viridis()

The \(\beta\) coefficients are not uncorrelated. We can also observe this in the estimated variance-covariance matrix from the regression,

cov2cor(vcov(mod)[1:6, 1:6])

##               (Intercept) ln_export_pop          age         age2
## (Intercept)    1.00000000  -0.296892726 -0.540483926  0.509289743
## ln_export_pop -0.29689273   1.000000000 -0.013298893  0.004992375
## age           -0.54048393  -0.013298893  1.000000000 -0.975296583
## age2           0.50928974   0.004992375 -0.975296583  1.000000000
## male          -0.07803927   0.029867740 -0.001162037 -0.025516600
## urban_dum     -0.05286134  -0.045585417  0.012435687 -0.018532057
##                       male   urban_dum
## (Intercept)   -0.078039269 -0.05286134
## ln_export_pop  0.029867740 -0.04558542
## age           -0.001162037  0.01243569
## age2          -0.025516600 -0.01853206
## male           1.000000000  0.01914562
## urban_dum      0.019145621  1.00000000

The confidence interval using the bootstrap standard errors is,

# OLS estimate
betahat <- coef(mod_bs)[["ln_export_pop"]]
z <- 1.96
filter(beta_bs, term == "ln_export_pop") %>%
  ungroup() %>%
  summarize(se = sd(estimate)) %>%
  mutate(conf.low = betahat - z * se,
         conf.high = betahat + z * se)

## Source: local data frame [1 x 3]
## 
##           se   conf.low  conf.high
##        (dbl)      (dbl)      (dbl)
## 1 0.04459353 -0.4833508 -0.3085442

The confidence interval using bootstrap quantiles is,

# OLS estimate
filter(beta_bs, term == "ln_export_pop") %>%
  ungroup() %>%
  summarize(conf.low = quantile(estimate, 0.025),
            conf.high = quantile(estimate, 0.975))

## Source: local data frame [1 x 2]
## 
##     conf.low  conf.high
##        (dbl)      (dbl)
## 1 -0.4772805 -0.3139828

In this case the two methods produce similar confidence intervals, since the sampling distribution of \(\hat\beta\) is symmetric. These confidence interval are also similar to those produced by OLS.

confint(mod_bs)

##                    2.5 %     97.5 %
## (Intercept)    1.7619152  1.7949429
## ln_export_pop -0.4825835 -0.3093115

This is somewhat disappointing. Bootstrapping in this manner does not help us produce confidence intervals accounting for the clustering that produces larger standard errors. The reason is that we are sampling from the “population”. In order to account for the correlation within clusters, we would need to account for the way that the sample was generated from the population. There are several ways to bootstrap with clustering. The most important part is to randomly sample ethnicities from the set of ethnicities.

See above.

The instructions given were too simplistic, and this requires more complicated code than I thought when writing it. The group_by argument for bootstrap() allows for sampling within the clusters, but we need to sample the clusters themselves.

ethnicities <- select(nunn, ethnicity) %>% unique()
iter <- 1024
bs_samples <- list()
for (i in seq_len(iter)) {
  # Resample ethnic groups
  ethnic_sample <- sample_frac(ethnicities, 1, replace = TRUE)
  # Merge with Nunn data to get a new dataset
  newdata <- left_join(ethnic_sample, nunn, by = "ethnicity")
  bs_samples[[i]] <- tidy(lm(trust_neighbors ~ ln_export_pop, data = newdata))
}
bs_samples <- bind_rows(bs_samples)

bs_samples %>%
  filter(term == "ln_export_pop")  %>%
  ungroup() %>%
  summarize(conf.low = quantile(estimate, 0.025),
            conf.high = quantile(estimate, 0.975),
            se = sd(estimate))

## Source: local data frame [1 x 3]
## 
##    conf.low conf.high        se
##       (dbl)     (dbl)     (dbl)
## 1 -0.807128 0.3850288 0.3199544

Now the confidence interval and standard errors are much larger than those reported by OLS, and the bootstrapping that did not account for clusters.

F-tests can only compare nested models on the same data. The models in Table 1 are not nested, because they are not subsets of each other. Each model uses a different measure of slave exports. The \(p\)-value of the test of a single coefficient and the equivalent \(F\)-test are the same. The \(F\)-statistic is the square of the p-value.

anova(mod_controls, mod_1_6)

## Analysis of Variance Table
## 
## Model 1: trust_neighbors ~ age + age2 + male + urban_dum + education + 
##     occupation + religion + living_conditions + district_ethnic_frac + 
##     frac_ethnicity_in_district + isocode
## Model 2: trust_neighbors ~ ln_export_pop + age + age2 + male + urban_dum + 
##     education + occupation + religion + living_conditions + district_ethnic_frac + 
##     frac_ethnicity_in_district + isocode
##   Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
## 1  17567 14888                                  
## 2  17566 14777  1     111.4 132.43 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tvalue <- tidy(mod_1_6) %>% filter(term == "ln_export_pop") %>%
  magrittr::extract2("statistic")
tvalue

## [1] -11.50777

tvalue ^ 2 

## [1] 132.4287