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

poisson_mle_sims_bias_comp.R

library(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, 2)  # true coefficients
lambda <- exp(beta[1] + beta[2]*x)  # implied mean
n_sims <- 1000  # 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])
}

se <- function(x) {
  sd(x)/sqrt(n_sims)
}

# plot expected values and cis
sims_df_sep <- data.frame(true_qi = exp(beta[1] + beta[2]*x0)*beta[2], 
                          mle_ev = apply(tau_hat_mle, 2, mean),
                          avg_ev = apply(tau_hat_avg, 2, mean),
                          med_ev = apply(tau_hat_median, 2, mean),
                          med_pre_ev = apply(tau_hat_median_pre, 2, mean),
                          mean_pre_ev = apply(tau_hat_mean_pre, 2, mean),
                          mle_se = apply(tau_hat_mle, 2, se),
                          avg_se = apply(tau_hat_avg, 2, se),
                          med_se = apply(tau_hat_median, 2, se),
                          med_pre_se = apply(tau_hat_median_pre, 2, se),
                          mean_pre_se = apply(tau_hat_mean_pre, 2, se))
sims_ev_df <- sims_df_sep %>%
  gather(method, ev, ends_with("_ev")) %>%
  select(-ends_with("_se")) %>%
  separate(method, into = c("method"), extra = "drop") %>%
  glimpse()
sims_se_df <- sims_df_sep %>%
  gather(method, se, ends_with("_se")) %>%
  select(-ends_with("_ev")) %>%
  separate(method, into = c("method"), extra = "drop") %>%
  glimpse()
sims_df <- left_join(sims_ev_df, sims_se_df) %>%
  glimpse()

ggplot(sims_df, aes(x = true_qi, y = ev, ymin = ev - 2*se, ymax = ev + 2*se, color = method)) + 
  geom_pointrange()

# t-test median and mle
p_value <- est <- lwr <- upr <- numeric(n_qi)
for (i in 1:n_qi) {
  print(i)
  t_test <- t.test(tau_hat_mle[, i] - tau_hat_median[, i])
  p_value[i] <- t_test$p.value
  est[i] <- t_test$est
  lwr[i] <- t_test$conf.int[1]
  upr[i] <- t_test$conf.int[2]
}
test_df <- data.frame(p_value, lwr, upr, true = exp(beta[1] + beta[2]*x0)*beta[2])
ggplot(test_df, aes(x = true, y = est, ymin = lwr, ymax = upr)) + 
  geom_pointrange()