Large-scale assessment data typically include numerous variables, often affected by missing values. Motivated by the challenges arising in this framework, we extend the knockoffs method for selecting predictors to settings with missing values. Our proposal relies on a preliminary phase of multiple imputation (MI) of missing values. Each imputed dataset is then processed using a suitable knockoff filter. We evaluate the performance of the proposed method through simulation studies, showing satisfactory results consistent with a recently advocated cutting-edge method. We apply the method to large-scale assessment data collected by INVALSI on test scores of Italian students in grade 5, including many background variables. This case study is challenging, as most predictors have unordered categories, a setting not considered by traditional knockoff methods. In addition, some of the key predictors are affected by missing values. Our proposal to implement the knockoffs method within an MI framework is feasible, flexible, and effective.

The article proposes a new approach to estimating the latent distribution of item response theory (IRT) using kernel density estimation (KDE), particularly the solve-the-equation (STE) algorithm developed by Sheather and Jones (1991). As with existing methods, the KDE method aims to estimate the latent distribution of IRT to reduce biases in parameter estimates when the normality assumption on the latent variable is violated. Simulation studies and an empirical example confirm the robustness of algorithmic convergence of the KDE approach, and show that the KDE approach yields parameter estimates that are more accurate than or comparable to existing methods. Unlike other approaches that require multiple model fits for smoothing parameter selection, KDE requires only a single model-fitting step, substantially reducing computation time. These findings highlight KDE as a practical and efficient method for estimating latent distributions in IRT.

Multivariate bounded discrete data arises in many fields. In the setting of dementia studies, such data are collected when individuals complete neuropsychological tests. We outline a modeling and inference procedure that can model the joint distribution conditional on baseline covariates, leveraging previous work on mixtures of experts and latent class models. Furthermore, we illustrate how the work can be extended when the outcome data are missing at random using a nested EM algorithm. The proposed model can incorporate covariate information and perform imputation and clustering. We apply our model to simulated data and an Alzheimer’s disease data set.

Several textbooks give the incorrect formula to calculate the standard error estimates for standardized regression coefficients. As a remedy, two analytic methods have been developed: (1) the delta method and (2) the covariance structure modeling method. However, neither method is applicable to compute the standard error estimates for unstandardized regression coefficients of products of Z-scores. In the literature, a nonparametric bootstrap procedure is advocated to test the significance of unstandardized regression coefficients of products of Z-scores. In this article, we propose a simple analytic approach that can produce the standard error estimates not only for standardized regression coefficients when interaction terms are not included, but also for unstandardized regression coefficients when interaction terms of Z-scores are included. Two numeric examples are used to compare our analytic approach with the existing methods, and simulation studies are conducted to further evaluate the performances of regular regression, our analytic approach, and the nonparametric bootstrap procedure at finite sample sizes. It is found that (1) regular regression performs well only when the variances of predictor variables are small, (2) our analytic approach performs well at the sample size of 200 or larger, and (3) the nonparametric bootstrap procedure performs (almost) perfectly in all conditions.

Change point analysis (CPA) detects structural shifts in a response sequence by partitioning it into segments with different statistical properties. This paper proposes three CPA approaches based on the Schwarz information criterion (SIC; hereafter SIC-CPA): response data only, response time (RT) data only, and the combination of response and RT data, to detect the prevalent test speededness in time-limit tests. To comprehensively investigate the efficiency and accuracy of the proposed approaches, six simulation studies were conducted under diverse conditions. Simulation results demonstrate that SIC-CPA can effectively enhance the power of change point detection and reduce Type I errors, while improving computational efficiency compared to the likelihood ratio and Wald tests. Moreover, the SIC-CPA combining response and RT data outperforms the SIC-CPA based solely on RTs, and the latter is substantially superior to the SIC-CPA based solely on responses. In addition, SIC-CPA accurately identifies two change points in RT patterns, corresponding to early warm-up and later test speededness. Using an iterative detect–clean–recalibrate procedure, SIC-CPA achieves more reliable Type I error control than likelihood ratio and Wald tests when item parameters are estimated from contaminated data. A real data analysis was conducted to show the application of the proposed approaches.

