Last updated: 2022-01-20

Checks: 4 2

Knit directory: Padgett-Dissertation/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity itโ€™s best to always run the code in an empty environment.

The command set.seed(20210401) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

The following chunks had caches available:
  • model1
  • study4-model1-ppd

To ensure reproducibility of the results, delete the cache directory study4_model1_results_cache and re-run the analysis. To have workflowr automatically delete the cache directory prior to building the file, set delete_cache = TRUE when running wflow_build() or wflow_publish().

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Tracking code development and connecting the code version to the results is critical for reproducibility. To start using Git, open the Terminal and type git init in your project directory.


This project is not being versioned with Git. To obtain the full reproducibility benefits of using workflowr, please see ?wflow_start.


# Load packages & utility functions
source("code/load_packages.R")
source("code/load_utility_functions.R")
# environment options
options(scipen = 999, digits=3)

Describing the Observed Data

# Load diffIRT package with data
library(diffIRT)
data("extraversion")
mydata <- na.omit(extraversion)

# separate data then recombine
d1 <- mydata %>%
  as.data.frame() %>%
  select(contains("X"))%>%
  mutate(id = 1:n()) %>%
  pivot_longer(
    cols=contains("X"),
    names_to = c("item"),
    values_to = "Response"
  ) %>%
  mutate(
    item = ifelse(nchar(item)==4,substr(item, 3,3),substr(item, 3,4))
  )
d2 <- mydata %>%
  as.data.frame() %>%
  select(contains("T"))%>%
  mutate(id = 1:n()) %>%
  pivot_longer(
    cols=starts_with("T"),
    names_to = c("item"),
    values_to = "Time"
  ) %>%
  mutate(
    item = ifelse(nchar(item)==4,substr(item, 3,3),substr(item, 3,4))
  )
dat <- left_join(d1, d2)
Joining, by = c("id", "item")
dat_sum <- dat %>%
  select(item, Response, Time) %>%
  group_by(item) %>%
  summarize(
    M1 = mean(Response, na.rm=T),
    Mt = mean(Time, na.rm=T),
    SDt = sd(Time, na.rm=T),
    Mlogt = mean(log(Time), na.rm=T),
  )

colnames(dat_sum) <-
  c(
    "Item",
    "Proportion Endorsed",
    "Mean Response Time",
    "SD Response Time",
    "Mean Log Response Time"
  )

kable(dat_sum, format = "html", digits = 3) %>%
  kable_styling(full_width = T)
Item Proportion Endorsed Mean Response Time SD Response Time Mean Log Response Time
1 0.739 1.488 0.805 0.288
10 0.866 0.979 0.520 -0.115
2 0.535 1.354 0.648 0.208
3 0.852 1.115 0.632 0.002
4 0.923 1.001 0.664 -0.114
5 0.542 1.301 0.706 0.163
6 0.901 1.255 0.682 0.119
7 0.944 1.143 0.546 0.054
8 0.965 1.067 0.575 -0.030
9 0.824 1.728 0.745 0.463
# covariance among items
kable(cov(mydata[,colnames(mydata) %like% "X"]), digits = 3) %>%
  kable_styling(full_width = T)
X[1] X[2] X[3] X[4] X[5] X[6] X[7] X[8] X[9] X[10]
X[1] 0.194 -0.001 0.039 0.029 0.000 0.002 0.014 0.005 0.011 0.015
X[2] -0.001 0.251 0.023 0.006 0.077 0.011 0.002 0.012 0.031 0.030
X[3] 0.039 0.023 0.127 0.038 0.024 0.028 0.020 0.016 0.016 0.051
X[4] 0.029 0.006 0.038 0.072 0.014 0.006 0.017 0.019 0.029 0.025
X[5] 0.000 0.077 0.024 0.014 0.250 0.004 0.017 0.005 0.032 0.031
X[6] 0.002 0.011 0.028 0.006 0.004 0.090 0.009 0.011 0.004 0.015
X[7] 0.014 0.002 0.020 0.017 0.017 0.009 0.054 0.019 0.004 0.007
X[8] 0.005 0.012 0.016 0.019 0.005 0.011 0.019 0.034 0.008 0.009
X[9] 0.011 0.031 0.016 0.029 0.032 0.004 0.004 0.008 0.146 0.033
X[10] 0.015 0.030 0.051 0.025 0.031 0.015 0.007 0.009 0.033 0.117
# correlation matrix
psych::polychoric(mydata[,colnames(mydata) %like% "X"])
Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was done
Call: psych::polychoric(x = mydata[, colnames(mydata) %like% "X"])
Polychoric correlations 
      X[1]  X[2]  X[3]  X[4]  X[5]  X[6]  X[7]  X[8]  X[9]  X[10]
