simPsyStudy

Simulation of Ordinal Responses for Psychometric Studies

simPsyStudy defines factorial simulation conditions and generates binary or ordinal item responses under common-factor and probit graded response model parameterizations. It records the parameters, seeds, datasets, and response summaries needed to reproduce each condition.

The multiple-condition interface lets you define named sets of theta values, loadings, and thresholds, then generate every combination with separate output folders and recorded seeds.

Development status: simPsyStudy continues development of the original AUTTT package. This is a development release; CRAN preparation is in progress.

Installation

Install the development version from GitHub:

# Install remotes once, if needed.
install.packages("remotes")

remotes::install_github("Boklauth/simPsyStudy", upgrade = "never")
library(simPsyStudy)

If you update an already loaded package, restart R before loading the newly installed version.

Current capabilities

Quick start: simulate multiple conditions

This example creates eight conditions: two sample sizes, two loading levels, and two threshold patterns. With two replications per condition, it produces 16 datasets.

1. Define the factors and theta values

There are six items and two factors. Items 1-3 measure Factor 1; items 4-6 measure Factor 2. Each theta matrix has one row per person and one column per factor.

library(simPsyStudy)

model <- list(1:3, 4:6)

set.seed(2026)
make_theta <- function(n) {
  factor1 <- rnorm(n)
  factor2 <- 0.4 * factor1 + sqrt(1 - 0.4^2) * rnorm(n)
  cbind(Factor1 = factor1, Factor2 = factor2)
}

theta_sets <- list(
  N100 = make_theta(100),
  N200 = make_theta(200)
)

This generates normal theta values with a population correlation of 0.4. Sample correlations will vary. These matrices are held fixed across the item-response replications.

2. Define loadings and thresholds

Supply one standardized FAOV (Factor Analysis with Ordinal Variables) loading per item. A threshold vector is shared by all items; alternatively, supply a matrix with one row per item. Two thresholds give three response categories.

loading_sets <- list(
  low = rep(0.4, 6),
  high = rep(0.8, 6)
)

threshold_sets <- list(
  symmetric = c(-0.5, 0.5),
  asymmetric = c(-1.2, 0.2)
)

Use FAOV thresholds, not already converted GRM intercepts. The function converts the loadings and thresholds internally.

3. Run the simulation

# Choose a new folder for each run. This path is relative to getwd().
output_dir <- file.path(getwd(), "simPsyStudy_example_run_01")

conditions <- simulate_conditions(
  model = model,
  theta_sets = theta_sets,
  loading_sets = loading_sets,
  threshold_sets = threshold_sets,
  replications = 2,
  output_dir = output_dir,
  seed = 1234,
  file_prefix = "demo"
)

conditions[, c("condition", "theta", "loading", "threshold", "status")]

The function creates the output directory and a subfolder for each condition. If the destination already contains matching condition folders or run metadata, it stops before overwriting them. To repeat the example, change simPsyStudy_example_run_01 to a new folder name.

4. Read a generated dataset

first_file <- file.path(
  output_dir, conditions$folder[1], "demo_grm_rep1.dat"
)

responses <- read.table(first_file, header = FALSE)
dim(responses)  # 100 persons, 6 items
head(responses)

The .dat files contain space-separated responses with no header or row names. This example uses categories 1, 2, and 3. Binary simulations use categories 0 and 1.

What is saved?

A failed condition stops the run and records the error in conditions.csv; completed outputs are preserved. Automatic resuming is not currently supported.

Response-generation methods: "U" and "N"

The single-condition function simdata_grm() provides two methods for drawing the random values used to assign item-response categories:

These options describe the response-generation draws, not the distribution of the supplied theta values. Both methods draw independently for every person and item and have the same target distribution. simulate_conditions() accepts either method and defaults to "U".

Relationship between simdata_faov() and simdata_grm()

The two functions provide equivalent probit response models when their parameters are placed on corresponding scales. In the FAOV formulation, the continuous response underlying item j is

\[ Y_j^* = \lambda_j\theta + \sqrt{1-\lambda_j^2}\epsilon_j, \qquad \epsilon_j \sim N(0,1), \]

and thresholds \(\tau_{jk}\) divide that response into ordered categories. The corresponding GRM cumulative probability is

\[ P(Y_j \geq k+1 \mid \theta) = \Phi(a_j\theta-d_{jk}), \]

with the conversion

\[ a_j = \frac{\lambda_j}{\sqrt{1-\lambda_j^2}}, \qquad d_{jk} = \frac{\tau_{jk}}{\sqrt{1-\lambda_j^2}}. \]

In R, the conversion is:

a <- loadings / sqrt(1 - loadings^2)
d <- sweep(thresholds, 1, sqrt(1 - loadings^2), "/")

Here, d is a positive boundary parameter. If the alternative intercept form \(\Phi(a\theta+\delta)\) is used, its intercept is \(\delta=-d\). The IRT difficulty is \(b=d/a=\tau/\lambda\) when the loading is nonzero.

The implementations are equivalent when:

• The same theta values are used in both functions.
• Theta is standardized and each factor has variance 1.
• Every item loads on exactly one factor.
• Item residuals are independent and have variance \(1-\lambda_j^2\).
• The GRM uses the probit link and parameters are converted as shown above.

