An Examination of the Linking Error Currently Used in PISA

ABSTRACT

Educational large-scale assessment (LSA) studies like the program for international student assessment (PISA) provide important information about trends in the performance of educational indicators in cognitive domains. The change in the country means in a cognitive domain like reading between two successive assessments is an example of a trend estimate. The uncertainty of trend estimates includes sampling and linking errors, which are regularly reported in the PISA study. This article focuses on the linking error that assesses the variability in trend estimation regarding the choice of items. Since PISA 2015, the linking error estimation method has changed. This article compares the statistical behavior and the concept of the new PISA linking error and the PISA linking error utilized until PISA 2012. It turns out that the newly proposed linking error provides no generalization of the old PISA linking error but rather reflects different aspects of uncertainties and model error.

KEYWORDS:

• Large-scale assessment
• trend estimation
• linking
• differential item functioning
• item parameter drift
• linking error
• PISA

Introduction

One major aim of international large-scale assessments (ILSAs) is to monitor changes in the levels of educational outcomes (e.g., student performance). Every three or four years, for example, the Programme for International Student Assessment (PISA) provides international comparisons of student performance in three content areas (reading, mathematics, and science; OECD, 2014). The successive assessment of these content domains makes it possible to estimate national trends within each participating country, which provides policymakers with important information for the evaluation of educational reforms. In order to be able to estimate these trends in student performance, repeated assessments need to be reported on a common scale that is comparable across time (Mazzeo & von Davier, 2014). To accomplish this task, a set of common items is repeatedly administered in each assessment, and linking methods are used to align the results from the different assessments on a common scale. When calculating trend estimates, however, it is important to identify the source and magnitude of errors affecting national trend estimation (Wu, 2010).

In order to estimate country trends in student achievement, the results from different assessments need to be linked so that the achievement scores in a respective cognitive domain can be directly compared. The general idea of a linking approach is to use a set of common items administered in more than one assessment in order to establish a common metric that makes it possible to compare the test results across the different assessments (Dorans et al., 2007; von Davier & Sinharay, 2014). Figure 1 illustrates a typical linking design used in two assessments of an ILSA study. In both assessments, a set of I₀ link items (also referred to as common or anchor items) is administered to a cohort of students (e.g., 15-year-old students in a country; OECD, 2014). In addition, I₁ and I₂ unique items are presented in only one of the two assessments. One advantage of including unique items in an ILSA is that they can be made publicly available for secondary analysis, and the item pool can be renewed in later assessments (Mazzeo & von Davier, 2014).

Figure 1. Linking in two PISA assessments.

In the present article, we are dealing with the quantification of uncertainty in trend estimates in country means. It has been argued that trend estimates are prone to sampling errors of persons and linking errors due to the selection of items in successive assessments (Wu, 2010). PISA changed the linking error estimation approach since PISA 2015 (OECD, 2014, 2018). There is a lack of research on how the newly proposed PISA linking error compares to the old linking error employed until PISA 2012. The present article has two main goals. First, we clarify with simulations and analytical derivations the error concept that is behind the new PISA linking error. It is shown that the new PISA linking error mainly quantifies error due to country differential item functioning. Second, we demonstrate in various settings that the new PISA linking error is not a generalization of the old linking error. This message is important for practitioners because official reporting indicates that the new linking error could replace the old linking error without essentially changing the interpretation of the linking error concept.

Scaling models for trend estimation

In this section, alternative item response models (i.e., scaling models) for trend estimation are introduced. In the one-parameter logistic (1PL; Rasch, 1960) model, the item response function for item i (i = 1, … ,I) in country c (c = 1, … ,C) at time t = 1,2 is given by

(1) $P_{\mathit{ict}}\left(\mathit{\theta }\right)=\mathrm{\Psi }\left(\mathit{\theta }-b_{\mathit{ict},\mathrm{model}}\right),\mathit{\theta }\sim \mathrm{N}\left({\mathit{\mu }}_{\mathit{ct}},{\sigma }_{ct}^2\right),\mathit{\hspace{0.17em}}\left(\mathit{t}=1,2\right).$(1)

The assumed item difficulties $b_{\mathit{ict},\mathrm{model}}$ can be time point-specific or country-specific. However, item parameters can also be assumed invariant across countries or both time points. Furthermore, the latent ability variable $\mathit{\theta }$ is normally distributed with mean ${\mathit{\mu }}_{\mathit{ct}}$ and standard deviation ${\sigma }_{ct}^{\hspace{0.17em}}$. In empirical applications like the PISA study, the assumed item difficulties in the scaling model will differ to some extent from true data-generating item difficulties $b_{\mathit{ict},\mathrm{true}}$. This difference can introduce additional variability or bias in country means and country standard deviations. The 1PL model for dichotomous data or the partial credit model (Masters & Wright, 1997) for polytomous data has been used in PISA until PISA 2012 (OECD, 2014).

Since PISA 2015, the 2PL model (Birnbaum, 1968) for dichotomous items or the generalized partial credit model (Muraki, 1997) for polytomous items has been utilized. The item response function for the 2PL model is given by

(2) $P_{\mathit{ict}}\left(\mathit{\theta }\right)=\mathrm{\Psi }\left(a_{\mathrm{ict},\mathrm{model}}\left(\mathit{\theta }-b_{\mathrm{ict},\mathrm{model}}\right)\right),\mathit{\theta }\sim \mathrm{N}\left({\mathit{\mu }}_{\mathit{ct}},{\sigma }_{ct}^2\right),\mathit{t}=1,2,$(2)

where $a_{\mathit{ict},\mathrm{model}}$ are item discrimination parameters.

The partial invariance assumption (Oliveri & von Davier, 2011, 2014; von Davier et al., 2019) relies on the idea that most parameters are invariant across countries and time points. That is, it holds in the 1PL model that $b_{\mathit{ict},\mathrm{model}}=b_i$ with common item difficulties $b_i$. In PISA, for most of the items, item discriminations are assumed to be invariant across countries and time in the 2PL model (Joo et al., 2021; OECD, 2018, 2020).

Item fit statistics are utilized to assess whether the invariance assumption of item parameters can be justified. The mean deviation (MD; OECD, 2018) statistic for item i in a particular group g is defined as

(3) $\mathit{M}{D}_{\mathit{ig}}=\int \left(P_i,\mathrm{obs}\left(\mathit{\theta }\right)-P_i,\mathrm{model}\left(\mathit{\theta }\right)\right)f_g\left(\mathit{\theta }\right)\mathit{d\theta },$(3)

where $P_{\mathit{i},\mathrm{model}}$ is the model-implied item response function, $P_{\mathit{i},\mathrm{obs}}$ is the empirical (i.e., observed) item response function, which estimates the true item response function, and $f_g$ denotes the empirical density in group g. The MD statistic in (3) can be interpreted as a weighted average of the difference $P_{\mathit{i},\mathrm{obs}}\left(\mathit{\theta }\right)-P_{\mathit{i},\mathrm{model}}\left(\mathit{\theta }\right)$ in item response functions. In the MD statistic, positive and negative deviations in $P_{\mathit{i},\mathrm{obs}}\left(\mathit{\theta }\right)-P_{\mathit{i},\mathrm{model}}\left(\mathit{\theta }\right)$ can cancel on average.