X[1]   1.00                                                      
X[2]  -0.01  1.00                                                
X[3]   0.45  0.24  1.00                                          
X[4]   0.50  0.11  0.70  1.00                                    
X[5]   0.00  0.46  0.26  0.23  1.00                              
X[6]   0.04  0.15  0.50  0.21  0.06  1.00                        
X[7]   0.32  0.05  0.52  0.58  0.36  0.32  1.00                  
X[8]   0.18  0.38  0.57  0.71  0.17  0.48  0.78  1.00            
X[9]   0.12  0.29  0.24  0.55  0.31  0.08  0.13  0.31  1.00      
X[10]  0.19  0.34  0.69  0.54  0.35  0.32  0.22  0.39  0.47  1.00

 with tau of 
           1
X[1]  -0.642
X[2]  -0.088
X[3]  -1.046
X[4]  -1.422
X[5]  -0.106
X[6]  -1.290
X[7]  -1.586
X[8]  -1.809
X[9]  -0.930
X[10] -1.109

Model 1: Traditional IFA

Model details

cat(read_file(paste0(w.d, "/code/study_4/model_1.txt")))
model {
### Model
  for(p in 1:N){
    for(i in 1:nit){
      # data model
      y[p,i] ~ dbern(pi[p,i,2])

      # LRV
      ystar[p,i] ~ dnorm(lambda[i]*eta[p], 1)

      # Pr(nu = 2)
      pi[p,i,2] = phi(ystar[p,i] - tau[i,1])
      # Pr(nu = 1)
      pi[p,i,1] = 1 - phi(ystar[p,i] - tau[i,1])

    }
  }
  ### Priors
  # person parameters
  for(p in 1:N){
    eta[p] ~ dnorm(0, 1) # latent ability
  }

  for(i in 1:nit){
    # Thresholds
    tau[i, 1] ~ dnorm(0.0,0.1)
    # loadings
    lambda[i] ~ dnorm(0, .44)T(0,)
    # LRV total variance
    # total variance = residual variance + fact. Var.
    theta[i] = 1 + pow(lambda[i],2)
    # standardized loading
    lambda.std[i] = lambda[i]/pow(theta[i],0.5)
  }

  # compute omega
  lambda_sum[1] = lambda[1]
  for(i in 2:nit){
    #lambda_sum (sum factor loadings)
    lambda_sum[i] = lambda_sum[i-1]+lambda[i]
  }
  reli.omega = (pow(lambda_sum[nit],2))/(pow(lambda_sum[nit],2)+nit)
}

Model results

# Save parameters
jags.params <- c("tau", "lambda", "theta", "reli.omega", "lambda.std")
# initial-values
jags.inits <- function(){
    list(
      "tau"=matrix(c(-0.64, -0.09, -1.05, -1.42, -0.11, -1.29, -1.59, -1.81, -0.93, -1.11),
                   ncol=1, nrow=10),
      "lambda"=rep(0.7,10),
      "eta"=rnorm(142),
      "ystar"=matrix(c(0.7*rep(rnorm(142),10)), ncol=10)
    )
  }
# data
jags.data <- list(y = mydata[,1:10],
               N = nrow(mydata),
               nit = 10)
