library(bayesplot)
library(coda)
library(ggplot2)
library(tidyverse)
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)
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
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)
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)
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() )
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