With identical theta values and seeds, simdata_faov() and simdata_grm(method = "N") generate identical responses because both use standard-normal response draws. With method = "U", the individual datasets generally differ, but they follow the same response distribution.

The functions are not directly equivalent when simdata_faov(theta_matrix = NULL) generates new theta values within each replication, when factor variances differ from 1, or when the intended model includes cross-loadings, correlated item residuals, or a logistic link. If factor variances are not 1 and the underlying item responses must remain standardized, residual variances need to be based on \(1-\operatorname{diag}(\Lambda\Phi\Lambda')\) instead of \(1-\lambda^2\).

Input requirements and current scope

• The condition interface supports one or more factors. Every item must belong to exactly one factor; item indices do not have to be consecutive within a factor.
• Theta matrices must have finite numeric values, at least two rows, and nonconstant columns.
• Loadings must be strictly between -1 and 1. Thresholds must be finite and strictly increasing within each item.
• Name every input set using letters, digits, underscores, or hyphens, starting with a letter or digit. Names must be unique within each list, ignoring case.
• The GRM interface uses P(Y[j] >= k + 1 | theta) = pnorm(a[j] * theta - d[j, k]). The condition interface converts standardized FAOV loadings and thresholds to this parameterization.
• Threshold sets vary fastest, then loading sets, then theta sets. For condition i and replication r, the response-generation seed is seed + (i - 1) * replications + r. The caller’s random-number state is restored afterward.
• simulate_conditions() and simulate_condition() call the same interface. The plural name is recommended because one call can generate many conditions.

Other simulation helpers include TSK() for thresholds and response-distribution summaries, matrix-conversion functions, and multivariate normal and correlated gamma theta generators. gamma_from_normal() retains the original rejection-sampling procedure, while gamma_from_normal2() uses a Gaussian-copula transformation.

Gamma-based skewed multivariate distributions

Both gamma functions are designed to generate positively skewed multivariate data, but they use different procedures and should not be expected to return identical samples or correlations.

Use the original rejection approach as follows:

correlation <- to_cormatrix(c(0.49, 0.74, 0.87), n_dim = 3)
covariance <- to_covmatrix(correlation, c(0.95, 0.98, 1.1))

theta_rejection <- gamma_from_normal(
  shape = 5,
  rate = 5,
  mean_vec = c(-0.1, 0, 0.1),
  cov_matrix = covariance,
  size = 300,
  c = 4,
  seed_num = 45679
)

Use the Gaussian-copula approach with:

theta_copula <- gamma_from_normal2(
  shape = 5,
  rate = 5,
  mean_vec = c(-0.1, 0, 0.1),
  cov_matrix = covariance,
  size = 300,
  seed_num = 45679
)

Both functions return X, scaled.X, marginal skewness and kurtosis, and sample correlation and covariance matrices. The covariance matrix controls the dependence of the normal values used by each procedure; the sample covariance of the final skewed variables can differ from the input covariance.

Other tools and help

TSK(): thresholds and distribution summaries

TSK() converts expected category proportions into normal-theory thresholds and summarizes the shape of the corresponding ordinal response distribution. Supply the planned sample size with n and a vector of category probabilities with res_prop. The probabilities must be positive and sum to 1.

response_probabilities <- c(0.04, 0.06, 0.11, 0.37, 0.42)

item_distribution <- TSK(
  n = 300,
  res_prop = response_probabilities
)

item_distribution$thresholds
item_distribution$skew
item_distribution$Kurt1

The returned list contains the original category probabilities, cumulative probabilities, thresholds, skewness, adjusted skewness, excess kurtosis (Kurt1), ordinary kurtosis (Kurt2), and a sample-size-adjusted excess-kurtosis estimate (Kurt3). The thresholds can be repeated across items or assembled into an item-by-threshold matrix for simdata_faov() or simulate_conditions().

create_theta_mvn(): multivariate normal theta values

create_theta_mvn() generates correlated multivariate-normal latent-trait values. Supply the sample size, one mean and standard deviation per factor, the pairwise factor correlations, and a seed. Pairwise correlations are ordered (1, 2), (1, 3), ..., (m - 1, m).

theta <- create_theta_mvn(
  size = 300,
  mean_vec = c(-0.1, 0, 0.1),
  sd_vec = c(0.95, 0.98, 1.10),
  ifcor_vec = c(0.30, 0.50, 0.60),
  seed_num = 45679
)

head(theta$X)
head(theta$scaled.X)
theta$cor_mat

The returned X matrix has the requested means, standard deviations, and population correlation structure subject to sampling variation. scaled.X standardizes every generated factor and is the recommended input when loadings are interpreted as standardized FAOV loadings. The function also returns marginal skewness and kurtosis and the sample correlation and covariance matrices.

For complete function documentation, run:

?simulate_conditions
?TSK
?create_theta_mvn

See NEWS.md for development changes. Report problems through the GitHub issue tracker, including a small reproducible example and the output of sessionInfo().

Citation

To obtain the package citation, run:

cite <- format(citation("simPsyStudy"), style = "text")
cat(gsub("_", "", cite, fixed = TRUE), "\n")

Klauth B (2026). simPsyStudy: Simulation of Ordinal Responses for Psychometric Studies. R package version 1.1.8, https://github.com/Boklauth/simPsyStudy.

Author and license

Developed by Bo Klauth. simPsyStudy is distributed under the MIT license.