Robust estimation of the latent trait in graded response models

TL;DR

This paper introduces a new method to improve the accuracy of psychological assessments by reducing the impact of careless or incorrect responses in surveys with multiple-choice items.

Contribution

A novel robust estimator for the graded response model is proposed to handle response disturbances in Likert-type items.

Findings

The robust estimator significantly reduces bias in latent trait estimates under various conditions.

Stable standard errors were observed, indicating reliable performance of the estimator.

Application to real data from the Big Five Inventory-2 demonstrates practical utility.

Abstract

Aberrant responses (e.g., careless responses, miskeyed items, etc.) often contaminate psychological assessments and surveys. Previous robust estimators for dichotomous IRT models have produced more accurate latent trait estimates with data containing response disturbances. However, for widely used Likert-type items with three or more response categories, a robust estimator for estimating latent traits does not exist. We propose a robust estimator for the graded response model (GRM) that can be applied to Likert-type items. Two weighting mechanisms for downweighting “suspicious” responses are considered: the Huber and the bisquare weight functions. Simulations reveal the estimator reduces bias for various test lengths, numbers of response categories, and types of response disturbances. The reduction in bias and stable standard errors suggests that the robust estimator for the GRM is effective in counteracting the harmful effects of response disturbances and providing more accurate scores on psychological assessments. The robust estimator is then applied to data from the Big Five Inventory-2 (Ober et al., 2021) to demonstrate its use. Potential applications and implications are discussed.

Figures

Weights as a function of the residual (*black* = bisquare weights, *grey* = Huber weights). *Solid lines* represent the tuning parameters used in the study (*B*=4.0 and *H*=1.0). *Dashed lines* represent weights when using higher tuning parameters (*B*=8.0 and *H*=2.0)

Comparing the MLE of each latent trait with its robust estimate. The figure on the *left* is the bisquare-weighted robust estimate and the figure on the right is the Huber-weighted robust estimate. *Points in black* represent the 20% of subjects with aberrant data, where 33% of their responses are oppositely coded. *Points in gray* represent subjects with “nonaberrant” data. Subjects were not plotted if their MLE or robust estimate did not converge

Full text

Introduction

Likert-type items are commonly used in psychological and educational assessments (exams, tests, surveys, etc.). Item response theory (IRT) models such as the graded response model (GRM) have been developed to model the probability that a subject responds in a certain category on these items (Samejima, 1969). Given a subject’s response vector and the item parameters, one can estimate their latent trait, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ,$$\end{document} via methods such as maximum likelihood (ML) estimation. However, responses to these items may not follow the assumed IRT model if the respondent exhibits one or more types of aberrant behaviors.

In the context of high-stakes testing, a variety of aberrant behaviors can plague a subject’s responses on such assessments, such as the six examples outlined by Meijer (1996). A subject may exhibit sleeping or warm up behavior, giving less effort to or being less focused on early items but paying greater attention to items appearing later in the test. Guessing occurs when a subject of low ability guesses blindly to items of higher difficulty or items towards the end of an assessment under time pressure. Cheating can occur when a subject with low ability copies answers from a subject of high ability. Plodding describes a subject who responds very slowly and carefully to each item. Alignment errors may occur when the response on the answer sheet does not match the intended response. Finally, some subjects of high ability may be extremely creative, reinterpreting easy items and responding incorrectly to them. In the low-stakes assessment context, aberrant behaviors such as careless or insufficient effort (C/IE) responding have also received a lot of attention, in which the subject responds with reduced attention to the item (Curran, 2016; Hong et al., 2020). C/IE responding can be nonrandom (e.g., endorsing the same response category) or random (Meade & Craig, 2012).

When a response vector contains even just some aberrant responses, traditional methods such as ML estimation result in biased estimates of the latent trait (De Ayala et al., 2000; Mislevy & Bock, 1982). Spuriously high or spuriously low assessment scores may yield undesirable consequences, such as giving opportunities to subjects who are not qualified, or failing to identify subjects who need further attention in clinical settings. Beyond latent trait estimates, other statistical properties of the test may be compromised, such as decreased scale reliability (Huang et al., 2012), poor factor model fit (Woods, 2006), lower predictive validity (Hong & Cheng, 2019a), and biased structural parameters (Oshima, 1994; Wise et al., 2004). If the aberrant subjects are known, steps can be taken (e.g., retesting or removing these subjects from the data) in order to eliminate the bias produced by the response disturbances. However, the identity of aberrant subjects is often unknown. Therefore, in order to minimize the bias, it is important to detect aberrant responses and then reduce or eliminate their negative impact.

Several methods have been developed in order to address aberrant responses to polytomous items. However, much of the previous literature focuses on fully removing cases marked as aberrant. Given the controversy of data removal, and to avoid losing information, we seek to downweight responses according to their degree of aberrancy instead of removing them. Hence, a robust estimator for the GRM will be derived to obtain more accurate latent trait estimates. With two examples of aberrant behavior, we aim to demonstrate that these robust estimates are less biased than the maximum likelihood estimates (MLE) in the presence of response disturbances.

Statistical methods to handle response disturbances

In the IRT literature, which focuses more on achievement or ability testing, a variety of methods have been developed as well to address the issue of aberrant responses. These methods generally handle response disturbances in one of two ways: modeling aberrant behavior or detecting aberrant behavior. In the former approach, traditional IRT models (e.g., the 1-PL, 2-PL, and 3-PL IRT models) are replaced with other models that accommodate the aberrant behavior. For example, speededness may be modeled as a reduction in ability, either abruptly (Yamamoto, 1982, 1995) or gradually (Goegebeur et al., 2008). Or, the dimension representing ability may be accompanied by another dimension to represent the aberrant behavior, such as an auxiliary dimension to represent speededness (Van Der Linden et al., 1999; Van Der Linden & Xiong, 2013), or an additional factor representing careless response styles (Wetzel & Carstensen, 2017). In some cases, careless response behaviors have been modeled through a sequential decision-making process (Böckenholt, 2017) and flexible item parameter estimates (Falk & Cai, 2016). When the data contain both normal and aberrant responses, mixture models can be used to account for a combination of “solution” or normal response behavior and aberrant response behaviors, including speeding (Bolt et al., 2002; Wang & Xu, 2015) or cheating (Wang et al., 2018).

However, modeling aberrant behaviors requires the researcher to know which aberrant behavior occurs in the data and by the mechanism in which it manifests, an unrealistic assumption in practice since the nature of the aberrant response is rarely known (Meade & Craig, 2012). On the other hand, if the solution or normal behavior is known, any type of behavior that deviate from the known pattern can be detected, regardless of the type or mechanism. This process is the approach to “detect” instead of “model” aberrant responses. Many statistical methods have been proposed to detect potential response disturbances in low-stakes assessments. For instance, Mahalanobis distance can detect aberrant response vectors whose patterns deviate from the centroid of data (Mahalanobis, 1936). Longstring analysis identifies individuals with the longest consecutive string of the same response, repeatedly endorsed (Huang et al., 2012; Costa Jr. & McCrae, 2008; Johnson, 2005). Previous articles provide a thorough review of best practices for using these and other methods (e.g., intra-individual response variability, psychometric synonyms, and psychometric antonyms) of aberrant response detection (Curran, 2016; Hong et al., 2020; Meade & Craig, 2012; Niessen et al., 2016). These methods typically do not involve latent variable modeling.

