A Bayesian Item Response Model for Examining Item Position Effects in Complex Survey Data

Abstract

A multidimensional Bayesian item response model is proposed for modeling item position effects. The first dimension corresponds to the ability that is to be measured; the second dimension represents a factor that allows for individual differences in item position effects called persistence. This model allows for nonlinear item position effects on the item side as well as on the person side. Moreover, a flexible loading structure on the two dimensions is allowed. A fully Bayesian estimation procedure is proposed, and its performance is investigated by a simulation study. Further, the model is applied to empirical data collected in the Programme for International Student Assessment 2000 in the reading domain. The additional value of the model’s extended flexibility compared to more restrictive models is shown. The findings show that the linear hypothesis of change in performance during a test does not hold in general.

Keywords

Bayesian IRT MCMC estimation complex survey data multidimensional IRT item position effects

Introduction

Educational testing generally aims to measure individual traits or states that represent abilities or competencies in various domains of education. In order to derive these measures, a sequence of items is usually administered to persons. In an ideal situation of educational testing, the behavior a person shows in responding to items within a test is influenced only by their trait or state, which is what the test intended to measure. However, in realistic situations of educational testing, response behavior is influenced by several other factors besides the trait or state of the testees. For instance, the behavior a person shows in responding to items within a test that consists of several items depends on the position within this test in which the item is administered to that person. Almost seven decades ago, Mollenkopf (1950) investigated the influence of item positions on the response behavior of persons in speed and power tests. He already found that items tended to be more difficult when the item appeared later in the test rather than in earlier item positions.

In the literature, several approaches based on item response theory (IRT) have been proposed for analyzing item position effects. Item position effects can be viewed from the item side and from the person side (De Boeck & Wilson, 2004; Hartig & Buchholz, 2012). Viewed from the item side, the position effects of items are modeled as additional effects on item parameters and, therefore, are constant across persons. In the Rasch model, item position effects can be analyzed using the linear logistic test model (Alexandrowicz & Matschinger, 2008; Hohensinn et al., 2008; Kubinger, 2008). Viewed from the person side, test takers may be affected differently by item position effects. To allow for individual differences, item position effects can be specified in such a way that the item position loads on a factor that can vary across persons (Debeer & Janssen, 2013; Hartig & Buchholz, 2012; Schweizer et al., 2009). The effect, that is the loading of the item position, is usually held constant across items. The factor that the item position loads on is often called persistence (e.g., Debeer & Janssen, 2013). It should be noted that, to analyze position effects on the item side, items have to be administered to the test takers in different orders (Hartig & Buchholz, 2012). However, to analyze position effects on the person side, tests with a fixed item order (i.e., using only one test booklet; see Schweizer et al., 2009) can be used and individual differences in persistence can also be identified. In this article, we extend the multidimensional model for item position effects proposed by Debeer and Janssen (2013) by allowing for a more flexible functional form of persistence and for a loading structure of the items.

The remainder of the article is organized as follows. First, the model is introduced and its main features are summarized. Second, the identifiability of the parameters under ideal conditions is discussed, that is, for data arising from a test design in which every item occurs in every position and in which there are substantial links of items in different positions within a booklet. Third, a Bayesian estimation procedure to fit the proposed model is described. Fourth, the statistical properties of the Bayesian parameter estimates are evaluated in terms of bias, root mean square error of approximation, and coverage rates, calculated based on the results of a simulation study. Fifth, the model is applied to data collected from the first seven booklets, including the first three block positions, of the Programme for International Student Assessment (PISA) 2000. Finally, further extensions of the model and prospective research directions are discussed.

Modeling Item Positions in an IRT Framework

Item Response Models for Item Position Effects

We assume a test consisting of items $i = 1, . . ., I$ organized in different booklets that are administered to persons $p = 1, . . ., N$ . The booklet design is supposed to form a block design, that is, a booklet consists of $t = 1, . . ., T$ sections, also referred to as blocks. Furthermore, a block consists of more than 1 item, and 1 item is assigned to exactly one block (an item always has the same position in its block). Note that an item of a booklet has a fixed position t, which is defined by the position of its block, and an item is administered to a person in one position only. Table 1 illustrates an example for such a booklet design.

Table 1.

Example for a Booklet Design Consisting of Four Booklets

Booklet	Block Position t
	1	2	3	4
1	B1	B2	B3	B4
2	B2	B3	B4	B1
3	B3	B4	B1	B2
4	B4	B1	B2	B3

Note. There are four blocks (B1, B2, B3, and B4), and no two blocks have one item in common.

Let $Y_{p i t}$ denote the binary random variable representing the response of person p to Item i when administered in position t. A model for analyzing item position effects can be formulated as follows:

Φ^{- 1} [P (Y_{p i t} = 1)] = a_{i 1} θ_{p 1} + a_{i 2} γ_{t} θ_{p 2} + b_{i} + β_{t},

where b_i is a fixed effect representing the item-specific difficulty. Note that we used the $probit$ function $Φ^{- 1} (x)$ as link function for computational reasons, where $Φ$ is the standard normal distribution function. The parameter $β_{t}$ represents the position-specific mean change in difficulty of position t (Trendtel & Robitzsch, 2018) and is allowed to deviate from linearity (see Le, 2009). $θ_{p 1}$ is the random effect representing the ability of person p. The random effect $θ_{p 2}$ captures the individual deviation of person p from the position-specific mean change in difficulty $β_{t}$ and is also referred to as the persistence of person p (e.g., Hartig & Buchholz, 2012). The parameters $a_{i 1}$ and $a_{i 2}$ are dimension-specific loadings. The parameters $a_{i 1}$ are also known as the item discrimination parameters, and $a_{i 2}$ can be interpreted as the intensity with which the performance in solving Item i is affected by the persistence $θ_{p 1}$ averaged over all positions. Additionally, $γ_{t}$ is a position-specific loading parameter and represents the intensity of the persistence’ influence at position t averaged over all items. The position effects $β_{t}$ and the loadings $a_{i 2}$ could also be specified as random effects (e.g., Fox & Verhagen, 2010), which could make the estimation easier.

The vector ${(θ_{p 1}, θ_{p 2})}^{T}$ follows a bivariate normal distribution with mean vector $μ_{θ_{1}, θ_{2}} = (μ_{θ_{1}}, μ_{θ_{2}})^{T} {= (0, 0)}^{T} = 0$ and covariance matrix

Σ_{θ_{1}, θ_{2}} = (\begin{matrix} σ_{θ_{1}}^{2} \\ σ_{θ_{1} θ_{2}} & σ_{θ_{2}}^{2} \end{matrix}),

where the parameter $σ_{θ_{1} θ_{2}}$ allows for dependencies between ability and persistence.

To ensure the identifiability of the model, $γ_{1}$ and $β_{1}$ were set to 0; $γ_{2}$ , $σ_{θ_{1}}$ , and $σ_{θ_{2}}$ were set to 1; and the loadings $a_{i 1}$ , $a_{i 2}$ , $γ_{t}$ are restricted to be greater than zero for all $i = 1, . . ., I$ . In doing so, Position 1 served as the reference position for the effects representing the mean change in difficulty, that is, $β_{t}$ . For the position-specific loadings, $γ_{t}$ , Position 2 served as the reference position.

The model in Equation 1 can be considered to be a generalized nonlinear mixed effects model (Debeer & Janssen, 2013; De Boeck & Wilson, 2004) with fixed item effects. However, a generalization to models with random item effects, resulting in a cross-classification, could easily be made (e.g., Fox & Verhagen, 2010). To estimate the parameters of such models becomes computationally demanding when using maximum likelihood methods. Therefore, we propose a Markov chain Monte Carlo (MCMC) approach that easily allows the model to be extended to even more complicated models.

The model in Equation 1 can be restricted in several ways, and three versions of these model restrictions are discussed below.

Regarding the item position effects, we could assume that these effects have a certain shape. For instance, we could assume that the effects can be described by a function $f (t)$ of position t. A function that often can be found in literature in the context of item position effects is the linear function (e.g., Debeer & Janssen, 2013; Hartig & Buchholz, 2012). This can be incorporated by restrictions regarding the parameters $γ_{t}$ and $β_{t}$ leading to

