Standard Error Estimation for Subpopulation Non-invariance

Abstract

Score linking is widely used to place scores from different assessments, or the same assessment under different conditions, onto a common scale. A central concern is whether the linking function is invariant across subpopulations, as violations may threaten fairness. However, evaluating subpopulation differences in linked scores is challenging because linking error is not independent of sampling and measurement error when the same data are used to estimate the linking function and to compare score distributions. We show that common approaches involving neglecting linking error or treating it as independent substantially overestimate the standard errors of subpopulation differences. We introduce new methods that account for linking error dependencies. Simulation results demonstrate the accuracy of the proposed methods, and a practical example with real data illustrates how improved standard error estimation enhances power for detecting subpopulation non-invariance.

Keywords

linking equating standard error subpopulation invariance mode bridging

Introduction

Alternate assessments are frequently linked to support score comparability across different administrations or across distinct assessments of the same construct (Battauz & Leôncio, 2023; Deng & Rios, 2022; Leôncio et al., 2023; Liu & Jurich, 2023; Moses, 2022; Wallmark et al., 2023). Changes in administration format, such as transitioning from paper-based to digital delivery, shifting from in-person testing to remote proctoring, or moving from non-adaptive to adaptive test designs, may also require linking across assessments (Anselmi et al., 2023; Castellano et al., 2023, 2024; Chao & Chen, 2023; Fishbein et al., 2018; Jewsbury et al., 2020; Jones et al., 2022; Tian & Choi, 2023; Weirich et al., 2024; Wyse, 2023). Linking may also be applied across different assessments, such as linking national to international tests, or state to national assessments (Lim & Sireci, 2017; Reardon et al., 2021). Looking ahead, emerging assessment paradigms featuring increased personalization, dynamic interactivity, and AI-driven functionalities will necessitate linking procedures for which the fairness and validity across diverse subpopulations should be empirically evaluated (Arslan et al., 2024; Bulut et al., 2024; Hao et al., 2024; Johnson, 2025; Sinharay & Johnson, 2024; Wu et al., 2025).

The gold standard data collection design for linking is the random-groups design, in which the examinees assessed on alternate assessments or administration conditions are randomly equivalent (Dorans, 2018; Kim et al., 2010). Depending on the study design, examinees or other sampling units such as schools may be assigned to multiple conditions, resulting in dependent samples. Alternatively, when no sampling units overlap across conditions, the samples may be independent. Assuming the assessments measure the same underlying construct, linking is used to align the score distributions across groups, typically by equating sample moments or percentiles of the score distributions (Kolen & Brennan, 2014).

Because linking aligns scores at the aggregate level, it is essential to evaluate score comparability for subpopulations to ensure that assessment or administration changes do not systematically benefit or disadvantage any major subpopulation. The extent to which the linking function holds across subpopulations, referred to as subpopulation invariance, can be empirically assessed to evaluate whether the link supports fair and valid comparisons for all groups (Brennan, 2008; Dorans, 2004; Dorans et al., 2008; Huggins, 2014; Liu & Dorans, 2013; Petersen, 2008; Sinharay & Haberman, 2014).

For illustrative purposes, consider scores from a new assessment that have been linearly linked to the metric of a prior assessment, or the target metric. Let ${\hat{s}}_{Y}$ denote the transformed scores on the new assessment, defined by the linear function

{\hat{s}}_{Y} = \hat{a} s_{θ} + \hat{b},

(1)

where

s_{θ}

are the original scores on the new assessment, and

\hat{a}

and

\hat{b}

are the estimated linking coefficients. Under the random-groups design, these coefficients are computed as

\hat{a} = \frac{{\hat{σ}}_{X}}{{\hat{σ}}_{θ}},

(2)

\hat{b} = {\hat{μ}}_{X} - \hat{a} {\hat{μ}}_{θ},

(3)

where

{\hat{μ}}_{X}

and

{\hat{σ}}_{X}

are the sample estimates of the mean and standard deviation for the prior assessment (in the target metric), and

{\hat{μ}}_{θ}

and

{\hat{σ}}_{θ}

are the corresponding estimates for the new assessment (in their original metric, the source metric).

Subpopulation invariance may be evaluated by comparing summary statistics computed for the same subpopulation across the prior and new assessments (e.g., Yin, Bezirhan, et al., 2023). Let

\hat{D} = \hat{X} - \hat{Y},

(4)

where

\hat{X} = f (s_{X})

denotes a summary statistic (e.g., a mean) computed from the prior assessment scores

s_{X}

, and

\hat{Y} = f ({\hat{s}}_{Y})

denotes the corresponding summary statistic computed from the linked scores

{\hat{s}}_{Y}

on the new assessment. For example,

\hat{X}

and

\hat{Y}

might represent the means of girls’ scores on the two assessments. Subpopulation invariance can be assessed by testing the null hypothesis

\hat{D} = 0

. A statistically significant result provides evidence of subpopulation non-invariance.

As the linking coefficients Equations (2) and (3) are estimated from sample data, they are subject to sampling variability. This introduces uncertainty into the linked scores ${\hat{s}}_{Y}$ and any statistics derived from them, such as $\hat{Y}$ . This additional source of uncertainty is commonly referred to as linking variance or equating error (A. A. von Davier, 2013). Accordingly, linking variance should be incorporated into the standard error when testing subpopulation invariance using the difference $\hat{D}$ . However, in practice, analyses often fail to account for linking variance in the standard error or omit reporting standard errors altogether (e.g., Huggins & Elbaum, 2013; Ikeda et al., 2025; Inal & Arıkan, 2017; von Davier, 2006; Yin, Fishbein, et al., 2023).

While most approaches for quantifying linking variance focus on individual scores (e.g., Moses & Zhang, 2011; A. A. von Davier, 2013), some methods have been developed to estimate linking variance for summary statistics. Li (2022) proposed a bootstrap estimator for the standard error attributed to equating error on means. Mazzeo et al. (2018) introduced a Monte Carlo estimator for a linking variance component on means and other aggregate statistics. Johnson (1998) derived a first-order delta method approximation for the linking variance on means, and Kawahashi et al. (2024) derived a second-order delta method approximation for equating error affecting both scores and means.

However, as noted by each original author (Johnson, 1998; Kawahashi et al., 2024; Li, 2022; Mazzeo et al., 2018), these methods do not account for dependencies between linking error with sampling and measurement errors. Such dependencies may be particularly consequential when evaluating subpopulation invariance. Specifically, the same scores used to estimate the linking function coefficients are also used to compute the subpopulation estimates for the new assessment, $\hat{Y}$ , introducing dependence between the linking error with sampling and measurement errors on $\hat{Y}$ . Furthermore, because $\hat{Y}$ is computed from scores transformed using estimated linking coefficients that were also based on scores from the prior assessment, $\hat{Y}$ is dependent on $\hat{X}$ , which should be accounted for when estimating the variance of $\hat{D}$ .

Less work has addressed the dependencies between linking error and other sources of error. Jewsbury (2024) examined this issue and noted that such dependencies are difficult to model because their structure varies depending on the type of score comparison being made. To facilitate methodological development, Jewsbury proposed classifying score comparisons into categories that share common dependency structures, referred to as scenarios. Two example scenarios include: (1) post-bridge to bridge: comparing an administration of the new assessment to the administration used to estimate the linking coefficients (i.e., the bridge administration) and (2) post-bridge to pre-bridge: comparing an administration of the new assessment after the bridge to an administration of the prior assessment before the bridge.

Jewsbury (2024) proposed a taxonomy of methods that address dependencies in score comparisons. Total-based methods account for all sources of uncertainty simultaneously, component-based methods decompose the variance into independent components corresponding to distinct sources (e.g., independent samples), and difference-based methods correct standard error estimates for linking error. Within the total-based and component-based categories, a further distinction is made between equation-based methods (i.e., analytically derived plug-in estimators) and resampling-based methods (e.g., Monte Carlo, jackknife, or bootstrap), while difference-based methods are always equation-based.

While Jewsbury (2024) derived methods within each class in the taxonomy for several scenarios, their focus was on score comparisons over time, particularly in the context of reporting trends in large-scale assessments. As such, their methods cannot be directly applied to subpopulation invariance, a distinct type of score comparison distinct from the scenarios they considered. However, the principles underlying their taxonomy of methods can be extended to derive applicable methods for subpopulation invariance, which is the focus of this paper.

Beyond issues of linking error and its dependencies, a third issue is whether the new and prior assessments were administered in dependent or independent samples. Samples are dependent if they share units (e.g., schools or individuals). Our proposed component-based method assumes independence, while our proposed total-based method accommodates dependencies. Applying the component-based method to dependent samples may introduce bias due to assumption violation, whereas using the total-based method with independent samples may be inefficient.

Present Paper

In this paper, we introduce novel methods for subpopulation invariance based on Jewsbury’s (2024) taxonomy: total-based resampling, total-based equation, component-based resampling, component-based equation, and a difference-based correction-term method. The three equation-based methods are derived analytically, while the two resampling-based methods are formulated using the jackknife resampling technique. For comparison, we include the conventional method, which ignores linking error, as well as the linking component method, which accounts for linking error but neglects linking error dependencies. These seven methods are detailed in the Theory and Methods section.

In the Simulation section, we evaluate the statistical properties of the seven methods through simulation. In the Practical Example section, we apply these methods to real data to demonstrate the impact of appropriately accounting for linking error and linking error dependencies on inferences.

Theory and Methods

Prior Methods

For comparison with the methods novel to this paper, we consider two alternate classes of methods: first, conventional methods that ignore linking error; and second, linking variance component methods that account for linking error but neglect dependencies between linking error and other sources.

Conventional Methods

A conventional method refers to any variance estimation approach that does not account for linking error, such as standard or default estimators not specifically adapted to incorporate linking variance. By analogy, we define conventional sources of error as those captured by conventional methods; namely, sampling and measurement error. For example, the conventional application of the jackknife (e.g., Wolter, 2007) to estimate the variance of $\hat{Y}$ may be written as

{var}_{Conventional} (\hat{Y}) = \frac{N - 1}{N} \sum_{r = 1}^{N} {\{f (a^{(0)} s_{θ}^{(r)} + b^{(0)}) - \frac{1}{N} \sum_{r^{'} = 1}^{N} f (a^{(0)} s_{θ}^{(r^{'})} + b^{(0)})\}}^{2},

(5)

where

a

and

b

are the slope and intercept linking coefficients, respectively, and

s_{θ}

