Replication material for

transformation_bias_blog.R

library(pacman)

p_load(MASS, tidyverse)

# set seed
set.seed(8764)

# set simulation paramters
n <- 100  # sample size
x <- rnorm(n)  # single explanatory variable
n_qi <- 100  # number of points at which to calculate the marginal effect
x0 <- seq(-3, 3, length.out = n_qi)  # points at which to calculate the marginal effect
beta <- c(-2, 1)  # true coefficients
lambda <- exp(beta[1] + beta[2]*x)  # implied mean
n_sims <- 10000  # number of mc simulations

tau_hat_mean_pre <- tau_hat_median_pre <- tau_hat_median <- tau_hat_mle <- tau_hat_avg <- matrix(NA, nrow = n_sims, ncol = n_qi)
for (i in 1:n_sims) {
    y <- rpois(n, lambda = lambda)
    fit <- glm(y ~ x, family = poisson)
    beta_hat <- coef(fit)
    Sigma_hat <- vcov(fit)
    beta_tilde <- MASS::mvrnorm(1000, mu = beta_hat, Sigma = Sigma_hat)
    tau_tilde <- t(exp(cbind(1, x0)%*%t(beta_tilde)))*beta_tilde[, 2]
    tau_hat_avg[i, ] <- apply(tau_tilde, 2, mean)
    tau_hat_mle[i, ] <- exp(beta_hat[1] + beta_hat[2]*x0)*beta_hat[2]
    tau_hat_median[i, ] <- apply(tau_tilde, 2, median)
    tau_hat_median_pre[i, ] <- exp(median(beta_tilde[, 1]) + median(beta_tilde[, 2])*x0)*median(beta_tilde[, 2])
    tau_hat_mean_pre[i, ] <- exp(mean(beta_tilde[, 1]) + mean(beta_tilde[, 2])*x0)*mean(beta_tilde[, 2])
}

sims_df_sep <- data.frame(true_qi = exp(beta[1] + beta[2]*x0)*beta[2], 
                          mle = apply(tau_hat_mle, 2, mean),
                          avg = apply(tau_hat_avg, 2, mean),
                          med = apply(tau_hat_median, 2, mean),
                          med_pre = apply(tau_hat_median_pre, 2, mean),
                          mean_pre = apply(tau_hat_mean_pre, 2, mean)) 
sims_df <- sims_df_sep %>%
    gather(method, ev, mle, avg, med, med_pre, mean_pre)

hist(with(sims_df_sep, (mle - true_qi) - (med - true_qi)))
hist(with(sims_df_sep, (mle - true_qi) - (med_pre - true_qi)))
hist(with(sims_df_sep, (mle - true_qi) - (mean_pre - true_qi)))

ggplot(sims_df, aes(x = true_qi, y = ev - true_qi, 
                    linetype = method, color = method)) + 
    geom_line() + 
    theme_bw() + 
    labs(title = "Bias in Estimates of Poisson Marginal Effects",
         x = "True Marginal Effect",
         y = "Bias in Estimates of Marginal Effect",
         linetype = "Method", colour = "Method")