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

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

library(simloglm)

# Clustered Standard Errors
# Function according to Arai:

clx <- function(fm, dfcw, cluster) {
  require(sandwich)
  require(lmtest)
  require(zoo)

  M <- length(unique(cluster)) # Number of Clusters
  N <- length(cluster) # Number of Observations

  dfc <-
    (M / (M - 1)) * ((N - 1) / (N - fm$rank)) # Adjust degrees of freedom

  u <- apply(sandwich::estfun(fm), 2, function(x)
    tapply(x, cluster, sum))

  vcovCL <- dfc * sandwich::sandwich (fm, meat = crossprod(u) / N) * dfcw

  #print(lmtest::coeftest(fm, vcovCL)) # Print Results to console

  return(vcovCL) # Return Variance-Covariance Matrix
}

df <- simloglm::shepherd_you_2020

df$ln_meansalary <- log(df$meansalary + 1)
df$nolobstaff <- df$numstaff - df$futurelob

df$femaleratio <- df$numfemale / df$numstaff


df$lnles <- log(df$les + 1)
df$lntotbill <- log(df$totbill + 1)
df$lnssbill <- log(df$ss_bills + 1)

# Define the Model Formulas

member_char1 <-
  c(
    "majority dwnom1 budget chair subchr seniority maj_leader min_leader power dem membecamelob female afam latino state_leg south_dem"
  )
member_char1 <- gsub(" ", " + ", member_char1)

member_char2 <-
  c("majority dwnom1 budget chair subchr seniority maj_leader min_leader power")
member_char2 <- gsub(" ", " + ", member_char2)

staff_char1 <-
  c("ln_meansalary nolobstaff futurelob cstaff_lob femaleratio")
staff_char1 <- gsub(" ", " + ", staff_char1)

ff <-
  paste0("lnles ~ ",
         staff_char1,
         " + ",
         member_char2,
         " + congress + icpsr")

df$congress <- as.factor(df$congress)
df$icpsr <- as.factor(df$icpsr)
# Run Model 4

m4 <- lm(ff, data = df)

vcovCL <-
  clx(m4, 1, m4$model$'icpsr') # Clustered Variance-Covariance Matrix

# Residualize Treatment for interpretation according to

est <- m4$model

## make an object that is the variable name of the treatment

x <- "futurelob"

## store the coefficient estimate

b <- m4$coefficients[x]
x.resid <-
  lm(
    futurelob ~ ln_meansalary + nolobstaff  + cstaff_lob + femaleratio + majority + dwnom1 + budget + chair + subchr + seniority + maj_leader + min_leader + power +congress + icpsr,
    data = est
  )$residuals

reg2 <- lm(est$lnles ~ x.resid)


y.resid <-
  lm(
    lnles ~ ln_meansalary + nolobstaff  + cstaff_lob + femaleratio + majority + dwnom1 + budget + chair + subchr + seniority + maj_leader + min_leader + power + congress + icpsr,
    data = est
  )$residuals

reg3 <- lm(y.resid ~ x.resid)

# Range of original treatment variable (futurelob)
range(est[, x])

# SD of original treatment variable (futurelob)
(sdx <- sd(est[, x]))

# SD of residualized Treatment
(sdxresid <- sd(x.resid))
round(sdxresid, 2)

# Set scenario for simulation

# Extract Model Frame
mf <- model.frame(ff, m4$model)

# Transform to Model Frame to Design Matrix
mm <- model.matrix(as.formula(ff), mf)

# Average Case Approach, Baseline all set to its mean
scenario <- apply(mm, 2, mean)

# Average Case Approach, 1 SD increase in futurelob
scenario_sd <- scenario
scenario_sd["futurelob"] <- scenario_sd["futurelob"] + sdxresid

# Combine the two scenarios in one matrix
scenario <- rbind(scenario, scenario_sd)

set.seed(20220608)
results_aca <-
  simloglm(
    m4,
    nsim_est = 10,
    X = scenario,
    fast = T
  )

Results with average case approach

