Compare bias in Poisson MLE and various Sim methods (10 August)

poisson_mle_sims_bias_comp.R

library(pacman)

p_load(MASS, tidyverse)

# set seed
set.seed(4321)

# 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)

par(mfrow = c(1, 3))
plot(sims_df_sep$true_qi, with(sims_df_sep, (med_pre - true_qi) - (mle - true_qi)), 
     ylab = "Bias of Median of Simulated Betas for QI vs. MLE", xlab = "True Value")
plot(sims_df_sep$true_qi, with(sims_df_sep, (mean_pre - true_qi) - (mle - true_qi)), 
     ylab = "Bias of Mean of Simulated Betas for QI vs. MLE", xlab = "True Value")
plot(sims_df_sep$true_qi, with(sims_df_sep, (med - true_qi) - (mle - true_qi)), 
     ylab = "Bias of Median of Simulated QI vs. MLE", xlab = "True Value")