Last updated: 2022-01-16

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.

Environment: empty

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.

Seed: set.seed(20210401)

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.

Session information: recorded

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

Cache: detected

The following chunks had caches available:

model2
study1-model2-ppd

To ensure reproducibility of the results, delete the cache directory study1_model2_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().

File paths: relative

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

Repository version: no version control

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)
# generate data for study 1
source("code/study_1/study_1_generate_data.R")

Simulated Data

# data parameters
paravec <- c(
  N = 500
  , J = 5 # N_items
  , C = 3 # N_cat
  , etaCor = .23
  , etasd1 = 1
  , etasd2 = sqrt(0.1)
  , lambda=0.7
  , nu=1.5
  , sigma.ei=0.25
  , rho1=0.1
)
# simulated then saved below
sim_tau <- matrix(
  c(-0.822, -0.751, -0.616, -0.392, -0.865,
    0.780, 0.882, 0.827, 1.030, 0.877),
  ncol=2, nrow=5
)
# Use parameters to simulate data
sim.data <- simulate_data_misclass(paravec, tau=sim_tau)

Describing the Observed (simulated) Data

d1 <- sim.data$Ysampled %>%
  as.data.frame() %>%
  select(contains("y")) %>%
  mutate(id = 1:n()) %>%
  pivot_longer(
    cols = contains("y"),
    names_to = c("item"),
    values_to = "Response"
  ) %>%
  mutate(item = ifelse(nchar(item) > 2, substr(item, 2, 3), substr(item, 2, 2)))
d2 <- sim.data$logt %>%
  as.data.frame() %>%
  select(contains("logt")) %>%
  mutate(id = 1:n()) %>%
  pivot_longer(
    cols = contains("logt"),
    names_to = c("item"),
    values_to = "Time"
  ) %>%
  mutate(item = ifelse(nchar(item) > 5, substr(item, 5, 6), substr(item, 5, 5)))
dat <- left_join(d1, d2)

Joining, by = c("id", "item")

dat_sum <- dat %>%
  select(item, Response, Time) %>%
  group_by(item) %>%
  summarize(
    p1 = table(Response)[1] / n(),
    p2 = table(Response)[2] / n(),
    p3 = table(Response)[3] / n(),
    M1 = mean(Response, na.rm = T),
    Mt = mean(Time, na.rm = T),
    SDt = sd(Time, na.rm = T)
  )

colnames(dat_sum) <-
  c(
    "Item",
    "Prop. R == 1",
    "Prop. R == 2",
    "Prop. R == 3",
    "Mean Response",
    "Mean Response Time",
    "SD Response Time"
  )
dat_sum$Item <- paste0("item_", 1:N_items)

kable(dat_sum, format = "html", digits = 3) %>%
  kable_styling(full_width = T)

Item	Prop. R == 1	Prop. R == 2	Prop. R == 3	Mean Response	Mean Response Time	SD Response Time
item_1	0.308	0.404	0.288	1.98	1.39	0.597
item_2	0.310	0.414	0.276	1.97	1.43	0.618
item_3	0.338	0.386	0.276	1.94	1.43	0.613
item_4	0.362	0.384	0.254	1.89	1.40	0.592
item_5	0.292	0.422	0.286	1.99	1.36	0.582

# covariance among items
cov(sim.data$Ysampled)

       y1     y2     y3     y4     y5
y1 0.5968 0.0634 0.0428 0.0640 0.0319
y2 0.0634 0.5860 0.0440 0.0364 0.0258
y3 0.0428 0.0440 0.6114 0.0394 0.0457
y4 0.0640 0.0364 0.0394 0.6055 0.0655
y5 0.0319 0.0258 0.0457 0.0655 0.5791

# correlation matrix
psych::polychoric(sim.data$Ysampled)

Call: psych::polychoric(x = sim.data$Ysampled)
Polychoric correlations 
   y1   y2   y3   y4   y5  
