Mediation analysis is one of the most popularly used methods in social sciences and related areas. To estimate the indirect effect, the least-squares regression is routinely applied, which is also the most efficient when the errors are normally distributed. In practice, however, real data sets are often non-normally distributed, either heavy-tailed or skewed, so that the least-squares estimators may behave very badly. To overcome this problem, we propose a robust M-estimation for the indirect effect via a general loss function, with a main focus on the Huber loss which is more slowly varying at large values than the squared loss. We further propose a data-driven procedure to select the optimal tuning constant by minimizing the asymptotic variance of the Huber estimator, which is more robust than the least-squares estimator facing outliers and non-normal data, and more efficient than the least-absolute-deviation estimator. Simulation studies compare the finite sample performance of the Huber loss with the existing competitors in terms of the mean square error, the type I error rate, and the statistical power. Finally, the usefulness of the proposed method is also illustrated using two real data examples.

We study how the Conditioning on Added Regression Predictions (CARP) statistics from different item pairs can be aggregated into a single overall test of monotone homogeneity. As a pairwise statistic, we use the mean conditional covariance (MCC) or its standardized value ($Z$). We use three different estimates of the covariance matrix of the pairwise test statistics: (1) the covariance matrix of the MCCs, based on the sample moments; (2) the covariance matrix of the MCCs or $Z$s, based on bootstrapping; and (3) the covariance matrix of the $Z$s, equated to the identity matrix. We consider various aggregation methods, including (a) the chi-bar-square statistic; (b) the preselected standardized partial sum of pairwise statistics; (c) the product of preselected $p$-values; (d) the minimum of preselected $p$-values; and (e–h) the same statistics, but now conditioned on post-selecting only the negative values in the test sample. We study the Type 1 error rate and power of the ensuing 20 tests based on simulations. The tests with the highest power among the tests that control the Type I error rate are based on $Z$-statistics with the identity matrix: the conditional likelihood ratio test, the conditionalized product of $p$-values, the conditionalized sum of Z-values, and the preselected product of $p$-values.

Nonparametric item response theory (IRT) models consist of assumptions that restrict the joint item-score distribution. These assumptions imply stochastic ordering properties that allow ordering of respondents and items using the simple sum score and item mean score, respectively, and imply observable data properties that are useful for investigating model fit. In this paper, we investigate these properties for two-level nonparametric IRT. We introduce four two-level nonparametric IRT models. Two models pertain to respondents nested in groups: The MHM-1, useful for ordering respondents and groups, and the DMM-1, useful for ordering respondents, groups, and items. Two models pertain to groups rated by multiple respondents: The MHM-2, useful for ordering groups, and the DMM-2, useful for ordering groups and items. We define the model assumptions, derive implied stochastic ordering properties, and derive observable data properties that are useful for model fit investigation. Relations between models and properties are also presented.

Multivariate analysis using graphical models is rapidly gaining ground in psychology. In particular, Markov random field (MRF) graphical models have become popular because their graph structure reflects the conditional associations between psychological variables. Despite the fact that most psychological variables are assessed on an ordinal scale, the analysis of MRFs for ordinal variables has received little attention in the psychometric literature. To fill this gap, we present an MRF for ordinal data that so far has not been considered in network psychometrics. We present statistical methodology to test the structure of the proposed MRF, which requires us to determine the plausibility of the opposing hypotheses of conditional dependence and independence. To this end, we develop a Bayesian approach using the inclusion Bayes factor to quantify the (lack of) evidence for a given edge. We use a Bayesian variable selection approach to model the inclusion and exclusion of edges in the network, and Bayesian model averaging to compare network structures with and without the given edge. We provide an implementation in the new R package bgms, evaluate its performance in simulations, and illustrate it with empirical data.

The basic local independence model (BLIM) is a probabilistic model developed in knowledge space theory (KST). Recently, Stefanutti, de Chiusole, et al. (2020, Psychometrika 85, 684–715) proposed the polytomous local independence model (PoLIM), which is an extension of the BLIM to items with more than two response alternatives (polytomous items). In a Commentary to this paper, Chiu et al. (2023, Psychometrika 88, 656–671) claimed that (i) the BLIM is just a deterministic input noisy AND-gate (DINA) model where every item has a single skill and, as a consequence of this, (ii) the “PoLIM is simply a paraphrase of a DINA model in cognitive diagnosis (CD) for polytomous items” (p. 656). This rejoinder shows that such statements are invalid and totally misleading. Its aim is to clarify the nature of the relationship between the BLIM and the DINA, as well as that between the PoLIM and the Polytomous DINA. It builds upon formal results by Heller, et al. (2015, Psychometrika 80(4), 995–1019) on the intimate relation between KST and CD notions, and shows that the BLIM/PoLIM may be conceived as marginal models for whole classes of CD models.

