Main simulation
The main simulation uses the PROsetta package.
The simulation consists of 8 conditions * 20 trials = 160 tasks.
# Code: Sangdon Lim
library(PROsetta)
library(mvnfast)
library(mirt)
library(stringr)
library(parallel)
library(foreach)
library(doParallel)
n_cores <- detectCores() - 2
cl <- makeCluster(n_cores)
registerDoParallel(cl)
source("sim_get_ipar.R")
source("sim_get_data.R")
source("sim_eval_table.R")
d <- loadData(
response = "dat_dePHQ9_validation.csv",
itemmap = "imap_dePHQ9.csv",
anchor = "anchor_dePHQ9.csv",
input_dir = "data"
)
lhs_scale <- 2
rhs_scale <- 1
c1 <- seq(1.0, 0.6, -.05)
c2 <- 1:20
tasks <- expand.grid(c1, c2)
n_tasks <- nrow(tasks)
out <- foreach(
idx_task = 1:n_tasks,
.packages = c("mvnfast", "PROsetta", "mirt", "stringr")) %dopar% {
# Fetch task information
true_corr <- tasks[idx_task, 1]
idx_trial <- tasks[idx_task, 2]
# Generate thetas and response data
set.seed(idx_trial)
mu_sigma <- list()
mu_sigma$mu <- rep(0, 2)
mu_sigma$sigma <- matrix(true_corr, 2, 2)
diag(mu_sigma$sigma) <- 1
true_theta <- gen_theta(mu_sigma)
X <- gen_data(true_theta, ipar_a, ipar_d)
dX <- d
dX@response <- X
dX@response[[dX@person_id]] <- 990000 + 1:dim(X)[1]
# Unidimensional
set.seed(1)
tmp_calib <- runLinking(dX, method = "FIXEDPAR")
tmp_rsss <- runRSSS(
dX, tmp_calib, min_score = 0,
min_theta = -4.0, max_theta = 4.0, inc = 0.05
)
table_UNI <- tmp_rsss[[as.character(lhs_scale)]]
table_UNI <- table_UNI[, c(1, 4:5)]
# Equipercentile
set.seed(1)
tmp <- runEquateObserved(
dX, smooth = "loglinear",
scale_from = lhs_scale,
scale_to = rhs_scale,
type_to = "theta",
rsss = tmp_rsss,
degrees = list(6, 1),
reps = 1000
)
table_EQP <- tmp$concordance
# Calibrated projection
if (true_corr != 1) {
set.seed(1)
tmp_calib <- runLinking(dX, method = "CP")
tmp_rsss <- runRSSS(
dX, tmp_calib, min_score = 0,
min_theta = -4.0, max_theta = 4.0, inc = 0.05
)
table_CP <- tmp_rsss[[as.character(lhs_scale)]]
} else {
table_CP <- NULL
}
# Evaluate
X_lhs <- getResponse(dX, lhs_scale)
table_out <- eval_table(
X_lhs, true_theta,
lhs_scale, rhs_scale,
table_UNI, table_EQP, table_CP
)
# Write output
fn <- sprintf("out_%s_%s.csv", true_corr, idx_trial)
write.csv(table_out, file.path("results", fn), row.names = FALSE)
NULL
}
stopCluster(cl)
