Nonignorable Consequences of (Partially) Ignoring Missing Item Responses: Students Omit (Constructed Response) Items Due to a Lack of Knowledge

Abstract

In recent literature, alternative models for handling missing item responses in large-scale assessments have been proposed. Based on simulations and arguments based on psychometric test theory, it is argued in this literature that missing item responses should never be scored as incorrect in scaling models but rather treated as ignorable or handled based on a model. The present article shows that these arguments have limited validity and illustrates the consequences in a country comparison using the PIRLS 2011 study. It is argued that students omit (constructed response) items because they do not know the correct item answer. A different treatment of missing item responses than scoring them as incorrect leads to significant changes in country rankings, which induces nonignorable consequences regarding the validity of the results. Additionally, two alternative item response models are proposed based on different assumptions for missing item responses. In the first pseudo-likelihood approach, missing item responses for a particular student are replaced by a score that ranges between zero and a model-implied probability computed based on the non-missing items. In the second approach, the probability of a missing item response is predicted by a latent response propensity variable and the item response itself. The models were applied to the PIRLS 2011 study, demonstrating that country comparisons change under different modeling assumptions for missing item responses.

1. Introduction

In educational large-scale assessment studies [1], such as the progress in international reading literacy study (PIRLS; [2]), the trends in international mathematics and science study (TIMSS; [3]), or the programme for international student assessment (PISA; [4]), students’ competencies are assessed using cognitive test items. However, students often do not respond to specific items, leading to missing item responses [5]. It is not obvious how item non-response [6,7] should be treated in the computation of values of competencies (i.e., values of the latent trait) in item response theory (IRT) models [8,9,10,11,12] that are used as scaling models.

While the treatment of missing data in statistical analyses in social sciences is now widely used [13,14], in recent literature, there has been criticism of conventional methods of treating missing item responses in item response models in large-scale assessment studies [5,15,16,17]. Typically, the treatment of item responses can be distinguished between the process of calibration (determination of item parameters) and scoring (determination of ability distributions) [18].

In PIRLS, TIMSS, and PISA (until PISA 2012), only missing item responses at the end of a test booklet (so-called not-reached items) are omitted (i.e., ignored) in calibration in order not to distort estimates of item parameters [2,18,19]. In the scoring step, missing item responses are evaluated as incorrect to compute values of the latent trait as plausible values [20,21]. Since PISA 2015, not-reached items have been ignored in calibration and scoring (although the proportion of missing items is used as a predictor in the latent regression model; see [17]).

However, some researchers [17,22,23,24,25,26] argued that missing item responses should never be treated as incorrect (i.e., in calibration, as well as in scoring). These researchers proposed alternative item response models to handle missing item responses properly. This article will first elaborate on the criticisms made in the literature [16,17,25] to treat missing item responses as incorrect and the proposed model-based alternatives. In contrast to the above-mentioned literature, we argue that only for reasons of validity and not for test-theoretical reasons should a particular method for handling (i.e., both in calibration and in scaling) be preferred. Since assumptions about the missingness of item responses are not empirically testable, we propose two alternative item response models that directly parametrize different assumptions of the missingness mechanism. Finally, we discuss potential psychometric consequences in large-scale assessment for the assessment of competencies with varying treatments for missing item responses.

2. Analysis of the Critique of Traditional Approaches to Handling Missing Item Responses

Typically, missing item responses are scored as incorrect when determining ability distributions in large-scale assessment studies [24]. In a series of recent publications, it was claimed that missing item responses should never be scored as incorrect. It is recommended that this traditional ad hoc method should never be used in large-scale assessment studies [5,17,23,24,25]. In the following, we discuss the arguments of these authors and conclude that they are based on assumptions and conclusions that are atypical in the application of item response models. Therefore, we find the recommendations in the above-mentioned literature somewhat doubtful.

2.1. Aleatoric and Epistemic Uncertainty

Useful for the following considerations is the distinction of Denoeux [27] in the quantification of uncertainty in data (the item responses) and corresponding IRT models. Denoeux distinguishes the concepts of probability and possibility. Uncertainty due to sampling or the modeling of existing data (a dataset containing all item responses) is denoted as aleatoric uncertainty [27]. Uncertainty (or fuzziness) for every single datum (i.e., each individual item response) is denoted as epistemic uncertainty [27]. Epistemic uncertainty (or systematic uncertainty) can be attributed to a lack of knowledge with respect to determining a single datum [28]. In contrast, aleatoric uncertainty refers to the notion of randomness. It describes the variability in the outcome of an experiment (which can be an educational large-scale assessment study) due to inherently random (i.e., sampling) effects. In [27], it is argued that aleatoric uncertainty can be best expressed by probabilities in a statistical model. In the case of educational large-scale assessment studies, this means that an IRT model for the dataset of item responses with fixed values is fitted. Each student’s individual item responses are viewed under the concept of epistemic uncertainty, which would result in fuzzy item responses [27,29]. Therefore, it is unclear how one should score missing item responses.

Nevertheless, in this interpretation, we would like to emphasize that it is vital that assumptions about the individual data points (i.e., the individual item responses) be strictly distinguished from the assumptions in a probability model, which is a statistical model for item responses of all students. To illustrate this conceptual distinction, we assume a distribution of the running times in a 50 m run for a defined population of Austrian students in the fourth grade. For this purpose, a sample of students is drawn. The associated uncertainty of sampling students from the population is associated with probabilities in the context of aleatoric uncertainty. However, the measurement of the running time of a student can never be precisely carried out (e.g., due to measurement uncertainty, the runtime is located in the interval between 8.75 and 8.85 s, which represents a fuzzy date). This kind of uncertainty corresponds to the concept of epistemic uncertainty.

In summary, we reiterate that assumptions about the emergence of data (aleatoric uncertainty) must be distinguished from assumptions about the modeling of existing data (epistemic uncertainty). It does not seem wise to mix up these two processes. Unfortunately, this is precisely the case in the reasoning of why missing item responses should not be scored as incorrect [15,23,24,25].