In PISA, the root mean square deviation (RMSD; OECD, 2018) statistic is frequently used and defined as

(4) $\mathrm{RMS}{\mathrm{D}}_{\mathit{ig}}=\sqrt{\int {({\mathit{P}}_{\mathit{i},\mathrm{obs}}(\mathit{\theta })-P_{\mathit{i},\mathrm{model}}(\mathit{\theta }))}^2{f}_g(\mathit{\theta })\mathrm{d}\mathit{\theta }}.$(4)

The RMSD statistic only takes non-negative values and quantifies the weighted variability of the difference $P_{\mathit{i},\mathrm{obs}}\left(\mathit{\theta }\right)-P_{\mathit{i},\mathrm{model}}\left(\mathit{\theta }\right)$. One can show that $\left|\mathrm{M}{\mathrm{D}}_{\mathit{ig}}\right|\le \mathrm{RMS}{\mathrm{D}}_{\mathit{ig}}$ (Robitzsch & Lüdtke, 2020). If the true item response function is given by $\mathrm{\Psi }\left(\mathit{\theta }-b_{\mathit{i},\mathrm{true}}\right)$ and there were no sampling errors, we would have $\mathrm{RMS}{\mathrm{D}}_{\mathit{ig}}=\left|\mathrm{M}{\mathrm{D}}_{\mathit{ig}}\right|$, which is the reason why we only rely on the MD statistic in this article. Computational details regarding the MD and RMSD statistics can be found in Köhler et al. (2020), Robitzsch (2022), or Tijmstra et al. (2020).

We now discuss the behavior of the MD statistic in a misspecified 1PL model. Assume that $b_{\mathit{i},\mathrm{model}}$ denote the assumed item difficulty used for scaling and $b_{\mathit{i},\mathrm{true}}$ is the true data-generating item difficulty. If f_g is the normal density with mean ${\mathit{\mu }}_g$ and standard deviation ${\mathit{\sigma }}_g$, it was shown by Robitzsch and Lüdtke (2020) that the MD statistic could be approximated by

(5) $\left|\mathit{M}{D}_{\mathit{ig}}\right|=\left|\mathrm{\Phi }\left(\frac{{\mathit{\mu }}_g-b_i,\mathrm{model}}{\sqrt{D^2+{\sigma }_g^2}}\right)-\mathrm{\Phi }\left(\frac{{\mathit{\mu }}_g-b_i,\mathrm{true}}{\sqrt{D^2+{\sigma }_g^2}}\right)\right|,\mathrm{where}\mathit{D}=\mathrm{1.701.}$(5)

We can apply a Taylor expansion with respect to item difficulties in (5), and we obtain

(6) $\left|\mathit{M}{D}_{\mathit{ig}}\right|=\frac{1}{\sqrt{D^2+{\sigma }_g^2}}\text{phii}\left(\frac{{\mathit{\mu }}_g-b_{\mathit{i},\mathrm{model}}}{\sqrt{D^2+{\sigma }_g^2}}\right)\left|b_{\mathit{i},\mathrm{model}}-b_{\mathit{i},\mathrm{true}}\right|.$(6)

By assuming ${\mathit{\sigma }}_g$ = 1, we can determine an upper bound of $\left|\mathrm{M}{\mathrm{D}}_{\mathit{ig}}\right|$ in (6) by determining ${\left(D^2+{\sigma }_g^2\right)}^{-1/2}$ and using $\mathit{\varphi }\left(\mathit{x}\right)\le \mathit{\varphi }\left(0\right)=0.399$ as

(7) $\left|{\mathrm{MD}}_{\mathit{ig}}\right|\le 0.202\cdot \left|b_{\mathit{i},\mathrm{model}}-b_{\mathit{i},\mathrm{true}}\right|.$(7)

For standard deviations ${\mathit{\sigma }}_g$ larger than 1, the multiplication factor in (7) is smaller than 0.202; for standard deviations smaller than 1, it is slightly larger than 0.202. The bound in (7) can only be attained if the group distribution density f_g matches the item difficulty, that is, it holds that ${\mathit{\mu }}_g-b_{\mathit{i},\mathrm{model}}=0$. Formulas (5), (6), and (7) can be used to the size of unmodelled DIF or IPD effects as the discrepancy $b_{\mathit{i},\mathrm{model}}-b_{\mathit{i},\mathrm{true}}$ into the size of the MD statistic. Hence, a cutoff value of the MD statistic of 0.20 requires at least a value of 1.00 of the deviation |$b_{\mathit{i},\mathrm{model}}-b_{\mathit{i},\mathrm{true}}|$. Thus, using a large cutoff MD cutoff of 0.20 is frequently very similar, if not identical, to the full invariance approach in which all parameters are assumed to be invariant across groups or countries because very large deviations $b_{\mathit{i},\mathrm{model}}-b_{\mathit{i},\mathrm{true}}$ rarely occur. It can be argued that model deviations that are represented by the difference $b_{\mathit{i},\mathrm{model}}-b_{\mathit{i},\mathrm{true}}$ become less important if items have extremely small or high item difficulties compared to the center of the θ distribution (see Tijmstra et al., 2020).

Linking errors for trend estimation

In this section, a variance component model for item difficulties in the 1PL model is presented that allows the derivation of linking errors. The model follows the framework presented by Robitzsch and Lüdtke (2019). The item response function in the 1PL model for item i in study t = 1,2 for country c is given by

(8) $\mathrm{logit}\mathrm{P}(X_{\mathit{ict}}=1|\mathit{\theta })=\mathit{\theta }-b_{\mathit{ict}},$(8)

where item difficulties $b_{\mathit{itc}}$ are centered for each country c in study t. The country mean ${\mathit{\mu }}_{\mathit{tc}}$ and the country standard deviation ${\mathit{\sigma }}_{\mathit{tc}}$ of country c can be estimated. In most applications, statistical inference (i.e., standard errors or uncertainty quantification in general) for ${\mathit{\mu }}_{\mathit{tc}}$ and ${\mathit{\sigma }}_{\mathit{tc}}$ should be conducted.

The consequences of the violation of measurement invariance (i.e., assuming item difficulties $b_{\mathit{i},\mathrm{model}}$ in the scaling model that differ from the data-generating item difficulties $b_{\mathit{i},\mathrm{true}}$) can be distinguished regarding three different aspects. First, items can function differently across assessments, known as item parameter drift (IPD; Hanson & Béguin, 2002; Meade et al., 2005). Second, items can function differently (differential item functioning; DIF) across countries at one time point, indicating that an item is relatively easier or more difficult for a specific country than at the international level (Camilli, 2006). These cross-national differences have been studied extensively and termed country DIF (Oliveri & von Davier, 2011, 2014). Third, country DIF can vary across time points, meaning that the relative difficulty changes across assessments (DIF × IPD; Carstensen, 2013).

We now specify a variance component model for item effects for national item parameters $b_{\mathit{itc}}$ that parametrizes DIF, IPD, and DIF × IPD (see Robitzsch & Lüdtke, 2019)

(9) $b_{\mathit{ict}}=b_i+u_{\mathit{ic}}+v_{\mathit{it}}+w_{\mathit{ict}}.$(9)

International item parameters ${\mathit{\beta }}_{\mathit{it}}$ that are used for linking the scales at the two time points (see the linking INT1$\leftrightarrow$INT2 in Figure 1) are given by

