Regularized Mixture Rasch Model

Abstract

The mixture Rasch model is a popular mixture model for analyzing multivariate binary data. The drawback of this model is that the number of estimated parameters substantially increases with an increasing number of latent classes, which, in turn, hinders the interpretability of model parameters. This article proposes regularized estimation of the mixture Rasch model that imposes some sparsity structure on class-specific item difficulties. We illustrate the feasibility of the proposed modeling approach by means of one simulation study and two simulated case studies.

1. Introduction

In education and psychological science, multivariate data of cognitive test items such as intelligence tests are frequently analyzed. The Rasch model (RM; [1,2,3,4]) is likely the most popular statistical model in applied research for analyzing a vector of random variables $\mathit{X}=(X_1,\dots ,X_I)$ of I dichotomous item responses (i.e., $X_i\in \{0,1\}$ for $i=1,\dots ,I$). The multivariate probability distribution $\mathrm{P}(\mathit{X}=\mathit{x})$ in the RM is given as

$$\mathrm{P}(\mathit{X}=\mathit{x})=\int \prod _{i=1}^I{P}_i(x_i,\theta ;b_i)\varphi (\theta ;\mu ,\sigma )\mathrm{d}\theta \mathrm{for}\mathit{x} =(x_1,\dots x_{\mathrm{I}}) \in {\{0,1\}}^I,$$

(1)

where $P_i(1,\theta ;b_i) = \mathrm{\Psi }(\theta - b_i)$ (which is also referred to as the item response function), $P_i(0,\theta ;b_i) = 1 - P_i(1,\theta ;b_i)$ and $\mathrm{\Psi }$ denotes the logistic distribution function. Moreover, $\varphi$ is the density function of the normal distribution with mean $\mu$ and standard deviation $\sigma$. The latent variable $\theta$ can be thought of as an underlying unidimensional factor that represents multivariate dependencies of the discrete vector $\mathit{X}$. Notably, the normal distribution assumption in the RM could be weakened [5]. The item difficulties $b_i$ represent a nonlinear transformation of proportion correct values of items $X_i$. Note that an identification constraint in (1) in the estimation of the RM must be imposed. Frequently, the mean $\mu$ is set to zero, or one fixes the mean of item difficulties to zero (i.e., ${\sum }_{i=1}^I{b}_i=0$).

The mixture Rasch model (MRM; [6,7,8]) models a heterogeneous distribution for $\mathit{X}$. In a nutshell, it is assumed that the RM in each of C latent classes and the marginal distribution can be interpreted as a mixture distribution [9]. The distribution of the MRM with C latent classes is given by

$$\mathrm{P}(\mathit{X}=\mathit{x})=\sum _{c=1}^C{p}_c\int \prod _{i=1}^I{P}_i(x_i,\theta ;b_{ic})\varphi (\theta ;{\mu }_c,{\sigma }_c)\mathrm{d}\theta ,$$

(2)

where the non-negative mixture probabilities $p_c$ ($c=1,\dots ,C$) add to one. The class-specific item difficulties $b_{ic}$ in (2) indicate the difficulty (i.e., some nonlinear transformation of proportion correct values) of item $X_i$ in latent class c. The distributional differences of latent classes are captured in the mean ${\mu }_c$ and standard deviation ${\sigma }_c$. The MRM can be interpreted as a model in which subjects are allocated into one of the C latent classes. The multivariate relationships in the vector $\mathit{X}$ of items can differ across latent classes.

Like in the RM defined in Equation (1), identification constraints in the MRM defined in Equation (2) are required [10]. One can fix all class means ${\mu }_c$ to zero or set the mean of item difficulties within each of the classes to zero (i.e., ${\sum }_{i=1}^I{b}_{ic}=0$ for all $c=1,\dots ,C$). The latter constraint has the advantage that differences between item parameters across latent classes can be interpreted.

After applying a standardization of class-specific item difficulties such as the above-mentioned mean centering, differences between class-specific item difficulties can be computed. The so-called latent differential item functioning (DIF; [11,12,13]) effects qualitatively describe the distinctive behavior of latent classes at the level of items [14]. Studying these latent DIF effects is an important exploratory step in understanding the differential performance of test takers on items [15].

The MRM has been extended to polytomous item responses [16,17] and more complex item response functions $P_i$ [18,19,20,21,22,23,24]. A disadvantage of the MRM in (2) is that all item difficulties are allowed to differ across classes. In empirical data, some parameters are likely to equal each other. This is the motivation for proposing a regularized mixture Rasch model (RMRM) that presupposes that only a subset of DIF effects differs from zero. Alternatively, to put it differently, subsets of class-specific item parameters are set to be equal to each other in model estimation. This property substantially eases the model interpretation in exploratory research.

The rest of the article is structured as follows. In Section 2, we present the estimation approach for the RMRM. In Section 3, we present a simulation study for two latent classes. In Section 4 and Section 5, we present two simulated case studies that involve two or three latent classes with a particular structure of DIF effects, respectively. Finally, the paper closes with a discussion in Section 6.

2. Regularized Mixture Models

In this section, we present the estimation of the RMRM. Regularized estimation recently became popular in psychometrics, such as item response modeling [25,26], structural equation modeling [27,28], and structured latent class analysis [29,30,31]. The MRM involves C latent classes. The allocation of persons (or subjects) to latent classes is unknown. If it were known, a multiple group RM with known (i.e., manifest) group allocation would result. The investigation of known demographic groups, such as gender or language groups, is an important topic in educational measurement. Moreover, regularization techniques were recently discussed for manifest DIF detection in the RM [32,33,34,35,36,37,38].

The main idea of using regularization techniques (see [39] for an overview) for the MRM is that by subtracting an appropriate penalty term from the log-likelihood function, some simplified structure on DIF effects is imposed. Let $\mathit{X}={\left(x_{pi}\right)}_{pi}$ denote the matrix of dichotomous item responses. The marginal log-likelihood function in the MRM is given by