model.fit <-  R2jags::jags(
  model = paste0(w.d, "/code/study_4/model_1.txt"),
  parameters.to.save = jags.params,
  inits = jags.inits,
  data = jags.data,
  n.chains = 4,
  n.burnin = 5000,
  n.iter = 10000
)
module glm loaded
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 1420
   Unobserved stochastic nodes: 1582
   Total graph size: 8742

Initializing model
print(model.fit, width=1000)
Inference for Bugs model at "C:/Users/noahp/Documents/GitHub/Padgett-Dissertation/code/study_4/model_1.txt", fit using jags,
 4 chains, each with 10000 iterations (first 5000 discarded), n.thin = 5
 n.sims = 4000 iterations saved
               mu.vect sd.vect    2.5%     25%     50%     75%   97.5% Rhat n.eff
lambda[1]        0.508   0.233   0.099   0.341   0.496   0.658   1.003 1.01   430
lambda[2]        0.657   0.251   0.215   0.480   0.643   0.813   1.184 1.00   650
lambda[3]        1.996   0.660   0.919   1.524   1.929   2.372   3.516 1.03   130
lambda[4]        1.657   0.595   0.690   1.225   1.589   2.033   3.010 1.02   160
lambda[5]        0.647   0.254   0.180   0.471   0.635   0.806   1.181 1.01   480
lambda[6]        0.828   0.346   0.228   0.586   0.804   1.034   1.598 1.01   350
lambda[7]        1.047   0.432   0.268   0.753   1.015   1.306   2.027 1.01   450
lambda[8]        1.533   0.612   0.564   1.104   1.449   1.864   3.019 1.00   760
lambda[9]        0.749   0.305   0.234   0.540   0.722   0.925   1.452 1.00   580
lambda[10]       1.567   0.513   0.742   1.205   1.504   1.848   2.779 1.01   310
lambda.std[1]    0.431   0.159   0.099   0.323   0.444   0.550   0.708 1.01   430
lambda.std[2]    0.525   0.145   0.210   0.432   0.541   0.631   0.764 1.00   700
lambda.std[3]    0.869   0.075   0.677   0.836   0.888   0.921   0.962 1.02   330
lambda.std[4]    0.823   0.099   0.568   0.775   0.846   0.897   0.949 1.01   210
lambda.std[5]    0.519   0.149   0.177   0.426   0.536   0.628   0.763 1.01   590
lambda.std[6]    0.601   0.160   0.223   0.506   0.627   0.719   0.848 1.01   400
lambda.std[7]    0.679   0.159   0.259   0.602   0.712   0.794   0.897 1.01   630
lambda.std[8]    0.797   0.119   0.491   0.741   0.823   0.881   0.949 1.02   540
lambda.std[9]    0.569   0.153   0.228   0.475   0.586   0.679   0.824 1.00   710
lambda.std[10]   0.816   0.089   0.596   0.769   0.833   0.879   0.941 1.01   350
reli.omega       0.922   0.022   0.875   0.911   0.926   0.937   0.953 1.01   370
tau[1,1]        -0.955   0.180  -1.331  -1.073  -0.950  -0.830  -0.622 1.00  1400
tau[2,1]        -0.132   0.167  -0.462  -0.242  -0.131  -0.019   0.198 1.00   910
tau[3,1]        -2.499   0.580  -3.905  -2.823  -2.425  -2.082  -1.613 1.03   140
tau[4,1]        -3.059   0.645  -4.570  -3.397  -2.975  -2.606  -2.052 1.02   200
tau[5,1]        -0.151   0.166  -0.481  -0.261  -0.150  -0.039   0.169 1.00  4000
tau[6,1]        -2.126   0.304  -2.803  -2.311  -2.100  -1.912  -1.616 1.01   510
tau[7,1]        -2.802   0.456  -3.846  -3.063  -2.748  -2.479  -2.053 1.00  1600
tau[8,1]        -3.767   0.812  -5.775  -4.192  -3.617  -3.195  -2.551 1.00  2200
tau[9,1]        -1.485   0.230  -1.971  -1.621  -1.474  -1.328  -1.080 1.01   470
tau[10,1]       -2.324   0.460  -3.404  -2.584  -2.260  -2.000  -1.601 1.01   610
theta[1]         1.313   0.275   1.010   1.116   1.246   1.433   2.006 1.00   620
theta[2]         1.494   0.367   1.046   1.230   1.413   1.662   2.401 1.00   650
theta[3]         5.421   2.964   1.845   3.322   4.721   6.626  13.360 1.03   100
theta[4]         4.098   2.268   1.476   2.501   3.524   5.133  10.063 1.02   160
theta[5]         1.483   0.367   1.032   1.222   1.403   1.650   2.395 1.01   340
theta[6]         1.806   0.679   1.052   1.344   1.646   2.069   3.552 1.01   360
theta[7]         2.283   1.040   1.072   1.567   2.030   2.705   5.110 1.01   340
theta[8]         3.724   2.237   1.318   2.220   3.100   4.475  10.116 1.00  1200
theta[9]         1.655   0.540   1.055   1.291   1.522   1.855   3.108 1.01   450
theta[10]        3.720   1.872   1.550   2.451   3.263   4.414   8.724 1.01   300
deviance       704.721  30.592 647.942 684.707 703.635 723.982 763.976 1.00  1600

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 467.4 and DIC = 1172.1
DIC is an estimate of expected predictive error (lower deviance is better).

