In martinmodrak/revize-rs: Code for analysis of revize-rs project data

library(rstan)
library(tidybayes)
library(tidyverse)
library(here)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
source(here("R","simulator.R"))
source(here("R","evaluation_tools.R"))
source(here("R","sampling_multi.R"))
source(here("R","sbc.R"))

model <- stan_model(here("stan","main_model.stan"))

N_questions = 100; N_fixed = 10
corr_matrix <- rlkjcorr(1, N_questions, eta = 0.1)
corr_matrix_chol <- t(chol(corr_matrix)) #Transpose to lower triangular to match Stan
#b_sd <- abs(rnorm(N_fixed, 0, 1))
b_sd <- rep(1, length.out = N_fixed)
mult_res <- (corr_matrix_chol %*% matrix(rnorm(N_fixed * N_questions, 0, 1), nrow = N_questions, ncol = N_fixed)) %*% diag(b_sd)
#corr_matrix - cor(t(mult_res))
mult_sd_diff <- b_sd - apply(mult_res, FUN = sd, MARGIN = 1)
mean(mult_sd_diff)
hist(mult_sd_diff)

mvn_res <- t(rmvnorm(N_fixed, sigma = corr_matrix)) #%*% diag(b_sd)
mvn_sd_diff <- b_sd - apply(mvn_res, FUN = sd, MARGIN = 1)
mean(mvn_sd_diff)
hist(mvn_sd_diff)
#corr_matrix - cor(t(mvn_res))
#tibble(mvn = as.numeric(mvn_res), mult = as.numeric(mult_res)) %>% ggplot(aes(mvn, mult)) + geom_point(alpha = 0.1) + geom_smooth(method = "lm")

ncat = 4
repeat {
  test_data <- simulate_data(N = 200, N_questions = 4, N_fixed = 2, ncat = 4, N_random_groups = c(2,4),  N_monotonic = 2, questions_correlation = TRUE)
  if(length(unique(test_data$observed$Y)) == ncat) {
    break;
  }
  cat("Reject\n")
}
fit <- sampling(model, data = test_data$observed, control = list(adapt_delta = 0.95))
check_all_diagnostics(fit)
evaluation_summary(fit, test_data$true)
#evaluation_summary(fit, test_data$computed)

#set.seed(6847567)
sbc_steps <- 10
sbc_results <- sbc(model, function() { simulate_data_reject(N = 200, N_questions = 4, N_fixed = 2, ncat = 4, N_random_groups = c(2,4),  N_monotonic = 2, questions_correlation = TRUE)}, N_steps = sbc_steps, control = list(adapt_delta = 0.9))
saveRDS(sbc_results, "sbc_100.rds")
#saveRDS(sbc_results, "sbc_no_mono_no_cor_no_random.rds")
pattern <- "zeta|q_corr|b_sd"
plot_sbc_params(sbc_results$params %>% filter(!grepl(pattern, param_name)), binwidth = 10)
plot_sbc_params(sbc_results$params %>% filter(grepl(pattern, param_name)), binwidth = 10)
summarise_sbc_diagnostics(sbc_results)

sbc_results$params %>% filter(grepl("^b\\[", param_name)) %>% sbc_power()
sbc_results$params %>% filter(grepl("b_monotonic", param_name)) %>% sbc_power()
sbc_results$params %>% filter(grepl("b_random", param_name)) %>% sbc_power()

sbc_orig <- readRDS("sbc_500.rds")
plot_sbc_params(sbc_orig$params %>% filter(!grepl("zeta|q_corr", param_name)))
plot_sbc_params(sbc_orig$params %>% filter(grepl("zeta|q_corr", param_name)))

model_code_order_test <- "
data {
  int N;
  int sigma;
}

parameters {
  ordered[N] x;
}

model {
  x ~ student_t(3, 0, sigma);
}
"
model_order_test <- stan_model(model_code = model_code_order_test)

N = 1
sigma = 3
sbc_order_results <- sbc(model, function() { list(observed = list(N = N, sigma = sigma), true = list(x = sort(rt(N, 3) * sigma)) }, N_steps = 1000)
#saveRDS(sbc_results, "sbc_500.rds")
saveRDS(sbc_results, "sbc_no_mono_no_cor_no_random.rds")
pattern <- "zeta|q_corr|b_sd"
plot_sbc_params(sbc_results$params %>% filter(!grepl(pattern, param_name)), binwidth = 5)