2.2. Reasoning Based on Foundations of Psychometric Test Theory

The criticism of treating missing item responses as incorrect is based on simulations and test-theoretical arguments [16,25], which we will discuss in detail below.

Most simulation studies convey that scoring of missing item responses as incorrect provides distorted item parameters and ability distributions [30]. In these simulation studies, a missing item response is mostly simulated due to a dependence on latent variables or person covariates but not depending on the item itself. Then, that the model would produce biased estimates when treating missing item responses as incorrect obviously (see [30,31]) follows because the simulation is based on a data-generating model that is not in correspondence with the analysis model (i.e., treating missing item responses as incorrect). In contrast, if a particular IRT model fits, and missing values are generated for all items for which the item was solved incorrectly, then scoring missing item responses not as incorrect will consequently provide distorted model parameters (see [30]). Whether simulation-based reasoning against the scoring of missing item responses as incorrect is valid depends on the plausibility of the data-generating model for the missingness process [31]. We argue that approaches that model missingness on an item independently of the unknown item response itself are implausible. Thus, simulation studies are not helpful in terms of justification for choosing an adequate treatment of missing item responses in real data [31,32].

In [33], simulation studies are classified into method-centred and neutral simulation studies. Method-centred simulation studies are included in research papers that propose a new statistical model. These kinds of simulation studies frequently simulate data that fit the proposed model. The criticized papers on missing item responses [15,17,24] are particular papers that include method-centred simulation studies. Neutral simulation studies are studies aiming at comparing different methods. They should not intentionally prefer one particular statistical approach. The issue with some method-centred simulation studies like [17,34] in the literature on missing item responses is not that they propose new models. However, they make recommendations on how to generally handle missing item responses based on generated datasets that were perfectly suited to their item response models implemented in their favorite toolbox. Therefore, we would like to emphasize that simulation conditions and studied models must be as wide as possible if there is a desire to derive recommendations for practitioners. Otherwise, it sounds more like an advertisement than serious research practice if researchers only focus on their selective toolbox.

In a reasoning based on psychometric test theory, the core of the criticism is that scoring missing item responses as incorrect would model missing item responses deterministically [16,25]. This modeling strategy would contradict model-implied item response probabilities in an IRT model, and it is therefore deemed invalid [17,25]. Let us denote the missingness indicator variable with $R_{pi}$, that for a person p an item i is observed (i.e., $R_{pi}=1$). Then, formally, for the incorrect scoring, a missing item response $Y_{pi}$ is set to zero; thus, $\mathsf{P}(Y_{pi}=1|R_{pi}=0)=0$, independent of the latent trait ${\mathsf{\theta }}_p$ of person p. This would lead to a contradiction because the prediction in the IRT model (e.g., in a two-parameter logistic (2PL; [35]) model) would result in a probability of $\mathsf{P}(Y_{pi}=1|{\mathsf{\theta }}_p)>0$, i.e., the data handling is deterministic, but the IRT model is probabilistic [16,25]. Hence, model assumptions in the IRT model would be violated. The contradiction can be tackled by using the different concepts of uncertainty according to [27] because in the reasoning of [16] or [17], epistemic and aleatory uncertainty are confounded. In this literature, the same probability model (i.e., the item response probability from an IRT model) for the occurrence of a single data point and modeling all item responses in the sample is employed.

If we embed our considerations now more firmly in concepts of psychometrics, it is noticeable that the rejection of treating missing item responses as incorrect relies on an intraindividual interpretation of the probabilities in the IRT model. This reasoning corresponds to a stochastic subject perspective (i.e., probabilities of the IRT model can be interpreted for each combination of a person p and an item i; see [36]), which is not typically employed in large-scale assessment studies [37]. If one uses a random sampling perspective [36,38], the probabilities are interpreted as the result of a sampling of persons, and item parameters, as well as a latent trait distribution for $\mathsf{\theta }$, are introduced as a representation of a high-dimensional contingency table for discrete multivariate observations $\mathit{Y}$. Then, formally, the IRT model for I discrete items $Y_1,\dots ,Y_I$ can be written as

$$\mathsf{P}(\mathit{Y}=(y_1,\dots ,y_I))=\int \prod _{i=1}^I\mathsf{P}(Y_i=y_i|\mathsf{\theta })f\left(\mathsf{\theta }\right)\mathrm{d}\mathsf{\theta },$$

(1)

where f denotes the density function of the latent variable (i.e., latent trait or (latent) ability) $\mathsf{\theta }$. In this notation, it becomes clear that in IRT models, a distribution is modeled for a sample of persons, and individuals are not represented (see also [30]). However, if one assumes the intraindividual perspective in the sense of the stochastic subject perspective, the reasoning is based on the assertion that the scoring of missing item responses as incorrect would lead to a deterministic item response probability of zero because the scoring as incorrect leads to a wrong item response with a probability of one (i.e., $\mathsf{P}(Y_{pi}=0|R_{pi}=0,{\mathsf{\theta }}_p)=1$). However, the scoring of an item response as incorrect is not related to the probabilistic modeling of item responses of a person. This is because the missing item responses could also have been generated by a deterministic process, which does not contradict the probabilistic modeling in the statistical model for all item responses of all persons. Therefore, one has to distinguish between a real response process—which can be either deterministic or probabilistic but certainly follows the given probabilistic specification of the IRT model in only rare cases—and the statistical model assumptions for the definition of ability as defining elements for the construction of an ability scale $\mathsf{\theta }$ (see [30]).