Posterior Distribution Summary

# extract for plotting
jags.mcmc <- as.mcmc(model.fit)
a <- colnames(as.data.frame(jags.mcmc[[1]]))
fit.mcmc <- data.frame(as.matrix(jags.mcmc, chains = T, iters = T))
colnames(fit.mcmc) <- c("chain", "iter", a)
fit.mcmc.ggs <- ggmcmc::ggs(jags.mcmc) # for GRB plot

Categroy Thresholds (\(\tau\))

# tau
bayesplot::mcmc_areas(fit.mcmc, regex_pars = "tau", prob = 0.8); ggsave("fig/study4_model1_tau_dens.pdf")

Saving 7 x 5 in image
bayesplot::mcmc_acf(fit.mcmc, regex_pars = "tau"); ggsave("fig/study4_model1_tau_acf.pdf")

Saving 7 x 5 in image
bayesplot::mcmc_trace(fit.mcmc, regex_pars = "tau"); ggsave("fig/study4_model1_tau_trace.pdf")

Saving 7 x 5 in image
ggmcmc::ggs_grb(fit.mcmc.ggs, family = "tau"); ggsave("fig/study4_model1_tau_grb.pdf")

Saving 7 x 5 in image

Factor Loadings (\(\lambda\))

bayesplot::mcmc_areas(fit.mcmc, regex_pars = "lambda", prob = 0.8)

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "lambda")

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "lambda")

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "lambda")

bayesplot::mcmc_areas(fit.mcmc, regex_pars = "lambda.std", prob = 0.8); ggsave("fig/study4_model1_lambda_dens.pdf")

Saving 7 x 5 in image
bayesplot::mcmc_acf(fit.mcmc, regex_pars = "lambda.std"); ggsave("fig/study4_model1_lambda_acf.pdf")

Saving 7 x 5 in image
bayesplot::mcmc_trace(fit.mcmc, regex_pars = "lambda.std"); ggsave("fig/study4_model1_lambda_trace.pdf")

Saving 7 x 5 in image
ggmcmc::ggs_grb(fit.mcmc.ggs, family = "lambda.std"); ggsave("fig/study4_model1_lambda_grb.pdf")

Saving 7 x 5 in image

Latent Response Total Variance (\(\theta\))

bayesplot::mcmc_areas(fit.mcmc, regex_pars = "theta", prob = 0.8); ggsave("fig/study4_model1_theta_dens.pdf")

Saving 7 x 5 in image
bayesplot::mcmc_acf(fit.mcmc, regex_pars = "theta"); ggsave("fig/study4_model1_theta_acf.pdf")

