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

View a summary of model developed here: https://rpubs.com/steveneheil/pwc-MCMC-CDAT

load("Out/summary1d.Rdata")
summary1d
##              Mean    SD Naive SE Time-series SE   2.5%    25%    50%    75%
## theta.1121 -0.117 0.277    0.012          0.022 -0.687 -0.287 -0.117  0.068
## theta.1154 -0.946 0.328    0.014          0.023 -1.637 -1.142 -0.939 -0.739
## theta.1222  0.128 0.279    0.012          0.023 -0.435 -0.062  0.137  0.321
## theta.1637  0.018 0.283    0.012          0.022 -0.577 -0.155  0.020  0.210
## theta.1995 -0.479 0.292    0.013          0.025 -1.063 -0.653 -0.481 -0.283
## theta.3369  1.722 0.460    0.020          0.030  0.958  1.383  1.701  2.022
## theta.4004  0.372 0.302    0.013          0.024 -0.280  0.189  0.379  0.569
## theta.4209 -0.232 0.288    0.012          0.024 -0.789 -0.430 -0.233 -0.034
## theta.4268  0.282 0.290    0.012          0.023 -0.300  0.098  0.296  0.485
## theta.4841 -0.556 0.293    0.013          0.024 -1.137 -0.732 -0.551 -0.367
## theta.5074 -0.145 0.279    0.012          0.024 -0.715 -0.326 -0.129  0.035
## theta.5495 -0.833 0.315    0.014          0.023 -1.465 -1.052 -0.828 -0.604
## theta.5710 -1.007 0.335    0.014          0.024 -1.650 -1.217 -1.004 -0.778
## theta.8294  0.159 0.288    0.012          0.024 -0.450 -0.024  0.165  0.355
## theta.8980  0.703 0.316    0.014          0.025  0.086  0.487  0.718  0.917
## theta.9899  0.521 0.303    0.013          0.022 -0.098  0.329  0.522  0.705
## alpha.1     2.696 0.731    0.031          0.037  1.444  2.207  2.663  3.131
## alpha.25    1.129 0.340    0.015          0.015  0.562  0.886  1.113  1.349
## alpha.26    1.881 0.554    0.024          0.028  1.039  1.481  1.807  2.175
## alpha.27    1.795 0.531    0.023          0.025  0.924  1.425  1.714  2.101
## alpha.28    0.985 0.339    0.015          0.017  0.431  0.755  0.929  1.211
## alpha.29    0.853 0.305    0.013          0.013  0.365  0.640  0.828  1.025
## alpha.30    0.703 0.260    0.011          0.012  0.249  0.516  0.683  0.848
## alpha.31    0.334 0.188    0.008          0.009  0.019  0.205  0.320  0.448
## alpha.34    1.048 0.313    0.013          0.012  0.519  0.815  1.022  1.249
## alpha.35    0.468 0.203    0.009          0.010  0.130  0.310  0.447  0.614
## alpha.36    1.004 0.354    0.015          0.015  0.428  0.737  0.977  1.209
## alpha.37    1.028 0.357    0.015          0.015  0.466  0.781  0.984  1.227
## alpha.38    0.618 0.222    0.010          0.011  0.257  0.465  0.594  0.756
##             97.5% Model
## theta.1121  0.422    1d
## theta.1154 -0.272    1d
## theta.1222  0.664    1d
## theta.1637  0.553    1d
## theta.1995  0.093    1d
## theta.3369  2.630    1d
## theta.4004  0.952    1d
## theta.4209  0.319    1d
## theta.4268  0.840    1d
## theta.4841  0.006    1d
## theta.5074  0.395    1d
## theta.5495 -0.218    1d
## theta.5710 -0.395    1d
## theta.8294  0.720    1d
## theta.8980  1.293    1d
## theta.9899  1.140    1d
## alpha.1     4.401    1d
## alpha.25    1.851    1d
## alpha.26    3.295    1d
## alpha.27    2.968    1d
## alpha.28    1.713    1d
## alpha.29    1.553    1d
## alpha.30    1.290    1d
## alpha.31    0.730    1d
## alpha.34    1.687    1d
## alpha.35    0.932    1d
## alpha.36    1.863    1d
## alpha.37    1.812    1d
## alpha.38    1.162    1d

Plot 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) %>%
  select(Model, Variable, ID, everything() )

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