2-Dimensional Model with Clusters

library(bayesplot)
library(coda)
library(ggplot2)
library(ggrepel)
library(tidyverse)

Load 2-dimensional DP prior model MCMC object

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)

Summary of MCMC chain output

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

Plot 2dDP Model posterior parameters from samples

Trace Plots

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

Density Plots

densplot(mcmc_list2dDP, right = FALSE) 

# Right = FALSE helps makes each cluster integer fall neatly to the right of its own bin in histogram

Autocorrelation Plots

Next the autocorrelation plots

autocorr.plot(mcmc_list2dDP, auto.layout = FALSE, lag.max=20)

Create a tidy data table

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%'))

Plot the posterior means of items

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