Saving 7 x 5 in image
bayesplot::mcmc_trace(fit.mcmc, regex_pars = "theta"); ggsave("fig/study4_model1_theta_trace.pdf")

Saving 7 x 5 in image
ggmcmc::ggs_grb(fit.mcmc.ggs, family = "theta"); ggsave("fig/study4_model1_theta_grb.pdf")

Saving 7 x 5 in image

Factor Reliability Omega (\(\omega\))

bayesplot::mcmc_areas(fit.mcmc, regex_pars = "reli.omega", prob = 0.8); ggsave("fig/study4_model1_omega_dens.pdf")

Saving 7 x 5 in image
bayesplot::mcmc_acf(fit.mcmc, regex_pars = "reli.omega"); ggsave("fig/study4_model1_omega_acf.pdf")

Saving 7 x 5 in image
bayesplot::mcmc_trace(fit.mcmc, regex_pars = "reli.omega"); ggsave("fig/study4_model1_omega_trace.pdf")

Saving 7 x 5 in image
ggmcmc::ggs_grb(fit.mcmc.ggs, family = "reli.omega"); ggsave("fig/study4_model1_omega_grb.pdf")

Saving 7 x 5 in image
# extract omega posterior for results comparison
extracted_omega <- data.frame(model_1 = fit.mcmc$reli.omega)
write.csv(x=extracted_omega, file=paste0(getwd(),"/data/study_4/extracted_omega_m1.csv"))

Posterior Predictive Distributions

# Posterior Predictive Check
Niter <- 200
model.fit$model$recompile()
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 1420
   Unobserved stochastic nodes: 1582
   Total graph size: 8742

Initializing model
fit.extra <- rjags::jags.samples(model.fit$model, variable.names = "pi", n.iter = Niter)
N <- model.fit$model$data()[["N"]]
nit <- 10
nchain=4
C <- 2
n <- i <- iter <- ppc.row.i <- 1
y.prob.ppc <- array(dim=c(Niter*nchain, nit, C))
for(chain in 1:nchain){
  for(iter in 1:Niter){
    # initialize simulated y for this iteration
    y <- matrix(nrow=N, ncol=nit)
    # loop over item
    for(i in 1:nit){
      # simulated data for item i for each person
      for(n in 1:N){
        y[n,i] <- sample(1:C, 1, prob = fit.extra$pi[n, i, 1:C, iter, chain])
      }
      # computer proportion of each response category
      for(c in 1:C){
        y.prob.ppc[ppc.row.i,i,c] <- sum(y[,i]==c)/N
      }
    }
    
    # update row of output
    ppc.row.i = ppc.row.i + 1
  }
}

yppcmat <- matrix(c(y.prob.ppc), ncol=1)
z <- expand.grid(1:(Niter*nchain), 1:nit, 1:C)
yppcmat <- data.frame(  iter = z[,1], nit=z[,2], C=z[,3], yppc = yppcmat)

ymat <- model.fit$model$data()[["y"]]
y.prob <- matrix(ncol=C, nrow=nit)
for(i in 1:nit){
  for(c in 1:C){
    y.prob[i,c] <- sum(ymat[,i]==c-1)/N
  }
}
yobsmat <- matrix(c(y.prob), ncol=1)
z <- expand.grid(1:nit, 1:C)
yobsmat <- data.frame(nit=z[,1], C=z[,2], yobs = yobsmat)
plot.ppc <- full_join(yppcmat, yobsmat)
Joining, by = c("nit", "C")
p <- plot.ppc %>%
  mutate(C    = as.factor(C),
         item = nit) %>%
  ggplot()+
  geom_boxplot(aes(x=C,y=y.prob.ppc), outlier.colour = NA)+
  geom_point(aes(x=C,y=yobs), color="red")+
  lims(y=c(0, 1))+
  labs(y="Posterior Predictive Category Proportion", x="Item Category")+
  facet_wrap(.~nit, nrow=1)+
  theme_bw()+
  theme(
    panel.grid = element_blank(),
    strip.background = element_rect(fill="white")
  )
p

