library(tidyverse)

setwd('/path/to/data')

t_results <- tibble(
  homophily_rate_citedf = c(
    0.116,
    0.081,
    0.078,
    0.212,
    0.133,
    0.124,
    0.068,
    0.16,
    0.214,
    0.085,
    0.102,
    0.099
  ),
  study = c(
    'Ferber (1986)',
    'Ferber (1988)',
    'Lutz (1990)',
    'Davenport and Snyder (1995)',
    'Mcelhinny, Hols, Holtzkener,\nUnger, and Hicks (2003)',
    'Håkanson (2005)',
    'Ferber and Brün (2011)',
    'Knobloch-Westerwick\nand Glynn (2013)',
    'Mitchell, Lange, and Brus (2013)',
    'Potthoff and Zimmermann (2017)',
    'Dion, Sumner,\nand Mitchell (2018)',
    'Ghiasi, Mongeon, Sugimoto,\nand Larivière (2018)'
  ),
  n_citation_links = c(
    2394,
    11669,
    10593,
    4951,
    16766,
    23483,
    3256,
    2958,
    3013,
    25853,
    30066,
    20395382
  ),
  citingm_citedf = c(
    94,
    516,
    1085,
    423,
    2113,
    3005,
    112,
    411,
    242,
    1900,
    NA,
    2711742
  ),
  citingm = c(
    1597,
    6100,
    6957,
    3443,
    9715,
    13139,
    1608,
    1956,
    2389,
    21779,
    NA,
    9682639
  ),
  citingf_citedf = c(
    139,
    920,
    851,
    505,
    2463,
    3646,
    227,
    365,
    197,
    700,
    NA,
    3885636
  ),
  citingf = c(
    797,
    5569,
    3636,
    1508,
    7051,
    10344,
    1648,
    987,
    624,
    4074,
    NA,
    10254123
  ),
  study_order = 1:12
)
t_results$study <- factor(t_results$study, levels = t_results$study[order(t_results$study_order)])

t_results$mciting_p <- t_results$citingm_citedf / t_results$citingm
t_results$mciting_var <- t_results$mciting_p * (1 - t_results$mciting_p)
t_results$fciting_p <- t_results$citingf_citedf/ t_results$citingf
t_results$fciting_var <- t_results$fciting_p * (1 - t_results$fciting_p)

t_results$threshold <- 
  qt(0.975, t_results$citingm + t_results$citingf - 2)

t_results$ci_l <- 
  (t_results$fciting_p - t_results$mciting_p) - 
  t_results$threshold * sqrt(t_results$mciting_var / t_results$citingm + t_results$fciting_var / t_results$citingf)
t_results$ci_u <- 
  (t_results$fciting_p - t_results$mciting_p) + 
  t_results$threshold * sqrt(t_results$mciting_var / t_results$citingm + t_results$fciting_var / t_results$citingf)

dion_data <- matrix(nrow = 5, ncol = 4)
colnames(dion_data) <- c('citingm', 'citingm_citedf', 'citingf', 'citingf_citedf')
dion_data[1, ] <- c(10489, 1178, 2762, 702)
dion_data[2, ] <- c(629, 313, 3982, 2500)
dion_data[3, ] <- c(4476, 246, 123, 10)
dion_data[4, ] <- c(8207, 329, 441, 51)
dion_data[5, ] <- c(2731, 214, 331, 51)
dion_data <- as_tibble(dion_data)

dion_data$mciting_p <- dion_data$citingm_citedf / dion_data$citingm
dion_data$fciting_p <- dion_data$citingf_citedf / dion_data$citingf
dion_data$mciting_var <- dion_data$mciting_p * (1 - dion_data$mciting_p)
dion_data$fciting_var <- dion_data$fciting_p * (1 - dion_data$fciting_p)
dion_data$homophily <- dion_data$fciting_p - dion_data$mciting_p
dion_data$total_links <- dion_data$citingm + dion_data$citingf
dion_data$weighted_hom <- dion_data$homophily * dion_data$total_links / sum(dion_data$total_links)
dion_within_var_f <- sum(dion_data$citingf * dion_data$fciting_p * (1 - dion_data$fciting_p)) / sum(dion_data$citingf)
dion_within_var_m <- sum(dion_data$citingm * dion_data$mciting_p * (1 - dion_data$mciting_p)) / sum(dion_data$citingm)
dion_between_var_f <- var(dion_data$fciting_p)
dion_between_var_m <- var(dion_data$mciting_p)
dion_var_f <- dion_within_var_f + dion_between_var_f
dion_var_m <- dion_within_var_m + dion_between_var_m

dion_threshold <- qt(0.975, sum(dion_data$total_links - 2))

dion_weighted_homophily <- sum(dion_data$weighted_hom)

dion_ci_l <- dion_weighted_homophily - dion_threshold * sqrt(dion_var_m / sum(dion_data$citingm) + dion_var_f / sum(dion_data$citingf))
dion_ci_u <- dion_weighted_homophily + dion_threshold * sqrt(dion_var_m / sum(dion_data$citingm) + dion_var_f / sum(dion_data$citingf))
t_results$ci_l[str_starts(t_results$study, 'Dion')] <- dion_ci_l
t_results$ci_u[str_starts(t_results$study, 'Dion')] <- dion_ci_u

t_results$point_size <- sqrt(t_results$n_citation_links)

meta_plot <- ggplot(data = t_results) +
  geom_point(aes(x = homophily_rate_citedf, y = study, size = point_size), shape = 19) +
  geom_linerange(aes(x = homophily_rate_citedf, y = study, xmin = ci_l, xmax = ci_u)) +
  scale_x_continuous(breaks = seq(0, .25, by = .05), limits = c(0, .255)) +
  theme_classic() +
  theme(
    axis.text.x = element_text(size = 10, colour = 'black'),
    axis.text.y = element_text(size = 10, colour = 'black'),
    axis.title.x = element_text(size = 12, colour = 'black'),
    axis.title.y = element_blank(),
    legend.position = 'none'
  ) +
  labs(x = 'Gender homophily rate') +
  geom_hline(yintercept = 1.5, colour = 'grey') +
  geom_hline(yintercept = 2.5, colour = 'grey') +
  geom_hline(yintercept = 3.5, colour = 'grey') +
  geom_hline(yintercept = 4.5, colour = 'grey') +
  geom_hline(yintercept = 5.5, colour = 'grey') +
  geom_hline(yintercept = 6.5, colour = 'grey') +
  geom_hline(yintercept = 7.5, colour = 'grey') +
  geom_hline(yintercept = 8.5, colour = 'grey') +
  geom_hline(yintercept = 9.5, colour = 'grey') +
  geom_hline(yintercept = 10.5, colour = 'grey') +
  geom_hline(yintercept = 11.5, colour = 'grey')

ggsave('meta_plot.png', plot = meta_plot, device = 'png', width = 160, height = 120, units = 'mm')