y1 1.00                    
y2 0.14 1.00               
y3 0.09 0.09 1.00          
y4 0.13 0.08 0.08 1.00     
y5 0.07 0.05 0.09 0.14 1.00

 with tau of 
       1    2
y1 -0.50 0.56
y2 -0.50 0.59
y3 -0.42 0.59
y4 -0.35 0.66
y5 -0.55 0.57

Model 2: IFA with RT

Model details

cat(read_file(paste0(w.d, "/code/study_1/model_2.txt")))

model {
### Model
  for(p in 1:N){
    for(i in 1:nit){
      # data model
      y[p,i] ~ dcat(pi[p,i, ])

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

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

      # log-RT model
      dev[p,i]<-lambda[i]*(eta[p] - (tau[i,1]+tau[i,2])/2)
      lrt[p,i] ~ dnorm(icept[i] - speed[p] - rho * abs(dev[p,i]), prec[i])

    }
  }
  ### Priors
  # person parameters
  for(p in 1:N){
    eta[p] ~ dnorm(0, 1) # latent ability
    speed[p]~dnorm(sigma.ts*eta[p],prec.s)  # latent speed
  }
  sigma.ts ~ dnorm(0, 0.1)
  prec.s ~ dgamma(.1,.1)
  for(i in 1:nit){
    # lrt parameters
    icept[i]~dnorm(0,.1)
    prec[i]~dgamma(.1,.1)
    # Thresholds
    tau[i, 1] ~ dnorm(0.0,0.1)
    tau[i, 2] ~ dnorm(0, 0.1)T(tau[i, 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)
  }
  rho~dnorm(0,.1)I(0,)

  # 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

This model is similar to the BL-IRT model for jointly modeling item responses and response times (Molenaar et al., 2015).

# Save parameters
jags.params <- c("tau",
                 "lambda","lambda.std",
                 "theta",
                 "icept",
                 "prec",
                 "prec.s",
                 "sigma.ts",
                 "rho",
                 "reli.omega")
# initial-values
jags.inits <- function(){
    list(
      "tau"=matrix(c(-0.822, -0.751, -0.616, -0.392, -0.865,
                     0.780, 0.882, 0.827, 1.030, 0.877), ncol=2, nrow=5),
      "lambda"=rep(0.7,5),
      "rho"=0.1,
      "icept"=rep(1.5, 5),
      "prec.s"=10,
      "prec"=rep(4, 5),
      "sigma.ts"=0.1,
      "eta"=sim.data$eta[,1,drop=T],
      "speed"=sim.data$eta[,2,drop=T],
      "ystar"=t(sim.data$ystar)
    )
  }

mydata <- list(
  y = sim.data$Ysampled,
  lrt = sim.data$logt,
  N = nrow(sim.data$Ysampled),
  nit = ncol(sim.data$Ysampled)
)

# Run model
# Model 2

model.fit <-  R2jags::jags(
  model = paste0(w.d, "/code/study_1/model_2.txt"),
  parameters.to.save = jags.params,
  inits = jags.inits,
  data = mydata,
  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: 5000
   Unobserved stochastic nodes: 3528
   Total graph size: 44073

Initializing model

print(model.fit, width=1000)

Inference for Bugs model at "C:/Users/noahp/Documents/GitHub/Padgett-Dissertation/code/study_1/model_2.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
icept[1]         1.542   0.061    1.433    1.499    1.540    1.583    1.667 1.01   180
icept[2]         1.552   0.063    1.438    1.508    1.547    1.591    1.687 1.03    99
icept[3]         1.655   0.082    1.492    1.599    1.655    1.711    1.810 1.01   200
icept[4]         1.571   0.063    1.458    1.529    1.569    1.612    1.696 1.01   190
icept[5]         1.569   0.078    1.423    1.515    1.569    1.623    1.722 1.02   140
lambda[1]        0.438   0.122    0.219    0.357    0.431    0.510    0.714 1.01   550
lambda[2]        0.318   0.109    0.094    0.247    0.318    0.392    0.533 1.04   240
lambda[3]        0.607   0.158    0.318    0.499    0.599    0.705    0.947 1.01   200
lambda[4]        0.448   0.124    0.232    0.363    0.441    0.521    0.706 1.01   670
lambda[5]        0.565   0.134    0.311    0.475    0.563    0.650    0.833 1.01   530
lambda.std[1]    0.395   0.092    0.214    0.336    0.396    0.455    0.581 1.01   500
lambda.std[2]    0.299   0.094    0.094    0.240    0.303    0.365    0.470 1.04   230
lambda.std[3]    0.509   0.098    0.303    0.447    0.514    0.576    0.688 1.01   260
lambda.std[4]    0.402   0.091    0.226    0.342    0.403    0.462    0.577 1.00   800
lambda.std[5]    0.484   0.089    0.297    0.429    0.490    0.545    0.640 1.01   550
prec[1]          3.982   0.297    3.432    3.779    3.969    4.177    4.590 1.00  2500
prec[2]          3.956   0.291    3.423    3.753    3.944    4.142    4.570 1.00  4000
prec[3]          4.174   0.354    3.553    3.932    4.153    4.385    4.929 1.00   820
prec[4]          4.122   0.310    3.560    3.910    4.113    4.317    4.761 1.00  3800
prec[5]          4.707   0.391    4.014    4.433    4.694    4.947    5.539 1.00   960
prec.s          10.391   1.473    8.063    9.361   10.213   11.203   13.802 1.02   180
reli.omega       0.527   0.057    0.402    0.492    0.532    0.567    0.623 1.01   240
rho              0.465   0.144    0.202    0.369    0.456    0.554    0.766 1.04    75
sigma.ts         0.077   0.030    0.017    0.057    0.077    0.097    0.134 1.00  3300
tau[1,1]        -0.751   0.088   -0.928   -0.811   -0.751   -0.692   -0.582 1.00  4000
tau[2,1]        -0.731   0.084   -0.899   -0.789   -0.729   -0.676   -0.568 1.00  1400
tau[3,1]        -0.674   0.092   -0.859   -0.736   -0.672   -0.612   -0.498 1.01   420
tau[4,1]        -0.497   0.083   -0.657   -0.554   -0.497   -0.441   -0.330 1.00  3800
tau[5,1]        -0.826   0.088   -0.999   -0.885   -0.827   -0.765   -0.657 1.00  3400
tau[1,2]         0.824   0.089    0.652    0.764    0.824    0.884    0.998 1.00  3500
tau[2,2]         0.853   0.085    0.690    0.796    0.852    0.910    1.018 1.00  2900
tau[3,2]         0.884   0.092    0.703    0.822    0.882    0.945    1.064 1.00  1100
tau[4,2]         1.016   0.092    0.837    0.954    1.015    1.079    1.196 1.00  1900
tau[5,2]         0.874   0.093    0.694    0.812    0.873    0.935    1.062 1.00   750
theta[1]         1.207   0.116    1.048    1.127    1.186    1.260    1.509 1.00  1100
theta[2]         1.113   0.072    1.009    1.061    1.101    1.153    1.284 1.00   640
theta[3]         1.393   0.202    1.101    1.249    1.359    1.497    1.896 1.03   120
theta[4]         1.216   0.122    1.054    1.132    1.194    1.271    1.499 1.02   330
theta[5]         1.337   0.155    1.097    1.226    1.317    1.423    1.694 1.00   530
deviance      7467.165  76.519 7320.231 7415.326 7467.329 7518.652 7618.684 1.00  2600

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 = 2926.5 and DIC = 10393.6
DIC is an estimate of expected predictive error (lower deviance is better).

Posterior Distribution Summary

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

# save posterior draws for later
write.csv(x=fit.mcmc, file=paste0(getwd(),"/data/study_1/posterior_draws_m2.csv"))

Categroy Thresholds (\(\tau\))

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

Saving 7 x 5 in image

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "tau"); ggsave("fig/study1_model2_tau_acf.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "tau"); ggsave("fig/study1_model2_tau_trace.pdf")

Saving 7 x 5 in image

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "tau"); ggsave("fig/study1_model2_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/study1_model2_lambda_dens.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "lambda.std"); ggsave("fig/study1_model2_lambda_acf.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "lambda.std"); ggsave("fig/study1_model2_lambda_trace.pdf")

Saving 7 x 5 in image

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "lambda.std"); ggsave("fig/study1_model2_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/study1_model2_theta_dens.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "theta"); ggsave("fig/study1_model2_theta_acf.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "theta"); ggsave("fig/study1_model2_theta_trace.pdf")

Saving 7 x 5 in image

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "theta"); ggsave("fig/study1_model2_theta_grb.pdf")

