Simulation study 1 repeating JGV2023

svp_01.R

styler::style_dir()

library(SimDesign)

design <- SimDesign::createDesign(
  lambda_32 = c(0, .2, .5, .7),
  n = c(seq(50, 150, 25), seq(200, 300, 50), c(400, 500, 1000)),
  group_iter = 1:10
)
# condition <- design[42, ]
# condition

generate <- function(condition, fixed_objects = NULL) {
  lambda_mat <- cbind(
    c(rep(.7, 3), rep(0, 3)),
    c(rep(0, 3), rep(.7, 3))
  )
  lambda_mat[3, 2] <- condition$lambda_32
  phi_mat <- matrix(c(1, .25, .25, 1), 2)
  sigma <- tcrossprod(tcrossprod(lambda_mat, phi_mat), lambda_mat) +
    .51 * diag(6)
  n_ast <- condition$n - 1
  s_mat <- rWishart(1, n_ast, sigma / n_ast)[, , 1]
  var_names <- paste0("x", seq_len(nrow(s_mat)))
  dat <- list(s_mat = s_mat, n = condition$n)
  dimnames(dat$s_mat) <- list(var_names, var_names)
  return(dat)
}
# dat <- generate(condition)
# dat

analyse <- function(condition, dat, fixed_objects = NULL) {
  model_syntax <- paste(
    paste0("F1 =~ ", paste0("x", 1:3, collapse = " + ")),
    paste0("F2 =~ ", paste0("x", 4:6, collapse = " + ")),
    sep = "\n"
  )
  mod <- minorbsem::minorbsem(
    model_syntax,
    sample_cov = dat$s_mat, sample_nobs = dat$n,
    method = "none", simple_struc = FALSE,
    warmup = 500, sampling = 500, show = FALSE
  )

  load_draws <- mod@stan_fit$draws(variables = paste0(
    "Load_mat[", c(4:6, 1:3), ",", c(rep(1, 3), rep(2, 3)), "]"
  ))
  load_rh_dt <- posterior::summarise_draws(load_draws, posterior::rhat)
  load_rh <- load_rh_dt$`posterior::rhat`
  names(load_rh) <- paste0(
    "l", gsub("Load_mat\\[|,|\\]", "", load_rh_dt$variable), "_rh"
  )
  load_rh <- load_rh[grepl("1_rh", names(load_rh))]

  load_mat <- posterior::as_draws_matrix(load_draws)
  load_mode <- apply(load_mat, 2, function(x) {
    x_dens <- density(x, bw = "SJ")
    m1 <- x_dens$x[which.max(x_dens$y)]
    m2 <- modeest::shorth(x)
    m3 <- modeest::hsm(x)
    m4 <- modeest::meanshift(x)
    c(m1, m2, m3, m4)
  })
  load_mode_vec <- as.vector(load_mode)
  names(load_mode_vec) <- paste(
    paste0(
      "l",
      rep(gsub("Load_mat\\[|,|\\]", "", colnames(load_mode)), rep(4, 6))
    ),
    paste0("m", rep(1:4, 6)),
    sep = "_"
  )

  maj_draws <- mod@stan_fit$draws(
    variables = c(
      "ppp", "Load_mat[3,2]", "Load_mat[1,2]", "Load_mat[2,2]", "phi_mat[2,1]",
      "gdp_loadings_complex", "sigma_loadings_complex"
    )
  )
  maj_sum <- posterior::summarise_draws(
    maj_draws,
    mean, median, sd, ~ quantile(.x, c(.025, .975)),
    posterior::rhat, posterior::ess_bulk, posterior::ess_tail
  )
  full_vec <- unlist(as.vector(maj_sum[, -1]))
  pars <- c("ppp", "l32", "l12", "l22", "p21", "dpa", "dps")
  stats <- c("mn", "md", "sd", "lo", "hi", "rh", "eb", "et")
  names(full_vec) <- paste(
    rep(pars, length(stats)),
    rep(stats, rep(length(pars), length(stats))),
    sep = "_"
  )
  full_vec <- full_vec[!is.na(full_vec)]
  ret_vec <- c(full_vec, load_mode_vec, load_rh)
  return(ret_vec)
}
# ana <- analyse(dat = dat)
# ana

summarise <- function(condition, results, fixed_objects = NULL) {
  colMeans(results, na.rm = TRUE)
}

res <- runSimulation(
  design = design, generate = generate,
  replications = 100,
  analyse = analyse, summarise = summarise, beep = TRUE,
  parallel = TRUE, ncores = 20, save = TRUE, save_results = TRUE,
  progress = TRUE, filename = "./svp_01"
)
res

svp_01_process.R

Sys.time()

