Week 4 Exercises

Complete solutions are available. Work each exercise before you open one.

Note: These exercises are tending toward “tedious” and “rote.” However, I want to make sure that we understand the implementation of the basic framework, so there’s a lot of tedious repetition in this homework.

Exercise 1 Berk (2010)

Read Berk (2010, pdf).

  1. Berk distinguishes Level I (descriptive), Level II (statistical inference), and Level III (causal inference) regression analyses. Summarize the features of each level in your own words. Explain why Berk argues that most regressions in criminology are Level I. Does his claim about criminology apply to political science as well?
  2. Select an empirical article you are already familiar with from your subfield that uses non-experimental (i.e., observational) data. Identify how the authors use regression. At what level does the authors claim to operate? Use direct quotes if possible. At what level does the analysis actually operate? Explain.
  3. Berk critiques observational regression treated as causal. He also lists common responses (e.g., “the assumptions are reasonable,” “we all have models,” “you have to make assumptions to make progress”). Pick one response and assess whether you find it persuasive. Defend your position. How might you design a study in your area that could move from Level I or II toward Level III without weak rhetorical defenses? Perhaps give an example.
  4. Berk concludes that, with rare exceptions, regression analyses of observational data should be treated as Level I. Do you agree? If you take Berk’s view seriously, what happens to publication standards, peer review, funding, and training?
Complete solution

No solution intended.

Exercise 2 Clark and Golder (2006), normal distribution

Clark and Golder (2006, pdf) fit the following regression model (that we’ve seen several times now).

  1. Use an R formula and model.matrix() to create the required design matrix X. Think carefully about the formula needed.
  2. Reproduce their estimates with the matrix multiplication \((X^\top X)^{-1} X^\top y\) or solve(t(X) %*% X) %*% t(X) %*% y.
  3. Reproduce their estimates using a normal linear model with optim(). You can easily adapt code like est_logit() from the notes, remembering that you must combine all parameters (i.e., the vector \(\beta\) and scalar \(\sigma\)) into a single vector \(\theta\) to give optim(). See the template below.
normal_ll <- function(par, y, X) {
  beta  <- par[1:ncol(X)]     # pull out the vector of betas
  sigma <- par[ncol(X) + 1]   # pull out the scalar sigma
  ... [and so on]
}

The least squares and ML estimates of the coefficients should match; for the normal linear model, the ML estimate is just the least squares solution.

# load Clark and Golder's data
cg <- crdata::cg2006  # from my data package
Hint 1

The table’s rows tell you which interactions the model includes. In an R formula, a*b adds a, b, and a:b. See Section 10.4 of the notes.

Hint 2

Section 11.0.1 of the notes computes solve(t(X)%*%X)%*%t(X)%*%y on real data. est_logit() in Section 11.1.2 is the optim() code to adapt. Notice that dnorm() needs a positive sd, so give \(\sigma\) a positive starting value.

Complete solution

Solution

# load data
cg <- crdata::cg2006
glimpse(cg)
Rows: 487
Columns: 8
$ country           <chr> "Argentina", "Argentina", "Argentina", "Argentina", …
$ year              <dbl> 1946, 1951, 1954, 1958, 1960, 1963, 1965, 1973, 1983…
$ average_magnitude <dbl> 10.53, 10.53, 4.56, 8.13, 4.17, 8.35, 4.17, 10.13, 1…
$ eneg              <dbl> 1.342102, 1.342102, 1.342102, 1.342102, 1.342102, 1.…
$ enep              <dbl> 5.750, 1.970, 1.930, 2.885, 5.485, 5.980, 5.155, 3.1…
$ upper_tier        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ en_pres           <dbl> 2.09, 1.96, 1.96, 2.65, 2.65, 3.90, 3.90, 2.66, 2.30…
$ proximity         <dbl> 1.00, 1.00, 0.20, 1.00, 0.20, 1.00, 0.33, 1.00, 1.00…
# create design matrix
f <- enep ~ eneg*log(average_magnitude) + eneg*upper_tier + en_pres*proximity
mf <- model.frame(f, data = cg)
X  <- model.matrix(f, data = mf)
y  <- model.response(mf)

# OLS
solve(t(X) %*% X) %*% t(X) %*% y
                                   [,1]
(Intercept)                  2.91570805
eneg                         0.11160359
log(average_magnitude)       0.07798753
upper_tier                  -0.05655491
en_pres                      0.26384754
proximity                   -3.09756561
eneg:log(average_magnitude)  0.26366122
eneg:upper_tier              0.05919037
en_pres:proximity            0.68317107
# MLE via optim (normal-errors)
normal_ll <- function(theta, y, X){
  k <- ncol(X)
  beta <- theta[1:k]
  sigma <- theta[k+1]
  mu <- X%*%beta
  sum(dnorm(y, mu, sigma, log = TRUE))
}

# use optim()
theta_start <- c(rep(0, ncol(X)), 2)
est <- optim(
  par     = theta_start,
  fn      = normal_ll,
  y       = y,
  X       = X,
  method  = "BFGS",
  control = list(fnscale = -1)
)
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced
est$par  # matches ols from above
 [1]  2.91570243  0.11161188  0.07798356 -0.05655386  0.26384097 -3.09755108
 [7]  0.26366069  0.05918970  0.68317320  1.31999212

Exercise 3 Matching formulas to linear predictors

The formulas below use the ZeligData::turnout dataset. race is a factor with levels others (the baseline) and white; age, educate, and income are numeric.

# load the turnout data
turnout <- ZeligData::turnout

The table below pairs fifteen linear predictors with the fifteen formulas that produce them, but formulas are scrambled.

