# Bayesian hierarchical model of reporting report_propensity, implemented in stan
# run on cluster with 20 cores (but this could be modified)

rm(list=ls())
library(rstan)
library(Hmisc)

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

setwd(path_to_data)

ncores = 20
# ncores = 4

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

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

D = readRDS('reporting_propensity_resampled.Rds')

D = data.frame(D)
 
these_predictor_names=c('participation_days_std2', 'participation_days_std2_sq','participation_days_std2_cu') 

test = make_stan_data(data=D, response_name='any_reports', predictor_names=these_predictor_names, group_names=c('user'))

test$for_stan$n_level_1

result = stanml(data=D, response_name='any_reports', predictor_names=these_predictor_names, group_names=c('user'), model_code=paste0(path_to_code, 'stan/reporting_propensity_model_logit_RI_1L.stan'), model_name='report_propensity', machine_name="cluster", save=TRUE, print=TRUE, plot=FALSE, traceplot=TRUE, dir=path_to_untracked_output, nchains=stan_settings$chains, iterations=stan_settings$iterations, warmup=stan_settings$warmup, include=FALSE, pars='eta') 

fit = result$fit
this_data = result$this_data

names(this_data)

x = extract(fit)
bmat = cbind(x$mu_a, x$b)
colnames(bmat) = c('mu_user_intercepts', these_predictor_names)
bs = sapply(1:ncol(bmat), function(i) mean(bmat[,i]))
names(bs) = c('mu_user_intercepts', these_predictor_names)
nsims = 400
bsims = bmat[sample(1:nrow(bmat), nsims),]

user_intercepts = x$a
user_intercepts = sapply(1:ncol(user_intercepts), function(i) mean(user_intercepts[,i]))

saveRDS(list(rp=D,this_data=this_data, bs=bs, bsims=bsims, user_intercepts=user_intercepts), file='reporting_propensity_sampled_stanfit.Rds')