(10) ${\mathit{\beta }}_{\mathit{it}}=b_i+v_{\mathit{it}},$(10)

where the item effects $u_{\mathit{ic}}$, $v_{\mathit{it}}$, and $w_{\mathit{ict}}$ have zero means and variances ${\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F},c}^2=\mathrm{Var}\left(u_{\mathit{ic}}\right)$, ${\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2=\mathrm{Var}\left(v_{\mathit{it}}\right)$, and ${\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F}\times \mathrm{I}\mathrm{P}\mathrm{D},c}^2=\mathrm{Var}\left(w_{\mathit{ict}}\right)$. Note that the variance components of DIF and DIF × IPD effects can be country-specific. In the following, the different approaches for assessing the linking error of the country mean trend estimate ${\widehat{\mathrm{\Delta \mu }}}_{\mathit{c}}={\hat{\mu }}_{\mathit{c}2}-{\hat{\mu }}_{\mathit{c}1}$ are discussed.

Old PISA linking error until 2012: assessment of IPD effects

Until PISA 2012, the old linking error only addresses the effects of IPD in linking items. The linking procedure in PISA used mean–mean linking (Kolen & Brennan, 2014; Monseur & Berezner, 2007), which equates the mean of international item parameters ${\mathit{\beta }}_{\mathit{it}}$ for t = 1 (i.e., T1) and t = 2 (i.e., T2). It should be emphasized that the item parameters ${\mathit{\beta }}_{\mathit{it}}$ are the same for each country c. The old PISA linking error LE that was operationally used until PISA 2012 is estimated as the empirical variance of the difference in the international parameters of the two successive time points

(11) $\widehat{L{E}^2}=\frac{1}{I_0-1}\sum _{\mathit{i}=1}^{I_0}{\left({\mathit{\beta }}_{\mathit{i}2}-{\mathit{\beta }}_{\mathit{i}1}\right)}^2$(11)

It can be shown that the expected value of the linking error estimate in (11) is given by (see Monseur & Berezner, 2007; Wu, 2010)

(12) $\mathit{L}{E}^2=\mathrm{E}\left(\widehat{L{E}^2}\right)=\frac{2}{I_0}{\sigma }_{IPD}^2$(12)

Thus, the old PISA linking error is a function of the effects of IPD (i.e., ${\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2$) and the number of link items (i.e., $I_0$). However, it does not take into account variations in trend estimation that are due to country DIF effects (i.e., ${\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F},c}^2$) and relative difficulty changes across assessments that are due to DIF×IPD effects (i.e., ${\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F}\times \mathrm{I}\mathrm{P}\mathrm{D},c}^2$). In former PISA studies, the linking error estimate in (11) has also been generalized to accommodate the testlet structure of items (Monseur & Berezner, 2007; Robitzsch & Lüdtke, 2019). The reason for preferring the more complex estimate is that IPD or DIF effects are often positively correlated for items within a testlet. Hence, ignoring the testlet structure will typically result in underestimated linking errors.

A linking error formula that accounts for DIF, IPD, and DIF×IPD effects

It has been shown in a simulation study by Sachse et al. (2016) that the old PISA linking error substantially underestimates the true linking error of country mean trend estimates. Robitzsch and Lüdtke (2019) analytically derived the linking error for the trend estimate in the country means for the 1PL model based on the variance component model defined in (9) and (10). The analytical derivation assumed mean–mean linking steps, although it can be expected that the linking errors based on (concurrent) marginal maximum likelihood estimation are very similar. The linking error of the trend estimate ${\widehat{\mathrm{\Delta }\mathit{\mu }}}_{\mathit{c}}={\hat{\mu }}_{\mathit{c}2}-{\hat{\mu }}_{\mathit{c}1}$ in the country means is given by

(13) ${\mathrm{LE}}^2=\mathrm{Var}({\widehat{\mathrm{\Delta }\mathit{\mu }}}_{\mathit{c}})=\frac{2}{I_0}{\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2+\frac{I_1+I_2}{(I_0+I_1)(I_0+I_2)}{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F},c}^2+\frac{2I_0+I_1+I_2}{(I_0+I_1)(I_0+I_2)}{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F}\times \mathrm{I}\mathrm{P}\mathrm{D},c}^2.$(13)

Importantly, only the first term in (13) (i.e., the variance due to IPD) is reflected in the old PISA linking error (see Equation (12)(12) $\mathit{L}{E}^2=\mathrm{E}\left(\widehat{L{E}^2}\right)=\frac{2}{I_0}{\sigma }_{IPD}^2$(12) ). In empirical applications, the DIF and DIF × IPD are also relevant contributions to the true linking error, which explains the reported underestimation of Sachse et al. (2016). It should be emphasized that the middle term in (13) that refers to the DIF variance component disappears if there are only link items, that is, I₁ = I₂ = 0. It can also be seen in the linking error formula in (13) that the linking error is a country-specific quantity. The variance components in (13) can be estimated for each country through linear mixed-effects models (Robitzsch & Lüdtke, 2019).

New PISA linking error since 2015: assessment of variability due to noninvariance

Since PISA 2015, a new linking error method for trend estimates in the country means has been established (OECD, 2018, 2020). The change in the linking error method was motivated by the assumption that the newly proposed linking error includes additional aspects in the uncertainty quantification of trend estimates (OECD, 2018). First, it assesses the effects of the choice of link items and unique items on trend estimates. In contrast, the old PISA linking error solely quantifies the IPD effects of link items. Second, it is argued that model error such as DIF that is not fully accounted for by modeling is more adequately reflected in the new linking error. Third, the new linking error assesses the effects of using different calibration samples for obtaining item parameters utilized in the scaling models. Fourth, changes in the scaling model, such as the switch from the 1PL to the 2PL model in PISA, were also claimed to be more adequately reflected in the new PISA linking error. Moreover, the new linking error is also intended to be applicable to the partial invariance approach that has been implemented since PISA 2015 (von Davier et al., 2019).

The linking error estimation relies on the recalibration (RC) approach in which a data set of the first time point is reestimated using a different set of item parameters (Martin et al., 2012). We formally denote by ${\hat{\mu }}_{\mathit{c}1}=\mathcal{C}\left({\mathbf{b}}_{\mathit{c}1},{\mathcal{D}}_{\mathit{c}1}\right)$ the original estimate of the country mean for country c at T1, where ${\mathbf{b}}_{\mathit{c}1}$ are item parameters for country c that were originally employed in the estimation of the country mean and ${\mathcal{D}}_{\mathit{c}1}$ denotes the dataset of item responses for country c at T1. The label $\mathcal{C}$ indicates the application of a calibration procedure. The recalibration typically involves item parameters that are used in the scaling models at T2 (OECD, 2018 see also Martin et al., 2012). Moreover, the country mean estimate under recalibration could also be based on only link items that appear at T1 and T2 (i.e., item set I₀ in Figure 1) by removing unique items at T1 from the estimation. We formally denote by ${\hat{\mu }}_{\mathit{c}1,\mathit{RC}}=\mathcal{C}\left({\mathbf{b}}_{\mathit{c}2},{\tilde{\mathcal{D}}}_{\mathit{c}1}\right)$ the country mean estimate at T1 under recalibration, where ${\mathbf{b}}_{\mathit{c}2}$ are the country-specific item parameters used at T2 and ${\tilde{\mathcal{D}}}_{\mathit{c}1}\subseteq {\mathcal{D}}_{\mathit{c}1}$ is a subset or a full set of item responses for country c at T1¹. The consequence of the choice of link items and unique items is also quantified in the difference ${\hat{\mu }}_{\mathit{c}1,\mathit{RC}}-{\hat{\mu }}_{\mathit{c}1}$. The new PISA linking error estimate is given by (see OECD, 2018, 2020)