Linear predictor \(\eta_i\) Formula
1 \(\beta_0 + \beta_1\text{age}_i\) A vote ~ age*educate
2 \(\beta_0 + \beta_1\log(\text{age}_i)\) B vote ~ race - 1
3 \(\beta_0 + \beta_1\sqrt{\text{income}_i}\) C vote ~ income + educate
4 \(\beta_0 + \beta_1\text{white}_i\), where \(\text{white}_i = 1\) for white respondents D vote ~ age:educate
5 a separate intercept for each race, and no common intercept E vote ~ sqrt(income)
6 \(\beta_0 + \beta_1(\text{income}_i + \text{educate}_i)\) F vote ~ race*educate
7 \(\beta_0 + \beta_1(\text{age}_i \times \text{educate}_i)\) G vote ~ age^2
8 \(\beta_0 + \beta_1\text{age}_i + \beta_2\text{educate}_i\) H vote ~ I(income + educate)
9 \(\beta_0 + \beta_1\text{income}_i + \beta_2\text{educate}_i\) I vote ~ race
10 \(\beta_0 + \beta_1\text{age}_i + \beta_2\text{age}_i^2\) J vote ~ (age + educate + income)^2
11 one intercept shared by both races, but a separate education slope for each K vote ~ age + I(age^2)
12 \(\beta_0 + \beta_1\text{age}_i + \beta_2\text{age}_i^2 + \beta_3\text{age}_i^3\) L vote ~ log(age)
13 \(\beta_0 + \beta_1\text{age}_i + \beta_2\text{educate}_i + \beta_3(\text{age}_i \times \text{educate}_i)\) M vote ~ race:educate
14 \(\beta_0 + \beta_1\text{white}_i + \beta_2\text{educate}_i + \beta_3(\text{white}_i \times \text{educate}_i)\) N vote ~ age + educate
15 all three of age, educate, and income, plus all three pairwise products O vote ~ poly(age, 3, raw = TRUE)
  1. Match each linear predictor to the one formula that produces it. Each formula is used exactly once. Work these out without running R.
  2. Now run model.matrix() on the ones you were unsure about and compare the column names to your answer. Note any you got wrong and explain what you expected instead.
  3. Several of these are worth a second look. For each group below, explain what R did and why.
    • G and J both use ^, and neither one produces a squared term. What does ^ do to the terms inside it, and why does age^2 collapse to a single age column? Which formula in the table is the only one that puts an actual square in the design matrix?
    • K and O both put a polynomial in age into the design matrix (one written out by hand and one built by poly()). What does raw = TRUE do? Refit O without it and compare the coefficients and the fitted values. Why does the fit not change?
    • D and M both use : without *. age:educate gives a product and no main effects, but race:educate gives two interaction columns where F gives one. Why does dropping the main effect change the number of interaction columns for the factor but not for the numeric variable?
    • C and H use the same two variables but produce a different number of columns. What restriction on C’s coefficients does H impose?
    • B and I describe the same two groups. What does each coefficient means? Confirm that the two models produce the same fitted values. Why can a model not include an intercept and an indicator for both races?
Hint 1

We need four operators:

  • + adds a column
  • : forms a product
  • * expands to both main effects and the product
  • - removes

Important: ^ is not arithmetic. It probably doesn’t work like you expect.

See Section 10.4 for * and :, Section 10.5.2 for ^ and I().

See Section 10.2 for factors and dummy variables.

Hint 2

For part 3, a categorical variable with \(m\) categories contributes \(m - 1\) indicators when the model has an intercept. See the “On the baseline” callout in Section 9.6 of the notes.

Complete solution

Solution

Part 1. 1-G, 2-L, 3-E, 4-I, 5-B, 6-H, 7-D, 8-N, 9-C, 10-K, 11-M, 12-O, 13-A, 14-F, 15-J

Part 2.

# load data
turnout <- ZeligData::turnout

# the fifteen formulas, labeled as in the table
formulas <- c(
  A = "vote ~ age*educate",
  B = "vote ~ race - 1",
  C = "vote ~ income + educate",
  D = "vote ~ age:educate",
  E = "vote ~ sqrt(income)",
  F = "vote ~ race*educate",
  G = "vote ~ age^2",
  H = "vote ~ I(income + educate)",
  I = "vote ~ race",
  J = "vote ~ (age + educate + income)^2",
  K = "vote ~ age + I(age^2)",
  L = "vote ~ log(age)",
  M = "vote ~ race:educate",
  N = "vote ~ age + educate",
  O = "vote ~ poly(age, 3, raw = TRUE)"
)

# print the column names of each design matrix
for (i in seq_along(formulas)) {
  X <- model.matrix(as.formula(formulas[i]), data = turnout)
  cat(names(formulas)[i], ": ", formulas[i], "\n",
      "   ", paste(colnames(X), collapse = " | "), "\n", sep = "")
}
A: vote ~ age*educate
   (Intercept) | age | educate | age:educate
B: vote ~ race - 1
   raceothers | racewhite
C: vote ~ income + educate
   (Intercept) | income | educate
D: vote ~ age:educate
   (Intercept) | age:educate
E: vote ~ sqrt(income)
   (Intercept) | sqrt(income)
F: vote ~ race*educate
   (Intercept) | racewhite | educate | racewhite:educate
G: vote ~ age^2
   (Intercept) | age
H: vote ~ I(income + educate)
   (Intercept) | I(income + educate)
I: vote ~ race
   (Intercept) | racewhite
J: vote ~ (age + educate + income)^2
   (Intercept) | age | educate | income | age:educate | age:income | educate:income
K: vote ~ age + I(age^2)
   (Intercept) | age | I(age^2)
L: vote ~ log(age)
   (Intercept) | log(age)
M: vote ~ race:educate
   (Intercept) | raceothers:educate | racewhite:educate
N: vote ~ age + educate
   (Intercept) | age | educate
