setwd("trust-jrssa")

library(tidyr)
library(dplyr)
library(brms)

options(mc.cores = parallel::detectCores())

country <- "ie"

nchains <- 8
niter <- 1000
nwarmup <- 500

sum_coding <- function(x, lvls = levels(x)) {
  # codes the first category with -1
  nlvls <- length(lvls)
  stopifnot(nlvls > 1)
  cont <- diag(nlvls)[, -nlvls, drop = FALSE]
  cont[nlvls, ] <- -1
  cont <- cont[c(nlvls, 1:(nlvls - 1)), , drop = FALSE]
  colnames(cont) <- lvls[-1]
  x <- factor(x, levels = lvls)
  contrasts(x) <- cont
  x
}


df <- readRDS("data/all_cat_27_edu3_subset_2_20211212.rds") %>%
  filter(t_year >= 1989, 
         t_cntry == toupper(country)) %>%
  mutate(caseid = row_number()) %>%
  select(caseid, t_project, t_wave, t_year, t_cntry, survey, 
         weight_total_despst, weight_total_des,
         sex, age_cat, educ3, thres, trust_parl, trust_jus, trust_polpart) %>%
  gather(inst, trust, c("trust_parl", "trust_jus", "trust_polpart")) %>%
  mutate(Intercept = 1,
         survey_inst = paste0(survey, inst, sep = "_"),
         age_cat = factor(age_cat),
         sex = factor(sex),
         educ3 = factor(educ3),
         project_f = sum_coding(factor(t_project)),
         inst_f = sum_coding(factor(inst)))
	
	
m2 <- brm(
  trust | thres(thres, gr = survey_inst) + weights(weight_total_despst, scale = TRUE) ~ 
    Intercept + inst_f*project_f +  (1 | caseid) + s(t_year) + 
    s(t_year, age_cat, bs="fs") + s(t_year, sex, bs="fs") + s(t_year, educ3, bs="fs"),
  data = df,
  family = cumulative("probit", threshold = "sum_to_zero"),
  prior = prior(normal(0, 1.5), class = "b") +
    prior(normal(0, 2), class = "sds") +
    prior(normal(0, 1.5), class = "Intercept"),
  chains = nchains,
  iter = niter,
  warmup = nwarmup,
  inits = 0,
  file = paste0("models/", country, "_despst")
)