are the scores for the new assessment in the original metric. The function

f (\cdot)

denotes a statistical summary function, such as the mean. A superscript of (0) denotes the sample realization in the original sample, while a superscript of

(r)

denotes the sample realization in jackknife resample

r \in {1, \dots, N}

, where

N

is the total sample size and the number of jackknife resamples.

Note that the conventional method refers to applying a standard estimation procedure to the transformed scores of the new assessment, ${\hat{s}}_{Y}$ , without considering that the scores have been linked. Because linking has not been considered, the same sample estimates of the linking coefficients that were used to transform the scores are held constant across all resamples. In Equation (5) we write $a^{(0)} s_{θ}^{(r)} + b^{(0)}$ instead of ${\hat{s}}_{Y}^{(r)}$ to emphasize that the linking coefficients remain constant across the resamples.

The conventional estimator of the variance of $\hat{D}$ assuming independent samples is

{var}_{Conventional, ind} (\hat{D}) = var (\hat{X}) + {var}_{Conventional} (\hat{Y}),

(6)

where

var (\hat{X})

is the conventional estimate of the variance of

\hat{X}

. The subscript is excluded as there is no linking error on

\hat{X}

, so the conventional estimate is correct.

The conventional estimator of the variance of $\hat{D}$ assuming dependent samples is

\begin{align} {var}_{Conventional, dep} (\hat{D}) = \\ \frac{N - 1}{N} \sum_{r = 1}^{N} {\{f (s_{X}^{(r)}) - f (a^{(0)} s_{θ}^{(r)} + b^{(0)}) - \frac{1}{N} \sum_{r^{'} = 1}^{N} [f (s_{X}^{(r^{'})}) - f (a^{(0)} s_{θ}^{(r^{'})} + b^{(0)})]\}}^{2}, \end{align}

(7)

where

s_{X}

are the scores for the prior assessment in the prior metric.

Linking Component Methods

Second, the linking component approach estimates the linking variance on a statistic separately to the conventional sources of sampling and measurement, defined as the linking (variance) component. By assuming that linking error is independent to conventional sources, the sum of the estimated linking component and the conventionally estimated variance provides an estimate of the total variance (Jewsbury, 2019).

As described earlier, prior research on linking variance almost exclusively assumes independency of linking error and other sources (e.g., Johnson, 1998; Li, 2022; Mazzeo et al., 2018; Moses & Zhang, 2011). While Li (2022) uses a bootstrap estimator and Mazzeo et al. (2018) use a Monte Carlo estimator, we implement a jackknife estimator for consistency with other methods. The jackknife implementation of the linking component approach can be written as

\begin{align} {var}_{Link component} (\hat{D}) = {var}_{Conventional} (\hat{D}) \\ + \frac{N - 1}{N} \sum_{r = 1}^{N} {\{f (a^{(r)} s_{θ}^{(0)} + b^{(r)}) - \frac{1}{N} \sum_{r^{'} = 1}^{N} f (a^{(r^{'})} s_{θ}^{(0)} + b^{(r^{'})})\}}^{2}, \end{align}

(8)

where

{var}_{Conventional} (\hat{D})

is calculated with Equation (6) if samples are independent and Equation (7) if samples are dependent.

Notice that the linking coefficients are re-estimated within each resample, whereas the untransformed score is held constant. Consequently, the linking component accounts for linking variance but not conventional sources of error, whereas the conventional estimator Equations (6) and (7) accounts for conventional sources of error, but not linking variance.

Equations for the linking component can be found in Jewsbury (2019; see Equations (42) and (48)).

New Methods

In this section, we introduce several new methods under the taxonomy of Jewsbury (2024) that account for linking dependencies with other sources and apply to tests of subpopulation invariance. The derivations follow a similar approach as Jewsbury (2024), but in contrast to their methods that concern differences over time, our methods apply to subpopulation invariance, that is, Equation (4). In addition to the derivations of the equation-based methods, we also denote how resampling-based methods may be implemented with the jackknife.

For greater generality, we derive the equations for composite scores. Composite scores are weighted sums of subscale scores and are used by some assessments such as NAEP (NCES, 2025d). Specifically, a composite score is defined as

s_{X} = \sum_{i} w_{i} s_{X_{i}},

(9)

s_{θ} = \sum_{i} w_{i} s_{θ_{i}},

(10)

where

s_{X_{i}}

and

s_{θ_{i}}

are the subscale

i

scores for the prior assessment and the new assessment prior to transformation, respectively.

w_{i}

is the subscale weight.

In the case of composite scores, we assume that the link has been applied to the subscores. That is,

{\hat{s}}_{Y_{i}} = {\hat{a}}_{i} s_{θ_{i}} + {\hat{b}}_{i},

(11)

where

{\hat{a}}_{i} = \frac{{\hat{σ}}_{X_{i}}}{{\hat{σ}}_{θ_{i}}},

(12)

{\hat{b}}_{i} = {\hat{μ}}_{X_{i}} - {\hat{a}}_{i} {\hat{μ}}_{θ_{i}},

(13)

and where, for subscale

i

{\hat{μ}}_{X_{i}}

, and

{\hat{σ}}_{X_{i}}

are the estimators of the population mean and standard deviation, respectively, for the previous assessment in the target metric. Similarly,

{\hat{μ}}_{θ_{i}}

and

{\hat{σ}}_{θ_{i}}

are the corresponding estimators for the new assessment in the source metric.

To avoid separate derivations for each type of statistic that might be evaluated for subpopulation invariance (e.g., mean and standard deviation), we follow Jewsbury (2024) in deriving error variances for the type U estimator, defined generally to include many types of statistical summary functions that satisfy the expression

\hat{Y} = \sum_{i} ω_{i} (\hat{a} {\hat{θ}}_{i} + m {\hat{b}}_{i}),

where

\hat{Y} = f ({\hat{s}}_{Y})

is the estimator based on the new assessment in the metric of the prior assessment and

{\hat{θ}}_{i} = f (s_{θ, i})

is the estimator based on the new assessment in the metric of the new assessment for subscale

i

The general form of the type U estimator allows us to derive variance expressions for many commonly used statistics, such as means and percentiles, without rederiving each separately. Specifically, $ω_{1} = 1$ and $ω_{> 1} = 0$ for type AB and A estimators, $m = 1$ for type AB and CAB estimators, and $m = 0$ for type A and CA estimators. Type AB estimators include means and percentiles, type A estimators include standard deviations and estimators equal to differences between type AB or type A estimators, type CAB are estimators of means or percentiles of composite scores, and type CA are estimators equal to differences between type CAB estimators.

In the simulation and practical example, we focus on means for a unidimensional score. As noted, this classifies as type AB, so $ω_{1} = 1$ , $ω_{> 1} = 0$ , and $m = 1$ . However, the equations below are written for greatest generality to apply to other types of statistics.

Total-Based Approaches

Total-based methods directly evaluate the variance of an estimator while simultaneously accounting for all sources of variance, including linking and conventional sources such as sampling and measurement. Total-based resampling methods involves a resampling method such as the bootstrap or the jackknife, and accounts for linking variance by re-estimating and re-applying the linking function within each resample (for more general discussion, see Jewsbury, 2019, 2024). Total-based equation methods use expressions to estimate the variance of the estimator, considering all sources of error.

Total-Based Resampling

The jackknife implementation of total-based resampling may be written as

\begin{align} {var}_{Total resampling} (\hat{D}) = \\ \frac{N - 1}{N} \sum_{r = 1}^{N} {\{f (s_{X}^{(r)}) - f (a^{(r)} s_{θ}^{(r)} + b^{(r)}) - \frac{1}{N} \sum_{r^{'} = 1}^{N} [f (s_{X}^{(r^{'})}) - f (a^{(r^{'})} s_{θ}^{(r^{'})} + b^{(r^{'})})]\}}^{2} . \end{align}

(14)

Notice that both the linking coefficients and the untransformed score vary across the resamples. This contrasts to the linking component approach where one is held fixed at a time while the other is allowed to vary; see Equations (6) and (8). Both linking and conventional sources of uncertainty are accounted for simultaneously with total-based resampling, inherently accounting for dependencies (Jewsbury, 2024).

Total-Based Equation

The estimator Equation (4) may also be written as

\begin{array}{l} \hat{D} & = \sum_{i} ω_{i} ({\hat{X}}_{i} - {\hat{Y}}_{i}), \\ = \sum_{i} ω_{i} {\hat{σ}}_{X_{i}} ({\hat{Z}}_{X_{i}} - {\hat{Z}}_{θ_{i}}), \end{array}

where,

\begin{array}{l} {\hat{Z}}_{θ_{i}} & = \frac{{\hat{θ}}_{i} - m {\hat{μ}}_{θ_{i}}}{{\hat{σ}}_{θ_{i}}}, \\ {\hat{Z}}_{X_{i}} & = \frac{{\hat{X}}_{i} - m {\hat{μ}}_{X_{i}}}{{\hat{σ}}_{X_{i}}}, \end{array}

The variance of $\hat{D}$ may be written as

\begin{array}{l} var (\hat{D}) = \sum_{i} \sum_{j} ω_{i} ω_{j} {[covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) + σ_{X_{i}} σ_{X_{j}}] covar ({\hat{Z}}_{X_{i}} - {\hat{Z}}_{θ_{i}}, {\hat{Z}}_{X_{j}} - {\hat{Z}}_{θ_{j}}) \\ + (Z_{X_{i}} - Z_{θ_{i}}) (Z_{X_{j}} - Z_{θ_{j}}) covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{X_{j}})}, \end{array}

where variables without hats are the expectation of associated estimator, for example,

σ_{X_{i}} = E [{\hat{σ}}_{X_{i}}]

Alternatively, the variance of $\hat{D}$ may be written as

\begin{align} {var}_{Total equation} (\hat{D}) = \sum_{i} \sum_{j} ω_{i} ω_{j} { \\ [covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) + σ_{X_{i}} σ_{X_{j}}] \cdot [covar ({\hat{Z}}_{θ_{i}}, {\hat{Z}}_{θ_{j}}) + covar ({\hat{Z}}_{X_{i}}, {\hat{Z}}_{X_{j}}) - 2 covar ({\hat{Z}}_{θ_{i}}, {\hat{Z}}_{X_{j}})] \\ + (Z_{θ_{i}} Z_{θ_{j}} + Z_{X_{i}} Z_{X_{j}} - 2 Z_{θ_{i}} Z_{X_{j}}) \cdot covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{X_{j}})} . \end{align}

(15)

Component-Based Approaches

Component-based methods use Taylor series-like approximations to estimate components of variance, where each component represents an independent source of variance (see Jewsbury, 2024). For differences on linked scores, component-based methods are theoretically appropriate when the new and prior assessments are administered to independent samples. In this case, a component may be estimated to account for the variance of $\hat{D}$ due to each sample separately.

Note that the difference between the linking component method (Li, 2022; Mazzeo et al., 2018) and Jewsbury’s (2024) component-based approach is that the linking component method defines a component associated with linking error, regardless of whether it is actually independent of other sources. In contrast, Jewsbury’s (2024) component based approach defines components for each independent source of error. In the context of subpopulation invariance, the linking component is the variance of the linking error on $\hat{Y}$ in Equation (4). The component-based approach does not make this decomposition, as linking error on $\hat{Y}$ is not independent of other sources of error on $\hat{Y}$ and $\hat{X}$ .

Applying the component-based approach to subpopulation invariance, when the two samples are independent, produces one component associated with the new assessment sample and one component associated with the prior assessment sample. The former represents both conventional error on $\hat{Y}$ and part of the linking error on $\hat{Y}$ . The latter represents both conventional error on $\hat{X}$ as well as the other part of linking error on $\hat{Y}$ . Linking error is decomposed across the two parts because it is a function of both samples.

Component-Based Resampling

The jackknife implementation of component-based resampling may be written as

\begin{align} {var}_{Component resampling} (\hat{D}) = \\ \frac{N - 1}{N} \sum_{r = 1}^{N} {f (s_{X}^{(r)}) - f (\frac{σ_{X}^{(r)}}{σ_{θ}^{(0)}} (s_{θ}^{(0)} - μ_{θ}^{(0)}) + μ_{X}^{(r)}) \\ - \frac{1}{N} \sum_{r^{'} = 1}^{N} [f (s_{X}^{(r^{'})}) - f (\frac{σ_{X}^{(r^{'})}}{σ_{θ}^{(0)}} (s_{θ}^{(0)} - μ_{θ}^{(0)}) + μ_{X}^{(r^{'})})]}^{2} \\ + \frac{N - 1}{N} \sum_{r = 1}^{N} {f (s_{X}^{(0)}) - f (\frac{σ_{X}^{(0)}}{σ_{θ}^{(r)}} (s_{θ}^{(r)} - μ_{θ}^{(r)}) + μ_{X}^{(0)}) \\ - \frac{1}{N} \sum_{r^{'} = 1}^{N} [f (s_{X}^{(0)}) - f (\frac{σ_{X}^{(0)}}{σ_{θ}^{(r^{'})}} (s_{θ}^{(r^{'})} - μ_{θ}^{(r^{'})}) + μ_{X}^{(0)})]}^{2} . \end{align}

(16)

Notice that each component is defined by holding constant all random variables associated with one sample, while re-estimating the random variables associated with the other sample. Consequently, component-based resampling is expected to be more efficient than total-based resampling (Equation (14)) when samples are independent, but biased when samples are dependent.

Component-Based Equation

Component estimators may be defined by the new assessment and the prior assessment separately, such that

\begin{array}{l} {\hat{D}}_{X} & = \sum_{i} ω_{i} {\hat{σ}}_{X_{i}} ({\hat{Z}}_{X_{i}} - Z_{θ_{i}}), \\ {\hat{D}}_{θ} & = \sum_{i} ω_{i} σ_{X_{i}} (Z_{X_{i}} - {\hat{Z}}_{θ_{i}}), \end{array}

where

\begin{array}{l} Z_{X_{i}} & = \frac{X_{i} - m μ_{X}}{σ_{X}}, \\ Z_{θ_{i}} & = \frac{θ_{i} - m μ_{θ}}{σ_{θ}} . \end{array}

The variances of the components are

\begin{array}{l} var ({\hat{D}}_{X}) = \sum_{i} \sum_{j} ω_{i} ω_{j} \{[covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) + σ_{X_{i}} σ_{X_{j}}] covar ({\hat{Z}}_{X_{i}}, {\hat{Z}}_{X_{j}}) \\ + (Z_{θ_{i}} Z_{θ_{j}} + Z_{X_{i}} Z_{X_{j}} - 2 Z_{θ_{i}} Z_{X_{j}}) covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) \}, \end{array}