O: vote ~ poly(age, 3, raw = TRUE)
   (Intercept) | poly(age, 3, raw = TRUE)1 | poly(age, 3, raw = TRUE)2 | poly(age, 3, raw = TRUE)3

The three that usually cause trouble:

  • G age^2 produces (Intercept) and age only, so it matches 1 rather than 10.
  • H I(income + educate) produces one coefficient on the sum; C income + educate produces two.
  • M race:educate produces raceothers:educate and racewhite:educate: a separate education slope for each race, with one shared intercept.

Part 3.

G and J. Inside a formula, ^2 means “interact all of these terms with each other, up to two at a time.” It doesn’t square anything. Interacting a variable with itself just gives that variable back, so G age^2 is just age. J (age + educate + income)^2 gives the three variables plus their three pairwise interactions. K is the only formula in the table that actually squares a variable. I() tells R to do the arithmetic rather than read ^ as a formula operator. See Section 10.5.2.

K and O. With raw = TRUE, poly() gives you the columns age, age^2, and age^3. Those are the columns you’d write by hand, so the coefficients are the same. Without raw = TRUE, R uses three new columns that are mixtures of those powers, built to be uncorrelated with each other. Any fit you can get with one set of columns you can get with the other, so the predicted probabilities are identical. Only the coefficients change, and they aren’t \(\beta_1\), \(\beta_2\), and \(\beta_3\) from the equation. See Section 10.3.1.

fit_raw  <- glm(vote ~ poly(age, 3, raw = TRUE), family = binomial, data = turnout)
fit_hand <- glm(vote ~ age + I(age^2) + I(age^3), family = binomial, data = turnout)
fit_orth <- glm(vote ~ poly(age, 3), family = binomial, data = turnout)

all.equal(unname(coef(fit_raw)), unname(coef(fit_hand)))  # same coefficients
[1] TRUE
rbind(raw = unname(coef(fit_raw)), orthogonal = unname(coef(fit_orth)))
                [,1]      [,2]          [,3]         [,4]
raw        -2.736302 0.1917304  -0.002614735 9.716534e-06
orthogonal  1.114008 8.3811586 -15.519527075 2.478206e+00
all.equal(fitted(fit_raw), fitted(fit_orth))              # same fit
[1] TRUE

D, M, and F. : adds the interaction. The two numeric variables have just one interaction between them, so D adds a single column. Howevver, R (re)codes a factor according to what else is in the model. M has no racewhite column, so R gives each race its own education slope and a shared intercept. F already has racewhite, so R adds only racewhite:educate, and that coefficient is the difference between the two slopes. See Section 10.4.

C and H. C gives income and education each its own coefficient. H forces them to share a coefficient (i.e., it adds the two variables together first, then estimates a single coefficient for the sum, so it estimates one fewer parameter). Note that the + inside I() is ordinary addition, not the formula operator.

B and I. These are the same model written two ways, in the sense that the predicted probabilities are identical. In I, the intercept is the log-odds of voting for others, and racewhite is how much higher whites are. In B, each coefficient is that group’s own log-odds. You can’t have the intercept and both indicators. The two indicators add up to a column of 1s, which is what the intercept already is, so this would create perfect collinearity. See Section 9.6.

fit_i <- glm(vote ~ race, family = binomial, data = turnout)
fit_b <- glm(vote ~ race - 1, family = binomial, data = turnout)

coef(fit_i)
(Intercept)   racewhite 
  0.5181383   0.6699191 
coef(fit_b)
raceothers  racewhite 
 0.5181383  1.1880573 
all.equal(fitted(fit_i), fitted(fit_b))
[1] TRUE

Exercise 4 Probit

The probit model uses the cdf of the normal distribution (e.g., often denoted as \(\Phi(\cdot)\); pnorm() in R) rather than the inverse logit.

  1. Adapt the code from the notes to fit a probit model (or see this Gist). Use the same ZeligData::turnout data and model specification vote ~ age + educate + income + race.
  2. Compare your probit coefficients to the logit coefficients from the notes. Notice that multiplying probit coefficients by 1.6 gives you roughly the logit coefficients.
  3. Adapt the code to compute the first difference and 90% confidence interval. Compare to the first difference and 90% confidence interval from the notes.
Hint

The only change from the logit is the inverse link: use pnorm() in place of plogis(), both in the log-likelihood and in the function that computes the first difference. You are adapting est_logit() in Section 11.1.2 and the first-difference code in Section 11.2.2 of the notes.

Complete solution
# load packages
library(tidyverse)
library(numDeriv)  

# bernoulli log-likelihood with probit inverse link
probit_ll <- function(beta, y, X) {
  p <- pnorm(X %*% beta)               # Φ(η)
  sum(dbinom(y, size = 1, prob = p, log = TRUE))
}

# function to fit model
est_probit <- function(f, data) {
  mf <- model.frame(f, data = data)
  X  <- model.matrix(f, data = mf)
  y  <- model.response(mf)

  beta_start <- rep(0, ncol(X))

  est <- optim(
    par     = beta_start,
    fn      = probit_ll,
    y       = y,
    X       = X,
    hessian = TRUE,
    method  = "BFGS",
    control = list(fnscale = -1)
  )

  if (est$convergence != 0) warning("Probit model did not converge.")

  list(beta_hat = est$par, var_hat = solve(-est$hessian))
}

# formula and data 
turnout <- ZeligData::turnout
f <- vote ~ age + educate + income + race

# fit model
fit_probit <- est_probit(f, data = turnout)