Saving 7 x 5 in image

Response Time Intercept (\(\beta_{lrt}\))

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

Saving 7 x 5 in image

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "icept"); ggsave("fig/study1_model2_icept_acf.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "icept"); ggsave("fig/study1_model2_icept_trace.pdf")

Saving 7 x 5 in image

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "icept"); ggsave("fig/study1_model2_icept_grb.pdf")

Saving 7 x 5 in image

Response Time Precision (\(\sigma_{lrt}\))

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

Saving 7 x 5 in image

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "prec"); ggsave("fig/study1_model2_prec_acf.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "prec"); ggsave("fig/study1_model2_prec_trace.pdf")

Saving 7 x 5 in image

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "prec"); ggsave("fig/study1_model2_prec_grb.pdf")

Saving 7 x 5 in image

Speed Factor Variance (\(\sigma_s\))

bayesplot::mcmc_areas(fit.mcmc, regex_pars = "prec.s", prob = 0.8); ggsave("fig/study1_model2_precs_dens.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "prec.s"); ggsave("fig/study1_model2_precs_acf.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "prec.s"); ggsave("fig/study1_model2_precs_trace.pdf")

Saving 7 x 5 in image

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "prec.s"); ggsave("fig/study1_model2_precs_grb.pdf")

Saving 7 x 5 in image