Polychoric correlation is often an important building block in the analysis of rating data, particularly for structural equation models. However, the commonly employed maximum likelihood (ML) estimator is highly susceptible to misspecification of the polychoric correlation model, for instance, through violations of latent normality assumptions. We propose a novel estimator that is designed to be robust against partial misspecification of the polychoric model, that is, when the model is misspecified for an unknown fraction of observations, such as careless respondents. To this end, the estimator minimizes a robust loss function based on the divergence between observed frequencies and theoretical frequencies implied by the polychoric model. In contrast to existing literature, our estimator makes no assumption on the type or degree of model misspecification. It furthermore generalizes ML estimation, is consistent as well as asymptotically normally distributed, and comes at no additional computational cost. We demonstrate the robustness and practical usefulness of our estimator in simulation studies and an empirical application on a Big Five administration. In the latter, the polychoric correlation estimates of our estimator and ML differ substantially, which, after further inspection, is likely due to the presence of careless respondents that the estimator helps identify.

Recent research shows that amortized variational inference (AVI) can be used to efficiently estimate high-dimensional latent variable models on large datasets. However, its use has remained limited to item response theory (IRT), and generalizing the approach to discrete latent variable models is not straightforward. We propose two ways to deal with this problem. In an initial simulation, we verify that these approaches can be used to estimate simple discrete latent variable models, such as latent class analysis and the generalized deterministic inputs, noisy and gate model. In these cases, AVI provides accurate parameter estimates, although the computational advantage over marginal maximum likelihood (MML) and standard variational inference (VI) is limited. We then apply the same approach to estimate mixture IRT models. In this case, AVI is computationally faster than MML estimation and standard VI. To demonstrate the practical applicability of our AVI approach, we use it to fit a seven-dimensional mixture IRT model to a narcissism inventory. Whereas quadrature-based methods cannot feasibly estimate models of this dimensionality, the efficient AVI approach even allows for computation of bootstrapped standard errors. We provide our code, along with an easy-to-use tool for fitting these models to new datasets.

Relationships among teachers are known to influence their teaching-related perceptions. We study whether and how teachers’ advising relationships (networks) are related to their perceptions of satisfaction, students, and influence over educational policies, recorded as their responses to a questionnaire (item responses). We propose a novel joint model of network and item responses (JNIRM) with correlated latent variables to understand these co-varying ties. This methodology allows the analyst to test and interpret the dependence between a network and item responses. Using the JNIRM, we discover that teachers’ advising relationships contribute to their perceptions of satisfaction and students more often than their perceptions of influence over educational policies. In addition, we observe that the complementarity principle applies in certain schools, where teachers tend to seek advice from those who are different from them. The JNIRM has the potential to improve the cost efficiency of network studies by integrating network and item response information during model estimation.

Time-varying coefficient modeling (TVCM), which represents regression coefficients as smooth functions of continuous time, provides a flexible framework for uncovering complex patterns of change in levels and associations in intensive longitudinal data. However, conventional TVCM remains limited to investigating directional effects across individuals. By introducing a TVCM formulation of the multivariate normal distribution, the present study extends TVCM to explore change in undirected associations (couplings) and variability, thereby broadening its utility for psychological research. We discuss three versions of this approach: an aggregate-level model and two hierarchical versions capturing interindividual differences in unfolding change, either via person-specific intercepts accounting for onset differences or through fully person-specific coefficient functions smoothed via partial pooling. To illustrate the proposed developments, we apply them to six weeks of intensive longitudinal data from 16 anxiety patients undergoing therapy and examine unfolding changes in the level and volatility of nervousness and threat monitoring, their coupling, as well as between-person heterogeneity in each of these. We further show how inspecting first-order derivatives of the coefficient functions supports identifying periods of stability and change. Finally, we discuss extensions incorporating person-level characteristics to explain heterogeneity in patterns of change and predict outcomes.

Cognitive diagnostic models (CDMs) provide fine-grained diagnostic feedback by modeling the relationship between latent attributes and item responses. Two key components required for CDM implementation are the Q-matrix, which links items to attributes, and the attribute hierarchy, which defines prerequisite relationships among attributes. In many practical settings, both structures are specified by experts based on cognitive theory. In this article, we propose a novel Bayesian estimation method that simultaneously learns the Q-matrix, the attribute hierarchy, and the parameters of the deterministic inputs, noisy “and” gate (DINA) model. We develop a Metropolis–Hastings within Gibbs algorithm and integrate a mini-batch strategy to improve computational efficiency. We conducted a series of simulation studies to evaluate the performance of the proposed algorithm under varying conditions, including sample size, test length, hierarchy structure, mini-batch size, and threshold settings. The results demonstrate strong recovery rates for both latent structures and item parameters, confirming the accuracy and robustness of our method. A real data application further illustrates the utility of the proposed framework in uncovering interpretable diagnostic structures. Our findings offer practical guidance for researchers seeking to implement CDMs when both the Q-matrix and attribute hierarchy are unknown.

