PWC 1-Dimensional Model CDAT
Steve Heil
2025-06-29
1-Dimensional Paired Comparison Model
library(bayesplot)
library(coda)
library(ggplot2)
library(tidyverse)Load 1-dimensional model MCMC object
load("Out/mcmc_list1d.Rdata")
# Create a list of parameters for each type of parameter
thetas <- grep("^theta", colnames(mcmc_list1d[[1]]), value = TRUE)
alphas <- grep("^alpha", colnames(mcmc_list1d[[1]]), value = TRUE)Summary of MCMC 1d with normal prior
load("Out/summary1d.Rdata")
summary1d## Mean SD Naive SE Time-series SE 2.5% 25% 50% 75%
## theta.1121 -0.177 0.277 0.003 0.008 -0.730 -0.359 -0.173 0.010
## theta.1154 -0.947 0.343 0.004 0.009 -1.650 -1.172 -0.931 -0.711
## theta.1222 0.208 0.278 0.003 0.008 -0.342 0.022 0.204 0.398
## theta.1637 -0.029 0.270 0.003 0.008 -0.571 -0.209 -0.028 0.153
## theta.1995 -0.401 0.285 0.003 0.008 -0.971 -0.591 -0.400 -0.206
## theta.3369 1.719 0.480 0.006 0.011 0.845 1.379 1.696 2.028
## theta.4004 0.348 0.293 0.003 0.008 -0.224 0.145 0.345 0.543
## theta.4209 -0.211 0.279 0.003 0.008 -0.767 -0.390 -0.206 -0.022
## theta.4268 0.204 0.280 0.003 0.008 -0.347 0.015 0.205 0.391
## theta.4841 -0.485 0.288 0.003 0.008 -1.066 -0.668 -0.479 -0.287
## theta.5074 -0.064 0.273 0.003 0.008 -0.607 -0.244 -0.064 0.120
## theta.5495 -0.721 0.321 0.004 0.008 -1.387 -0.931 -0.713 -0.499
## theta.5710 -0.958 0.342 0.004 0.009 -1.664 -1.178 -0.945 -0.723
## theta.8294 0.290 0.279 0.003 0.008 -0.266 0.101 0.290 0.475
## theta.8980 0.803 0.333 0.004 0.009 0.167 0.572 0.793 1.021
## theta.9899 0.405 0.296 0.003 0.008 -0.177 0.203 0.402 0.609
## alpha.25 1.106 0.393 0.005 0.006 0.492 0.830 1.053 1.327
## alpha.26 2.027 0.626 0.007 0.011 1.033 1.576 1.959 2.384
## alpha.27 1.946 0.582 0.007 0.010 1.010 1.522 1.880 2.292
## alpha.28 1.114 0.416 0.005 0.006 0.452 0.816 1.058 1.349
## alpha.29 0.860 0.333 0.004 0.005 0.323 0.625 0.827 1.055
## alpha.30 0.691 0.284 0.003 0.004 0.229 0.494 0.655 0.854
## alpha.31 0.356 0.206 0.002 0.003 -0.008 0.215 0.337 0.482
## alpha.34 0.993 0.344 0.004 0.005 0.440 0.747 0.949 1.193
## alpha.35 0.515 0.235 0.003 0.003 0.126 0.351 0.490 0.650
## alpha.36 1.097 0.405 0.005 0.006 0.461 0.809 1.038 1.326
## alpha.37 1.018 0.389 0.004 0.006 0.420 0.745 0.962 1.231
## alpha.38 0.656 0.265 0.003 0.004 0.230 0.471 0.626 0.800
## 97.5% Model
## theta.1121 0.356 1d
## theta.1154 -0.310 1d
## theta.1222 0.763 1d
## theta.1637 0.496 1d
## theta.1995 0.144 1d
## theta.3369 2.742 1d
## theta.4004 0.935 1d
## theta.4209 0.328 1d
## theta.4268 0.750 1d
## theta.4841 0.066 1d
## theta.5074 0.465 1d
## theta.5495 -0.103 1d
## theta.5710 -0.317 1d
## theta.8294 0.843 1d
## theta.8980 1.484 1d
## theta.9899 0.985 1d
## alpha.25 2.028 1d
## alpha.26 3.446 1d
## alpha.27 3.289 1d
## alpha.28 2.065 1d
## alpha.29 1.609 1d
## alpha.30 1.348 1d
## alpha.31 0.794 1d
## alpha.34 1.781 1d
## alpha.35 1.047 1d
## alpha.36 2.025 1d
## alpha.37 1.947 1d
## alpha.38 1.288 1dPlot 1d Model posterior parameters from samples
Before plotting, set separate parameters
# Create a list of parameters for each type of parameter ahead
thetas <- grep("^theta", colnames(mcmc_list1d[[1]]), value = TRUE)
alphas <- grep("^alpha", colnames(mcmc_list1d[[1]]), value = TRUE)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_list1d, pars = thetas, facet_args = list(ncol = 4, scales = "fixed"))
Next the alphas
# Trace plot only for theta parameters
mcmc_trace(mcmc_list1d, pars = alphas, facet_args = list(ncol = 4, scales = "fixed"))
Next density plots of all posterior parameter distributions #### Density Plots
densplot(mcmc_list1d) 
Autocorrelation Plots
Next the autocorrelation plots
autocorr.plot(mcmc_list1d[[1]], auto.layout = FALSE, lag.max=20)
Create a tidy data table
summary1d <- summary1d %>%
mutate(
Variable = str_extract(rownames(.), "^[^\\.]+"),
ID = str_extract(rownames(.), "(?<=\\.)\\d+")
) %>%
`rownames<-`(NULL) %>%
dplyr::select(Model, Variable, ID, Mean, SD, 'Time-series SE', `50%`, `2.5%`, `97.5%` )Plot the posterior means of items
summary1d %>%
filter(Variable == "theta") %>%
arrange(`50%`) %>%
mutate(ID = factor(ID, levels = ID)) %>%
ggplot(aes(x = ID, y = `50%`)) +
geom_point() +
geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`), width = 0.2) +
labs(
title = "1-D Posterior Median Probit and 95% Credible Intervals for Theta Parameters",
x = "Item ID",
y = "Median (with 95% CI)"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
summary1d %>%
filter(Variable == "alpha") %>%
arrange(`50%`) %>%
mutate(ID = factor(ID, levels = ID)) %>%
ggplot(aes(x = ID, y = `50%`)) +
geom_point() +
geom_errorbar(aes(ymin = `2.5%`, ymax = `97.5%`), width = 0.2) +
labs(
title = "Posterior Median Probits and 95% Credible Intervals for Alpha Parameters",
x = "Rater ID",
y = "Median alpha (with 95% Credible Interval)"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
\(\alpha\) inversely related to the latent variance:
Large \(\alpha\) ⇒ smaller variance ⇒ more decisive preferences
Small \(\alpha\) ⇒ larger variance ⇒ more randomness in judgments