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

library(pwr)
library(here)
library(rstan)
library(tidyverse)
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"))

A quick frequentist estimate

power = 0.8
alpha = 0.05
power_vals <- tibble(n_predictors = c(10, 20, 30, 40, 100)) %>%
  #crossing(tibble(f2 = c(0.02, 0.15, 0.35))) %>%
  crossing(tibble(sample_size = seq(from = 200, to = 1500, length.out = 100) %>% round())) %>%
  rowwise() %>%
  mutate(f2 = pwr.f2.test(u = n_predictors, v = sample_size - n_predictors - 1, sig.level = alpha, power = power)$f2) %>%
  ungroup()


power_vals %>%
  mutate(n_predictors = as.factor(n_predictors)) %>%
  ggplot(aes(x = sample_size, y = f2, color = n_predictors, group = n_predictors )) +
  geom_line() +
  geom_hline(yintercept = 0.02) +
  geom_hline(yintercept = 0.15) +
  geom_hline(yintercept = 0.35) +
  scale_y_log10()

Using the SBC code and full simulations:

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

base_random_groups <- c(3,14,3,15,8,3,8,3,2,7,4,9,5,9,6)
#settings <- tibble(N = c(50, 75, 100)) %>% crossing(subset_questions =  c(11)) %>% 
settings <- tibble(N = c(500, 750, 1000)) %>% #crossing(subset_questions =  c(11, 22)) %>% 
  crossing(tibble(N_fixed = c(4, 6), N_monotonic = c(4, 6), N_random_groups = list(base_random_groups, c(base_random_groups, 6, 6, 6, 5, 6, 8, 8, 6)))) %>%
  mutate(N_monotonic_cat = 6, 
         #N_questions = 22, questions_correlation = TRUE, 
         N_questions = 1, questions_correlation = FALSE, 
         ncat = 7)


settings <- settings %>% rbind(tibble(N = c(250)) %>% 
  crossing(tibble(N_fixed = c(4, 6), N_monotonic = c(4, 6), N_random_groups = list(base_random_groups, c(base_random_groups, 6, 6, 6, 5, 6, 8, 8, 6)))) %>%
  mutate(N_monotonic_cat = 6, 
         #N_questions = 22, questions_correlation = TRUE, 
         N_questions = 1, questions_correlation = FALSE, 
         ncat = 7))

settings

sbc_steps = 50
sbc_results <- list()
#for(i in 1:nrow(settings)) {
for(i in rev(1:nrow(settings))) {
  result_file <- paste0("power_",i,".rds")
  if(file.exists(result_file)) {
    sbc_results[[i]] <- readRDS(result_file)
  } else {
    params <- as.list(settings[i, ])
    params$N_random_groups <- settings$N_random_groups[[i]]
    sbc_results[[i]] <- sbc(model, function() { do.call(simulate_data_reject, params)}, N_steps = sbc_steps, control = list(adapt_delta = 0.9))
    saveRDS(sbc_results[[i]], result_file)
  }
}

all_params_power <- sbc_results %>% imap(function(x, i) { x$params %>% 
    mutate(group = paste0("N = ",settings$N[i],", fixed=", settings$N_fixed[i])) }
    ) %>% do.call(rbind, .) %>% sbc_power()

sbc_power_plot(all_params_power %>% filter(grepl("^b\\[", param_name)), group = group)

sbc_power_plot(all_params_power %>% filter(grepl("^b_random\\[", param_name)), group = group)

sbc_power_plot(all_params_power %>% filter(grepl("^b_monotonic\\[", param_name)), group = group)

martinmodrak/revize-rs documentation built on March 9, 2021, 5:30 a.m.