# compare probit coefficients to logit coefficients from notes
fit_probit$beta_hat
[1] -1.76361803  0.01652753  0.10425348  0.09632643  0.16270805
fit_probit$beta_hat*1.6 # similar to logit
[1] -2.82178884  0.02644405  0.16680557  0.15412229  0.26033288
# build chosen covariate rows
X_lo <- cbind(
  "constant" = 1,
  "age"      = quantile(turnout$age, probs = 0.25),
  "educate"  = median(turnout$educate),
  "income"   = median(turnout$income),
  "white"    = 1
)
X_hi <- X_lo
X_hi[, "age"] <- quantile(turnout$age, probs = 0.75)

# first difference (for probit)
fd_fn_probit <- function(beta, hi, lo) {
  pnorm(hi %*% beta) - pnorm(lo %*% beta)
}

# first difference and se
fd_hat_probit <- fd_fn_probit(fit_probit$beta_hat, X_hi, X_lo)
grad_probit <- grad(
  func = fd_fn_probit,
  x    = fit_probit$beta_hat,
  hi   = X_hi,
  lo   = X_lo
)
se_fd_hat_probit <- sqrt(grad_probit %*% fit_probit$var_hat %*% grad_probit)

# estimate
fd_hat_probit
         [,1]
25% 0.1421895
# 90% ci
fd_hat_probit - 1.64 * se_fd_hat_probit # lower
         [,1]
25% 0.1141743
fd_hat_probit + 1.64 * se_fd_hat_probit # upper
         [,1]
25% 0.1702048

Exercise 5 Russett and Oneal (2001)

This exercise uses data from Russett and Oneal (2001). The model the probability of a dispute in a dyad-year as a function of several predictors. A quick summary of their book is here. For our purposes, we’ll focus on dem.lo, which is the lower of the two Polity scores in the dyad. For more details on the variables, see ?crdata::ro2001 or Table 1 on p. 277 of Oneal and Russet (1997, pdf).

# load russett and oneal data from {crdata}
ro <- crdata::ro2001

Fit the logit model dispute ~ allies + lcaprat2 + contiguity + dem.lo + logdstab + power.

# reproduce their simplest logistic regression model p. 314, table A3.1
f <- dispute ~ allies + lcaprat2 + 
  contiguity + dem.lo + logdstab + power
fit <- glm(f, family = "binomial", data = ro)

# extract coefficient estimates
beta_hat <- coef(fit)

# extract covariance estimates
v_hat <- vcov(fit)

We can make a summary of coefficient estimates with modelsummary(), but it isn’t very informative about the quantity of interest. We want to know how the probability of a dispute is changing!

library(modelsummary)
modelsummary(fit, stars = TRUE, gof_map = NA)
Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.
(1)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) -0.468*
(0.186)
alliesNot Allies 0.694***
(0.063)
lcaprat2 -0.264***
(0.019)
contiguityNoncontiguous -0.980***
(0.073)
dem.lo -0.091***
(0.005)
logdstab -0.306***
(0.026)
powerMinor Powers -0.467***
(0.070)

Use beta_hat and v_hat from above with the invariance property and the delta method (via numDeriv::grad() and for loops) to compute the following quantities of interest and their 90% confidence interval.

  1. Plot the estimated probability (i.e., expected value) of a dispute as dem.lo varies from -10 to 10 (in small steps, of maybe one unit or 0.1 units). Fix all other numeric variables at their medians and alliesNot Allies, contiguityNoncontiguous, and powerMinor Powers at 0. Interpret. Compare the information in this plot to the information in the table above. (Aside: How can these probabilities vary nonlinearly if the dem.lo variable is included in the linear predictor linearly?)
  2. What is the change in the probability of a dispute (i.e., “first difference”; \(\hat{\pi}_{max} - \hat{\pi}_{min}\)) when dem.lo changes from its minimum to its maximum, fixing the other variables at the values described above? Notice that you can “kinda see” this value in the plot, but computing it directly gives you the exact value and the SE. The SE cannot be inferred from the plot.1
Hint 1

Section 11.2.2 of the notes loops over the rows of X_c and computes the expected value and its SE at each one, then computes a first difference the same way. That code is the template for both parts here.

Hint 2

For the aside: the linear predictor is linear in dem.lo, but the inverse logit is not. See Section 11.2.1 of the notes.

Complete solution

Solution

# load packages
library(tidyverse)

# load data
ro <- crdata::ro2001
glimpse(ro)
Rows: 39,996
Columns: 12
$ stateaname <chr> "United States", "United States", "United States", "United …
$ statebname <chr> "Canada", "Canada", "Canada", "Canada", "Canada", "Canada",…
$ statea     <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,…
$ stateb     <dbl> 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20,…
$ year       <dbl> 1920, 1921, 1922, 1923, 1924, 1925, 1926, 1927, 1928, 1929,…
$ dispute    <fct> No Dispute, No Dispute, No Dispute, No Dispute, No Dispute,…
$ allies     <fct> Not Allies, Not Allies, Not Allies, Not Allies, Not Allies,…
$ lcaprat2   <dbl> 3.241306, 3.171194, 3.247957, 3.180491, 3.191799, 3.190288,…
$ contiguity <fct> Contiguous, Contiguous, Contiguous, Contiguous, Contiguous,…
$ dem.lo     <int> 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, …
$ logdstab   <dbl> 5.820108, 5.820108, 5.820108, 5.820108, 5.820108, 5.820108,…
$ power      <fct> At Least One Great Power, At Least One Great Power, At Leas…
# reproduce their simplest logistic regression model p. 314, table A3.1
f <- dispute ~ allies + lcaprat2 + 
  contiguity + dem.lo + logdstab + power
fit <- glm(f, family = "binomial", data = ro)
arm::display(fit, digits = 4)
glm(formula = f, family = "binomial", data = ro)
                        coef.est coef.se