$$l(\mathit{b},\mathbf{\gamma };\mathit{X})=\sum _{p=1}^{N}log\left(\sum _{c=1}^C{p}_c\int \prod _{i=1}^I{P}_i(x_{pi},\theta ;b_{ic})\varphi (\theta ;{\mu }_c,{\sigma }_c)\mathrm{d}\theta \right).$$

(3)

In practice, the integration in (3) can be substituted by a summation and evaluating $\theta$ at a finite grid ${\theta }_t$ for $t=1,\dots ,T$:

$$l(\mathit{b},\mathbf{\gamma };\mathit{X})=\sum _{p=1}^{N}log\left(\sum _{c=1}^C{p}_c\sum _{t=1}^T\prod _{i=1}^I{P}_i(x_{pi},{\theta }_t;b_{ic})\omega ({\theta }_t;{\mu }_c,{\sigma }_c)\right),$$

(4)

where $\omega$ is a discrete analog of the normal density. The latent class probabilities $p_c$ can be represented by logistically transformed parameters $q_c$:

$$p_c=\frac{exp\left(q_c\right)}{{\displaystyle \sum _{d=1}^{C}exp\left(q_d\right)}}\mathrm{for}c=1,\dots ,C,$$

(5)

and one sets $q_1=0$.

Regularization techniques use penalty functions to control the variability in subsets of model parameters. For a scalar parameter x, the lasso penalty is defined as

$${𝒫}_{\mathrm{Lasso}}(x,\lambda )=\lambda \left|x\right|,$$

(6)

where $\lambda$ is a non-negative regularization parameter. It is known that the lasso penalty induces bias in estimated parameters. To circumvent this issue, the smoothly clipped absolute deviation (SCAD; [40]) penalty has been proposed. It is defined by