Survey questionnaires are commonly used by psychologists and social scientists to measure various latent traits of study subjects. Various causal inference methods such as the potential outcome framework and structural equation models have been used to infer causal effects. However, the majority of these methods assume the knowledge of true causal structure, which is unknown for many applications in psychological and social sciences. This calls for alternative causal approaches for analyzing such questionnaire data. Bayesian networks are a promising option as they do not require causal structure to be known a priori but learn it objectively from data. Although we have seen some recent successes of using Bayesian networks to discover causality for psychological questionnaire data, their techniques tend to suffer from causal non-identifiability with observational data. In this paper, we propose the use of a state-of-the-art Bayesian network that is proven to be fully identifiable for observational ordinal data. We develop a causal structure learning algorithm based on an asymptotically justified BIC score function, a hill-climbing search strategy, and the bootstrapping technique, which is able to not only identify a unique causal structure but also quantify the associated uncertainty. Using simulation studies, we demonstrate the power of the proposed learning algorithm by comparing it with alternative Bayesian network methods. For illustration, we consider a dataset from a psychological study of the functional relationships among the symptoms of obsessive-compulsive disorder and depression. Without any prior knowledge, the proposed algorithm reveals some plausible causal relationships. This paper is accompanied by a user-friendly open-source R package OrdCD on CRAN.

In modeling change over time, developmental theories often emphasize meaningful quantities like peaks, inflections, timing, and tempo. However, longitudinal analyses typically rely on simple polynomial models that estimate powered terms of time in a linear, additive form which are disconnected from these meaningful quantities. While these linear parameterizations are computationally efficient and produce stable results, the quantities estimated in these models are difficult to directly connect to theoretical hypotheses. To address this disconnect between estimation and theory development, I propose several approaches for linear estimation with nonlinear inference (LENI), a framework that transforms results from stable, easily-estimated linear models into nonlinear estimates which align with theoretical quantities of interest through a set of principled transformation functions. I first outline derivations for the interpretable nonlinear parameters, and transform the results of the corresponding linear model—including fixed and random effects as well as conditional covariates effects —into the results we would have obtained by fitting a nonlinear version of the model. I conclude by summarizing a linearized structural equation model approach which can flexibly accommodate any known nonlinear target function within a linearly-estimable framework. I conclude with recommendations for applied researchers and directions for fruitful future work in this area.

In this paper, we propose a novel and highly effective variational Bayesian expectation maximization-maximization (VBEM-M) inference method for log-linear cognitive diagnostic model (CDM). In the implementation of the variational Bayesian approach for the saturated log-linear CDM, the conditional variational posteriors of the parameters that need to be derived are in the same distributional family as the priors, the VBEM-M algorithm overcomes this problem. Our algorithm can directly estimate the item parameters and the latent attribute-mastery pattern simultaneously. In contrast, Yamaguchi and Okada’s (2020a) variational Bayesian algorithm requires a transformation step to obtain the item parameters for the log-linear cognitive diagnostic model (LCDM). We conducted multiple simulation studies to assess the performance of the VBEM-M algorithm in terms of parameter recovery, execution time, and convergence rate. Furthermore, we conducted a series of comparative studies on the accuracy of parameter estimation for the DINA model and the saturated LCDM, focusing on the VBEM-M, VB, expectation-maximization, and Markov chain Monte Carlo algorithms. The results indicated that our method can obtain more stable and accurate estimates, especially for the small sample sizes. Finally, we demonstrated the utility of the proposed algorithm using two real datasets.

Multidimensional item response theory (MIRT) offers psychometric models for various data settings, most popularly for dichotomous and polytomous data. Less attention has been devoted to count responses. A recent growth in interest in count item response models (CIRM)—perhaps sparked by increased occurrence of psychometric count data, e.g., in the form of process data, clinical symptom frequency, number of ideas or errors in cognitive ability assessment—has focused on unidimensional models. Some recent unidimensional CIRMs rely on the Conway–Maxwell–Poisson distribution as the conditional response distribution which allows conditionally over-, under-, and equidispersed responses. In this article, we generalize to the multidimensional case, introducing the Multidimensional Two-Parameter Conway–Maxwell–Poisson Model (M2PCMPM). Using the expectation-maximization (EM) algorithm, we develop marginal maximum likelihood estimation methods, primarily for exploratory M2PCMPMs. The resulting discrimination matrices are rotationally indeterminate. Recently, regularization of the discrimination matrix has been used to obtain a simple structure (i.e., a sparse solution) for dichotomous and polytomous data. For count data, we also (1) rotate or (2) regularize the discrimination matrix. We develop an EM algorithm with lasso ($\ell _1$) regularization for the M2PCMPM and compare (1) and (2) in a simulation study. We illustrate the proposed model with an empirical example using intelligence test data.