(Intercept)             -0.4684   0.1864
alliesNot Allies         0.6938   0.0626
lcaprat2                -0.2643   0.0186
contiguityNoncontiguous -0.9803   0.0730
dem.lo                  -0.0908   0.0049
logdstab                -0.3057   0.0260
powerMinor Powers       -0.4674   0.0703
---
  n = 39996, k = 7
  residual deviance = 13506.2, null deviance = 15166.4 (difference = 1660.2)
# extract coefficient estimates
beta_hat <- coef(fit)

# extract covariance estimates
v_hat <- vcov(fit)

# ---- helpers (copied from notes) ----
ev_fn <- function(beta, X) {
  plogis(X%*%beta)
}

fd_fn <- function(beta, hi, lo) {
  plogis(hi%*%beta) - plogis(lo%*%beta)
}


# ---- part 2 ----

# create chosen values for X
X_c <- cbind(
  "constant"   = 1, # intercept
  "allies"     = 0, 
  "lcaprat2"   = median(ro$lcaprat2, na.rm = TRUE),
  "contiguity" = 0,
  "dem.lo"     = -10:10, 
  "logdstab"   = median(ro$logdstab, na.rm = TRUE),
  "power"      = 0
)

# containers for estimated quantities of interest and ses
ev_hat <- numeric(nrow(X_c))
se_ev_hat <- numeric(nrow(X_c))

# loop over each row of X_c and compute qi and se
for (i in 1:nrow(X_c)) {
  # for the ith row of X...
  
  # invariance property
  ev_hat[i] <- ev_fn(beta_hat, X_c[i, ])
  
  # delta method
  grad <- grad(
    func = ev_fn, # what function are we taking the derivative of?
    x = beta_hat, # what variable(s) are we taking the derivative w.r.t.?
    X = X_c[i, ])  # what other values are needed?
  se_ev_hat[i] <- sqrt(grad %*% v_hat %*% grad)
}

# put X_c, qi estimates, and se estimates in data frame
qi <- cbind(X_c, ev_hat, se_ev_hat) |>
  data.frame() |>
  glimpse()
Rows: 21
Columns: 9
$ constant   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ allies     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ lcaprat2   <dbl> 2.830705, 2.830705, 2.830705, 2.830705, 2.830705, 2.830705,…
$ contiguity <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ dem.lo     <dbl> -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, …
$ logdstab   <dbl> 8.042244, 8.042244, 8.042244, 8.042244, 8.042244, 8.042244,…
$ power      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ ev_hat     <dbl> 0.05913318, 0.05427729, 0.04979905, 0.04567245, 0.04187273,…
$ se_ev_hat  <dbl> 0.004693969, 0.004266219, 0.003886446, 0.003549778, 0.00325…
# plot
ggplot(qi, aes(x = dem.lo, y = ev_hat, 
               ymin = ev_hat - 1.64*se_ev_hat, 
               ymax = ev_hat + 1.64*se_ev_hat)) + 
  geom_ribbon() + 
  geom_line() +
  labs(x = "dem.lo", y = "Pr(dispute)")

# ---- part 3 ----

# make X_lo
X_lo <- cbind(
  "constant"   = 1, # intercept
  "allies"     = 0, 
  "lcaprat2"   = median(ro$lcaprat2, na.rm = TRUE),
  "contiguity" = 0,
  "dem.lo"     = -10, 
  "logdstab"   = median(ro$logdstab, na.rm = TRUE),
  "power"      = 0
)

# make X_hi by modifying the relevant value of X_lo
X_hi <- X_lo
X_hi[, "dem.lo"] <- 10

# invariance property
fd_hat <- fd_fn(beta_hat, X_hi, X_lo)

# delta method
grad <- grad(
  func = fd_fn, 
  x = beta_hat, 
  hi = X_hi,
  lo = X_lo)  
se_fd_hat <- sqrt(grad %*% v_hat %*% grad)

# estimated fd
fd_hat
            [,1]
[1,] -0.04901941
# estimated se
se_fd_hat
            [,1]
[1,] 0.004235749
# 90% ci
fd_hat - 1.64*se_fd_hat  # lower
            [,1]
[1,] -0.05596604
fd_hat + 1.64*se_fd_hat  # upper
            [,1]
[1,] -0.04207278

Exercise 6 Holland (2015)

Holland (2015, pdf) writes:

My first hypothesis is that [the number of] enforcement operations [operations] drop off with the fraction of poor residents [lower] in an electoral district. So district poverty should be a negative and significant predictor of enforcement, but only in politically decentralized cities [Lima and Santiago]. Poverty should have no relationship with enforcement in politically centralized cities [Bogota] once one controls for the number of vendors (p. 362)”

We can use optim() to fit a Poisson regression model [Gist] and a negative binomial regression model to Holland’s data for Santiago. The logic of the Poisson regression model follows the logic of logistic regression: choose a distribution, choose an inverse link function, and fit the model with maximum likelihood.

The Poisson model. The outcome operations is a count (0, 1, 2, …). The Poisson distribution is a natural starting point for counts. It has one parameter \(\lambda\), which is both the mean and the variance. In Week 2, we fit a Poisson model without covariates, so every observation had the same \(\lambda\). Now we want \(\lambda_i\) to depend on covariates, just as \(\pi_i\) does in the logit model. But \(\lambda_i\) must be positive, while \(X_i\beta\) can be any real number. Thus we need an inverse link function that maps the real numbers onto the positive numbers. The usual choice is the exponential function, which makes the model

\[y_i \sim \text{Poisson}(\lambda_i) \text{, where } \lambda_i = \exp(X_i\beta).\]

Compare this to the logit model, where \(\pi_i = \operatorname{logit}^{-1}(X_i\beta)\). (Only the distribution and the inverse link function change.)