ggsave(filename = "fig/study4_model1_ppc_y.pdf",plot=p,width = 6, height=3,units="in")
ggsave(filename = "fig/study4_model1_ppc_y.png",plot=p,width = 6, height=3,units="in")
ggsave(filename = "fig/study4_model1_ppc_y.eps",plot=p,width = 6, height=3,units="in")

Manuscript Table and Figures

Table

# print to xtable
print(
  xtable(
    model.fit$BUGSoutput$summary,
    caption = c("study4 Model 1 posterior distribution summary")
    ,align = "lrrrrrrrrr"
  ),
  include.rownames=T,
  booktabs=T
)
% latex table generated in R 4.0.5 by xtable 1.8-4 package
% Thu Jan 20 13:28:36 2022
\begin{table}[ht]
\centering
\begin{tabular}{lrrrrrrrrr}
  \toprule
 & mean & sd & 2.5\% & 25\% & 50\% & 75\% & 97.5\% & Rhat & n.eff \\ 
  \midrule
deviance & 704.72 & 30.59 & 647.94 & 684.71 & 703.64 & 723.98 & 763.98 & 1.00 & 1600.00 \\ 
  lambda[1] & 0.51 & 0.23 & 0.10 & 0.34 & 0.50 & 0.66 & 1.00 & 1.01 & 430.00 \\ 
  lambda[2] & 0.66 & 0.25 & 0.21 & 0.48 & 0.64 & 0.81 & 1.18 & 1.00 & 650.00 \\ 
  lambda[3] & 2.00 & 0.66 & 0.92 & 1.52 & 1.93 & 2.37 & 3.52 & 1.03 & 130.00 \\ 
  lambda[4] & 1.66 & 0.60 & 0.69 & 1.23 & 1.59 & 2.03 & 3.01 & 1.02 & 160.00 \\ 
  lambda[5] & 0.65 & 0.25 & 0.18 & 0.47 & 0.63 & 0.81 & 1.18 & 1.01 & 480.00 \\ 
  lambda[6] & 0.83 & 0.35 & 0.23 & 0.59 & 0.80 & 1.03 & 1.60 & 1.01 & 350.00 \\ 
  lambda[7] & 1.05 & 0.43 & 0.27 & 0.75 & 1.01 & 1.31 & 2.03 & 1.01 & 450.00 \\ 
  lambda[8] & 1.53 & 0.61 & 0.56 & 1.10 & 1.45 & 1.86 & 3.02 & 1.00 & 760.00 \\ 
  lambda[9] & 0.75 & 0.31 & 0.23 & 0.54 & 0.72 & 0.92 & 1.45 & 1.00 & 580.00 \\ 
  lambda[10] & 1.57 & 0.51 & 0.74 & 1.20 & 1.50 & 1.85 & 2.78 & 1.01 & 310.00 \\ 
  lambda.std[1] & 0.43 & 0.16 & 0.10 & 0.32 & 0.44 & 0.55 & 0.71 & 1.01 & 430.00 \\ 
  lambda.std[2] & 0.53 & 0.15 & 0.21 & 0.43 & 0.54 & 0.63 & 0.76 & 1.00 & 700.00 \\ 
  lambda.std[3] & 0.87 & 0.08 & 0.68 & 0.84 & 0.89 & 0.92 & 0.96 & 1.02 & 330.00 \\ 
  lambda.std[4] & 0.82 & 0.10 & 0.57 & 0.77 & 0.85 & 0.90 & 0.95 & 1.01 & 210.00 \\ 
  lambda.std[5] & 0.52 & 0.15 & 0.18 & 0.43 & 0.54 & 0.63 & 0.76 & 1.01 & 590.00 \\ 
  lambda.std[6] & 0.60 & 0.16 & 0.22 & 0.51 & 0.63 & 0.72 & 0.85 & 1.01 & 400.00 \\ 
  lambda.std[7] & 0.68 & 0.16 & 0.26 & 0.60 & 0.71 & 0.79 & 0.90 & 1.01 & 630.00 \\ 
  lambda.std[8] & 0.80 & 0.12 & 0.49 & 0.74 & 0.82 & 0.88 & 0.95 & 1.02 & 540.00 \\ 
  lambda.std[9] & 0.57 & 0.15 & 0.23 & 0.47 & 0.59 & 0.68 & 0.82 & 1.00 & 710.00 \\ 
  lambda.std[10] & 0.82 & 0.09 & 0.60 & 0.77 & 0.83 & 0.88 & 0.94 & 1.01 & 350.00 \\ 
  reli.omega & 0.92 & 0.02 & 0.88 & 0.91 & 0.93 & 0.94 & 0.95 & 1.01 & 370.00 \\ 
  tau[1,1] & -0.96 & 0.18 & -1.33 & -1.07 & -0.95 & -0.83 & -0.62 & 1.00 & 1400.00 \\ 
  tau[2,1] & -0.13 & 0.17 & -0.46 & -0.24 & -0.13 & -0.02 & 0.20 & 1.00 & 910.00 \\ 
  tau[3,1] & -2.50 & 0.58 & -3.90 & -2.82 & -2.43 & -2.08 & -1.61 & 1.03 & 140.00 \\ 
  tau[4,1] & -3.06 & 0.64 & -4.57 & -3.40 & -2.98 & -2.61 & -2.05 & 1.02 & 200.00 \\ 
  tau[5,1] & -0.15 & 0.17 & -0.48 & -0.26 & -0.15 & -0.04 & 0.17 & 1.00 & 4000.00 \\ 
  tau[6,1] & -2.13 & 0.30 & -2.80 & -2.31 & -2.10 & -1.91 & -1.62 & 1.01 & 510.00 \\ 
  tau[7,1] & -2.80 & 0.46 & -3.85 & -3.06 & -2.75 & -2.48 & -2.05 & 1.00 & 1600.00 \\ 
  tau[8,1] & -3.77 & 0.81 & -5.77 & -4.19 & -3.62 & -3.19 & -2.55 & 1.00 & 2200.00 \\ 
  tau[9,1] & -1.48 & 0.23 & -1.97 & -1.62 & -1.47 & -1.33 & -1.08 & 1.01 & 470.00 \\ 
  tau[10,1] & -2.32 & 0.46 & -3.40 & -2.58 & -2.26 & -2.00 & -1.60 & 1.01 & 610.00 \\ 
  theta[1] & 1.31 & 0.28 & 1.01 & 1.12 & 1.25 & 1.43 & 2.01 & 1.00 & 620.00 \\ 
  theta[2] & 1.49 & 0.37 & 1.05 & 1.23 & 1.41 & 1.66 & 2.40 & 1.00 & 650.00 \\ 
  theta[3] & 5.42 & 2.96 & 1.85 & 3.32 & 4.72 & 6.63 & 13.36 & 1.03 & 100.00 \\ 
  theta[4] & 4.10 & 2.27 & 1.48 & 2.50 & 3.52 & 5.13 & 10.06 & 1.02 & 160.00 \\ 
  theta[5] & 1.48 & 0.37 & 1.03 & 1.22 & 1.40 & 1.65 & 2.39 & 1.01 & 340.00 \\ 
  theta[6] & 1.81 & 0.68 & 1.05 & 1.34 & 1.65 & 2.07 & 3.55 & 1.01 & 360.00 \\ 
  theta[7] & 2.28 & 1.04 & 1.07 & 1.57 & 2.03 & 2.70 & 5.11 & 1.01 & 340.00 \\ 
  theta[8] & 3.72 & 2.24 & 1.32 & 2.22 & 3.10 & 4.48 & 10.12 & 1.00 & 1200.00 \\ 
  theta[9] & 1.65 & 0.54 & 1.05 & 1.29 & 1.52 & 1.85 & 3.11 & 1.01 & 450.00 \\ 
  theta[10] & 3.72 & 1.87 & 1.55 & 2.45 & 3.26 & 4.41 & 8.72 & 1.01 & 300.00 \\ 
   \bottomrule
