Identifying Local Dependence With a Score Test Statistic Based on the Bifactor Logistic Model

Abstract

Local dependence (LD) refers to the violation of the local independence assumption of most item response models. Statistics that indicate LD between a pair of items on a test or questionnaire that is being fitted with an item response model can play a useful diagnostic role in applications of item response theory. In this article, a new score test statistic, S_b, based on the bifactor logistic model is described. To compare the performance of S_b with the score test statistic (S_t) based on the threshold shift model, and the LD X² statistic, data were simulated under locally independent, bifactor, and threshold shift conditions. The results summarize the null distributions of all three diagnostic statistics, and their power for various degrees of bifactor and threshold shift LD. Future research directions are discussed, including the straightforward generalization of S_b for polytomous item response models, and the challenges involved in the corresponding generalizations of S_t and LD X².

Keywords

item response theory local dependence bifactor model score test

Item response theory (IRT) provides a collection of latent variable models and statistical procedures (including parameter estimation techniques, model evaluation methods and diagnostics) used for item analysis and test scoring (Thissen & Steinberg, 2009). One basic assumption of IRT models is local independence (LI), which requires that responses to different items are independent conditional on the latent variable of interest,¹ θ. Formally, LI implies that the probability of observing responses x = (x₁, . . . , x_j)′ to all J items given θ is the product of each item’s trace line function T_j(x_j|θ) for response x_j:

π (x | θ) = \prod_{j = 1}^{J} T_{j} (x_{j} | θ) .

Equation 1 formalizes the strong version of LI (McDonald, 1982); it serves as the basis for both estimating item parameters and scoring, where the likelihood function of the overall response pattern for a given θ is usually treated as though the contribution of each item is independent. Although strong LI is important, there is no reliable statistical procedure to test it, because it requires no conditional dependency among item responses in marginal subtables of any order. Consequently, only violations of pairwise independence in certain parametric forms, which are introduced in the next section, are tested in practice.

When pairwise local dependence (LD) is of interest, it is natural to distinguish between the positive and negative valence of LD. Positive LD indicates that positive association remains between the pair after the interdependency attributable to the latent variable(s) is taken into account. One frequently cited example of positive LD involves redundant items, that is, (nearly) the same question asked more than once in a test. Meanwhile, negative LD indicates negative conditional association, which is in general less common and more difficult to explain.

Existing Methods to Identify LD Models

Models

Two parametric models are considered in the current study for the purpose of identifying LD: the bifactor model (Gibbons & Hedeker, 1992) and a threshold shift model (Glas & Suárez Falcón, 2003). Figure 1 illustrates the two models as path diagrams.²

Figure 1.

Path diagrams: left, bifactor model; right, threshold shift model

The bifactor model is a special multidimensional factor analysis model (Gibbons & Hedeker, 1992; for the more general two-tier model, see Cai, 2010). The specification of the item bifactor model implies that an extra latent variable θ₂ produces residual covariances for a subset of items when the data are fitted with a unidimensional model. A bifactor logistic model for dichotomous item pair p and q with responses x_p, x_q = {0, 1} is as follows:

{\begin{matrix} T_{p} (x_{p} | θ_{1}, θ_{2}) = \frac{1}{1 + e x p [{(- 1)}^{x_{p}} (a_{p} θ_{1} + a_{p q} θ_{2} + c_{p})]} \\ T_{q} (x_{q} | θ_{1}, θ_{2}) = \frac{1}{1 + e x p [{(- 1)}^{x_{q}} (a_{q} θ_{1} \pm a_{p q} θ_{2} + c_{q})]} \end{matrix},

in which a_p and a_q are the slope parameters on the primary latent dimension θ₁, c_p and c_q are intercept parameters, and a_pq is the slope parameter on the secondary dimension θ₂, representing LD.

To identify the model when the LD subset comprises only two items, some constraint must be imposed on the two secondary dimension slope parameters; a convenient restriction imposed in Equation 2 is equality of the absolute value of bifactor slopes. The positivity/negativity of the a_pqθ₂ term in the second formula is determined by positive/negative LD to be modeled: A positive sign implies positive LD, and a negative sign implies negative LD. The locally independent two-parameter logistic (2PL) model is a special case of bifactor model obtained by setting a_pq = 0 in Equation 2.

The secondary factor loading can be understood analogously as representing error covariance in the linear-normal structural equation model: It can be proved (see Appendix A) that the LD bifactor model is equivalent to an error-covariance model if the indicator variables are normally distributed and the link function is identity. Another parameterization of a restricted bifactor model was given by Bradlow, Wainer, and Wang (1999), in which bifactor slopes are fixed to be equal to their corresponding primary slope, and the variance of the secondary factor is estimated.

The development of a threshold shift model may be traced to work by Kelderman (1984) and Jannarone (1986) on the Rasch model. It also appeared in work by Hoskens and De Boeck (1997), where it was called an “ordered-constant” interaction model, among four types of item response models using additional pairwise interaction parameters. Glas and Suárez Falcón (2003) derived a score test for this model’s interaction term (see also van der Linden & Glas, 2010).

In the 2PL case, for example, the trace line function of the second item in a locally dependent item pair p and q (p < q) may be written as follows:

T_{q} (x_{q} | θ_{1}) = \frac{1}{1 + e x p [{(- 1)}^{x_{q}} (a_{q} θ_{1} + c_{q} + δ_{p q} x_{p})]},

where δ_pq can be considered a threshold shift for the second response x_q when the first item response x_p is positive. In other words, endorsing the first item makes the second item less or more likely to be endorsed depending on δ_pq being negative or positive. δ_pq is also an ANOVA-like interaction term when the log odds of the pairwise response pattern is considered (Hoskens & De Boeck, 1997). Notice that if δ_pq = 0, the threshold shift model reduces to the locally independent model.

In Appendix A, it is proved that for the linear-normal structural model, the threshold shift model with the shift parameter taken as a regression coefficient is equivalent to the error-covariance model, and thus to the bifactor model. However, in IRT where the indicator variables are Bernoulli distributed and the link function is logistic, the bifactor model and the threshold shift model are in general not exactly equivalent (with an exception of a degenerate 2PL case).

Apart from these models, there are other limited or full information parametric models for LD. For example, both LISREL and Mplus can include error covariances in a factor analytic model based on the polychoric correlation matrix of ordered categorical indicators (Muthén, 1978; see Wirth and Edwards, 2007, for a review). Braeken, Tuerlinckx, and De Boeck (2007) proposed to estimate the dependency parameter of the bivariate logistic distribution constructed with copula functions. Recently, Kim, De Ayala, Ferdous, and Nering (2011) compared the performance of 10 existing indices of LD under various simulated conditions. However, Kim et al.’s set of statistics did not overlap with those in the current study,³ so no attempt is made to compare the results of this study with theirs. Interested readers are referred to their article for more information.

Statistical Procedures

Score test: Let η be the vector of parameters for some item response model; also define the partition of η = (η₀, η₁), which leaves out a single parameter η₁. Consider the following hypothesis testing problem with nuisance parameters η₀:

H_{0} : η_{1} = 0 vs . H_{1} : η_{1} \neq 0,

which is a special type of nested-model comparison in which the alternative model is reduced to the null after fixing one parameter to zero. The connection of this formal statistical test with the bifactor LD model can be established by setting η₁ = a_pq and η₀ = (a, c) in which a includes all primary dimension slopes and c all intercepts; then it becomes a test for the violation of pairwise LI triggered by a secondary latent variable specific to item p and q. Similarly, identifying η₁ = δ_pq and η₀ in the same fashion gives a parallel test based on the threshold shift model.

The score test (Lehmann, 1999, also known as Lagrange multiplier test) statistic for testing Problem 4 is defined as

S^{*} = \frac{1}{N} g ({\tilde{η}}_{0}, 0)' ℐ {({\tilde{η}}_{0}, 0)}^{- 1} g ({\tilde{η}}_{0}, 0),

in which g(η₀, η₁) is the score function (i.e., gradient vector of log-likelihood with respect to η) for the alternative model and ℐ(η₀, η₁) is the Fisher information matrix for one observation. Both are evaluated at $({\tilde{η}}_{0}, 0)$ , which is the constrained maximum likelihood estimate (MLE) under the null hypothesis. For the aforementioned tests of LD based on the bifactor model and the threshold shift model, the constrained MLE is the collection of parameter estimates obtained with locally independent model.

To compute the LD score test statistics, partial derivatives of log-likelihood under each LD model will be derived in the next section. When LI is true (i.e., under H₀), it can be proved (e.g., Lehmann, 1999) that

S * \overset{d}{\to} X_{1}^{2} .

One advantage of the score test over other asymptotic tests, such as the Wald test and the likelihood ratio test, is that the computation of the test statistic does not rely on the MLE under the alternative, which would require model refitting. In LD detection, the hypothesis is usually tested for all pairs of items; the process of obtaining MLE for each pairwise LD model would be very time-consuming for long tests. Therefore, the score test is more efficient computationally.

Residual measures

Summary statistics of bivariate marginal residuals also provide evidence about the assumption of LI. One popular choice in this category is the LD X² (Chen & Thissen, 1997) statistic. Taking the dichotomous case as an example, for each pair of items the 2× 2 tables as shown in Table 1 can be constructed for observed and expected probabilities; ps denote the observed cell proportions and πs the expected ones.

Table 1.

Bivariate Marginal Tables for Observed Proportions $p_{x_{p} x_{q}}$ and Modeled Proportions $π_{x_{p} x_{q}}$ for Items p and q (x_p, x_q = {0, 1})

		Item q
		0	1
Item p	0	p ₀₀	p ₀₁
	1	p ₁₀	p ₁₁

		Item q
		0	1

Item p	0	π₀₀	π₀₁
	1	π₁₀	π₁₁

Expected cell proportions can be calculated under the locally independent item response model:

π_{x_{p} x_{q}} = \int_{- \infty}^{\infty} T_{p} (x_{p} | θ) T_{q} | (x_{q} | θ) ϕ (θ) d θ,

in which ϕ(θ) is the probability density function, assumed to be N(0, 1). The LD X² statistic is the Pearson’s X², a quadratic form of the discrepancy of these observed and expected marginal probabilities (i.e., residuals):

X^{2} = N \sum_{x_{p} = 0}^{1} \sum_{x_{q} = 0}^{1} \frac{{(p_{x_{p} x_{q}} - π_{x_{p} x_{q}})}^{2}}{π_{x_{p} x_{q}}} .

The asymptotic null distribution for X² is difficult to compute. Chen and Thissen suggested the use of $χ_{1}^{2}$ as an approximation based on their simulation results. They argued that the degree of freedom is one for the test of independence in the 2 × 2 marginal table, whereas estimating slope parameters from the relationships among items could be regarded as imposing fractional loss of the one degree of freedom.

Other residual diagnostics have been proposed. For example, Yen’s (1984) Q₃ is defined as the sample Pearson correlation between paired residuals. However, it has the same problem as LD X², that the asymptotic distribution of the residual correlation is not clear under the null hypothesis.⁴ Another residual diagnostic, the DIMTEST T statistic, was proposed by Stout (1987) for testing unidimensionality. As compared with Yen’s Q₃, DIMTEST is a nonparametric procedure testing whether the average residual covariance in a given item subset is significantly larger than zero. Since all the procedures investigated in the current study are parametric, further details of the DIMTEST statistic are not included here.

The research described here investigates the feasibility and usefulness of a score test based on the bifactor model. Whether the performance of the LD score test depends on the choice of alternative models is also examined by comparing this statistic with the one based on the threshold shift alternative (Glas & Suárez Falcón, 2003) in various simulated conditions. LD X² is also involved in the simulation study, which is intended to verify and extend the results of Chen and Thissen (1997). Because LD X² statistics are reported by the computer software IRTPRO (Cai, Thissen, & du Toit, 2011) as item-pair level diagnostics, their performance would be of interest to many substantive users.

Theory of a Score Test

Derivatives in General Form

Let the response for each item be dichotomous (i.e., x_ij is 0 or 1). Bock (2010) provided, for the multidimensional dichotomous IRT model, a general form of the first derivatives of the marginal log-likelihood function log L(η, x) with respect to an item parameter ηs out of η:

\begin{matrix} \frac{\partial l o g L (η, x)}{\partial η_{s}} = \sum_{i = 1}^{N} \frac{1}{π (x_{i})} \frac{\partial π (x_{i})}{\partial η_{s}} \\ = \sum_{i = 1}^{N} \frac{1}{π (x_{i})} \int_{θ} \frac{\partial π (x_{i} | θ)}{\partial η_{s}} ϕ (θ) d θ \\ = \sum_{i = 1}^{N} \frac{1}{π (x_{i})} \int_{θ} [\sum_{j = 1}^{J} \frac{x_{i j}}{T_{j} (x_{i j} | θ)} \frac{\partial T_{j} (x_{i j} | θ)}{\partial η_{s}}] π (x_{i} | θ) ϕ (θ) d θ, \end{matrix}

where $π (x_{i} | θ) = Π_{j = 1}^{J} T_{j} (x_{ij} | θ)$ , $π (x_{i}) = \int_{θ} π (x_{i} | θ) ϕ (θ) d θ$ , x_ij is subject i’s response to the jth item, and ϕ(θ) is the multivariate normal probability density function of latent variables θ with prespecified parameter values. Equation 9 gives the expression for elements in the gradient vector g.⁵ The computation of Fisher information involves summing over all possible response patterns

ℐ_{s t} = \sum_{{all possible x_{i}' s}} \frac{1}{π (x_{i})} \frac{\partial π (x_{i})}{\partial η_{s}} \frac{\partial π (x_{i})}{\partial η_{t}} .

In practice, however, this computation may become impossible for long tests because the number of operations to be performed increases as an exponential function of test length. One computationally feasible consistent estimate is the cross-product (XPD) approximation (Bock & Lieberman, 1970; Kendall & Stuart, 1961):

{({\hat{ℐ}}_{XPD})}_{s t} = \sum_{{observed x_{i}' s}} \frac{p (x_{i})}{π {(x_{i})}^{2}} \frac{\partial π (x_{i})}{\partial η_{s}} \frac{\partial π (x_{i})}{\partial η_{t}},

in which the sum is taken over all observed response patterns. Because XPD approximation is a consistent estimate of Fisher information, by Slutsky’s theorem (e.g., Lehmann, 1999) and Equation 6, an asymptotically equivalent version of score test statistic can be defined as

S = \frac{1}{N} g' {\hat{ℐ}}_{XPD}^{- 1} g \overset{d}{\to} χ_{1}^{2} .

In Equations 9 and 11, only the first derivatives of trace line with respect to the parameters of the alternative model are involved when calculating the XPD estimates of information.

Derivatives for Multidimensional Logistic Model

Considering the bifactor LD model (and then 2PL model) as a special case of the multidimensional logistic model, the derivatives of the trace line function of the latter are needed. The trace line function has the following form:

T_{j} (x_{j} | θ) = \frac{1}{1 + e x p [{(- 1)}^{x_{j}} (a_{j}' θ + c_{j})]},

in which a_j is a vector of slope parameters for the vector of latent variables θ. Equation 2 is a special case of Equation 13. The first derivative with respect to slope and intercept parameters are given by Equations 14 and 15, respectively:

\frac{\partial T_{j} (x_{j} | θ)}{\partial a_{j}} = θ T_{j} (x_{j} | θ) [x_{j} - T_{j} (1 | θ)],

\frac{\partial T_{j} (x_{j} | θ)}{\partial c_{j}} = T_{j} (x_{j} | θ) [x_{j} - T_{j} (1 | θ)] .

Derivatives for Threshold Shift Model

In the threshold shift LD model, the trace line for the second item in an item pair is given by Equation 3, and the corresponding first derivative with respect to δ_pq is

\frac{\partial T_{q} (x_{q} | θ_{1}, x_{p})}{\partial δ_{p q}} = x_{p} T_{q} (x_{q} | θ_{1}) [x_{q} - T_{q} (1 | θ_{1})] .

The derivatives of item qs trace line with respect to its slope and intercept are in the same form of Equations 14 and 15, except that the multidimensional logistic trace line function is replaced by Equation 3. Partial derivatives of trace lines are inserted in Equations 9 and 11 to yield the gradient and XPD information, which are further used in Equation 12 to compute the score test statistic, S_t, based on the threshold shift model (Glas & Suárez Falcón, 2003).

A Special Issue With the Score Test for Bifactor Slopes

It was mentioned earlier that the bifactor LD model reduces to the locally independent 2PL model when the secondary slope is zero; however, the test statistic cannot be computed at exactly a_pq = 0. A detailed proof is provided in Appendix B that the score function is zero when a_pq is set to zero, which indicates a local maximum or saddle point. Figure 2 provides a visualization of this point by plotting the profile log-likelihood⁶ as a function of a_pq. The data used to compute this profile log-likelihood were generated from a bifactor positive LD model. When the true model is used (i.e., positive LD model, solid line in Figure 2), a_pq = 0 is a saddle point; when the wrong model is used (i.e., negative LD model, dashed line of Figure 2), it is a local maximum.

Figure 2.

Profile log-likelihood for bifactor model as a function of secondary slope

As a result, the null hypothesis a_pq = ε is used, with ε some small value (0.0001 in this research), and the score test statistic based on the bifactor model is computed as

S_{b} = \frac{1}{N} g (\tilde{η}, ε)' {\hat{ℐ}}_{XPD} {(\tilde{η}, ε)}^{- 1} g (\tilde{η}, ε) .

Due to the symmetry (shown as the dotted line in Figure 2) of the likelihood function with respect to a_pq = 0, it is only necessary to consider a_pq > 0 when computing the score statistics for both positive and negative LD models. After the two statistics associated with these two models are obtained, the sign of score function at a_pq = ε is checked to determine which model to keep: The model with a positive score function value is retained; if both models produce negative values, then the LI model is taken as approximately correct.

Method

Simulation Design

To investigate the performance of LD statistics, dichotomous item response data were generated using R (R Development Core Team, 2010) from the locally independent 2PL model (null conditions), the threshold shift model (TS conditions), and the bifactor model (BF conditions), respectively.

For the null case, data were generated with 10, 25, and 50 items for the sample sizes of 200, 500, and 1,000 (see Table 2). Item parameters were sampled from the same distributions as Chen and Thissen (1997) used, which were believed to resemble the empirical distributions seen in practice. Specifically, the slopes were drawn from log-normal distribution loga ~ N(0, 0.5), and the thresholds (b = −c/a) were drawn from N(0, 1:5). The data generation procedure was replicated 1,000 times for each condition; LD statistics were extracted only for the first item pair of each replication.

Table 2.

Simulation Conditions

Condition	No. of items	Sample size	Strength of LD
Null	10	200, 500, 1,000	—
	25	200, 500, 1,000	—
	50	200, 500, 1,000	—
TS	10	200, 500, 1,000	δ_pq = 0.16, 0.34, 0.66
	25	200, 500, 1,000	δ_pq = 0.16, 0.34, 0.66
	50	200, 500, 1,000	δ_pq = 0.16, 0.34, 0.66
BF	10	200, 500, 1,000	δ_pq = 0.54, 0.96, 1.67
	25	200, 500, 1,000	δ_pq = 0.54, 0.96, 1.67
	50	200, 500, 1,000	δ_pq = 0.54, 0.96, 1.67

Note: TS = threshold shift model; BF = bifactor model. Design: No. of Items × Sample Sizes (× Strength).

To evaluate the power of statistics when LD is present, data were generated from both the threshold shift model and the bifactor model. Apart from sample size and number of items, there was also variation of the strength of LD, as shown in Table 2. The association parameters were generated on the factor loading scale, and then transformed to produce a_pq and δ_pq parameters. Specifically, secondary factor loadings λ_pq were generated from N(μ_λ, 0.01), where μ_λ = 0.3, 0.5, 0.7, indicating low, moderate, and strong LD. To obtain valid values for loadings and control the dispersion of item parameters, all these sampling distributions were truncated to the interval μ_λ − 0,2, μ_λ + 0,2.

The bifactor slopes can be transformed from loadings using the following exact formula (see Wirth & Edwards, 2007):

a_{p q} = \frac{1.702 λ_{p q}}{\sqrt{1 - λ_{p q}^{2}}} .

The transformation between δ_pq and θ_pq is not exact because no apparently equivalent factor analysis model is defined for the threshold shift model. One approximation is to use the result proved (in Appendix A) with continuous indicators:

λ_{p q} = \sqrt{δ_{p q} (1 - λ_{p}^{2})},

where λ_p is the primary slope for the first item in each LD pair.

Although negative LD does not have the same substantive interpretation as positive LD, neither of them is mathematically more important. Thus, under each simulated condition, 500 negative LD pairs were generated out of the total 1,000 replications by switching the sign of a_pq and δ_pq to negative in the trace line of item q.

A modified version of the computational engine of the software IRTPRO was used for estimating item parameters and computing LD statistics.

Evaluation of LD Statistics

The distribution of the bifactor score test statistic S_b, the threshold shift score test statistic S_t, and LD X² under the null hypothesis can be obtained by pooling over all replications within each cell. For the locally independent conditions, the empirical distributions are compared to the $χ_{1}^{2}$ distribution⁷ by means of quantiles. This shows the Type I error rates and empirical p values of the LD statistics.

The receiver operating characteristic (ROC) curves of all three statistics for all conditions are presented to provide information about the power of these tests. The horizontal axis of each ROC curve represents the false positive rate, or the alpha level in the setting of hypothesis testing; the vertical axis represents the true positive rate, which is the power of the statistics.

Results

Null Conditions

Several important quantiles (i.e., 0.25, 0.5, 0.75, 0.9, 0.95, and 0.99⁸) of the empirical distributions are tabulated for each condition of the simulation in Table 3; the corresponding quantiles of the $χ_{1}^{2}$ distribution are included in the footnote of the table.

Table 3.

Empirical Quantiles of the Null Distribution of the Three Statistics

	25%	50%	75%	90%	95%	99%
10 items
N = 200
S_b	0.13	0.49	1.48	3.28	4.95	10.63
S_t	0.13	0.50	1.45	3.38	5.16	29.86
X²	0.06	0.25	0.71	1.46	2.21	4.10
N = 500
S_b	0.11	0.46	1.29	2.56	3.73	8.96
S_t	0.10	0.45	1.31	2.67	3.88	8.28
X²	0.05	0.26	0.71	1.52	2.26	4.69
N = 1,000
S_b	0.11	0.46	1.34	2.85	4.03	7.46
S_t	0.11	0.48	1.33	2.94	4.16	7.35
X²	0.06	0.27	0.75	1.77	2.39	4.03
25 items
N = 200
S_b	0.14	0.65	2.16	4.32	6.12	19.84
S_t	0.15	0.65	2.10	4.28	6.52	45.26
X²	0.08	0.32	1.01	2.34	2.89	4.08
N = 500
S_b	0.10	0.49	1.61	3.46	5.28	11.27
S_t	0.11	0.50	1.68	3.43	5.16	11.44
X²	0.07	0.31	0.96	2.08	2.83	4.78
N = 1,000
S_b	0.09	0.48	1.44	2.82	4.37	6.79
S_t	0.09	0.48	1.37	2.87	4.46	7.11
X²	0.07	0.32	0.93	1.83	2.67	4.35
50 items
N = 200
S_b	0.19	0.84	2.88	6.53	10.47	325.45
S_t	0.19	0.83	2.83	7.12	10.84	240.26
X²	0.11	0.39	1.07	2.29	3.35	42.08
N = 500
S_b	0.11	0.54	1.70	3.46	5.26	10.56
S_t	0.12	0.58	1.70	3.48	5.38	11.85
X²	0.11	0.39	1.06	1.99	2.89	5.58
N = 1,000
S_b	0.11	0.53	1.58	3.34	4.80	8.79
S_t	0.13	0.52	1.58	3.34	4.87	10.57
X²	0.14	0.40	1.05	2.03	3.12	5.61

Note: LD = local dependence. The corresponding quantiles for $χ_{1}^{2}$ is 0.10, 0.45, 1.32, 2.71, 3.84, and 6.63.

A general trend is that both S_b and S_t tend to be liberal, whereas LD X² is conservative if the $χ_{1}^{2}$ cutoffs are used. Both increasing test length and decreasing sample size result in larger statistics, which makes the distribution of LD X² closer to $χ_{1}^{2}$ , but which further exacerbates the liberality of both score test statistics. The empirical quantiles of S_b and S_t are nearly identical; that can be seen more clearly in empirical cumulative density function (CDF) plots in Figure 3.

Figure 3.

Empirical CDF plots under the null conditions

From Figure 3, it can be concluded that $χ_{1}^{2}$ serves as a good approximation of the null distribution for the score statistics S_b and S_t in a short test (i.e., 10 items) with medium (i.e., N = 500) or large (i.e., N = 1, 000) samples, as well as a medium length test (i.e., 25 items) with large samples. Otherwise, both score statistics are liberal. In all the conditions of this study, the empirical distributions of LD X² are always stochastically smaller than $χ_{1}^{2}$ , which indicates conservativeness.

Threshold Shift Conditions

Figures 4 to 6 show the ROC curves for the three statistics with data generated by the threshold shift model.

Figure 4.

ROC curves for weak TS condition (δ _pq = 0.16)

Figure 5.

ROC curves for medium TS condition (δ _pq = 0.34)

Figure 6.

ROC curves for strong TS condition (δ _pq = 0.66)

In Figure 4, the results for weak threshold shift of only 0.16 standard units, show that there is very little difference between the empirical curves and the diagonal line of the ROC plots (i.e., random rejection of the null hypothesis) for any of the statistics, which indicates low power. Even though some power is gained with increasing sample size, it is no more than 0.2 with 1,000 observations if 0.05 is used as the nominal level.

Figure 5 shows that when δ _pq = 0.34, there is, again, more power to detect LD in larger samples; the power is affected very little by test length. To illustrate, choosing 0.05 as the nominal level and only considering 50-item tests, the power of bifactor statistics is 0.11 for sample size 200, 0.28 for sample size 500, and 0.43 for sample size 1,000. The empirical ROC curves for all three statistics are very similar.

Figure 6 shows that as the threshold shift parameter increases to 0.66, the power of all statistics becomes higher (i.e., greater than 0.6) if sample size is 500 or 1,000, but remains low/moderate (i.e., about 0.25) for small samples.

To summarize, the ROC curves of all three statistics are roughly identical, which reflects their similar performance in all the TS conditions. The power increases when sample size increases, but remains about the same when test length changes.

Bifactor Conditions

Figures 7 to 9 show the ROC curves for the three statistics with data generated by bifactor LD model.

Figure 7.

ROC curves for weak BF condition (a_pq = 0.54)

Figure 8.

ROC curves for medium BF condition (a_pq = 0.96)

Figure 9.

ROC curves for strong BF condition (a_pq = 1.67)

Figure 7 shows that when the impact of secondary factor is weak (a_pq = 0.54; error covariance is 0.09), the power is relatively low or moderate for all combinations of test lengths and sample sizes; the power at α = 0.05 ranges from (approximately) 0.1 to 0.3. However, the pattern of results is similar: Power increases as the number of observations increases and is not influenced very much by the number of items; all three statistics have nearly the same ROC curves.

Figure 8 shows that the power at nominal level 0.05 falls between 0.4 and 0.8 for all the statistics when a_pq = 0.96 on average (error covariance is 0.25), which is high compared with the medium TS conditions. Again, all three statistics yield similar results.

Figure 9 shows that when the residual association is strong (i.e., a_pq = 1.67, error covariance is 0.49), the power for all statistics is high for all conditions; especially with large samples, all three procedures have nearly perfect power to detect LD, which is shown by their extremely steep ROC curves.

In all, the results for BF conditions are similar to those obtained for TS conditions; the only difference is that all statistics seem to be more powerful. That may be attributed to the approximate formula (i.e., Equation 19) used to transform loadings to threshold shift parameters, which may lead to the unmatched LD levels between the threshold shift and bifactor generating models.

Discussion

For locally independent data, the empirical distributions of both score test statistics closely resemble $χ_{1}^{2}$ , especially when sample size is large and the number of items is small (which means fewer parameters are estimated). However, $χ_{1}^{2}$ as the null distribution should be used with caution in smaller samples and/or with longer tests, because the results suggest that both score test statistics tend to be liberal. In contrast, LD X² is conservative if treated as $χ_{1}^{2}$ , which is consistent with the results found by Chen and Thissen (1997); however, it might be a good alternative for small samples and/or long tests. Research is needed to provide a better approximation of LD X²’s asymptotic distribution.

With respect to power, all three statistics exhibit similar patterns in terms of ROCs, whether a secondary factor or a threshold shift is the source of LD. Both score test statistics are very close numerically. It may suggest that the two mathematically distinct LD models are not empirically distinguishable; however, further investigation is required before any conclusion can be reached.

In summary, all three statistics considered in this study provide useful diagnostic information about LD. LD X² is the easiest to compute; its power is comparable with the score statistics if its null distribution is known. In practice, however, unless resampling methods (e.g., jackknife or bootstrap) are invoked, it is in general not easy to approximate its null distribution. Both score statistics behave similarly in all the conditions considered in the present study; however, the computation of the bifactor statistic S_b requires numerically integrating over one more dimension, which makes it computationally more expensive than the threshold shift statistic S_t.

In future research, it would be interesting to generalize both score test statistics to polytomous IRT models, for example, the graded response model (Samejima, 1969). The generalization of bifactor statistic is straightforward; that involves adding another dimension to the original model. In contrast, because graded response model has more than one threshold, there are complications in generalizing the threshold shift statistic. Suppose the item pair has K response categories; one possibility is to incorporate different shift parameters for each of the K − 1 nonzero responses of the first item on each of the K − 1 thresholds of the second item, which leads the number of additional parameters to increase from 1 to (K − 1)². It might be unwieldy to explain each of these shift parameters; nevertheless, it is acceptable to compute a multivariate score test statistic (i.e., with limiting distribution, $χ_{{(K - 1)}^{2}}^{2}$ ) without obtaining actual parameter estimates. Other possible extensions involve imposing constraints on the (K − 1)² shift parameters; for example, shift parameters on the same threshold of the second item may be constrained to be equal, which results in a K − 1 dimensional score test statistic. Further simulation studies are needed to evaluate these possible alternatives. Extension of the LD X² statistic to the polytomous case is less ambiguous; see Maydeu-Olivares and Liu (2012) for some simulation results if a different reference distribution is used instead of the one for the independence model. Further simulation studies are necessary to compare alternative reference distributions.

As one of the reviewers pointed out, further simulation studies considering multiple LD items should be done to investigate the performance of the proposed LD indices under more realistic circumstances. Following the paradigm of Kim et al. (2011) in this situation, the false-detection rates with respect to locally independent pairs in the test should also be incorporated as a criterion of evaluation. Nevertheless, it should be noticed that introducing multiple LD pairs (or subsets) blurs the boundary between the primary latent variable representing ability or trait and the secondary latent variables representing LD. Caution is required when interpreting simulation results of this kind.

Conclusions that can be drawn from the results of the current study are that (a) LD X² is the easiest to compute, performs well, and has been implemented in commercial software; (b) score test statistic S_b based on the threshold shift model is also easy to compute, and works better than LD X² with larger samples and shorter tests; (c) score test statistic S_t based on the bifactor model behaves as well as S_t—however, S_t is still preferred for 2PL model due to its computational efficiency; and (d) despite being distinct statistical models with different likelihood functions, the bifactor model and the threshold shift model are essentially indistinguishable as in asymptotic tests.

Footnotes

Appendix A

Appendix B

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.

Notes

References

Bock

R. D.

(2010). Notes on parameter estimation for polytomous item responses models. Unpublished manuscript.

Bock

R. D.

Lieberman

(1970). Fitting a response model for N dichotomously scored items. Psychometrika, 35, 179-197.

Bradlow

Wainer

Wang

(1999). A Bayesian random effects model for testlets. Psychometrika, 64, 153-168.

Braeken

Tuerlinckx

De Boeck

(2007). Copula functions for residual dependency. Psychometrika, 72, 393-411.

Cai

(2010). A two-tier full-information item factor analysis model with applications. Psychometrika, 75, 581-612.

Cai

Thissen

du Toit

S. H. C.

(2011). IRTPRO for Windows. Lincolnwood, IL: Scientific Software International.

Chen

W.-H.

Thissen

(1997). Local dependence indexes for item pairs using item response theory. Journal of Educational and Behavioral Statistics, 22, 265-289.

Gibbons

Hedeker

(1992). Full-information item bi-factor analysis. Psychometrika, 57, 423-436.

Glas

C. A.

Suárez Falcón

J. C.

(2003). A comparison of item-fit statistics for the three-parameter logistic model. Applied Psychological Measurement, 27, 87-106.

10.

Hoskens

De Boeck

(1997). A parametric model for local dependence among test items. Psychological Methods, 2, 261-277.

11.

(2001). Testing for local dependency in dichotomous and polytomous item response models. Psychometrika, 66, 109-132.

12.

Jannarone

(1986). Conjunctive item response theory kernels. Psychometrika, 51, 357-373.

13.

Kelderman

(1984). Loglinear Rasch model tests. Psychometrika, 49, 223-245.

14.

Kendall

Stuart

(1961). The advanced theory of statistics. Vol. II. London, England: Griffin.

15.

Kim

De Ayala

R. J.

Ferdous

A. A.

Nering

M. L.

(2011).The comparative performance of conditional independence indices. Applied Psychological Measurement, 35, 447-471.

16.

Lehmann

(1999). Elements of large-sample theory. New York, NY: Springer.

17.

Maydeu-Olivares

Liu

(2012). Item diagnostics in multivariate discrete data. Unpublished manuscript.

18.

McDonald

R. P.

(1982). Linear versus nonlinear models in item response theory. Applied Psychological Measurement, 6, 379-396.

19.

Muthén

(1978). Contributions to factor analysis of dichotomous variables. Psychometrika, 43, 551-560.

20.

R Development Core Team. (2010). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing.

21.

Samejima

(1969). Estimation of ability using a response pattern of graded scores (Monograph No. 17). Psychometrika.

22.

Stout

(1987). A nonparametric approach for assessing latent trait unidimensionality. Psychometrika, 52, 589-617.

23.

Thissen

Steinberg

(2009). Item response theory. In Millsap

Maydeu-Olivares

(Eds.), The SAGE handbook of quantitative methods in psychology (pp. 148-177). London, England: Sage.

24.

van der Linden

Glas

(2010). Statistical tests of conditional independence between responses and/or response times on test items. Psychometrika, 75, 120-139.

25.

Wirth

Edwards

M. C.

(2007). Item factor analysis: Current approaches and future directions. Psychological Methods, 12, 58-79.

26.

Yen

W. M.

(1984). Effects of local item dependence on the fit and equating performance of the three-parameter logistic model. Applied Psychological Measurement, 8, 125-145.