Benchmarking Firth's estimators

benchmark-firth.R

# install crdata package to get weisiger2014 data set
remotes::install_github("carlislerainey/crdata")

# load packages
library(tidyverse)
library(brglm2)
library(logistf)
library(microbenchmark)

# load data
weis <- crdata::weisiger2014

# rescale weisiger2014 explanatory variables using arm::rescale()
rs_weis <- weis %>%
  mutate(across(polity_conq:coord, arm::rescale)) 

# create functions to fit models
f <- resist ~ polity_conq + lndist + terrain + soldperterr + gdppc2 + coord
f1 <- function() {
  glm(f, data = rs_weis, family = "binomial")
}
f2 <- function() {
  glm(f, data = rs_weis, family = "binomial", method = brglmFit)
}
f3 <- function() {
  logistf(f, data = rs_weis)
}
f4 <- function() {
  logistf(f, data = rs_weis, pl = FALSE)
}

# results don't exactly match (different convergence criteria?)
cbind(coef(f2()), coef(f4()))

# do benchmarking
bm <- microbenchmark("regular glm()" = f1(), 
               "brglm2" = f2(), 
               "logistf (default)" = f3(),
               "logistf (w/ pl = FALSE)" = f4(),
               times = 100) 

print(bm)

# plot times
bm %>%
  group_by(expr) %>%
  summarize(avg_time = mean(time)*10e-5) %>%  # convert to milliseconds
  ggplot(aes(x = fct_rev(expr), y = avg_time)) + 
  geom_col() + 
  labs(x = "Method", 
       y = "Avg. Time (in milliseconds)") + 
  coord_flip()