## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE
)


## ----setup--------------------------------------------------------------------
library(tibble)
library(dplyr)
library(tidyr)
library(ggplot2)
library(RStanTVA)
library(brms)

## -----------------------------------------------------------------------------
rstan_options(threads_per_chain = parallel::detectCores())

## ----meanprobitk--------------------------------------------------------------
meanprobitK <- Vectorize(function(mK, sigmaK, nS) {
  p <- pnorm(1:nS, mK, sigmaK)
  lps <- c(0, p)
  ups <- c(p, 1)
  sum(0:nS * (ups - lps))
}, c("mK","sigmaK"))

## ----load-data----------------------------------------------------------------
data("tva_recovery")
head(tva_recovery)

## ----load-true-vals-----------------------------------------------------------
data("tva_recovery_true_params")
str(tva_recovery_true_params)

## ----narrow-true-vals---------------------------------------------------------
true_params <- tva_recovery_true_params$coef_subject

true_params_narrow <- true_params %>% mutate(`italic(w)[lambda]` = w[,2], `italic(K)` = meanprobitK(mK,sK,6))  %>% select(-w,-mK,-sK) %>% rename(`italic(t)[0]` = mu0, `sigma[0]` = sigma0, `italic(C)` = C) %>% pivot_longer(-c(subject,condition), names_to = "param", values_to = "true_value")
head(true_params_narrow)

true_params_narrow2 <- bind_rows(
  true_params_narrow %>% filter(condition == "low" & param != "italic(C)") %>% select(-condition),
  true_params_narrow %>% filter(param == "italic(C)") %>% transmute(subject, param = sprintf("%s[%s]", param, condition), true_value)
)
head(true_params_narrow2)

## ----fit_mle_nested_reg, include=FALSE----------------------------------------
fit_mle_nested_reg <- readRDS("tva_recovery_cache/fit_mle_nested_reg.rds")

## ----fit_bayes_fixed, include=FALSE-------------------------------------------
fit_bayesian_nested <- readRDS("tva_recovery_cache/fit_bayesian_nested.rds")

## ----fit_bayesian_hierarchical_globals, include=FALSE-------------------------
fit_bayesian_hierarchical_globals <- readRDS("tva_recovery_cache/fit_bayesian_hierarchical_globals.rds")

## ----plot-bayes, fig.width = 8, fig.height = 5, dpi = 100---------------------
figdat <- bind_rows(
  fit_mle_nested_reg %>% add_column(method = "MLRN"),
  fit_bayesian_hierarchical_globals %>% add_column(method = "HB"),
  fit_bayesian_nested %>% add_column(method = "NB")
) %>%
  mutate(method = factor(method, levels = c("MLRN","NB","HB")))

fig <- ggplot(figdat) +
  theme_minimal() +
  theme(text = element_text(family = "sans", size = 10)) +
  labs(x = "True value", y = "Recovered value") +
  facet_wrap(method~param, ncol = 7, scales = "free", labeller = function(x) label_parsed(list(sprintf("%s~(\"%s\")", x$param, x$method)))) +
  geom_abline(linetype = "dashed") +
  geom_linerange(aes(x=true_value,ymin=cri[,1],ymax=cri[,2]), color = "gray50", linewidth = 0.2) +
  geom_point(aes(x=true_value,y=fitted_value), size = 0.5, color = "blue") +
  geom_point(aes(x=vx,y=vy), color = "transparent", figdat %>% group_by(param) %>% reframe(vx = range(true_value), vy = range(c(fitted_value,cri), na.rm = TRUE)))

print(fig)