This article investigates the problem of automatically flagging test takers who exhibit atypical responses or behaviors for further review by human experts. The objective is to develop a selection policy that maximizes the expected number of test takers correctly identified as warranting additional scrutiny while maintaining a manageable volume of reviews per test administration. The selection procedure should learn from the outcomes of the expert reviews. Since typically only a fraction of test takers are reviewed, this leads to a semi-supervised learning problem. The latter is formalized in a Bayesian setting, and the corresponding optimal selection policy is derived. Since calculating the policy and the underlying posterior distributions is computationally infeasible, a variational approximation and three heuristic selection policies are proposed. These policies are informed by properties of the optimal policy and correspond to different exploration/exploitation trade-offs. The performance of the approximate policies is assessed via numerical experiments using both synthetic and real-world data and is compared with procedures based on off-the-shelf algorithms as well as theoretical performance bounds.

While the joint modeling of item responses and response times (RTs) has received considerable attention, most existing approaches remain limited to dichotomous items and are not applicable to assessments involving polytomous or mixed-format items. To address this limitation, this article proposes a novel joint modeling framework for graded item responses and RTs. Specifically, we develop a conditional RT model given item responses and integrate it with a marginal response model based on Samejima’s graded response model, yielding a conditional joint model for graded item responses and RTs. The model is then embedded within a two-level hierarchical framework to account for the relationship between ability and speed at the population level. A key methodological contribution is the development of a stochastic approximation EM (SAEM) algorithm for estimating the proposed model, which efficiently computes its marginal maximum likelihood estimates. Simulation studies demonstrate the accurate parameter recovery of the SAEM algorithm and indicate that the proposed model outperforms the hierarchical model assuming conditional independence across various testing conditions. Finally, an empirical analysis using data from the 2022 Programme for International Student Assessment illustrates the effectiveness of the graded response–response time model in large-scale assessments.

Latent class models (LCMs) are widely used for identifying unobserved subgroups from multivariate categorical data in social sciences, with binary data as a particularly popular example. However, accurately recovering individual latent class memberships remains challenging, especially when handling high-dimensional datasets with many items. This work proposes a novel two-stage algorithm for LCMs suited for high-dimensional binary responses. Our method first initializes latent class assignments by an easy-to-implement spectral clustering algorithm, and then refines these assignments with a one-step likelihood-based update. This approach combines the computational efficiency of spectral clustering with the improved statistical accuracy of likelihood-based estimation. We establish theoretical guarantees showing that this method is minimax-optimal for latent class recovery in the statistical decision theory sense. The method also leads to exact clustering of subjects with high probability under mild conditions. As a byproduct, we propose a computationally efficient consistent estimator for the number of latent classes. Extensive experiments on both simulated data and real data validate our theoretical results and demonstrate our method’s superior performance over alternative methods.

Learning diagnosis is essential for effective education, with formative assessments shown to significantly enhance academic performance. Diagnostic classification models (DCMs) have been developed to assess students’ learning status and provide remedial instruction. However, the impact of mastery or non-mastery of specific attributes on long-term learning development remains uncertain. If certain non-mastered attributes hinder the growth of mathematical ability, early intervention becomes essential. In this study, we developed a random-effects DCM for multilevel growth curves (RDC–MGC) model to identify the specific effects of attribute mastery on individual-level mathematics ability growth. The simulation studies showed that the Bayesian estimation procedure provided appropriate parameter recovery and coverage probabilities, whereas ignoring the multilevel structure resulted in biased parameter estimates. The model was applied to arithmetic test data from second- to sixth-grade elementary school students. Diagnosis was conducted in the second grade, and the effects of mastery on mathematics ability growth from the third to sixth grades were assessed. The results showed that attribute mastery in second grade was associated with both the intercept and slope of individual ability growth, suggesting the potential importance of early-stage diagnostic information for understanding later mathematical development. Potential extensions of the proposed RDC–MGC model are also discussed.