There are also many statistical methods to detect aberrant responses that are based on a latent variable model and residuals, which quantify how responses deviate from the main model for normal or solution behavior. Such residuals include person-fit statistics that are used to detect misfit between a subject’s test performance and their true latent trait. Typically, person-fit statistics compare an observed item response with an expected item response for a test item (Meijer & Sijtsma, 2001). Various aberrant behaviors can be identified, given that they deviate from the expected behavior of the subject according to the assumed model. Working under the IRT framework, we can leverage the item properties by considering the “appropriateness” of the subject response respective to the IRT model. A person-fit statistic that indicates that the response is atypical given the IRT model suggests that the response may be aberrant. A variety of person-fit statistics have been developed to detect unusual response vectors, such as appropriateness measures (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_3$$\end{document} and M), residual-based statistics (e.g., U and W), and standardized extended caution indices such as ECI1 and ECI2 (Drasgow et al., 1985, 1987, 1991; Levine & Rubin, 1979; Levine & Drasgow, 1988; Molenaar & Hoijtink, 1990; Tatsuoka, 1984; Trabin & Weiss, 1983; Wright & Stone, 1979; Wright & Masters, 1982). Karabatsos (2003) demonstrates the relative power of these and other person-fit statistics in detecting various aberrant responses across different test designs. A commonly used IRT-based person-fit statistic is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _z,$$\end{document} which compares the standardized log-likelihood of a response vector against the standard normal distribution (Drasgow et al., 1985). Previous studies have shown the promising ability of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _z$$\end{document} to detect C/IE and other types of aberrant responses (Hong et al., 2020; Meijer & Sijtsma, 2001; Niessen et al., 2016).

Responses that are flagged as aberrant by these statistical methods are often removed from the data set. Fully removing cases of aberrant data has been practiced for many decades. For instance, Cronbach (1950) argues that full removal is more appropriate than considering the aberrant responses as valid as other responses. Full removal continues to be a conventional albeit controversial method for handling response disturbances (Hong et al., 2020; Meijer & Sijtsma, 2001).

More recently, other methods stemming from the statistical quality control literature have also been developed to detect aberrant response patterns, such as change-point analyses (CPA) and cumulative sum (CUSUM) chart. In addition to indicating the presence of aberrant responses, these methods also try to determine when a certain response pattern starts deviating from the expected response pattern (Shao et al., 2016; Sinharay, 2016). For example, one CPA procedure detects if back random responding (BRR, which refers to random responding towards the end of an assessment, potentially due to fatigue or loss of interest) occurs and, if yes, when it begins (Yu & Cheng, 2019). Compared to the aforementioned methods, CPA enables partial removal of data, instead of complete removal. For example, when BRR is detected and a change point is identified, all responses after the change point would be excluded from the latent trait estimation.

Likewise, CUSUM procedures, proposed by Page (1954), identify fluctuations in the mean of the variable of interest. Instead of identifying a sharp change in responding during a test, CUSUM can detect a change in the quality of responses based on gradual shifts in the estimated latent trait. For instance, Meijer (2002) demonstrates, with an empirical certification test, how a positive CUSUM trend can indicate “warm-up” behavior (i.e., the respondent answers more items correctly as the test continues) while a negative CUSUM trend can indicate fatigue or guessing due to time constraints. Yu and Cheng (2022) outline a number of CUSUM statistics used to detect aberrant response patterns such as those developed by Armstrong and Shi (2009), Bradlow et al. (1998) and Van Krimpen-Stoop and Meijer (2000, 2001).

However, both CPA and CUSUM procedures still make a strong assumption about respondent behavior and are inflexible. For instance, a respondent could resort to random responding when they encounter a set of challenging questions but later returns to normal responding. Removing all responses beyond a certain point would not be appropriate in that case.

Robust estimation is a procedure where all data are leveraged to the extent possible while counteracting aberrant behavior. The type of robust estimation applied in this paper follows Huber’s M-estimation, a procedure in which data points further away from the centroid of the data receive a smaller weight (Huber, 1981). Accordingly, extreme cases have a smaller impact.

The idea of robust estimation has been applied to counteract aberrant responses. An advantage to using robust estimation is that it downweights data according to the degree to which each observation is aberrant. A more “suspicious” observation receives less weight, such that it will contribute less to the overall parameter estimation than a less “suspicious” observation. Thus, robust estimation accounts for partially aberrant responses, instead of removing data according to a binary standard. In the context of latent variable modeling, robust estimation has been employed in a number of applications. For example, in the context of IRT models, Hong and Cheng (2019b) proposed a robust marginal maximum likelihood estimator for estimating item parameters (or structural parameters) for Likert-scale items that follow the GRM.

This paper will focus on the maximum likelihood estimate (MLE) of the latent trait in the IRT framework. Some approaches to counteract bias in the MLE include Warm (1989) weighted maximum likelihood estimation (WLE). Although similar to the robust MLE in that the likelihood function is adjusted, WLE addresses bias inherent in the MLE (e.g., due to short test length) rather than bias due to aberrant responding. To address bias in the latent trait due to aberrant responding, one proposed solution is to weigh item responses linearly based on their position on a test. Items appearing later on an assessment are weighed less in ability estimation than the earlier items when performance decline is of concern, and vice versa when performance increase is of concern (Wise et al., 2023). However, this method solely pertains to strict performance decline/incline, and cannot be generalized to a wider range of response disturbances. The weights are also arbitrarily determined by item position, instead of being data-driven.

More general and systematic robust ML estimators for unidimensional IRT (UIRT) and multidimensional IRT (MIRT) models with dichotomous data have been developed for estimating the latent trait (Filonczuk et al., 2022; Mislevy & Bock, 1982; Schuster & Yuan, 2011). These robust estimators are designed to mitigate the effects of response disturbances, on the basis that outlying cases indicate more aberrant behavior. Item responses that poorly fit the IRT model are downweighted such that they contribute less to the estimation of the respondent’s ability. The estimators are able to target a variety of response disturbances and are not restricted to specific aberrant behavior. Both Huber and bisquare weight functions (Huber, 1981; Mosteller & Tukey, 1977) are utilized, which are described later in more detail.

However, these robust IRT developments have been limited to dichotomous data. No robust estimator has been developed for estimating latent traits from polytomously scored Likert-scale items. Given the wide application of Likert-scale items in psychological and educational assessments, we propose a robust estimator using the GRM.

Methods

According to the GRM, the probability that a subject i responds in or above a category k for item j is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P^*_{jk}(\theta _i) = P(X_{ij} \ge k| \theta _i) = \frac{1}{1+ e^{-Da_j (\theta _i-b_{jk})}}, \end{aligned}$$\end{document}

Embretson and Reise (2000), where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_j$$\end{document} is the item discrimination parameter applied to all category boundary functions for item j. D is a constant fixed to 1.7, commonly applied in order to scale the item parameters to that of a normal ogive model. Suppose there are K categories and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K-1$$\end{document} threshold parameters ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{j,1}; ... b_{j,K-1}$$\end{document} ), where the location parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_{jk}$$\end{document} separates response category k and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k+1$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=1,..K-1$$\end{document} ). The probability of endorsing exactly category k is therefore:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P_{jk}(\theta ) = P^*_{j,k}(\theta ) - P^*_{j,k+1}(\theta ), \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P^*_{j1}(\theta ) \equiv 1.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P^*_{jK}(\theta ) \equiv 0.0.$$\end{document}

The probability that a response falls in a category k for an individual with latent trait \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _i$$\end{document} is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P(X_{ij}=k|\theta _i) = P_{jk}(\theta _i) = \prod ^K_{k=1} P_{jk}(\theta _i)^{u_{ijk}}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{ijk}$$\end{document} is an indicator variable:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} u_{ijk} = {\left\{ \begin{array}{ll} 1 & \text {if } X_{ij} = k; \\ 0 & \text {otherwise.} \end{array}\right. } \end{aligned}$$\end{document}

We assume the number of categories remains constant for all J items on the assessment. In ML estimation, we estimate the latent trait by first finding the likelihood of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _i$$\end{document} given the response vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{x}_i = (x_{i1}, x_{i2}, ..., x_{iJ}),$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L(\theta _i|\varvec{x}_i) = \prod ^J_{j=1}\prod ^K_{k=1} P_{jk}(\theta _i)^{u_{ijk}}. \end{aligned}$$\end{document}

The log-likelihood,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {log}L(\theta _i|\varvec{x}_i) = \sum ^J_{j=1} \sum ^K_{k=1} u_{ijk}\text {log}P_{jk}(\theta _i), \end{aligned}$$\end{document}

is maximized by setting its first derivative equal to 0 and solving for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _i$$\end{document} .

In robust ML estimation, the likelihood of the item response function is weighted by a function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_{ij}$$\end{document} , resulting in the weighted likelihood \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$WL(\theta _i|\varvec{x}_i)$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} WL(\theta _i|\varvec{x}_i) = \prod ^J_{j=1}\left( \prod ^K_{k=1} P_{jk}(\theta _i)^{u_{ijk}}\right) ^{w_{ij}} \end{aligned}$$\end{document}