\end{tabular}
\caption{study4 Model 1 posterior distribution summary} 
\end{table}

Figure


sessionInfo()
R version 4.0.5 (2021-03-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22000)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] car_3.0-10           carData_3.0-4        mvtnorm_1.1-1       
 [4] LaplacesDemon_16.1.4 runjags_2.2.0-2      lme4_1.1-26         
 [7] Matrix_1.3-2         sirt_3.9-4           R2jags_0.6-1        
[10] rjags_4-12           eRm_1.0-2            diffIRT_1.5         
[13] statmod_1.4.35       xtable_1.8-4         kableExtra_1.3.4    
[16] lavaan_0.6-7         polycor_0.7-10       bayesplot_1.8.0     
[19] ggmcmc_1.5.1.1       coda_0.19-4          data.table_1.14.0   
[22] patchwork_1.1.1      forcats_0.5.1        stringr_1.4.0       
[25] dplyr_1.0.5          purrr_0.3.4          readr_1.4.0         
[28] tidyr_1.1.3          tibble_3.1.0         ggplot2_3.3.5       
[31] tidyverse_1.3.0      workflowr_1.6.2     

loaded via a namespace (and not attached):
 [1] minqa_1.2.4        TAM_3.5-19         colorspace_2.0-0   rio_0.5.26        
 [5] ellipsis_0.3.1     ggridges_0.5.3     rprojroot_2.0.2    fs_1.5.0          
 [9] rstudioapi_0.13    farver_2.1.0       fansi_0.4.2        lubridate_1.7.10  
