In mneunhoe/simloglm: Simulate Correct Quantities of Interest from Linear Models with Logged Dependent Variables

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(simloglm)

df <- simloglm::hollibaugh_rothenberg_2018

# Data frame with the full complement of control variables.
ols.frame <- na.omit(subset(
  df,
  select = c(
    logdist,
    signedagency,
    sameside,
    dim1,
    dim2,
    dim1_alt,
    dim2_alt,
    dim1_nofixed,
    dim2_nofixed,
    pres_fil_dist,
    term_timediffyear,
    toptier,
    midtier,
    first.hundred,
    first.term,
    pres.approval,
    pres,
    nom_fil_dist
  )
))


m4 <- lm(
  logdist ~ signedagency * sameside + dim1 + dim2 +
    term_timediffyear + toptier + midtier +
    pres_fil_dist + first.hundred +
    first.term + pres.approval ,
  data = ols.frame
)



n_range <- 101

sd_range <- seq(-2, 2, length.out = n_range)

range_ind <-
  seq(
    min(df$dim1, na.rm = TRUE),
    max(df$dim1, na.rm = TRUE),
    length.out = n_range
  )

range_div <-
  seq(
    min(df$pres_fil_dist, na.rm = TRUE),
    max(df$pres_fil_dist, na.rm = TRUE),
    length.out = n_range
  )

# First we calculate the QoI with the original scenario of the authors

results_ind_aca <-
  simloglm(
    m4,
    scenario = list(
      signedagency = mean(m4$model$signedagency),
      sameside = mean(m4$model$sameside),
      dim1 = range_ind,
      dim2 = mean(m4$model$dim2),
      term_timediffyear = mean(m4$model$term_timediffyear),
      toptier = 0,
      midtier = 0,
      pres_fil_dist = mean(m4$model$pres_fil_dist),
      first.hundred = 0,
      first.term = 0,
      pres.approval = mean(m4$model$pres.approval)
    )
  )

# Second we choose an arbitrary extreme scenario.

results_ind_aca_extreme <-
  simloglm(
    m4,
    scenario = list(
      signedagency = min(m4$model$signedagency),
      sameside = max(m4$model$sameside),
      dim1 = range_ind,
      dim2 = min(m4$model$dim2),
      term_timediffyear = max(m4$model$term_timediffyear),
      toptier = 0,
      midtier = 0,
      pres_fil_dist = max(m4$model$pres_fil_dist),
      first.hundred = 0,
      first.term = 1,
      pres.approval = max(m4$model$pres.approval)
    )
  )

cols <- viridis::viridis(2)

m <- cbind(1, 2)
layout(m,
       widths = c(1, 1))

summary_mean <- get_summary(results_ind_aca, which_qoi = "mean")
summary_median <- get_summary(results_ind_aca, which_qoi = "median")


fd_extreme <- get_first_difference(results_ind_aca_extreme, which_qoi = "median", which_scenarios = c(1, 101))
fd <- get_first_difference(results_ind_aca, which_qoi = "median", which_scenarios = c(1, 101))


par(mfrow = c(1,2))

# Plot 1: Agency Decision-Maker Independence
par(mar = c(5,     # bottom (5)
            9.7,   # left   (4)
            2,     # top    (4)
            0.5),    # right  (2)),
    xpd = T)


hist(
  df$dim1,
  breaks = 20,
  xaxt = "n",
  yaxt = "n",
  main = "",
  ylab = "",
  xlab = "",
  xlim = c(-1, 2.5),
  ylim = c(0, 10000),
  border = adjustcolor("grey50", alpha = 0),
  col = adjustcolor("grey50", alpha = 0.5)
)
par(new = TRUE)  # new plot
plot(
  range_ind,
  summary_median$point_estimate,
  type = "n",
  ylim = c(0.1, 2),
  main = "",
  xlab = "",
  ylab = "",
  yaxt = "n",
  xaxt = "n",
  xlim = c(-1, 2.5),
  bty = "n"
)
box(lty = "solid", col = "grey50")
segments(y0 = seq(0.1, 2, 0.1), x0 = -1.14, x1 = 2.64,  lty = "dashed", col = "grey90")

axis(
  2,
  at = seq(0.1, 2, 0.1),
  las = 2,
  col.axis = "grey50",
  col = "grey50"
)
axis(1,
     at = seq(-1, 2.5, 0.5),
     col.axis = "grey50",
     col = "grey50")
mtext(2,
      text = "Nominee-President Ideological Divergence",
      line = 4.1)
mtext(3,
      text = expression(bold("Agency Decision-Maker Independence")),
      line = 0.4)
mtext(1,
      text = "Value of Independent Variable",
      line = 3.1)
text(x = 2.4,
     y = 0.25,
     labels = "Conditional\nMedian",
     srt = 90,
     cex = 0.8)
text(x = 2.4,
     y = 0.7,
     labels = "Conditional\nMean",
     srt = 90,
     cex = 0.8)


# Mean
polygon(
  c(rev(range_ind), range_ind),
  c(rev(summary_mean$quantiles[1, ]), summary_mean$quantiles[2, ]),
  col = adjustcolor(cols[2], 0.6),
  border = NA
)
lines(range_ind, summary_mean$point_estimate)

# Median

polygon(
  c(rev(range_ind), range_ind),
  c(rev(summary_median$quantiles[1, ]), summary_median$quantiles[2, ]),
  col = adjustcolor(cols[1], 0.6),
  border = NA
)
lines(range_ind, summary_median$point_estimate)


# Plot 2: President-Filibuster Divergence
par(mar = c(5,       # bottom (5)
            0.5,     # left   (4)
            2,       # top    (4)
            9.7),    # right  (2)
    xpd = T)



plot(fd$fd, 1, ylim = c(0.5, 2.5), xlim = c(-0.25, 0.8), yaxt = "n", xaxt = "n", xlab = "", bty = "n", ylab = "", type = "n")
box(lty = "solid", col = "grey50")

axis(1,
     at = seq(0, 0.5, 0.5),
     col.axis = "grey50",
     col = "grey50")
mtext(1,
      text = "First Difference of Medians",
      line = 3.1)
mtext(3,
      text = expression(bold("First Differences")),
      line = 0.4)
segments(x0 = seq(-0.5, 0.75, 0.25), y0 = 0.43, y1 = 2.57,  lty = "dashed", col = "grey90")
segments(x0 = 0, y0 = 0.43, y1 = 2.57,  lty = "dashed", col = "black")
segments(x0 = fd$quantiles[1], x1 = fd$quantiles[2], y0 = 1, lwd = 6, lend = 1, col = viridis::viridis(2)[1])
points(x = fd$fd, y = 1, cex = 2, pch = 124, col = "white")
segments(x0 = fd_extreme$quantiles[1], x1 = fd_extreme$quantiles[2], y0 = 2, lwd = 6, lend = 1, col = viridis::viridis(2)[1])
points(fd_extreme$fd, 2, cex = 2, pch = 124, col = "white")
text(x = fd$quantiles[2], y = 1, labels = "Average Case Scenario", pos = 4)
text(x = fd_extreme$quantiles[2], y = 2, labels = "Alternative Scenario", pos = 4)

mneunhoe/simloglm documentation built on June 14, 2022, 5:18 p.m.