Factor Covariance (\(\sigma_{ts}\))

bayesplot::mcmc_areas(fit.mcmc, regex_pars = "sigma.ts", prob = 0.8); ggsave("fig/study1_model2_sigmats_dens.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "sigma.ts"); ggsave("fig/study1_model2_sigmats_acf.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "sigma.ts"); ggsave("fig/study1_model2_sigmats_trace.pdf")

Saving 7 x 5 in image

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "sigma.ts"); ggsave("fig/study1_model2_sigmats_grb.pdf")

Saving 7 x 5 in image

PID (\(\rho\))

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

Saving 7 x 5 in image

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "rho"); ggsave("fig/study1_model2_rho_acf.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "rho"); ggsave("fig/study1_model2_rho_trace.pdf")

Saving 7 x 5 in image

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "rho"); ggsave("fig/study1_model2_rho_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/study1_model2_omega_dens.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_acf(fit.mcmc, regex_pars = "reli.omega"); ggsave("fig/study1_model2_omega_acf.pdf")

Saving 7 x 5 in image

bayesplot::mcmc_trace(fit.mcmc, regex_pars = "reli.omega"); ggsave("fig/study1_model2_omega_trace.pdf")

Saving 7 x 5 in image

ggmcmc::ggs_grb(fit.mcmc.ggs, family = "reli.omega"); ggsave("fig/study1_model2_omega_grb.pdf")

Saving 7 x 5 in image

# extract omega posterior for results comparison
extracted_omega <- data.frame(model_2 = fit.mcmc$reli.omega)
write.csv(x=extracted_omega, file=paste0(getwd(),"/data/study_1/extracted_omega_m2.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: 5000
   Unobserved stochastic nodes: 3528
   Total graph size: 44073

Initializing model

fit.extra <- rjags::jags.samples(model.fit$model, variable.names = "pi", n.iter = Niter)

NOTE: Stopping adaptation

N <- model.fit$model$data()[[1]]
nit <- 5
nchain=4
C <- 3
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()[[4]]
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)/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, 0.67))+
  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/study1_model2_ppc_y.pdf",plot=p,width = 6, height=4,units="in")