Φ^{- 1} [P (Y_{p i t} = 1)] = a_{i 1} θ_{p 1} + a_{i 2} (t - 1) θ_{p 2} + b_{i} + (t - 1) β .

The parameter $β$ can be considered to be the average change in performance from one item position to the next and $θ_{p 2}$ can be considered to be the individual deviation from this average as before, while the position-specific loading of position t is defined by $(t - 1)$ .

Further, the model in Equation 1 can be restricted to have constant item-specific loadings. The loading parameters are set to $a_{i 1} = 1$ and $a_{i 2} = 1$ for all $i = 1, . . ., I$ . This leads to

Φ^{- 1} [P (Y_{p i t} = 1)] = θ_{p 1} + γ_{t} θ_{p 2} + b_{i} + β_{t} .

In this model, $σ_{θ_{1}}$ , $σ_{θ_{2}}$ , and $σ_{θ_{1} θ_{2}}$ are free parameters.

Finally, the position effects can be modeled on the item side only. This means that the second person factor $θ_{p 2}$ is left out of the model, which leads to

Φ^{- 1} [P (Y_{p i t} = 1)] = a_{i 1} θ_{p 1} + b_{i} + β_{t} .

In this model, the person factor $θ_{p 1}$ follows a univariate normal distribution with mean of 0 and variance of 1. Obviously, some of the aforementioned restrictions could be combined leading to further versions of the model. In the next section, we discuss the identifiability of the model described in Equation 1.

Identification of the Parameters

We show that the parameters of the statistical model in Equation 1 can be identified after integrating out the random effects. The main principle of statistical identification is that unknown parameters can be uniquely estimated based on observed data only (see San Martin et al., 2013) and not by another source of information (like prior distributions, see San Martin & González, 2010). This means that parameters can be determined if the probability distribution of the observed data is fully known (i.e., there is no sampling error).

In the following, we assume that the model restrictions described in the previous section are made, that is, $σ_{θ_{1}}^{2} = σ_{θ_{2}}^{2} = γ_{2} = 1$ , $γ_{1} = β_{1} = μ_{θ_{1}} = μ_{θ_{2}} = 0$ , and $γ_{t}, a_{i 1}, a_{i 2} > 0$ for all $t > 2$ and all $i = 1, . . ., I$ .