# Summarize results
# First calculate the offset to fully retransform the values
results_aca_retransformed <- lapply(results_aca, function(x){ if(!("list" %in% class(x))) {
x - 1
}
})
class(results_aca_retransformed) <- "simloglm"
# First Difference as difference between expected values on the original scale
fd_mean_aca_max <- get_first_difference(results_aca_retransformed, which_qoi = "mean", alpha = 0)
fd_mean_aca_99 <- get_first_difference(results_aca_retransformed, which_qoi = "mean", alpha = 0.01)
fd_mean_aca_95 <- get_first_difference(results_aca_retransformed, which_qoi = "mean", alpha = 0.05)
fd_mean_aca_80 <- get_first_difference(results_aca_retransformed, which_qoi = "mean", alpha = 0.2)
fd_median_aca_max <- get_first_difference(results_aca_retransformed, which_qoi = "median", alpha = 0)
fd_median_aca_99 <- get_first_difference(results_aca_retransformed, which_qoi = "median", alpha = 0.01)
fd_median_aca_95 <- get_first_difference(results_aca_retransformed, which_qoi = "median", alpha = 0.05)
fd_median_aca_80 <- get_first_difference(results_aca_retransformed, which_qoi = "median", alpha = 0.2)
# Ratio instead of fd
# The minus 1 offset matters here (doesn't for the fd above)
percent_increase_aca_max <- get_ratio(results_aca_retransformed, which_qoi = "median", alpha = 0)
percent_increase_aca_99 <- get_ratio(results_aca_retransformed, which_qoi = "median", alpha = 0.01)
percent_increase_aca_95 <- get_ratio(results_aca_retransformed, which_qoi = "median", alpha = 0.05)
percent_increase_aca_80 <- get_ratio(results_aca_retransformed, which_qoi = "median", alpha = 0.2)
par(
mar = c(4.1, 3.1, 0.1, 0.1),
lend = 1,
las = 1,
mfrow = c(3, 1)
)
plot(
1,
1,
ylim = c(0, 2),
pch = "|",
col = adjustcolor(viridis::viridis(2)[1], 0.2),
type = "n",
bty = "n",
ylab = "",
yaxt = "n",
xaxt = "n",
xlab = "",
xlim = pretty(c(min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)))[c(1, length(pretty(c(
min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)
))))]
)
title("",
xlab = "First Difference (Mean) in Legislative Effectiveness Scores (LES)",
ylab = "",
col.lab = "grey")
axis(
1,
at = c(min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)),
labels = round(c(min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)), 2),
lwd = 0,
lwd.ticks = 1,
col = "grey",
col.axis = "grey"
)
axis(
1,
at = pretty(c(min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)))[c(2:(length(pretty(c(
min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)
))) - 1))],
lwd = 0.5,
lwd.ticks = 1,
col = "grey",
col.axis = "grey"
)
axis(
2,
at = 1,
labels = "(A)",
lwd = 1,
lwd.ticks = 0,
col = "grey",
col.axis = "grey"
)
segments(
x0 = fd_mean_aca_80$quantiles[1,],
y0 = 1,
x1 = fd_mean_aca_80$quantiles[2,],
col = adjustcolor(viridis::viridis(2)[1], 1),
lwd = 8,
lend = 1
)
segments(
x0 = fd_mean_aca_95$quantiles[1,],
y0 = 1,
x1 = fd_mean_aca_95$quantiles[2,],
col = adjustcolor(viridis::viridis(2)[1], 1),
lwd = 4,
lend = 1
)
segments(
x0 = fd_mean_aca_99$quantiles[1,],
y0 = 1,
x1 = fd_mean_aca_99$quantiles[2,],
col = adjustcolor(viridis::viridis(2)[1], 1),
lwd = 2,
lend = 1
)
points(fd_mean_aca_99$fd, 1, pch = 124, col = viridis::viridis(2)[2])
legend(
"top",
legend = c("Point Estimate", "80% CI", "95% CI", "99% CI"),
pch = c(124, NA, NA, NA),
lty = c(NA, "solid", "solid", "solid"),
lwd = c(NA, 8, 4, 2),
col = c(
viridis::viridis(2)[2],
viridis::viridis(2)[1],
viridis::viridis(2)[1],
viridis::viridis(2)[1]
),
bty = "n",
ncol = 5,
text.col = "grey"
)
plot(
1,
1,
ylim = c(0, 2),
pch = "|",
col = adjustcolor(viridis::viridis(2)[1], 0.2),
type = "n",
bty = "n",
ylab = "",
yaxt = "n",
xaxt = "n",
xlab = "",
xlim = pretty(c(min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)))[c(1, length(pretty(c(
min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)
))))]
)
title("",
xlab = "First Difference (Median) in Legislative Effectiveness Scores (LES)",
ylab = "",
col.lab = "grey")
axis(
1,
at = c(min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)),
labels = round(c(min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)), 2),
lwd = 0,
lwd.ticks = 1,
col = "grey",
col.axis = "grey"
)
axis(
1,
at = pretty(c(min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)))[c(2:(length(pretty(c(
min(fd_mean_aca_max$quantiles), max(fd_mean_aca_max$quantiles)
))) - 1))],
lwd = 0.5,
lwd.ticks = 1,
col = "grey",
col.axis = "grey"
)
axis(
2,
at = 1,
labels = "(B)",
lwd = 1,
lwd.ticks = 0,
col = "grey",
col.axis = "grey"
)
segments(
x0 = fd_median_aca_80$quantiles[1,],
y0 = 1,
x1 = fd_median_aca_80$quantiles[2,],
col = adjustcolor(viridis::viridis(2)[1], 1),
lwd = 8,
lend = 1
)
segments(
x0 = fd_median_aca_95$quantiles[1,],
y0 = 1,
x1 = fd_median_aca_95$quantiles[2,],
col = adjustcolor(viridis::viridis(2)[1], 1),
lwd = 4,
lend = 1
)
segments(
x0 = fd_median_aca_99$quantiles[1,],
y0 = 1,
x1 = fd_median_aca_99$quantiles[2,],
col = adjustcolor(viridis::viridis(2)[1], 1),
lwd = 2,
lend = 1
)
points(fd_median_aca_99$fd, 1, pch = 124, col = viridis::viridis(2)[2])
legend(
"top",
legend = c("Point Estimate", "80% CI", "95% CI", "99% CI"),
pch = c(124, NA, NA, NA),
lty = c(NA, "solid", "solid", "solid"),
lwd = c(NA, 8, 4, 2),
col = c(
viridis::viridis(2)[2],
viridis::viridis(2)[1],
viridis::viridis(2)[1],
viridis::viridis(2)[1]
),
bty = "n",
ncol = 5,
text.col = "grey"
)
plot(
1,
1,
ylim = c(0, 2),
pch = "|",
col = adjustcolor(viridis::viridis(2)[1], 0.2),
type = "n",
bty = "n",
ylab = "",
yaxt = "n",
xaxt = "n",
xlab = "",
xlim = pretty(c(min(percent_increase_aca_max$quantiles), max(percent_increase_aca_max$quantiles))-1)[c(1, length(pretty(c(min(
percent_increase_aca_max$quantiles
), max(
percent_increase_aca_max$quantiles
))-1)))]
)
title("",
xlab = "% Change Legislative Effectiveness Scores (LES)",
ylab = "",
col.lab = "grey")
axis(
1,
at = c(min(percent_increase_aca_max$quantiles), max(percent_increase_aca_max$quantiles))-1,
labels = round(c(min(percent_increase_aca_max$quantiles), max(percent_increase_aca_max$quantiles))-1, 2) * 100,
lwd = 0,
lwd.ticks = 1,
col = "grey",
col.axis = "grey"
)
axis(
1,
at = pretty(c(min(percent_increase_aca_max$quantiles), max(percent_increase_aca_max$quantiles))-1)[c(2:(length(pretty(c(
min(percent_increase_aca_max$quantiles), max(percent_increase_aca_max$quantiles)
)-1)) -1))],
labels =  pretty(c(min(percent_increase_aca_max$quantiles), max(percent_increase_aca_max$quantiles))-1)[c(2:(length(pretty(c(
min(percent_increase_aca_max$quantiles), max(percent_increase_aca_max$quantiles)
)-1))-1 ))] * 100,
lwd = 0.5,
lwd.ticks = 1,
col = "grey",
col.axis = "grey"
)
axis(
2,
at = 1,
labels = "(C)",
lwd = 1,
lwd.ticks = 0,
col = "grey",
col.axis = "grey"
)
segments(
x0 = percent_increase_aca_80$quantiles[1,]-1,
y0 = 1,
x1 = percent_increase_aca_80$quantiles[2,]-1,
col = adjustcolor(viridis::viridis(2)[1], 1),
lwd = 8,
lend = 1
)
segments(
x0 = percent_increase_aca_95$quantiles[1,]-1,
y0 = 1,
x1 = percent_increase_aca_95$quantiles[2,]-1,
col = adjustcolor(viridis::viridis(2)[1], 1),
lwd = 4,
lend = 1
)
segments(
x0 = percent_increase_aca_99$quantiles[1,]-1,
y0 = 1,
x1 = percent_increase_aca_99$quantiles[2,]-1,
col = adjustcolor(viridis::viridis(2)[1], 1),
lwd = 2,
lend = 1
)
points(percent_increase_aca_99$ratio-1, 1, pch = 124, col = viridis::viridis(2)[2])

# Average case result
percent_increase_aca_95

# One sided test
mean(results_aca_retransformed$median[,2,1]/results_aca_retransformed$median[,1,1] -1  > 0)

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