The negative binomial model. But the Poisson model makes a strong assumption. Because \(\lambda_i\) is both the mean and the variance, the model assumes \(\text{Var}(y_i) = E(y_i)\). Count data often have a variance much larger than the mean (we call this “overdispersion”). The negative binomial distribution relaxes this assumption. It has two parameters, a mean \(\mu\) and a dispersion parameter \(\theta\) (and its variance is \(\mu + \mu^2/\theta\)). The variance is always larger than the mean, and it grows as \(\theta\) shrinks. (As \(\theta \to \infty\), the negative binomial becomes the Poisson.) We model the mean exactly as in the Poisson model and treat \(\theta\) as a single constant to estimate. The model is

\[y_i \sim \text{NegBin}(\mu_i, \theta) \text{, where } \mu_i = \exp(X_i\beta).\]

Notice that \(\theta\) plays the same role as \(\sigma^2\) in the normal linear model. It does not depend on covariates, and it does not appear in the linear predictor. It’s a second parameter that describes the spread of the outcome, and we estimate it along with \(\beta\).

Model Distribution Modeled parameter Other parameters
Normal linear \(y_i \sim N(\mu_i, \sigma^2)\) \(\mu_i = X_i\beta\) \(\sigma^2\)
Logit \(y_i \sim \text{Bernoulli}(\pi_i)\) \(\pi_i = \operatorname{logit}^{-1}(X_i\beta)\) none
Poisson \(y_i \sim \text{Poisson}(\lambda_i)\) \(\lambda_i = \exp(X_i\beta)\) none
Negative binomial \(y_i \sim \text{NegBin}(\mu_i, \theta)\) \(\mu_i = \exp(X_i\beta)\) \(\theta\)

In R, glm(f, family = poisson) fits the Poisson model and MASS::glm.nb(f) fits the negative binomial model. Both use the log link by default, so both have \(\exp(X_i\beta)\) as the inverse link function.

  1. Fit the same model below using glm() and glm.nb(). Compare the coefficients and standard errors. Use {modelsummary} to make a side-by-side table of coefficient and SE estimates. Explain the differences in the estimates across models.
  2. Compute a relevant quantity of interest. Explain the differences in the estimates across models.
  3. Simulate the five fake data sets from the two fitted models (i.e., the predictive distribution). Compare the five data sets to the observed data. Note: Now that we have covariates, each observation has it’s own parameters, to make sure to that \(\lambda\)/lambda and \(\mu\)/mu are specific to each observation for the Poisson and NB models, respectively. You can do this by giving rpois() a vector lambda of length length(y) or nrow(X) rather than a scalar, for example. \(\theta\)/size is also needed to simulate fake data from the negative binomial distribution, but it’s constant.
# load data
holland <- crdata::holland2015 |>
  filter(city == "santiago")

# formula corresponds to model 1 for each city in holland (2015) table 2
f <- operations ~ lower + vendors + budget + population
Hint 1

Section 3.1 of the notes simulates five fake data sets from a Poisson model and compares them to the observed data. That code is the template for part 3.

Hint 2

predict(fit, type = "response") returns each observation’s \(\lambda\) or \(\mu\). glm.nb() stores \(\hat{\theta}\) in fit$theta.

Complete solution

Solution

# load packages
library(tidyverse)
library(numDeriv)
library(modelsummary)

# ---- part 1 ----
# (holland and f come from the chunk above)

# fit the poisson and negative binomial models
fit_pois <- glm(f, family = poisson, data = holland)
fit_nb <- MASS::glm.nb(f, data = holland)

# side-by-side table of coefficient and se estimates
modelsummary(list("Poisson" = fit_pois, "negative binomial" = fit_nb), 
             stars = TRUE, gof_map = NA)
Poisson negative binomial
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 2.619*** 2.706
(0.519) (1.850)
lower -0.038*** -0.046+
(0.009) (0.027)
vendors -0.221+ -0.188
(0.122) (0.225)
budget -0.001 -0.000
(0.000) (0.002)
population 0.015+ 0.021
(0.009) (0.027) 
# the negative binomial's extra parameter
fit_nb$theta
[1] 0.3675202
# the mean and variance of the observed counts
mean(holland$operations)
[1] 2.705882
var(holland$operations)
[1] 24.39572
  • The coefficient estimates are similar. Both models use the inverse link \(e^{X_i\beta}\), so both describe the average count the same way.
  • The standard errors are about two to four times larger for the negative binomial.
  • The Poisson model assumes \(\text{Var}(y_i) = E(y_i)\). The observed counts have a mean of about 2.7 and a variance of about 24, so the Poisson standard errors are much too small.
  • \(\hat{\theta} \approx 0.37\) is the extra parameter that lets the variance exceed the mean. Smaller values allow more extra variance.
# ---- part 2 ----

# function to compute the first difference
fd_fn <- function(beta, hi, lo) {
  exp(hi%*%beta) - exp(lo%*%beta)
}

# make X_lo
X_lo <- cbind(
  "constant"   = 1, # intercept
  "lower"      = quantile(holland$lower, probs = 0.25), 
  "vendors"    = median(holland$vendors),
  "budget"     = median(holland$budget),
  "population" = median(holland$population)
)

# make X_hi by modifying the relevant value of X_lo
X_hi <- X_lo
X_hi[, "lower"] <- quantile(holland$lower, probs = 0.75) 

# function to compute the first difference and its 90% ci for a fitted model
fd_ci <- function(fit, hi, lo) {
  
  # pull the estimates out of the fitted model
  beta_hat <- coef(fit)
  v_hat <- vcov(fit)
  
  # invariance property
  fd_hat <- fd_fn(beta_hat, hi, lo)
  
  # delta method
  grad_fd <- grad(
    func = fd_fn, 
    x = beta_hat, 
    hi = hi,
    lo = lo)  
  se_fd_hat <- sqrt(grad_fd %*% v_hat %*% grad_fd)
  
  # estimate and 90% ci
  c(fd_hat = fd_hat, 
    lwr = fd_hat - 1.64*se_fd_hat, 
    upr = fd_hat + 1.64*se_fd_hat)
}