Let the index $i (t), j (t^{'})$ denote the case that Item i and Item j are observed in positions t and $t^{'}$ , respectively. We show that all unknown model parameters can be calculated as functions of the observed univariate and bivariate distributions¹ $P (Y_{i (t)})$ and $P (Y_{i (t)}, Y_{j (t^{'})})$ . A proportion $p_{i (t)}$ can be obtained as the success probability of $Y_{i (t)}$ . A tetrachoric correlation $ρ_{i (t), j (t^{'})}$ can be computed from the bivariate distribution $P (Y_{i (t)}, Y_{j (t^{'})})$ , which is the basis of model identification (Muthén, 1979). The tetrachoric correlation can be written as functions of the item loadings $a_{i 1}$ , $a_{i 2}$ , and $γ_{t}$ and the correlation $σ_{θ_{1} θ_{2}}$ (Renard et al., 2004).

The identification proof consists of six steps. In each step, it is shown that a parameter or a group of parameters can be calculated as a function of the observed data distribution.

First, we show how the item loadings $a_{i 1}$ can be identified from the observed data. Using data for Items i and j in the first position, it holds for the tetrachoric correlation that $ρ_{i (1), j (1)} = a_{i 1} a_{j 1} {(a_{i 1}^{2} + 1)}^{- 1 / 2} {(a_{j 1}^{2} + 1)}^{- 1 / 2}$ (Renard et al., 2004). It follows that $v_{i j h} = ρ_{i (1), j (1)} ρ_{i (1), h (1)} ρ_{j (1), h (1)}^{- 1} = a_{i 1}^{2} {(a_{i 1}^{2} + 1)}^{- 1}$ . Hence, $a_{i 1} = v_{i j h}^{1 / 2} {(1 - v_{i j h})}^{- 1 / 2}$ and $a_{i 1}$ can be calculated as a function of the observed data (i.e., a function of $v_{i h j}$ ).

Second, we show how the item loadings $a_{i 2}$ can be identified using data for Item i in the second position and for Item j in the first position; it holds that $ρ_{j (1) i (2)} = a_{j 1} (a_{i 1} + a_{i 2} σ_{θ_{1} θ_{2}}) (a_{j 1}^{2} + {1)}^{- 1 / 2} {(a_{i 1}^{2} + 2 a_{i 1} a_{i 2} σ θ_{1} θ_{2} + a_{i 2}^{2} + 1)}^{- 1 / 2}$ . This equation has the unknowns $a_{i 2}$ and $σ_{θ_{1} θ_{2}}$ , whereas $ρ_{j (1) i (2)}$ is observed and $a_{i 1}$ and a_j are identified in the first step. Squaring the equation leads to a solution of a quadratic equation for $a_{i 2}$ as a function of $σ_{θ_{1} θ_{2}}$ . Hence, $a_{i 2}$ can be identified as long as $σ_{θ_{1} θ_{2}}$ can be identified.

Third, we identify the correlation $σ_{θ_{1} θ_{2}}$ . Using Items i and j both in the second position, it holds that $ρ_{i (2) j (2)} = (a_{i 1} a_{j} + a_{i 1} a_{j 2} σ_{θ_{1} θ_{2}} + a_{j 1} a_{i 2} σ_{θ_{1} θ_{2}} + a_{i 2} a_{j 2}) (a_{i 1}^{2} + 2 a_{i 1} a_{i 2} σ θ_{1} θ_{2} + a_{i 2}^{2} + {1)}^{- 1 / 2} {(a_{j 1}^{2} + 2 a_{j 1} a_{j 2} σ_{θ_{1} θ_{2}} + a_{j 2}^{2} + 1)}^{- 1 / 2} .$

This equation has the only unknown $σ_{θ_{1} θ_{2}}$ because $a_{i 2}$ and $a_{j 2}$ are functions in $σ_{θ_{1} θ_{2}}$ . It can be shown that this equation has one and only solution, $σ_{θ_{1} θ_{2}}$ , and, hence, this parameter is identified from the observed data.

Fourth, we show how to identify $γ_{t}$ ( $t \geq 3$ ). Using data for Item i in position t and an Item j in position 1, it holds that $ρ_{j (1) i (t)} = a_{j 1} (a_{i 1} + a_{i 2} γ_{t} σ_{θ_{1} θ_{2}}) (a_{j 1}^{2} + {1)}^{- 1 / 2} {(a_{i 1}^{2} + 2 a_{i 1} a_{i 2} γ_{t} σ θ_{1} θ_{2} + a_{i 2}^{2} γ_{t}^{2} + 1)}^{- 1 / 2}$ . This equation can be solved with respect to the unknown $γ_{t}$ as all other parameters in the equation have already been identified.

Fifth, we show how to identify item intercepts b_i . Using data from Item i in the first position, it holds that $p_{i (1)} = Φ [b_{i} {(a_{i 1}^{2} + 1)}^{- 1 / 2}]$ (Renard et al., 2004). Hence, $b_{i} = (a_{i 1}^{2} + {1)}^{1 / 2} Φ^{- 1} (p_{i (1)})$ can be identified from the data because $a_{i 1}$ has already been identified from data.

Sixth, we show how to identify the position parameters $β_{t}$ ( $t \geq 2$ ). Using data from Item i in position t, it holds that $p_{i (t)} = Φ [(b_{i} + β_{t}) (a_{i 1}^{2} + 2 a_{i 1} a_{i 2} γ_{t} σ_{θ_{1} θ_{2}} + a_{i 2}^{2} σ_{θ_{1} θ_{2}}^{2} + {1)}^{- 1 / 2}]$ . The only unknown in this equation is the parameter $β_{t}$ for which the equation has a unique solution.

A detailed analysis of the requirements of the six identification steps reveals that specific test designs are needed to identify all of the parameters. Such sufficient requirements can formulated as follows: All booklets consist of blocks containing more than 1 item. The booklets are designed such that the total number of block positions is greater than 2. For every item exist two booklets k and l such that the item appears in booklet k at Block Position 1 and in booklet l at Block Position 2. For every block position $t > 1$ exists a booklet containing an item pair $(i, j)$ for which holds that Item i is administered at Block Position 1 and Item j at Block Position t.

The booklets we selected in the empirical example from the PISA 2000 design (see Table 2) fulfill these requirements, and therefore, it is ensured that all of the univariate and bivariate distributions needed in the six steps could be estimated.²

Table 2.

Booklet Design of PISA 2000 Assessments and the Position-Specific Relative Frequencies of Correct Responses Regarding Reading Items and the First Three Block Positions Only

Block Position	Booklet									Percentage of Correct Responses
Block Position	1	2	3	4	5	6	7	8	9	Percentage of Correct Responses
1	R1	R2	R3	R4	R5	R6	R7	M	S	67.4
2	R2	R3	R4	R5	R6	R7	R1	S	M	61.9
3	R4	R5	R6	R7	R1	R2	R3	R8	R9	64.6
4	M	S	M	S	M	S	R8	R9	R8

Note. R1–R9 are clusters consisting of reading items. M and S represent clusters consisting of mathematics and science items, respectively. Cells displayed in gray were not considered for the analyses.

MCMC Estimation Procedure

In the following, a full Bayesian procedure to estimate the proposed model is described, namely, an MCMC procedure using the Gibbs sampler (Gelfand & Smith, 1990). To make the estimation procedure applicable to stratified clustered samples, we extended the usual MCMC procedure by incorporating sampling weights (e.g., Goldstein, 2011). The procedure was implemented in R (R Core Team, 2018). MCMC procedures are often a facilitation of estimation (see Gelman et al., 2013; Hoff, 2009; Levy & Mislevy, 2016), and, in general, researchers can implement MCMC methods relatively easily.

This is one reason for why the MCMC approach has found widespread use in the estimation of item response models (e.g., Fox, 2010; Rupp et al., 2004). Using the Gibbs sampling approach, a complicated high-dimensional estimation problem can be reduced to stepwise evaluations of lower dimensional posterior distributions (which are often unidimensional). However, compared to maximum likelihood estimation, the Bayesian approach additionally requires the specification of prior distributions (with densities denoted by $π$ ) for all parameters (Hoff, 2009).

Our Gibbs sampler employed augmented latent data z which can be interpreted as continuous versions of corresponding dichotomous item responses y (Albert & Chib, 1993). The use of augmented data simplifies the Gibbs sampling steps (Levy & Mislevy, 2016). By incorporating the augmented data, the joint posterior density of the parameters given the data is specified as

\begin{array}{l} p (a_{1}, b, a_{2}, γ, β, θ_{1}, θ_{2}, z, Σ_{θ_{1} {,θ}_{2}} | y) \\ \propto & [\prod_{p, i (t)} P (Y_{p i t} = y_{p i t} | Z_{p i t}) P (Z_{p i t} | θ_{p 1}, θ_{p 2}, a_{i 1}, a_{i 2}, b_{i}, γ_{t}, β_{t})] \times \\ [\prod_{p} p (θ_{p 1}, θ_{p 2} | Σ_{θ_{1} {,θ}_{2}})] \times p (a_{1}, b) p (a_{2}) p (γ) p (β) p (Σ_{θ_{1} {,θ}_{2}}), \end{array}

where the proportional sign is used because a normalizing constant is left out of the equation.

For notational reasons, let $ω$ denote the vector consisting of all model parameters except $Σ_{θ_{1} {,θ}_{2}}$ with density function $π (ω)$ of the corresponding joint prior distribution. Using this notation and integrating out the augmented data as well as the factors $θ_{1}$ and $θ_{2}$ in Equation 6, the joint posterior simplifies to

p (ω | y) \propto [\prod_{p} \int P (y_{p} | z_{p}) P (z_{p} | θ_{p 1}, θ_{p 2}, ω) P (θ_{p 1}, θ_{p 2} | Σ_{θ_{1}, θ_{2}}) d z_{p} d θ_{p 1} d θ_{p 2}] \times π (ω) π (Σ_{θ_{1}, θ_{2}}) .

It is easily seen that the posterior consists of a likelihood part involving data from all persons and a joint prior part $π (ω) π (Σ_{θ_{1} {,θ}_{2}})$ .

Conventional statistical inference is based on the assumption that all elements of a sample are independent. However, especially in educational testing, samples often exhibit a stratified clustered structure, and, consequently, the assumption of independence is violated. Ignoring such a sample structure might lead to underestimated standard errors and might threaten the validity of the inferences drawn from those samples (Rust, 2013). To meet this challenging sample structure, sample weights in conjunction with appropriate resampling methods are usually applied. In order to be in line with these techniques, as a first step, we incorporated sample weights into our estimation procedure (for a description of the procedures regarding inference and model assessment see the next section). Therefore, we weighted the likelihood part that led to a pseudo-likelihood part and an associated pseudo-posterior density (e.g., Goldstein, 2011; Savitsky & Toth, 2016) in Equation 7:

\begin{array}{l} p (ω | y) \propto & [\prod_{p} {(\int P (y_{p} | z_{p}) P (z_{p} | θ_{p 1}, θ_{p 2}, ω) P (θ_{p 1}, θ_{p 2} | Σ_{θ_{1}, θ_{2}}) d z_{p} d θ_{p 1} d θ_{p 2})}^{w_{p}}] \times \\ π (ω) π (Σ_{θ_{1}, θ_{2}}), \end{array}

where w_p is the sample weight of person p. To be consistent with usual likelihood expressions, we set $\sum_{p} w_{p} = N$ .

To apply the Bayesian procedure, the parameter vector is divided into components. Then, for each component, the respective sampling step from the corresponding conditional posterior distribution can be derived by using common Bayesian techniques (Hoff, 2009). A conditional distribution of parameters related to one component is a distribution of parameters conditional on the sampled values for all other parameters (related to all other components). Therefore, the procedure can be described by several estimation steps with one step for each component. To incorporate the sample weights mentioned above into the MCMC steps, let W be the diagonal matrix with entries w_p on the diagonal.

For notational reasons, we define

τ_{p i t} = a_{i 1} θ_{p 1} + a_{i 2} γ_{t} θ_{p 2} + b_{i} + β_{t} .

Step 1

In the first step, the data y are augmented with latent data Z. The $Z_{p i t}$ are independent and normally distributed with mean $τ_{p i t}$ and standard deviation equal to 1, truncated at the left by 0 if $y_{p i t} = 1$ and truncated at the right by 0 otherwise (Albert & Chib, 1993), that is:

Z_{p i t} | Y_{p i t} \sim {\begin{matrix} N (τ_{p i t},1) truncated at the left by 0 if & y_{p i t} = 1 \\ N (τ_{p i t},1) truncated at the right by 0 if & y_{p i t} = 0 \end{matrix} .

For the sake of better readability, we omitted the indices p, i, and t in some cases and used matrix or vector notation. Then, vectors and matrices consist only of components of appropriate combinations of persons, items, and item positions, that is, combinations that are represented in the data.

Step 2

In this step, the item parameters $a_{i 1}$ and b_i are drawn. These parameters are modeled as fixed effects with improper prior distributions (e.g., mean of 0 and infinite variance). Defining $Z_{i}^{*}$ with components $Z_{p i t}^{*} = Z_{p i t} - a_{i 2} γ_{t} θ_{p 2} - β_{t} = a_{i 1} θ_{p 1} + b_{i} + ∊_{p t}$ , the components of parameter vector ${(a_{i 1}, b_{i})}^{T}$ can be considered to be coefficients of the normal regression of $Z_{i}^{*}$ on $X_{i} = (θ_{1}_{i}, 1)$ where $θ_{1}_{i}$ is the vector with appropriate components $θ_{p 1}$ (with person p responded to Item i) and $1$ is the vector with components equal to 1 with appropriate dimension. The conditional distribution can be written as a bivariate normal distribution:

(a_{i 1}, b_{i}) \sim N ({(X_{i}^{T} W_{i} X_{i})}^{- 1} X_{i}^{T} W_{i} Z_{i}^{*}, (X_{i}^{T} W_{i} X_{i})^{- 1}),

where W _i is the corresponding submatrix of W with respective sampling weights on the diagonal.

Step 3

The parameters $a_{i 2}$ are modeled as fixed effects with prior distribution $N (μ_{a_{2}}, σ_{a_{2}}^{2})$ . After implementing the Gibbs sampler, in some cases, the estimation of the $a_{i 2}$ parameters showed some instabilities. To improve the stability of the estimation, we allowed informative prior distributions to be specified for the $a_{i 2}$ parameters. To draw $a_{i 2}$ , we defined $Z_{i}^{*}$ with components $Z_{p i t}^{*} = Z_{p i t} - a_{i 1} θ_{p 1} - b_{i} - β_{t} = a_{i 2} γ_{t} θ_{p 2} + ∊_{p t}$ . The parameter $a_{i 2}$ can be considered to be the coefficient of the normal regression of $Z_{i}^{*}$ on X _i with respective components $γ_{t} θ_{p 2}$ and corresponding matrix W _i .

Incorporating the prior distribution, the conditional distribution can be written as

a_{i 2} \sim N (m_{i} v_{i}, v_{i}),

with $m_{i} = (X_{i}^{T} W_{i} X_{i})^{- 1} X_{i}^{T} W_{i} Z_{i}^{*} + μ_{a_{2}} / σ_{a_{2}}^{2}$ , $v_{i} = {({(X_{i}^{T} W_{i} X_{i})}^{- 1} + 1 / σ_{a_{2}}^{2})}^{- 1}$ . A noninformative prior distribution is specified by setting $μ_{a_{2}} = 0$ and setting $σ_{a_{2}}$ to a very large number.

Step 4

The parameters $β_{t}$ are modeled as fixed effect with noninformative prior distribution. To draw $β = (β_{1}, . . ., β_{T})^{T}$ , we define $Z_{i}^{*}$ with components $Z_{p i t}^{*} = Z_{p i t} - a_{i 1} θ_{p 1} - a_{i 2} γ_{t} θ_{p 2} - b_{i} = β_{t} + ∊_{p i}$ with corresponding matrix W _i . A parameter $β_{t}$ can be viewed as the coefficient of the normal regression of $Z_{i}^{*}$ on $1$ . The conditional distribution can be written as

β_{t} \sim N ({(1^{T} W_{t} 1)}^{- 1} 1^{T} W_{t} Z_{t}^{*}, (1^{T} W_{t} 1)^{- 1}) .

For reasons of identifiability, $β_{1}$ was set to 0.

Step 5

The parameters $γ_{t}$ are modeled as fixed effects with noninformative prior distribution. To draw $γ = (γ_{1}, . . ., γ_{T})^{T}$ , we define $Z_{i}^{*}$ with components $Z_{p i t}^{*} = Z_{p i t} - a_{i 1} θ_{p 1} - (b_{i} + β_{t}) = a_{i 2} γ_{t} θ_{p 2} + ∊_{p i}$ and corresponding matrix W _t. The parameter $γ_{t}$ can be viewed as the coefficient of the normal regression of $Z_{i}^{*}$ on X _t where X _t has appropriate components $g_{i} θ_{p 2}$ . The conditional distribution can be written as

γ_{t} \sim N ({(X_{t}^{T} W_{t} X_{t})}^{- 1} X_{t}^{T} W_{t} Z_{i}^{*}, (X_{t}^{T} W_{t} X_{t})^{- 1}) I (γ_{t} > 0),

where I denotes the indicator function. For reasons of identifiability, $γ_{1}$ and $γ_{2}$ were set to 0 and 1, respectively.

Step 6

To draw $θ_{1}$ and $θ_{2}$ , we consider the vector $θ_{p} = (θ_{p 1}, θ_{p 2})^{T}$ . Equation 9 can be written as

τ_{p i t} = (A_{p i t}^{(1)}, A_{p i t}^{(2)}) θ_{p} + b_{i} + β_{t} + ∊_{i t},

where $A_{p i t}^{(1)} = a_{i 1}$ and $A_{p i t}^{(2)} = a_{i 2} γ_{t}$ . Let A _p be the matrix consisting of the column vectors $A_{p}^{(1)}$ and $A_{p}^{(2)}$ with appropriate components $A_{p i t}^{(1)}$ and $A_{p i t}^{(2)}$ , respectively. Further, let $τ_{p}$ and $β_{p}$ vectors consist of appropriate components $τ_{p i t}$ and $b_{i} + β_{t}$ , respectively. Then, Equation 15 can be written in matrix notation as

τ_{p} = A_{p} θ_{p} + β_{p} .

Let z _p and $∊_{p}$ be vectors consisting of appropriate components $z_{p i t}$ and $∊_{i t}$ , respectively, and $Z_{p}^{*} = z_{p} - β_{p}$ . Then, $Z_{p}^{*} = A_{p} θ_{p} + ∊_{p}$ and the conditional distribution is

θ_{p 1} | Z_{p}^{*} \sim N (D_{p} d_{p}, D_{p}),

where $D_{p}^{- 1} = A_{p}^{T} A_{p} + Σ_{θ_{1} {,θ}_{2}}^{- 1}$ and $d_{p} = A_{p}^{T} Z_{p}^{*}$ .

Step 7

The conjugate prior distribution of $Σ_{θ_{1} {,θ}_{2}}$ is a normal-inverse-Wishart distribution (e.g., Béguin & Glas, 2001) and is chosen to be noninformative for our purposes. Let ${\bar{θ}}_{1} = \sum_{p} w_{p} θ_{p 1} / N$ and ${\bar{θ}}_{2} = \sum_{p} w_{p} θ_{p 2} / N$ denote the weighted mean of $θ_{p 1}$ and $θ_{p 2}$ , respectively. Further, let S be the matrix with entries $S_{11} = \sum_{p} w_{p} {(θ_{p 1} - {\bar{θ}}_{1})}^{2}, S_{22} = \sum_{p} w_{p} {(θ_{p 2} - {\bar{θ}}_{2})}^{2}$ and $S_{12} = S_{21} = \sum_{p} w_{p} (θ_{p 1} - {\bar{θ}}_{1}) (θ_{p 2} - {\bar{θ}}_{2})$ . The conditional distribution is a normal-inverse-Wishart distribution with N degrees of freedom and a scale matrix S, that is

Σ_{θ_{1} {,θ}_{2}} \sim W^{- 1} (S, N) .

This MCMC sampling procedure provides approximations of the posterior distributions of the model parameters. The mode or mean of the corresponding marginal distributions serves as the point estimate for the respective parameters. Statistical inference regarding these estimates when accounting for sampling weights can be considered as a Bayesian version of pseudo-likelihood estimation (Rabe-Hesketh & Skrondal, 2006).

Inference and Model Assessment for Complex Samples

To calculate the likelihood of the model parameters, we follow Merkle et al. (2018) and use the marginal likelihood to derive the model selection criteria because we want to infer to a population beyond the specific sample considered. The correlation $σ_{θ_{1} θ_{2}}$ is the only component of $Σ_{θ_{1} {,θ}_{2}}$ that needs to be estimated, so we define $δ = (ω, σ_{θ_{1} θ_{2}})^{T}$ as the vector consisting of all model parameters to be estimated.

We consider the likelihood part with sampling weights from Equation 8 in which augmented data are integrated out:

L (y | δ) = \prod_{p} {(\int P (y_{p} | θ_{p 1}, θ_{p 2}, ω) P (θ_{p 1}, θ_{p 2} | Σ_{θ_{1}, θ_{2}}) d θ_{p 1} d θ_{p 2})}^{w_{p}},

which can be seen as a Bayesian version of the pseudo-likelihood (e.g., Rabe-Hesketh & Skrondal, 2006). As $θ_{p 1}$ and $θ_{p 2}$ have the same joint prior distribution for all persons p, the likelihood can be simplified to

L (y | δ) = \prod_{p} {(\int P (y_{p} | θ_{1}, θ_{2}, ω) P (θ_{1}, θ_{2} | Σ_{θ_{1}, θ_{2}}) d θ_{1} d θ_{2})}^{w_{p}},

which is a weighted version of the marginal likelihood described in Merkle et al. (2018) with corresponding log pseudo-likelihood

l (y | δ) = \sum_{p} log {(\int P (y_{p} | θ_{1}, θ_{2}, ω) P (θ_{1}, θ_{2} | Σ_{θ_{1}, θ_{2}}) d θ_{1} d θ_{2})}^{w_{p}} = \sum_{p} w_{p} l_{p} (y_{p} | δ),

where $l_{p} (y_{p} | δ)$ is the marginal log likelihood regarding person p.

The weighted analogue of the Fisher information $I (δ)$ is calculated by

I (δ) = - \sum_{p} w_{p} \frac{\partial^{2} l_{p} (y_{p} | δ)}{\partial δ \partial δ^{T}},

and can be approximated numerically using the central difference approximation (e.g., Abramowitz & Stegun, 1972).

Let $\hat{δ}$ be a consistent estimator of $δ$ . To derive standard errors and correlation regarding parameter estimates, we estimate the asymptotic covariance matrix $V (\hat{δ})$ using a set of replicate weights (Kolenikov, 2010). For replication $r = 1, . . ., R$ let ${\hat{δ}}^{(r)}$ be the estimate calculated for the $r th$ replication with corresponding replication weight $w_{p}^{(r)}$ . The asymptotic covariance matrix can be calculated by

V (\hat{δ}) = A \sum_{r = 1}^{R} ({\hat{δ}}^{(r)} - \hat{δ}) ({\hat{δ}}^{(r)} - \hat{δ})^{T},

where A is a scaling constant. For instance, in PISA 2000, the balanced repeated replication was used as the replication method with $R = 80$ replicates and a scaling constant $A = .05$ (Organization for Economic Cooperation and Development [OECD], 2002).

A weighted version of the design effect $Δ$ can then be calculated by (Lumley & Scott, 2015)

Δ = tr (I (\hat{δ}) V (\hat{δ})) / ν,

where $ν$ is the number of estimated parameters.

To assess the models fitted in the empirical example, we used a design-based version of the Akaike information criterion (dAIC; Lumley & Scott, 2015) to account for the complex structure of samples in large-scale assessments. Using the expected a posteriori (EAP) $\hat{δ}$ as a consistent estimator for $δ$ , the dAIC can be calculated by (Lumley & Scott, 2015)

d A I C = - 2 l (y | \hat{δ}) + 2 ν Δ .

Simulation Study

Design and Method

To investigate the performance of the implemented Gibbs sampler, a simulation study was conducted. The model in Equation 1 was used as the data-generating model. As the true population parameters, we chose one fixed set of parameters for 20 items in booklets with four block positions for the data-generating model. Each block consisted of 5 items, and the underlying booklet design (see Table 1) formed a generalized Youden design (e.g., Kiefer, 1975).

The latent factors $θ_{1}$ and $θ_{2}$ were simulated according to a bivariate normal distribution with a mean of 0 and variances of 1. The correlation between $θ_{1}$ and $θ_{2}$ was set to $- .5$ , representing a realistic negatively correlated trait and persistence (e.g., Debeer & Janssen, 2013). Previous studies have shown that in general, the variance of $θ_{1}$ is greater than the variance of $θ_{2}$ (e.g., Hartig & Buchholz, 2012). Consequently, the loadings $a_{i 1}$ and $a_{i 2}$ were chosen within the range of 0.5 to 1.5 and 0.25 to 0.75, respectively. The loading parameters were combined with difficulty parameters b_i within the range of −1.282 to 1.282, corresponding to the $10 th$ and $90 th$ percentiles from the standard normal distribution, respectively. To achieve a high level of diversity of the combinations of the parameter values, the parameters $a_{i 1}, a_{i 2}$ , and b_i were assigned to the single items separately and randomly for each replication. This can be interpreted as a form of item sampling. According to the model, $γ_{1}$ and $β_{1}$ were set to 0 and $γ_{2}$ was set to 1. Further, we set $γ_{3} = 0.3$ , $γ_{4} = 1.2$ , $β_{2} = 0.6$ , $β_{3} = - 0.1$ , and $β_{4} = - 0.3$ .

The simulations were conducted for sample sizes of $500$ , 1,500, and 2,500. This resulted in three settings, each consisting of the same item parameters (with different combinations) and position parameters but with different sample sizes. For each setting, 1,000 replications were conducted.

Noninformative priors were specified for all parameters. We used Gibbs sampling with a burn-in period of 2,500 iterations and a chain of 10,000 iterations. Effective sample sizes were calculated to ensure a sufficiently small MCMC error (Hoff, 2009). In all conditions and for all parameters, the average effective sample size was larger than 200.

The statistical behavior was evaluated in terms of bias, root mean square error (RMSE), and coverage rates. To investigate the frequentist behavior of the parameter estimates, the respective Bayesian credibility interval (BCI) was calculated by using the $5 th$ and $95 th$ percentiles of the corresponding posterior distribution for each replication. Coverage was given if the population parameter $λ$ was included in the BCI. The coverage rate is the relative frequency coverage that occurred across all replications. To save space, we summarize the results for the item-specific parameters only in an aggregate form. A detailed description of the results regarding the item-specific parameters can be found in Appendix Tables A1 through A3 in the online version of the journal.

Results

Bias

In the left part of Table 3, the observed bias for the item parameters is displayed in aggregate form. The estimates for the sample size of 2,500, for the loadings $a_{i 1}$ and $a_{i 2}$ , as well as for the item difficulties b_i were very close to the true values with averaged bias ranging from $- 0.000$ to $0.011$ . For the sample size of $500$ , the $a_{i 1}$ and $a_{i 2}$ parameters showed higher averaged bias of $0.054$ and $0.076$ , respectively, with a tendency to be overestimated; the b_i parameters showed small bias of $0.002$ . The bias of the single-item parameters is shown separately in Table 3. While the bias of the $β$ parameters was small for all sample sizes, the $γ$ parameters tended to be biased toward zero. The parameter $γ_{4}$ showed crucial bias with a value of $- 0.144$ when estimated based on a sample size of $500$ . However, for a sample size of 2,500, the bias was with a value of $- 0.029$ within an acceptable range. The correlation between $θ_{1}$ and $θ_{2}$ seemed to be slightly biased toward zero for small sample sizes. However, all parameters can be viewed as being approximately unbiased for large sample sizes.

Table 3.

Summary of the Results of the Simulation Study

Parameter	True	Bias			RMSE			COV90
		Sample Size (N)			Sample Size (N)			Sample Size (N)
		500	1,500	2,500	500	1,500	2,500	500	1,500	2,500
$a_{i 1}$	1.000	.054	.017	.010	.147	.076	.059	.892	.902	.902
$a_{i 2}$	0.500	.076	.020	.011	.207	.107	.081	.890	.894	.897
b_i	0.000	.002	.001	0.000	.089	.049	.048	.935	.935	.936
$γ_{3}$	0.300	−.030	−.019	−.014	.101	.070	.052	.880	.866	.897
$γ_{4}$	1.200	−.144	−.054	−.029	.194	.095	.071	.693	.824	.860
$β_{2}$	0.600	.020	.005	.005	.051	.027	.022	.906	.935	.912
$β_{3}$	−0.100	−.003	−.001	.001	.044	.025	.019	.907	.904	.897
$β_{4}$	−0.300	−.006	−.004	−.001	.046	.026	.020	.940	.932	.942
$σ_{θ_{1} θ_{2}}$	−0.500	.049	.019	.015	.087	.047	.035	.881	.910	.914

Note. For the data-generating model, one fixed set of parameters was used and sample sizes varied from 500 to 2,500. $N =$ sample size; bias $=$ bias of the expected a posteriori (EAP); RMSE $=$ root mean square error of the EAP. COV90 $=$ coverage rate based on the respective 90% Bayesian credibility interval (BCI; 5th and 95th percentile of the posterior distribution). For the item-specific parameters $a_{i 1}$ , $a_{i 2}$ , and b_i , the respective statistics were averaged over all items, and the respective means are displayed here. The parameters $γ_{1}$ , $γ_{2}$ , and $β_{1}$ were fixed for identification purposes during the estimation procedure and are therefore not listed here.

RMSE

The summaries regarding the RMSE are displayed in the middle part of Table 3. For all parameters, the overall accuracy increased with increasing sample size. For the sample size of $500$ , the estimates of the item parameters showed RMSE within the range of $0.072$ and $0.261$ , where the b_i parameters showed the highest accuracy $[0.072, 0.127]$ and the $a_{i 2}$ parameters the lowest $[0.155, 0.261]$ . The maximum RMSE regarding the item parameters was $0.091$ for the item parameters when the sample size was 2,500. The $γ$ and $β$ parameters showed an RMSE (see Table 3) ranging from $0.044$ to $0.194$ for a sample size of $N = 500$ and ranging from $0.019$ to $0.071$ for a sample size of $N =$ 2,500. The overall accuracy of the estimate of the correlation between $θ_{1}$ and $θ_{2}$ also increased to an acceptable level of $0.035$ for the sample size of $N =$ 2,500 compared to an RMSE of $0.087$ for the sample size of $N = 500$ .

Coverage Rate

To evaluate the accuracy of the standard errors of the estimates, coverage rates were calculated based on the 90% BCI. Again, we briefly summarize the results for the a _i1, a _i2, and b _i parameters. The estimates of all item parameters showed rates very close to or above the nominal level. Additionally, the sample size did not seem to have an influence on the coverage rates (range for sample size of $500$ : $[86.5, 96.0]$ ; range for sample size of 2,500: $[86.6, 92.5]$ ). The coverage rates of the loadings were very close to the nominal level while the rates regarding the b_i parameters were consistently higher than the nominal level for sample sizes of 1,500 and 2,500.

The coverage rates for the remaining parameters are shown in Table 3. All parameters except for $γ_{4}$ showed coverage rates very close to or above the nominal level, respectively. The parameter $γ_{4}$ showed a poor coverage rate for sample size of $500$ , while the coverage rate for sample size of 2,500 was slightly below the nominal level.

Some coverage rates showed a tendency to be overestimated. This might be traced back to the fact that we only used noninformative priors in this simulation study. The use of noninformative priors lead to posterior distributions that represent higher variability. Sampling parameters from informative prior, which are used in the Gibbs sampler, can lead to more accurate coverage rates in simulations.

Discussion

The statistical properties of the proposed model regarding bias, RMSE, and coverage rates were acceptable for sample sizes of $N =$ 2,500. For a sample size of $N =$ 1,500 and smaller, some parameter estimates should be handled with care. The precision of some parameters could be lower than sufficient. For instance, the tendency to overestimate the loading parameters for smaller sample sizes might be traced back to the use of priors that only allows positives values. This leads to a positive skewness of the respective posterior distributions. One could conclude that in these situations, the mode of the posterior distribution (MAP) can be a better estimator compared to the EAP. Indeed, in this simulation study for the sample size of 500, the MAP of the loading parameters showed a smaller bias ranging from 0.002 to 0.086 (bias of the EAP ranged from 0.011 to 0.120). However, for the sample size of 2,500, the bias of the two estimators was of comparable magnitudes: bias of MAP ranging from −0.001 to 0.015 and the bias EAP ranging from 0.001 to 0.020.

Altogether, for sample sizes $N \geq$ 2,500, the simulation study showed that the estimation procedure produces acceptable estimates of the parameters. However, for smaller sample sizes, we recommend using reasonable informative priors to obtain more accurate and precise estimates.

Empirical Example: PISA 2000 Reading Data

To illustrate the proposed approach, we applied the model in Equation 1 and restricted versions of it to a subsample of the published data from the PISA 2000 reading assessment (OECD, 2002).

Data

The sample we analyzed comprised 2,853 students from New Zealand. Note that the sample size differs from that reported in OECD (2002) because students to whom Booklets 8 and 9 were administered were not included in our analyses (see below). As booklets were assigned to students randomly, this subsample can still be considered to be a representative sample of 15-year-old students from New Zealand.

In PISA 2000, items were administered in nine different booklets. Because not all booklets contained reading items in the first three block positions, we confined the set of booklets and items considered in our analyses in two ways: We considered (1) the first three block positions instead of all four block positions in the PISA test and (2) the first seven booklets only. The arrangement of these selected booklets (see Table 2) forms a balanced incomplete block design (BIBD; Frey et al., 2009), that is, every block (or cluster) of items appears in as many booklets as every other block, and any pair of blocks appears in as many booklets as any other pair of blocks. Furthermore, this subdesign forms a generalized Youden design (e.g., Kiefer, 1975), consisting of seven different blocks, which means that, beyond the key constraints of the BIBD, every block is assigned to every block position exactly once. This confinement provides empirical data that make it possible to illustrate the identifiability of the model parameters. Originally, 17 items were polytomously scored with three categories (incorrect, partially correct, and full credit). These items were dichotomized by scoring the full credit as correct and partially correct and incorrect as incorrect. This resulted in data from 96 dichotomously scored items.³ All missing responses coded as omitted or not reached were treated as incorrect responses. All other missing responses were treated as missing by design.

Further, to specify the aforementioned prior distribution for the $a_{i 2}$ parameters, we used the results provided by Debeer et al. (2014) for PISA 2009. The authors analyzed the PISA 2009 reading data from 65 countries, including New Zealand, which is the country we are focusing on in the present article. For their analyses, they fitted—among others—a model for linear position effects similar to the model in Equation 3 to the country-specific data. For New Zealand, they presented the standardized estimates for the standard deviation of $θ_{2}$ , denoted by ${\hat{σ}}^{prior} = 0.17$ . It should be noted that the estimates provided by Debeer et al. (2014) were calculated based on the logit metric. Therefore, we transferred the respective estimate to the probit metric using the factor $1 / 1.7$ (e.g., Lord & Novick, 1968). For New Zealand, we specified the prior distribution of all parameters $a_{i 2}$ as a normal distribution with a mean of $0.1$ and a standard deviation of $0.2$ for all models. For a similar sample size and regarding the item loadings, Patz and Junker (1999, p. 163) report for a standard deviation of $0.5$ as a prior specification that the estimates were not overly sensitive. Therefore, we assume our choice to be a weakly informative prior,⁴ which, however, sufficiently stabilized estimation.

Data analysis

To demonstrate the superiority of the proposed model over more restricted ones, we applied the models described in Equation 1 (M1), Equation 3 (M2), Equation 4 (M3), and Equation 5 (M4) to the data. After inspecting the block position-specific percentage of correct responses (see right part of Table 2), we decided to fit a further model (M2b) with the restrictions than the linear one from Equation 3:

Model M2b : γ {= (0, 1, 1)}^{T}, β = (0, β, β) .

Model M2b assumes that the effects of Position 2 are the same as those of Position 3. One has to note that, in PISA 2000, there was a 15-minute break between Block 2 and Block 3. This might be (at least partly) the reason for why the pattern of the position-specific percentage of correct responses differs from a linear form.

For all models, a Gibbs sampler was applied, while for M2, M2b, M3, and M4, the procedure described in the previous section was adapted. The Gibbs sampling had a burn-in period of 25,000 iterations, and the whole chain was 50,000 iterations long. All reported estimators are the respective means of the posterior distribution (EAP).

To assess the convergence of the MCMC algorithm, we inspected the trace plots of every estimated parameter and the respective autocorrelation function with an increasing lag between the successive draws. All models met the respective diagnostic criteria. Further, effective sample sizes and the proportional scale reduction factor (PSR; Gelman & Rubin, 1992) were calculated for each model parameter. To calculate the PSR, the respective chains were cut in half. Only two estimated parameters (of 872 in total) showed effective sample sizes smaller than 400 but with a minimum of 329. And only one estimated parameter showed (with a value of $1.0662$ ) a PSR above the typically recommended threshold of $1.05$ . Altogether, one can assume sufficient convergence for all fitted models.

To be able to compare the standard deviations of $θ_{1}$ and $θ_{2}$ of the models, some normalization restrictions would have to be imposed on the a _i1 and a _i2. The parameters $σ_{θ_{1}}$ and $σ_{θ_{2}}$ could then be free parameters, and the model would remain identifiable. Another way would be to consider the respective mean of the parameters $a_{i 1}$ and $a_{i 2}$ . Of course, these two values are conceptually not exactly the same as the respective standard deviations. The generalized Youden design and the randomized assignment of students to booklets lead to approximately equally weighted contributions of all items to the respective factor, and, therefore, the respective means of $a_{i 1}$ and $a_{i 2}$ can be considered to be good approximations. However, one has to note that the $a_{i 2}$ are weighted by $γ_{t}$ which is except for Position 1 different from 1 in general and the approximation corresponds to Position 1. To provide an impression of the influence of the different modeling approaches of item position effects on the mean of the item difficulties b_i , the mean of b_i is reported for all models. The statistical significance of the parameters was assessed by calculating the respective standard errors using the appropriate resampling method (Kolenikov, 2010; OECD, 2002). We report significance starting from a nominal α level of $.05$ . The models were compared using a dAIC (see previous section), where a smaller dAIC indicates a better fit to the data.

Results

The results of the analysis regarding the distribution of $θ_{1}$ and $θ_{2}$ , the item position effects, and the respective model fit are shown in the first four columns, in the middle three columns, and in the last three columns of Table 4, respectively.

Table 4.

Summary of the Results of the Analyses Based on the PISA 2000 Reading Data of New Zealand, Confined to the First Three Blocks of Booklets 1–7

Model	$M (a_{i 1})$	$M (a_{i 2})$	$M (b_{i})$	$σ_{θ_{1} θ_{2}}$	$γ_{3}$	$β_{2}$	$β_{3}$	dAIC	$ν$	$Δ$
M1	.743*	.369*	.668*	.250*	0.038	−.162*	−.107*	114,953.8	292	2.208
M2	.744*	.150*	.601*	.129*	2.000^†	−.056*	−.112*	116,140.0	290	3.099
M2b	.732*	.215*	.644*	.205*	1.000^†	−.146*	−.146*	114,958.4	290	1.588
M3	.702*^‡	.346*^‡	.623*	.133*	0.016	−.190*	.105*	115,693.4	102	2.499
M4	.765*		.651*			−.215*	−.109*	115,207.3	194	2.215

Note. $M (a_{i 1}) =$ mean of the loadings parameters $a_{i 1}$ ; $M (a_{i 2}) =$ mean of the loading parameters $a_{i 2}$ ; $M (b_{i}) =$ mean of the item difficulties b_i ; $σ_{θ_{1} θ_{2}} =$ correlation between the ability $θ_{1}$ and the persistence $θ_{2}$ ; dAIC $=$ design-based version of the Akaike information criterion; $ν =$ number of parameters to be estimated; $Δ =$ estimated design effect. $γ_{3}$ is the position-specific parameter of Position 3. $β_{2}$ and $β_{3}$ are the position-specific average changes in difficulty of Position 2 and 3, respectively. $β_{1}$ was set to 0 for all models. $γ_{1}$ was set to 0 for all models except for Model M4. $γ_{2}$ was set to 1 for all models except for Model M4. Values marked with an asterisk are significantly different from zero on a 5% α level. Values marked with a dagger indicate that the respective parameter is fixed due to model restrictions. Values marked with a double dagger indicate that the respective value is the estimated standard deviation of $θ_{1}$ or $θ_{2}$ , respectively, in Model M3.

Considering the mean of the loading parameters $M (a_{i 1})$ , the specification of restrictions did not seem to have any impact on the ability variance for each country. In contrast, the mean of the loadings on $θ_{2}$ seemed to strongly depend on the model restrictions. This is not surprising because the restrictions regarding parameter $γ_{3}$ have a direct influence on the loading structure on $θ_{2}$ related to Position 3. The highest variance was observed for the less restrictive model, M1, followed by M2b. The lowest modeled variance was observed for model M2, in which linear position effects were assumed.

A positive correlation between $θ_{1}$ and $θ_{2}$ here means that students with a higher ability level show lower declines in performance over the course of the test and students with a lower ability level show higher declines. The correlation ranged from $0.129$ for M2 to $0.250$ for M1, where all correlations were significantly different from zero.

To compare the mean changes in difficulty from one position to the other, we used $M (a_{i 1})$ as a proxy of the standard deviation of $θ_{1}$ and approximated standardized effects by dividing the respective estimate by $M (a_{i 1})$ . The standardized mean changes of difficulty from Position 1 to Position 2 ranged from $- 0.075$ $S D$ for M2 to $- 0.281$ $S D$ for M4. Considering M1, the parameter value can be interpreted to show that, on average, there was a decline in performance from Position 1 to Position 2 as high as one would expect between two students whose SD in ability is about two fifths apart.

The standardized mean change in difficulty from Position 1 to Position 3 ranged from $- 0.144$ for M1 to $- 0.199$ $S D$ for M2b.

The estimate of $γ_{3}$ for M1 was not significantly different from zero. The very small value of $γ_{3}$ indicates that the modeled persistence had a very weak (if any) influence on the performance in Block Position 3. In contrast, since $γ_{2}$ is fixed to 1, one can assume that there were systematic individual differences in position-specific performances regarding Position 2, only, which is reflected in $M (a_{i 2}) = 0.369$ . Regarding Position 3, it seems that the persistence did not have any influence on the position-specific performance. Taking the 15-minute break into account, this might be interpreted as a kind of resetting persistence in that it functions in Position 3 as it does in Position 1. This represents a contrast to the findings regarding the mean change in difficulty from Position 1 to Position 3, where a significant and substantial increase in difficulty was observed for M1.

Inspecting the dAIC, M1 is preferable even though the difference from M2b was not very high (see Table 4). This small difference in the model-fit statistics of these two models is at least partly due to the much smaller design effect of M2b. Another noticeable finding is that Model M2 shows the worst fit even if its structure regarding item loadings is more flexible than the structure of Model M3. Also the unidimensional Model M4 shows a better fit than Model M2.

Discussion

Altogether, we can summarize that there were substantial item position effects. Further, the assumption that there are linear decreases in the mean changes in performance during a test did not hold in general. This might be traced back, in our case, to the relatively long break of 15 minutes between Blocks 2 and 3. Finally, there are indications of a factor on the person side that can be interpreted as persistence. In general, the degree to which the person’s persistence influences their response behavior cannot be modeled linearly.

Discussion

We introduced a flexible two-dimensional model for modeling nonlinear position effects on the item side as well as on the person side. Further, the proposed model allows for a flexible loading structure of the items on the two modeled dimensions. A fully Bayesian estimation procedure of the model was proposed, and its performance was demonstrated on simulated data. Further, we applied our approach to reading data collected in PISA 2000 and showed the additional value of the extended flexibility compared to more restrictive models.

The algorithm was implemented in R (R Core Team, 2018). Of course, the usage of other software such as WinBUGS (Lunn et al., 2000) or JAGS (Plummer, 2017) and Stan (Carpenter et al., 2017) was also considered. Handling complex survey data in MCMC estimation as a general-purpose Bayesian software requires the function of weighing individual posterior contributions (Savitsky & Toth, 2016). This seems to be possible in Stan (https://groups.google.com/forum/#!msg/stan-dev/5pJdH72hoM8/GLW_mTeaObA), but it is not easily applicable in WinBUGS or JAGS.

In the following, we discuss the limitations of the proposed model and suggest further developments. Our proposed IRT model for position effects is limited to dichotomous item responses. However, an extension allowing for ordinal response data is possible. For an Item i with C ordered response categories, in the data augmentation step $C - 1$ threshold parameters have to be introduced that divide the continuum in C intervals. For an observed category k, the auxiliary variable $Z_{i k}$ is then drawn from a normal distribution truncated to the interval representing category k. To achieve better convergence properties in the case of ordinal responses, the substitution of some Gibbs sampling steps with Metropolis–Hastings sampling steps has been proposed (e.g., Fox, 2010, p. 83).

Further, the model could be extended by introducing the residual item-specific position effects of item intercepts as random effects. This strategy can be extended by including random effects for item-specific position effects for item loadings (Le, 2009).

The central assumption underlying the proposed model was a bivariate normal distribution for ability and persistence. Consequently, persistence takes negative values as well as positive values that decrease and increase the success probability of a correct item response, respectively. In contrast, if persistence is intended to represent performance decline only, restricting persistence to negative values could be reasonable. Hence, the persistence component in the bivariate distribution could be assumed to follow a log-normal distribution multiplied by −1, which only admits negative values. This model is related to a performance-decline item response model (Jin & Wang, 2014). Although the model is formulated for investigating declines in each item in a fixed-order test, the model could also be applied to the situation of position effects in this article, and the decline of a student can occur at the beginning of the test or in the second or third positions.

Supplemental Material

Supplemental Material, 01_JEBS-2019-MAN-0038.R2_Appendix - A Bayesian Item Response Model for Examining Item Position Effects in Complex Survey Data

Supplemental Material, 01_JEBS-2019-MAN-0038.R2_Appendix for A Bayesian Item Response Model for Examining Item Position Effects in Complex Survey Data by Matthias Trendtel and Alexander Robitzsch in Journal of Educational and Behavioral Statistics

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Matthias Trendtel

Notes

References

Abramowitz

Stegun

I. A.

(1972). Handbook of mathematical functions: With formulas, graphs, and mathematical tables (Vol. 55). Dover Publications.

Albert

J. H.

Chib

(1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association, 88(422), 669–679.

Alexandrowicz

Matschinger

(2008). Estimation of item location effects by means of the generalized logistic regression model: A simulation study and an application. Psychology Science Quarterly, 50(1), 64–74.

Béguin

A. A.

Glas

C. A.

(2001). MCMC estimation and some model-fit analysis of multidimensional IRT models. Psychometrika, 66(4), 541–561.

Carpenter

Gelman

Hoffman

M. D.

Lee

Goodrich

Betancourt

Brubaker

M. A.

Guo

Riddell

(2017). Stan: A probabilistic programming language. Journal of Statistical Software, 76(1), 1–32.

Debeer

Buchholz

Hartig

Janssen

(2014). Student, school, and country differences in sustained test-taking effort in the 2009 PISA reading assessment. Journal of Educational and Behavioral Statistics, 39(6), 502–523.

Debeer

Janssen

(2013). Modeling item-position effects within an IRT framework. Journal of Educational Measurement, 50(2), 164–185.

De Boeck

Wilson

(2004). A framework for item response models. In De Boeck

Wilson

(Eds.), Explanatory item response models (pp. 3–41). Springer.

Fox

J.-P.

(2010). Bayesian item response modeling: Theory and applications. Springer.

10.

Fox

J.-P.

Verhagen

A. J.

(2010). Random item effects modeling for cross-national survey data. In Davidov

Schmidt

Billiet

(Eds.), Cross-cultural analysis: Methods and applications (pp. 467–488). Routledge Academic.

11.

Frey

Hartig

Rupp

A. A.

(2009). An NCME instructional module on booklet designs in large-scale assessments of student achievement: Theory and practice. Educational Measurement: Issues and Practice, 28(3), 39–53.

12.

Gelfand

A. E.

Smith

A. F.

(1990). Sampling-based approaches to calculating marginal densities. Journal of the American Statistical Association, 85(410), 398–409.

13.

Gelman

Carlin

J. B.

Stern

H. S.

Dunson

D. B.

Vehtari

Rubin

D. B.

(2013). Bayesian data analysis (3rd ed.). CRC Press.

14.

Gelman

Rubin

D. B.

(1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7, 457–472.

15.

Goldstein

(2011). Multilevel statistical models (4th ed.). John Wiley.

16.

Hartig

Buchholz

(2012). A multilevel item response model for item position effects and individual persistence. Psychological Test and Assessment Modeling, 54(4), 418–431.

17.

Hoff

P. D.

(2009). A first course in Bayesian statistical methods. Springer.

18.

Hohensinn

Kubinger

K. D.

Reif

Holocher-Ertl

Khorramdel

Frebort

(2008). Examining item-position effects in large-scale assessment using the linear logistic test model. Psychology Science Quarterly, 50(3), 391.

19.

Jin

K.-Y.

Wang

W.-C.

(2014). Item response theory models for performance decline during testing. Journal of Educational Measurement, 51(2), 178–200.

20.

Kiefer

(1975). Balanced block designs and generalized Youden designs, I. Construction (patchwork). The Annals of Statistics, 3(1), 109–118.

21.

Kolenikov

(2010). Resampling variance estimation for complex survey data. Stata Journal, 10(2), 165–199.

22.

Kubinger

K. D.

(2008). On the revival of the Rasch model-based LLTM: From constructing tests using item generating rules to measuring item administration effects. Psychology Science Quarterly, 50(3), 311–327.

23.

L. T.

(2009). Effects on item positions on their difficulty and discrimination: A study in PISA science data across test language and countries. In Shigemasu

Okada

Imaizumi

Hoshino

(Eds.), New trends in psychometrics (pp. 207–214). Universal.

24.

Levy

Mislevy

R. J.

(2016). Bayesian psychometric modeling. CRC Press.

25.

Lord

F. M.

Novick

M. R.

(1968). Statistical theories of mental test scores. Addison-Wesley.

26.

Lumley

Scott

(2015). AIC and BIC for modeling with complex survey data. Journal of Survey Statistics and Methodology, 3(1), 1–18.

27.

Lunn

D. J.

Thomas

Best

Spiegelhalter

(2000). WinBUGS-a Bayesian modelling framework: Concepts, structure, and extensibility. Statistics and Computing, 10(4), 325–337.

28.

Merkle

Furr

Rabe-Hesketh

(2018). Bayesian model assessment: Use of conditional vs marginal likelihoods. ArXiv preprint. https://arxiv.org/abs/1802.04452

29.

Mollenkopf

W. G.

(1950). An experimental study of the effects on item analysis data of changing item placement and test-time limit. Psychometrika, 15, 291–315.

30.

Muthén

(1979). A structural probit model with latent variables. Journal of the American Statistical Association, 74(368), 807–811.

31.

Organization for Economic Cooperation and Development. (2002). PISA 2000 Technical Report. OECD Publishing.

32.

Patz

R. J.

Junker

B. W.

(1999). A straightforward approach to Markov chain Monte Carlo methods for item response models. Journal of Educational and Behavioural Statistics, 24(2), 146–178.

33.

Plummer

(2017). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. Version 4.3.0. http://mcmc-jags.sourceforge.net

34.

Rabe-Hesketh

Skrondal

(2006). Multilevel modelling of complex survey data. Journal of the Royal Statistical Society: Series A (Statistics in Society), 169(4), 805–827.

35.

R Core Team. (2018). R: A language and environment for statistical computing [Computer software manual]. R Core Team. https://www.R-project.org/

36.

Renard

Molenberghs

Geys

(2004). A pairwise likelihood approach to estimation in multilevel probit models. Computational Statistics & Data Analysis, 44(4), 649–667.

37.

Rupp

A. A.

Dey

D. K.

Zumbo

B. D.

(2004). To Bayes or not to Bayes, from whether to when: Applications of Bayesian methodology to modeling. Structural Equation Modeling, 11(3), 424–451.

38.

Rust

(2013). Sampling, weighting, and variance estimation in international large-scale assessments. In Rutkowski

von Davier

Rutkowski

(Eds.), Handbook of international large-scale assessment: Background, technical issues, and methods of data analysis (pp. 117–153). CRC Press.

39.

San Martin

González

(2010). Bayesian identifiability: Contributions to an inconclusive debate. Chilean Journal of Statistics, 1(2), 69–91.

40.

San Martin

Rolin

J.-M.

Castro

L. M.

(2013). Identification of the 1PL model with guessing parameter: Parametric and semi-parametric results. Psychometrika, 78(2), 341–379.

41.

Savitsky

T. D.

Toth

(2016). Bayesian estimation under informative sampling. Electronic Journal of Statistics, 10(1), 1677–1708.

42.

Schweizer

Schreiner

Gold

(2009). The confirmatory investigation of APM items with loadings as a function of the position and easiness of items: A two-dimensional model of APM. Psychology Science Quarterly, 51(1), 47–64.

43.

Trendtel

Robitzsch

(2018). Modeling item position effects with a Bayesian item response model applied to PISA 2009–2015 data. Psychological Test and Assessment Modeling, 60(2), 241–263.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.20 MB

0.00 MB