We consider a test which allows students to attempt a multiple-choice question multiple times (i.e., multiple attempts). The most extreme form of multiple attempts is the answer-until-correct (AUC) procedure. Previous research has demonstrated that multiple-attempt procedures such as AUC could potentially enhance learning and increase reliability. However, for multiple-choice items, guessing is still non-ignorable. Traditional models of sequential item response theory (SIRT) could model multiple-attempt procedures but fail to take guessing into account. The purpose of this study is to propose SIRT models for multiple-choice, multiple-attempt items (SIRT-MM). First, we defined a family of SIRT-MM models to account for the idiosyncrasies of items, answer options, and examinee behavior. We also explained how these models could improve person parameter estimates by taking into account partial (mis)information of examinees. Second, we conducted model comparisons between the SIRT-MM models, the graded response model, and the nominal response model. Third, we discussed how the item and person parameters can be estimated, and evaluated item and person parameter recovery of SIRT-MM models under different conditions through a simulation study. Finally, we applied the SIRT-MM models to a real dataset and demonstrated their utility through model selection, person parameter recovery, and information functions.

This study extends the loss function-based parameter estimation method for diagnostic classification models proposed by Ma, de la Torre, et al. (2023, Psychometrika) to consider prior knowledge and uncertainty of sampling. To this end, we integrate the loss function-based estimation method with the generalized Bayesian method. We establish the consistency of attribute mastery patterns of the proposed generalized Bayesian method. The proposed generalized Bayesian method is compared in a simulation study and found to be superior to the previous nonparametric diagnostic classification method—a special case of the loss function-based method. Moreover, the proposed method is applied to real data and compared with previous parametric and nonparametric estimation methods. Finally, practical guidelines for the proposed method and future research directions are discussed.

For classroom teaching and learning, classifying students’ skills into more than two categories (e.g., no, basic, and advanced masteries) is more instructionally relevant. Such classifications require polytomous attributes, for which most existing cognitive diagnosis models (CDMs) are inapplicable. This paper proposes the saturated polytomous cognitive diagnosis model (sp-CDM), a general model that subsumes existing CDMs for polytomous attributes as special cases. The generalization is shown by mathematically illustrating the relationships between the proposed and existing CDMs. Moreover, algorithms to estimate the proposed model is proposed. A simulation study is conducted to evaluate the parameter recovery of the sp-CDM using the proposed estimation algorithms, as well as to illustrate the consequences of improperly fitting constrained or unnecessarily complex polytomous-attribute CDMs. A real-data example involving polytomous attributes is presented to demonstrate the practical utility of the proposed model.

The family of polychoric models (PM) categories ordinal data with latent multivariate normal variables. This modeling framework is commonly used to study the association between ordinal variables, often leading to a polychoric correlation model (PCM). Moreover, PM subsumes several well-known psychometric models, such as the structural equation modeling (SEM) with ordinal data. That said, the identifiability of PM has not been addressed in the literature. Meanwhile, in recent years researchers have suggested that the latent variables underlying PM could be generalized to the family of elliptical distributions, such as the multivariate logistic and t distributions. This article concerns the identifiability of PM and PCM with latent elliptical distributions, for which we show that PM is not identifiable and PCM is identifiable. In particular, we prove the identifiability of the polychoric t correlation model based on the copula representation. We then move on to find the set of identifiability constraints of PM through an “equivalence-classes approach of identifiability,” and demonstrate its use in two applications: one concerns the identifiability of PM on Likert scales and on comparative judgment, and the other concerns the identifiability of ordinal SEM and item factor analysis. Possible implications induced by these identifiability constraints are discussed.

When analyzing scaling conditions in latent variable structural equation models (SEMs) with continuous observed variables, analysts scaling a latent variable typically set the factor loading of one indicator to one and either set its intercept to zero or the mean of its latent variable to zero. When binary and ordinal observed variables are part of SEMs, the identification and scaling choices are more varied and multifaceted. Longitudinal data further complicate this. In SEM software, such as lavaan and Mplus, fixing the underlying variables’ variances or the error variances to one are two primary scaling conventions. As demonstrated in this paper, choosing between these constraints can significantly impact longitudinal analysis, affecting model fit, degrees of freedom, and assumptions about the dynamic process and error structure. We explore alternative parameterizations and conditions of model equivalence with categorical repeated measures.