# the change in expected operations as lower moves from its 25th to its 75th percentile
round(fd_ci(fit_pois, X_hi, X_lo), 2)
fd_hat    lwr    upr 
 -1.98  -2.86  -1.10 
round(fd_ci(fit_nb, X_hi, X_lo), 2)
fd_hat    lwr    upr 
 -2.13  -5.14   0.87
  • The estimates nearly agree at about two fewer operations in both models.
  • The negative binomial’s 90% confidence interval is about three times wider and includes zero.
  • Notice that the two models give about the same estimate, but different conclusions about how precise that estimate is.
# ---- part 3 ----

# each observation has its own lambda and mu, but theta is constant
lambda_hat <- predict(fit_pois, type = "response")
mu_hat <- predict(fit_nb, type = "response")
theta_hat <- fit_nb$theta

# create observed data set
observed_data <- tibble(operations = holland$operations, 
                        type = "observed")

# simulate five fake data sets from each model
set.seed(1234)
sim_list <- list()
for (i in 1:5) {
  sim_list[[i]] <- bind_rows(
    tibble(operations = rpois(nrow(holland), lambda = lambda_hat),
           type = paste0("Poisson #", i)),
    tibble(operations = rnbinom(nrow(holland), mu = mu_hat, size = theta_hat),
           type = paste0("negative binomial #", i))
  )
}

# combine the fake and observed data sets
gg_data <- bind_rows(sim_list) |>
  bind_rows(observed_data) |>
  glimpse()
Rows: 374
Columns: 2
$ operations <dbl> 0, 1, 1, 2, 4, 2, 0, 1, 4, 1, 1, 6, 8, 6, 0, 2, 2, 3, 5, 0,…
$ type       <chr> "Poisson #1", "Poisson #1", "Poisson #1", "Poisson #1", "Po…
# plot the observed and fake data sets
ggplot(gg_data, aes(x = operations)) + 
  geom_histogram() + 
  facet_wrap(vars(type))

# count the zeros and find the largest count in each data set
gg_data |>
  group_by(type) |>
  summarize(zeros = sum(operations == 0), 
            sd = sd(operations), 
            max = max(operations)) |>
  print(n = 11)
# A tibble: 11 × 4
   type                 zeros    sd   max
   <chr>                <int> <dbl> <dbl>
 1 Poisson #1               9  2.29     8
 2 Poisson #2               4  2.50     9
 3 Poisson #3               7  2.53     9
 4 Poisson #4               7  3.06    10
 5 Poisson #5               7  2.58     9
 6 negative binomial #1    18 23.7    127
 7 negative binomial #2    16  3.62    18
 8 negative binomial #3    16 10.1     32
 9 negative binomial #4    17  9.26    47
10 negative binomial #5    20  4.04    19
11 observed                18  4.94    16
  • 18 of the 34 districts have zero operations. The observed SD is about 4.9 and the largest count is 16.
  • The Poisson data sets have too few zeros (about four to nine) and too little spread (SD of about 2.3 to 3.1). None of them produces a count larger than 10.
  • The negative binomial data sets get the number of zeros about right (16 to 20) and do produce large counts. But the largest counts are sometimes much larger than anything observed, such as the 127 in the first data set.

Exercise 7 Clark and Golder (2006); t distribution

The code below reproduces the coefficient estimates for the normal linear model from Clark and Golder (2006) as shown in an earlier exercise. Re-estimate their model using a location-scale t distribution (see Exercise 11 for Week 2) modeling the mean \(\mu = X\beta\). Use ML to estimate the coefficients \(\beta\), the scale \(\sigma\), and the df \(\nu\). Comment on the differences between the coefficient estimates and standard errors using the normal model and the t. What does the ML estimate \(\hat{\nu}\) of the degrees of freedom parameter \(\nu\) tell us about Clark and Golder’s data?

# load packages
library(tinytable)

# load data
cg <- crdata::cg2006

# create design matrix
f <- enep ~ eneg*log(average_magnitude) + eneg*upper_tier + en_pres*proximity

# normal model
fit <- glm(f, data = cg, family = gaussian)

For convenience, here is how me might make a quick table using tinytable::tt().

# make a table
data.frame(term = names(coef(fit)),
           beta_hat = coef(fit),
           se = sqrt(diag(vcov(fit)))) |>
  tt(digits = 2, escape = TRUE) |>
  theme_latex(placement = "H")
term beta_hat se
(Intercept) 2.916 0.176
eneg 0.112 0.071
log(average_magnitude) 0.078 0.116
upper_tier -0.057 0.02
en_pres 0.264 0.064
proximity -3.098 0.352
eneg:log(average_magnitude) 0.264 0.067
eneg:upper_tier 0.059 0.014
en_pres:proximity 0.683 0.137
Hint 1

The log-likelihood has the same shape as the one you wrote for Exercise 2, with \(\nu\) as a third parameter in par.

Hint 2
  • The R function metRology::dt.scaled(..., log = TRUE) computes the required log-likelihood for each observation.
  • You need to estimate the \(\sigma =\)sd and \(\nu =\) df parameters, but they are fixed.
  • Model the mean as a with the usual linear predictor \(X_i\beta\). No inverse link function is needed, because mean is unbounded.
  • When creating the log-likelihood function, the first argument must be a single parameter vector that include all parameters stacked into a single vector. I often call this par. Then extract par into beta, sigma, and nu immediately inside the function.
Complete solution

Solution.

# assuming code from question already included

# create design matrix
mf <- model.frame(f, data = cg)
X  <- model.matrix(f, data = mf)
y  <- model.response(mf)

