# 心理学データ解析応用/伴走サイトコード -----------------------------------------------------
#  Programmed by kosugitti
#  Licence ; Creative Commons BY-SA license (CC BY-SA) version 4.0
## Lesson 7. Alternatives for ANOVA #2. Rstan Version

# 準備 ----------------------------------------------------------------------

rm(list = ls())
library(tidyverse)
library(bayesplot)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library(MASS)
## MAP関数
map_estimation <- function(z) {
  density(z)$x[which.max(density(z)$y)]
}
## MCMCサンプルをデータフレームにする関数
MCMCtoDF <- function(fit) {
  fit %>%
    rstan::extract() %>%
    as.data.frame() %>%
    tibble::as_tibble() %>%
    tibble::rowid_to_column("iter") %>%
    dplyr::select(-lp__) %>%
    tidyr::pivot_longer(-iter) -> MCMCsample
  return(MCMCsample)
}

## MCMCデータフレームを要約する関数
MCMCsummary <- function(MCMCsample) {
  MCMCsample %>%
    dplyr::group_by(name) %>%
    dplyr::summarise(
      EAP = mean(value),
      MED = median(value),
      MAP = map_estimation(value),
      SD = sd(value),
      L95 = quantile(value, prob = 0.025),
      U95 = quantile(value, prob = 0.975),
      Rhat = posterior::rhat_basic(value)
    ) %>%
    mutate(across(where(is.numeric), ~ num(., digits = 3)))
}


# リバース・リバースエンジニアリング -------------------------------------------------------

mu <- c(50, 50)
sd1 <- 10
sd2 <- 5
rho <- 0.7
SIG <- matrix(ncol = 2, nrow = 2)
SIG[1, 1] <- sd1 * sd1
SIG[2, 2] <- sd2 * sd2
SIG[1, 2] <- sd1 * sd2 * rho
SIG[2, 1] <- sd1 * sd2 * rho
N <- 100
X <- mvrnorm(N, mu, SIG)

dataSet <- list(N = N, X = X)

model <- rstan::stan_model("rstan/paired_t.stan")
fit <- rstan::sampling(model,
  data = dataSet,
  chains = 4,
  iter = 6000,
  warmup = 1000
)

fit %>%
  MCMCtoDF() %>%
  MCMCsummary()


# tidyなデータ ----------------------------------------------------------------

X <- mvrnorm(N, mu, SIG)
tidy_data <- X %>%
  as.data.frame() %>%
  rename(pre = V1, post = V2) %>%
  as_tibble() %>%
  rowid_to_column("ID") %>%
  pivot_longer(-ID) %>%
  mutate(cond = if_else(name == "pre", 1, 2))

dataSet <- list(
  L = NROW(tidy_data),
  N = max(tidy_data$ID),
  IDindex = tidy_data$ID,
  Condition = tidy_data$cond,
  val = tidy_data$value
)

model <- rstan::stan_model("rstan/paired_t_tidy.stan")
fit <- rstan::sampling(model,
  data = dataSet,
  chains = 4,
  iter = 6000,
  warmup = 1000
)

fit %>%
  MCMCtoDF() %>%
  MCMCsummary()


# 課題 --------------------------------------------------------------------

dat_raw <- data.frame(
  ID = 1:4,
  period1 = c(10, 9, 4, 7),
  period2 = c(5, 4, 2, 3),
  period3 = c(9, 5, 3, 5)
)

tidy_dat <- dat_raw %>%
  pivot_longer(-ID) %>%
  mutate(name = as.factor(name) %>% fct_relevel("period1", "period2")) %>%
  mutate(cond = as.numeric(name))