Furthermore, we argue, in the following, that statistical reasons from the stochastic subject perspective invalidate the criticism of scoring missing item responses as incorrect. When modeling intraindividual distributions in an IRT model for a person, item parameters are defined as fixed and known parameters, and the latent trait ${\mathsf{\theta }}_p$ is required as a fixed effect that is assigned to a person. Then, the assumption of local stochastic independence refers to conditional independence of item responses for a fixed person p based on a set of items, and it is therefore not empirically falsifiable but a means for identifying the ability ${\mathsf{\theta }}_p$ of this person. Hence, the ability ${\mathsf{\theta }}_p$ is only defined by the specification of likelihood and item responses $\mathit{Y}$. In the Rasch model [39], for a fixed person p and item i, the probability of a correct item response is expressed by $\mathsf{P}(Y_{pi}=1|{\mathsf{\theta }}_p)=\mathsf{\Psi }({\mathsf{\theta }}_p-b_i)$ (where $\mathsf{\Psi }$ denotes the logistic function), which is used in the likelihood estimation for the distribution of the trait ${\mathsf{\theta }}_p$. Whether the value zero or one for a single observation ($x_{pi}$ fits) in the Rasch model is empirically undecidable, so scoring missing item responses as incorrect does not necessarily violate the IRT model [30].

It becomes clear that we prefer the random sampling perspective because the item response probabilities $P_{pi}$ from an IRT model do not describe how likely a particular student p is to correctly solve item i. In contrast, we believe that the probability is likely deterministic at the item level. Randomness in IRT models (and, hence, probabilities between zero and one) only emerge because we summarize interindividual differences of multivariate item response data. For example, a particular written arithmetic item can be solved by a student if she or he is capable of applying the necessary rules required for written arithmetic. Hence, whether students correctly solve this item relates to opportunities to learn and whether items are instructionally sensitive [40,41]. As a consequence, specific item effects are regarded as construct-relevant and not as unsystematic errors [42,43].

We would like to emphasize that our arguments only apply to missing item responses for items that were administered to students in a test booklet. In large-scale assessment studies such as PIRLS, planned missingness designs (i.e., balanced incomplete block designs) are utilized in which only a subset of items is administered to students in a test booklet [44]. Such missing item responses should—of course—not be scored as incorrect because students never saw these items. Hence, lack of knowledge cannot be the reason for planned missing data. However, in the case of administered items, we cannot rule out the possibility that students do not provide an item response because they do not know the solution to this item.

To sum up, in our opinion, the reasonings given in the literature [16,17,25] for not scoring missing item responses as incorrect that are based on simulation studies or psychometric test-theory-based reasoning do not seem to be typical for applications in large-scale assessment studies. However, we note that for reasons of validity, quite different scoring methods than incorrect ones could be adequate in applications. Central to the argumentation of the authors’ criticism is that modeling the response processes of students in a cross-sectional measurement (interindividual perspective) is erroneously confounded with response processes associated with a single person (intraindividual perspective; see [45]).

3. Model-Based Treatment of Missing Item Responses

In this section, different model-based approaches for the treatment of missing item responses are discussed. In the following, we only consider IRT models for dichotomous data in Rasch models [39,46]. However, our considerations are also applicable to more general model classes, such as 2PL models or IRT models for polytomous data [47].

Among the model-based procedures, IRT models are often used for distinguishing between ignorable and nonignorable item responses [48]. By ignoring missing item responses for students in the likelihood function, missing item responses are omitted in the estimation. If one ignores missing item responses in the estimation, it can be shown that the missing item responses can be imputed (i.e., scored) under this assumption with a probability of $P_{pi}\left({\mathsf{\theta }}_{p,M}\right)$ (i.e., scored) as correct. Here, ${\mathsf{\theta }}_{p,M}$ denotes the ability of person p that is being calculated exclusively with the observed item responses. This practically means that in the IRT model, observed item responses can be used to impute missing item responses. Since the probability $P_{pi}\left({\mathsf{\theta }}_{p,M}\right)$ is always larger than zero, ignoring missing item responses always leads to larger estimated trait values ${\mathsf{\theta }}_p$ than scoring them as incorrect. Ignoring missing item responses does not mean that the item nonresponse is independent of the latent trait (i.e., the missingness is not missing completely at random (MCAR); see [13]). The missing process is characterized by the fact that the whole of the information about the ability ${\mathsf{\theta }}_p$ is already recoverable from the observed item responses, i.e., the missing data are missing at random (MAR; [13]). Under this assumption, the proportion of missing item responses might be larger for students with lower abilities than for students with higher abilities.

The treatment of missing item responses as incorrect can be considered one extremum, with ignoring them as the opposing pole. In the first case, there is a risk of underestimating the ability of a person. In the second case, there is a risk of overestimating the ability of a person.

As an alternative to ignoring missing item responses, multidimensional IRT models for nonignorable item responses have been proposed [34,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63]. The two-dimensional model proposed in [48] introduces a latent individual response propensity (response tendency) $\mathsf{\xi }$ (besides the latent ability $\mathsf{\theta }$), which predicts missingness in item responses. For the response indicators $R_{pi}$, a Rasch model is assumed.