(Hu, 1994). The resulting weighted log-likelihood (Eq. 7) demonstrates how each item has a unique weighted contribution when determining the overall score:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {log}WL(\theta _i|\varvec{x}_i) = \sum ^J_{j=1} w_{ij} \sum ^K_{k=1} u_{ijk}\text {log}P_{jk}(\theta _i). \end{aligned}$$\end{document}

Assigning an item a weight of zero is the equivalent of removing the response from estimation. Setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_j$$\end{document} equal for all J items for a subject gives a weighted log-likelihood equivalent to that in ML estimation (Eq. 5) for finding the most likely \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} .

The Newton–Raphson method is then used to estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} iteratively with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \hat{\theta }_{i,(t+1)} = \hat{\theta }_{i,(t)}- \frac{\partial \text {log} WL(\hat{\theta }_{i,(t)})/ \partial \theta _i}{\partial ^2 \text {log} WL(\hat{\theta }_{i,(t)})/ \partial \theta _i^2}, \end{aligned}$$\end{document}

where t is the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^{th}$$\end{document} iteration. The first and second derivatives of the weighted log-likelihood with respect to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _i$$\end{document} are found as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial \text {log}WL}{\partial \theta _i}= & \sum ^J_{j=1} w_{ij}\sum ^K_{k=1} u_{ijk} a_j \nonumber \\ & \times \left( \frac{P^*_{j,k-1}(\theta _i)Q^*_{j,k-1}(\theta _i) - P^*_{j,k}(\theta _i)Q^*_{j,k}(\theta _i)}{P_{j,k}} \right) ,\nonumber \\ \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\partial ^2 \text {log} WL}{\partial \theta _i^2}= & \sum ^J_{j=1} w_{ij}\sum ^K_{k=1} u_{ijk} a_j^2\nonumber \\ & \times \left( \frac{P^*_{j,k-1}(\theta _i)Q^*_{j,k-1}(\theta _i)(Q^*_{j,k-1}(\theta _i)-P^*_{j,k-1}(\theta _i))}{P_{jk}(\theta _i)} \right. \nonumber \\ & \left. -\frac{P^*_{j,k}(\theta _i) Q^*_{j,k}(\theta _i)(Q^*_{j,k}(\theta _i)-P^*_{j,k}(\theta _i))}{P_{jk}(\theta _i)} \right. \\ & \left. -\frac{(P^*_{j,k-1}(\theta _i) Q^*_{j,k-1}(\theta _i) -P^*_{j,k}(\theta _i)Q^*_{j,k}(\theta _i))^2}{P_{j,k}^2(\theta _i)} \right) .\nonumber \end{aligned}$$\end{document}

As suggested by Schuster and Yuan (2011), the MLE is used as a starting value for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{i,(t=0)}$$\end{document} in the robustified Newton–Raphson algorithm. In our simulation and applied analyses, convergence was deemed reached when the log-likelihoods between consecutive iterations of the Newton-Raphson algorithm fell within a difference of 0.01. If a trait estimate converged outside the interval [-3.0, 3.0], the estimate was fixed to the endpoint closer in value (i.e., an estimate less than -3.0 was replaced with -3.0, and an estimate greater than 3.0 was replaced with 3.0), consistent with common practice in IRT research (Patton et al., 2014). Both of these practices were used for both ML and robust estimation in our simulations and applied example.

Standard errors of each estimate were obtained using the Huber-White sandwich estimator, a technique for estimating the variance of the MLE under a misspecified model (Huber, 1967; White, 1980). In the case of aberrant responses, the main IRT model is misspecified for them. Hence, the Huber–White SE is more appropriate than the usual Fisher-information-based SE for the MLE. The standard error (SE) of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_i$$\end{document} for subject i is estimated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$SE_{\hat{\theta }_i}=\sqrt{\hat{V}_i}$$\end{document} where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{V}_i=(-A_i)^{-1}B_i(-A_i)^{-1}.$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_i$$\end{document} is the second derivative of the weighted log-likelihood with respect to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} (Eq. 10) evaluated at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_i:$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} A_i=\frac{\partial ^2 \text {log} L}{\partial \theta ^2} \bigg | _{\hat{\theta }_i} . \end{aligned}$$\end{document}

Using Eq. 9, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_i$$\end{document} is estimated as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} B_i=\left[ \frac{\partial \text {log}WL}{\partial \theta } \bigg | _{\hat{\theta }_i} \right] ^2. \end{aligned}$$\end{document}

If the model is correct, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-A$$\end{document} should equal B,  and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{V}$$\end{document} simplifies to the inverse of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-A$$\end{document} . 95% confidence intervals were calculated with the following expressions

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left[ \hat{\theta }_i- 1.96*SE_{\hat{\theta }_i}, \hat{\theta }_i+ 1.96*SE_{\hat{\theta }_i} \right] . \end{aligned}$$\end{document}

Coverage rates were then computed as the proportion of subjects whose true \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} lays within the 95% confidence interval, over all replications.

Weighting mechanisms

It is essential to have a systematic and data-driven way to properly assign weights. Oftentimes the weight is a function of a residual \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{ij}$$\end{document} , which indicates the “deviation” of an observation from normal or expected behavior. In the IRT literature, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{ij}$$\end{document} is largely determined by the misfit between observed and expected (or model-based) responses, assuming that the IRT model fits well regular response data. The goal of the weighting function is to downweight contributions of items with large residuals. Previous residuals for robust estimation of ability with dichotomous data have drawn on item information to detect misfitting item responses when little information is to be expected from a subject on an item (Filonczuk et al., 2022; Mislevy & Bock, 1982; Schuster & Yuan, 2011). Previous residuals for responses on the GRM have been defined as a function of the difference between the observed response \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_{ij}$$\end{document} and expected score \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ E[X_{ij}|\hat{\theta }_i] = \sum _{k=1}^KkP_{j,k}(\hat{\theta }_i)$$\end{document} , such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{ij} = X_{ij} - E[X_{ij}|\hat{\theta }_i]$$\end{document} used in CUSUM procedures to evaluate person-fit (Van Krimpen-Stoop & Meijer, 2002). Although a variety of residuals could be used, we propose using the standardized residual

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} r_{ij} = \frac{X_{ij} - E[X_{ij}|\hat{\theta }_i]}{\sigma _{X_{ij}}}. \end{aligned}$$\end{document}

The residual is standardized by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{X_{ij}}$$\end{document} where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sigma _{X_{ij}}^2 = \sum _{k=1}^K (X_{ijk}-E[X_{ij}|\hat{\theta }_i])^2P_{jk}(\hat{\theta }_i). \end{aligned}$$\end{document}

The resulting standardized residual follow asymptotically the standard normal distribution, which allows for easy interpretation.Fig. 1. Weights as a function of the residual (black = bisquare weights, grey = Huber weights). Solid lines represent the tuning parameters used in the study (B=4.0 and H=1.0). Dashed lines represent weights when using higher tuning parameters (B=8.0 and H=2.0)