ggsave(filename = "fig/study1_model2_ppc_y.png",plot=p,width = 6, height=4,units="in")
ggsave(filename = "fig/study1_model2_ppc_y.eps",plot=p,width = 6, height=4,units="in")

Manuscript Table and Figures

Table

# print to xtable
print(
  xtable(
    model.fit$BUGSoutput$summary,
    caption = c("study1 Model 2 posterior distribution summary")
    ,align = "lrrrrrrrrr"
  ),
  include.rownames=T,
  booktabs=T
)

% latex table generated in R 4.0.5 by xtable 1.8-4 package
% Sun Jan 16 14:07:35 2022
\begin{table}[ht]
\centering
\begin{tabular}{lrrrrrrrrr}
  \toprule
 & mean & sd & 2.5\% & 25\% & 50\% & 75\% & 97.5\% & Rhat & n.eff \\ 
  \midrule
deviance & 7467.17 & 76.52 & 7320.23 & 7415.33 & 7467.33 & 7518.65 & 7618.68 & 1.00 & 2600.00 \\ 
  icept[1] & 1.54 & 0.06 & 1.43 & 1.50 & 1.54 & 1.58 & 1.67 & 1.01 & 180.00 \\ 
  icept[2] & 1.55 & 0.06 & 1.44 & 1.51 & 1.55 & 1.59 & 1.69 & 1.03 & 99.00 \\ 
  icept[3] & 1.65 & 0.08 & 1.49 & 1.60 & 1.66 & 1.71 & 1.81 & 1.01 & 200.00 \\ 
  icept[4] & 1.57 & 0.06 & 1.46 & 1.53 & 1.57 & 1.61 & 1.70 & 1.01 & 190.00 \\ 
  icept[5] & 1.57 & 0.08 & 1.42 & 1.52 & 1.57 & 1.62 & 1.72 & 1.02 & 140.00 \\ 
  lambda[1] & 0.44 & 0.12 & 0.22 & 0.36 & 0.43 & 0.51 & 0.71 & 1.01 & 550.00 \\ 
  lambda[2] & 0.32 & 0.11 & 0.09 & 0.25 & 0.32 & 0.39 & 0.53 & 1.04 & 240.00 \\ 
  lambda[3] & 0.61 & 0.16 & 0.32 & 0.50 & 0.60 & 0.70 & 0.95 & 1.02 & 200.00 \\ 
  lambda[4] & 0.45 & 0.12 & 0.23 & 0.36 & 0.44 & 0.52 & 0.71 & 1.01 & 670.00 \\ 
  lambda[5] & 0.57 & 0.13 & 0.31 & 0.48 & 0.56 & 0.65 & 0.83 & 1.01 & 530.00 \\ 
  lambda.std[1] & 0.40 & 0.09 & 0.21 & 0.34 & 0.40 & 0.45 & 0.58 & 1.01 & 500.00 \\ 
  lambda.std[2] & 0.30 & 0.09 & 0.09 & 0.24 & 0.30 & 0.36 & 0.47 & 1.04 & 230.00 \\ 
  lambda.std[3] & 0.51 & 0.10 & 0.30 & 0.45 & 0.51 & 0.58 & 0.69 & 1.01 & 260.00 \\ 
  lambda.std[4] & 0.40 & 0.09 & 0.23 & 0.34 & 0.40 & 0.46 & 0.58 & 1.01 & 800.00 \\ 
  lambda.std[5] & 0.48 & 0.09 & 0.30 & 0.43 & 0.49 & 0.55 & 0.64 & 1.01 & 550.00 \\ 
  prec[1] & 3.98 & 0.30 & 3.43 & 3.78 & 3.97 & 4.18 & 4.59 & 1.00 & 2500.00 \\ 
  prec[2] & 3.96 & 0.29 & 3.42 & 3.75 & 3.94 & 4.14 & 4.57 & 1.00 & 4000.00 \\ 
  prec[3] & 4.17 & 0.35 & 3.55 & 3.93 & 4.15 & 4.39 & 4.93 & 1.00 & 820.00 \\ 
  prec[4] & 4.12 & 0.31 & 3.56 & 3.91 & 4.11 & 4.32 & 4.76 & 1.00 & 3800.00 \\ 
  prec[5] & 4.71 & 0.39 & 4.01 & 4.43 & 4.69 & 4.95 & 5.54 & 1.00 & 960.00 \\ 
  prec.s & 10.39 & 1.47 & 8.06 & 9.36 & 10.21 & 11.20 & 13.80 & 1.02 & 180.00 \\ 
  reli.omega & 0.53 & 0.06 & 0.40 & 0.49 & 0.53 & 0.57 & 0.62 & 1.01 & 240.00 \\ 
  rho & 0.47 & 0.14 & 0.20 & 0.37 & 0.46 & 0.55 & 0.77 & 1.04 & 75.00 \\ 
  sigma.ts & 0.08 & 0.03 & 0.02 & 0.06 & 0.08 & 0.10 & 0.13 & 1.00 & 3300.00 \\ 
  tau[1,1] & -0.75 & 0.09 & -0.93 & -0.81 & -0.75 & -0.69 & -0.58 & 1.00 & 4000.00 \\ 
  tau[2,1] & -0.73 & 0.08 & -0.90 & -0.79 & -0.73 & -0.68 & -0.57 & 1.00 & 1400.00 \\ 
  tau[3,1] & -0.67 & 0.09 & -0.86 & -0.74 & -0.67 & -0.61 & -0.50 & 1.01 & 420.00 \\ 
  tau[4,1] & -0.50 & 0.08 & -0.66 & -0.55 & -0.50 & -0.44 & -0.33 & 1.00 & 3800.00 \\ 
  tau[5,1] & -0.83 & 0.09 & -1.00 & -0.89 & -0.83 & -0.76 & -0.66 & 1.00 & 3400.00 \\ 
  tau[1,2] & 0.82 & 0.09 & 0.65 & 0.76 & 0.82 & 0.88 & 1.00 & 1.00 & 3500.00 \\ 
  tau[2,2] & 0.85 & 0.08 & 0.69 & 0.80 & 0.85 & 0.91 & 1.02 & 1.00 & 2900.00 \\ 
  tau[3,2] & 0.88 & 0.09 & 0.70 & 0.82 & 0.88 & 0.95 & 1.06 & 1.00 & 1100.00 \\ 
  tau[4,2] & 1.02 & 0.09 & 0.84 & 0.95 & 1.01 & 1.08 & 1.20 & 1.00 & 1900.00 \\ 
  tau[5,2] & 0.87 & 0.09 & 0.69 & 0.81 & 0.87 & 0.93 & 1.06 & 1.00 & 750.00 \\ 
  theta[1] & 1.21 & 0.12 & 1.05 & 1.13 & 1.19 & 1.26 & 1.51 & 1.01 & 1100.00 \\ 
  theta[2] & 1.11 & 0.07 & 1.01 & 1.06 & 1.10 & 1.15 & 1.28 & 1.00 & 640.00 \\ 
  theta[3] & 1.39 & 0.20 & 1.10 & 1.25 & 1.36 & 1.50 & 1.90 & 1.03 & 120.00 \\ 
  theta[4] & 1.22 & 0.12 & 1.05 & 1.13 & 1.19 & 1.27 & 1.50 & 1.02 & 330.00 \\ 
  theta[5] & 1.34 & 0.16 & 1.10 & 1.23 & 1.32 & 1.42 & 1.69 & 1.01 & 530.00 \\ 
   \bottomrule
