library(bayesplot)
library(coda)
library(ggplot2)
library(ggrepel)
library(tidyverse)
load("Out/mcmc_list2d.Rdata")
# Create a list of parameters for each type of parameter
thetas <- grep("^theta", colnames(mcmc_list2d[[1]]), value = TRUE)
gammas <- grep("^gamma", colnames(mcmc_list2d[[1]]), value = TRUE)
View a summary of model developed here: https://rpubs.com/steveneheil/pwc-MCMC-CDAT
load("Out/summary2d.Rdata")
summary2d
## Mean SD Naive SE Time-series SE 2.5% 25% 50% 75%
## theta1.1121 1.243 0.401 0.007 0.011 0.475 0.969 1.241 1.519
## theta1.1154 0.007 0.387 0.007 0.010 -0.743 -0.258 0.005 0.265
## theta1.1222 0.613 0.362 0.007 0.009 -0.091 0.366 0.613 0.859
## theta1.1637 0.561 0.362 0.007 0.010 -0.118 0.309 0.555 0.805
## theta1.1995 0.872 0.407 0.007 0.011 0.101 0.594 0.860 1.150
## theta1.3369 1.631 0.463 0.008 0.010 0.743 1.317 1.620 1.936
## theta1.4004 -0.854 0.399 0.007 0.010 -1.671 -1.118 -0.861 -0.584
## theta1.4209 -1.234 0.391 0.007 0.010 -2.025 -1.491 -1.231 -0.965
## theta1.4268 -0.894 0.398 0.007 0.010 -1.697 -1.158 -0.887 -0.619
## theta1.4841 -0.912 0.383 0.007 0.010 -1.671 -1.167 -0.910 -0.645
## theta1.5074 1.154 0.386 0.007 0.010 0.429 0.880 1.150 1.420
## theta1.5495 -2.243 0.446 0.008 0.012 -3.162 -2.538 -2.234 -1.930
## theta1.5710 -0.873 0.395 0.007 0.009 -1.644 -1.140 -0.862 -0.606
## theta1.8294 -0.116 0.369 0.007 0.008 -0.844 -0.366 -0.115 0.137
## theta1.8980 0.679 0.384 0.007 0.009 -0.062 0.421 0.674 0.942
## theta1.9899 0.357 0.376 0.007 0.009 -0.410 0.098 0.363 0.624
## theta2.1121 -1.224 0.365 0.007 0.010 -1.940 -1.465 -1.217 -0.971
## theta2.1154 -1.546 0.359 0.007 0.009 -2.264 -1.781 -1.549 -1.305
## theta2.1222 -0.160 0.328 0.006 0.009 -0.788 -0.379 -0.154 0.068
## theta2.1637 -0.565 0.330 0.006 0.009 -1.220 -0.780 -0.565 -0.335
## theta2.1995 -1.239 0.363 0.007 0.010 -1.969 -1.476 -1.239 -0.998
## theta2.3369 1.183 0.394 0.007 0.010 0.389 0.919 1.182 1.454
## theta2.4004 1.316 0.354 0.006 0.009 0.660 1.069 1.299 1.554
## theta2.4209 0.764 0.334 0.006 0.009 0.132 0.532 0.758 0.980
## theta2.4268 1.122 0.345 0.006 0.008 0.470 0.888 1.114 1.339
## theta2.4841 -0.262 0.329 0.006 0.008 -0.892 -0.475 -0.270 -0.057
## theta2.5074 -1.056 0.358 0.007 0.010 -1.801 -1.289 -1.051 -0.806
## theta2.5495 0.579 0.332 0.006 0.009 0.046 0.327 0.549 0.793
## theta2.5710 -0.845 0.345 0.006 0.009 -1.505 -1.071 -0.848 -0.622
## theta2.8294 0.753 0.324 0.006 0.009 0.137 0.537 0.755 0.971
## theta2.8980 0.874 0.346 0.006 0.008 0.194 0.648 0.863 1.103
## theta2.9899 0.581 0.342 0.006 0.009 -0.106 0.354 0.589 0.806
## gamma.25 0.729 0.137 0.002 0.004 0.450 0.636 0.731 0.820
## gamma.26 0.832 0.127 0.002 0.004 0.582 0.747 0.833 0.918
## gamma.27 0.842 0.126 0.002 0.004 0.594 0.755 0.845 0.928
## gamma.28 0.463 0.140 0.003 0.004 0.190 0.365 0.463 0.562
## gamma.29 0.917 0.142 0.003 0.004 0.639 0.821 0.916 1.011
## gamma.30 0.327 0.148 0.003 0.004 0.054 0.221 0.321 0.426
## gamma.31 0.086 0.072 0.001 0.002 0.002 0.030 0.068 0.124
## gamma.34 0.934 0.147 0.003 0.005 0.629 0.839 0.938 1.032
## gamma.35 1.327 0.141 0.003 0.003 1.023 1.236 1.339 1.432
## gamma.36 1.309 0.125 0.002 0.003 1.048 1.227 1.312 1.400
## gamma.37 1.323 0.119 0.002 0.003 1.075 1.247 1.327 1.409
## gamma.38 1.442 0.091 0.002 0.002 1.236 1.386 1.457 1.514
## 97.5% Model
## theta1.1121 2.023 2d
## theta1.1154 0.785 2d
## theta1.1222 1.312 2d
## theta1.1637 1.255 2d
## theta1.1995 1.675 2d
## theta1.3369 2.558 2d
## theta1.4004 -0.082 2d
## theta1.4209 -0.472 2d
## theta1.4268 -0.152 2d
## theta1.4841 -0.185 2d
## theta1.5074 1.917 2d
## theta1.5495 -1.398 2d
## theta1.5710 -0.116 2d
## theta1.8294 0.596 2d
## theta1.8980 1.443 2d
## theta1.9899 1.067 2d
## theta2.1121 -0.533 2d
## theta2.1154 -0.845 2d
## theta2.1222 0.458 2d
## theta2.1637 0.060 2d
## theta2.1995 -0.527 2d
## theta2.3369 1.932 2d
## theta2.4004 2.034 2d
## theta2.4209 1.427 2d
## theta2.4268 1.844 2d
## theta2.4841 0.399 2d
## theta2.5074 -0.365 2d
## theta2.5495 1.333 2d
## theta2.5710 -0.164 2d
## theta2.8294 1.396 2d
## theta2.8980 1.570 2d
## theta2.9899 1.256 2d
## gamma.25 0.992 2d
## gamma.26 1.081 2d
## gamma.27 1.083 2d
## gamma.28 0.728 2d
## gamma.29 1.193 2d
## gamma.30 0.631 2d
## gamma.31 0.268 2d
## gamma.34 1.215 2d
## gamma.35 1.551 2d
## gamma.36 1.534 2d
## gamma.37 1.531 2d
## gamma.38 1.564 2d
First traceplots of theta parameter iterations. This bayesplot function has a default fixed scale for x axis, which enables better comparison between parameters than the traceplot function in coda. So fixed scales work best with same parameters.
Start with thetas
# Trace plot only for theta parameters
mcmc_trace(mcmc_list2d, pars = thetas, facet_args = list(ncol = 2, scales = "fixed"))
Next the gammas with their own fixed y axis.
# Trace plot only for gamma parameters
mcmc_trace(mcmc_list2d, pars = gammas, facet_args = list(ncol = 2, scales = "fixed"))
Next density plots of all posterior parameter distributions
densplot(mcmc_list2d)
Next the autocorrelation plots
autocorr.plot(mcmc_list2d, auto.layout = FALSE, lag.max=20)
### Create a tidy data table
summary2d <- summary2d %>%
mutate(
Variable = str_extract(rownames(.), "^[^\\.]+"),
ID = str_extract(rownames(.), "(?<=\\.)\\d+")
) %>%
`rownames<-`(NULL) %>%
dplyr::select(Model, Variable, ID, everything() )
summary2dWideThetas <- summary2d %>%
filter(Variable == "theta1" | Variable == "theta2") %>%
dplyr::select(ID, Variable, `50%`, `2.5%`, `97.5%`) %>%
pivot_wider(names_from = Variable, values_from = c('50%', '2.5%', '97.5%'))
ggplot(summary2dWideThetas, aes(x = `50%_theta1`, y = `50%_theta2`, label = ID)) +
# Add horizontal error bars (theta1 CI)
geom_errorbarh(aes(xmin = `2.5%_theta1`, xmax = `97.5%_theta1`), height = 0.05, color = "gray60", linewidth = 0.3) +
# Add vertical error bars (theta2 CI)
geom_errorbar(aes(ymin = `2.5%_theta2`, ymax = `97.5%_theta2`), width = 0.05, color = "gray60", linewidth = 0.3) +
# Add point and label
geom_point() +
geom_text(vjust = -0.5, size = 3) +
labs(
title = "2D Model Posterior Medians and 95% Credible Intervals",
x = "Theta 1 Creative Elaboration & Fluency",
y = "Theta 2 Skill with Pictorial Convention"
) +
theme_minimal()
Plot the median gamma for each rater
gammaAngles <- summary2d %>%
filter(Variable == "gamma") %>%
dplyr::select(ID, gamma = `50%`) %>%
mutate(
x = cos(gamma), # x component
y = sin(gamma) # y component
)
ggplot(gammaAngles, aes(x = 0, y = 0, xend = x, yend = y)) +
geom_segment(arrow = arrow(length = unit(0.15, "cm")),
color = "black", linewidth = 0.7) +
coord_fixed() + # 1:1 aspect ratio
labs(
title = "2D Model Median Gamma Rater Preference",
x = "Theta 1 Creative Elaboration & Fluency",
y = "Theta 2 Skill with Pictorial Convention"
) +
# Smart label positioning
geom_text_repel(
data = gammaAngles,
aes(x = x, y = y, label = ID),
size = 3,
segment.color = NA,
nudge_x = 0.1, nudge_y = 0.1,
max.overlaps = Inf) +
theme_minimal()