var ({\hat{D}}_{θ}) = \sum_{i} \sum_{j} ω_{i} ω_{j} σ_{X_{i}} σ_{X_{j}} covar ({\hat{Z}}_{θ_{i}}, {\hat{Z}}_{θ_{j}}) .

The component-based approach assumes the total variance is equal to the sum of the variance of the components. That is,

{var}_{Component equation} (\hat{D}) = var ({\hat{D}}_{X}) + var ({\hat{D}}_{θ}) .

(17)

Difference-Based Approaches: The Correction-Term Method

Difference-based methods aim to estimate the difference between the correct variance and the variance returned by a conventional method that neglects linking variance (Jewsbury, 2024, 2025). The primary motivation for developing difference-based approaches is that implementing them may be more convenient and less computationally intensive than other methods, especially if the conventionally estimated variance is already available. Consequently, we develop a difference-based method, the correction-term method for subpopulation invariance, that can be calculated without estimating any additional variances. Necessarily, this means that the correction-term method may be more approximate than other methods.

The estimator assumed by conventional methods is

\begin{array}{l} {\hat{D}}_{conv} & = \sum_{i} ω_{i} [{\hat{X}}_{i} - (a_{i} {\hat{θ}}_{i} + b_{i})], \\ = \sum_{i} w_{i} [{\hat{X}}_{i} - (\frac{σ_{X_{i}}}{σ_{θ_{i}}} ({\hat{θ}}_{i} - μ_{θ_{i}}) + μ_{X_{i}})] . \end{array}

Therefore, the variance of the estimator assumed by conventional methods is

var ({\hat{D}}_{conv}) = \sum_{i} \sum_{j} ω_{i} ω_{j} [covar ({\hat{X}}_{i}, {\hat{X}}_{j}) + a_{i} a_{j} covar ({\hat{θ}}_{i}, {\hat{θ}}_{j}) - 2 a_{i} covar ({\hat{θ}}_{i}, {\hat{X}}_{j})] .

The variance of the estimator may also be written as