Using data from the National Longitudinal Survey of Youth 1997, we empirically explore how different parameterizations lead to varying conclusions in longitudinal categorical analysis. More specifically, we provide insights into the specifications of the autoregressive latent trajectory model and its special cases—the linear growth curve and first-order autoregressive models—for categorical repeated measures. These findings have broader implications for a wide range of longitudinal models.

In redundancy analysis (RA), the redundancy variates are interpreted in terms of the predictor variables that have the prominent redundancy loadings. Israels (1986) advocated the rotation of redundancy loadings to facilitate the interpretation of the rotated redundancy variates. In this paper, the purpose is to obtain the standard error estimates for rotated redundancy loadings that can facilitate the interpretation of the rotated redundancy variates. To this end, we modify the original RA-L model (Gu et al., 2023) and specify two modified RA-L models for orthogonal and oblique rotations, separately. On the basis of the modified RA-L models, we describe the infinitesimal jackknife (IJ) method that can produce the standard error estimates for rotated RA estimates. A simulation study is conducted to validate the standard error estimates from the IJ method, and two real examples are used to demonstrate the use of the standard error estimates for rotated redundancy loadings. Finally, we summarize the paper and provide additional remarks regarding the rotation methods and the use of numeric derivatives in the implementation of the IJ method.

Extended redundancy analysis (ERA) is a statistical approach to component-based multivariate regression modeling that explores interrelationships among multiple sets of while incorporating regression with a data-reduction technique. The extant models that utilize ERA have assumed the outcome variables with the same data type. Also, ERA models focused on estimating direct pathways only without explicitly addressing mediation effects. In this paper, ERA is extended to handle multiple mediators and mixed types of outcome variables by adopting a Bayesian framework, taking into account correlation structure among all of the outcome variables. The proposed method develops an algorithm that derives the joint posterior distribution of parameters using a Markov chain Monte Carlo algorithm. Simulations and an empirical dataset are provided to illustrate the usefulness of the proposed method.

Data harmonization is an emerging approach to strategically combining data from multiple independent studies, enabling addressing new research questions that are not answerable by a single contributing study. A fundamental psychometric challenge for data harmonization is to create commensurate measures for the constructs of interest across studies. In this study, we focus on a regularized explanatory multidimensional item response theory model (re-MIRT) for establishing measurement equivalence across instruments and studies, where regularization enables the detection of items that violate measurement invariance, also known as differential item functioning (DIF). Because the MIRT model is computationally demanding, we leverage the recently developed Gaussian Variational Expectation–Maximization (GVEM) algorithm to speed up the computation. In particular, the GVEM algorithm is extended to a more complicated and improved multi-group version with categorical covariates and Lasso penalty for re-MIRT, namely, the importance weighted GVEM with one additional maximization step (IW-GVEMM). This study aims to provide empirical evidence to support feasible uses of IW-GVEMM for re-MIRT DIF detection, providing a useful tool for integrative data analysis. Our results show that IW-GVEMM accurately estimates the model, detects DIF items, and finds a more reasonable number of DIF items in a real world dataset. The proposed method has been integrated into R package VEMIRT (https://map-lab-uw.github.io/VEMIRT).

We develop a fully nonparametric, easy-to-use, and powerful test for the missing completely at random (MCAR) assumption on the missingness mechanism of a dataset. The test compares distributions of different missing patterns on random projections in the variable space of the data. The distributional differences are measured with the Kullback-Leibler Divergence, using probability Random Forests (Malley et al., 2011). We thus refer to it as “Projected Kullback–Leibler MCAR” (PKLM) test. The use of random projections makes it applicable even if very few or no fully observed observations are available or if the number of dimensions is large. An efficient permutation approach guarantees the level for any finite sample size, resolving a major shortcoming of most other available tests. Moreover, the test can be used on both discrete and continuous data. We show empirically on a range of simulated data distributions and real datasets that our test has consistently high power and is able to avoid inflated type-I errors. Finally, we provide an R-package PKLMtest with an implementation of our test.

Despite the versatility of generalized linear mixed models in handling complex experimental designs, they often suffer from misspecification and convergence problems. This makes inference on the values of coefficients problematic. In addition, the researcher’s choice of random and fixed effects directly affects statistical inference correctness. To address these challenges, we propose a robust extension of the “two-stage summary statistics” approach using sign-flipping transformations of the score statistic in the second stage. Our approach efficiently handles within-variance structure and heteroscedasticity, ensuring accurate regression coefficient testing for 2-level hierarchical data structures. The approach is illustrated by analyzing the reduction of health issues over time for newly adopted children. The model is characterized by a binomial response with unbalanced frequencies and several categorical and continuous predictors. The proposed approach efficiently deals with critical problems related to longitudinal nonlinear models, surpassing common statistical approaches such as generalized estimating equations and generalized linear mixed models.