# Resampling reporting propensity data
# run on laptop

rm(list=ls())
library(rstan)
library(Hmisc)
library(data.table)
library(grid)
library(gridExtra)


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

setwd(path_to_data)
 
report_propensity = readRDS('report_propensity.Rds')
users = readRDS('user_reporting.Rds')

describe(users$n_reports)

# taking just those report propensities from users who sent at any report
users$any_reports = users$n_reports > 0
reporters = users$user_UUID[users$any_reports]
report_propensity = report_propensity[report_propensity$user %in% reporters,]

report_propensity$any_reports[is.na(report_propensity$any_reports)] = FALSE

describe(report_propensity$any_reports)


rp_orig = report_propensity

# removing extreme values
report_propensity = report_propensity[which(report_propensity$participation_days <= as.numeric(quantile(report_propensity$participation_days, .95)) & report_propensity$season_day > as.numeric(quantile(report_propensity$season_day, .05)) & report_propensity$season_day < as.numeric(quantile(report_propensity$season_day, .95))),]


# sampling weights for resampling
participation_days_sampling_weights = data.table('participation_days' = as.integer(names(table(report_propensity$participation_days))), 'participation_day_weight' = nrow(report_propensity)/as.numeric(table(report_propensity$participation_days)))

report_propensity = merge(report_propensity, participation_days_sampling_weights, by='participation_days', all.x=T, all.y=F)

samp_size = 10000 # sample size for resampling

# resampling
rp = report_propensity[sample(1:nrow(report_propensity), samp_size, prob=report_propensity$participation_day_weight, replace=F),]


ggplot(data=report_propensity, aes(x=participation_days)) + geom_histogram(binwidth = 1, colour="#cccccc") + xlab("Total Participation Time (days)") + ylab("Participants (N)") + theme_bw() + geom_histogram(data=rp, aes(x=participation_days), binwidth = 1, colour="black")
ggsave(paste0(path_to_figures, 'imbalance.png'), width=7, height=6, dpi=1000)
ggsave(paste0(path_to_figures, 'imbalance.pdf'), width=7, height=6)

mean(rp$any_reports)
mean(report_propensity$any_reports)


rp$participation_days_std2 = std2(rp$participation_days)
rp$participation_days_std2_sq = rp$participation_days_std2^2
rp$participation_days_std2_cu = rp$participation_days_std2^3

saveRDS(rp, file='reporting_propensity_resampled.Rds')

