---
title: "modelChecking"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
##-- Load packages --##
```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(bayesplot)
library(tidybayes)
library(rstan)
options(scipen = 50)
options(max.print = 1000000)
memory.limit(size = 24000)
```
# ==============================================================================
# M3a
# ==============================================================================
##-- Load the model --##
```{r}
load("C:/PhD project.Chapter3.Study1/Data/reinforcement learning models/kang.V4/m3a.Rdata")
```
##-- Check all parameters --##
```{r}
m3a@model_pars
# Parameters of interest: 'mu_xi', 'mu_ep', 'mu_b', 'mu_pi', 'mu_rho', 'sigma'
```
##-- Select parameters of interest and check their corresponding Rhat values and the effective sample size --##
```{r}
m3a.parametersOfInterest <- summary(m3a,
pars = c("mu_xi", "mu_ep", "mu_b",
"mu_pi", "mu_rho", 'sigma'),
probs = c(.025, .975))$summary
```
```{r}
print(m3a.parametersOfInterest)
```
##-- Check the traceplot --##
```{r}
rstan::stan_trace(m3a,
pars = c("mu_xi", "mu_ep", "mu_b",
"mu_pi", "mu_rho"),
inc_warmup = FALSE,
#nrow = 2,
ncol = 3)
```
```{r}
rstan::stan_trace(m3a,
pars = c('sigma'),
inc_warmup = FALSE,
#nrow = 2,
ncol = 3)
```
##-- trankplot --##
```{r message=FALSE}
rstan::extract(m3a,
pars = c("mu_xi", "mu_ep", "mu_b",
"mu_pi", "mu_rho"),
inc_warmup = FALSE, permuted = FALSE) %>%
bayesplot::mcmc_rank_overlay() +
scale_color_manual(values = c("#80A0C7", "#B1934A", "#A65141", "#EEDA9D")) +
scale_x_continuous(breaks = 0:16 * 1e3,
labels = c(0, str_c(1:16, "k"))) +
theme(legend.position = "bottom") +
facet_wrap(~parameter,
labeller = label_parsed,
# nrow = 2,
ncol = 3) +
coord_cartesian(ylim = c(100, NA))
```
```{r message=FALSE}
rstan::extract(m3a,
pars = c("sigma[1]", "sigma[2]", "sigma[3]",
"sigma[4]", "sigma[5]"),
inc_warmup = FALSE, permuted = FALSE) %>%
bayesplot::mcmc_rank_overlay() +
scale_color_manual(values = c("#80A0C7", "#B1934A", "#A65141", "#EEDA9D")) +
scale_x_continuous(breaks = 0:16 * 1e3,
labels = c(0, str_c(1:16, "k"))) +
theme(legend.position = "bottom") +
facet_wrap(~parameter,
labeller = label_parsed,
# nrow = 2,
ncol = 3) +
coord_cartesian(ylim = c(100, NA))
```
##-- Check the autocorrelation plot of the mcmc chains --##
```{r}
posterior.m3a <- as.array(m3a)
```
```{r}
bayesplot::mcmc_acf(posterior.m3a,
pars = c("mu_xi", "mu_ep", "mu_b",
"mu_pi", "mu_rho"),
lags = 20)
bayesplot::mcmc_acf(posterior.m3a,
pars = c("sigma[1]", "sigma[2]", "sigma[3]",
"sigma[4]", "sigma[5]"),
lags = 20)
```