meta-regression-sundell-2014.R

library(pacman)
p_load(tidyverse, janitor, brms)
tibble::tribble(~Year, ~OR, ~z, 
  1765, 1.80, 1.69,
  1775, 2.41, 2.86,
  1790, 2.36, 2.67,
  1805, 2.58, 3.07,
  1820, 1.83, 2.63,
  1835, 1.55, 1.48,
  1850, 1.06, 0.14,
  1866, 1.77, 1.32,
  1881, 1.60, 1.11,
  1895, 0.85, 0.23,
  1910, 1.27, 0.33,
  1925, 2.01, 1.96,
  1940, 1.02, 0.07,
  1955, 1.36, 1.52,
  1970, 0.96, 0.26) %>% 
  mutate(logOR= log(OR), 
         se = abs(logOR/z)) -> out2

brm(formula=logOR|se(se)~gp(Year), data=out2, control=list(adapt_delta=.999, max_treedepth=20), 
    warmup=1500, iter=8000, file="brm_gp") -> brm_gp

brm_gp %>% conditional_effects() %>% .[[1]] -> estimated_effects

out2 %>% mutate(conf_low = (logOR-se*1.96) %>% exp,
                conf_high = (logOR+se*1.96) %>% exp) -> out3

estimated_effects %>% 
  as_tibble %>% 
  ggplot(aes(Year, y=estimate__ %>% exp, 
             ymin=lower__ %>% exp, 
             ymax=upper__ %>% exp))+
  geom_line()+
  geom_ribbon(alpha=0.3)+
  geom_point(data=out3, 
             inherit.aes = F, 
             mapping = aes(x=Year, y=OR))+
  geom_errorbar(data=out3,
                aes(x=Year, ymin=conf_low, ymax=conf_high), 
                inherit.aes = F, 
                width=0.3)+
  geom_hline(yintercept = 1, linetype="dashed")+
  coord_trans(y="log10")+
  labs(y="Odds Ratio", x = "Year", 
       title="Meta-Regression of Estimates of a Nepotism Measure for Swedish Central Public Administration",
       subtitle="Estimates from \"Nepotism and Meritocracy\", Anders Sundell (2014)")+
  scale_x_continuous(breaks=seq(1740, 1980, 20))+
  scale_y_continuous(breaks=c(0.25, 0.5,1,2,3,4,5))