# Bayesian hierarchical model of reporting report_propensity, implemented in stan

rm(list=ls())
library(rstan)
library(Hmisc)
library(rstanarm)
library(loo)
library(data.table)

setwd('~/research/gTiger/Article_citsci/open_code')
source('settings/parameters.r')
source('settings/paths.r')
source('settings/functions.r')

setwd(path_to_data)

ncores = 4

rstan_options(auto_write = TRUE)
options(mc.cores = ncores)

stan_settings = list(iterations=1000, chains=options()$mc.cores, warmup=500)

invstd2 = function(x, d) (x/(2*sd(d, na.rm=T)))

D = readRDS('reporting_propensity_resampled.Rds')

D$trials = 1

M00 = stan_glmer(cbind(any_reports, trials-any_reports)~1 + (1 | user),data=D, family=binomial(link = "logit"), prior = student_t(df = 1, location=0, scale=2.5), prior_intercept = student_t(df = 1, location=0, scale=2.5), chains = ncores, cores = ncores)


M0 = stan_glmer(cbind(any_reports, trials-any_reports)~participation_days_std2 + (1 | user),data=D, family=binomial(link = "logit"), prior = student_t(df = 1, location=0, scale=2.5), prior_intercept = student_t(df = 1, location=0, scale=2.5), chains = ncores, cores = ncores)

M1 = stan_glmer(cbind(any_reports, trials-any_reports)~participation_days_std2 + participation_days_std2_sq + (1 | user), data=D, family=binomial(link = "logit"), prior = student_t(df = 1, location=0, scale=2.5), prior_intercept = student_t(df = 1, location=0, scale=2.5), chains = ncores, cores = ncores)

M2 = stan_glmer(cbind(any_reports, trials-any_reports)~participation_days_std2 + participation_days_std2_sq + participation_days_std2_cu + (1 | user), data=D, family=binomial(link = "logit"), prior = student_t(df = 1, location=0, scale=2.5), prior_intercept = student_t(df = 1, location=0, scale=2.5), chains = ncores, cores = ncores)

M00_loo = loo(M00) 
M0_loo = loo(M0) 
M1_loo = loo(M1) 
M2_loo = loo(M2) 

M00_waic = waic(M00) 
M0_waic = waic(M0) 
M1_waic = waic(M1) 
M2_waic = waic(M2)

M0_loo
M1_loo
M2_loo
compare(M2_loo, M1_loo, M0_loo, M00_loo)

compare(M2_waic, M1_waic, M0_waic, M00_waic)

t_coef_M00 = coef(M00)
t_coef_M0 = coef(M0)
t_coef_M1 = coef(M1)
t_coef_M2 = coef(M2)

t_se_M00 = se(M00)
t_se_M0 = se(M0)
t_se_M1 = se(M1)
t_se_M2 = se(M2)

t_waics = list(M00_waic,M0_waic, M1_waic, M2_waic)
t_loos = list(M00_loo,M0_loo, M1_loo, M2_loo)
t_ses = list(t_se_M00, t_se_M0, t_se_M1, t_se_M2)
t_coefs = list(t_coef_M00, t_coef_M0, t_coef_M1, t_coef_M2)

names(t_waics) = c("M00","M0", "M1", "M2")
names(t_loos) = c("M00","M0", "M1", "M2")
names(t_ses) = c("M00","M0", "M1", "M2")
names(t_coefs) = c("M00","M0", "M1", "M2")

for(i in 1:length(t_waics)){
  t_waics[[i]]$elpd_waic = round(t_waics[[i]]$elpd_waic,0)
  t_waics[[i]]$se_elpd_waic = round(t_waics[[i]]$se_elpd_waic,0)
}
for(i in 1:length(t_loos)){
  t_loos[[i]]$elpd_loo = round(t_loos[[i]]$elpd_loo,0)
  t_loos[[i]]$se_elpd_loo = round(t_loos[[i]]$se_elpd_loo,0)
}

t_coef_M0$user$participation_days_std2[1]


t_models_table = data.table("p"=c("participation time (days)","", "participation time (days) sq.","", "participation time (days) cu","", "ELPD", "", "N"),
                            M00=c(
                              NA,
                              NA,
                              NA, 
                              NA, 
                              NA, 
                              NA, 
                              t_loos[["M00"]]$elpd_loo,t_loos[["M00"]]$se_elpd_loo, nrow(D)), 
                            M0=c(
                              invstd2(t_coef_M0$user$participation_days_std2[1], D$participation_days),
                              invstd2(t_se_M0['participation_days_std2'],D$participation_days),
                              NA, 
                              NA, 
                              NA, 
                              NA, 
                              t_loos[["M0"]]$elpd_loo,t_loos[["M0"]]$se_elpd_loo, nrow(D)), 
                            M1=c(
                              invstd2(t_coef_M1$user$participation_days_std2[1], D$participation_days),
                              invstd2(t_se_M1['participation_days_std2'],D$participation_days),
                              invstd2(t_coef_M1$user$participation_days_std2_sq[1], D$participation_days),
                              invstd2(t_se_M1['participation_days_std2_sq'],D$participation_days),
                              NA, 
                              NA, 
                              t_loos[["M1"]]$elpd_loo,t_loos[["M1"]]$se_elpd_loo, nrow(D)), 
                            M2=c(
                              invstd2(t_coef_M2$user$participation_days_std2[1], D$participation_days),
                              invstd2(t_se_M2['participation_days_std2'],D$participation_days),
                              invstd2(t_coef_M2$user$participation_days_std2_sq[1], D$participation_days),
                              invstd2(t_se_M2['participation_days_std2_sq'],D$participation_days),
                              invstd2(t_coef_M2$user$participation_days_std2_cu[1], D$participation_days),
                              invstd2(t_se_M2['participation_days_std2_cu'],D$participation_days),                              t_loos[["M2"]]$elpd_loo,t_loos[["M2"]]$se_elpd_loo, nrow(D)))




saveRDS(list(M00, M0, M1, M2), file='untracked/reporting_propensity_model_comparison_models.Rds')

saveRDS(t_models_table, file='reporting_propensity_model_comparison_table.Rds')

saveRDS(list(t_coefs, t_ses, t_loos, t_waics), file='reporting_propensity_model_comparison_summaries.Rds')