Given properly defined residuals that characterize how “suspicious” or misfitting a response is, a weight can be assigned to each response. In this study, two commonly used weight functions are considered: bisquare (Mosteller & Tukey, 1977),

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} w(r_j) = {\left\{ \begin{array}{ll} [1-(r_{ij}/B)^2]^2 & \text { if } |r_{ij}| \le B \\ 0 & \text { if } |r_{ij}| > B \end{array}\right. } \end{aligned}$$\end{document}

and Huber (Huber, 1981),

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} w(r_{ij}) = {\left\{ \begin{array}{ll} 1 & \text { if } |r_{ij}| \le H \\ H/|r_{ij}| & \text { if } |r_{ij}| > H, \end{array}\right. } \end{aligned}$$\end{document}

where B and H are tuning parameters that regulate the proportion of item responses and the extent to which they are downweighted. Prior studies have incorporated these weight functions when robustly estimating latent traits, so both are considered in our study (Mislevy & Bock, 1982; Schuster & Yuan, 2011). Figure 1 illustrates the weight functions across different residuals with different values for the tuning parameters, B and H. A greater absolute residual indicates a greater misfit. Therefore, we wish to downweight these cases more. As the residual increases in absolute value, less weight is given, while residuals close to 0.0 are given weights close to 1.0. Tuning parameters H and B regulate the proportion and degree of the item responses being downweighted. A higher tuning parameter leads to less downweighting. For example, when we use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B = 8.0$$\end{document} in comparison to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=4.0$$\end{document} (the dashed black line in comparison to the solid black line in Fig. 1), item responses with residuals greater in magnitude receive higher weights; thus, such observations are downweighted less. The same effect occurs when we increase the Huber tuning parameter from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=1.0$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=2.0$$\end{document} , shifting from the grey solid line to the grey dashed line.

To sum up, our robust estimation procedure operates as follows. First, we find the MLE for each respondent. Second, we calculate the residuals on all items for each respondent, identifying the degree of mismatch between the expected and observed item responses. Third, a weighting mechanism is chosen, and a weight is assigned to each item response as a function of the residual. Fourth, the robust estimate for each individual is updated incorporating these weights. Steps two through four are repeated, calculating the residual as a function of the latest updated latent trait estimate, until the Newton–Raphson algorithm converges.

Simulation study

While this estimation procedure is capable of counteracting a wide variety of response disturbances, our main simulation examines two examples of commonly observed aberrant responses in low-stakes assessments: back random responding (BRR) and overlooking reverse-coded statements. BRR occurs when respondents initially respond to items according to their true latent trait, then switch to a careless response pattern, which is continued for the remaining items. Previous studies have reported that careless responding primarily occurs on the latter items of an assessment, particularly when motivation is low for taking the assessment (Clark et al., 2003). Because of its prevalence, we chose to model BRR as a response disturbance in our simulations.

Furthermore, reverse-coded items are worded in a sentiment opposite from the scale (i.e., negatively worded), but are assumed to measure the same underlying construct (Hughes, 2009). An example item is “I feel unhappy” in a scale that measures happiness. A greater score on a reverse-coded item indicates a lower level of the latent trait measured by the item, and vice versa. They are often included to provoke the subject to pay more attention and process item content individually instead of forming response sets (i.e., responding according to general attitudes towards the test as a whole; Weems et al., 2003). In addition, adding reverse-coded items may help avoid biased responding that can result if all items are worded in a socially desirable direction (Weems et al., 2006). When including these items, it is assumed that subjects will respond to negatively worded items in a manner opposite to the positively worded items. However, subjects and researchers may both violate this assumption when encountering reverse-coded items. First, subjects may be inattentive and overlook reverse-coded items, responding erroneously (Hughes, 2009). Reading and interpreting reverse-coded items may require more mental effort and be more difficult for subjects, resulting in an aberrant response (Marsh, 1986; Williams & Swanson, 2001) as well. Second, researchers who perform secondary data analysis may be unaware of negatively worded items, due to carelessness on their end or unfamiliarity with the data set. Then, these items may not be properly reverse-coded prior to analysis. These aberrant behaviors violate the IRT model, since the response does not reflect the subject’s true latent trait when calculating the score. Thus, we use robust estimation in lieu of ML estimation to provide more accurate latent trait scores.

In our simulation study, we varied the following properties of a test: the number of categories on each Likert-type item, the test length, the severity, or proportion of items that are disturbed in each response vector, and the type of aberrant response: BRR or overlooked reverse coding.

Likert-type scales are commonly developed with either three, five, or seven categories, where each number of categories has its respective advantages (Aybek & Toraman, 2022). Therefore, we simulated tests with three, five, and seven categories, keeping the number of categories constant for all items within each test (e.g., a test could not have some items with three categories and some items with five categories).

Tests were generated with lengths of 10, 30, 50, and 70 items. Although multidimensional psychological inventories can include hundreds of items, such as the 567-item MMPI-2 (Sellbom & Anderson, 2013), we are concerned with unidimensional scales, so in this study the maximum test length was set at 70.

The robust estimation is most successful when the subject exhibits some behavior true to their latent trait and not all responses are aberrant, otherwise complete data removal would be sufficient. Previous studies have set the severity, or proportion of aberrant items, between 10%-50% (Meade & Craig, 2012). Therefore, we examined two severity conditions: 10% and 30%.

In effect, under back random responding, random guessing began after 90% or 70% (severities of 0.10 and 0.30, respectively) of the test. At this change point, each response category had an equal probability (e.g., 1/K) of being selected according to the discrete uniform distribution U[1, K], for an item with K categories.

In the neglected reverse coded simulations (i.e., subjects are inattentive to the negatively worded items), 10% or 30% of randomly selected items from each test were designated as negatively worded. Aberrant responses to these items were generated as follows. First, the item response \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_j$$\end{document} was generated based on the subject’s probabilities of endorsing each category according to their true ability. The response was then disturbed by being reversed in value, such that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} u^*_{ij} = {\left\{ \begin{array}{ll} u_{ij} & \text {if subject}\, i \, \text {is attentive to item } j \\ K+1- u_{ij} & \text {if subject}\, i \, \text {is inattentive to item } j \end{array}\right. } \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^*_{ij}$$\end{document} is the new response to replace \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{ij}$$\end{document} . For example, on a 5-pt Likert scale ranging from 1 to 5, consider a subject with a high latent trait who is likely to respond 1 to the negatively worded item. After reverse coding, the response will be coded as 5, reflecting the subject’s true latent trait. However, if the subject fails to recognize it’s a negatively worded item, their response would be 5. After reverse coding, the response would be 1, which indicates a low latent trait and not their true latent trait.

In total, we generated a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3\times 4\times 2\times 2$$\end{document} simulation design, incorporating three levels of item response categories, four test lengths, two severity levels, and two types of aberrant responses. For item parameter generation, we follow Dodd et al. (1989), who suggest generating \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j$$\end{document} s at random uniformly across the latent trait continuum, while varying the characteristics of the category threshold parameters within each item. In our study, we simulated \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_j \sim Unif(-2.5, 2.5)$$\end{document} to generate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K-1$$\end{document} values, one value for each threshold parameter, for each item j. Within each item, these values were sorted in ascending order such that the highest value corresponded with the greatest threshold. Item discrimination parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_i$$\end{document} s have been previously simulated according to Unif(0.90, 2.15) (Dodd et al., 1989). Therefore, our items used values from [0.90, 2.15] with equal increments between each value; e.g., for a 10-item assessment, the 10 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_i$$\end{document} s would be (0.900, 0.943, 0.986, ..., 2.106, 2.150). All item parameters were treated as fixed across all replications for each of the 24 conditions.

True \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} s from [-2.0, 2.0] in increments of 0.25 were used to generate data from each of the tests. Six hundred replications of each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} were used, forming 600 response vectors given the item parameters based on the following steps. Probabilities of responding in each category were generated for each of these true \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} s according to the item parameters and the item response function (Eq. 2). Polytomous data was generated through sampling, using these probabilities as the probability of landing in each category.

