--- 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) ```