(14) $\mathit{L}{E}^2=\frac{1}{C}\sum _{\mathit{c}=1}^C{\left({\hat{\mu }}_{\mathit{c}1,\mathit{RC}}-{\hat{\mu }}_{\mathit{c}1}\right)}^2.$(14)

It summarizes the differences in the country means due to different assumed item parameters and used items across all countries. Hence, the linking error defined in (14) applies to all countries, and no country-specific linking error estimates can be assessed using the new PISA linking error method. Notably, PISA replaces the estimated variance (or the standard deviation, respectively) in (14) with a robust standard deviation estimator (OECD, 2018). In the rest of this article, we only consider the nonrobust variant in (14) because using a robust estimator in the simulation studies and analytical derivations does not change the expected values of linking errors. In contrast, robust standard deviation estimates are only expected to substantially increase the efficiency of linking error estimates if the differences ${\hat{\mu }}_{\mathit{c}1,\mathit{RC}}-{\hat{\mu }}_{\mathit{c}1}$ strongly deviate from normality.

Assessing the total error for statistical inference of trend estimates

In all PISA studies, statistical inference for trend estimates in the country means included two aspects of errors, namely, the sampling error SE and the linking error LE (OECD, 2014, 2018; Wu, 2010). The total error TE is used to assess statistical uncertainty in trend estimates using the following formula (OECD, 2014):

(15) $\mathit{TE}=\sqrt{\mathit{S}{E}^2+\mathit{L}{E}^2}.$(15)

Note that (15) presupposes that the two error component estimates are uncorrelated. However, the linking error can also be affected by sampling error. Additional variability in the linking error due to sampling error can be eliminated by analytical or resampling approaches (Robitzsch, 2021; Robitzsch & Lüdtke, 2019).

In this article, we ignore the impact of sampling error in linking error estimation and apply the total error estimation formula (15). Note that PISA also ignores potential biases in linking error estimates due to sampling error, which is the reason why we also follow the approaches implemented in operational practice in this article.

Purpose

To our knowledge, there is currently no published literature on the performance of the newly proposed PISA linking error. Therefore, this article aims at understanding the behavior of the new PISA linking error compared to the old linking error used in PISA until 2012. Consequently, we restrict ourselves to the simpler 1PL model instead of the 2PL model because it was claimed that the new PISA linking error is a generalization of the old one. It has been shown in the Section “A Linking Error Formula That Accounts for DIF, IPD, and DIF × IPD Effects” that the true linking error for the country mean trend estimate depends on the three variance components DIF, IPD, DIF × IPD and the numbers of link and unique items. In the next two studies, we restrict the investigation of the new PISA linking error by studying the effects of IPD and DIF separately. The advantage of this approach is that the behavior of the newly proposed linking error can be understood in idealized settings. We believe that it is advantageous to try to understand the behavior of a new estimator in clearly defined but idealized settings instead of investigating complex settings that mimic real data applications because the latter confounds several aspects. In our view, it is frequently impossible to get a deep insight into the performance of a statistical approach by studying complex simulation findings. In Study 1, we investigate the behavior of the new linking error under IPD with only link items utilizing a simulation study and analytical arguments. In Study 2, the behavior of the new linking error is investigated under DIF in a test design with link items and unique items utilizing a simulation study and analytical arguments.

Study 1: item parameter drift in minor-minor design

The old PISA linking error only addressed IPD in link items. Study 1 investigates the performance of the new PISA linking error in minor–minor designs in PISA based on the 1PL scaling model in which only link items and no unique items exist (see Figure 2). Only IPD and no DIF effects were considered. First, a simulation study was carried out to study the behavior of the new linking error empirically. Afterward, analytical arguments are provided to deepen the understanding of the empirical findings.

Figure 2. Minor-minor design.

Method

Our simulation study comprised C = 30 countries at two time points. We took the country means and standard deviations from the PISA trend study between 2006 and 2009 for reading reported by Robitzsch and Lüdtke (2019). The country means at T1 (i.e., PISA 2006) had an average of 0.00 (SD = 0.21, Min = −0.32, Max = 0.46). The country means at T2 (i.e., PISA 2009) were slightly larger on average (M = 0.05, SD = 0.25, Min = −0.59, Max = 0.52). The trend estimates in the country means, defined as the difference between the means at the two time points, had an average of 0.05 (SD = 0.17, Min = −0.53, Max = 0.34). We assumed equal country standard deviations across the two assessments for all countries (M = 1.18, SD = 0.09, Min = −0.99, Max = 1.40). The data-generating parameters can be found at https://osf.io/vq5pn/.

All countries had the same sample size N. The sample sizes were varied in the simulation with factor levels N = 1000, 2500, and 5000. Item responses were simulated from the 1PL model. In total, I₀ = 30 link items were chosen. Common item difficulties b_i were defined equispaced in the interval [−2.0,2.0]. IPD effects $v_{\mathit{it}}$ were added to the common item difficulties (see Equation (9)(9) $b_{\mathit{ict}}=b_i+u_{\mathit{ic}}+v_{\mathit{it}}+w_{\mathit{ict}}.$(9) ). A mixture distribution was chosen to simulate IPD effects. The motivation for choosing a mixture distribution was to simultaneously accommodate unsystematically and normally distributed IPD effects as well as large IPD effects that could resemble a data-generating model fulfilling the partial invariance assumption. The data-generating model for IPD effects $v_{\mathit{it}}$ consists of a subset of items ${\mathcal{J}}_{\mathrm{random}}$ that have unsystematic (and small) IPD effects (random IPD) and a distinct set ${\mathcal{J}}_{\mathrm{bias}}$ of biased (and large) IPD effects (biased IPD; see Robitzsch & Lüdtke, 2022b, for a similar distinction),