folder <- "~/Dropbox/small_var_prior/"

styler::style_file(paste0(folder, "svp_01_process.R"))

library(ggplot2)
theme_set(theme_bw())
library(scales)
library(data.table)
library(directlabels)

cb_palette <- c(
  "#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00",
  "#CC79A7"
)
cbb_palette <- c(
  "#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00",
  "#CC79A7"
)

# import simulation results ----
results_folder <- paste0(folder, "svp_01-results_COE-STD-149221/")
sim_results <- rbindlist(lapply(
  list.files(results_folder), function(file_name) {
    dat <- as.data.frame(readRDS(
      paste0(results_folder, file_name)
    )$results)
    design <- as.integer(regmatches(file_name, regexpr("\\d+", file_name)))
    dat$design <- design
    dat
  }
), fill = TRUE, use.names = TRUE)
sim_results[, iter := seq_len(.N)]
sim_results[]
sim_results[, .N, list(design, iter)][order(design)]
sim_results[]
colnames(sim_results)

sim_results <- melt.data.table(sim_results, id.vars = c("design", "iter"))
sim_results[]

design <- SimDesign::createDesign(
  lambda_32 = c(0, .2, .5, .7),
  n = c(seq(50, 150, 25), seq(200, 300, 50), c(400, 500, 1000)),
  group_iter = 1:10
)
design$design <- seq_len(nrow(design))
sim_results <- merge(
  sim_results, design[, which(colnames(design) != "group_iter")],
  all.x = TRUE
)
sim_results[]
sim_results <- sim_results[complete.cases(sim_results)]
sim_results[]

# clean up simulation results ----
str(sim_results)
sim_results[, variable_t := as.character(variable)]
sim_results[, .N, list(variable_t, variable)]
sim_results[
  , quantity := substr(
    variable_t, nchar(variable_t) - 1, nchar(variable_t)
  )
]
sim_results[, .N, quantity]
sim_results[, param := substr(variable_t, 1, 3)]
sim_results[, .N, param][order(param)]
sim_results[, pop := fcase(
  param %in% paste0("l", c(1:2, 4:6), c(rep(2, 2), rep(1, 3))), 0,
  param %in% c("l32"), lambda_32,
  param %in% c("p21"), .25,
  default = NA
)]
sim_results[]
sim_results[
  , rh_check := value[quantity == "rh"], list(design, iter, param)
]
sim_results[]
sim_results[, mean(rh_check > 1.1), list(param, lambda_32)]

# create and plot metrics ----
mo_dt <- sim_results[
  quantity %in% paste0("m", 1:4) | param == "p21" & quantity == "mn"
]
mo_dt[, mean(rh_check > 1.1)]
mo_dt <- mo_dt[rh_check < 1.1]
mo_dt[]

rb_dt <- mo_dt[
  ,
  .(est = mean(value), est_95 = quantile(abs(value), .95)),
  list(
    lambda_32, n,
    param_i = ifelse(
      param == "l32", 1, ifelse(param %in% c("l12", "l22"), 2,
        ifelse(param == "p21", 4, 3)
      )
    ),
    pop, quantity
  )
]
rb_dt[]
rb_dt[, param_f := factor(param_i, 1:4, c("L3,2", "L1:2,2", "L4:6,1", "P2,1"))]
rb_dt[, lambda_32_f := factor(
  lambda_32, c(0, .2, .5, .7), paste0("L3,2=", c(0, .2, .5, .7))
)]
rb_dt[]

ggplot(
  rb_dt[quantity == "m3" | param_i == 4],
  aes(n, est, group = param_f, shape = param_f, linetype = param_f)
) +
  geom_line(linewidth = .4) +
  geom_hline(yintercept = .01, linetype = 1, linewidth = .05) +
  geom_hline(
    aes(yintercept = pop * .9),
    linetype = 1, linewidth = .05,
    data = rb_dt[quantity == "m3" & param_i == 1]
  ) +
  geom_hline(
    aes(yintercept = pop * 1.1),
    linetype = 1, linewidth = .05,
    data = rb_dt[quantity == "m3" & param_i == 1]
  ) +
  annotate(
    "text",
    x = 140, y = .01, label = .01, hjust = 1, vjust = .02, size = 2, alpha = .7
  ) +
  scale_linetype_manual(values = c(1, 3, 5, 6)) +
  scale_x_continuous(breaks = c(50, 300, 500, 1000)) +
  scale_y_continuous(breaks = c(0, .2, .25, .3, .5, .7)) +
  scale_alpha(range = c(.25, 1)) +
  facet_grid(. ~ lambda_32_f, scales = "free") +
  guides(alpha = "none") +
  theme_bw() +
  theme(
    legend.position = "top",
    panel.grid.major.x = element_line(linewidth = .05, colour = 1),
    panel.grid.major.y = element_line(linewidth = .05, colour = 1),
    panel.grid.minor.x = element_blank(),
    panel.grid.minor.y = element_blank(),
    panel.spacing.x = unit(.5, "cm"),
    strip.background = element_blank(),
    panel.border = element_blank(),
    axis.ticks.y = element_blank(),
    axis.line.x = element_line(linewidth = .25),
    axis.title.y = element_blank()
  ) +
  labs(
    title = "Mean of posterior mode",
    linetype = "Which parameter(s)?", x = "Sample size"
  )