Three trait estimates were obtained: the MLE, the bisquare-weighted robust estimate, and the Huber-weighted robust estimate. The same response vector was used for all three estimators. Initial values of zero were used in the Newton–Raphson algorithm to estimate the MLE. Following past practices in robust IRT literature, the MLE was then used as the initial value to estimate the robust estimate, using zero in cases where the MLE did not converge. The average bias and MSE were evaluated for the three latent trait estimates. If the estimate did not converge within 30 iterations or approached infinity, the estimate was considered missing and removed from our calculation of bias and MSE. If more than 30% of the estimates for a trait level on a certain test design did not converge, we did not report the bias and MSE for the trait level on that test design.Table 1. Bias using maximum likelihood (ML), biweight (BI, B = 4), and Huber (HU, H = 1) estimation of the latent trait when there is back random responding in the data at a severity of 10%. (Nonconvergence rates higher than 30% displayed in parentheses) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} K**J−2−1.5−1−0.500.511.52ML310(79.5%)(43.3%)0.180.08−0.05−0.18−0.28(67.5%)(89.5%)300.420.370.280.170.06−0.01−0.09−0.15−0.21500.230.120.02−0.04−0.08−0.14−0.24−0.30−0.31700.310.230.150.07−0.01−0.08−0.16−0.19−0.25510(66.8%)0.15−0.03−0.17−0.27−0.30−0.24(35.5%)(80.8%)30(33.7%)0.180.120.03−0.02−0.07−0.10−0.17−0.2250(48.8%)0.160.100.050.00−0.05−0.11−0.17−0.31700.240.170.130.090.04−0.02−0.09−0.15(33.2%)710(78.7%)0.290.200.170.110.01−0.06−0.20(44.0%)30(89.8%)0.150.050.02−0.01−0.06−0.09−0.15(70.3%)50(80.3%)0.110.070.04−0.01−0.04−0.09−0.15(84.0%)70(91.0%)0.140.100.070.03−0.02−0.06−0.13(88.0%) HU310−0.08−0.050.130.10−0.03−0.12−0.080.040.0530−0.100.020.070.050.02−0.01−0.02−0.030.0750−0.08−0.04−0.08−0.020.010.02−0.07−0.110.0670−0.010.01−0.020.010.040.01−0.02−0.040.025100.01−0.04−0.05−0.10−0.13−0.12−0.040.090.2130−0.060.040.03−0.04−0.02−0.030.020.050.0750−0.040.030.020.010.030.00−0.03−0.010.02700.000.030.010.030.050.00−0.04−0.040.02710−0.150.020.090.100.090.02−0.04−0.08−0.0630−0.040.010.030.020.00−0.03−0.03−0.010.0250−0.020.010.030.01−0.01−0.03−0.02−0.030.0070−0.040.040.050.020.02−0.01−0.03−0.04−0.02BI310−0.33−0.120.120.09−0.02−0.12−0.080.090.1130−0.13−0.030.000.020.00−0.030.000.000.0950−0.13−0.09−0.10−0.020.030.05−0.04−0.060.0770−0.05−0.04−0.060.010.030.020.010.000.05510−0.09−0.05−0.04−0.08−0.10−0.09−0.020.07(35.0%)30−0.110.000.00−0.05−0.030.000.040.070.0850−0.07−0.010.000.000.020.01−0.010.020.0870−0.04−0.01−0.020.020.030.01−0.02−0.010.08710−0.29−0.060.060.080.070.01−0.03−0.05−0.0330−0.08−0.010.020.010.00−0.01−0.010.020.0550−0.05−0.010.010.01−0.02−0.010.00−0.010.0370−0.070.000.020.010.010.00−0.02−0.010.01Table 2Bias using maximum likelihood (ML), biweight (BI, B = 4), and Huber (HU, H = 1) estimation of the latent trait when there is back random responding in the data at a severity of 30%. (Nonconvergence rates higher than 30% displayed in parentheses) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} K**J−2−1.5−1−0.500.511.52ML310(58.2%)(30.2%)0.620.500.290.01−0.26(48.0%)(65.7%)300.950.710.440.220.01−0.19−0.36−0.55−0.71500.820.570.350.16−0.01−0.21−0.41−0.58−0.74700.710.500.300.12−0.08−0.29−0.49−0.66−0.78510(38.3%)0.350.07−0.13−0.25−0.37−0.47−0.70(33.8%)300.720.530.330.13−0.07−0.23−0.41−0.55−0.74500.630.450.260.10−0.05−0.18−0.34−0.52−0.75700.620.480.340.240.11−0.06−0.21−0.36−0.587100.940.700.560.430.310.150.00−0.30−0.63300.600.370.210.10−0.05−0.19−0.36−0.50−0.65500.580.440.310.170.02−0.12−0.26−0.43−0.60700.590.440.320.180.02−0.12−0.25−0.39−0.57 HU310−0.08−0.070.160.250.190.140.11−0.09−0.36300.010.090.180.13−0.01−0.10−0.12−0.14−0.07500.120.090.100.050.03−0.03−0.14−0.19−0.11700.060.070.040.070.01−0.09−0.15−0.15−0.095100.12−0.01−0.11−0.14−0.15−0.12−0.04−0.13−0.20300.200.250.180.05−0.08−0.16−0.15−0.05−0.03500.130.140.080.01−0.02−0.04−0.10−0.16−0.19700.090.110.090.090.090.00−0.07−0.11−0.187100.080.130.200.190.200.150.05−0.10−0.28300.140.120.070.05−0.01−0.08−0.15−0.17−0.13500.140.160.130.090.00−0.05−0.08−0.16−0.18700.150.160.140.06−0.01−0.06−0.08−0.14−0.17BI310−0.41−0.180.140.250.170.150.15−0.02−0.2630−0.14−0.120.040.08−0.03−0.10−0.05−0.020.0050−0.05−0.050.010.000.030.01−0.07−0.060.0170−0.08−0.07−0.050.040.02−0.06−0.04−0.010.035100.00−0.06−0.15−0.15−0.13−0.080.02−0.09−0.19300.040.110.110.03−0.07−0.11−0.070.060.07500.000.010.010.00−0.020.00−0.03−0.06−0.0470−0.04−0.020.020.050.060.02−0.02−0.01−0.02710−0.11−0.020.080.110.160.140.04−0.07−0.23300.030.030.020.020.00−0.02−0.06−0.06−0.01500.030.060.050.05−0.01−0.01−0.01−0.06−0.05700.020.040.060.03−0.02−0.02−0.03−0.05−0.05Table 3Sandwich standard error using maximum likelihood (ML), biweight (BI, B = 4), and Huber (HU, H = 1) estimation of the latent trait when there is back random responding in the data at a severity of 10%. (Nonconvergence rates higher than 30% displayed in parentheses) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} K**J−2−1.5−1−0.500.511.52ML310(79.5%)(43.3%)0.710.760.770.720.59(67.5%)(89.5%)300.620.470.320.210.210.280.340.340.36500.390.300.220.170.180.230.260.290.32700.280.230.180.140.140.190.250.280.29510(66.8%)0.470.420.400.420.460.46(35.5%)(80.8%)30(33.7%)0.310.230.170.170.210.280.380.5350(48.8%)0.240.170.130.150.180.200.250.34700.270.230.180.130.120.130.140.18(33.2%)710(78.7%)0.590.440.340.330.360.400.46(44.0%)30(89.8%)0.250.200.170.170.190.230.30(70.3%)50(80.3%)0.220.190.170.150.150.170.21(84.0%)70(91.0%)0.190.150.120.110.120.140.18(88.0%) HU3101.110.740.790.951.010.860.600.581.14300.960.570.390.250.220.300.350.360.50500.620.350.250.180.190.250.280.310.39700.350.260.200.150.150.210.270.290.345100.740.530.460.420.400.370.400.762.51300.640.350.250.180.170.210.300.450.69500.390.260.180.140.160.190.210.290.51700.310.240.180.130.130.130.150.200.327101.650.670.440.350.350.380.430.540.75300.460.280.210.180.170.190.240.350.53500.370.240.200.180.160.150.170.220.34700.330.200.160.130.120.130.150.200.29BI3102.661.260.840.951.000.870.610.651.29301.080.610.420.260.220.300.360.370.53500.750.370.260.180.200.250.290.310.39700.380.270.210.150.150.220.280.300.375101.190.610.470.420.400.370.440.76(35.0%)300.720.370.260.180.170.220.310.470.70500.420.260.190.140.160.190.210.310.57700.330.250.190.130.130.140.160.210.357102.490.890.470.350.350.380.440.560.76300.550.290.220.180.180.200.250.360.56500.390.240.200.180.160.150.170.230.35700.350.210.160.130.120.130.150.200.32Table 4Sandwich standard error using maximum likelihood (ML), biweight (BI, B = 4), and Huber (HU, H = 1) estimation of the latent trait when there is back random responding in the data at a severity of 30%. (Nonconvergence rates higher than 30% displayed in parentheses) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} K**J−2−1.5−1−0.500.511.52ML310(58.2%)(30.2%)0.730.720.690.640.57(48.0%)(65.7%)300.440.350.270.220.210.250.300.330.34500.260.210.180.170.190.220.250.270.29700.230.200.170.140.140.160.200.240.27510(38.3%)0.450.430.420.430.450.460.46(33.8%)300.300.250.200.170.170.200.240.290.36500.240.190.150.130.150.180.190.210.23700.230.190.150.130.120.130.140.160.197100.540.450.380.350.350.380.410.450.50300.260.220.190.170.170.180.200.240.31500.220.200.180.160.150.150.160.180.21700.190.170.140.120.110.120.130.150.18 HU3101.060.750.670.770.800.700.540.490.60300.910.560.360.240.220.270.330.350.47500.460.300.200.170.200.240.270.300.36700.340.240.190.150.150.190.250.280.325100.770.570.480.440.420.400.430.611.33300.460.300.230.180.170.190.250.380.63500.320.230.170.130.150.180.200.250.38700.290.230.170.120.130.140.150.200.297101.140.600.390.330.350.400.470.580.73300.400.270.210.170.170.190.220.290.45500.300.210.190.170.160.160.170.210.30700.270.190.150.130.120.120.140.180.25BI3102.741.310.780.760.810.720.560.560.79301.220.710.430.250.220.280.340.370.50500.620.340.230.180.200.250.280.310.39700.430.280.210.150.150.200.260.290.355101.280.650.510.460.430.410.470.651.27300.590.330.240.180.170.200.260.440.71500.390.250.180.140.150.180.210.280.47700.340.250.180.130.130.140.160.220.377101.720.850.450.340.350.400.480.610.80300.490.290.220.180.180.200.240.330.53500.350.230.200.180.160.160.180.230.35700.320.210.160.130.120.130.150.190.29

