
# Star Wars at Central Banks Script - Apply p-curve method 
# Last Updated: 24 February 2020

library(tidyverse)
library(poibin)

# create list of sub-group data frames
df_list <- list(full_brodeur = full_brodeur_data, full_cb = full_cb_data, exploratory_cb = exploratory_cb_data, 
                datadriven_cb = datadriven_cb_data, total_less_exp_and_dd = total_cb_less_exp_and_dd)

# create new variables
raw_chi <- map(df_list, function(x) x[["t_stat_raw"]]^2) # z-score^2 ~ chi-squared with 1 df (we assume t-tests have enough df to approximate z-scores)
names(raw_chi) <- names(df_list)

sample_counts <- map(raw_chi, length) # number of test statistics included in each sub group
dof <- map(sample_counts, function(x) rep(1, times = x)) # some of the cho square functions used later require degrees of freedom inputs (always 1 for z scores)
ncp33 <- map(dof, function(x) map(x, function(y) unlist(get_chi_ncp(df = y, power = 1/3))))

# compute p-values
pvalues <- map(raw_chi, function(x) as.numeric(unlist(1 - pchisq(x, df = 1))))

# create p-curve graph data
full_curve <- map(pvalues, function (x) get_pcurve_data(p = data.frame(pvalues = x), alpha = 0.05))
triple_star_curve <- map(pvalues, function (x) get_pcurve_data(p = data.frame(pvalues = x), alpha = 0.01))

# compute pp-values
pp_full <- map(pvalues, function (x) get_pp_values(p = x, alpha = 0.05))
pp_half <- map(pvalues, function (x) get_pp_values(p = x, alpha = 0.025))
pp_triple_star <- map(pvalues, function (x) get_pp_values(p = x, alpha = 0.01))

# run meta-analysis statistical tests for p-curve inferences
#  Stouffer method (prefered)
stouffer_full_chi <- map(pp_full, stouffer)
stouffer_full_pval <- map(stouffer_full_chi, pnorm)

stouffer_half_chi <- map(pp_half, stouffer)
stouffer_half_pval  <- map(stouffer_half_chi, pnorm)

stouffer_triplestar_chi <- map(pp_triple_star, stouffer)
stouffer_triplestar_pval <- map(stouffer_triplestar_chi, pnorm)

#  Fisher method (robustness check)
fisher_full <- map(pp_full, function (x) 1 - pchisq(fisher(x), df = sum(!is.na(x))))
fisher_half <- map(pp_half, function (x) 1 - pchisq(fisher(x), df = sum(!is.na(x))))
fisher_triplestar <- map(pp_triple_star, function (x) 1 - pchisq(fisher(x), df = sum(!is.na(x))))
