PWC 2-Dimensional Model CDAT
Steve Heil
2025-06-15
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.247 0.401 0.016 0.017 0.559 0.955 1.244 1.523
## theta1.1154 -0.019 0.373 0.015 0.018 -0.713 -0.257 -0.030 0.223
## theta1.1222 0.604 0.354 0.014 0.015 -0.107 0.358 0.616 0.859
## theta1.1637 0.606 0.362 0.015 0.016 -0.107 0.362 0.608 0.847
## theta1.1995 0.849 0.406 0.017 0.020 0.088 0.559 0.845 1.120
## theta1.3369 1.686 0.441 0.018 0.021 0.901 1.360 1.682 1.977
## theta1.4004 -0.821 0.395 0.016 0.016 -1.620 -1.097 -0.810 -0.549
## theta1.4209 -1.271 0.407 0.017 0.019 -2.023 -1.547 -1.261 -1.010
## theta1.4268 -0.850 0.396 0.016 0.018 -1.692 -1.101 -0.812 -0.558
## theta1.4841 -0.954 0.374 0.015 0.015 -1.684 -1.194 -0.949 -0.713
## theta1.5074 1.134 0.389 0.016 0.017 0.431 0.857 1.137 1.391
## theta1.5495 -2.323 0.427 0.017 0.017 -3.159 -2.598 -2.317 -2.044
## theta1.5710 -0.952 0.386 0.016 0.016 -1.733 -1.235 -0.943 -0.724
## theta1.8294 -0.138 0.356 0.015 0.016 -0.851 -0.374 -0.136 0.104
## theta1.8980 0.685 0.362 0.015 0.015 -0.006 0.443 0.702 0.935
## theta1.9899 0.419 0.360 0.015 0.016 -0.280 0.186 0.417 0.662
## theta2.1121 -1.168 0.364 0.015 0.016 -1.933 -1.396 -1.153 -0.918
## theta2.1154 -1.536 0.351 0.014 0.017 -2.231 -1.778 -1.531 -1.309
## theta2.1222 -0.160 0.333 0.014 0.015 -0.833 -0.376 -0.168 0.061
## theta2.1637 -0.514 0.322 0.013 0.014 -1.140 -0.735 -0.507 -0.302
## theta2.1995 -1.265 0.351 0.014 0.015 -1.910 -1.528 -1.254 -1.032
## theta2.3369 1.246 0.374 0.015 0.018 0.486 1.014 1.250 1.492
## theta2.4004 1.380 0.333 0.014 0.015 0.755 1.142 1.381 1.581
## theta2.4209 0.779 0.341 0.014 0.014 0.203 0.523 0.754 1.015
## theta2.4268 1.206 0.350 0.014 0.016 0.554 0.970 1.205 1.430
## theta2.4841 -0.250 0.324 0.013 0.013 -0.875 -0.472 -0.255 -0.048
## theta2.5074 -1.068 0.342 0.014 0.017 -1.719 -1.289 -1.067 -0.847
## theta2.5495 0.575 0.333 0.014 0.015 0.073 0.326 0.542 0.793
## theta2.5710 -0.834 0.328 0.013 0.015 -1.516 -1.037 -0.837 -0.637
## theta2.8294 0.728 0.318 0.013 0.014 0.148 0.509 0.727 0.939
## theta2.8980 0.893 0.343 0.014 0.018 0.235 0.650 0.922 1.124
## theta2.9899 0.665 0.331 0.014 0.014 -0.032 0.464 0.666 0.892
## gamma.1 0.872 0.126 0.005 0.006 0.622 0.794 0.872 0.964
## gamma.25 0.725 0.140 0.006 0.008 0.438 0.633 0.726 0.829
## gamma.26 0.830 0.121 0.005 0.006 0.588 0.747 0.828 0.917
## gamma.27 0.837 0.121 0.005 0.006 0.591 0.756 0.837 0.928
## gamma.28 0.454 0.135 0.005 0.007 0.197 0.367 0.457 0.544
## gamma.29 0.917 0.137 0.006 0.007 0.648 0.826 0.920 1.009
## gamma.30 0.317 0.150 0.006 0.007 0.032 0.209 0.317 0.418
## gamma.31 0.076 0.071 0.003 0.003 0.002 0.024 0.050 0.112
## gamma.34 0.944 0.143 0.006 0.007 0.677 0.849 0.942 1.041
## gamma.35 1.341 0.136 0.006 0.005 1.050 1.258 1.349 1.444
## gamma.36 1.313 0.128 0.005 0.006 1.049 1.235 1.325 1.407
## gamma.37 1.326 0.116 0.005 0.005 1.068 1.250 1.332 1.402
## gamma.38 1.450 0.086 0.004 0.004 1.243 1.405 1.465 1.517
## 97.5% Model
## theta1.1121 2.095 2d
## theta1.1154 0.758 2d
## theta1.1222 1.212 2d
## theta1.1637 1.294 2d
## theta1.1995 1.621 2d
## theta1.3369 2.587 2d
## theta1.4004 -0.108 2d
## theta1.4209 -0.468 2d
## theta1.4268 -0.180 2d
## theta1.4841 -0.245 2d
## theta1.5074 1.896 2d
## theta1.5495 -1.540 2d
## theta1.5710 -0.175 2d
## theta1.8294 0.514 2d
## theta1.8980 1.369 2d
## theta1.9899 1.102 2d
## theta2.1121 -0.481 2d
## theta2.1154 -0.823 2d
## theta2.1222 0.486 2d
## theta2.1637 0.141 2d
## theta2.1995 -0.595 2d
## theta2.3369 1.964 2d
## theta2.4004 2.030 2d
## theta2.4209 1.443 2d
## theta2.4268 1.912 2d
## theta2.4841 0.404 2d
## theta2.5074 -0.348 2d
## theta2.5495 1.325 2d
## theta2.5710 -0.165 2d
## theta2.8294 1.367 2d
## theta2.8980 1.565 2d
## theta2.9899 1.278 2d
## gamma.1 1.104 2d
## gamma.25 0.987 2d
## gamma.26 1.062 2d
## gamma.27 1.058 2d
## gamma.28 0.722 2d
## gamma.29 1.191 2d
## gamma.30 0.609 2d
## gamma.31 0.261 2d
## gamma.34 1.203 2d
## gamma.35 1.552 2d
## gamma.36 1.538 2d
## gamma.37 1.538 2d
## gamma.38 1.564 2dPlot 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) %>%
select(Model, Variable, ID, everything() )
summary2dWideThetas <- summary2d %>%
filter(Variable == "theta1" | Variable == "theta2") %>%
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 = "Posterior Medians of and 95% Credible Intervals",
x = "Theta 1",
y = "Theta 2"
) +
theme_minimal()
Plot the median gamma for each rater
gammaAngles <- summary2d %>%
filter(Variable == "gamma") %>%
select(ID, gamma = `50%`) %>%
mutate(
x = cos(gamma) * 2, # x component
y = sin(gamma)* 2 # 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 = "Rater Preference Vectors (Gamma Angles)",
x = "Theta 1",
y = "Theta 2"
) +
# 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()