library(bayesplot)
library(coda)
library(ggplot2)
library(ggrepel)
library(tidyverse)
load("Out/mcmc_list2dDP.Rdata")
thetas <- grep("^theta", colnames(mcmc_list2dDP[[1]]), value = TRUE)
gammas <- grep("^gamma", colnames(mcmc_list2dDP[[1]]), value = TRUE)
clusters <- grep("^cluster", colnames(mcmc_list2dDP[[1]]), value = TRUE)
View a summary of model developed here: https://rpubs.com/steveneheil/pwc-MCMC-CDAT
load("Out/summary2dDP.Rdata")
summary2dDP
## Mean SD Naive SE Time-series SE 2.5% 25% 50% 75%
## theta1.1121 1.227 0.422 0.005 0.008 0.442 0.938 1.212 1.503
## theta1.1154 0.053 0.406 0.005 0.007 -0.724 -0.233 0.049 0.321
## theta1.1222 0.676 0.361 0.005 0.006 0.000 0.431 0.665 0.907
## theta1.1637 0.590 0.364 0.005 0.006 -0.111 0.342 0.581 0.832
## theta1.1995 0.862 0.433 0.006 0.008 0.060 0.562 0.843 1.142
## theta1.3369 1.539 0.451 0.006 0.007 0.654 1.234 1.538 1.839
## theta1.4004 -0.857 0.412 0.005 0.007 -1.698 -1.132 -0.847 -0.571
## theta1.4209 -1.300 0.415 0.005 0.007 -2.142 -1.575 -1.285 -1.015
## theta1.4268 -0.902 0.431 0.006 0.007 -1.769 -1.188 -0.885 -0.598
## theta1.4841 -0.894 0.387 0.005 0.006 -1.656 -1.152 -0.889 -0.630
## theta1.5074 1.222 0.404 0.005 0.007 0.468 0.943 1.208 1.483
## theta1.5495 -2.279 0.465 0.006 0.008 -3.229 -2.590 -2.263 -1.952
## theta1.5710 -0.807 0.400 0.005 0.006 -1.569 -1.081 -0.805 -0.538
## theta1.8294 -0.154 0.368 0.005 0.006 -0.884 -0.397 -0.152 0.095
## theta1.8980 0.655 0.389 0.005 0.006 -0.124 0.403 0.652 0.914
## theta1.9899 0.385 0.381 0.005 0.006 -0.368 0.134 0.380 0.643
## theta2.1121 -1.119 0.360 0.005 0.007 -1.848 -1.357 -1.107 -0.874
## theta2.1154 -1.546 0.360 0.005 0.007 -2.261 -1.783 -1.538 -1.305
## theta2.1222 -0.188 0.331 0.004 0.007 -0.852 -0.406 -0.185 0.039
## theta2.1637 -0.558 0.330 0.004 0.006 -1.221 -0.771 -0.561 -0.334
## theta2.1995 -1.173 0.356 0.005 0.007 -1.892 -1.407 -1.164 -0.930
## theta2.3369 1.270 0.401 0.005 0.007 0.468 1.005 1.272 1.538
## theta2.4004 1.215 0.356 0.005 0.007 0.538 0.970 1.208 1.451
## theta2.4209 0.713 0.357 0.005 0.007 0.021 0.471 0.709 0.942
## theta2.4268 1.028 0.354 0.005 0.006 0.345 0.786 1.023 1.268
## theta2.4841 -0.324 0.341 0.004 0.007 -0.994 -0.553 -0.325 -0.099
## theta2.5074 -1.035 0.365 0.005 0.007 -1.788 -1.276 -1.024 -0.781
## theta2.5495 0.455 0.409 0.005 0.008 -0.274 0.170 0.439 0.719
## theta2.5710 -0.936 0.357 0.005 0.007 -1.629 -1.174 -0.939 -0.702
## theta2.8294 0.747 0.333 0.004 0.007 0.075 0.523 0.747 0.966
## theta2.8980 0.885 0.341 0.004 0.007 0.208 0.659 0.887 1.110
## theta2.9899 0.549 0.336 0.004 0.006 -0.115 0.322 0.558 0.774
## gamma.25 0.859 0.153 0.002 0.003 0.527 0.760 0.873 0.970
## gamma.26 0.887 0.137 0.002 0.003 0.601 0.795 0.894 0.986
## gamma.27 0.888 0.137 0.002 0.003 0.602 0.795 0.895 0.987
## gamma.28 0.455 0.200 0.003 0.004 0.050 0.329 0.465 0.590
## gamma.29 0.897 0.139 0.002 0.003 0.612 0.803 0.904 0.996
## gamma.30 0.349 0.207 0.003 0.004 0.018 0.166 0.358 0.508
## gamma.31 0.129 0.100 0.001 0.001 0.005 0.050 0.107 0.183
## gamma.34 0.900 0.144 0.002 0.003 0.614 0.803 0.904 0.997
## gamma.35 1.424 0.115 0.001 0.002 1.137 1.365 1.448 1.514
## gamma.36 1.423 0.112 0.001 0.002 1.146 1.361 1.446 1.511
## gamma.37 1.423 0.113 0.001 0.002 1.138 1.361 1.447 1.512
## gamma.38 1.437 0.102 0.001 0.002 1.186 1.380 1.457 1.518
## cluster.25 2.514 1.516 0.020 0.037 1.000 1.000 2.000 3.000
## cluster.26 2.432 1.470 0.019 0.039 1.000 1.000 2.000 3.000
## cluster.27 2.450 1.480 0.019 0.038 1.000 1.000 2.000 3.000
## cluster.28 3.299 1.593 0.021 0.027 1.000 2.000 3.000 5.000
## cluster.29 2.475 1.491 0.019 0.037 1.000 1.000 2.000 3.000
## cluster.30 3.258 1.570 0.020 0.029 1.000 2.000 3.000 4.000
## cluster.31 3.579 1.604 0.021 0.030 1.000 2.000 4.000 5.000
## cluster.34 2.472 1.493 0.019 0.038 1.000 1.000 2.000 3.000
## cluster.35 2.675 1.501 0.019 0.039 1.000 1.000 2.000 4.000
## cluster.36 2.670 1.499 0.019 0.039 1.000 1.000 2.000 4.000
## cluster.37 2.673 1.502 0.019 0.037 1.000 1.000 2.000 4.000
## cluster.38 2.680 1.506 0.019 0.038 1.000 1.000 2.000 4.000
## n.clusters 4.476 0.782 0.010 0.012 3.000 4.000 4.000 5.000
## alpha 1.060 0.605 0.008 0.009 0.247 0.618 0.936 1.372
## 97.5% Model
## theta1.1121 2.093 2dDP
## theta1.1154 0.882 2dDP
## theta1.1222 1.402 2dDP
## theta1.1637 1.316 2dDP
## theta1.1995 1.746 2dDP
## theta1.3369 2.422 2dDP
## theta1.4004 -0.090 2dDP
## theta1.4209 -0.512 2dDP
## theta1.4268 -0.100 2dDP
## theta1.4841 -0.134 2dDP
## theta1.5074 2.048 2dDP
## theta1.5495 -1.404 2dDP
## theta1.5710 -0.014 2dDP
## theta1.8294 0.557 2dDP
## theta1.8980 1.417 2dDP
## theta1.9899 1.134 2dDP
## theta2.1121 -0.443 2dDP
## theta2.1154 -0.850 2dDP
## theta2.1222 0.459 2dDP
## theta2.1637 0.075 2dDP
## theta2.1995 -0.499 2dDP
## theta2.3369 2.038 2dDP
## theta2.4004 1.927 2dDP
## theta2.4209 1.443 2dDP
## theta2.4268 1.726 2dDP
## theta2.4841 0.354 2dDP
## theta2.5074 -0.352 2dDP
## theta2.5495 1.320 2dDP
## theta2.5710 -0.216 2dDP
## theta2.8294 1.408 2dDP
## theta2.8980 1.553 2dDP
## theta2.9899 1.213 2dDP
## gamma.25 1.115 2dDP
## gamma.26 1.130 2dDP
## gamma.27 1.131 2dDP
## gamma.28 0.842 2dDP
## gamma.29 1.143 2dDP
## gamma.30 0.732 2dDP
## gamma.31 0.372 2dDP
## gamma.34 1.158 2dDP
## gamma.35 1.564 2dDP
## gamma.36 1.565 2dDP
## gamma.37 1.564 2dDP
## gamma.38 1.565 2dDP
## cluster.25 6.000 2dDP
## cluster.26 6.000 2dDP
## cluster.27 6.000 2dDP
## cluster.28 6.000 2dDP
## cluster.29 6.000 2dDP
## cluster.30 6.000 2dDP
## cluster.31 6.000 2dDP
## cluster.34 6.000 2dDP
## cluster.35 6.000 2dDP
## cluster.36 6.000 2dDP
## cluster.37 6.000 2dDP
## cluster.38 6.000 2dDP
## n.clusters 6.000 2dDP
## alpha 2.560 2dDP
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_list2dDP, pars = thetas, facet_args = list(ncol = 2, scales = "fixed"))
Next the gammas with their own fixed y axis.
# Trace plot only for gamma
mcmc_trace(mcmc_list2dDP, pars = gammas, facet_args = list(ncol = 2, scales = "fixed"))
Start with finally the clusters with their own fixed y axis.
# Trace plot only for cluster parameters
mcmc_trace(mcmc_list2dDP, pars = clusters, facet_args = list(ncol = 2, scales = "fixed"))
Next density plots of all posterior parameter distributions
densplot(mcmc_list2dDP, right = FALSE)
# Right = FALSE helps makes each cluster integer fall neatly to the right of its own bin in histogram
Next the autocorrelation plots
autocorr.plot(mcmc_list2dDP, auto.layout = FALSE, lag.max=20)
summary2dDP.1 <- summary2dDP %>%
mutate(
Variable = str_extract(rownames(.), "^[^\\.]+"),
ID = str_extract(rownames(.), "(?<=\\.)\\d+")
) %>%
`rownames<-`(NULL) %>%
dplyr::select(Model, Variable, ID, everything() )
summary2dDPWideThetas <- summary2dDP.1 %>%
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(summary2dDPWideThetas, 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 Clustering Model Posterior Medians and 95% Credible Intervals",
x = "Theta 1 Creative Elaboration & Fluency",
y = "Theta 2 Skill with Pictorial Convention"
) +
theme_minimal()
Plot the gamma median for each rater.
gammaAngles <- summary2dDP.1 %>%
filter(Variable == "gamma") %>%
pivot_wider(
id_cols = "ID",
names_from = c("Variable"),
values_from = c("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 Clustering 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,
max.iter = 1000,
segment.color = NA,
nudge_x = 0.1, nudge_y = 0.1,
max.overlaps = Inf) +
theme_minimal()