Tuning parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=4.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=1.0$$\end{document} were used for the bisquare and Huber weight functions, respectively, in all simulations. These values were chosen through a series of pilot simulations in which B and H were each varied, and the optimal result was the value that produced the most reduced bias and MSE across the test conditions in our main simulations. For brevity, we omit these results from the paper but can provide them upon request. Although we used a limited number of conditions to come up with these values, the parameter values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=4.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=1.0$$\end{document} are also consistent with the tuning parameters suggested in previous robust estimation of latent traits (Mislevy & Bock, 1982; Schuster & Yuan, 2011). However, for data with low degrees of aberrant responses (e.g., a severity of 10%) generated from tests with few categories (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=3$$\end{document} ), greater tuning parameters, such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=2.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=8.0$$\end{document} , tend to optimize the bias and MSE. Using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=1.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=4.0$$\end{document} in this situation results in a minor overcorrection of the data (i.e., for a MLE with a negative bias, the robust estimate may produce positive bias). Since increasing the tuning parameter decreases the amount of downweighting applied to the data, these results appear consistent: with fewer response disturbances, less downweighting should be applied to the data.

Results

Similar patterns were observed for both back random responding and neglected reverse coding, so our findings are generalized to both types of aberrant responding. Due to a large number of simulation conditions, tables from the simulations with neglected reverse coding were omitted from this paper but are available on OSF (https://osf.io/jku4d/).

When an estimation method led to 30% or greater nonconvergence for a certain true trait and test design, the bias and average standard error were not reported. Instead, the nonconvergence rate is reported in parentheses with the results (Tables 1, 2, 3 and 4). This incident occurred predominantly in ML estimation for extreme traits (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta =\pm 1.5$$\end{document} or even more extreme). However, the robust procedures both tend to result in much lower nonconvergence rates, suggesting that, when the ML estimation fails in convergence, robust estimation may be able to provide a converging estimate.

Overall, the bisquare- and Huber-weighted robust estimates resulted in less bias compared to the MLE. Tables 3 and 4 display the average bias for each of the three estimates (ML, Huber-weighted, and bisquare-weighted estimation) over each test design. For conciseness, a subset of true latent traits from the full simulation is shown ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta = $$\end{document} -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, and 2.0). Although both estimators were typically effective in reducing the bias, the bisquare weight function tended to produce the smaller bias of the two robust methods. For example, when the severity is 0.10, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=70$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=3$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta =1.5$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{ML}$$\end{document} gives a bias of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-0.19$$\end{document} , while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{H}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{B}$$\end{document} give \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-0.04$$\end{document} and 0.00, respectively. While the bisquare procedure does not always eliminate the bias, it typically is closer to zero than the Huber-weighted estimate. When K increases, the MLE tends to produce smaller bias, but not to the extent of elimination. Likewise, the Huber-and bisquare-weighted estimates provide smaller bias as K increases, to the point where the bisquare-weighted estimate almost reaches a negligible bias when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=7$$\end{document} .

Furthermore, when the severity is lower (Table 1), the bias is less in the MLE compared to instances of greater severity (Table 2). This pattern aligns with our expectations that, as more aberrant responses are present, the estimates will produce more bias. As the severity, number of categories, and test length all decrease, a slight overcorrection tends to appear in the bias by the robust estimators. For example, when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta =-2.00$$\end{document} , severity is 0.10, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=3$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=30$$\end{document} , the bias for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{MLE}$$\end{document} is positive, 0.40, while that of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{H}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{B}$$\end{document} is -0.13 and -0.16, respectively. Beyond \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=10$$\end{document} (with greater K and severity) the test length does not have much of an impact on the bias.

Therefore, the robust estimation works best with assessments with more than ten items and/or more than three categories. For such assessments with fewer items and categories, we recommend using robust estimation with caution, though the overcorrection that occurs in shorter tests can be mitigated with a larger tuning parameter.

Average sandwich standard errors for the latent trait estimates can be found in Tables 3 and 4. Standard errors are slightly larger for the robust estimates compared to MLE in some cases of extreme latent traits (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm 2.0$$\end{document} ) or cases of short test lengths (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=10$$\end{document} ), but is otherwise negligible. Previous studies have reported increases in standard error consistent with these results when using robust estimation in general (Carroll & Pederson, 1993) and more precisely in item response models (Schuster & Yuan, 2011). When evaluating estimation procedures, the bias-variance trade-off must be considered: reducing the bias risks an increase in variance. However, the standard error does not increase substantially with the robust estimation in our study, compared to that of the MLE. Therefore, we achieve a desirable bias-variance balance, and the total error is not inflated with robust estimation.