# t log-likelihood
t_ll <- function(par, y, X){
  
  # pull out parameters from par for easy reading
  k <- ncol(X)
  beta <- par[1:k]
  sigma <- par[k+1]  # 2nd to last
  nu <- par[k+2] # last
  
  #
  mu <- X%*%beta
  sum(metRology::dt.scaled(y, 
                           mean = mu, 
                           sd = sigma, 
                           df = nu, 
                           log = TRUE))
}

# use optim()
par_start <- c(rep(0, ncol(X)), 2, 10)
est <- optim(
  par     = par_start,
  fn      = t_ll,
  y       = y,
  X       = X,
  method  = "BFGS",
  control = list(fnscale = -1), 
  hessian = TRUE
)
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in stats::dt((x - mean)/sd, df, ncp = ncp, log = TRUE): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in stats::dt((x - mean)/sd, df, ncp = ncp, log = TRUE): NaNs produced
Warning in stats::dt((x - mean)/sd, df, ncp = ncp, log = TRUE): NaNs produced
Warning in log(sd): NaNs produced
Warning in stats::dt((x - mean)/sd, df, ncp = ncp, log = TRUE): NaNs produced
Warning in stats::dt((x - mean)/sd, df, ncp = ncp, log = TRUE): NaNs produced
Warning in stats::dt((x - mean)/sd, df, ncp = ncp, log = TRUE): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in log(sd): NaNs produced
Warning in stats::dt((x - mean)/sd, df, ncp = ncp, log = TRUE): NaNs produced
Warning in stats::dt((x - mean)/sd, df, ncp = ncp, log = TRUE): NaNs produced
Warning in log(sd): NaNs produced
data.frame("term" = c(colnames(X), "sigma", "nu"),
      "beta_hat" = est$par, 
      "se" = sqrt(diag(solve(-est$hessian))))  |>
  tt(digits = 2, escape = TRUE) |>
  theme_latex(placement = "H")
term beta_hat se
(Intercept) 2.535 0.16
eneg 0.211 0.083
log(average_magnitude) 0.4 0.095
upper_tier -0.046 0.017
en_pres 0.193 0.059
proximity -2.92 0.229
eneg:log(average_magnitude) 0.023 0.057
eneg:upper_tier 0.05 0.013
en_pres:proximity 0.707 0.114
sigma 0.715 0.05
nu 2.243 0.314

The estimate of the degrees of freedom parameter is 2.3 with a standard error of about 0.3. This means that the distribution of the effective number of electoral parties deviates substantially from normal (when df > 10, the t starts to look somewhat normal). \(\nu \approx 2.5\) suggests that observed values often fall several standard deviations away from the mean. This has a nice substantive interpretation. While there seems to be some predictability to ENEP, some cases are highly unusual, which suggests “something different” might be going on in those cases.

Exercise 8 Connections

We have covered a wide range of tools. Use the list below as a menu and sketch the connections between two or more ideas. Represent the connections in any format you prefer; a diagram with arrows might make sense to you or you might prefer short written explanations. Simply noting that two things are connected is insufficient—most things are related in this course. Focus on why ideas are connected or how they fit together.

Hint: Definitely describe the connections between maximum likelihood estimates, Fisher information, the invariance property, and the delta method. Try for 2-3 other collections of 3 to 8 things that fit together as well.

  • Mathematical Tools: log, first derivative, second derivative, gradient, Hessian, matrix multiplication, pdf/pmf, expected value, variance
  • Engines: OLS, maximum likelihood, method of moments; also likelihood function, log-likelihood function
  • Distributions: Bernoulli, exponential, Weibull, log-normal, Poisson, negative binomial, t, normal, beta, uniform
  • Confidence intervals: parametric bootstrap, Fisher information, delta method, Wald method
  • Hypothesis Tests: nothing so far
  • Quantities of Interest: invariance property, odds, exponential rate/mean, beta shape parameters to mean and variance, expected value, first difference
  • Estimate Evaluation: predictive distribution
  • Estimator Evaluation: sampling distribution, bias, consistency, coverage, asymptotic distribution of estimator, Monte Carlo simulation
  • R: optim(), numDeriv()::grad(), %*%, density + distribution + quantile + random number generator functions for various distributions
  • Data Sets: toothpaste cap, ZeligData::turnout, Clark and Golder’s (2006) crdata::cg2006, Lahman’s batting_average, WDI, Herron’s hockey, Holland’s (2015) crdata::holland2015, the UF Election Lab turnout rates
Complete solution

No solution intended.

References

Berk, Richard. 2010. “What You Can and Can’t Properly Do with Regression.” Journal of Quantitative Criminology 26 (4): 481–87. https://doi.org/10.1007/s10940-010-9116-4.
Clark, William Roberts, and Matt Golder. 2006. “Rehabilitating Duverger’s Theory.” Comparative Political Studies 39 (6): 679–708. https://doi.org/10.1177/0010414005278420.
Holland, Alisha C. 2015. “The Distributive Politics of Enforcement.” American Journal of Political Science 59 (2): 357–71. https://doi.org/10.1111/ajps.12125.
Oneal, John R., and Bruce M. Russet. 1997. “The Classical Liberals Were Right: Democracy, Interdependence, and Conflict, 1950-1985.” International Studies Quarterly 41 (2): 267–94. https://doi.org/10.1111/1468-2478.00042.
Russett, Bruce, and John R. Oneal. 2001. Triangulating Peace: Democracy, Interdependence, and International Organizations. New York: W. W. Norton & Company.

Footnotes

  1. I usually recommend computing a first difference as a change from the 25th to the 75th percentile, rather than the min-to the max. The min-max comparison can be sensitive to model specification and represent usual, unrealistic scenarios.↩︎