$$\mathsf{P}(R_{pi}=1|{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)=\mathsf{\Psi }({\mathsf{\xi }}_p-{\mathsf{\beta }}_i)$$

(2)

The probability of responding to an item, therefore, depends on the response tendency ${\mathsf{\xi }}_p$ and parameter ${\mathsf{\beta }}_i$. Whether an item is missing or not does not depend on the unknown item response itself. This means that $Y_{pi}$ and $R_{pi}$ are conditionally independent, given ${\mathsf{\theta }}_p$ and ${\mathsf{\xi }}_p$. The complete information for model parameters can be extracted from observed item responses and response indicators by specifying a particular parametric model. Formally, missing data are missing not at random (MNAR). However, missingness can be traced back to a further latent variable. This kind of missingness process is also denoted as latent ignorability ([64]; see also [31,65,66,67]).

Note that in a joint model for item responses $X_{pi}$ and response indicators $R_{pi}$, a bivariate distribution for the latent variable $(\mathsf{\theta },\mathsf{\xi })$ is estimated [15,48,68]. Alternative model-based procedures generalize the two-dimensional IRT model of Holman and Glas (see [25,69]) or use mixture distribution approaches [70,71]. It should be emphasized that for an individual person, the maximum likelihood estimate for the ability ${\mathsf{\theta }}_p$ does not exist; only the EAP estimate can be used [47]. As a consequence, a correlation of the ability $\mathsf{\theta }$ and the response tendency $\mathsf{\xi }$ has to be known to estimate individual ability values.

If the response tendency ${\mathsf{\xi }}_p$ can be expressed as ${\mathsf{\xi }}_p=\mathsf{\rho }{\mathsf{\theta }}_p+{\epsilon }_p$, one can show that missing item responses are approximately scored as $P_{pi}\left({\mathsf{\theta }}_{p,M}\right)-\mathsf{\rho }$ (see [69]). In contrast to the model with ignorable item responses, the ability estimates are therefore adjusted by a constant $\mathsf{\rho }$. Only if the response tendency and the ability are uncorrelated (i.e., $\mathsf{\rho }=0$) do both models lead to the same ability estimates.

4. Two Alternative Item Response Models for Nonignorable Item Responses: Approaches for a Sensitivity Analysis

In this section, two alternative model-based approaches for nonignorable missing item responses are discussed. In the first model, missing item responses are scored as partially correct [72]. In the second IRT model, the two-dimensional IRT model proposed in [48] is extended for the potential dependence of the missingness of an item from the unknown item response itself. In this model, the treatment of missing item responses as incorrect turns out to be a particular case. The two proposed models can be seen as a sensitivity analysis in which different assumptions about the missingness process can be investigated (e.g., [73]).

4.1. Pseudo-likelihood Approach for Partially Correct Scoring of Missing Item Responses

In the pseudo-likelihood approach of Lord [72], missing item responses are scored as partially correct. The pseudo-likelihood function $L_p$ (more precisely, the pseudo-log-likelihood function) for person p is defined as

$$log{L}_p=\sum _{i=1}^I\left[w_{pi}log{P}_{pi}+(1-w_{pi})log(1-P_{pi})\right],$$

(3)

where $P_{pi}$ is the probability of a correct answer to item i for person p, and $w_{pi}$ is the score of person p on item i. In the case of non-missing item responses, the score $w_{pi}$ is equal to either one (for a correct answer) or zero (for an incorrect answer). In [72], it is argued that in the case of missing item responses for multiple-choice items with M alternatives, the scores $w_{pi}$ can be set to $1/M$, and the resulting pseudo-likelihood function (3) is maximized. As argued above, by ignoring missing item responses, the values of the missing responses are be scored by $w_{pi}=P_{pi}\left({\mathsf{\theta }}_{p,M}\right)$, where ${\mathsf{\theta }}_{p,M}$ is the ability estimate for a person p based on the non-missing item responses. If the missing item is treated as incorrect, a score of $w_{pi}=0$ would be chosen. These should be regarded as two extreme cases of missing item treatment responses. Interim assumptions can be expressed by a sensitivity parameter $\mathsf{\rho }=0$ (scoring missings responses as incorrect) and $\mathsf{\rho }=1$ (treating missing responses as ignorable) that can be examined in the context of a sensitivity analysis. This approach defines the score as $w_{pi}=\mathsf{\rho }{P}_{pi}\left({\mathsf{\theta }}_{p,M}\right)$, given a preliminary ability estimate ${\mathsf{\theta }}_{p,M}$. Then, ability estimates are derived as a function of the parameter $\mathsf{\rho }$. Note that one can also use a marginal pseudo-likelihood function in (3) in which a distribution for ability values $\mathsf{\theta }$ is assumed. This technique of sensitivity analysis is frequently used for imputations of nonignorable missing data [73,74]. It should be noted that the parameter $\mathsf{\rho }$ itself is not estimable but must be fixed in the estimation. This is also supported by the fact that the maximization of the pseudo-likelihood function with respect to $\mathsf{\rho }$ can only result in $\mathsf{\rho }=0$ or $\mathsf{\rho }=1$.

The likelihood function (3) can be written as

$$L_p=\prod _{i=1}^I{L}_{pi}=\prod _{i=1}^I\left\{P_{pi}^{w_{pi}}{(1-P_{pi})}^{1-w_{pi}}\right\}.$$

(4)

Consequently, the likelihood for person p on item i contributes with $L_{pi}=P_{pi}^{w_{pi}}{(1-P_{pi})}^{1-w_{pi}}$. The missing item responses can therefore be regarded as fuzzy data with values of zero and one and an associated membership function $m_{pi}\left(0\right)=1-w_{pi}$ or $m_{pi}\left(1\right)=w_{pi}$ [29]. The multiplicative term $L_{pi}$ is also denoted as partial membership (see [75]). Alternatively, for fuzzy data, an additive likelihood term according to $L_{pi}=w_{pi}{P}_{pi}+(1-w_{pi})(1-P_{pi})$ has been proposed (see [29]), which is also denoted as mixed membership ([75]).

When estimating the IRT model according to the pseudo-likelihood approach, scores $w_{pi}$ are required in Equation (3). The concrete implementation of my proposed pseudo-likelihood approach proceeds as follows. The scores are computed for missing item responses based on an initial first scaling step. In this first step, missing item responses are treated as ignorable variables, and individual personal ability estimates are obtained (e.g., weighted likelihood estimates, see [76]), given the observed data and predetermined item difficulties. Alternatively, plausible values can be drawn from the individual posterior distribution [20]. The pseudo-likelihood approach (3) with fixed weights $w_{pi}$ can be estimated with the usual EM algorithm (see, e.g., [21,77]). In this approach, the M step remains unchanged compared to evaluations in usual IRT models. In the evaluation of the individual likelihood and the expected counts, the pseudo-likelihood is used instead of an ordinary likelihood.

4.2. Modeling the Missing Response Process

In the second IRT model, missing item responses are handled similarly to the two-dimensional model by Holman and Glas ([48]; see also [15]). Let us again denote $\mathsf{\Psi }$ as the logistic function so that in the Rasch model [47], the probability of a correct item response is expressed as

$$\mathsf{P}(Y_{pi}=1|{\mathsf{\theta }}_p)=\mathsf{\Psi }({\mathsf{\theta }}_p-b_i).$$

(5)

The probability for a missing item response of item i conditional on the response propensity $\mathsf{\xi }$ may now also be affected by the unknown item response $Y_i$ itself (see [6,7,78]). Hence, we define

$$\mathsf{P}(R_{pi}=1|Y_{pi}=k,{\mathsf{\xi }}_p)=\mathsf{\Psi }({\mathsf{\xi }}_p-{\mathsf{\beta }}_i-k\mathsf{\delta })$$

(6)

for $k=0,1$. The probability of a non-missing item response (i.e., $R_{pi}=1$) depending on the unknown item response $Y_i$ can be modeled with $\mathsf{\delta }\ne 0$. If $\mathsf{\delta }=0$ is chosen, the IRT model for nonignorable item responses in [48] is obtained. It has been shown that the parameter $\mathsf{\delta }$ can be identified from data ([31,79,80,81]; see also [82] for a similar model). Scaling results can also be studied in a sensitivity analysis as a function of a sequence of $\mathsf{\delta }$ values. For very small values of $\mathsf{\delta }$ (e.g., $\mathsf{\delta }=-10$), one obtains $\mathsf{P}(R_{pi}=1|Y_{pi}=1,{\mathsf{\xi }}_p)=1$ and, subsequently, $\mathsf{P}(R_{pi}=0|Y_{pi}=1,{\mathsf{\xi }}_p)=0$. This means that all students who know the item would provide a correct non-missing item response with a probability of 1. Using the Bayes theorem, it follows that $\mathsf{P}(Y_{pi}=1|R_{pi}=0;{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)=0$. Hence, all missing item responses are evaluated as incorrect for a very small value of $\mathsf{\delta }$. As in the pseudo-likelihood approach, one can use the proposed model to investigate the extremes of ignoring missing item responses but allowing the dependence on the latent response propensity ($\mathsf{\delta }=0$) and scoring missing item responses as incorrect ($\mathsf{\delta }=-10$).

In the proposed IRT model, item difficulties $b_i$, item parameter ${\mathsf{\beta }}_i$ for response tendency, and the bivariate distribution of $(\mathsf{\theta },\mathsf{\xi })$ are estimated. In the calculation of the likelihood for item i, the probabilities of three disjunctive events—$\mathsf{P}(R_{pi}=1,Y_{pi}=0|{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)$, $\mathsf{P}(R_{pi}=1,Y_{pi}=1|{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)$, and $\mathsf{P}(R_{pi}=0|{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)$—are used. The probability of a missing item response $\mathsf{P}(R_{pi}=0|{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)$ can be computed as the total probability:

$$\mathsf{P}(R_{pi}=0|{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)=\mathsf{P}(R_{pi}=0|Y_{pi}=0;{\mathsf{\xi }}_p)\mathsf{P}(Y_{pi}=0|{\mathsf{\theta }}_p)+\mathsf{P}(R_{pi}=0|Y_{pi}=1;{\mathsf{\xi }}_p)\mathsf{P}(Y_{pi}=1|{\mathsf{\theta }}_p),$$

(7)

  which can be calculated based on the item parameters and response models (5) and (6). Hence, missing item responses $Y_{pi}$ are integrated analogously to missing items in large-scale assessments studies by employing the integrated marginal maximum likelihood estimation method ([21]; see also [6,83]).

The proposed IRT model is based on the evaluation of item response probabilities conditional on the latent variables $\mathsf{\theta }$ and $\mathsf{\xi }$. For observed correct item responses, the probability is expressed as

$$\mathsf{P}(Y_{pi}=1,R_{pi}=1|{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)=\mathsf{P}(R_{pi}=1|Y_{pi}=1,{\mathsf{\xi }}_p)\mathsf{P}(Y_{pi}=1|{\mathsf{\theta }}_p)=\mathsf{\Psi }({\mathsf{\xi }}_p-{\mathsf{\beta }}_i-\mathsf{\delta })\mathsf{\Psi }({\mathsf{\theta }}_p-b_i).$$

(8)

For observed incorrect item responses, one obtains

$$\mathsf{P}(Y_{pi}=0,R_{pi}=1|{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)=\mathsf{P}(R_{pi}=1|Y_{pi}=0,{\mathsf{\xi }}_p)\mathsf{P}(Y_{pi}=0|{\mathsf{\theta }}_p)=\mathsf{\Psi }({\mathsf{\xi }}_p-{\mathsf{\beta }}_i)\mathsf{\Psi }(-{\mathsf{\theta }}_p+b_i)$$

(9)

using $1-\mathsf{\Psi }\left(x\right)=\mathsf{\Psi }(-x)$. For missing item responses, according to Equation (7), one obtains

$$\mathsf{P}(R_{pi}=0|{\mathsf{\theta }}_p,{\mathsf{\xi }}_p)=\mathsf{\Psi }(-{\mathsf{\xi }}_p+{\mathsf{\beta }}_i)\mathsf{\Psi }(-{\mathsf{\theta }}_p+b_i)+\mathsf{\Psi }(-{\mathsf{\xi }}_p+{\mathsf{\beta }}_i+\mathsf{\delta })\mathsf{\Psi }({\mathsf{\theta }}_p-b_i).$$

(10)

The two-dimensional IRT model can be estimated using an EM algorithm. A two-dimensional distribution for $(\mathsf{\theta },\mathsf{\xi })$ is estimated. Frequently, the bivariate normal distribution is utilized. However, this assumption can be weakened [68].

5. Comparison of Four Countries in PIRLS 2011

In the following section, a comparison of reading literacy performance in PIRLS 2011 is conducted using a test booklet for four selected countries: Austria (AUT), Germany (GER), France (FRA), and the Netherlands (NLD). The dependence of the country means is investigated for different treatments of the missing item responses.

5.1. Data

In the following analysis, item responses of booklet 13 in PIRLS 2011 (i.e., the “PIRLS Reader”) consisting of 35 items (15 multiple-choice items with four response alternatives; 20 constructed response items) were used. For this booklet, item responses of 968 Austrian, 809 German, 901 French, and 802 Dutch students were available. For simplicity, all polytomous items were dichotomized, and only the highest scores were recoded as correct. The dataset has been made available as data.pirlsmissing in the R [84] package sirt [85]. Student sampling weights were taken into account in the analyses.

Descriptive analyses show that the average proportion of missing item responses varied considerably between items and countries (AUT: 11.2%; DEU: 7.9%; FRA: 13.6%; NLD: 2.7%). Missing proportions for constructed response items were much higher than for multiple-choice items (e.g., for AUT, constructed response: 17.7%; multiple-choice: 2.6%). When scoring missing item responses as incorrect, Austrian students provided, on average, 55.1% correct item responses and slightly surpassed France (53.7%) but achieved significantly lower results than German students (63.0%) and Dutch students (64.4%).

5.2. Analysis

For the country comparison, six different IRT models (referred to as models M1, M2, M3, M4, M5, and M6) were specified. Student weights were taken into account in marginal maximum likelihood estimation. For a simplified description of the results, country means in each IRT model were transformed so that Austrian (AUT) students had a country mean of 500 and a country standard deviation of 100. The trait distributions for Germany (DEU), France (FRA), and the Netherlands (NLD) were subjected to the same transformation so that country comparisons in all IRT models should always be interpreted relative to the performance of Austria. This choice of standardization was made because Austria was considered a reference country in the analysis, and the other three countries should be compared to Austria (see [86,87]). However, there could also be good reasons for using a fixed mean and a fixed standard deviation in a total population comprising all countries (see [31,88]).

In model M1 (that corresponds to the official procedure in PIRLS), a unidimensional Rasch model with four groups (i.e., the four countries AUT, GER, FRA, and NLD) was specified in which missing item responses were scored as incorrect. In model M2, missing item responses were ignored, which means that they were omitted in the computation of the marginal likelihood. As an extension to model M2, model M3 was a two-dimensional Rasch model with two latent variables of ability $\mathsf{\theta }$ and the response tendency $\mathsf{\xi }$. In model M4, a pseudo-likelihood approach of the Rasch model was specified in which missing item responses for multiple-choice items with four answer alternatives were scored with the score $w_{pi}=0.25$. In contrast, missing item responses for constructed response items were scored as incorrect (i.e., with a score of $w_{pi}=0$). In models M1, M2, M3, and M4, item parameters and the distribution parameters for the trait(s) in the four countries were concurrently estimated.

In models M5 and M6, in a first analysis, missing item responses were scored as incorrect, and common item parameters for all four countries were obtained, which were subsequently fixed in a second step. In model M5, the pseudo-likelihood approach was used for the Rasch model, in which missing item responses were scored as $w_{pi}=\mathsf{\rho }{P}_{pi}\left({\mathsf{\theta }}_{p,M}\right)$, where ${\mathsf{\theta }}_{p,M}$ is the personal ability, which was computed based on the non-missing item responses. The personal ability ${\mathsf{\theta }}_{p,M}$ was a random draw from the individual posterior distribution. The country means were calculated depending on the sensitivity parameter $\mathsf{\rho }=0,0.01,\dots ,0.99,1$. In model M6, the two-dimensional model M3 is extended by modeling the missingness process dependence of the absence on the item response itself. This model was calculated depending on the sensitivity parameter $\mathsf{\delta }=-10,-9.5,\dots ,-0.5,0$.

Additionally, we empirically determined the imposed missingness mechanism defined by Equation (6) for all four countries. In this model, the correlation $\mathsf{\rho }$ between the latent trait $\mathsf{\theta }$ and the response propensity $\mathsf{\xi }$ can be fixed to zero or can be estimated (i.e., $\mathsf{\rho }\ne 0$). Moreover, the parameter $\mathsf{\delta }$ can be fixed to 0 (i.e., ignoring the missing item response), fixed to $-10$ (i.e., treating the missing item response as incorrect), or be estimated from the data assuming a common $\mathsf{\delta }$ parameter across items. Hence, five different model specifications of the two-dimensional IRT model result. The analysis models were defined as N1: $\mathsf{\delta }=0$, $\mathsf{\rho }=0$; N2: $\mathsf{\delta }=-10$, $\mathsf{\rho }=0$; N3: $\mathsf{\delta }$ estimated, $\mathsf{\rho }=0$; N4: $\mathsf{\delta }=0$, $\mathsf{\rho }$ estimated; N5: $\mathsf{\delta }=-10$, $\mathsf{\rho }$ estimated; and N6: $\mathsf{\delta }$ estimated, $\mathsf{\rho }$ estimated. The models were compared using the Bayesian information criterion (BIC). Because the models were separately estimated for each country, the distributions of the latent trait $\mathsf{\theta }$ were subsequently brought onto a common metric across countries using the mean–mean linking method [89].

For all item response models, the country means of Austria, Germany, France, and the Netherlands were compared. The entire data preparation and estimation of the models was conducted in the R software [84]. For estimation of models M1, M2, and M3, the R package TAM [90] was used. The models M4, M5, M6, and N1 to N6 were estimated with the rasch.mml2() function of the R package sirt [85].

5.3. Results

In

Table 1. Country means for Austria (AUT), Germany (GER), France (FRA), and the Netherlands (NLD) in PIRLS 2011 for different treatments of missing item responses.

Model	AUT	GER	FRA	NLD
M1: missing = incorrect	500	537.5	488.7	540.3
M2: missing = ignorable	500	534.2	492.4	523.4
M3: 2-dim. model	500	534.9	492.5	524.8
M4: pseudo-likelihood  (for multiple-choice items)	500	537.6	489.4	539.5
M5: pseudo-likelihood
$\mathsf{\rho }=0\phantom{.3}$	500	537.3	488.9	539.9
$\mathsf{\rho }=0.3$	500	537.0	490.1	535.9
$\mathsf{\rho }=0.7$	500	535.9	491.8	529.5
$\mathsf{\rho }=1\phantom{.3}$	500	534.6	493.1	524.0
M6: 2-dim. model
$\mathsf{\delta }=-10\phantom{.}$	500	538.0	489.1	540.7
$\mathsf{\delta }=-1.5$	500	535.9	490.6	532.4
$\mathsf{\delta }=-0.5$	500	535.1	491.5	528.0
$\mathsf{\delta }=0\phantom{-.3}$	500	534.6	492.1	525.7

Note. The mean and the standard deviation for Austria (AUT) were fixed at 500 and 100 in all item response models, respectively.

, the country means for all models (M1 to M6) are shown. It is noteworthy that the two extreme treatments of missing item responses in model M1 (missing item responses were scored as incorrect) and model M2 (missing item responses were ignored) resulted in slight differences for Germany and France in comparison to the country mean of Austria (GER: M1: 537.5, M2: 534.2; FRA: M1: 488.7, M2: 492.4), while the different model specification led to significant differences for the Netherlands (NLD: M1: 540.3, M2: 523.4). These findings can be explained by the fact that the proportion of missing item responses for Dutch students (2.7%) was substantially lower than that for Austrian students (11.2%).

In models M5 and M6, the sensitivity parameters $\mathsf{\rho }$ and $\mathsf{\delta }$ are varied within the framework of sensitivity analysis that refers to different assumptions of the missing process for item responses. In practice, the findings presented in Table 1 span the range of the country means between the extremes of treating missing responses as incorrect (model M1) and treating them as ignorable (model M2).

In Figure 1, the country means are shown as a function of the sensitivity parameters in models M5 and M6. The ranks of Germany and the Netherlands switch depending on the parameters $\mathsf{\rho }$ and $\mathsf{\delta }$. Notably, the country means are continuous and monotone functions of the sensitivity parameters $\mathsf{\rho }$ and $\mathsf{\delta }$, respectively.

In

Table 2. Model comparison based on the Bayesian information criterion (BIC) for Austria (AUT), Germany (GER), France (FRA), and the Netherlands (NLD) in PIRLS 2011.

Model	AUT	DEU	FRA	NLD
N1: $\mathsf{\delta }=0$, $\mathsf{\rho }=0$	47,741	36,366	45,029	33,142
N2: $\mathsf{\delta }=-10$, $\mathsf{\rho }=0$	47,827	36,414	45,263	33,144
N3: $\mathsf{\delta }$ estimated, $\mathsf{\rho }=0$	47,722	36,365	45,028	33,130
N4: $\mathsf{\delta }=0$, $\mathsf{\rho }$ estimated	47,677	36,285	44,888	33,127
N5: $\mathsf{\delta }=-10$, $\mathsf{\rho }$ estimated	47,790	36,355	45120	33134
N6: $\mathsf{\delta }$ estimated, $\mathsf{\rho }$ estimated	47,666	36,288	44,887	33,120

Note. The parameter $\mathsf{\rho }$ is the correlation between the latent trait $\mathsf{\theta }$ and the latent response propensity $\mathsf{\xi }$. The parameter $\mathsf{\delta }$ indicates the dependency of the item response indicator $R_{pi}$ from the (unobserved) item response $X_{pi}$ (see Equation (6)). The model with the lowest BIC is printed in bold.

, the BIC is displayed for all five estimated models (N1, …, N5) and the four countries (AUT, DEU, FRA, and NLD). It turned out that for three of the four countries (i.e., AUT, FRA, and NLD), the most general model was preferred that estimated both missingness parameters $\mathsf{\rho }$ and $\mathsf{\delta }$. Hence, the latent response propensity is not fully sufficient to explain missing item responses. The true but unobserved item response itself also governs the occurrence of missing item responses. Only for DEU, model N3, which includes a correlation $\mathsf{\rho }$ between the latent trait and the latent response propensity but does not allow an influence of the item response on the response indicator, was the best-fitting model in terms of the BIC.

From model N6, country-specific missingness parameters can be extracted. The latent trait $\mathsf{\theta }$ and the latent response $\mathsf{\xi }$ moderately correlated positively (AUT: $\mathsf{\rho }=0.31$; DEU: $\mathsf{\rho }=0.40$; FRA: $\mathsf{\rho }=0.47$; NLD: $\mathsf{\rho }=0.23$). The parameter $\mathsf{\delta }$ was smaller than zero, indicating that students who omitted items did so because they did not know the item (AUT: $\mathsf{\delta }=-0.96$; DEU: $\mathsf{\delta }=-0.75$; FRA: $\mathsf{\delta }=-0.57$; NLD: $\mathsf{\delta }=-2.02$). Interestingly, the $\delta$ parameter was slightly more negative in model N3, in which the correlation $\mathsf{\rho }$ between the trait and the response propensity was set to zero (AUT: $\mathsf{\delta }=-1.20$; DEU: $\mathsf{\delta }=-1.08$; FRA: $\mathsf{\delta }=-0.58$; NLD: $\mathsf{\delta }=-2.39$). Importantly, these findings demonstrate that the response mechanisms differ across countries.

In

Table 3. Country means for Austria (AUT), Germany (GER), France (FRA), and the Netherlands (NLD) in PIRLS 2011 for different model specifications of the missingness mechanism.

Model	AUT	DEU	FRA	NLD
N1: $\mathsf{\delta }=0$, $\mathsf{\rho }=0$	500	535.4	494.1	526.9
N2: $\mathsf{\delta }=-10$, $\mathsf{\rho }=0$	500	537.6	489.8	542.0
N3: $\mathsf{\delta }$ estimated, $\mathsf{\rho }=0$	500	537.1	497.9	530.9
N4: $\mathsf{\delta }=0$, $\mathsf{\rho }$ estimated	500	535.1	493.3	526.9
N5: $\mathsf{\delta }=-10$, $\mathsf{\rho }$ estimated	500	537.6	489.9	542.0
N6: $\mathsf{\delta }$ estimated, $\mathsf{\rho }$ estimated	500	537.4	495.8	530.4

Note. The mean and the standard deviation for Austria (AUT) were fixed at 500 and 100 in all item response models, respectively.

, country means for the five different model specifications (N1 to N5) are presented. The largest relative changes across models occur for NLD because this country had the lowest missing item response rates. The models do not substantially differ regarding the question of whether the correlation $\mathsf{\rho }$ between the trait $\mathsf{\theta }$ and the response propensity $\mathsf{\theta }$ is estimated or fixed to zero.

6. Discussion

In contrast to statements in the literature [17,24], our analyses of PIRLS 2011 showed that ignoring missing item responses has an impact on the central results of a study. Country means varied considerably in different treatments of missing item responses, which raises the question of the choice of a valid scaling approach.

We argued that the criticism against the scoring of missing item responses as incorrect based on test-theoretical arguments [17,25] utilizes an intraindividual interpretation of item response probabilities. In this reasoning, the concepts of aleatoric uncertainty (the statistical model) and epistemic uncertainty (concerning the assessment of data) are erroneously confounded. It was also argued that simulation studies to choose the appropriate way to handle missing item responses are irrelevant to choosing a scaling method (see also [30]). Since defining an adequate scoring rule for missing responses is not empirically determinable, two different IRT models were proposed, which served as a sensitivity analysis of different assumptions about missing item responses.

The first alternative one-dimensional IRT model (model M5) is based on a pseudo-likelihood approach in which the missing item responses can take values between zero and one. However, this approach has the disadvantage that for item responses, weights $w_{pi}$ have to be computed. In the second model (model M6), the additional latent response tendency is introduced, which can be correlated with ability. The missingness of an item response is allowed to be dependent on the unknown item response itself.

Typically, the sensitivity analyses performed in this paper are preferred to more complex model specifications, such as pattern mixture models [34]. When imputing missing data, a range of analysis results can be obtained by varying plausible assumptions on the missingness [74]. If a prior distribution for a sensitivity parameter characterizing the assumptions of the missingness process is specified, then a multimodel inference analysis can be conducted that weighs different theoretical assumptions about the missingness process [91].

Within the scope of our analysis, we restricted ourselves to the Rasch model. However, the analysis can be extended to more complex IRT models such as the 2PL, the 3PL, or the generalized partial credit model. Instead of using a common sensitivity parameter for all items, this parameter can be separately estimated for each individual item or for individual item groups (such as multiple-choice items and constructed response items; see [31,80]). If a country-specific missingness process is estimated (i.e., by assuming country-specific $\delta$ parameters), the comparability of the country means obtained from such IRT models could be strongly questioned.

As an alternative to model-based treatments of missing item responses, missing values on items can be multiply imputed [92,93,94]. The advantage of imputation-based procedures might be that more complicated (e.g., multidimensional) IRT models might be used for imputation, while the analysis model is simpler and typically mis-specified (see [31,95]).

Undoubtedly, the best strategy for handling missing data would be to avoid their appearance [96]. In computerized administered tests, missing item responses can be avoided by forcing students to respond to an item. However, careless or unmotivated responses cannot be excluded by such modifications. Students who do not know a constructed response item will unlikely provide a valid item response if forced by the computer.

We would like to emphasize that we think that the choice of how to score missing item responses should mainly be made under from a validity perspective [97,98]. That is, a particular missing data treatment implies a scoring rule for obtaining ability estimates. The validity of the scoring rule must be discussed regarding the use of the score and its interpretations [99]. In this sense, we think that scoring missing item responses as incorrect is the most valid strategy for country comparisons in official reporting of large-scale educational assessment studies such as PIRLS but may be inappropriate for research questions involving relations of student abilities and questionnaire scales such as motivation or self-concept [32]. Hence, there cannot be a general optimal rule for how to handle missing item responses in all potential applications of a dataset. Our motivation for writing this article was that we find scoring strategies inappropriate that do not score missing item responses as incorrect (see [5]) in studies that involve reporting and comparison of several states or countries.

As pointed out by an anonymous reviewer, missing item responses should be classified into omitted and not-reached item responses [100]. Our arguments mainly apply to omitted item responses within a test. For speed tests, it might be preferable not to score not-reached item responses as incorrect. However, large-scale assessment studies such PIRLS are not strongly speeded, and we also think that not-reached items should be scored as incorrect [32]. However, the recently introduced use of multistage adaptive testing in educational assessment studies complicates the meaning of not-reached items across students [101].

The consequences of treating missing item responses as ignorable or incorrect can also be significant for scaling in longitudinal studies if missing proportions vary by grade, school track, or time point. In IRT models that assume ignorable missing item responses, the meaning of the identified ability remains unclear. If an ability in a particular domain is just defined by the correctly solved items in a test, in our view, a non-circular and more valid definition of ability is obtained. Consequently, it can be concluded that missing item responses cannot be ignorable or partially ignorable.

7. Conclusions

In this article, we argued that treating missing item responses as ignorable or partially ignorable in IRT models involving a latent response propensity does not necessarily lead to valid substantive conclusions. To handle such nonignorable item responses, we proposed two types of statistical sensitivity analyses that include scoring missing item responses as incorrect and the treatment as ignorable as two extreme cases. By doing so, a range of plausible parameter estimates can be obtained. The two types of sensitivity analyses were illustrated using PIRLS 2011 with four countries. It has been shown that assumptions about the missing data mechanism did impact the country means and the ranking, as well as differences between countries.