Our findings are consistent with the robust literature. Overall, the robust estimates resulted in less bias along with comparable standard errors, compared to the MLEs. We find that the bisquare-weighted estimate performs better than the Huber-weighted estimate in most cases. This finding contradicts the findings of Schuster and Yuan (2011), which suggested the Huber-weighted estimate outperforms the bisquare-weighted estimate due to a high incidence of infinite estimates with the bisquare weight. Because our study pertains to Likert-type data, as opposed to dichotomous data studied by Schuster and Yuan (2011), less nonconvergence occurs since it is less likely that all responses fall in the lowest/highest response category (similar to all zeros or all 1’s in dichotomous data). Therefore, the issue of infinite estimates does not impact our bisquare estimation here in polytomous items as much as in dichotomous items, and the bisquare weight indeed tends to outperform the Huber weight in the robust estimation.

The Huber weight outperforms the bisquare weight in very few instances, such as when there is both low severity and smaller K. Here, the bisquare weight function results in a greater degree of overcorrection than the Huber weight function. Moreover, the Huber weight tends to result in a smaller standard error in short test lengths with more extreme traits. In all other cases, the bisquare weight tends to result in less bias and standard errors comparable to the Huber weight function.Table 5. Coverage rates for the sandwich standard error using maximum likelihood (ML), biweight (BI, B = 4), and Huber (HU, H = 1) estimation of the latent trait when there is back random responding in the data at a severity of 10%. (Nonconvergence rates higher than 30% displayed in parentheses) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} K**J−2−1.5−1−0.500.511.52ML310(79.5%)(43.3%)0.970.950.910.920.98(67.5%)(89.5%)300.950.890.840.840.990.940.950.980.98500.970.960.970.980.910.890.880.880.92700.860.870.870.940.960.920.900.950.96510(66.8%)0.971.000.990.920.920.96(35.5%)(80.8%)30(33.7%)0.950.920.950.920.940.960.960.9850(48.8%)0.950.900.940.950.970.970.970.95700.980.970.900.870.910.930.900.92(33.2%)710(78.7%)0.980.950.940.980.970.970.96(44.0%)30(89.8%)0.960.930.920.920.930.940.97(70.3%)50(80.3%)0.990.980.960.980.970.950.95(84.0%)70(91.0%)0.960.910.880.930.940.950.96(88.0%) HU3101.001.001.000.990.980.971.001.000.98301.000.990.960.910.970.950.971.001.00501.001.000.990.960.960.980.990.981.00701.000.990.970.940.960.991.001.001.005100.980.990.990.980.930.940.940.980.99301.000.990.980.970.920.960.980.991.00501.001.000.980.950.970.980.990.991.00701.001.000.990.940.930.950.960.971.007101.000.980.960.950.960.970.980.980.99301.000.990.960.960.970.970.980.991.00501.001.000.990.980.980.970.980.991.00701.000.990.980.960.950.970.980.991.00 BI3101.001.001.001.000.990.961.001.000.98301.000.990.970.940.970.940.961.001.00501.001.000.990.950.960.981.000.991.00701.001.000.980.940.960.981.001.001.005100.980.980.990.970.920.940.950.99(35.0%)301.000.990.980.960.900.960.980.991.00501.001.000.970.960.970.981.000.991.00701.001.000.980.950.920.940.970.991.007101.000.990.960.960.960.960.980.970.99301.000.990.980.950.960.960.990.991.00501.001.000.990.990.970.960.980.991.00701.001.000.980.960.950.970.990.991.00Table 6Coverage rates for the sandwich standard error using maximum likelihood (ML), biweight (BI, B = 4), and Huber (HU, H = 1) estimation of the latent trait when there is back random responding in the data at a severity of 30%. (Nonconvergence rates higher than 30% displayed in parentheses) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} K**J−2−1.5−1−0.500.511.52ML310(58.2%)(30.2%)0.840.760.770.850.97(48.0%)(65.7%)300.430.420.500.710.940.750.680.630.49500.160.280.490.780.880.770.620.450.29700.170.310.540.770.850.530.340.290.20510(38.3%)0.790.950.910.800.760.800.75(33.8%)300.370.440.540.780.840.680.580.490.43500.260.390.550.800.800.750.580.300.09700.220.340.410.500.680.820.600.390.177100.630.600.570.690.800.910.950.870.72300.360.590.710.720.780.690.540.450.44500.240.430.550.730.850.780.590.370.19700.130.300.390.600.790.720.530.320.13 HU3100.980.980.950.900.840.760.981.000.87300.980.850.790.800.920.820.870.950.94500.950.880.820.830.880.910.900.880.93700.950.920.870.800.870.870.850.940.955100.940.960.960.910.830.870.930.940.94300.890.790.780.870.830.730.800.940.98500.950.900.880.860.800.870.880.900.90700.980.940.880.800.770.900.870.910.917100.950.910.820.860.820.890.970.950.89300.950.890.870.820.840.850.820.890.96500.930.870.850.870.920.860.880.860.92700.960.860.790.850.830.810.860.880.92 BI3100.990.980.940.880.840.781.001.000.92300.990.930.830.850.910.790.850.970.97500.980.940.870.850.850.900.930.930.98700.990.980.890.790.870.870.930.980.995100.950.950.940.890.830.880.960.940.94300.950.880.800.880.800.770.860.970.99500.980.960.910.850.780.850.930.950.97701.000.980.930.830.780.880.910.950.977100.980.950.880.880.830.870.950.940.89300.990.930.880.810.800.850.910.950.98500.980.950.910.880.890.840.910.930.97700.990.960.870.860.800.820.920.940.98

Fig. 2. Comparing the MLE of each latent trait with its robust estimate. The figure on the left is the bisquare-weighted robust estimate and the figure on the right is the Huber-weighted robust estimate. Points in black represent the 20% of subjects with aberrant data, where 33% of their responses are oppositely coded. Points in gray represent subjects with “nonaberrant” data. Subjects were not plotted if their MLE or robust estimate did not converge

Tables 5 and 6 display the coverage rates for BRR at severities of 0.1 and 0.3, respectively. Across all conditions, the coverage rates for the robust estimates are greater than those of the MLE. While the coverage rates for the MLE decrease with increasing severity, those for the robust estimate remain relatively similar. Thus, for higher severities, the advantage of the robust estimation on the coverage rates is more prominent, as it is able to mitigate the effects of increased aberrant behavior.

Applied analysis

We applied our robust estimation procedure to data from the Big Five Inventory-2 (BFI-2) administered to adolescents ages 14 to 17 (Ober et al., 2021). The sample size of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=838$$\end{document} contained no missing responses. The scale consists of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=60$$\end{document} items on a five-point Likert scale, with 12 items loading onto each of the five personality factors. However, because we are concerned with the unidimensional GRM, we focused on one personality dimension, neuroticism, which is measured by 12 items. Six of these items (50%) were negatively worded such that they should be reverse coded before analysis. Item parameters were estimated using the properly reverse coded data for all 838 subjects.

160 respondents were randomly chosen to be aberrant (20% prevalence). For these subjects, aberrant responses were created such that four negatively worded items had responses opposite to the true intended responses (33% severity). In doing so, we tried to mimic a scenario where subjects were being inattentive and overlooked the negative wording, for example, missing a “NOT” in the question. Discrepancies between the MLE ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{MLE}$$\end{document} ) and robust estimates ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{B}$$\end{document} when the bisquare weight is used and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{H}$$\end{document} when the Huber weight is used) were compared between the aberrant and the nonaberrant, or undisturbed, responders. It is important to note that there may be other response disturbances that occur in the data (e.g., other types of careless responding), because they frequently plague real data. Therefore, our “nonaberrant” group data may not be entirely undisturbed. However, it is still presumably much less undisturbed than our artificially aberrant data. Likewise, the aberrant group may contain additional response disturbances beyond those artificially created.