In educational testing, inferences of ability have been mainly based on item responses, while the time taken to complete an item is often ignored. To better infer the ability, a new class of state space models, which conjointly model response time with time series of dichotomous responses, is developed. Simulations for the proposed models demonstrate that the biases of ability estimation are reduced as well as the precisions of ability estimation are improved. An empirical study is conducted using EdSphere datasets, where the two competing relationships (i.e., monotone and inverted U-shape) for the distance between ability and difficulty are investigated in modeling response times. The results of model comparison support that the inverted U-shape relationship better captures the behaviors and psychology of examinees in exams for EdSphere datasets.

Violations of the assumption of local independence are a fundamental issue in item response theory as they threaten model validity and bias the parameter estimates. For such a reason, a plethora of tests and approaches has been devised in the last 40 years to detect or to model such violations. Nonetheless, local dependence (LD) remains an open problem, with somewhat blurred boundaries due to the lack of a general framework for dealing with the different notions of dependence that have been suggested in the literature. The present contribution has a two-fold aim: On the one hand, to review and collect some of the approaches available in the literature; on the other hand, by following a unified perspective on assessment models introduced by Noventa et al. (2024, Journal of Mathematical Psychology 122, 102872) to suggest a possible systematization of some existing and some new approaches to LD. As a result, deterministic and probabilistic modeling mechanisms of LD are formalized and discussed.

In a Bayesian modeling approach for differential item functioning (DIF), the dependence structure of (non-)uniform DIF is represented by a structured covariance matrix. Item–group interactions (uniform DIF) but also item-specific person–group interactions (non-uniform DIF) are represented by additional correlations in (latent) item responses. DIF in discriminations and difficulties is modeled simultaneously across items and multiple groups in the covariance matrix, making it possible to examine (non-)uniform DIF without needing anchor item(s) or multiple-step procedures. The modeling framework is very efficient, avoids the computation of any interaction parameter, and requires only a single covariance parameter for the DIF assessment of an item parameter for any number of groups. This supports a simultaneous non-uniform DIF analysis of all items, even for small sample sizes. The proposed DIF procedure is applied to PISA data, where the advantages of the method are illustrated and compared to a multiple-group IRT analysis.

In latent space item response models (LSIRMs), subjects and items are embedded in a low-dimensional Euclidean latent space. As such, interactions among persons and/or items can be revealed that are unmodeled in conventional item response theory models. Current estimation approach for LSIRMs is a fully Bayesian procedure with Markov Chain Monte Carlo, which is, while practical, computationally challenging, hampering applied researchers to use the models in a wide range of settings. Therefore, we propose an LSIRM based on two variants of regularized joint maximum likelihood (JML) estimation: penalized JML and constrained JML. Owing to the absence of integrals in the likelihood, the JML methods allow for various models to be fit in limited amount of time. This computational speed facilitates a practical extension of LSIRMs to ordinal data, and the possibility to select the dimensionality of the latent space using cross-validation. In this study, we derive the two JML approaches and address different issues that arise when using maximum likelihood to estimate the LSIRM. We present a simulation study demonstrating acceptable parameter recovery and adequate performance of the cross-validation procedure. In addition, we estimate different binary and ordinal LSIRMs on real datasets pertaining to deductive reasoning and personality. All methods are implemented in R package ‘LSMjml’ which is available from CRAN.

ShinyItemAnalysis (SIA) is an R package and shiny application for an interactive presentation of psychometric methods and analysis of multi-item measurements in psychology, education, and social sciences in general. In this article, we present a new feature introduced in the recent version of the package, called “SIA modules,” which allows researchers and practitioners to offer new analytical methods for broader use via add-on extensions. SIA modules are designed to integrate with and build upon the SIA interactive application, enabling them to leverage the existing infrastructure for tasks, such as data uploading and processing. They can access and further use a range of outputs from various analyses, including models and datasets. Because SIA modules come in R packages (or extend existing ones), they can be bundled with their datasets, utilize object-oriented systems, or even compiled code. We illustrate the concepts using sample modules from the newly introduced SIAmodules package and other packages. After providing a general overview of building Shiny applications, we describe how to develop the SIA add-on modules with the support of the new SIAtools package. Finally, we discuss the possibilities of future development and emphasize the importance of freely available, interactive psychometric software for disseminating methodological innovations.