\end{tabular}
\caption{study1 Model 2 posterior distribution summary} 
\end{table}

Figure

plot.dat <- fit.mcmc %>%
  select(!c("iter", "deviance", "reli.omega", paste0("lambda[",1:5,"]")))%>%
  pivot_longer(
    cols= !c("chain"),
    names_to="variable",
    values_to="value"
  )

meas.var <- c(
        "lambda.std[1]", "theta[1]", "tau[1,1]", "tau[1,2]",
        "lambda.std[2]", "theta[2]", "tau[2,1]", "tau[2,2]",
        "lambda.std[3]", "theta[3]", "tau[3,1]", "tau[3,2]",
        "lambda.std[4]", "theta[4]", "tau[4,1]", "tau[4,2]",
        "lambda.std[5]", "theta[5]", "tau[5,1]", "tau[5,2]"
      )
plot.dat1 <- plot.dat %>%
  filter(variable %in% meas.var) %>%
  mutate(
    variable = factor(
      variable,
      levels = meas.var, ordered = T
    )
  )

spd.var <- c(
        "icept[1]", "icept[2]", "icept[3]", "icept[4]", "icept[5]", 
        "prec[1]", "prec[2]", "prec[3]", "prec[4]", "prec[5]",
        "rho", "prec.s", "sigma.ts"
      )