[13] xml2_1.3.2         codetools_0.2-18   splines_4.0.5      mnormt_2.0.2      
[17] knitr_1.31         jsonlite_1.7.2     nloptr_1.2.2.2     broom_0.7.5       
[21] dbplyr_2.1.0       compiler_4.0.5     httr_1.4.2         backports_1.2.1   
[25] assertthat_0.2.1   cli_2.3.1          later_1.1.0.1      htmltools_0.5.1.1 
[29] tools_4.0.5        gtable_0.3.0       glue_1.4.2         reshape2_1.4.4    
[33] Rcpp_1.0.7         cellranger_1.1.0   jquerylib_0.1.3    vctrs_0.3.6       
[37] svglite_2.0.0      nlme_3.1-152       psych_2.0.12       xfun_0.21         
[41] ps_1.6.0           openxlsx_4.2.3     rvest_1.0.0        lifecycle_1.0.0   
[45] MASS_7.3-53.1      scales_1.1.1       ragg_1.1.1         hms_1.0.0         
[49] promises_1.2.0.1   parallel_4.0.5     RColorBrewer_1.1-2 curl_4.3          
[53] yaml_2.2.1         sass_0.3.1         reshape_0.8.8      stringi_1.5.3     
[57] highr_0.8          zip_2.1.1          boot_1.3-27        rlang_0.4.10      
[61] pkgconfig_2.0.3    systemfonts_1.0.1  evaluate_0.14      lattice_0.20-41   
[65] labeling_0.4.2     tidyselect_1.1.0   GGally_2.1.1       plyr_1.8.6        
[69] magrittr_2.0.1     R6_2.5.0           generics_0.1.0     DBI_1.1.1         
[73] foreign_0.8-81     pillar_1.5.1       haven_2.3.1        withr_2.4.1       
[77] abind_1.4-5        modelr_0.1.8       crayon_1.4.1       utf8_1.1.4        
[81] tmvnsim_1.0-2      rmarkdown_2.7      grid_4.0.5         readxl_1.3.1      
[85] CDM_7.5-15         pbivnorm_0.6.0     git2r_0.28.0       reprex_1.0.0      
[89] digest_0.6.27      webshot_0.5.2      httpuv_1.5.5       textshaping_0.3.1 
[93] stats4_4.0.5       munsell_0.5.0      viridisLite_0.3.0  bslib_0.2.4       
[97] R2WinBUGS_2.1-21