ggsave("svp_01_01_post_mode.png", width = 6.5, height = 3.5)

l32_first <- mo_dt[
  param != "p21",
  .(sum(param[which.max(abs(value))] == "l32")),
  list(lambda_32, n, quantity, iter)
][, .(count = sum(V1), .N), list(lambda_32, n)]
l32_first[]
l32_first[, lambda_32_f := factor(
  lambda_32, c(0, .2, .5, .7), paste0("L3,2=", c(0, .2, .5, .7))
)]

ggplot(
  l32_first,
  aes(
    n, count / N,
    shape = lambda_32_f, linetype = lambda_32_f
  )
) +
  scale_linetype_manual(values = c(3, 2, 5, 1)) +
  scale_shape_manual(values = as.character(c(0, 2, 5, 7))) +
  geom_point(size = 2) +
  geom_line(linewidth = .3) +
  geom_dl(aes(label = lambda_32_f), method = list("first.points", cex = .8)) +
  scale_x_continuous(
    breaks = c(50, 150, 300, 500, 1000), limits = c(-50, 1000)
  ) +
  scale_y_continuous(labels = percent_format()) +
  theme(
    legend.position = "top",
    panel.grid.major.x = element_line(linewidth = .05, colour = 1),
    panel.grid.major.y = element_line(linewidth = .05, colour = 1),
    panel.grid.minor.x = element_blank(),
    panel.spacing.x = unit(.5, "cm"),
    strip.background = element_blank(),
    panel.border = element_blank(),
    axis.ticks.y = element_blank(),
    axis.line.x = element_line(linewidth = .25),
    axis.title.y = element_blank()
  ) +
  guides(shape = "none", linetype = "none") +
  labs(
    title = "How often is L3,2 selected as the largest loading?",
    linetype = "", shape = "", x = "Sample size"
  )
ggsave("svp_01_02_largest_mode.png", width = 6.5, height = 3.5)

max_tab <- dcast.data.table(mo_dt[
  quantity == "m3",
  .(max(abs(value))),
  list(lambda_32, n, quantity, iter, param == "l32")
], lambda_32 + n + iter ~ param, value.var = "V1")
max_tab[]
max_tab[, gap := `TRUE` - `FALSE`]
max_tab[]
max_tab <- max_tab[, .(mean_gap = mean(gap, na.rm = TRUE)), list(lambda_32, n)]
max_tab[]
max_tab[, lambda_32_f := factor(
  lambda_32, c(0, .2, .5, .7), paste0("L3,2=", c(0, .2, .5, .7))
)]

ggplot(
  max_tab,
  aes(
    n, mean_gap,
    shape = lambda_32_f, linetype = lambda_32_f
  )
) +
  scale_linetype_manual(values = c(3, 2, 5, 1)) +
  scale_shape_manual(values = as.character(c(0, 2, 5, 7))) +
  geom_point(size = 2) +
  geom_line(linewidth = .3) +
  geom_dl(aes(label = lambda_32_f), method = list("last.points", cex = .8)) +
  scale_x_continuous(
    breaks = c(50, 150, 300, 500, 1000), limits = c(50, 1100)
  ) +
  scale_y_continuous(breaks = c(0, .2, .5, .7)) +
  theme_bw() +
  theme(
    legend.position = "top",
    panel.grid.major.x = element_line(linewidth = .05, colour = 1),
    panel.grid.major.y = element_line(linewidth = .05, colour = 1),
    panel.grid.minor.x = element_blank(),
    panel.spacing.x = unit(.5, "cm"),
    strip.background = element_blank(),
    panel.border = element_blank(),
    axis.ticks.y = element_blank(),
    axis.line.x = element_line(linewidth = .25),
    axis.title.y = element_blank()
  ) +
  guides(shape = "none", linetype = "none") +
  labs(
    title = "Average gap between L3,2 and the largest of remaining loadings",
    linetype = "", x = "Sample size"
  )
ggsave("svp_01_03_mean_mode_gap.png", width = 6.5, height = 3.5)