2-Dimensinoal Paired Comparison Model

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

Load 2-dimensional model MCMC object

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)

Summary of MCMC 2d with normal prior

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

Plot 2d 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_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

Density Plots

densplot(mcmc_list2d) 

Autocorrelation Plots

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

Plot the posterior means of items

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