plot.dat2 <- plot.dat %>%
  filter(variable %in% spd.var) %>%
  mutate(
    variable = factor(
      variable,
      levels = spd.var, ordered = T
    )
  )


p1 <- ggplot(plot.dat1, aes(x=value, group=variable))+
  geom_density(adjust=2)+
  facet_wrap(variable~., scales="free_y", ncol=4) +
  labs(x="Magnitude of Parameter",
       y="Posterior Density")+
  theme_bw()+
  theme(
    panel.grid = element_blank(),
    strip.background = element_rect(fill="white"),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank() 
  )
p1

p2 <- ggplot(plot.dat2, aes(x=value, group=variable))+
  geom_density(adjust=2)+
  facet_wrap(variable~., scales="free", ncol=5) +
  labs(x="Magnitude of Parameter",
       y="Posterior Density")+
  theme_bw()+
  theme(
    panel.grid = element_blank(),
    strip.background = element_rect(fill="white"),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank() ,
    axis.text.x = element_text(size=8, angle=90, hjust=1, vjust=0.50)
  )
p2

# all as one
plot.dat <- fit.mcmc %>%
  select(!c("iter", "deviance", "reli.omega", paste0("lambda[",1:5,"]")))%>%
  pivot_longer(
    cols= !c("chain"),
    names_to="variable",
    values_to="value"
  ) %>%
  mutate(
    variable = factor(
      variable,
      # 33
      # 10x3 + 3 === horizontal page
      levels = c(
        paste0("lambda.std[",1:5,"]"), paste0("theta[",1:5,"]"),
        paste0("tau[",1:5,",1]"), paste0("tau[",1:5,",2]"),
        paste0("icept[",1:5,"]"), paste0("prec[",1:5,"]"),
        "rho", "prec.s","sigma.ts"
      ), ordered = T
    )
  )
p <- ggplot(plot.dat, aes(x=value, group=variable))+
  geom_density(adjust=2)+
  facet_wrap(variable~., scales="free", ncol=5) +
  labs(x="Magnitude of Parameter",
       y="Posterior Density")+
  theme_bw()+
  theme(
    panel.grid = element_blank(),
    strip.background = element_rect(fill="white"),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank() ,
    axis.text.x = element_text(size=7)
  )
p

ggsave(filename = "fig/study1_model2_posterior_dist.pdf",plot=p,width = 10, height=7,units="in")
ggsave(filename = "fig/study1_model2_posterior_dist.png",plot=p,width = 10, height=7,units="in")
ggsave(filename = "fig/study1_model2_posterior_dist.eps",plot=p,width = 10, height=7,units="in")

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

Simulation Study 1

Model 2 Results

R. Noah Padgett

2022-01-10