\begin{array}{l} var ({\hat{D}}_{conv}) = \sum_{i} \sum_{j} ω_{i} ω_{j} {[covar ({\hat{σ}}_{x_{i}}, {\hat{σ}}_{x_{j}}) + σ_{x_{i}} σ_{x_{j}}] covar ({\hat{Z}}_{X_{i}}, {\hat{Z}}_{X_{j}}) + Z_{X_{i}} Z_{X_{j}} covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) \\ + covar ({\hat{μ}}_{X_{i}}, {\hat{μ}}_{X_{j}}) + [a_{i} a_{j} covar ({\hat{σ}}_{θ_{i}}, {\hat{σ}}_{θ_{j}}) + a_{i} a_{j} σ_{θ_{i}} σ_{θ_{j}}] covar ({\hat{Z}}_{θ_{i}}, {\hat{Z}}_{θ_{j}}) \\ + Z_{θ_{i}} Z_{θ_{j}} a_{i} a_{j} covar ({\hat{σ}}_{θ_{i}}, {\hat{σ}}_{θ_{j}}) + a_{i} a_{j} covar ({\hat{μ}}_{θ_{i}}, {\hat{μ}}_{θ_{j}}) - 2 a_{j} σ_{X_{i}} σ_{θ_{j}} covar ({\hat{Z}}_{X_{i}}, {\hat{Z}}_{θ_{j}}) \\ - 2 a_{j} Z_{X_{i}} Z_{θ_{j}} covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{θ_{j}}) - 2 a_{j} covar ({\hat{μ}}_{X_{i}}, {\hat{μ}}_{θ_{j}}) + 2 Z_{X_{j}} covar ({\hat{μ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) + 2 a_{i} a_{j} Z_{θ_{j}} covar ({\hat{μ}}_{θ_{i}}, {\hat{σ}}_{θ_{j}}) \\ - 2 a_{j} Z_{X_{j}} covar ({\hat{μ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) - 2 a_{j} Z_{θ_{j}} covar ({\hat{μ}}_{X_{i}}, {\hat{σ}}_{θ_{j}}) + 2 Z_{X_{j}} σ_{X_{i}} covar ({\hat{Z}}_{X_{i}}, {\hat{σ}}_{X_{j}}) + 2 a_{j} Z_{X_{j}} σ_{X_{i}} covar ({\hat{Z}}_{θ_{i}}, {\hat{σ}}_{θ_{j}}) \\ + 2 σ_{X_{i}} covar ({\hat{Z}}_{X_{i}}, {\hat{μ}}_{X_{j}}) + 2 a_{j} σ_{X_{i}} covar ({\hat{Z}}_{θ_{i}}, {\hat{μ}}_{θ_{j}}) - 2 σ_{X i} σ_{X_{j}} covar ({\hat{Z}}_{X_{i}}, {\hat{Z}}_{θ_{j}}) - 2 σ_{X_{i}} Z_{X_{j}} covar ({\hat{Z}}_{θ_{i}}, {\hat{σ}}_{X_{j}}) \\ - 2 σ_{X_{i}} Z_{θ_{j}} covar ({\hat{Z}}_{X_{i}}, {\hat{σ}}_{θ_{j}}) - 2 a_{j} σ_{X_{i}} covar ({\hat{Z}}_{X_{i}}, {\hat{μ}}_{θ_{j}}) - 2 σ_{X_{i}} covar ({\hat{Z}}_{θ_{i}}, {\hat{μ}}_{X_{j}})} . \end{array}

The difference between the variance of the correct estimator and the conventional approach is

\begin{align} var (\hat{D}) - var ({\hat{D}}_{conv}) = \sum_{i} \sum_{j} ω_{i} ω_{j} {covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) [covar ({\hat{Z}}_{θ_{i}}, {\hat{Z}}_{θ_{j}}) - 2 covar ({\hat{Z}}_{θ_{i}}, {\hat{Z}}_{X_{j}})] \\ + [Z_{θ_{i}} Z_{θ_{j}} - 2 Z_{θ_{i}} Z_{X_{j}}] covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) - covar ({\hat{μ}}_{X_{i}}, {\hat{μ}}_{X_{j}}) - a_{i} a_{j} covar ({\hat{σ}}_{θ_{i}}, {\hat{σ}}_{θ_{j}}) covar ({\hat{Z}}_{θ_{i}}, {\hat{Z}}_{θ_{j}}) \\ - a_{i} a_{j} Z_{θ_{i}} Z_{θ_{j}} covar ({\hat{σ}}_{θ_{i}}, {\hat{σ}}_{θ_{j}}) - a_{i} a_{j} covar ({\hat{μ}}_{θ_{i}}, {\hat{μ}}_{θ_{j}}) + 2 a_{j} Z_{X_{i}} Z_{θ_{j}} covar ({\hat{σ}}_{X_{i}}, {\hat{σ}}_{θ_{j}}) + 2 a_{j} covar ({\hat{μ}}_{X_{i}}, {\hat{μ}}_{θ_{j}}) \\ - 2 Z_{X_{j}} covar ({\hat{μ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) - 2 a_{i} a_{j} Z_{θ_{j}} covar ({\hat{μ}}_{θ_{i}}, {\hat{σ}}_{θ_{j}}) + 2 a_{j} Z_{X_{j}} covar ({\hat{μ}}_{X_{i}}, {\hat{σ}}_{X_{j}}) + 2 a_{j} Z_{θ_{j}} covar ({\hat{μ}}_{X_{i}}, {\hat{σ}}_{θ_{j}}) \\ - 2 a_{j} Z_{X_{j}} σ_{X_{i}} covar ({\hat{Z}}_{θ_{i}}, {\hat{σ}}_{θ_{j}}) - 2 σ_{X_{i}} covar ({\hat{Z}}_{X_{i}}, {\hat{μ}}_{X_{j}}) - 2 a_{j} σ_{X_{i}} covar ({\hat{Z}}_{θ_{i}}, {\hat{μ}}_{θ_{j}}) + 2 σ_{X_{i}} Z_{θ_{j}} covar ({\hat{Z}}_{X_{i}}, {\hat{σ}}_{θ_{j}}) \\ + 2 a_{j} σ_{X_{i}} covar ({\hat{Z}}_{X_{i}}, {\hat{μ}}_{θ_{j}}) + 2 σ_{X_{i}} covar ({\hat{Z}}_{θ_{i}}, {\hat{μ}}_{X_{j}})} . \end{align}

(18)

The simplified quadratic formula or correction-term method approximates the difference in Equation (18) as a quadratic function of the statistic of interest, with coefficients that can be pre-calculated to apply to any statistic of the same type and any scenario or type of score comparison (here, subpopulation invariance). That is,

Correction = A_{Y} Y^{2} + A_{Y, X} Y X + B_{Y} Y + B_{X} X + C .

(19)

To approximate Equation (18) in quadratic form, covariances involving the standardized scores are dropped because these terms are not invariant across all statistics. Furthermore, composite scores are approximated as if the composite scores were linked directly instead of the subscales being linked.

The parameters for the resulting correction term are

\begin{array}{l} A_{Y} & = d_{1}, \\ A_{X, Y} & = d_{2}, \\ B_{Y} & = d_{3} - m μ_{X} (2 d_{1} + d_{2}), \\ B_{X} & = d_{4} - m μ_{X} d_{2}, \\ C & = 2 a covar ({\hat{μ}}_{X}, {\hat{μ}}_{θ}) - var ({\hat{μ}}_{X}) - a^{2} var ({\hat{μ}}_{θ}) + m μ_{X}^{2} (d_{1} + d_{2}) - m μ_{X} (d_{3} + d_{4}), \end{array}

where

\begin{array}{l} d_{1} & = \frac{1}{σ_{X}^{2}} [var ({\hat{σ}}_{X}) - a^{2} var ({\hat{σ}}_{θ})], \\ d_{2} & = \frac{2}{σ_{X}^{2}} [a covar ({\hat{σ}}_{X}, {\hat{σ}}_{θ}) - var ({\hat{σ}}_{X})], \\ d_{3} & = \frac{2}{σ_{X}} [a covar ({\hat{μ}}_{X}, {\hat{σ}}_{θ}) - a^{2} covar ({\hat{μ}}_{θ}, {\hat{σ}}_{θ})], \\ d_{4} & = \frac{2}{σ_{X}} [a covar ({\hat{μ}}_{θ}, {\hat{σ}}_{X}) - covar ({\hat{μ}}_{X}, {\hat{σ}}_{X})] . \end{array}

The correction-term method in full is

{var}_{Correction - term} (\hat{D}) = {var}_{Conventional} (\hat{D}) + Correction,

(20)

where

{var}_{Conventional} (\hat{D})

is calculated with Equation (6) if samples are independent and Equation (7) if samples are dependent.

For reference, we define the complete correction-term method as a variant of the correction-term approach that retains the covariance terms involving ${\hat{Z}}_{X}$ and ${\hat{Z}}_{θ}$ . When these terms are included, the coefficients become specific to the focal statistic, underming the main advantage of the correction-term method. That is, to enable analysts to incorporate linking variance without requiring estimation of additional (co)variances. Nevertheless, the complete correction-term method can be used as a reference to evaluate the accuracy of the correction-term method approximations.

The complete correction-term method is defined with the same expressions as the correction-term method except that,

d_{3} = \frac{2}{σ_{X}} [a^{2} covar ({\hat{μ}}_{θ}, {\hat{σ}}_{θ}) - a covar ({\hat{μ}}_{X}, {\hat{σ}}_{θ})] + 2 covar ({\hat{Z}}_{X}, {\hat{σ}}_{θ}) - 2 a covar ({\hat{Z}}_{θ}, {\hat{σ}}_{θ}),

(21)

and

\begin{align} C = 2 a covar ({\hat{μ}}_{X}, {\hat{μ}}_{θ}) - var ({\hat{μ}}_{X}) - a^{2} var ({\hat{μ}}_{θ}) + m μ_{X}^{2} (d_{1} + d_{2}) - m μ_{X} (d_{3} + d_{4}) \\ + 2 a {\hat{σ}}_{X} covar ({\hat{Z}}_{X}, {\hat{μ}}_{θ}) + 2 {\hat{σ}}_{X} covar ({\hat{Z}}_{θ}, {\hat{μ}}_{X}) - 2 σ_{X} covar ({\hat{Z}}_{X}, {\hat{μ}}_{X}) - 2 a σ_{X} covar ({\hat{Z}}_{θ}, {\hat{μ}}_{θ}) . \end{align}

(22)

Simulation

The primary purpose of the simulation was to evaluate the statistical performance of the proposed and existing variance estimation methods in terms of bias, root mean square error (RMSE), and confidence interval coverage probability.

A secondary purpose was to evaluate the impact of the use of theoretically invalid methods. Specifically, (1) What is the impact of ignoring linking error by using a conventional variance estimator? (2) What is the impact of accounting for linking error but ignoring linking error dependencies, as in the linking component method? (3) What is the impact of assuming independence samples when samples are dependent?

The simulation followed a random-groups data collection design, in which a new and prior assessment were administered to randomly equivalent samples (Dorans, 2018). To induce non-normality, scores were generated using a linear regression model with dichotomous predictors, analogous to the statistical method used to generate the scores in large-scale assessments (Jewsbury, 2023).

Scores from the new administration were linearly transformed to the prior metric using Equations (1), (2), and (3). Subpopulation invariance was assessed by estimating the focal group mean difference across the two assessments with Equation (4), and several variance estimation methods were applied to estimate the associated standard error.

Each condition was replicated $R = 10,000$ times. The simulation was conducted in R 4.4.2 with base functions (e.g., rnorm; R Core Team, 2024).

Method

Conditions

The simulation varied the following conditions:

$•$ Focal group percentage sample size, $n \in {5 %, 10 %, 25 %, 50 %, 75 %, 90 %, 95 %}$ . This range was chosen based on the expectation that accounting for linking error becomes more important for groups that make up a larger proportion of the sample (Jewsbury, 2019).

$•$ Effect size of the subpopulation non-invariance, $d \in {0, . 05, . 10, . 25}$ . The effect size was varied to evaluate performance across different degrees of subpopulation non-invariance. The value $d = 0$ represents subpopulation invariance; values $d \in {. 05, . 10}$ reflect empirically observed values; and $d = . 25$ represents a larger effect than generally found empirically (Jewsbury et al., 2020; NCES, 2025a, 2025c, 2025b; Yin, Bezirhan, et al., 2023). Effect sizes were defined using the standard deviation of each sample individually (not pooled).

$•$ Independent versus dependent samples. The two samples were simulated as either dependent or independent. A central focus of this condition was the comparison between total-based methods, which account for dependencies, and component-based methods, which assume independence.

$•$ Normal versus non-normal scores. To evaluate robustness to score distribution, non-normal scores were simulated using a regression on dichotomous predictors, approximating latent regression models in large-scale assessments (Jewsbury, 2023; von Davier et al., 2006). Non-normal sores were generated using parameters from Jewsbury (2024) to facilitate comparison with prior work. For comparison, scores that were normally distributed conditional on the focal group were also simulated.

Data Generating Model

To simulate scores for two assessment samples, we generated the score vector $s_{i}^{⊤} = [\begin{matrix} s_{X i} & s_{θ i} \end{matrix}] \in R^{1 \times 2}$ for individual $i \in {1, \dots, N = 1000}$ from the model

s_{i} = [\begin{matrix} G_{X i}^{⊤} γ_{X} + τ_{X} \\ G_{Y i}^{⊤} γ_{Y} + τ_{Y} \end{matrix}] + ε_{i}, ε_{i} \sim N_{2} (0, Σ),

where

G_{X i}, G_{Y i} \in {0,1}^{10}

are group membership vectors for individual

i

in samples

X

and

Y

, respectively, and

γ_{X}, γ_{Y} \in R^{10}

are the corresponding regression coefficient vectors. The intercept vector is

τ^{⊤} = [\begin{matrix} τ_{X} & τ_{θ} \end{matrix}] = [\begin{matrix} 0 & 0.3 \end{matrix}]

. The residual covariance matrix is defined as

Σ = [\begin{matrix} σ_{X}^{2} & ρ σ_{X} σ_{Y} \\ ρ σ_{X} σ_{Y} & σ_{Y}^{2} \end{matrix}],

with the variances

σ_{X}^{2}

and

σ_{Y}^{2}

were set to ensure unit marginal variances of

s_{X}

and

s_{θ}

. The residual correlation

ρ

was 0.5 (dependent samples) or 0 (independent samples).

Letting $j \in {X, Y}$ index the assessment sample and $g \in {1, \dots, 10}$ index groups, group memberships were assigned independently across groups according to probabilities $p_{g}$ . In the dependent samples condition, the same group membership draw $q_{ijg} \sim Bernoulli (p_{g})$ was used for both samples; in the independent samples condition, independent draws were made for each sample. For non-focal groups $g \in {2, \dots, 10}$ , group sizes were drawn from $p_{g} \sim U (0.45, 0.55)$ and held constant across replications. The regression coefficients for these non-focal groups were generated under subpopulation invariance and held constant across replications, with $γ_{g X} = γ_{g Y} \sim N (0, η^{2})$ , where $η = 0$ (normal) or $η = 0.15$ (non-normal).

For the focal group $g = 1$ , the coefficients were generated under subpopulation non-invariance. Specifically, the effect size $d$ that took values in ${0, . 05, . 10, . 25, . 50}$ was split across samples such that $γ_{1 X} = - d / 2$ and $γ_{1 Y} = + d / 2$ . The focal group proportion $p_{1}$ varied across conditions, taking values in ${. 95, . 90, . 75, . 50, . 25, . 10, . 05}$ .

Linking coefficients were estimated with Equations (2) and (3) and applied to transform the sample $Y$ scores, $s_{θ}$ , using Equation (1). Subpopulation non-invariance of the focal group mean was estimated with Equation (4), and standard errors were obtained from multiple variance estimation methods.

Variance Estimation Methods

Two existing methods were implemented: (1) the conventional jackknife method that ignores linking error and (2) the linking component jackknife method that accounts for linking error but ignores linking error dependencies. In the independent samples condition, the conventional method was implemented with Equation (6) and the linking component method was implemented with Equations (6) and (8). In the dependent samples condition, the conventional method was implemented with Equation (7) and the linking component method was implemented with Equations (7) and (8).

The five novel methods that account for both linking error and linking error dependencies were also implemented. The total-based resampling method was implemented using the jackknife in Equation (14) for both independent and dependent sample conditions. The total-based equation was implemented using Equation (15). Although the cross-sample covariances in the total-based equation could have been set to zero in the independent samples condition, we estimated all covariances in both independent and dependent sample conditions to align with the total-based resampling method. Similarly, the component-based resampling method was implemented using the jackknife in Equation (16) and the component-based equation method using Equation (17) for both independent and dependent sample conditions.

In the independent samples condition, the correction-term method was implemented with Equations (6) and (20), and the cross-sample covariances were set to zero when computing the correction-term coefficients. In the dependent samples condition, the correction-term method was implemented with Equations (7) and (20).

Outcomes

We evaluated the percentage bias, percentage RMSE and 95% confidence interval coverage probability for all methods. The percentage bias was calculated as

Percentage bias = \frac{100}{\sqrt{var (\hat{y})}} \frac{1}{R} \sum_{r = 1}^{R} (\sqrt{{\hat{var}}_{r} (\hat{y})} - \sqrt{var (\hat{y})}),

(23)

where

\sqrt{{\hat{var}}_{r} (\hat{y})}

is the estimated standard error for the focal method in replication

r \in {1, \dots, R}

, and

\sqrt{var (\hat{y})}

is the actual standard deviation of the estimated mean difference across the

R

replications.

Percentage RMSE was calculated as

Percentage RMSE = \frac{100}{\sqrt{var (\hat{y})}} \sqrt{\frac{1}{R} \sum_{r = 1}^{R} {(\sqrt{{\hat{var}}_{r} (\hat{y})} - \sqrt{var (\hat{y})})}^{2}} .

(24)

The confidence interval coverage probability was calculated as

Coverage = \frac{100}{R} \sum_{r = 1}^{R} I ({\hat{y}}_{r} - t \sqrt{{\hat{var}}_{r} (\hat{y})} \leq y \leq {\hat{y}}_{r} + t \sqrt{{\hat{var}}_{r} (\hat{y})}),

(25)

where

I (\cdot)

is the indicator function, equal to 1 if the condition holds and 0 otherwise,

{\hat{y}}_{s}

is the estimated mean difference in replication

r

y

is the true mean difference, and

t

is the critical

t

-statistic corresponding to a 95% confidence level.

y

was calculated analytically from the data generating model parameters.

Results and Discussion

The pattern of results was consistent across score normality and degree of subpopulation non-invariance. Therefore, we illustrate the results with a representative condition: non-normal scores and degree of subpopulation non-invariance $d = . 10$ . A full table of results is provided in the Online Supplemental Materials.

Because the existing methods exhibited substantially higher bias, inflating the y-axis range, differences among the proposed methods are difficult to discern when plotted with the existing methods. For this reason, we include two figures: Figure 1 shows only the proposed methods; Figure 2 includes all methods.

Validity of Proposed Methods

The new total- and component-based methods performed well (Figure 1). As expected, total-based methods showed negligible bias and confidence interval coverage close to the nominal $95 %$ across all conditions. Component-based methods, which assume sample independence, also showed good bias and coverage when applied to independent samples. These methods slightly outperformed total-based methods in RMSE under independence, demonstrating the efficiency gain from correctly assuming independence. The correction-term method, which is more approximate, slightly underperformed the others in RMSE, particularly for dependent samples.

Notably, the corresponding resampling and equation pairs performed very similarly with largely overlapping results in the plots (Figure 1). Consequently, the equations can be used to understand and predict the behavior of the resampling methods.

Impact of the Use of Theoretically Invalid Methods

Methods that do not properly account for the linking error, linking error dependencies, or sample dependence yielded poor results (Figure 2). The conventional estimator, which ignores linking error, showed consistently high positive bias and inflated RMSE, especially when the focal group made up a larger portion of the sample.

The linking component estimator, which includes a linking error term but ignores its dependence with the conventional estimate, had even worse bias (Figure 2). Since this estimator adds a positive component to an already positively biased estimate (see Equation (8)), the bias is necessarily increased.

When applied to dependent samples, component-based methods (which assume independence) also showed non-negligible biases (Figure 1). This finding supports the recommendation to use total-based methods for dependent samples. When samples were independent, both total- and component-based methods were approximately unbiased, with the latter showing slightly lower RMSE as previously noted.

Practical Example

We conducted an empirical analysis using data from the 2021 Progress in International Reading Literacy Study (PIRLS; von Davier et al., 2023). PIRLS was selected because it represents one of the most methodologically complex assessments, involving sampling weights, multistage sampling, and plausible values. In addition, its international design enabled replication of the analysis across multiple country samples.

This analysis served two purposes. First, it illustrates the application of the proposed methods to real data and demonstrates how power to detect subpopulation non-invariance improves when standard errors are correctly estimated. Second, it enables an empirical evaluation of the methods by assessing the plausibility of the results and the consistency of the methods in a complex real-world assessment setting.

PIRLS 2021 was the first administration in which countries could choose to administer a computer-based version of the assessment, departing from the prior paper-based assessment. Countries opting into the digital assessment were also assessed with the paper assessment in a randomly equivalent sample, enabling the metrics of the new (digital) and prior (paper) assessments to be linked (Almaskut et al., 2023).

We deviated from the official PIRLS procedure, which linked a pooled sample across all the participating countries (Yin, Fishbein, et al., 2023), and instead performed the linking within each country separately. This allowed us to evaluate statistical power more directly by replicating equivalent analyses across multiple country samples.

Within each country, we examined subpopulation invariance of the mean reading score conditional on student’s responses to the yes-no question, “Do you have any of these things in your home? Your own computer or tablet” (see IEA, 2020 for the specific question). This item was chosen as a likely candidate to exhibit subpopulation non-invariance, and in directions that can be reasonably anticipated (i.e., better performance on digital for students with a computer or tablet, better performance on paper for students without).

The analysis was conducted in R 4.4.2 with base functions (R Core Team, 2024). The code is provided in the Online Supplemental Materials.

Method

Data

The R data for 2021 PIRLS was downloaded from the IEA Data Repository (Fishbein et al., 2024). We identified 26 pairs of matching paper and digital datasets. To simplify the results, we included only samples that represented entire countries, excluding Belgium (Flemish language) and Moscow City.¹ The paired data for the United States, the 27th participant in the digital administration, was absent from the IEA Data Repository, possibly due to a decision not to report digital scores (NCES, 2025e).

Sample sizes for the digital assessments ranged from 3,030 to 8,551, with an outlier of 27,448 (United Arab Emirates). Paper-based samples ranged from 835 to 3,207 participants. The paper and digital samples were designed to be randomly equivalent. In some countries the samples were dependent (Yin, Bezirhan, et al., 2023). However, it is unclear whether the sample dependence could be modeled because the jackknife sampling zones may not have been structured to reflect this dependency (Foy & Almaskut, 2023). Therefore, consistent with the approach used in the PIRLS technical documentation (Yin, Fishbein, et al., 2023), we treated the samples as independent for analysis.

Analysis

For each country, we first unapplied the official PIRLS linking transformation from the digital plausible values using the published linking coefficients (Yin, Bezirhan, et al., 2023). We then estimated country-specific linking coefficients using Equations (2) and (3), and applied them to the digital plausible values via Equation (1).

PIRLS employs standard large-scale assessment methodology, including case weights, replicate weights, and plausible values. Due to its complexity, we refer readers to detailed treatments elsewhere of the basic methodology (e.g., Jewsbury et al., 2024; M. von Davier et al., 2023). Only key details and extensions of the basic methodology will be described here.

Estimates of population parameters, including the population mean and standard deviation estimates in Equations (2) and (3), were obtained using the standard methodology used in PIRLS (M. von Davier et al., 2023). That is, the parameters were estimated using the case weights² and for each plausible value, with results averaged across the plausible values.

Following official procedures, we constructed 250 replicate weights per country for the paired jackknife method (Foy & Almaskut, 2023). Within-imputation variances were estimated using these replicate weights for each plausible value and averaged. Between-imputation variances were estimated using standard procedures, and total variance was obtained via Rubin’s rules (Rubin, 1987).

Degrees of freedom for the subpopulation mean estimates $\hat{X}$ and $\hat{Y}$ were calculated using the method of Barnard and Rubin (1999), and where the maximum number of jackknife zones was used as the complete-case degrees of freedom. For the difference in means Equation (4), degrees of freedom were estimated via the Welch–Satterthwaite approximation (Welch, 1938). These were used to compute the critical $t$ -value for constructing confidence intervals.

We estimated the subpopulation non-invariance of the mean reading score within levels of the variable “[do you have in your home] your own computer or tablet”³ using Equation (4). The levels of the variables were labeled yes, no, and missing (specifically, omitted or invalid). Standard errors, confidence intervals, and $t$ -tests for the null hypothesis $\hat{D} = 0$ were computed using multiple variance estimation methods.

Variance estimation methods. We implemented the two existing and five novel variance estimation methods previously introduced and evaluated in the simulation.

Applying the resampling methods (total and component resampling) was complex due to plausible values, which requires two variance estimators: one for within-imputation variance and one for between-imputation variance (Rubin, 1987). We adapted both the within- and between-imputation variance estimators to account for linking variance under the total and component approaches. Full details on these adaptations are provided in the Appendix.

Equation-based methods (total equation, component equation, and correction-term) were implementing by plugging-in estimated parameters and (co)variances into Equations (15), (17), and (20), respectively. The parameters and (co)variances were estimated using standard PIRLS estimation procedures (Foy & Almaskut, 2023). For the total equation and correction-term methods, the cross-sample covariances were set to zero. Consequently all methods with the exception of the total resampling method assumed sample independence. We also implemented the complete correction-term method using Equations (21) and (22) as a reference for the correction-term method.

Results and Discussion

Impact of Accounting for Linking Variance

Given the scope of the analysis—applying seven variance estimation methods across three subpopulations in 24 countries—the full set of results is provided in the Online Supplemental Materials. In the main text, Figure 3 illustrates the increase in power obtained by using valid standard error estimators versus ignoring linking error when testing for subpopulation invariance. To reduce visual complexity, only the component-based methods are shown alongside the conventional method; these were selected based on their superior performance under sample independence in the simulation study. Figure 4 displays a heatmap of significance probability values across countries and subpopulations to assess the empirical consistency of all of the proposed methods.

Consistent with the simulation results, the conventional method produced larger confidence intervals relative to the component resampling and component equation methods. This overestimation was more pronounced in the larger subgroup, Own computer at home: yes (mean proportion = 63.8%, SD = 12.8%), compared to the no subgroup (mean proportion = 33.7%, SD = 12.4%). As expected, all three methods yielded similar estimates for the very small missing subgroup (mean proportion = 1.4%, SD = 1.0%).

While truth is unknown with real data, we can consider the results in terms of plausibility. According to the law of total variance, the variance of $\hat{D}$ across countries should equal the sum of between-country variance in the true values $D$ ⁴ and the mean of the within-country estimator variances $var (\hat{D})$ . For the yes group, the observed variance of $\hat{D}$ across countries was 8.6, whereas the mean of the conventionally estimated variances was 33.7, implying an impossible negative between-country variance in $D$ . In contrast, the component resampling and equation methods produced mean estimated variances of 6.5 and 6.6, respectively, implying a plausible and positive between-country variance.

Figure 3 visually illustrates this discrepancy for the yes group. With correct standard errors, the cross-country spread in $\hat{D}$ should be comparable or larger than the width of the confidence intervals. Instead, the point estimates cluster near the same value, with the average $\hat{D} = - 3.0$ close to the center of the intervals in most countries. Moreover, despite wide intervals spanning both positive and negative values, $\hat{D}$ is negative in 22 out of 24 countries, suggesting greater precision than the conventional standard errors reflect.

The results indicate a meaningful increase in power when standard errors are correctly estimated (Figure 3). With the conventional method, none of the yes group comparisons reached statistical significance. However, using the component-based methods, $\hat{D}$ was significant in six countries. All significant effects were negative, indicating that students with their own computer at home tended to perform better on the digital assessment relative to the paper-based version.

For the no group, one country showed significant $\hat{D}$ values across all three methods, and an additional country was significant only when using component-based methods (Figure 3). All significant estimates were positive, suggesting better performance on the paper assessment among students without their own computer.

For the missing subgroup, two countries showed significant results across all three methods, and a third country was significant only under the component-based methods. Again, all significant estimates were positive, suggesting students who did not answer the question performed better on the paper version. Potentially, students who struggle with cognitive assessments administered in digital format may also struggle to respond to contextual questions administered in digital format.

Consistency of the Proposed Methods

Figure 4 displays a heatmap of significance probability values, illustrating the empirical consistency of the proposed methods. The four total and component-based methods produced identical significance outcomes in all but one instance (see Online Supplemental Materials). The only exception was a single result from the total resampling method, and some differences were expected as all other methods assumed sample independence. However, even in that case, the difference in $p$ -values was minimal (.005). In particular, the consistency between the resampling and equation methods, despite differing theoretical foundations, provides support for the methods in real data.

The correction-term method showed meaningful discrepancies from the other methods (Figure 4), and produced negative variance estimates for 11 countries in the yes group and 1 country in the no group. The poor performance can be attributed to the additional approximations of the correction-term approach where statistic-specific covariances are dropped, as the complete correction term method including the statistic-specific covariances was highly aligned with the other methods.

The underperformance of the correction-term method was more pronounced for the yes group, likely because of the magnitude of the correction. For the United Arab Emirates yes group, the conventional estimate of $var (\hat{D})$ was 76.2, while the component resampling estimate was 3.8. For the correction-term method to be consistent with the component resampling method, the correction must be $(76.2 - 3.8) / 3.8 = 19.1$ times larger than the component resampling method estimate. Consequently, only a small error in estimating the correction translates into a large error on the variance estimate.

Despite the apparent sensitivity of the correction-term method to approximations, the degree of underperformance of the method was not observed in the simulation, even with large correction magnitudes (e.g., 90% and 95% focal group conditions). The discrepancy may reflect the limited degrees of freedom of the PIRLS variance estimators. The paired jackknife procedure applied to imputed data has degrees of freedom lower than the number of jackknife sampling zones (Barnard & Rubin, 1999), which can be as few as 11 (Foy & Almaskut, 2023). The between-imputation variance estimator with five plausible values has only four degrees of freedom (Casella & Berger, 2002).

Discussion

This paper introduced five new methods for evaluating subpopulation invariance within the broader taxonomy of linking variance methods proposed by Jewsbury (2024). These methods were validated through simulation studies, where they demonstrated substantially improved statistical properties compared to existing approaches. In the context of subpopulation invariance, properly accounting for linking variance leads to reduced standard errors and increased statistical power. This increase in power was illustrated through an empirical application using data from the 2021 PIRLS assessment.

In random-groups designs, population invariance is ensured by construction via equating the marginal score distributions across groups (see Equation (1)). Accounting for linking variance in the estimation of standard errors enhances power to detect subpopulation non-invariance by correctly conditioning on the fact that population-level invariance has already been imposed by design (Equation (15)).

The results underscore the critical impact of correct standard error estimation. Using conventional methods that ignore linking variance led to dramatic overestimation of standard errors, resulting in severe loss of power, as seen in both simulation and the PIRLS example. Methods that include linking error but ignore its dependencies with other error sources, as most existing linking methods do (e.g., Johnson, 1998; Kawahashi et al., 2024; Li, 2022; Mazzeo et al., 2018), only exacerbated the overestimation. Similarly, applying methods that assume sample independence to dependent samples (e.g., using component-based methods where total-based methods are appropriate) resulted in notable bias.

Across both the simulation and the empirical analysis, resampling-based and equation-based methods yielded highly consistent results. Given their distinct theoretical foundations, this agreement offers additional validity evidence. Resampling methods, while more computationally intensive, are more flexible and can in principle accommodate non-linear linking functions (e.g., equipercentile) or more complex statistics (e.g., the standardized Root Mean Square Difference from Dorans & Holland, 2000). In contrast, equation-based methods offer greater transparency and interpretability, but are currently limited to linear linking and simpler summary statistics (e.g., means, standard deviations, and percentiles).

The distinction between our total-based and component-based methods centers on the treatment of sample dependence. Total-based methods explicitly accommodate dependence, while component-based methods assume independence. This was reflected in the simulation results: total-based methods outperformed under dependence, whereas component-based methods were slightly more efficient (lower RMSE) under independence.

One unexpected finding was the relative underperformance of the correction-term method, particularly in the empirical analysis. Although this method is intended to offer a computationally efficient approximation, prior work (Jewsbury, 2024) found that correction-term approaches performs comparably to less approximate methods in the context of trend comparisons. The discrepancy may stem from the relative magnitude of the correction term: in trend analyses, the correction tends to be positive and can be small relative to the total variance, whereas in subpopulation invariance analyses, the correction is negative and can exceed the total variance by several multiples (see e.g., Figures 2 and 3). Under such conditions, even modest proportional error in estimating the correction can lead to substantial mis-estimation of the standard error. Accordingly, use of the correction-term method should be accompanied by sensitivity checks against more robust methods.

As demonstrated in the practical example, the proposed methods are applicable even in the highly complex context of large-scale educational assessments. The PIRLS example involved case weights, multistage sampling, and plausible values. In particular, plausible values require that variance be decomposed into within-imputation and between-imputation variances (Jewsbury et al., 2024). Both PIRLS’ paired jackknife estimator for the within-imputation variance and the Monte Carlo estimator for the between-imputation variance were extended to account for linking variance under both the total and component resampling approaches (see Appendix).

While this paper focuses on linking with random-groups data collection designs (Dorans, 2018), other types of linking exist. The methods can also apply to pseudo-equivalent groups linking after the adjustment for non-equivalence has been applied (e.g., Haberman, 2015). However, linking variance with a non-equivalent anchor test design (NEAT) is more distinct (Robitzsch & Lüdtke, 2019; 2024; Xu & von Davier, 2010). For NEAT designs, the double jackknife may be the most analogous to the methods considered here, especially the total resampling method (Haberman et al., 2009).

In sum, this paper demonstrates the critical importance of properly accounting for linking variance when evaluating subpopulation invariance. We proposed several new methods that address this issue, validated them through simulation, and illustrated their utility through an empirical application.

Supplemental Material

Supplemental Material - Standard Error Estimation for Subpopulation Non-invariance

Supplemental Material for Standard Error Estimation for Subpopulation Non-invariance by Paul A. Jewsbury in Applied Psychological Measurement

Footnotes

ORCID iD

Paul A. Jewsbury

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The study was funded by ETS.

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.

Data Availability Statement

The PIRLS data used in the paper is available from the IEA Data Repository, .

Supplemental Material

Supplemental material for this article is available online.

Notes

Appendix

Applying the methods to large-scale assessments such as PIRLS involving case weights, complex samples, and plausible values adds additional complexity. The proposed equation-based estimators (total equation, component equation, and correction term) can be estimated by applying the standard procedures recommended by the large-scale assessment to estimate each (co)variance in the equation (Jewsbury et al., 2024). However, applying our proposed resampling methods is less straightforward. In this Appendix we describe how the resampling-based methods can be applied through extending the existing variance estimators used in large-scale assessments.

Large-scale assessments using plausible values follow Rubin’s rules to estimate parameters and their variances (Rubin, 1987). That is, (26)

\bar{Q} = \frac{1}{M} \sum_{l = 1}^{M} {\hat{Q}}^{(l)},

(27)

var (\bar{Q}) = \bar{U} + (1 + \frac{1}{M}) B,

where

\bar{Q}

is the statistic of interest,

{\hat{Q}}^{(l)}

is the statistic of interest based on plausible value

l \in {1, \dots, M}

, and (28)

\bar{U} = \frac{1}{M} \sum_{l = 1}^{M} U^{(l)},

(29)

B = \frac{1}{M - 1} \sum_{l = 1}^{M} {({\hat{Q}}^{(l)} - \bar{Q})}^{2} .

In the context of subpopulation invariance, the statistic of interest $Q$ is defined as the difference between two subpopulation summary statistics, estimated from two linked samples of scores (see Equation (4)). Specifically, (30)

{\hat{Q}}^{(l)} = {\hat{D}}^{(l)} = f (v_{X}^{(l)}, w_{X}) - f (v_{Y}^{(l)}, w_{Y}),

where

f (\cdot, \cdot)

denotes a weighted summary function (e.g., a weighted mean), with the first argument representing the scores and the second the corresponding weights. Here,

v_{X}^{(l)}

and

v_{Y}^{(l)}

are the

l

th sets of plausible values for the prior and new assessment samples, respectively. The new sample scores

v_{Y}^{(l)}

have been transformed to the metric of the prior assessment. The weights

w_{X}

and

w_{Y}

are the sampling weights for the prior and new assessment samples, respectively.

In the practical example, we adapt PIRLS’ variance estimators to incorporate linking variance. As is typical in large-scale assessments, PIRLS uses a resampling-based method to estimate the within-imputation variance. These estimators typically take the form: (31)

U^{(l)} = c \sum_{r = 1}^{R} {(δ_{U}^{Method} (l, r) - {\bar{δ}}_{U}^{Method} (l))}^{2},

where

c

is a design-specific scaling constant,

R

is the number of replicates,

δ_{U}^{Method} (l, r)

is the estimate of the statistic of interest based on the

l

th plausible value and the

r

th replicate weight, and

{\bar{δ}}_{U}^{Method} (l)

denotes the average of these estimates across replicates for the

l

th plausible value.

In PIRLS, $c = \frac{1}{2}$ . Note that large-scale assessments may instead use the estimate based on the original sampling weights ${\hat{Q}}^{(l)}$ rather than the average across replicates ${\bar{δ}}_{U} (l)$ in Equation (31). In this paper, we consistently use the average.

To incorporate linking variance into the within-imputation estimator (Equation (31)), we modify how the replicate estimate $δ$ is calculated. The superscript Method in $δ_{U}^{Method}$ indicates that the specific calculation depends on the approach used (total or component-based).

To fully account for linking variance, the between-imputation variance estimator must also be adapted. To express this, we re-write the standard formula as (32)

B = \frac{1}{M - 1} \sum_{l = 1}^{M} {(δ_{B}^{Method} (l) - {\bar{δ}}_{B}^{Method})}^{2},

where

δ_{B}^{Method} (l)

denotes the estimate of the statistic of interest based on the

l

th plausible value and the sampling weights, and

{\bar{δ}}_{B}^{Method}

denotes the average of these estimates across plausible values.

References

Almaskut

LaRoche

Foy

(2023). Sample design in PIRLS 2021. In von Davier

Mullis

I. V. S.

Fishbein

Foy

(Eds.), Methods and procedures: PIRLS 2021 technical report (pp. 3.1–3.32): TIMSS & PIRLS International Study Center. https://doi.org/10.6017/lse.tpisc.tr2103.kb9560

Anselmi

Robusto

Cristante

(2023). Enhancing computerized adaptive testing with batteries of unidimensional tests. Applied Psychological Measurement, 47(3), 167–182. https://doi.org/10.1177/01466216231165301

Arel-Bundock

Enevoldsen

Yetman

(2018). countrycode: An R package to convert country names and country codes. Journal of Open Source Software, 3(28), 848. https://doi.org/10.21105/joss.00848

Arslan

Lehman

Tenison

Sparks

J. R.

López

A. A.

Zapata-Rivera

(2024). Opportunities and challenges of using generative AI to personalize educational assessment. Frontiers in Artificial Intelligence, 7, Article 1460651. https://doi.org/10.3389/frai.2024.1460651

Barnard

Rubin

D. B.

(1999). Small sample degrees of freedom with multiple imputation. Biometrika, 86, 948–955. https://doi.org/10.1093/biomet/86.4.948

Battauz

Leôncio

(2023). A likelihood approach to item response theory equating of multiple forms. Applied Psychological Measurement, 47(3), 200–220. https://doi.org/10.1177/01466216231151702

Brennan

R. L.

(2008). A discussion of population invariance. Applied Psychological Measurement, 32(1), 102–114. https://doi.org/10.1177/0146621607311582

Bulut

Beiting-Parrish

Casabianca

J. M.

Slater

S. C.

Jiao

Song

Morilova

(2024). The rise of artificial intelligence in educational measurement: Opportunities and ethical challenges. Chinese/English Journal of Educational Measurement and Evaluation, 5(3), 3. https://doi.org/10.59863/MIQL7785

Casella

Berger

R. L.

(2002). Statistical inference (2nd ed.): Duxbury.

10.

Castellano

K. E.

Johnson

M. S.

Lawless

(2024). Are large-scale test scores comparable for at-home versus test center testing? Applied Psychological Measurement, 48(6), 295–299. https://doi.org/10.1177/01466216241253795

11.

Castellano

K. E.

Sinharay

Hao

(2023). An investigation into the impact of test session disruptions for at-home test administrations. Applied Psychological Measurement, 47(1), 76–82. https://doi.org/10.1177/01466216221128011

12.

Chao

H. Y.

Chen

J. H.

(2023). Controlling the minimum item exposure rate in computerized adaptive testing: A two-stage sympson–hetter procedure. Applied Psychological Measurement, 47(7-8), 460–477. https://doi.org/10.1177/01466216231209756

13.

Deng

Rios

J. A.

(2022). Investigating the effect of differential rapid guessing on population invariance in equating. Applied Psychological Measurement, 46(7), 589–604. https://doi.org/10.1177/01466216221108991

14.

Dorans

N. J.

(2004). Using subpopulation invariance to assess test score equity. Journal of Educational Measurement, 41(1), 43–68. https://doi.org/10.1111/j.1745-3984.2004.tb01158.x

15.

Dorans

N. J.

(2018). Scores, scales, and score linking. In Irwing

Booth

Hughes

D. J.

(Eds). The wiley handbook of psychometric testing: A multidisciplinary reference on survey, scale and test development. (V.II, pp. 573–606): Wiley. https://doi.org/10.1002/9781118489772.ch19

16.

Dorans

N. J.

Holland

P. W.

(2000). Population invariance and the equatability of tests: Basic theory and the linear case. Journal of Educational Measurement, 37(4), 281–306. https://doi.org/10.1111/j.1745-3984.2000.tb01088.x

17.

Dorans

N. J.

Liu

Hammond

(2008). Anchor test type and population invariance: An exploration across subpopulations and test administrations. Applied Psychological Measurement, 32(1), 81–97. https://doi.org/10.1177/0146621607311580

18.

Fishbein

Martin

M. O.

Mullis

I. V.

Foy

(2018). The TIMSS 2019 item equivalence study: Examining mode effects for computer-based assessment and implications for measuring trends. Large-scale Assessments in Education, 6(1), 1–23. https://doi.org/10.1186/s40536-018-0064-z

19.

Fishbein

Yin

Foy

(2024). PIRLS 2021 user guide for the international database (2nd ed). TIMSS & PIRLS International Study Center. https://pirls2021.org/data

20.

Foy

Almaskut

(2023). Estimating standard errors in the PIRLS 2021 results. In von Davier

Mullis

I. V. S.

Fishbein

Foy

(Eds.), Methods and procedures: PIRLS 2021 technical report (pp. 13.1–13.31): TIMSS & PIRLS International Study Center. https://doi.org/10.6017/lse.tpisc.tr2103.kb1576

21.

Haberman

S. J.

(2015). Pseudo-equivalent groups and linking. Journal of Educational and Behavioral Statistics, 40(3), 254–273. https://doi.org/10.3102/1076998615574772

22.

Haberman

S. J.

Lee

Y.-H.

Qian

(2009). Jackknifing techniques for evaluation of equating accuracy (Research Report No. RR-09-39): Educational Testing Service. https://doi.org/10.1002/j.2333-8504.2009.tb02196.x

23.

Hao

von Davier

A. A.

Yaneva

Lottridge

von Davier

Harris

D. J.

(2024). Transforming assessment: The impacts and implications of large language models and generative AI. Educational Measurement: Issues and Practice, 43(2), 16–29. https://doi.org/10.1111/emip.12602

24.

Huggins

A. C.

(2014). The effect of differential item functioning in anchor items on population invariance of equating. Educational and Psychological Measurement, 74(4), 627–658. https://doi.org/10.1177/0013164413506222

25.

Huggins

A. C.

Elbaum

(2013). Test accommodations and equating invariance on a fifth-grade science exam. Educational Assessment, 18(1), 49–72. https://doi.org/10.1080/10627197.2013.761536

26.

IEA . (2020). Progress in international reading literacy study: Student questionnaire. Retrieved from. https://pirls2021.org/wp-content/uploads/2023/05/P21_StudentQuestionnaire.pdf

27.

Ikeda

Clark

Papageorgiou

Ohta

Blackhurst

Bruce

(2025). Aligning scores of language proficiency tests: A score concordance study between ielts academic and toefl iBT (research report No. RR-25-02): Educational Testing Service.

28.

Inal

Arıkan

Ç. A.

(2017). An investigation of group invariance in test equating according to gender. Journal of Measurement and Evaluation in Education and Psychology, 8(1), 128–145. https://doi.org/10.21031/epod.301326

29.

Jewsbury

P. A.

(2019). Error variance in common population linking bridge studies (Research Report No. RR-19-42): Educational Testing Service.

30.

Jewsbury

P. A.

(2023). Educational surveys: Methodological foundations. In Tierney

Rizvi

Ercikan

(Eds.), International encyclopedia of education: Quantitative research/educational measurement (4th ed.). Elsevier.

31.

Jewsbury

P. A.

(2024). Generally applicable variance estimation methods for common-population linking. Journal of Educational and Behavioral Statistics. https://doi.org/10.3102/10769986241263976

32.

Jewsbury

P. A.

Finnegan

R. J.

Jia

Rust

K. F.

Burg

(2020). 2017 NAEP transition to digitally-based assessment in mathematics and reading at grades 4 and 8: Mode evaluation study. Retrieved from. https://nces.ed.gov/nationsreportcard/subject/publications/main2020/pdf/transitional_whitepaper.pdf

33.

Jewsbury

P. A.

Jia

Gonzalez

(2024). Considerations for the use of plausible values in large-scale assessments. Large-scale Assessments in Education, 12(24), 1–27. https://doi.org/10.1186/s40536-024-00213-y

34.

Jewsbury

P. A.

(2025). Linking error on achievement levels accounting for dependencies and complex sampling. Journal of Educational Measurement. https://doi.org/10.1111/jedm.12439

35.

Johnson

E. G.

(1998). Linking the national assessment of educational progress (NAEP) and the third international mathematics and science study (TIMSS): A technical report (NCES 98-499). Government Printing Office.

36.

Johnson

M. S.

(2025). Responsible AI for measurement and learning: Principles and practices. ETS. (Research Report No. RR-25-03). https://www.ets.org/Media/Research/pdf/RR-25-03.pdf

37.

Jones

Tong

Liu

Borglum

Primoli

(2022). Score comparability between online proctored and in-person credentialing exams. Journal of Educational Measurement, 59(2), 180–207. https://doi.org/10.1111/jedm.12320

38.

Kawahashi

Miyazaki

Kubo

(2024). Asymptotic standard errors of equating coefficients using second-order delta method for non-parametric ability distribution [Poster presentation]. International Meeting of the Psychometric Society, Prague, Czechia.

39.

Kim

Walker

M. E.

McHale

(2010). Comparisons among designs for equating mixed-format tests in large-scale assessments. Journal of Educational Measurement, 47(1), 36–53. https://doi.org/10.1111/j.1745-3984.2009.00098.x

40.

Kolen

M. J.

Brennan

R. L.

(2014). Test equating, scaling, and linking: Methods and practices (3rd ed.): Springer-Verlag.

41.

Leôncio

Wiberg

Battauz

(2023). Evaluating equating transformations in IRT observed-score and kernel equating methods. Applied Psychological Measurement, 47(2), 123–140. https://doi.org/10.1177/01466216221124087

42.

(2022). Assessing the impact of equating error on group means and group mean differences. Journal of Educational Measurement, 59(1), 62–79. https://doi.org/10.1111/jedm.12311

43.

Lim

Sireci

S. G.

(2017). Linking TIMSS and NAEP assessments to evaluate international trends in achievement. Education Policy Analysis Archives, 25(11), 1–21. https://doi.org/10.14507/epaa.25.2682

44.

Liu

Jurich

(2023). Outlier detection using t-test in Rasch IRT equating under NEAT design. Applied Psychological Measurement, 47(1), 34–47. https://doi.org/10.1177/01466216221124045

45.

Liu

Dorans

N. J.

(2013). Assessing a critical aspect of construct continuity when test specifications change or test forms deviate from specifications. Educational Measurement: Issues and Practice, 32(1), 15–22. https://doi.org/10.1111/emip.12001

46.

Mazzeo

Donoghue

J. R.

Liu

(2018). Estimating standard errors for NAEP that incorporate random-groups linking error for the transition from paper-based to digital-based assessments [Paper presentation]. In Annual Meeting of the National Council on Measurement in Education. New York, NY, United States.

47.

Moses

(2022). Linking and comparability across conditions of measurement: Established frameworks and proposed updates. Journal of Educational Measurement, 59(2), 231–250. https://doi.org/10.1111/jedm.12322

48.

Moses

Zhang

(2011). Standard errors of equating differences: Prior developments, extensions, and simulations. Journal of Educational and Behavioral Statistics, 36(6), 779–803. https://doi.org/10.3102/1076998610396892

49.

NCES . (2025a). 2018 NAEP transition to DBA and model evaluation for the U.S. History, geography, and civics assessments at grade 8. Retrieved from. https://www.nationsreportcard.gov/civics/supporting_files/2018_naep_social_sciences_dba_transition.docx

50.

NCES . (2025b). 2019 NAEP transition to DBA and mode evaluation for the mathematics and reading assessments at grade 12. Retrieved from. https://www.nationsreportcard.gov/mathematics/supportive_files/2019_rm_g12_dba_transition.pdf

51.

NCES . (2025c). 2019 NAEP transition to DBA and mode evaluation for the science assessments at grades 4, 8, and 12. Retrieved from. https://www.nationsreportcard.gov/science/supporting_files/NAEP_2019_Science_Transition_to_DBA_and_Mode_Evaluation_Study.docx

52.

NCES . (2025d). NAEP technical documentation: Definition of composite scales. Retrieved from. https://nces.ed.gov/nationsreportcard/tdw/analysis/scaling_determination_composite.aspx

53.

NCES . (2025e). Progress in international reading literacy study: Methodology and technical notes. Retrieved from. https://nces.ed.gov/surveys/pirls/pirls2021/technotes.asp

54.

Petersen

N. S.

(2008). A discussion of population invariance of equating. Applied Psychological Measurement, 32(1), 98–101. https://doi.org/10.1177/0146621607311581

55.

R Core Team . (2024). R: A language and environment for statistical computing: R Foundation for Statistical Computing.

56.

Reardon

S. F.

Kalogrides

A. D.

(2021). Validation methods for aggregate-level test scale linking: A case study mapping school district test score distributions to a common scale. Journal of Educational and Behavioral Statistics, 46(2), 138–167. https://doi.org/10.3102/1076998619874089

57.

Robitzsch

Lüdtke

(2019). Linking errors in international large-scale assessments: Calculation of standard errors for trend estimation. Assessment in Education: Principles, Policy & Practice, 26(4), 444–465. https://doi.org/10.1080/0969594x.2018.1433633

58.

Robitzsch

Lüdtke

(2024). An examination of the linking error currently used in PISA. Measurement: Interdisciplinary Research and Perspectives, 22(1), 61–77. https://doi.org/10.1080/15366367.2023.2198915

59.

Rubin

D. B.

(1987). Multiple imputation for nonresponse in surveys. John Wiley & Sons.

60.

Sinharay

Haberman

S. J.

(2014). An empirical investigation of population invariance in the value of subscores. International Journal of Testing, 14(1), 22–48. https://doi.org/10.1080/15305058.2013.822712

61.

Sinharay

Johnson

M. S.

(2024). Computation and accuracy evaluation of comparable scores on culturally responsive assessments. Journal of Educational Measurement, 61(1), 5–46. https://doi.org/10.1111/jedm.12381

62.

Tian

Choi

(2023). The impact of item model parameter variations on person parameter estimation in computerized adaptive testing with automatically generated items. Applied Psychological Measurement, 47(4), 275–290. https://doi.org/10.1177/01466216231165313

63.

von Davier

A. A.

(2006). Population invariance of test equating and linking: Theory extension and applications across exams (Research Report No. RR-06-31). Educational Testing Service. https://doi.org/10.1002/j.2333-8504.2006.tb02037.x

64.

von Davier

A. A.

(2013). Observed-score equating: An overview. Psychometrika, 78(4), 605–623. https://doi.org/10.1007/s11336-013-9319-3

65.

von Davier

Mullis

I. V. S.

Fishbein

Foy

(Eds.), (2023). Methods and procedures: PIRLS 2021 technical report. TIMSS & PIRLS International Study Center. https://pirls2021.org/methods

66.

von Davier

Sinharay

Oranje

Beaton

(2006). The statistical procedures used in National Assessment of Educational Progress: Recent developments and future directions. In Rao

C. R.

Sinharay

(Eds.) Handbook of statistics, (pp. 1039–1055): Elsevier.

67.

Wallmark

Josefsson

Wiberg

(2023). Efficiency analysis of item response theory Kernel equating for mixed-format tests. Applied Psychological Measurement, 47(7-8), 496–512. https://doi.org/10.1177/01466216231209757

68.

Weirich

Sachse

K. A.

Henschel

Schnitzler

(2024). Comparing test-taking effort between paper-based and computer-based tests. Applied Psychological Measurement, 48(1-2), 3–17. https://doi.org/10.1177/01466216241227535

69.

Welch

B. L.

(1938). The significance of the difference between two means when the population variances are unequal. Biometrika, 29(3/4), 350–362. https://doi.org/10.1093/biomet/29.3-4.350

70.

Wolter

K. M.

(2007). Introduction to variance estimation (2nd ed.) Springer.

71.

Kim

S. Y.

Westine

Boyer

(2025). IRT observed-score equating for rater-mediated assessments using a hierarchical rater model. Journal of Educational Measurement, 62(1), 145–171. https://doi.org/10.1111/jedm.12425

72.

Wyse

A. E.

(2023). Two statistics for measuring the score comparability of computerized adaptive tests. Applied Psychological Measurement, 47(7-8), 513–525. https://doi.org/10.1177/01466216231209749

73.

Von Davier

(2010). Linking errors in trend estimation. In: Large-scale surveys: A case study (research report No. RR-10-10). Educational Testing Service.

74.

Yin

Bezirhan

Fishbein

Foy

(2023). Implementing the PIRLS 2021 achievement scaling methodology. In von Davier

Mullis

I. V. S.

Fishbein

Foy

(Eds.) Methods and procedures: PIRLS 2021 technical report (pp. 11.1–11.82): TIMSS & PIRLS International Study Center. https://doi.org/10.6017/lse.tpisc.tr2103.kb3067

75.

Yin

Fishbein

Bezirhan

Foy

von Davier

(2023). Examining country-level differences between digitalPIRLS data and bridge data. In von Davier

Mullis

I. V. S.

Fishbein

Foy

(Eds.) Methods and procedures: PIRLS 2021 technical report (pp. 12.1–12.33): TIMSS & PIRLS International Study Center. https://doi.org/10.6017/lse.tpisc.tr2103.kb9281

Supplementary Material

Please find the following supplemental material available below.

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

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

14.92 MB

0.00 MB