Study 4: Extroversion Data Analysis
Model 4 Results
R. Noah Padgett
2022-01-17
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:
- model4
- study4-model4-ppd
To ensure reproducibility of the results, delete the cache directory study4_model4_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 doneCall: 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.109Model 4: Full IFA with Misclassification
Model details
cat(read_file(paste0(w.d, "/code/study_4/model_4.txt")))model {
### Model
for(p in 1:N){
for(i in 1:nit){
# data model
y[p,i] ~ dbern(omega[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])
# log-RT model
dev[p,i]<-lambda[i]*(eta[p] - tau[i,1])
mu.lrt[p,i] <- icept[i] - speed[p] - rho * abs(dev[p,i])
lrt[p,i] ~ dnorm(mu.lrt[p,i], prec[i])
# MISCLASSIFICATION MODEL
for(c in 1:ncat){
# generate informative prior for misclassificaiton
# parameters based on RT
for(ct in 1:ncat){
alpha[p,i,ct,c] <- ifelse(c == ct,
ilogit(lrt[p,i]),
(1/(ncat-1))*(1-ilogit(lrt[p,i]))
)
}
# sample misclassification parameters using the informative priors
gamma[p,i,c,1:ncat] ~ ddirch(alpha[p,i,c,1:ncat])
# observed category prob (Pr(y=c))
omega[p,i, c] = gamma[p,i,c,1]*pi[p,i,1] +
gamma[p,i,c,2]*pi[p,i,2]
}
}
}
### 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)
# transformations
sigma.t = pow(prec.s, -1) + pow(sigma.ts, 2) # speed variance
cor.ts = sigma.ts/(pow(sigma.t,0.5)) # LV correlation
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)
# 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
# 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.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),
"speed"=rnorm(142),
"ystar"=matrix(c(0.7*rep(rnorm(142),10)), ncol=10),
"rho"=0.1,
"icept"=rep(0, 10),
"prec.s"=10,
"prec"=rep(4, 10),
"sigma.ts"=0.1
)
}
jags.data <- list(
y = mydata[,1:10],
lrt = mydata[,11:20],
N = nrow(mydata),
nit = 10,
ncat = 2
)
# Run model
model.fit <- R2jags::jags(
model = paste0(w.d, "/code/study_4/model_4.txt"),
parameters.to.save = jags.params,
inits = jags.inits,
data = jags.data,
n.chains = 4,
n.burnin = 5000,
n.iter = 10000
)module glm loadedCompiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 2840
Unobserved stochastic nodes: 4587
Total graph size: 43537
Initializing modelprint(model.fit, width=1000)Inference for Bugs model at "C:/Users/noahp/Documents/GitHub/Padgett-Dissertation/code/study_4/model_4.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.688 0.158 1.447 1.577 1.662 1.773 2.071 1.02 180
icept[2] 1.463 0.088 1.313 1.401 1.455 1.515 1.665 1.04 75
icept[3] 1.607 0.307 1.159 1.379 1.560 1.778 2.341 1.07 43
icept[4] 1.270 0.280 0.962 1.092 1.193 1.364 2.026 1.13 30
icept[5] 1.376 0.088 1.232 1.318 1.365 1.420 1.596 1.00 710
icept[6] 1.454 0.207 1.208 1.307 1.393 1.547 1.998 1.06 51
icept[7] 1.295 0.151 1.101 1.191 1.261 1.365 1.679 1.02 210
icept[8] 1.270 0.229 1.030 1.126 1.203 1.326 1.912 1.04 180
icept[9] 1.933 0.166 1.700 1.817 1.899 2.012 2.364 1.01 430
icept[10] 1.359 0.233 1.001 1.183 1.331 1.504 1.883 1.05 67
lambda[1] 1.744 0.712 0.528 1.244 1.658 2.191 3.361 1.06 55
lambda[2] 2.095 0.763 0.778 1.548 2.013 2.580 3.837 1.21 18
lambda[3] 2.322 0.686 1.103 1.837 2.290 2.727 3.825 1.05 81
lambda[4] 1.212 0.791 0.072 0.590 1.105 1.719 2.899 1.12 28
lambda[5] 1.268 0.686 0.204 0.783 1.161 1.670 2.803 1.06 58
lambda[6] 0.696 0.493 0.037 0.315 0.608 0.982 1.880 1.04 78
lambda[7] 0.580 0.414 0.024 0.250 0.503 0.816 1.544 1.00 1000
lambda[8] 0.642 0.461 0.029 0.279 0.564 0.901 1.748 1.02 140
lambda[9] 1.346 0.554 0.392 0.953 1.281 1.692 2.608 1.07 45
lambda[10] 1.840 0.613 0.719 1.415 1.789 2.258 3.119 1.09 43
lambda.std[1] 0.823 0.126 0.467 0.779 0.856 0.910 0.958 1.09 63
lambda.std[2] 0.871 0.093 0.614 0.840 0.896 0.932 0.968 1.18 27
lambda.std[3] 0.900 0.060 0.741 0.878 0.916 0.939 0.967 1.08 78
lambda.std[4] 0.661 0.252 0.072 0.508 0.741 0.864 0.945 1.09 41
lambda.std[5] 0.709 0.195 0.200 0.617 0.758 0.858 0.942 1.05 84
lambda.std[6] 0.495 0.245 0.037 0.300 0.519 0.701 0.883 1.03 100
lambda.std[7] 0.441 0.239 0.024 0.242 0.450 0.632 0.839 1.00 1300
lambda.std[8] 0.469 0.245 0.029 0.269 0.491 0.669 0.868 1.02 160
lambda.std[9] 0.756 0.145 0.365 0.690 0.788 0.861 0.934 1.08 75
lambda.std[10] 0.848 0.101 0.584 0.817 0.873 0.914 0.952 1.15 61
prec[1] 1.776 0.226 1.368 1.616 1.758 1.921 2.253 1.00 3900
prec[2] 3.527 0.473 2.667 3.198 3.495 3.830 4.526 1.01 470
prec[3] 4.260 0.599 3.198 3.837 4.206 4.635 5.568 1.01 240
prec[4] 2.488 0.316 1.910 2.268 2.472 2.693 3.141 1.01 500
prec[5] 2.857 0.371 2.180 2.598 2.844 3.094 3.615 1.00 2400
prec[6] 3.027 0.390 2.315 2.754 3.010 3.277 3.828 1.00 4000
prec[7] 4.923 0.655 3.740 4.467 4.883 5.334 6.331 1.00 4000
prec[8] 3.910 0.500 2.994 3.560 3.886 4.219 4.954 1.00 4000
prec[9] 2.574 0.321 1.987 2.352 2.560 2.776 3.250 1.00 3700
prec[10] 6.633 0.968 4.970 5.943 6.575 7.266 8.716 1.00 730
prec.s 11.807 2.654 7.811 9.955 11.420 13.130 18.470 1.01 350
reli.omega 0.946 0.018 0.902 0.937 0.949 0.958 0.972 1.18 22
rho 0.057 0.032 0.008 0.034 0.055 0.076 0.125 1.02 280
sigma.ts 0.114 0.058 0.002 0.075 0.114 0.152 0.233 1.01 400
tau[1,1] -1.965 0.636 -3.447 -2.313 -1.872 -1.523 -0.986 1.01 310
tau[2,1] -0.209 0.391 -0.971 -0.454 -0.220 0.034 0.615 1.04 67
tau[3,1] -3.769 0.984 -5.877 -4.449 -3.666 -3.021 -2.132 1.05 66
tau[4,1] -3.905 1.370 -7.450 -4.608 -3.585 -2.913 -2.074 1.22 18
tau[5,1] -0.444 0.427 -1.398 -0.707 -0.406 -0.145 0.303 1.02 120
tau[6,1] -4.665 1.594 -8.511 -5.586 -4.358 -3.443 -2.457 1.12 26
tau[7,1] -4.699 1.449 -8.235 -5.493 -4.442 -3.643 -2.652 1.01 250
tau[8,1] -5.223 1.560 -9.009 -6.090 -4.952 -4.075 -2.904 1.00 1100
tau[9,1] -2.542 0.637 -4.011 -2.931 -2.452 -2.079 -1.549 1.02 180
tau[10,1] -3.605 0.861 -5.509 -4.150 -3.517 -2.983 -2.141 1.01 350
theta[1] 4.549 2.806 1.279 2.548 3.750 5.799 12.298 1.05 60
theta[2] 5.970 3.562 1.605 3.395 5.054 7.657 15.721 1.23 16
theta[3] 6.862 3.419 2.216 4.376 6.242 8.436 15.627 1.04 85
theta[4] 3.094 2.552 1.005 1.348 2.221 3.956 9.401 1.18 20
theta[5] 3.078 2.229 1.042 1.613 2.348 3.789 8.859 1.07 47
theta[6] 1.727 0.997 1.001 1.099 1.369 1.965 4.535 1.07 50
theta[7] 1.508 0.658 1.001 1.062 1.253 1.665 3.384 1.01 590
theta[8] 1.624 0.839 1.001 1.078 1.318 1.811 4.057 1.04 110
theta[9] 3.118 1.706 1.153 1.909 2.641 3.864 7.799 1.08 35
theta[10] 4.760 2.390 1.517 3.003 4.201 6.097 10.729 1.07 42
deviance 3210.160 45.020 3123.530 3179.362 3209.934 3240.679 3299.242 1.02 150
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 = 993.5 and DIC = 4203.7
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_4/posterior_draws_m4.csv"))Categroy Thresholds (\(\tau\))
# tau
bayesplot::mcmc_areas(fit.mcmc, regex_pars = "tau", prob = 0.8); ggsave("fig/study4_model4_tau_dens.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_acf(fit.mcmc, regex_pars = "tau"); ggsave("fig/study4_model4_tau_acf.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_trace(fit.mcmc, regex_pars = "tau"); ggsave("fig/study4_model4_tau_trace.pdf")
Saving 7 x 5 in imageggmcmc::ggs_grb(fit.mcmc.ggs, family = "tau"); ggsave("fig/study4_model4_tau_grb.pdf")
Saving 7 x 5 in imageFactor 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_model4_lambda_dens.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_acf(fit.mcmc, regex_pars = "lambda.std"); ggsave("fig/study4_model4_lambda_acf.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_trace(fit.mcmc, regex_pars = "lambda.std"); ggsave("fig/study4_model4_lambda_trace.pdf")
Saving 7 x 5 in imageggmcmc::ggs_grb(fit.mcmc.ggs, family = "lambda.std"); ggsave("fig/study4_model4_lambda_grb.pdf")
Saving 7 x 5 in imageLatent Response Total Variance (\(\theta\))
bayesplot::mcmc_areas(fit.mcmc, regex_pars = "theta", prob = 0.8); ggsave("fig/study4_model4_theta_dens.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_acf(fit.mcmc, regex_pars = "theta"); ggsave("fig/study4_model4_theta_acf.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_trace(fit.mcmc, regex_pars = "theta"); ggsave("fig/study4_model4_theta_trace.pdf")
Saving 7 x 5 in imageggmcmc::ggs_grb(fit.mcmc.ggs, family = "theta"); ggsave("fig/study4_model4_theta_grb.pdf")
Saving 7 x 5 in imageResponse Time Intercept (\(\beta_{lrt}\))
bayesplot::mcmc_areas(fit.mcmc, regex_pars = "icept", prob = 0.8); ggsave("fig/study4_model4_icept_dens.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_acf(fit.mcmc, regex_pars = "icept"); ggsave("fig/study4_model4_icept_acf.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_trace(fit.mcmc, regex_pars = "icept"); ggsave("fig/study4_model4_icept_trace.pdf")
Saving 7 x 5 in imageggmcmc::ggs_grb(fit.mcmc.ggs, family = "icept"); ggsave("fig/study4_model4_icept_grb.pdf")
Saving 7 x 5 in imageResponse Time Precision (\(\sigma_{lrt}\))
bayesplot::mcmc_areas(fit.mcmc, pars = paste0("prec[",1:5,"]"), prob = 0.8); ggsave("fig/study4_model4_prec_dens.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_acf(fit.mcmc, pars = paste0("prec[",1:5,"]")); ggsave("fig/study4_model4_prec_acf.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_trace(fit.mcmc, pars = paste0("prec[",1:5,"]")); ggsave("fig/study4_model4_prec_trace.pdf")
Saving 7 x 5 in imageggmcmc::ggs_grb(fit.mcmc.ggs, family = "prec"); ggsave("fig/study4_model4_prec_grb.pdf")
Saving 7 x 5 in imageSpeed Factor Variance (\(\sigma_s\))
bayesplot::mcmc_areas(fit.mcmc, regex_pars = "prec.s", prob = 0.8); ggsave("fig/study4_model4_precs_dens.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_acf(fit.mcmc, regex_pars = "prec.s"); ggsave("fig/study4_model4_precs_acf.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_trace(fit.mcmc, regex_pars = "prec.s"); ggsave("fig/study4_model4_precs_trace.pdf")
Saving 7 x 5 in imageggmcmc::ggs_grb(fit.mcmc.ggs, family = "prec.s"); ggsave("fig/study4_model4_precs_grb.pdf")
Saving 7 x 5 in imageFactor Covariance (\(\sigma_{ts}\))
bayesplot::mcmc_areas(fit.mcmc, regex_pars = "sigma.ts", prob = 0.8); ggsave("fig/study4_model4_sigmats_dens.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_acf(fit.mcmc, regex_pars = "sigma.ts"); ggsave("fig/study4_model4_sigmats_acf.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_trace(fit.mcmc, regex_pars = "sigma.ts"); ggsave("fig/study4_model4_sigmats_trace.pdf")
Saving 7 x 5 in imageggmcmc::ggs_grb(fit.mcmc.ggs, family = "sigma.ts"); ggsave("fig/study4_model4_sigmats_grb.pdf")
Saving 7 x 5 in imagePID (\(\rho\))
bayesplot::mcmc_areas(fit.mcmc, regex_pars = "rho", prob = 0.8); ggsave("fig/study4_model4_rho_dens.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_acf(fit.mcmc, regex_pars = "rho"); ggsave("fig/study4_model4_rho_acf.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_trace(fit.mcmc, regex_pars = "rho"); ggsave("fig/study4_model4_rho_trace.pdf")
Saving 7 x 5 in imageggmcmc::ggs_grb(fit.mcmc.ggs, family = "rho"); ggsave("fig/study4_model4_rho_grb.pdf")
Saving 7 x 5 in imageFactor Reliability Omega (\(\omega\))
bayesplot::mcmc_areas(fit.mcmc, regex_pars = "reli.omega", prob = 0.8); ggsave("fig/study4_model4_omega_dens.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_acf(fit.mcmc, regex_pars = "reli.omega"); ggsave("fig/study4_model4_omega_acf.pdf")
Saving 7 x 5 in imagebayesplot::mcmc_trace(fit.mcmc, regex_pars = "reli.omega"); ggsave("fig/study4_model4_omega_trace.pdf")
Saving 7 x 5 in imageggmcmc::ggs_grb(fit.mcmc.ggs, family = "reli.omega"); ggsave("fig/study4_model4_omega_grb.pdf")
Saving 7 x 5 in image# extract omega posterior for results comparison
extracted_omega <- data.frame(model_4 = fit.mcmc$reli.omega)
write.csv(x=extracted_omega, file=paste0(getwd(),"/data/study_4/extracted_omega_m4.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: 2840
Unobserved stochastic nodes: 4587
Total graph size: 43537
Initializing modelfit.extra <- rjags::jags.samples(model.fit$model, variable.names = "omega", n.iter = Niter)NOTE: Stopping adaptationN <- 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$omega[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_model4_ppc_y.pdf",plot=p,width = 6, height=3,units="in")
ggsave(filename = "fig/study4_model4_ppc_y.png",plot=p,width = 6, height=3,units="in")
ggsave(filename = "fig/study4_model4_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 4 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 14:07:25 2022
\begin{table}[ht]
\centering
\begin{tabular}{lrrrrrrrrr}
\toprule
& mean & sd & 2.5\% & 25\% & 50\% & 75\% & 97.5\% & Rhat & n.eff \\
\midrule
deviance & 3210.16 & 45.02 & 3123.53 & 3179.36 & 3209.93 & 3240.68 & 3299.24 & 1.02 & 150.00 \\
icept[1] & 1.69 & 0.16 & 1.45 & 1.58 & 1.66 & 1.77 & 2.07 & 1.02 & 180.00 \\
icept[2] & 1.46 & 0.09 & 1.31 & 1.40 & 1.45 & 1.51 & 1.67 & 1.04 & 75.00 \\
icept[3] & 1.61 & 0.31 & 1.16 & 1.38 & 1.56 & 1.78 & 2.34 & 1.07 & 43.00 \\
icept[4] & 1.27 & 0.28 & 0.96 & 1.09 & 1.19 & 1.36 & 2.03 & 1.13 & 30.00 \\
icept[5] & 1.38 & 0.09 & 1.23 & 1.32 & 1.37 & 1.42 & 1.60 & 1.00 & 710.00 \\
icept[6] & 1.45 & 0.21 & 1.21 & 1.31 & 1.39 & 1.55 & 2.00 & 1.06 & 51.00 \\
icept[7] & 1.30 & 0.15 & 1.10 & 1.19 & 1.26 & 1.37 & 1.68 & 1.02 & 210.00 \\
icept[8] & 1.27 & 0.23 & 1.03 & 1.13 & 1.20 & 1.33 & 1.91 & 1.04 & 180.00 \\
icept[9] & 1.93 & 0.17 & 1.70 & 1.82 & 1.90 & 2.01 & 2.36 & 1.01 & 430.00 \\
icept[10] & 1.36 & 0.23 & 1.00 & 1.18 & 1.33 & 1.50 & 1.88 & 1.05 & 67.00 \\
lambda[1] & 1.74 & 0.71 & 0.53 & 1.24 & 1.66 & 2.19 & 3.36 & 1.06 & 55.00 \\
lambda[2] & 2.09 & 0.76 & 0.78 & 1.55 & 2.01 & 2.58 & 3.84 & 1.21 & 18.00 \\
lambda[3] & 2.32 & 0.69 & 1.10 & 1.84 & 2.29 & 2.73 & 3.82 & 1.05 & 81.00 \\
lambda[4] & 1.21 & 0.79 & 0.07 & 0.59 & 1.10 & 1.72 & 2.90 & 1.12 & 28.00 \\
lambda[5] & 1.27 & 0.69 & 0.20 & 0.78 & 1.16 & 1.67 & 2.80 & 1.06 & 58.00 \\
lambda[6] & 0.70 & 0.49 & 0.04 & 0.31 & 0.61 & 0.98 & 1.88 & 1.04 & 78.00 \\
lambda[7] & 0.58 & 0.41 & 0.02 & 0.25 & 0.50 & 0.82 & 1.54 & 1.00 & 1000.00 \\
lambda[8] & 0.64 & 0.46 & 0.03 & 0.28 & 0.56 & 0.90 & 1.75 & 1.02 & 140.00 \\
lambda[9] & 1.35 & 0.55 & 0.39 & 0.95 & 1.28 & 1.69 & 2.61 & 1.08 & 45.00 \\
lambda[10] & 1.84 & 0.61 & 0.72 & 1.42 & 1.79 & 2.26 & 3.12 & 1.09 & 43.00 \\
lambda.std[1] & 0.82 & 0.13 & 0.47 & 0.78 & 0.86 & 0.91 & 0.96 & 1.09 & 63.00 \\
lambda.std[2] & 0.87 & 0.09 & 0.61 & 0.84 & 0.90 & 0.93 & 0.97 & 1.18 & 27.00 \\
lambda.std[3] & 0.90 & 0.06 & 0.74 & 0.88 & 0.92 & 0.94 & 0.97 & 1.08 & 78.00 \\
lambda.std[4] & 0.66 & 0.25 & 0.07 & 0.51 & 0.74 & 0.86 & 0.95 & 1.09 & 41.00 \\
lambda.std[5] & 0.71 & 0.20 & 0.20 & 0.62 & 0.76 & 0.86 & 0.94 & 1.05 & 84.00 \\
lambda.std[6] & 0.50 & 0.25 & 0.04 & 0.30 & 0.52 & 0.70 & 0.88 & 1.03 & 100.00 \\
lambda.std[7] & 0.44 & 0.24 & 0.02 & 0.24 & 0.45 & 0.63 & 0.84 & 1.00 & 1300.00 \\
lambda.std[8] & 0.47 & 0.25 & 0.03 & 0.27 & 0.49 & 0.67 & 0.87 & 1.02 & 160.00 \\
lambda.std[9] & 0.76 & 0.14 & 0.36 & 0.69 & 0.79 & 0.86 & 0.93 & 1.08 & 75.00 \\
lambda.std[10] & 0.85 & 0.10 & 0.58 & 0.82 & 0.87 & 0.91 & 0.95 & 1.15 & 61.00 \\
prec[1] & 1.78 & 0.23 & 1.37 & 1.62 & 1.76 & 1.92 & 2.25 & 1.00 & 3900.00 \\
prec[2] & 3.53 & 0.47 & 2.67 & 3.20 & 3.49 & 3.83 & 4.53 & 1.01 & 470.00 \\
prec[3] & 4.26 & 0.60 & 3.20 & 3.84 & 4.21 & 4.64 & 5.57 & 1.01 & 240.00 \\
prec[4] & 2.49 & 0.32 & 1.91 & 2.27 & 2.47 & 2.69 & 3.14 & 1.01 & 500.00 \\
prec[5] & 2.86 & 0.37 & 2.18 & 2.60 & 2.84 & 3.09 & 3.62 & 1.00 & 2400.00 \\
prec[6] & 3.03 & 0.39 & 2.31 & 2.75 & 3.01 & 3.28 & 3.83 & 1.00 & 4000.00 \\
prec[7] & 4.92 & 0.65 & 3.74 & 4.47 & 4.88 & 5.33 & 6.33 & 1.00 & 4000.00 \\
prec[8] & 3.91 & 0.50 & 2.99 & 3.56 & 3.89 & 4.22 & 4.95 & 1.00 & 4000.00 \\
prec[9] & 2.57 & 0.32 & 1.99 & 2.35 & 2.56 & 2.78 & 3.25 & 1.00 & 3700.00 \\
prec[10] & 6.63 & 0.97 & 4.97 & 5.94 & 6.58 & 7.27 & 8.72 & 1.00 & 730.00 \\
prec.s & 11.81 & 2.65 & 7.81 & 9.96 & 11.42 & 13.13 & 18.47 & 1.01 & 350.00 \\
reli.omega & 0.95 & 0.02 & 0.90 & 0.94 & 0.95 & 0.96 & 0.97 & 1.18 & 22.00 \\
rho & 0.06 & 0.03 & 0.01 & 0.03 & 0.05 & 0.08 & 0.12 & 1.02 & 280.00 \\
sigma.ts & 0.11 & 0.06 & 0.00 & 0.07 & 0.11 & 0.15 & 0.23 & 1.01 & 400.00 \\
tau[1,1] & -1.97 & 0.64 & -3.45 & -2.31 & -1.87 & -1.52 & -0.99 & 1.01 & 310.00 \\
tau[2,1] & -0.21 & 0.39 & -0.97 & -0.45 & -0.22 & 0.03 & 0.61 & 1.04 & 67.00 \\
tau[3,1] & -3.77 & 0.98 & -5.88 & -4.45 & -3.67 & -3.02 & -2.13 & 1.05 & 66.00 \\
tau[4,1] & -3.91 & 1.37 & -7.45 & -4.61 & -3.58 & -2.91 & -2.07 & 1.22 & 18.00 \\
tau[5,1] & -0.44 & 0.43 & -1.40 & -0.71 & -0.41 & -0.14 & 0.30 & 1.02 & 120.00 \\
tau[6,1] & -4.67 & 1.59 & -8.51 & -5.59 & -4.36 & -3.44 & -2.46 & 1.12 & 26.00 \\
tau[7,1] & -4.70 & 1.45 & -8.24 & -5.49 & -4.44 & -3.64 & -2.65 & 1.01 & 250.00 \\
tau[8,1] & -5.22 & 1.56 & -9.01 & -6.09 & -4.95 & -4.08 & -2.90 & 1.00 & 1100.00 \\
tau[9,1] & -2.54 & 0.64 & -4.01 & -2.93 & -2.45 & -2.08 & -1.55 & 1.02 & 180.00 \\
tau[10,1] & -3.61 & 0.86 & -5.51 & -4.15 & -3.52 & -2.98 & -2.14 & 1.01 & 350.00 \\
theta[1] & 4.55 & 2.81 & 1.28 & 2.55 & 3.75 & 5.80 & 12.30 & 1.05 & 60.00 \\
theta[2] & 5.97 & 3.56 & 1.60 & 3.40 & 5.05 & 7.66 & 15.72 & 1.22 & 16.00 \\
theta[3] & 6.86 & 3.42 & 2.22 & 4.38 & 6.24 & 8.44 & 15.63 & 1.04 & 85.00 \\
theta[4] & 3.09 & 2.55 & 1.01 & 1.35 & 2.22 & 3.96 & 9.40 & 1.18 & 20.00 \\
theta[5] & 3.08 & 2.23 & 1.04 & 1.61 & 2.35 & 3.79 & 8.86 & 1.07 & 47.00 \\
theta[6] & 1.73 & 1.00 & 1.00 & 1.10 & 1.37 & 1.96 & 4.53 & 1.08 & 50.00 \\
theta[7] & 1.51 & 0.66 & 1.00 & 1.06 & 1.25 & 1.67 & 3.38 & 1.01 & 590.00 \\
theta[8] & 1.62 & 0.84 & 1.00 & 1.08 & 1.32 & 1.81 & 4.06 & 1.04 & 110.00 \\
theta[9] & 3.12 & 1.71 & 1.15 & 1.91 & 2.64 & 3.86 & 7.80 & 1.08 & 35.00 \\
theta[10] & 4.76 & 2.39 & 1.52 & 3.00 & 4.20 & 6.10 & 10.73 & 1.07 & 42.00 \\
\bottomrule
\end{tabular}
\caption{study4 Model 4 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