$${𝒫}_{\mathrm{SCAD}}(x,\lambda )=\left\{\begin{array}{cc}\lambda \left|x\right|\hfill & \mathrm{if}\left|x\right|<\lambda \hfill \\ -(x^2-{2a\lambda \left|x\right|}^2+{\lambda }^2{)\left(2(a-1)\right)}^{-1}\hfill & \mathrm{if}\lambda \le \left|x\right|\le a\lambda \hfill \\ (a+1){\lambda }^2\hfill & \mathrm{if}\left|x\right|>a\lambda \hfill \end{array}\right.$$

(7)

with $a=3.7$.

2.1. Two Alternative Approaches to Regularizing the Mixture Rasch Model

We can estimate the RMRM in two variants of applying regularization. In the first approach, we use overidentified item parameters $b_{ic}$ and the fused regularization technique [41,42]. Let $\mathit{b}$ denote the vector of all class-specific item parameters and $\mathbf{\gamma }$ the vector of all distribution parameters. The following estimation function H (i.e., the negative of a regularized likelihood function) is minimized

$$H(\mathit{b},\mathbf{\gamma };\mathit{X})=-l(\mathit{b},\mathbf{\gamma };\mathit{X})+N\sum _{i=1}^I\sum _{c=1}^{C-1}\sum _{c^{\prime }=c+1}^C{𝒫}_{\mathrm{SCAD}}(b_{ic}-b_{i{c}^{\prime }},\lambda ).$$

(8)

Note that fused regularization (8) penalizes the presence of many nonvanishing item parameter differences $b_{ic}-b_{i{c}^{\prime }}$. With a regularization parameter $\lambda =0$, differences $b_{ic}-b_{i{c}^{\prime }}$ in item difficulties are unpenalized. With increasing values of $\lambda$, the penalty contribution in the estimation function H becomes larger. Eventually, with sufficiently large $\lambda$ values, item difficulties $b_{ic}$ and $b_{i{c}^{\prime }}$ are fused; that is, they receive the same estimate.

Moreover, note that sample size N is multiplied with the penalty function in (8). We prefer this choice because optimal values of the regularization parameter $\lambda$ are less sample size dependent in this case. Moreover, optimal $\lambda$ values can be more easily compared across different sample sizes.

It should be noted that in model estimation, the regularization parameter $\lambda$ in (8) has to be fixed. In practice, the regularization parameter $\lambda$ also has to be estimated. Hence, the minimization is performed on a grid of $\lambda$ values (e.g., $\lambda =0.01,0.02,\dots ,0.50$), and the model that is selected is optimal with respect to some criterion. Typical criteria are the cross-validated log-likelihood, the Akaike information criterion (AIC), and the Bayesian information criterion (BIC) [39]. See [43] for model selection for the (nonregularized) MRM.

The second estimation approach relies on the ordinary regularization of latent DIF effects. The latent DIF effects are included by using an overidentified model with common item parameters $b_{i0}$ and latent DIF effects $e_{ic}$ by relying on the decomposition.

$$b_{ic}=b_{i0}+e_{ic}.$$

(9)

Note that the difference in item difficulties for classes c and $c^{\prime }$ are given as

$$b_{ic}-b_{i{c}^{\prime }}=e_{ic}-e_{i{c}^{\prime }}$$

(10)

Hence, latent DIF effects quantify differences between item difficulties across latent classes after introducing an implicit identification constraint for determining means ${\mu }_c$ of latent classes $c=1,\dots ,C$ (while fixing ${\mu }_1=0$). Using latent DIF effects in the second approach (9) instead of considering the regularization of differences in item difficulties in (8) possesses advantages if the focus of the analysis lies in the detection and assessment of latent DIF effects.

The estimation function based on the decomposition (9) is defined by

$$H({\mathit{b}}_0,\mathit{e},\mathbf{\gamma };\mathit{X})=-l({\mathit{b}}_0,\mathit{e},\mathbf{\gamma };\mathit{X})+N\sum _{i=1}^I\sum _{c=1}^C{𝒫}_{\mathrm{SCAD}}(e_{ic},\lambda ),$$

(11)

where ${\mathit{b}}_0$ denotes the vector of all common item parameters $b_{i0}$.

The special case of two latent classes in the MRM requires further attention. In this case, only one DIF effect $e_i$ must be included in the model by relying on the decomposition

$$\begin{array}{ccc}{b}_{i1}& =& b_{i0}-e_i/2\\ b_{i2}& =& b_{i0}+e_i/2\end{array}.$$

(12)

In general, fused regularization will impose a bit more structured solutions for more than three latent classes if there are clusters of latent classes with the same DIF effect at the level of each item. In contrast, SCAD regularization (11) only presupposes one item-specific cluster of latent classes with zero DIF effects, and all other DIF effects differ from zero but do not potentially merge into another cluster of latent classes with the same DIF effect. Whether the more general structure of fused regularization is advantageous in Rasch mixture with at least four classes is an empirical question in concrete applications.

2.2. Estimation

The regularized likelihood functions can be optimized using marginal maximum likelihood estimation and the expectation maximization (EM) algorithm [26,31,44]. The EM algorithm alternates between the E-step and the M-step. The E-step computation is identical to the estimation in nonregularized item response models. In the M-step, the maximization of the regularized expected log-likelihood function involving expected counts is carried out. The difference in regularized estimation is that the optimization function becomes nondifferentiable because the SCAD penalty is nondifferentiable. The optimization of nondifferentiable optimization can be performed using gradient descent [39] approaches or by replacing the nondifferentiable optimization functions with differentiable approximating functions [31,42,45,46]. In our experience, the latter approach is quite satisfactory in applications.

As usually encountered in mixture models, the maximum likelihood optimization function is often prone to local optima. Hence, it is recommended to estimate the RMRM with a sufficiently large number of random starting values to ensure that the estimated solution corresponds to the global optimum of the likelihood function (see [47]).

The sketched EM algorithm can be practically implemented in the general estimation function xxirt() in the R package sirt [48]. This function is used in the simulation study and the two case studies in this paper.

2.3. Computation of Standard Errors

The computation of standard errors in regularized ML estimation is an active area of research [39]. In the simulation and case studies in this article, standard errors are computed based on nonparametric bootstrap [49]. The estimated model parameter of interest $\mathbf{\gamma }$ depends on a data-driven determined regularization parameter ${\widehat{\lambda }}_{\mathrm{opt}}$ that is determined by the AIC or the BIC criterion.

The bootstrap can either determine the optimal regularization parameter or one applies regularized ML using the optimal ${\widehat{\lambda }}_{\mathrm{opt}}$ parameter obtained from the original sample. Typically, the former introduces additional variability. In a preliminary analysis in Simulated Case Study 2, it turned out that the average chosen $\lambda$ parameter in bootstrap samples was substantially larger than the regularization parameter ${\widehat{\lambda }}_{\mathrm{opt}}$ from the original sample. For this reason, we only report results of standard errors in bootstrap samples that use the fixed regularization parameter ${\widehat{\lambda }}_{\mathrm{opt}}$.

Furthermore, it is vital to implement a test of statistical significance for regularized latent DIF effects $e_i$ or differences in class-specific item difficulties. It has been suggested to report the proportion of bootstrap samples $p_{\mathrm{boot}}$ in which a regularized DIF effect was estimated equal to zero [39]. Values of $p_{\mathrm{boot}}$ that are sufficiently close to zero indicate latent DIF effects $e_i$ that significantly differ from zero.

3. Simulation Study 1: Simulation Study Involving Two Latent Classes

In this section, results from a simulation study of an RMRM with two latent classes are presented.

3.1. Method

The simulated datasets consisted of $I=20$ items with two latent classes that followed an MRM. The class-specific item difficulties $b_{ic}$ were decomposed into common item difficulties $b_{i0}$ and DIF effects $e_i$ according to Equation (12). The common item difficulties of the 20 items had equidistant values between $-2.0$ and $2.0$.

Four out of twenty items had DIF effects that differed from zero. Items 6, 8, and 17 had a positive DIF effect $\delta$, while item 11 had a negative DIF effect $-\delta$. In the simulation study, the size of the DIF effect $\delta$ was either 0.5 or 1.0.

For identification, the mean ${\mu }_1$ of the first latent class was set to zero. The standard deviation ${\sigma }_1$ of the first latent class was set to 1.0. For the second latent class, ${\mu }_2=0.5$ and ${\sigma }_2=0.8$ were chosen throughout the simulation. The class probabilities were fixed to $p_1=0.7$ and $p_2=0.3$.

Moreover, we varied the sample size N in the simulation. We chose sample sizes of 1000, 2500, and 5000 to cover a range of moderate to large sample sizes.

To avoid label switching issues in estimating the RMRM, we utilized a weak prior distribution on the logistically transformed the second latent class probability $q_2$ (i.e., $p_2=\mathrm{\Psi }\left(q_2\right)$). The prior $\pi \left(q_2\right)$ was chosen as the normal distribution $\mathrm{N}(-0.7,0.4)$, meaning that the second class was the smaller one to avoid label switching issues when estimating the mixture Rasch model. The model parameters ${\mathit{b}}_0$ (i.e., common item difficulties $b_{i0}$ for $i=1,\dots ,I$), $\mathit{e}$ (i.e., all DIF effects $e_i$ for $i=1,\dots ,I$), and $\mathbf{\gamma }$ (i.e., ${\sigma }_1$, ${\mu }_2$, ${\sigma }_2$, and $q_2$) were obtained by minimizing the penalized likelihood function:

$$H({\mathit{b}}_0,\mathit{e},\mathbf{\gamma };\mathit{X})=-l({\mathit{b}}_0,\mathit{e},\mathbf{\gamma };\mathit{X})+N\sum _{i=1}^I{𝒫}_{\mathrm{SCAD}}(e_i,\lambda )-\pi \left(q_2\right).$$

(13)

In total, 5000 replications were simulated in the 2 (DIF effects) × 3 (sample size) $=6$ conditions.

The following values of the regularization parameter $\lambda$ were chosen in a decreasing order while using the obtained estimates from the previous estimation as starting values: 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.48, 0.46, 0.44, 0.42, 0.40, 0.38, 0.36, 0.34, 0.32, 0.30, 0.29, 0.28, 0.27, 0.26, 0.25, 0.24, 0.23, 0.22, 0.21, 0.20, 0.19, 0.18, 0.17, 0.16, 0.15, 0.14, 0.13, 0.12, 0.11, 0.10, 0.09, 0.08, 0.07, 0.06, 0.05, 0.04, 0.03, 0.02, 0.01, 0.009, 0.008, 0.007, 0.006, and 0.005. DIF effect estimates $e_i$ were considered as zero if their absolute values did not exceed 0.02. In the Results Section 3.2, we report model estimates for fixed $\lambda$ values of 0.05, 0.10, and 0.15, as well as parameter estimates resulting from models with minimum AIC or BIC values. The performance of parameter estimates was assessed by bias and root mean square error (RMSE).

The whole simulation was conducted in the R software [50] using the xxirt() function in the R package sirt [48]. The code for the data simulation and model estimation can be found at https://osf.io/wrs5k/ (accessed on 30 September 2022).

3.2. Results

In

Table 1. Simulation Study 1: Average number of detected DIF effects $e_i$.

		Choice of $\mathit{\lambda }$
$\left|\mathbf{DIF}\right|$	N	AIC	BIC	0.05$\phantom{\mathit{u}}$	0.1	0.15
0.5	1000	6.4	0.4	12.5	9.1	5.3
	2500	5.3	0.2	8.6	5.7	3.4
	5000	5.1	0.3	7.7	4.4	2.0
1	1000	6.9	0.9	13.0	9.4	5.6
	2500	6.7	2.5	9.8	6.4	4.4
	5000	6.3	3.9	7.9	4.9	4.0

Note. |DIF| = absolute value of DIF effects e_i; N = sample size; λ = regularization parameter. Note that 4 out of 20 items had DIF effects different from zero.

, the average number of detected DIF effects (i.e., estimated to be different from zero) is presented. Four out of twenty items had DIF effects different from zero. Interestingly, the number of detected DIF effects was substantially underestimated if the BIC was used as the regularization parameter selection criterion, except in the condition of large DIF effects of 1 (i.e., $\left|\mathrm{DIF}\right|=1$) and a large sample size of $N=5000$ with an average of 3.9 detected DIF effects. In particular, model selection based on BIC had worse performance in the case of a small DIF effect of $\left|\mathrm{DIF}\right|=0.5$. In contrast, the model estimated with the AIC selection criterion had, on average, 5 to 7 detected DIF effects, which was slightly higher than the number of true DIF effects.

If the RMRM is estimated with a fixed regularization parameter $\lambda$, the average number of detected DIF effects decreases with increasing sample sizes. Overall, model selection based on AIC might be preferred over BIC if it is more critical not to detect true DIF effects.

Table 2. Simulation Study 1: Average type I error rates for items with no DIF effects and average power rates for items with DIF effects.

		Type I Error Rate for Non-DIF Effects ${\mathit{e}}_{\mathit{i}}$					Power Rate for DIF Effects ${\mathit{e}}_{\mathit{i}}$
		Choice of $\mathit{\lambda }$					Choice of $\mathit{\lambda }$
$\left|\mathbf{DIF}\right|$	N	AIC	BIC	0.05	0.1	0.15	AIC	BIC	0.05	0.1	0.15
0.5	1000	30.7	1.6	61.1	44.1	25.3	38.3	2.7	69.0	52.3	31.1
	2500	21.3	0.5	37.1	23.2	13.3	46.8	2.7	67.4	48.6	31.2
	5000	17.5	0.4	29.6	14.8	$\phantom{1}$5.9	58.5	$\phantom{1}$5.4	74.4	52.0	27.1
1	1000	25.6	1.7	58.5	38.4	19.7	70.0	16.5	90.1	80.5	61.6
	2500	17.6	1.0	36.4	16.1	5.8	95.9	58.0	98.9	95.3	87.5
	5000	14.6	0.6	24.5	5.5	0.8	99.9	95.4	100.0	99.6	96.0

Note. |DIF| = absolute value of DIF effects e_i; N = sample size; λ = regularization parameter.

presents type I error for non-DIF effects $e_i$ (i.e., $e_i$ was zero in the simulated data) and power rates for DIF effects $e_i$ (i.e., $e_i$ had values different from zero). It turned out that type I error rates were relatively high if the AIC was used as the model selection criterion ($\mathrm{Min}=14.6$, $\mathrm{Max}=30.7$). Type I error rates for the model selection based on BIC ranged between 0.4 and 1.7. These low type I error rates for the BIC criterion come at the price that BIC has very low power for detecting DIF effects if the true DIF effect is small (i.e., $\left|\mathrm{DIF}\right|=0.5$) or in not too large sample sizes (i.e., $N=1000$). Interestingly, for the smallest sample size of $N=1000$ and small DIF effects, type I error rates and power rates were very close to each other (i.e., based on AIC: type I error rate was 30.7 and the power rate was 38.3), but power rates improve in large sample sizes or for large DIF effects.

Figure 1 displays the optimal regularization parameter ${\lambda }_{\mathrm{opt}}$ as a function of sample size, the size of DIF effects, and the chosen information criterion AIC or BIC. It can be seen that ${\lambda }_{\mathrm{opt}}$ values were generally smaller when based on AIC instead of BIC. Moreover, the optimal regularization parameter decreases with a larger sample size. It is evident that there is substantial variability in the estimated ${\lambda }_{\mathrm{opt}}$ values in repeated samples. It can be seen that the largest ${\lambda }_{\mathrm{opt}}$ value was frequently obtained for BIC (i.e., for a small DIF effect of $\left|\mathrm{DIF}\right|=0.5$ or $N=1000$). In this case, all latent DIF effects were regularized.

Table 3. Simulation Study 1: Average absolute bias (bias) and root mean square error of model parameters.

			Bias					RMSE
			Choice of $\mathit{\lambda }$					Choice of $\mathit{\lambda }$
Par	$\left|\mathbf{DIF}\right|$	N	AIC	BIC	0.05	0.1	0.15	AIC	BIC	0.05	0.1	0.15
$b_i$	0.5	1000	0.04	0.05	0.03	0.03	0.04	0.17	0.17	0.16	0.16	0.17
		2500	0.02	0.02	0.02	0.02	0.02	0.13	0.13	0.13	0.13	0.13
		5000	0.03	0.03	0.03	0.03	0.03	0.10	0.10	0.10	0.10	0.10
	1	1000	0.04	0.04	0.04	0.04	0.04	0.16	0.17	0.16	0.16	0.17
		2500	0.04	0.05	0.04	0.04	0.04	0.10	0.11	0.10	0.10	0.10
		5000	0.04	0.04	0.04	0.04	0.04	0.08	0.07	0.08	0.08	0.07
${\mu }_2$	0.5	1000	0.15	0.19	0.12	0.13	0.16	0.34	0.37	0.30	0.31	0.34
		2500	0.09	0.11	0.08	0.09	0.10	0.28	0.29	0.27	0.27	0.28
		5000	0.05	0.07	0.05	0.05	0.07	0.22	0.23	0.21	0.21	0.22
	1	1000	0.09	0.13	0.07	0.08	0.10	0.29	0.33	0.26	0.27	0.30
		2500	0.02	0.03	0.02	0.02	0.02	0.17	0.18	0.16	0.17	0.17
		5000	0.01	0.01	0.01	0.01	0.01	0.11	0.11	0.11	0.11	0.11
${\sigma }_1$	0.5	1000	0.00	0.01	0.00	0.00	0.01	0.10	0.10	0.09	0.09	0.10
		2500	0.00	0.01	0.00	0.00	0.01	0.07	0.07	0.07	0.07	0.07
		5000	0.00	0.00	0.00	0.00	0.00	0.05	0.05	0.05	0.05	0.05
	1	1000	0.01	0.02	0.00	0.00	0.01	0.09	0.10	0.09	0.09	0.09
		2500	0.00	0.00	0.00	0.00	0.00	0.05	0.06	0.05	0.05	0.05
		5000	0.01	0.01	0.01	0.01	0.00	0.04	0.04	0.04	0.04	0.04
${\sigma }_2$	0.5	1000	0.06	0.05	0.06	0.06	0.06	0.16	0.16	0.15	0.15	0.16
		2500	0.05	0.03	0.05	0.05	0.05	0.11	0.11	0.11	0.11	0.11
		5000	0.03	0.03	0.04	0.04	0.03	0.08	0.08	0.08	0.08	0.08
	1	1000	0.04	0.03	0.04	0.04	0.04	0.14	0.15	0.14	0.14	0.14
		2500	0.03	0.02	0.03	0.03	0.03	0.08	0.09	0.08	0.08	0.08
		5000	0.03	0.02	0.03	0.03	0.02	0.06	0.06	0.06	0.06	0.06
$p_2$	0.5	1000	0.03	0.03	0.03	0.03	0.03	0.04	0.03	0.04	0.04	0.04
		2500	0.03	0.03	0.03	0.03	0.03	0.04	0.03	0.04	0.04	0.04
		5000	0.03	0.03	0.03	0.03	0.03	0.03	0.03	0.03	0.03	0.04
	1	1000	0.03	0.03	0.03	0.03	0.03	0.04	0.04	0.05	0.05	0.04
		2500	0.02	0.02	0.02	0.02	0.02	0.04	0.04	0.04	0.04	0.04
		5000	0.01	0.01	0.01	0.01	0.01	0.04	0.03	0.03	0.03	0.03
$e_i$ (no DIF)	0.5	1000	0.04	0.01	0.03	0.04	0.05	0.57	0.20	0.60	0.59	0.56
		2500	0.03	0.00	0.02	0.03	0.03	0.36	0.10	0.38	0.37	0.34
		5000	0.02	0.00	0.02	0.02	0.02	0.25	0.07	0.27	0.25	0.21
	1	1000	0.05	0.01	0.04	0.04	0.05	0.46	0.19	0.50	0.48	0.45
		2500	0.01	0.00	0.01	0.02	0.01	0.23	0.09	0.26	0.23	0.18
		5000	0.01	0.00	0.01	0.01	0.00	0.15	0.05	0.17	0.11	0.05
$e_i$ (DIF)	0.5	1000	0.22	0.46	0.14	0.17	0.25	0.62	0.54	0.59	0.60	0.63
		2500	0.15	0.47	0.09	0.14	0.22	0.47	0.51	0.42	0.47	0.51
		5000	0.13	0.45	0.09	0.15	0.27	0.40	0.50	0.35	0.41	0.48
	1	1000	0.26	0.77	0.19	0.22	0.31	0.67	0.95	0.56	0.61	0.72
		2500	0.08	0.37	0.07	0.09	0.12	0.37	0.68	0.34	0.37	0.44
		5000	0.06	0.08	0.06	0.06	0.07	0.23	0.29	0.23	0.23	0.28

Note. Par = parameter group; |DIF| = absolute value of DIF effects e_i; N = sample size; λ = regularization parameter.

displays average absolute biases and RMSE for different parameters or averaged across groups of parameters. In general, bias and RMSE were reduced in larger samples and were smaller in the presence of large DIF effects than for small DIF effects. Interestingly, and in coherence with a statement in [51], bias and RMSE for model parameters can be smaller for a fixed regularization parameter (i.e., for $\lambda =0.10$) compared with model selection based on AIC or BIC. The property of thresholding parameter estimates to zero is helpful in parameter selection (i.e., detecting DIF effects) but has disadvantages for statistical frequentist properties of bias and RMSE. It remains to be investigated whether a noticeable increase in type I error rates for a fixed regularization parameter $\lambda$ is of concern in applications of the RMRM. Overall, bias and RMSE were smaller if the model selection was carried out based on AIC instead of BIC.

Table 3 only contains a number of selected values of the regularization parameter $\lambda$. In Figure 2, the RMSE of parameters ${\mu }_2$, ${\sigma }_2$, and $p_2$ are displayed and compared with the RMSE based on the optimal regularization parameter obtained from AIC or BIC. It can be seen that small fixed $\lambda$ values were competitive with optimal regularization parameters in terms of RMSE. The situation slightly differs for the $p_2$ parameter. For moderate sample sizes $N=1000$ or $N=2500$, very large $\lambda$ values near to one led to the lowest RMSE values.

Finally, Figure 3 presents the RMSE for parameter groups $b_i$ (parameters “b”), latent DIF effects $e_i$ with a true value of zero (i.e., non-DIF effects; parameters “e_nodif”), and latent DIF effects $e_i$ with a true value different from zero (parameters “e_dif”) for selected values of the regularization parameter $\lambda$. It can be beneficial for item difficulties $b_i$ in terms of RMSE if small fixed $\lambda$ values are chosen. Obviously, using a large $\lambda$ value for non-DIF effects is advantageous because these parameters would be correctly regularized. However, a large fixed $\lambda$ value comes at the price of not detecting true DIF effects. To sum up, these findings illustrate that choosing a fixed $\lambda$ value could outperform AIC- or BIC-based regularized estimation if RMSE were the statistical criterion that would drive the estimator choice.

4. Simulated Case Study 2: Illustrative Example with a Nonspeeded and a Speeded Latent Class

In a linear fixed administered test, there is a fixed order of test items. Frequently, items at later test positions are prone to position effects; that is, they are difficultly compared with the situation if they were administered at earlier test positions. Similarly, test takers can show a performance decline [52,53,54]. This means that persons show lower performance at the end of the test compared with the beginning of the test. Importantly, the extent of performance decline can vary across persons [55,56].

Performance decline can occur if the test is speeded; that is, not all test takers reach the end of the test due to low item processing, limited testing time, or a lack of motivation. MRMs have been proposed for handling speededness effects [57]. Bolt [57] proposed to use a MRM with two latent classes. The first class refers to the nonspeeded class of test takers, while the second class refers to the speeded class of test takers. The speeded class is typically characterized by increased item difficulties for items at the end of the test [57]. In this simulated case study, we assume two latent classes in the MRM, where the class-specific item difficulties are modeled as

$$\begin{array}{ccc}{b}_{i1}& =& \hfill b_{i0}\hfill \\ b_{i2}& =& \hfill b_{i0}+e_i\hfill \end{array}.$$

(14)

In the simulated dataset, we used item parameters adapted from [57]. The item difficulties are shown in Figure 4 and numerically presented in

Table 4. Simulated Case Study 2: True and estimated item parameters.

	${\mathit{b}}_{\mathit{i}}$			${\mathit{e}}_{\mathit{i}}$
Item	True	Est	SE	True	Est	${\mathit{p}}_{\mathbf{boot}}$
1	−1.4	−1.31	0.16	$\phantom{-}$0.0	$\phantom{-}$0.00	0.74
2	−0.9	−0.85	0.15	$\phantom{-}$0.0	$\phantom{-}$0.00	0.81
3	−1.6	−1.59	0.16	$\phantom{-}$0.0	$\phantom{-}$0.00	0.82
4	−1.1	−1.02	0.17	$\phantom{-}$0.0	$\phantom{-}$0.00	0.77
5	$\phantom{-}$0.3	$\phantom{-}$0.32	0.20	$\phantom{-}$0.0	$\phantom{-}$0.00	0.59
6	$\phantom{-}$0.4	$\phantom{-}$0.44	0.17	$\phantom{-}$0.0	$\phantom{-}$0.00	0.77
7	$\phantom{-}$0.4	$\phantom{-}$0.50	0.15	$\phantom{-}$0.0	$\phantom{-}$0.00	0.86
8	$\phantom{-}$0.9	$\phantom{-}$0.95	0.15	$\phantom{-}$0.0	$\phantom{-}$0.00	0.83
9	$\phantom{-}$0.5	$\phantom{-}$0.56	0.21	$\phantom{-}$0.0	$\phantom{-}$0.00	0.63
10	$\phantom{-}$0.5	$\phantom{-}$0.58	0.18	$\phantom{-}$0.0	$\phantom{-}$0.00	0.81
11	$\phantom{-}$0.9	$\phantom{-}$0.94	0.15	$\phantom{-}$0.0	$\phantom{-}$0.00	0.84
12	$\phantom{-}$0.4	$\phantom{-}$0.56	0.29	$\phantom{-}$0.0	−0.34	0.18
13	−1.6	−1.68	0.17	$\phantom{-}$0.0	$\phantom{-}$0.25	0.65
14	−0.6	−0.75	0.27	$\phantom{-}$0.0	$\phantom{-}$0.49	0.20
15	−0.6	−0.54	0.17	$\phantom{-}$0.0	$\phantom{-}$0.00	0.85
16	$\phantom{-}$0.9	$\phantom{-}$1.01	0.20	$\phantom{-}$0.0	$\phantom{-}$0.00	0.60
17	$\phantom{-}$0.4	$\phantom{-}$0.53	0.20	$\phantom{-}$0.0	−0.27	0.72
18	$\phantom{-}$0.9	$\phantom{-}$1.04	0.22	$\phantom{-}$0.0	−0.24	0.36
19	$\phantom{-}$0.5	$\phantom{-}$0.59	0.16	$\phantom{-}$0.1	$\phantom{-}$0.00	0.65
20	−0.1	$\phantom{-}$0.03	0.15	$\phantom{-}$0.3	$\phantom{-}$0.00	0.87
21	−1.9	−1.75	0.18	$\phantom{-}$0.5	$\phantom{-}$0.00	0.85
22	$\phantom{-}$0.3	$\phantom{-}$0.22	0.18	$\phantom{-}$0.4	$\phantom{-}$0.43	0.67
23	−0.9	−0.80	0.25	$\phantom{-}$0.8	$\phantom{-}$0.40	0.35
24	$\phantom{-}$0.0	$\phantom{-}$0.01	0.22	$\phantom{-}$0.7	$\phantom{-}$0.29	0.23
25	−1.2	−1.42	0.27	$\phantom{-}$0.8	$\phantom{-}$1.03	0.23
26	−0.2	−0.22	0.23	$\phantom{-}$0.6	$\phantom{-}$0.57	0.31

Note. b_i = item difficulty; e_i = DIF effect; True = true item parameters; Est = estimated item parameters; SE = standard error estimated by nonparametric bootstrap; p_boot = bootstrap probability of obtaining an estimate equal to zero.

. In total, there are 26 test items. Only items 19 to 26 were prone to speeededness effects and had DIF effects $e_i$ larger than zero, while items 1 to 18 had equal item difficulties in the two latent classes (i.e., they had no DIF effects). The nonspeeded class had a class probability of $p_1=0.75$, and the speeded class had a probability $p_2=0.25$. The means of the two classes in the MRM were ${\mu }_1=0$ and ${\mu }_2=-0.4$, respectively. Hence, the speeded class had a lower ability on average. Moreover, the standard deviations were set to ${\sigma }_1=1.1$ and ${\sigma }_2=1.4$, respectively.

A dataset of a sample size $N=6000$ was generated. We estimated an RMRM using the parameterization (14). We used the identification constraint ${\widehat{\mu }}_1=0$. The regularization parameter $\lambda$ was specified on an equidistant grid of values between 0.50 and 0.01 with decrements of 0.01. Replication material and the dataset can be found at https://osf.io/wrs5k/ (accessed on 30 September 2022).

For illustrating standard error computation, we used a nonparametric bootstrap with 500 bootstrap samples. We determined the standard error by using the robust scale parameter median absolute deviation (MAD implemented with the R function stats::mad(); [58]) of bootstrap parameter estimates to diminish the potential effect of outliers.

Figure 5 displays the AIC as a function of increasing values of the regularization parameter $\lambda$. It turned out that $\lambda =0.06$ provided the least AIC value. Therefore, we report results based on this regularization parameter.

The regularization paths for DIF effects $e_i$ are displayed in Figure 6. With increasing values of $\lambda$, fewer DIF effects were estimated as nonzero. For example, for $\lambda =0.20$, only one estimated DIF effect differed from zero.

The standard deviation of the first class was estimated as ${\widehat{\sigma }}_1=0.94$ ($\mathrm{SE}=0.10$), which somehow differed from the true value ${\sigma }_1=1.1$. The second speeded latent class had the following estimated parameters, which closely resembled the data-generating parameters: ${\widehat{p}}_2=0.36$ ($\mathrm{SE}=0.16$, true: $p_2=0.25$), ${\widehat{\mu }}_2=-0.15$ ($\mathrm{SE}=0.33$, true: ${\mu }_2=-0.4$), and ${\widehat{\sigma }}_2=1.38$ ($\mathrm{SE}=0.18$, true: ${\sigma }_2=1.40$).

The estimated item parameters are shown in Table 4. It can be seen that for the 8 DIF items, 5 DIF effects $e_i$ were correctly estimated as different from zero, while 3 DIF items had estimated DIF effects of 0. Notably, 5 non-DIF items (i.e., items 12, 13, 14, 17, and 18) had estimated DIF effects different from zero. Overall, the estimated common item difficulties ${\widehat{b}}_{0i}$ were close to the data-generating values. In accordance with findings in Simulation Study 1, the detection of DIF effects based on the BIC was less satisfactory than based on the AIC. Based on this illustrative study, it turned out that bootstrap probabilities $p_{\mathrm{boot}}$ (see Section 2.3) were substantially larger than 0.05 for true latent DIF effects that were estimated as different from zero.

5. Simulated Case Study 3: Illustrative Example Involving Three Latent Classes

In Simulated Case Study 3, we simulate data from an MRM with three latent classes. Bolt [59] presented an application in which sparse DIF effects occur. Figure 7 shows the data-generating item difficulties for the simulated dataset in this Simulated Case Study 3 that were adapted from [59]. It can be seen that many of the class-specific item difficulties were equal to each other. RMRM can be used to effectively estimate Rasch mixtures under some sparsity assumptions of DIF effects.

The simulated dataset had a sample size $N=5000$ and $I=19$ items. The data-generating item parameters can be found in

Table 5. Simulated Case Study 3: True and estimated item parameters.

	True			Est Fused Reg ${\mathit{b}}_{\mathit{ic}}$			Est Reg ${\mathit{e}}_{\mathit{ic}}$
Item	$\phantom{-}$${\mathit{b}}_{\mathit{i}1}$	$\phantom{-}$${\mathit{b}}_{\mathit{i}2}$	$\phantom{-}$${\mathit{b}}_{\mathit{i}3}$	$\phantom{-}$${\mathit{b}}_{\mathit{i}1}$	$\phantom{-}$${\mathit{b}}_{\mathit{i}2}$	$\phantom{-}$${\mathit{b}}_{\mathit{i}3}$	$\phantom{-}$${\mathit{b}}_{\mathit{i}1}$	$\phantom{-}$${\mathit{b}}_{\mathit{i}2}$	$\phantom{-}$${\mathit{b}}_{\mathit{i}3}$
1	−0.7	1.4	−0.7	−0.74	1.47	−0.74	−0.72	1.44	−0.72
2	−0.7	1.4	−0.7	−0.70	1.44	−0.70	−0.68	1.41	−0.68
3	−2.5	1.1	−2.5	−2.62	1.13	−2.19	−2.61	1.10	−2.18
4	−1.3	1.3	−1.3	−1.29	1.33	−1.29	−1.27	1.30	−1.27
5	−1.3	1.0	−1.3	−1.32	1.07	−1.32	−1.30	1.05	−1.30
6	1.1	1.1	−0.6	1.17	1.17	−0.57	1.17	1.17	−0.56
7	1.0	−1.1	−0.6	1.03	−1.21	−0.64	1.04	−1.24	−0.62
8	0.0	−1.8	−1.2	−0.05	−1.73	−1.34	−0.05	−1.76	−1.31
9	0.4	−1.2	−1.2	0.42	−1.16	−1.16	0.42	−1.16	−1.16
10	1.5	0.3	0.3	1.45	0.30	0.30	1.45	0.29	0.29
11	2.3	3.7	3.7	2.34	3.79	3.79	2.34	3.81	3.81
12	−0.9	−1.5	1.0	−0.84	−1.41	1.05	−0.83	−1.43	1.08
13	−0.9	−1.5	1.0	−0.92	−1.42	1.15	−0.91	−1.44	1.17
14	−1.2	−1.2	−1.2	−1.16	−1.16	−1.16	−1.15	−1.15	−1.15
15	−1.8	0.0	3.3	−1.70	0.19	2.89	−1.69	0.16	3.06
16	0.0	0.0	3.3	0.07	0.07	3.21	0.07	0.07	3.24
17	1.0	1.0	3.4	1.00	1.00	3.20	0.99	0.99	3.22
18	0.0	−1.2	−1.2	0.02	−1.26	−1.26	0.02	−1.26	−1.26
19	1.4	−0.4	−0.4	1.49	−0.48	−0.48	1.49	−0.63	−0.33

Note. True = true item parameters; Est Fused Reg b_ic = Estimated item parameters using fused regularization for item difficulties b_ic; Est Reg e_ic = Estimated item parameters using SCAD regularization for DIF effects e_ic; b_ic = item difficulty of item i in class c. Correctly detected DIF effects are printed in black bold font. Incorrectly detected DIF effects are printed in red bold font.

. The class-specific distribution parameters were ${\mu }_1=0$, ${\sigma }_1=1$, $p_1=0.45$ for Class 1, ${\mu }_2=0.8$, ${\sigma }_2=0.7$, $p_2=0.35$ for Class 2, and ${\mu }_3=-0.5$, ${\sigma }_3=1.2$, $p_3=0.2$ for Class 3.

We estimated the RMRM in two variants. First, we applied fused regularization to item parameter differences $b_{ic}-b_{i{c}^{\prime }}$ (see Equation (8)). Second, we used SCAD regularization for class-specific DIF effects $e_{ic}$ based on the decomposition (9) and the regularized likelihood function (11). We used the identification constraint ${\widehat{\mu }}_1=0$ in model estimation. Replication material can be found at https://osf.io/wrs5k/ (accessed on 30 September 2022). We did not carry out a bootstrap to compute standard errors because it would require some computational effort, and our primary interests were only interpretational purposes.

The optimal regularization parameter $\lambda$ was chosen using the least AIC value. It was $\lambda =0.07$ for the two estimation approaches. Overall, it turned out that the estimated model parameters were very close in the two estimation approaches. The estimated distribution parameters for fused regularization were ${\widehat{\mu }}_1=0$, ${\widehat{\sigma }}_1=1.03$, ${\widehat{p}}_1=0.53$ for Class 1, ${\widehat{\mu }}_2=0.83$, ${\widehat{\sigma }}_2=0.73$, ${\widehat{p}}_2=0.33$ for Class 2, and ${\widehat{\mu }}_3=-0.56$, ${\widehat{\sigma }}_3=1.12$, ${\widehat{p}}_3=0.14$ for Class 3. The estimated distribution parameters for SCAD regularization for DIF effects $e_{ic}$ were ${\widehat{\mu }}_1=0$, ${\widehat{\sigma }}_1=1.03$, ${\widehat{p}}_1=0.53$ for Class 1, ${\widehat{\mu }}_2=0.81$, ${\widehat{\sigma }}_2=0.74$, ${\widehat{p}}_2=0.33$ for Class 2, and ${\widehat{\mu }}_3=-0.54$, ${\widehat{\sigma }}_3=1.12$, ${\widehat{p}}_3=0.14$ for Class 3.

The estimated class-specific item difficulties are displayed in Table 5. It is evident that the pattern of true DIF effects was perfectly detected by fused and SCAD regularization. Moreover, one item (item 3) or two items (items 3 and 19) were detected to possess additional DIF effects that were not simulated for fused and SCAD regularization, respectively. Overall, the item parameter differences between the two estimation approaches were negligible.

6. Discussion

In this article, we proposed a regularized estimation approach of the mixture Rasch model. By putting a regularization penalty on differences in class-specific item difficulties or on latent DIF effects, the interpretability of latent classes in the mixture Rasch model is substantially eased. The regularization technique enables the automatic detection of latent DIF effects and provides a parsimonious model selection.

In the simulation study that involves two latent classes, model selection based on AIC tended to outperform model selection based on BIC. With AIC, there is a tendency to overestimate the number of DIF effects. At the same time, model selection based on BIC substantially underestimates the number of DIF effects. Which of the two criteria should be used in practice depends on the choice of how large type I error rates for non-DIF effects should be tolerated while guaranteeing sufficiently large power rates for the detection of DIF effects. In our view, AIC should be preferred because BIC would result in too many true DIF effects that remain undetected.

We presented two case studies to illustrate the potential of regularized mixture Rasch models. With sufficiently large sample sizes and using AIC for model selection, we successfully recovered the data-generating structure of DIF effects. Our observation that BIC should not be universally preferred over AIC was also confirmed in other research on the mixture Rasch model [43]. Moreover, we could also replicate this finding for other classes of item response models in our research [60].

We limited our simulation study only to sample sizes larger than 1000. Much smaller sample sizes might be interesting in applied research. However, we think that the maximum likelihood estimation of mixture models should involve large sample sizes (say, at least larger than 500) to ensure a sufficiently stable estimation of model parameters. Investigating the limits of applying the regularized mixture Rasch model might be an interesting topic of future research.

The computation of standard errors by nonparametric bootstrap has only been illustrated in Simulated Case Study 2. In future research, different standard error computation methods for estimating regularized mixture Rasch models might be investigated.

As with any newly proposed statistical technique, the future will tell whether the regularization approach can prove helpful in empirical applications. We think that this technique provides a means for obtaining more interpretable and less variable class-specific item parameter estimates. Likely, the regularization approach can also be applied to other classes of mixture latent variable models, such as the two- or three-parameter mixture logistic item response model or factor mixture models.

In conclusion, we believe that regularized mixture Rasch models can be used in exploratory analysis in the same way as nonregularized mixture Rasch models. We recognize the primary potential of regularization in obtaining more structured (and more stable) results if the true class-specific item difficulties follow a sparsity assumption. This assumption might not be realistic in all applications. However, one can at least include the regularized mixture Rasch model in the researcher’s toolbox for analyzing dichotomous item responses.