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

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

library(simloglm)

# set seed
set.seed(220608)

# set parameters
n <- 10
n_mc_sims <- 2000  # monte carlo simulations to obtain the sampling distribution Set to 10000 for results in paper

X_c <- 20
b_cons <- 2.5
b_educ <- 0.1
sigma <- 1.0  # should update this so it's robust to mixing up var and sd
educ <- seq(from = 10, to = 20, length.out = n)


# do simulation
sim_list <- vector("list", length = n_mc_sims)
#pb <- txtProgressBar(min = 0, max = n_sims, style = 3)
for (i in 1:n_mc_sims) {
  if (i %% 10 == 0) print(i)
  # find maximum likelihood estimate of quantity of interest
  e <- rnorm(n, mean = 0, sd = sigma)
  inc <- exp(b_cons + b_educ*educ + e)
  m <- lm(log(inc) ~ educ)

  sim_list[[i]] <- simloglm(m, nsim_est = 2000, scenario = list(educ = c(1,20)))

}


coverage_fct <- function(x, true_value){
  mean(x[1,] < true_value & x[2,] > true_value)

}

mean_summary <- lapply(sim_list, function(x) get_summary(x, which_qoi = "mean"))
median_summary <- lapply(sim_list, function(x) get_summary(x, which_qoi = "median"))
mean_fd_summary <- lapply(sim_list, function(x) get_first_difference(x, which_qoi = "mean"))
median_fd_summary <- lapply(sim_list, function(x) get_first_difference(x, which_qoi = "median"))
mean_ratio_summary <- lapply(sim_list, function(x) get_ratio(x, which_qoi = "mean"))
median_ratio_summary <- lapply(sim_list, function(x) get_ratio(x, which_qoi = "median"))


# Coverage scenario educ = 1, mean
cov_e1 <- coverage_fct(sapply(mean_summary, function(x) x$quantiles[,1] ), true_value = exp(b_cons+b_educ*1+1/2))
# Coverage scenario educ = 20, mean
cov_e20 <- coverage_fct(sapply(mean_summary, function(x) x$quantiles[,2] ), true_value = exp(b_cons+b_educ*20+1/2))
# Coverage scenario educ = 1, median
cov_m1 <- coverage_fct(sapply(median_summary, function(x) x$quantiles[,1] ), true_value = exp(b_cons+b_educ*1))
# Coverage scenario educ = 20, mean
cov_m20 <- coverage_fct(sapply(median_summary, function(x) x$quantiles[,2] ), true_value = exp(b_cons+b_educ*20))
# Coverage FD, mean
cov_fd_e <- coverage_fct(sapply(mean_fd_summary, function(x) x$quantiles ), true_value = exp(b_cons+b_educ*20 + 1/2) - exp(b_cons+b_educ*1 + 1/2))
# Coverage FD, median
cov_fd_m <- coverage_fct(sapply(median_fd_summary, function(x) x$quantiles ), true_value = exp(b_cons+b_educ*20 ) - exp(b_cons+b_educ*1 ))
# Coverage rate ratio, mean
cov_r_e <- coverage_fct(sapply(mean_ratio_summary, function(x) x$quantiles ), true_value = exp(b_cons+b_educ*20 + 1/2 ) / exp(b_cons+b_educ*1 + 1/2 ))
# Coverage rate ratio, median
cov_r_m <- coverage_fct(sapply(median_ratio_summary, function(x) x$quantiles ), true_value = exp(b_cons+b_educ*20 ) / exp(b_cons+b_educ*1 ))


out_tab <-
  data.frame("Estimand" = c("$\\text{Med}(y|\\text{edu}=1)$",
                            "$\\text{Med}(y|\\text{edu}=20)$",
                            "$E(y|\\text{edu}=1)$",
                            "$E(y|\\text{edu}=20)$",
                            "$\\text{Med}(y|\\text{edu}=20) - \\text{Med}(y|\\text{edu}=1)$",
                            "$E(y|\\text{edu}=20) - E(y|\\text{edu}=1)$",
                            "$\\text{Med}(y|\\text{edu}=20) / \\text{Med}(y|\\text{edu}=1)$",
                            "$E(y|\\text{edu}=20) / E(y|\\text{edu}=1)$"),
             "Coverage" =
               c(cov_m1,
                 cov_m20,
                 cov_e1,
                 cov_e20,
                 cov_fd_m,
                 cov_fd_e,
                 cov_r_m,
                 cov_r_e))

library(xtable)

print(xtable::xtable(out_tab,
                     digits = 3,
                     caption = "Monte Carlo Results for the coverage rate of simulated 95\\%-Confidence Intervals.",
                     label = "tab:mc_res"),
      include.rownames = F,
      booktabs = T,
      sanitize.text.function=function(x){x})

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