Focusing first on the nonaberrant group in Fig. 2 (gray data points), the robust estimates are approximately the same as the MLEs for both the Huber- and bisquare-weighted cases, as most gray points fall close to the identity line, especially for the Huber-weighted cases. However, when we consider the aberrant respondents (black points), a different pattern emerges. For aberrant respondents with small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{ML}$$\end{document} , the corresponding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_H$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_B$$\end{document} tend to fall below the identity line; whereas aberrant respondents whose \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{ML}$$\end{document} are large tend to have \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_H$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_B$$\end{document} that fall above the identity line. In other words, the robust estimate for a subject tends to be pulled to the extremes (e.g., -2.0 or 2.0), while their MLE tends to be approximately average (i.e., near 0.0). We expect to see this pattern when a respondent misses some negatively worded items, as the mix of some high responses with some low responses will bias the MLE towards 0.Table 7ML, Huber-weighted, and bisquare-weighted estimates ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{MLE}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_H$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_B$$\end{document} , respectively) of neuroticism for three subjects alongside their responses. Note: Aberrant (mistreated reverse coded) items are marked with “R”. Values marked with “ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${-}$$\end{document} ” did not convergeID \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{MLE}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_H$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_B$$\end{document} 1R2R345R6R789101112212−0.18−2.37−2.885512552114213431.051.632.87244422552255650−2.853.00144521554555

Table 7 provides the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }$$\end{document} s for three example aberrant subjects based on their responses for the 12 items labeled 1 through 12, with “R” indicating one of the four negatively worded items that a respondent might miss. Consider, for example, the response vector for Subject 212. Their responses to regular items mostly consist of 1s and 2s, indicating low neuroticism. Neglecting negatively worded items, their responses to these items remained low, resulting in 5s on these items after reverse coding, which misled the MLE to produce higher estimates. The MLE is biased towards zero, as the aberrant items are increasing the score: a mix of high and low responses leads to an average trait estimate. However, when robust estimation is employed, the trait estimate is much lower ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{H} = -2.37$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{B} = -2.88$$\end{document} ), indicative of the subject’s actual low neuroticism. This estimate more closely reflects the subject’s true latent trait as suggested by the undisturbed responses. The discrepancy between the robust estimates and the MLE is much greater for Subject 212 in Fig. 2a and b than the non-aberrant subjects in this sample, as it falls further away from the identity line.

Likewise, the response vector for Subject 343 in Table 7 consists of mostly 4s and 5s, indicating a high level of neuroticism. The aberrant responses on the four negatively worded items then wrongfully indicate a lower level of neuroticism, and reduce the MLE. However, the robust procedure properly recovers the trait estimate to the high level of neuroticism ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{H} = 1.63$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{B} = 2.87$$\end{document} ). Subject 343 falls above the identity line in Fig. 2a and b, demonstrating how their robust estimate increased from the MLE.

Moreover, the trait estimates for Subject 650 demonstrate a common phenomenon in this robust estimation. Oftentimes, the MLE for a subject will not converge while the robust estimates will converge. The robust estimates ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{H} = 2.85$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{B} = 3.00$$\end{document} ) then reflect the high ability as suggested by the undisturbed responses. Thus, robust estimation can be used to overcome nonconverging MLEs. Because the MLE did not converge, this subject was neither plotted in Fig. 2a nor Fig. 2b. Note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_{B}$$\end{document} did converge to a value greater than 3.00, so the estimate was truncated to 3.00.

Discussion

We have proposed a general robust estimation procedure for estimating latent traits with the GRM. We have found that the robust estimation proposed in this paper succeeds in reducing bias for data with various types of disturbances, without substantially increasing the standard error. The estimation is effective at various test lengths and number of response categories.

Both the Huber and bisquare weight functions are effective in reducing the bias in most cases, while providing the same or comparable standard error with respect to the MLE. The choice of either weighting mechanism is left to the researcher. Schuster and Yuan (2011) suggest that the Huber weight should be used over the bisquare weight in cases where nonconvergence may be an issue. Potential nonconverging cases may be identified by examining the data for response vectors which contain all, or predominantly all, responses in the same response category. However, given no concern about nonconvergence, the bisquare weight should be used if reduced bias is preferred. Indeed, in our simulations, we find that bias tends to be slightly reduced for the bisquare weight over the Huber weight, except in some conditions.

Although most differences between the results of the Huber and bisquare weight mechanisms are trivial, the decision of using one over the other may also be influenced by some observable characteristics of the data (e.g., short test lengths and/or few response categories). For data contaminated with a small proportion of aberrant responses (i.e., low severity), in cases of smaller K, the bisquare weight function may lead to a greater degree of overcorrection. In this case, the Huber weight function may be preferred. This discrepancy is negligible for greater K, where the bisquare weight function tends to reduce the bias more than the Huber weight function. Typically, the bisquare and Huber weights tend to produce very similar average standard errors. However, the Huber weight tends to give a smaller standard error than the bisquare weight for short test lengths ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J=10$$\end{document} ), especially for more extreme traits (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta = \pm 2.0, \theta = \pm 1.5$$\end{document} , etc.).

Although we recommend using the tuning parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=1.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=4.0$$\end{document} for the Huber and bisquare weight functions, respectively, it is ultimately up to the researcher to choose a tuning parameter. For a greater degree of downweighting, the tuning parameter can be decreased, while a lesser degree of downweighting can be accomplished with a greater tuning parameter. When the number of categories is low (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=3$$\end{document} ), the robust estimation is more sensitive to the choice of tuning parameters. For high degrees of aberrant behavior (e.g., 30% severity), using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=1.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=4.0$$\end{document} operates well; however, for lesser degrees of aberrant behavior (e.g., 10% severity), a greater tuning parameter such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=2.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=8.0$$\end{document} may be more optimal in order to avoid overcorrection when items have few categories. However, in our preliminary analysis of H and B, we found that, generally, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H=1.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=4.0$$\end{document} are the most robust across different test conditions and response behaviors. These values are consistent with the recommendations in the literature as well. Therefore, they are suggested as a rule of thumb.

While the study has demonstrated the advantage of the robust estimator for polytomous items, it can be extended in multiple ways. First, the study is limited to the specific type of residual and weight functions we chose. Different residuals and weighting mechanisms with the same objectives of detecting misfit and downweighting large misfit can be substituted. Future studies may compare this robust estimation with other weighting mechanisms that target aberrant responding, such as those proposed by Wise et al. (2023). In addition, it may be of interest to compare the robust estimation in this analysis to the WLE instead of the MLE, particularly when the test is short, such as ten items (Warm, 1989). Although bias due to aberrant responding is typically more severe, the MLE may contain bias due to short test lengths, and the WLE may be a stronger benchmark in such conditions.

Second, this study is limited by the assumption that the item parameters are known. This assumption is practical given that many IRT-based testing programs use item parameters calibrated from pretesting. However, the true item parameters are never exactly known and are not exempt from misestimations due to aberrant data in the calibration sample (Patton et al., 2019; Tsutakawa & Johnson, 1990). Future studies should investigate these limitations with regards to robust estimation.

Lastly, the analyses in this study are limited to unidimensional data. Future investigations with this estimator could extend to multidimensional polytomous models, such as the multidimensional graded response model (Muraki & Carlson, 1995). It is also important to consider mixed-format assessments, where dichotomous and polytomous items both exist (Hong et al., 2021). A general solution that subsumes multiple item types and underlying latent dimensions will be very desirable. Moreover, robust estimation has been applied to improve estimates of subjects’ working speed under the response time model (Hong et al., 2021), and both item responses and response times information have been used in person fit analysis (Gorney et al., 2024) and detection of rapid guessing behavior (Lu et al., 2020). It is therefore potentially plausible and useful to develop robust estimators for joint models of item responses and response times.