(16) $v_{\mathit{it}}\sim (1-{\mathit{\pi }}_{\mathrm{bias}})\mathrm{N}(0,{\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2)+{\mathit{\pi }}_{\mathrm{bias}}{F}_{\mathrm{bias}}.$(16)

Hence, unsystematic IPD is reflected in the normal distribution mixture component $\mathrm{N}\left(0,{\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2\right)$. Outlying IPD effects are represented by a contaminating distribution $F_{\mathrm{bias}}$ that has thicker tails than $\mathrm{N}\left(0,{\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2\right)$. In order to avoid confounding properties of biased estimation with the assessment of variability in the linking error, we assume that IPD effects v_it have zero means. Hence, the mean of the distribution $F_{\mathrm{bias}}$ is zero (i.e., $\int \mathit{x}\mathrm{d}{F}_{\mathrm{bias}}\left(\mathit{x}\right)=0$).

In the simulation study, the proportion ${\mathit{\pi }}_{\mathrm{bias}}$ of biased IPD items was either 0 or 0.15. The distribution of biased IPD effects was chosen as a uniform distribution on [−0.8,−0.5]$$\cup$$[0.5,0.8]. In addition, the standard deviation σ_IPD of IPD effects was either 0.00 or 0.30. In total, there are four conditions in which random IPD and biased IPD effects were crossed (i.e., no IPD with ${\mathit{\pi }}_{\mathrm{bias}}$ = 0 and σ_IPD = 0.00; IPD random with ${\mathit{\pi }}_{\mathrm{bias}}$ = 0 and σ_IPD = 0.30; IPD biased items with ${\mathit{\pi }}_{\mathrm{bias}}$ = 0.15 and σ_IPD = 0.00; IPD random and biased items with ${\mathit{\pi }}_{\mathrm{bias}}$ = 0.15 and σ_IPD = 0.30).

International item parameters at T1 were obtained by applying the 1PL model to the pooled sample comprising all students in all C = 30 countries. In the next step, the pooled sample of all students was scaled using item difficulties fixed to the obtained values from T1. It was assessed by means of the MD statistic whether item parameters at T2 can be regarded as invariant. Items that showed a misfit received a unique item parameter at T2 in the scaling model, while all other items had item difficulties fixed to the corresponding estimate at T1. Four cutoff values for the absolute value of the MD statistic were used in this simulation study: 0.20, 0.12, 0.08, and 0.05. The condition of the MD cutoff of 0.20 mimics the full invariance condition because items were only flagged in rare cases. Note that in this simulation study, we only assessed item misfit at the international level (i.e., the pooled sample) and not at the country level for each of the two time points because only IPD and no DIF effects were simulated.

The standard error referring to person sampling was estimated as ${\hat{\sigma }}_{\mathit{ct}}/\sqrt{\mathit{N}}$, where ${\hat{\sigma }}_{\mathit{ct}}$ is the estimated standard deviation in the scaling model². The new PISA linking error utilizes the used item parameters at T2 for recalibrating data at T1. With a larger MD cutoff value, more item parameters in the recalibration were equal to those used in the original calibration at T1. The total error (TE) was assessed with the formula (15). We computed the mean and the standard deviation of the country mean trend estimates as well as the average of the estimated linking error using the new PISA method. Means and standard deviations (SDs) were averaged across the results of all 30 countries.

In each of the conditions, 300 replications were carried out. The statistical software R was used for implementing the entire simulation. The R package TAM (Robitzsch et al., 2022) was employed for estimating the multiple-group 1PL models and the computation of the MD statistic.

Results

Across all conditions, the means of country mean trend estimates were essentially zero. Hence, only the assessment of variability (i.e., the SD) and the estimation of TEs are interesting in this simulation study.

In Table 1, the empirical standard deviation and average estimated total errors are compared as a function of sample size for different MD cutoff values of the partial invariance scaling approach. In the absence of IPD (i.e., the “No IPD” condition), the average of estimated TEs was slightly smaller than SDs of trend estimates. Due to the very large sample sizes of the pooled sample comprising all countries, item parameters were estimated practically without sampling errors, and no items were removed under PI for all cutoff values 0.20, 0.12, 0.08, and 0.05. Hence, no differences in SDs and estimated TEs between the methods occurred.

Table 1. Simulation study 1: standard deviations and average-estimated total errors for country mean trend estimates in the presence of item parameter drift (σ_IPD = 0.30) and no country differential item functioning (σ_DIF = 0.00, σ_DIF×IPD = 0.00) for a minor–minor design, I₀ = 30 link items, and C = 30 countries as a function of sample size per country (N).

	N = 1000		N = 2000		N = 5000
	SD	TE	SD	TE	SD	TE
	No IPD
PI MD .20	0.057	0.053	0.040	0.037	0.025	0.024
PI MD .12	0.057	0.053	0.040	0.037	0.025	0.024
PI MD .08	0.057	0.053	0.040	0.037	0.025	0.024
PI MD .05	0.057	0.053	0.040	0.037	0.025	0.024
	IPD Random
PI MD .20	0.102	0.055	0.088	0.041	0.085	0.028
PI MD .12	0.105	0.067	0.096	0.054	0.091	0.044
PI MD .08	0.110	0.069	0.098	0.057	0.093	0.046
PI MD .05	0.109	0.066	0.095	0.053	0.092	0.043
	IPD Biased Items
PI MD .20	0.086	0.056	0.075	0.041	0.070	0.028
PI MD .12	0.093	0.061	0.084	0.049	0.081	0.039
PI MD .08	0.068	0.067	0.055	0.056	0.045	0.047
PI MD .05	0.057	0.079	0.041	0.068	0.026	0.060
	IPD Random and Biased Items
PI MD .20	0.117	0.062	0.114	0.047	0.108	0.038
PI MD .12	0.124	0.075	0.122	0.064	0.115	0.054
PI MD .08	0.124	0.075	0.121	0.062	0.113	0.054
PI MD .05	0.122	0.071	0.115	0.056	0.112	0.049

Note. PI MD = scaling under partial invariance assumption using the MD statistic with a particular threshold.

Under the condition “IPD Random,” only a few items were removed with an MD cutoff value of 0.20. In this case, the total error mostly reflects sampling error and not linking error. It can be seen that the SD of trend estimates was smallest with the MD cutoff of 0.20. Due to normally distributed IPD effects, some items with large IPD effects were treated as outliers. It is known from robust statistics that the removal of outliers from a normal distribution results in less efficient estimates. This is also demonstrated in the larger standard deviations with a smaller MD cutoff value. As can be seen, the variability in trend estimates is grossly underestimated using the newly proposed PISA linking error. Hence, if there is only random and unsystematic IPD, the newly proposed PISA linking error is not a generalization of the old PISA linking error, which can be unbiasedly estimated (see Monseur & Berezner, 2007; Robitzsch & Lüdtke, 2019)³.

In the condition that only involved biased IPD items (i.e., condition “IPD Biased Items”), it turned out that specifying unique item parameters for biased items increased the precision in the country mean trend estimates. For example, with N = 5000, the standard deviation decreases from 0.070 with an MD cutoff of 0.20 to 0.026 with an MD cutoff of 0.05. Notably, the smaller variability in trend estimates using the partial invariance approach with a more strict MD cutoff value was not reflected in the PISA linking error. The PISA linking error was low when the SD was high and the other way around. Hence, the new PISA linking error does not reflect variability in trend estimates under repeated sampling but rather the number of items that were assumed noninvariant. Consequently, the error concept behind the new PISA linking error since 2015 differs from the old PISA linking error until PISA 2012.

In the presence of both random IPD and IPD due to biased items (i.e., condition “IPD Random and Biased Items”), it was observed that the linking error strongly underestimates the true variability in trend estimates in the country means. It can be concluded that the currently used linking error showed a very distinct behavior from the linking error previously used in PISA. While the old linking error correctly quantifies the variability in trend estimates due to IPD variability (Monseur & Berezner, 2007), the newly proposed linking error quantifies the extent of items that are assumed noninvariant. Therefore, the recently used new PISA linking error does not adequately capture variability in trend estimates due to IPD but constitutes a summary measure of the extent of noninvariant items. It was argued that the latter concept refers to the comparability of scaling results across groups (Joo et al., 2021; von Davier et al., 2019). In this sense, the new PISA linking error contributes to the TE to the extent of limited comparability.

Analytical findings

We now try to get further insights into the behavior of the new PISA linking error by providing analytical arguments. The derivations rely on the substitution of marginal maximum likelihood estimation of the 1PL model with unweighted least squares estimation (Cai & Moustaki, 2018; Robitzsch & Lüdtke, 2020). Hence, the findings are only approximately generalizable to operational practice that uses maximum likelihood estimation. The estimated country mean ${\hat{\mu }}_{\mathit{ct}}$ for country c at time point t in infinite sample sizes of persons can be written as

(17) ${\hat{\mu }}_{\mathit{ct}}={\mathit{\mu }}_{\mathit{ct}}-\frac{1}{I}{\sum }_{i=1}^I\left(b_{\mathit{ict},\mathrm{true}}-b_{\mathit{ict},\mathrm{model}}\right),$(17)

where $b_{\mathit{ict},\mathrm{model}}$ is the assumed item difficulty in the scaling model and $b_{\mathit{ict},\mathrm{true}}$ refers to the true data-generating item difficulty. Note that additional variability is introduced in (17) if the differences $b_{\mathit{ict},\mathrm{true}}-b_{\mathit{ict},\mathrm{model}}$ differ from zero. Linking errors quantify the variability of these differences.

We now derive the expected value of the new PISA linking error. Like in the simulation study, we rely on the mixture distribution for IPD effects (see (16)). Assume that there is no random IPD ${\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2$, but there only exist biased IPD effects. Furthermore, assume that these biased items are detected in the partial invariance approach using appropriate cutoff values for the MD statistic. The set of all items is partitioned into a set ${\mathcal{J}}_{\mathrm{random}}$ that has no IPD effects and a set ${\mathcal{J}}_{\mathrm{bias}}$ that has (large) IPD effects. By defining the variance in $F_{\mathrm{bias}}$ as ${\kappa }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2$ (i.e., $\int x^2\mathrm{d}{F}_{\mathrm{bias}}\left(\mathit{x}\right)={\kappa }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2)$, we get the true linking error due to IPD

(18) ${\mathrm{LE}}^2=\frac{2{\mathit{\pi }}_{\mathrm{bias}}}{\frac{}{}}{I}_0{\kappa }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2.$(18)

We now evaluate the new PISA linking error. The difference between the recalibrated country mean ${\hat{\mu }}_{\mathit{c}1,\mathit{RC}}$ and the original country mean ${\hat{\mu }}_{\mathit{c}1}$ can be calculated as

(19) ${\hat{\mu }}_{\mathit{c}1,\mathrm{RC}}{-\hat{\mu }}_{\mathit{c}1}=\frac{1}{I_o}\sum _{\mathit{i}\in {\mathit{J}}_{\mathrm{random}}}\left(v_{\mathit{i}2}-v_{\mathit{i}1}\right)+\frac{1}{I_0}\sum _{\mathit{i}\in {\mathit{J}}_{\mathrm{bias}}}\left(v_{\mathit{i}2}-v_{\mathit{i}1}\right).$(19)

The first term on the right-hand side in (19) equals zero because only large DIF effects exist. The second term refers to the set of biased items ${\mathcal{J}}_{\mathrm{bias}}$ that differ from zero, and we get

(20) $\mathrm{E}{\left({\hat{\mu }}_{\mathit{c}1,\mathit{RC}}-{\hat{\mu }}_{\mathit{c}1}\right)}^2=\frac{2{\mathit{\pi }}_{\mathrm{bias}}}{I_0}{\kappa }_{IPD}^2.$(20)

Hence, the result in (20) coincides with the true linking error computed in (18). As shown in this simulation study, the new PISA linking error does not assess the variability in trend estimates of the partial invariance approach in a repeated sampling of items. Instead, it quantifies the extent of comparability that coincides with the old PISA linking error in the particular data-generating model of partial invariance. Note that this derivation rests on the strong assumption that IPD effects fulfill the partial invariance assumption, and the noninvariant items were all correctly flagged by the item fit statistic.

We now apply mean–mean linking of the item parameters at the international level obtained from separate scaling models at T1 and T2. This approach was operationally implemented in PISA until PISA 2012 (OECD, 2014). Again, we assume infinite sample sizes of persons in order only to address the consequences of the linking error. At T1, estimated item difficulties at the international level are $b_i+v_{\mathit{i}1}$. At T2, item difficulties are used that are, on average, the same as those at T1 because the mean–mean linking method is utilized. Hence, for scaling, one uses item difficulties $b_i+v_{\mathit{i}2}-I_0^{-1}{\sum }_{\mathit{i}=1}^{I_0}\left(v_{\mathit{i}2}-v_{\mathit{i}1}\right)$. Hence, we get using (17)

(21) ${\hat{\mathit{\mu }}}_{\mathit{c}1,\mathrm{RC}}^{\mathit{\hspace{0.17em}}}-{\hat{\mathit{\mu }}}_{\mathit{c}1}^{\mathit{\hspace{0.17em}}}=\frac{1}{{\mathit{I}}_0}\sum _{\mathit{i}=1}^{{\mathit{I}}_0}\left[{\mathit{b}}_{\mathit{i}}+{\mathit{v}}_{\mathit{i}1}-\left({\mathit{b}}_{\mathit{i}}+{\mathit{v}}_{\mathit{i}2}-I_0^{-1}\sum _{\mathit{i}=1}^{{\mathit{I}}_0}\left({\mathit{v}}_{\mathit{i}2}-{\mathit{v}}_{\mathit{i}1}\right)\right)\right]-\frac{1}{{\mathit{I}}_0}\sum _{\mathit{i}=1}^{{\mathit{I}}_0}\left[{\mathit{b}}_{\mathit{i}}+{\mathit{v}}_{\mathit{i}1}-\left({\mathit{b}}_{\mathit{i}}+{\mathit{v}}_{\mathit{i}1}\right)\right]=0.$(21)

Consequently, when using mean–mean linking and unweighted least squares estimation instead of maximum likelihood estimation as simplifying assumptions, we expect that recalibrated country means coincide with original country means. With maximum likelihood estimation and mean–mean linking, our empirical experience is that the difference between both estimates also turned out very close to zero. Hence, the new PISA linking error is approximately zero under mean–mean linking, while the true linking error can substantially differ from zero.

Study 2: cross-sectional differential item functioning

In the second study, the impact of cross-sectional country DIF on the new PISA linking error is investigated. It is known that country DIF introduces additional variability in trend estimates if there are unique items in two successive assessments. Only cross-sectional DIF with a common SD σ_DIF across countries and no IPD effects are considered in the simulation study and the analytical derivations.

Method

In this study, we consider a major-minor test design (see Figure 3) in which there are I₀ link items, I₁ unique items at T1, and no unique items at T2. We used the same data-generating parameters for country means and standard deviations (see https://osf.io/vq5pn/). In total, I₀ = 30 link items and I₁ = 30 unique items were chosen. Both item sets had common equispaced item difficulties on the interval [−2,2]. DIF effects were added to common item difficulties. The same mixture distribution for DIF effects as in Simulation Study 1 was utilized. Random DIF was simulated with an SD σ_DIF of either 0 or 0.30. Again, biased DIF effects were uniformly drawn from [−0.8, −0.5]$$\cup$$[0.5,0.8] with a proportion ${\mathit{\pi }}_{\mathrm{bias}}$ of either 0 or 0.15. Item misfit was detected at the country level separately at T1 and T2 using the four choices of cutoff values 0.20, 0.12, 0.08, and 0.05. Like in Simulation Study 1, there are four types of DIF effects in which the occurrence of random DIF and biased DIF was crossed.

Figure 3. Major-minor design.

First, item responses were scaled at T1 with the 1PL model using a pooled sample comprising all students. Next, MD item fit statistics were assessed for each country at T1. If the absolute value of the MD statistic exceeded the prespecified cutoff, a unique item parameter for this country was introduced. All remaining items had fixed item difficulties that equaled the obtained item difficulties at the international metric. Second, item responses at T2 were separately scaled with the 1PL model using the pooled sample. The two assessments were linked onto a common international metric using the mean–mean linking approach. Like at T1, the linked item parameters at T2 were utilized for assessing the item fit at the country level. Those items whose MD statistic exceeded the prespecified cutoff value received unique item parameters. In assessing the new PISA linking error, the recalibrated country means only relied on the 30 link items. The 30 unique items at T1 were omitted in the computation of the recalibrated country mean.

Like in Simulation Study 1, the mean and the SD of trend estimates and the average TE were computed. In total, 300 replications were conducted, and the R software (R Core Team, 2022) and the R package TAM (Robitzsch et al., 2022) were utilized.

Results

In Table 2, empirical standard deviations and the average estimated linking errors are presented as a function of sample size N. The condition in which DIF was absent (i.e., condition “No DIF”) showed the same pattern as the one observed under the absence of IPD (see Table 1).

Table 2. Simulation study 2: standard deviations and average estimated total errors for country mean trend estimates in the presence of country differential item functioning (σ_DIF = 0.50), and no item parameter drift (σ_IPD = 0.00, σ_DIF×IPD = 0.00) for a Major-minor design, I₀ = 30 link items, I₁ = 30 unique items at T1, and C = 30 countries as a function of sample size per country (N).

	N = 1000		N = 2000		N = 5000
	SD	TE	SD	TE	SD	TE
	No DIF
PI MD .20	0.056	0.054	0.039	0.038	0.025	0.024
PI MD .12	0.056	0.054	0.039	0.038	0.025	0.024
PI MD .08	0.056	0.054	0.039	0.038	0.025	0.024
PI MD .05	0.056	0.054	0.039	0.038	0.025	0.024
	DIF Random
PI MD .20	0.080	0.079	0.070	0.069	0.063	0.062
PI MD .12	0.073	0.072	0.063	0.062	0.054	0.054
PI MD .08	0.070	0.067	0.058	0.056	0.048	0.047
PI MD .05	0.068	0.063	0.054	0.050	0.043	0.041
	DIF Biased Items
PI MD .20	0.069	0.068	0.057	0.056	0.048	0.047
PI MD .12	0.065	0.064	0.052	0.051	0.042	0.042
PI MD .08	0.059	0.057	0.043	0.042	0.031	0.030
PI MD .05	0.057	0.055	0.041	0.039	0.027	0.026
	DIF Random and Biased Items
PI MD .20	0.083	0.084	0.075	0.075	0.069	0.069
PI MD .12	0.079	0.078	0.070	0.069	0.063	0.062
PI MD .08	0.076	0.072	0.065	0.062	0.057	0.055
PI MD .05	0.073	0.066	0.060	0.055	0.049	0.046

PI MD = scaling under partial invariance assumption using the MD statistic with a particular threshold.

In the presence of random cross-sectional DIF (i.e., condition “DIF Random”), it can be seen that it is useful to flag items with large DIF effects and assign them unique item parameters to reduce the variability in trend estimates. Note that this property only applies to major-minor designs if unique items exist. Interestingly, the newly proposed PISA linking error approximately reflects the observed variability in trend estimates due to country DIF in major-minor designs.

If there is only DIF due to biased items (i.e., condition “DIF Biased Items”), the new PISA linking error even better captures the variability due to DIF. Items with large cross-sectional DIF effects are essentially removed from the analysis, which resulted in decreased variability of trend estimates.

Finally, if both random DIF and biased DIF were present (i.e., condition “DIF Random and Biased Items”), using more strict MD cutoff values again reduced the variability in trend estimates. Also, the new PISA linking error successfully assessed the observed variability in the country mean trend estimates.

It should be emphasized that the successful behavior of the new PISA linking error depends on the decision that the recalibrated country means only involved link items and no unique items. Linking errors were underestimated if unique items were also used in the recalibration approach.

Analytical findings

We now want to validate the findings of the simulation study with analytical results. We assume that there only exist country DIF and no IPD effects. The difference between the original and the recalibrated country mean can be determined as

(22) ${\hat{\mathit{\mu }}}_{\mathit{c}1,\mathbf{RC}}^{\mathit{\hspace{0.17em}}}-{\hat{\mathit{\mu }}}_{\mathit{c}1}^{\mathit{\hspace{0.17em}}}=\frac{1}{{\mathit{I}}_0}{\sum }_{\mathit{i}=1}^{{\mathit{I}}_0}{\mathit{u}}_{\mathit{ic}}-\frac{1}{{\mathit{I}}_0+{\mathit{I}}_1}{\sum }_{\mathit{i}=1}^{{\mathit{I}}_0+{\mathit{I}}_1}{\mathit{u}}_{\mathit{ic}}=\frac{{\mathit{I}}_1}{{\mathit{I}}_0\left({\mathit{I}}_0+{\mathit{I}}_1\right)}{\sum }_{\mathit{i}=1}^{{\mathit{I}}_0}{\mathit{u}}_{\mathit{ic}}-\frac{1}{{\mathit{I}}_0\left({\mathit{I}}_0+{\mathit{I}}_1\right)}{\sum }_{\mathit{i}={\mathit{I}}_0+1}^{{\mathit{I}}_0+{\mathit{I}}_1}{\mathit{u}}_{\mathit{ic}}.$(22)

By taking the variance in (22), we get

(23) $\mathrm{E}{\left({\hat{\mathit{\mu }}}_{\mathit{c}1,\mathbf{RC}}^{\mathit{\hspace{0.17em}}}-{\hat{\mathit{\mu }}}_{\mathit{c}1}^{\mathit{\hspace{0.17em}}}\right)}^2=\frac{{\mathit{I}}_1}{{\mathit{I}}_0\left({\mathit{I}}_0+{\mathit{I}}_1\right)}{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F},c}^2.$(23)

The expected value of the square of the linking error estimate is given by

(24) $\mathrm{E}\left(\mathrm{L}{\mathrm{E}}^2\right)=\frac{1}{\mathit{C}}{\sum }_{c=1}^C\mathrm{E}{\left({\hat{\mathit{\mu }}}_{\mathit{c}1,\mathrm{RC}}^{\mathit{\hspace{0.17em}}}-{\hat{\mathit{\mu }}}_{\mathit{c}1}^{\mathit{\hspace{0.17em}}}\right)}^2=\frac{{\mathit{I}}_1}{{\mathit{I}}_0\left({\mathit{I}}_0+{\mathit{I}}_1\right)}{\bar{\mathit{\sigma }}}_{\mathrm{DIF}}^2,$(24)

where ${\bar{\sigma }}_{\mathrm{DIF}}^2=C^{-1}{\sum }_{c=1}^C{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F},c}^2$ denotes the average DIF variance across countries. The quantity in (24) corresponds to the linking error obtained in the variance component in Robitzsch and Lüdtke (2019) (see Equation (13)(13) ${\mathrm{LE}}^2=\mathrm{Var}({\widehat{\mathrm{\Delta }\mathit{\mu }}}_{\mathit{c}})=\frac{2}{I_0}{\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2+\frac{I_1+I_2}{(I_0+I_1)(I_0+I_2)}{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F},c}^2+\frac{2I_0+I_1+I_2}{(I_0+I_1)(I_0+I_2)}{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F}\times \mathrm{I}\mathrm{P}\mathrm{D},c}^2.$(13) ). If I₂ = 0 and no IPD and DIF × IPD variance components exist, the DIF variance components are equal across items. Thus, the new PISA linking error captures variability in trend estimates due to country DIF. Notably, the DIF variance component is not reflected in the old PISA linking error, so including this property is an advantage of the new PISA linking error.

The analytical finding in (24) can also be used to investigate the adequacy of the new PISA linking error in general test designs in which there exist unique items at T1 and T2 (i.e., I₁ > 0 and I₂ > 0). The multiplication factor of the DIF variance of the new PISA linking error in (24) is given by $I_1{I}_0^{-1}{\left(I_0+I_1\right)}^{-1}$. The true linking error (see Equation (13)(13) ${\mathrm{LE}}^2=\mathrm{Var}({\widehat{\mathrm{\Delta }\mathit{\mu }}}_{\mathit{c}})=\frac{2}{I_0}{\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2+\frac{I_1+I_2}{(I_0+I_1)(I_0+I_2)}{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F},c}^2+\frac{2I_0+I_1+I_2}{(I_0+I_1)(I_0+I_2)}{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F}\times \mathrm{I}\mathrm{P}\mathrm{D},c}^2.$(13) ) has a multiplication factor of $\left(I_1+I_2\right){\left(I_0+I_1\right)}^{-1}{\left(I_0+I_2\right)}^{-1}$. The multiplication factor of the difference between the new PISA linking error and the true linking error is given by

(25) $\frac{I_2(I_0-I_1)}{I_0(I_0+I_1)(I_0+I_2).}$(25)

Hence, the new PISA linking error correctly captures DIF effects if there are no unique items at T2 (i.e., I₂ = 0) or the number of link items I₀ equals the number of unique items I₁ at T1.

Discussion

In this article, we studied the behavior of the newly proposed linking error using two simulation studies and analytical arguments that rely on the 1PL model. It turned out that it does not satisfactorily assess variability in trend estimates due to IPD effects. Hence, the new linking error does not provide a generalization of the old PISA linking error that (only) assesses variability due to IPD. In the case of IPD, the only exception in which the new PISA linking error would prove useful might be the rare constellation in which IPD effects exactly follow the partial invariance assumption, and the items with IPD effects would be correctly flagged by item fit statistics.

Interestingly, the new PISA linking error satisfactorily captures variability in trend estimates due to country DIF effects. This finding was confirmed in the simulation study and an analytical derivation. It should be emphasized that the procedure should be adapted to minor-major test designs in which unique items only appear at T2. In this case, the linking error must be based on evaluating recalibrated and originally estimated country means at T2 instead of T1.

If both IPD and DIF are present, it can be concluded that the new PISA linking error results in biased estimates of the variability. However, in our experience, DIF effects are much more important than IPD effects in empirical studies like PISA, which, fortunately, will probably not result in strongly biased linking error estimates using the new method. However, practitioners and applied researchers should be aware that the newly proposed linking error depends on the extent of assumed noninvariant items. Hence, the estimated linking error depends on the arbitrarily chosen cutoff statistic for assessing item fit (see Robitzsch, 2022). Consequently, the new PISA linking error might be regarded as a measure of comparability instead of the variability of scale scores in trend estimates.

We want to reiterate that the new PISA linking error estimates an error of zero if all items were assumed invariant. It is hard to understand why a linking error should be exactly zero by assuming invariance. In contrast, we are convinced that linking errors also exist in scaling models that assume invariant item parameters. We would like to point out that the extent of linking errors is unrelated to whether one uses concurrent calibration, mean–mean linking, or scaling under partial invariance. Linking errors refer to variability regarding item choice. It is not essential whether deviations are modeled using the partial invariance approach or mean–mean linking or remain unmodelled in misspecified concurrent calibration.

The recalibration approach is computationally feasible for trend estimates in the country means. However, reporting linking errors for trend estimates in other parameters of interest, such as standard deviations, percentiles, or proportions, could also be warranted. In this case, it might be more general to apply resampling methods of item groups such as jackknife (Monseur & Berezner, 2007; Robitzsch & Lüdtke, 2022a; Xu & von Davier, 2010) or balanced half-sampling (Robitzsch, 2021). Plausible values could be drawn for each subsample of items in the resampling approach. This approach would provide a simple computational device for calculating linking errors for arbitrary parameters of interest. We have argued in Robitzsch and Lüdtke (2019) that linking errors should also be reported for cross-sectional statistics because country DIF also impacts cross-sectional quantities such as means, standard deviations, or percentiles. We find it inconsistent that linking errors are only (or primarily) reported for trend estimates in country means for PISA but are ignored for all other trend estimate parameters and cross-sectional descriptive statistics.

Notably, we restricted ourselves to analytical and simulation evidence of the different performances of the two linking errors based on the 1PL model. However, the conceptual differences also transfer to the 2PL model. In Robitzsch (2023), a linking error similar, but more intricate, to Equation (13)(13) ${\mathrm{LE}}^2=\mathrm{Var}({\widehat{\mathrm{\Delta }\mathit{\mu }}}_{\mathit{c}})=\frac{2}{I_0}{\sigma }_{\mathrm{I}\mathrm{P}\mathrm{D}}^2+\frac{I_1+I_2}{(I_0+I_1)(I_0+I_2)}{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F},c}^2+\frac{2I_0+I_1+I_2}{(I_0+I_1)(I_0+I_2)}{\sigma }_{\mathrm{D}\mathrm{I}\mathrm{F}\times \mathrm{I}\mathrm{P}\mathrm{D},c}^2.$(13) for the trend estimate with log-mean-mean linking based on the 2PL model was obtained utilizing the M-estimation theory. Further simulation research can be carried out to demonstrate that the newly proposed linking error in PISA also substantially differs from the 2PL linking error in Robitzsch (2023) that includes DIF in item intercepts and item slopes (i.e., in situations of uniform or nonuniform DIF). Nevertheless, we believe that linking errors for any statistic of interest should be more generally obtained through resampling methods regarding items or item groups (i.e., testlets).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1. If all items were used for computing the recalibrated country mean, linking errors would become very small if the number of link items is much smaller than the number of unique items at T1. The reason is that the same item parameters were used for unique items. Hence, only link items can have different item parameters in the recalibrated country mean, which, in turn, can only impact the estimated linking error.

2. Using ${\hat{\sigma }}_{\mathit{ct}}/\sqrt{\mathit{N}}$ instead of the log-likelihood based standard error results in a slight underestimate of the true standard error regarding person sampling. However, this article mainly addresses the assessment of linking error and the degree of biased standard errors vanishes in a large sample sizes N of persons per country.

3. We do not present the old PISA linking error in this simulation.