3.1 Model specification

At the first level, the conditional distributions of the observed item responses and response times given the latent traits are specified. The choice of the measurement model is flexible. Here, we assume that the probability that examinee p answers item i correctly, given latent ability $\theta _p$ follows the 2PL model:

(7) $$ \begin{align} P(Y_{pi} = 1 | \theta_p, a_i, b_i) = \frac{1}{1 + \text{exp}\big(-a_i (\theta_p - b_i)\big)}. \end{align} $$

And the observed response time $T_{pi}$ of person p on item i follows a log-normal distribution with an item time intensity parameter $\beta _i$ and time discrimination $\alpha _i$ , conditional on latent speed $\tau _p$ :

(8) $$ \begin{align} T_{pi} \sim \text{log} \mathcal{N}(\beta_i - \tau_p, \alpha_i^{-2}), \end{align} $$

with density

(9) $$ \begin{align} f(t_{pi}) = \frac{\alpha_i}{t_{pi} \sqrt{2 \pi}} \text{exp} \Big\{ -\frac{1}{2} \big[ \alpha_i (\text{log} t_{pi} - (\beta_i - \tau_p)) \big]^2 \Big\}. \end{align} $$

Item responses and response times are assumed to be independent conditional on the latent traits $(\theta _p, \tau _p)$ .

At the second level, we assume that the marginal distributions of the latent ability and speed factors follow normal distributions:

(10) $$ \begin{align} \begin{aligned} \theta_p &\sim \mathcal{N} (\mu_\theta, \sigma^2_\theta), \\ \tau_p &\sim \mathcal{N} (\mu_\tau, \sigma^2_\tau). \end{aligned} \end{align} $$

For identifiability of the model, we set $\mu _\theta = \mu _\tau = 0$ and $\sigma ^2_\theta = 1$ . The joint distribution of $(\theta _p, \tau _p)$ is defined via a copula function C with a dependence parameter $\delta $ to model the asymmetric dependence between the latent ability and speed. Let $F_\theta $ and $F_\tau $ denote the CDFs of $\theta $ and $\tau $ . The joint CDF of $\theta $ and $\tau $ is

(11) $$ \begin{align} F_{(\theta, \tau)} = C(F_\theta, F_\tau; \delta). \end{align} $$

For instance, if we employ a Gumbel copula to model upper-tail dependence, its copula function in the bivariate case (Table 1) is

(12) $$ \begin{align} C(u, v; \delta) = \text{exp} \{-([-\text{log} u]^\delta + [-\text{log} v]^\delta)^{1 / \delta}\}, \quad u, v \in [0, 1], 1 \leq \delta < \infty. \end{align} $$

Setting, $x = -\text {log} u, y = -\text {log} v$ , the copula density is

(13) $$ \begin{align} c(u, v; \delta) = \text{exp} \{-[x^\delta + y^\delta]^{1 / \delta} \} [(x^\delta + y^\delta)^{1 / \delta} + \delta - 1] [x^\delta + y^\delta]^{1 / \delta - 2} (xy)^{\delta - 1} (uv)^{-1}. \end{align} $$

From Equations (4) and (5), the joint PDF of $(\theta , \tau )$ is

(14) $$ \begin{align} f_{\theta, \tau} (\theta, \tau \mid \sigma^2_{\tau}, \delta) = c(F_{\theta}(\theta), F_{\tau}(\tau) \mid \delta) f_{\theta}(\theta) f_{\tau}(\tau \mid \sigma^2_\tau), \end{align} $$

where $f_{\theta }$ and $f_{\tau }$ denote the respective normal density functions.

Note that the traditional hierarchical model with multivariate normality assumption can be seen as a special case of the proposed framework when the marginal distributions of the latent traits are normal, and their joint distribution is defined via a Gaussian copula. In fact, this framework is general and supports a wide range of copula and marginal specifications. The Gumbel copula could be replaced with the Joe copula to model even stronger upper-tail dependence as detailed in Appendix B and implemented in this study. Alternatively, the Student’s t copula could reflect symmetric tail dependence at both extremes. Similarly, the latent trait marginals are not restricted to a normal distribution. One might assume, for instance, a gamma distribution for speed to account for skewness (Maris, Reference Maris1993), or any other form that best fits the data.

An EM algorithm (Dempster et al., Reference Dempster, Laird and Rubin1977) with stochastic approximation (Robbins & Monro, Reference Robbins and Monro1951) for maximum marginal likelihood estimation is used. Details of the estimation algorithm and inference are described in Appendix A.

3.2 Rotated copulas

The standard Gumbel and Joe copulas only model positive association with upper-tail dependence. However, in practice, the association between latent traits could reflect a negative association (e.g., speed–accuracy trade-off) or lower-tail dependence, or both. In order to model different types of dependence structure, we use rotated copulas (Nelsen, Reference Nelsen2006) as illustrated in Figure 2. For instance, rotating a Gumbel copula by $90^\circ $ creates a negative association with lower-tail dependence, while a $270^\circ $ rotation also creates a negative association but with upper-tail dependence.

Figure 2

Bivariate distributions of latent traits under rotated copulas.

Note: The marginal distributions of $\theta $ and $\tau $ are standard normal, and the joint distributions are defined using rotated copulas.

Figure 2 Long description

A multi-panel scatter plot arrangement with two rows and three columns. All plots share identical axes: the x-axis represents theta and the y-axis represents tau, both ranging from negative 4 to 4.

Top Row: Gumbel Copulas

* Top-Left: 90 degree rotated Gumbel. Data points show a strong negative correlation, forming a dense cluster that slopes downward from the top-left quadrant toward the bottom-right.

* Top-Middle: 180 degree rotated Gumbel. Data points show a strong positive correlation, sloping upward from the bottom-left toward the top-right.

* Top-Right: 270 degree rotated Gumbel. Data points show a negative correlation similar to the 90 degree version but with a slightly different tail density.

Bottom Row: Joe Copulas

* Bottom-Left: 90 degree rotated Joe. Data points show a negative correlation with a distinct asymmetric tail, where the density is higher and more concentrated in the upper-left area.

* Bottom-Middle: 180 degree rotated Joe. Data points show a positive correlation with a sharp, narrow tail extending from the bottom-left, widening as it moves toward the top-right.

* Bottom-Right: 270 degree rotated Joe. Data points show a negative correlation with a sharp tail extending from the bottom-right, widening toward the top-left.

Rotated copulas can be implemented by transforming the arguments u and v within the original copula functions, typically by replacing them with $1 - u$ or $1 - v$ . For instance, a $270^\circ $ -rotated copula, $C_{270}(u, v)$ , is derived from its original copula function $C(u, v)$ as follows (Nelsen, Reference Nelsen2006):

(15) $$ \begin{align} C_{270}(u, v) = u - C(u, 1 - v). \end{align} $$

The corresponding copula density, $c_{270}(u, v),$ is then given by

(16) $$ \begin{align} c_{270} (u, v) = \frac{\partial^2}{\partial u \partial v} C_{270} (u, v) = \frac{\partial^2}{\partial u \partial v} [u - C(u, 1 - v)] = c(u, 1 - v). \end{align} $$

Similarly, the $90^\circ $ - and $180^\circ $ -rotated copula functions are

(17) $$ \begin{align} \begin{aligned} C_{90}(u, v) &= v - C(1 - u, v), \\ C_{180}(u, v) &= u + v - 1 + C(1 - u, 1 - v), \end{aligned} \end{align} $$

with their respective densities being

(18) $$ \begin{align} \begin{aligned} c_{90} (u, v) &= c(1 - u, v), \\ c_{180} (u, v) &= c(1 - u, 1 - v). \end{aligned} \end{align} $$

The gradients of the joint log-likelihood for the EM algorithm are derived for each rotated copula, as detailed in Appendix C. For model selection, likelihood-based criteria such as Akaike’s information criterion (AIC; Akaike, Reference Akaike1998) or the Bayesian information criterion (BIC; Schwarz, Reference Schwarz1978) can be used. However, because each bivariate copula in this framework has the same number of parameters, comparing their marginal likelihoods is sufficient to find the best-fitting model.

3.3 Functional complexity of bivariate copulas

Although the number of unknown parameters is equal between the MVN and the copula models, the complexity of the models can be explained beyond counting the number of parameters. This section compares the functional complexities of the bivariate copulas considered in this study. Specifically, we compare the bivariate Gaussian, Gumbel, and Joe copulas using the functional complexity penalty derived from the minimum description length (MDL; Hansen & Yu, Reference Hansen and Yu2001; Rissanen, Reference Rissanen2002). MDL is defined as

(19) $$ \begin{align} \text{MDL} = -\ln (f(y \mid \hat{\Delta})) + (k/2) \ln(n / 2\pi) + \ln \int_\Delta \sqrt{\det (I(\Delta))} d\Delta, \end{align} $$

where $\Delta $ is the model parameter vector, $\hat \Delta $ represents the parameter values maximizing the likelihood, and $I(\Delta )$ is the Fisher information matrix. MDL accounts for three dimensions of model complexity: the number of parameters, sample size, and functional form (Pitt & Myung, Reference Pitt and Myung2002). The last term in Equation (19) quantifies this functional complexity. To compare the MVN and copula models, we calculate this term for the bivariate Gaussian, Gumbel, and Joe copulas with uniform margins. Because evaluating the full joint model for speed and accuracy is computationally intensive, we focus on the bivariate copulas, as they represent the only structural differences between the full joint models.

For a bivariate copula with copula density c and dependence parameter $\delta $ , the Fisher information is

(20) $$ \begin{align} I(\delta) = E\left[ \left(\frac{\partial}{\partial \delta} \ln c(u, v; \delta)\right)^2\right] = \int_0^1 \int_0^1 \left(\frac{\partial}{\partial \delta} \ln c(u, v; \delta)\right)^2 c(u, v; \delta) du dv. \end{align} $$

For each bivariate copula, we calculated the complexity term:

(21) $$ \begin{align} \ln \int_\delta \sqrt{\det (I(\delta))} d\delta. \end{align} $$

The integration domain depends on the parameter space of each copula: $\rho \in (-1, 1)$ for the Gaussian, and $\delta \in [1, \infty )$ for the Gumbel and Joe copulas. Analytical expressions for the Fisher information of bivariate Gaussian and Gumbel copulas are available (Schepsmeier & Stöber, Reference Schepsmeier and Stöber2014) but not for the Joe copula. In this analysis, the integrals in Equations (20) and (21) were approximated using numerical integration to maintain consistency across copulas. As a result, the functional complexity for the bivariate Gaussian copula was 1.99, while the complexities for the Gumbel and Joe copulas were 1.02 and 1.11, respectively. These results indicate that the functional form of the bivariate Gaussian copula is approximately twice as complex. The result is expected, as the bivariate Gaussian copula allows for both negative and positive dependence, whereas the Gumbel and Joe copulas only allow for positive dependence, doubling the fitting propensity of the Gaussian copula. Although the proposed copula-based approach involves a model selection step among rotated copulas, the final model used for parameter estimation employs a single bivariate copula. Admittedly, the process of evaluating multiple rotated models to determine the appropriate dependence direction increases the overall modeling complexity. However, because the evaluated copulas have a lower functional complexity than the Gaussian copula, any observed improvements in model fit of a final selected copula model over the MVN model cannot be merely a byproduct of functional flexibility. Rather, such improvements would indicate that the selected copula captures the true underlying dependence structure.

4.1 Simulation design

In $R = 500$ replications, the true latent ability $\theta _p$ and speed $\tau _p$ for N persons were randomly generated from a Gumbel or Joe copula (both with normal margins) or an MVN distribution. The marginal distribution of $\theta _p$ was a standard normal, and the marginal distribution of $\tau _p$ was a normal with mean $0$ and variance $\sigma ^2_\tau = 0.25$ as in Guo et al. (Reference Guo, Xu, Ying and Zhang2022). Two levels of latent trait dependence were considered, corresponding to Kendall’s $\tau $ = $0.3$ and $0.5$ , representing moderate and strong dependence, respectively. The transformations to the copula dependence parameters (Joe, Reference Joe1997, Reference Joe2014) are as follows. For the Gumbel copula,

(22) $$ \begin{align} \delta = -1 / (\text{Kendall's}~\tau - 1). \end{align} $$

That is, Kendall’s $\tau $ values of $0.3$ and $0.5$ correspond to dependence parameter values $\delta = 1.43$ and $\delta = 2$ , respectively. For the Joe copula, there is no closed-form expression for $\delta (\tau )$ . Instead, one finds $\delta $ numerically by solving

(23) $$ \begin{align} \text{Kendall's}~\tau = 1 + 2(2 - \delta)^{-1}[\text{digamma}(2) - \text{digamma}(2 / \delta + 1)], \end{align} $$

where $\text {digamma}(z) = \frac {d}{dz} \text {log} \Gamma (z)$ . This is implemented in the copula R package (Hofert et al., Reference Hofert, Kojadinovic, Maechler and Yan2025). Kendall’s $\tau $ values of $0.3$ and $0.5$ correspond to $\delta = 1.77$ and $\delta = 2.86$ , respectively, under the Joe copula. For the Pearson correlation $\rho $ in a bivariate normal distribution, we have

(24) $$ \begin{align} \rho = \text{sin}\left(\frac{\pi}{2} \text{Kendall's}~\tau\right). \end{align} $$

Figure 3 displays the bivariate distributions of the latent traits randomly generated under each dependence condition.

Figure 3

Bivariate distributions of latent ability and speed under different dependence conditions.

Figure 3 Long description

The figure consists of six scatter plots arranged in a 2 by 3 grid.

Row 1: Kendall's tau equals 0.3.

* Top-left panel (Gumbel): Shows a moderately dispersed cloud of points with a slight upward linear trend and some clustering in the upper-right quadrant.

* Top-center panel (Joe): Displays a similar upward trend but with a more pronounced concentration of points in the upper-right corner, indicating upper tail dependence.

* Top-right panel (M V N): Shows a symmetric, elliptical cloud of points centered at 0, 0 with a moderate positive correlation.

Row 2: Kendall's tau equals 0.5.

* Bottom-left panel (Gumbel): The upward linear trend is tighter than the tau equals 0.3 version, with a distinct asymmetric stretch toward the top-right.

* Bottom-center panel (Joe): Shows a very sharp, narrow concentration of points in the upper-right quadrant, illustrating strong upper tail dependence.

* Bottom-right panel (M V N): Displays a tighter, symmetric elliptical cloud compared to the top-right panel, indicating a stronger positive correlation.

All plots share the same axes: the horizontal x-axis represents theta ranging from negative 3 to 3, and the vertical y-axis represents tau ranging from negative 1.5 to 1.5.

Item responses and response times on I items, conditional on the latent traits, were randomly generated from a 2PL model in Equation (7) and a log-normal model in Equation (8), respectively. The number of persons was $N = 500$ or $1,000$ , and the number of items was $I = 10$ or $20$ . The true item parameter values were set according to the parameter estimates from the empirical data on the 2018 Programme for International Student Assessment (PISA) science assessment reported in Guo et al. (Reference Guo, Xu, Ying and Zhang2022) and were held fixed across replications (Table 2). In the $I = 10$ condition, the first 10 items were used. In summary, the simulation design resulted in 24 conditions, including three dependence types (Gumbel, Joe, and MVN), two dependence strengths (Kendall’s $\tau = 0.3, 0.5)$ , two sample sizes ( $N = 500, 1,000$ ), and two test lengths ( $I = 10, 20$ ).

Table 2

True item parameter values in the simulation

Table 2 Long description

The table consists of five columns: Item, a, b, alpha, and beta. There are 20 numbered items.

* Item 1: a 0.76, b 0.43, alpha 2.87, beta 0.29.

* Item 2: a 0.18, b 0.05, alpha 2.41, beta 0.00.

* Item 3: a 1.10, b minus 0.98, alpha 2.37, beta 0.81.

* Item 4: a 1.00, b minus 2.86, alpha 2.59, beta 0.64.

* Item 5: a 0.82, b 0.08, alpha 2.40, beta 0.77.

* Item 6: a 1.40, b minus 0.85, alpha 3.11, beta 0.72.

* Item 7: a 0.70, b minus 0.47, alpha 2.61, beta 0.52.

* Item 8: a 1.47, b 0.17, alpha 2.23, beta 0.12.

* Item 9: a 1.05, b minus 0.69, alpha 2.14, beta 0.32.

* Item 10: a 1.56, b minus 0.62, alpha 2.21, beta 0.52.

* Item 11: a 1.66, b minus 3.27, alpha 2.69, beta 0.10.

* Item 12: a 1.11, b minus 0.58, alpha 2.81, beta 0.27.

* Item 13: a 1.15, b minus 1.72, alpha 2.76, beta 0.20.

* Item 14: a 1.24, b 0.71, alpha 2.22, beta 0.77.

* Item 15: a 1.18, b 1.89, alpha 2.25, beta 0.40.

* Item 16: a 0.92, b minus 2.51, alpha 2.93, beta minus 0.32.

* Item 17: a 0.70, b minus 0.04, alpha 3.12, beta 0.27.

* Item 18: a 1.03, b minus 0.22, alpha 2.43, beta 0.10.

* Item 19: a 0.96, b minus 0.78, alpha 1.96, beta minus 0.60.

* Item 20: a 1.32, b minus 1.10, alpha 2.54, beta 0.37.

Three models were fit to the data in each replication: the conventional hierarchical model with a multivariate normality assumption (MVN), the Gumbel copula model, and the Joe copula model. The EM algorithm with stochastic approximation, as described in the previous section, was used to obtain the MML parameter estimates for the copula-based models. For the conventional model, parameter estimates were obtained using an EM algorithm described in Glas and van der Linden (Reference Glas and van der Linden2010). Model convergence was monitored using the change in average individual marginal log-likelihood ( $< 10^{-4}$ ). The standard errors of the parameter estimates were obtained from the approximate information matrix. The 95% confidence intervals were then constructed as

(25) $$ \begin{align} \hat{\boldsymbol{\eta}} \pm 1.96 \cdot \text{SE}(\hat{\boldsymbol{\eta}}), \end{align} $$

where $\hat {\boldsymbol {\eta }}$ denotes the vector of model parameter estimates and $\text {SE}(\hat {\boldsymbol {\eta }})$ represents their corresponding standard errors. Model selection was evaluated by comparing the marginal log-likelihoods of the fitted models. For each condition, the proportion of replications in which the proposed model yielded the highest marginal log-likelihood over the other models was computed. As a baseline, EAP estimates of latent ability were also obtained from the independence 2PL model (Indep) without response time data. All simulations were conducted in R, using the copula and mirt (Chalmers, Reference Chalmers2012) R packages, and custom functions.

The bias, mean absolute deviance (MAD), and RMSE for the item parameters, structural parameters, and ability estimates were calculated as follows:

(26) $$ \begin{align} \begin{aligned} \text{Bias}(\eta) &= \frac{1}{R} \sum_{r = 1}^R (\hat{\eta}^{(r)} - \eta), \\ \text{MAD}(\eta) &= \frac{1}{R} \sum_{r = 1}^R (|\hat{\eta}^{(r)} - \eta|), \\ \text{RMSE}(\eta) &= \sqrt{\frac{1}{R} \sum_{r = 1}^R (\hat{\eta}^{(r)} - \eta)^2}, \end{aligned} \end{align} $$

where $\hat {\eta }^{(r)}$ is the estimate at the r-th iteration, and $\eta $ is the true parameter value. To evaluate the accuracy of standard error estimates, the ratio between the average standard error estimates and the empirical standard deviation of the estimates was computed. For a given parameter $\eta $ , this ratio is defined as

(27) $$ \begin{align} \text{SE Ratio}(\hat{\eta}) = \frac{\frac{1}{R} \sum_{r=1}^R \text{SE}(\hat{\eta}^{(r)})}{\sqrt{\frac{1}{R} \sum_{r=1}^R \left(\hat{\eta}^{(r)} - \bar{\hat{\eta}}\right)^2}}, \end{align} $$

where $\text {SE}(\hat {\eta }^{(r)})$ denotes the standard error estimate in the r-th replication, and $\bar {\hat {\eta }} = \frac {1}{R} \sum _{r=1}^R \hat {\eta }^{(r)}$ is the average estimate across replications. The 95% coverage rate was computed as the proportion of replications in which the true parameter value fell within the 95% confidence interval:

(28) $$ \begin{align} \text{95% C.R.}(\eta) = \frac{1}{R} \sum_{r=1}^{R} \mathbb{I} \left\{ \eta \in \left[ \hat{\eta}^{(r)} - 1.96 \cdot \text{SE}(\hat{\eta}^{(r)}), ~ \hat{\eta}^{(r)} + 1.96 \cdot \text{SE}(\hat{\eta}^{(r)}) \right] \right\}, \end{align} $$

where $\mathbb {I}$ is the indicator function.

In general, when the true latent traits are generated under strong upper-tail dependence (i.e., Gumbel or Joe copula), the copula models are expected to outperform the conventional model. Conversely, when the data are generated under the MVN distribution, the conventional model is expected to yield better performance than the copula models. The marginal likelihood-based model selection is expected to perform well in identifying the true data-generating models.

4.2 Simulation results

4.2.1 Model selection

Table 3 shows the proportion of replications in which each of the three fitted models yielded the highest marginal log-likelihood under each true data-generating condition. The average marginal log-likelihood values are presented in Appendix D (Table D1). The likelihood-based model selection proved effective at correctly identifying the true model. For data generated under the Gumbel and Joe copulas, selection accuracy ranged from 0.69 to 0.98 and 0.84 to 0.98, respectively. When data were generated under the MVN structure, the MVN model was chosen in 0.87–0.98 of the replications. Across all conditions, model selection accuracy improved with larger sample sizes and longer tests.

Table 3

Marginal likelihood-based model selection result

Table 3 Long description

The table consists of 12 columns and 8 data rows.

Columns 1 to 3 define the parameters:

- N (sample size): 500 or 1,000.

- I (test length): 10 or 20.

- tau sub K (Kendall's rank coefficient): 0.3 or 0.5.

Columns 4 to 12 are grouped under three ‘True’ model conditions, each containing three sub-columns for fitted models (Gumbel, M V N, Joe):

1. True: Gumbel condition:

- For N=500, I=10, tau=0.3: Gumbel 0.69, M V N 0.16, Joe 0.15.

- For N=1,000, I=10, tau=0.3: Gumbel 0.89, M V N 0.06, Joe 0.05.

- For N=500, I=20, tau=0.3: Gumbel 0.77, M V N 0.10, Joe 0.13.

- For N=1,000, I=20, tau=0.3: Gumbel 0.94, M V N 0.02, Joe 0.04.

- For N=500, I=10, tau=0.5: Gumbel 0.77, M V N 0.20, Joe 0.03.

- For N=1,000, I=10, tau=0.5: Gumbel 0.92, M V N 0.07, Joe 0.00.

- For N=500, I=20, tau=0.5: Gumbel 0.90, M V N 0.08, Joe 0.02.

- For N=1,000, I=20, tau=0.5: Gumbel 0.98, M V N 0.01, Joe 0.01.

2. True: M V N condition:

- Across all rows, the M V N fitted model is correctly selected with high proportions ranging from 0.87 to 0.98. Joe is consistently 0.00, and Gumbel ranges from 0.02 to 0.13.

3. True: Joe condition:

- Across all rows, the Joe fitted model is correctly selected with proportions ranging from 0.81 to 0.98. M V N is consistently 0.00, and Gumbel ranges from 0.02 to 0.19.

General trend: Model selection accuracy increases as sample size N and test length I increase.

Note: N is the sample size, and I is the test length. $\tau _K$ denotes Kendall’s $\tau $ rank coefficient. The main column headings indicate the true dependence structure used to generate the data. The subheadings indicate the fitted models. Values represent the proportion of replications in which a given fitted model was selected as the best-fitting model based on the marginal likelihood.

4.2.2 Model parameter recovery

Table 4 presents the parameter recovery results for the Gumbel copula model, estimated using the EM algorithm with stochastic approximation, when the data were generated from the Gumbel copula dependence structure. Note that the item parameter recovery results presented here are averaged across all items. The full, item-level parameter recovery tables can be found in Appendix D. The proposed estimation algorithm accurately recovered the true parameter values. Across all conditions, bias for item parameters ranged from $-$ 0.036 to 0.001 (MAD: 0.016–0.226), while the structural parameter bias ranged from $-$ 0.136 to 0.000 (MAD: 0.010–0.146). The RMSE of parameter estimates consistently decreased as sample size N increased. Note that the relatively large RMSE for the item difficulty parameter $\mathbf {b}$ was largely driven by items with low discrimination (Item 2, $a_i = 0.18$ ) or extreme difficulty (Item 11, $b_i = -3.27$ ), as shown in Table D5. The 95% coverage rates ranged from 0.72 to 0.98, and the SE ratios ranged from 0.86 to 1.41, indicating that the estimated standard errors generally approximated the empirical standard errors well. The item-level parameter recovery results were highly similar between the Gumbel model and the MVN model, showing the choice of dependence model in the trait level has minimal impact on the item parameter estimates (Appendix D, Tables D5 and D6). The parameter recovery for the Joe copula model was generally as accurate as the Gumbel model, except for the dependence parameter, $\delta $ , which was underestimated when dependence was strong (Appendix D, Table D2).

Table 4

Parameter recovery under the Gumbel copula dependence condition

Table 4 Long description

The table is organized into four horizontal blocks, each sharing the same ‘Condition’ columns on the left: N (sample size 500 or 1,000), I (test length 10 or 20), and delta (1.43 or 2).

* The first block displays Bias and M A D (Mean Absolute Deviation) for parameters a, b, alpha, beta, delta, and sigma sub tau. Bias values are generally low, ranging from minus 0.136 to 0.022. M A D values range from 0.009 to 0.226.

* The second block displays R M S E (Root Mean Square Error) and 95% coverage rate. R M S E values range from 0.012 to 1.633. The 95% coverage rates are mostly between 0.88 and 0.98, with a notable dip to 0.72 for delta when N is 1,000, I is 10, and delta is 2.

* The third block displays the S E (Standard Error) ratio for the same parameters. Most ratios are close to 1.0, ranging from 0.86 to 1.56.

Note: a, b, alpha, and beta represent item parameters, while delta and sigma sub tau represent structural parameters.

Note: N represents the sample size. I represents the test length. $\boldsymbol {a}, \boldsymbol {b}, \boldsymbol {\alpha }$ , and $\boldsymbol {\beta }$ represent item parameters, and $\delta $ and $\sigma _\tau $ represent structural parameters.

4.2.3 Latent ability estimation accuracy

Table 5 presents the accuracy of the latent ability estimates under Gumbel copula dependence conditions. The correctly specified Gumbel copula model yielded the lowest RMSE and the highest Pearson correlation with the true ability across all conditions. While its bias was slightly higher than misspecified models in some conditions (e.g., $\delta = 2$ ), its lower RMSE and higher correlation indicate more accurate ability estimation. This general finding was consistent across simulations where the model that matched the true dependence structure consistently outperformed the others in overall accuracy (Table D3).

Table 5

Accuracy of ability estimates under the Gumbel copula dependence condition

Table 5 Long description

The table is structured with 15 columns. The primary headers are Condition (subdivided into N, I, and delta), Bias of theta (subdivided by Gumbel, Joe, M V N, and Indep models), R M S E of theta (subdivided by the same four models), and Correlation of theta-hat and theta (subdivided by the same four models).

Key data points include:

* For N = 500, I = 10, delta = 1.43: Bias ranges from 0.005 to 0.007; R M S E ranges from 0.564 to 0.598; Correlation ranges from 0.800 to 0.825.

* For N = 1,000, I = 10, delta = 1.43: Bias ranges from 0.000 to 0.003; R M S E ranges from 0.561 to 0.596; Correlation ranges from 0.803 to 0.828.

* For N = 500, I = 20, delta = 1.43: Bias ranges from 0.004 to 0.005; R M S E ranges from 0.459 to 0.480; Correlation ranges from 0.877 to 0.888.

* For N = 1,000, I = 20, delta = 1.43: Bias is 0.000 or 0.001; R M S E ranges from 0.457 to 0.477; Correlation ranges from 0.879 to 0.890.

* For N = 500, I = 10, delta = 2: Bias ranges from 0.005 to 0.020; R M S E ranges from 0.509 to 0.597; Correlation ranges from 0.801 to 0.860.

* For N = 1,000, I = 10, delta = 2: Bias ranges from 0.001 to 0.017; R M S E ranges from 0.507 to 0.596; Correlation ranges from 0.803 to 0.862.

* For N = 500, I = 20, delta = 2: Bias ranges from 0.005 to 0.011; R M S E ranges from 0.424 to 0.480; Correlation ranges from 0.877 to 0.906.

* For N = 1,000, I = 20, delta = 2: Bias ranges from 0.001 to 0.007; R M S E ranges from 0.421 to 0.477; Correlation ranges from 0.879 to 0.907.

Across all conditions, the Gumbel model generally shows the lowest R M S E and highest correlation compared to the Indep baseline.

Note: N represents the sample size. I represents the test length. The column headers (Gumbel, Joe, MVN, and Indep) indicate the estimating model, while Gumbel was the data-generating model. Indep represents the baseline independence 2PL model without response time.

In addition to evaluating the overall accuracy of ability estimates, the RMSE and bias of $\theta $ conditional on the true ability level were computed to assess how estimation accuracy varies across the ability range. The conditional RMSE and bias values were calculated using equally spaced bins of width $0.5$ , covering the range of $\theta $ from $-4$ to $4$ . Under the Gumbel copula dependence condition (Figures 4 and 5), the copula models’ advantage over both the conventional model and the independence model is most prominent in the upper tail of the $\theta $ distribution. This effect is even stronger when the dependence between latent traits increases (Kendall’s $\tau = 0.5$ ). When the true data structure exhibits strong upper-tail dependence, response-time data provide more information about high-ability individuals and less information about low-ability individuals. The copula-based models borrow more information from response times for upper-tail abilities while borrowing less for lower-tail abilities, which aligns better with the true data structure in the simulation. In contrast, the conventional model’s MVN assumption implies asymptotic tail independence and therefore borrows less information from the response time data in the upper tail of ability. This discrepancy explains the performance gap in the upper tail between these models. The conditional RMSE shows a similar pattern when data are generated under the Joe copula dependence structure (Figure D1). Latent speed estimation accuracy was highly similar across the different dependence models, with negligible differences in bias, RMSE, and correlation (Table D4 and Figure D2).

Figure 4

Conditional RMSE of ability estimates under the Gumbel copula dependence condition.

Figure 4 Long description

A multi-panel figure with four line graphs arranged in a two by two grid. All graphs share the same axes: the x-axis represents theta ranging from minus 4 to 4, and the y-axis represents R M S E ranging from 0.0 to 2.0. Each graph contains four lines: a solid black line for gumbel, a dashed red line for m v n, a dotted green line for indep, and a dash-dot blue line for joe.

Top-left panel (N = 500, I = 10, delta = 1.4): All lines follow a U-shaped curve, overlapping closely between theta values of minus 2 and 2 with a minimum R M S E near 0.5 at theta equals 0. For theta greater than 2, the indep line rises sharply toward 2.0, followed by m v n, while gumbel and joe remain lower and closely aligned.

Top-right panel (N = 1000, I = 10, delta = 1.4): The trends are nearly identical to the top-left panel, showing that doubling the sample size N has minimal impact on the R M S E distribution for these parameters.

Bottom-left panel (N = 500, I = 10, delta = 2): The U-shape is maintained, but the divergence at higher theta values is more pronounced. The indep line reaches 2.0 much faster (around theta equals 3). The gumbel and joe lines show a slight peak around theta equals 3 before dipping slightly at theta equals 4.

Bottom-right panel (N = 1000, I = 10, delta = 2): Similar to the bottom-left, the indep line shows the highest error at positive theta. The gumbel and joe models continue to show the lowest R M S E across the entire theta range compared to the independence and multivariate normal models.

Figure 5

Conditional bias of ability estimates under the Gumbel copula dependence condition.

Figure 5 Long description

A four-panel line graph display. Each panel shares the same axes: the x-axis represents theta ranging from minus 4 to 4, and the y-axis represents Bias ranging from minus 2 to 2. A horizontal dashed gray line marks zero bias.

Legend in the top-left panel identifies four models:

- gumbel: solid black line

- m v n: dashed pink line

- indep: dotted green line

- joe: dash-dot blue line

Panel 1 (Top-Left): N = 500, I = 10, delta = 1.4. All lines show a downward trend. They overlap closely from theta = minus 4 to 0. For theta > 0, the indep line drops most sharply, followed by m v n. The gumbel and joe lines level off near Bias = minus 1.

Panel 2 (Top-Right): N = 1000, I = 10, delta = 1.4. The trends are nearly identical to Panel 1, indicating stability with increased sample size.

Panel 3 (Bottom-Left): N = 500, I = 10, delta = 2. The divergence between models for theta > 0 is more pronounced. The gumbel and joe lines level off earlier and higher (around Bias = minus 0.5) compared to the delta = 1.4 panels.

Panel 4 (Bottom-Right): N = 1000, I = 10, delta = 2. Similar to Panel 3, showing that higher delta values increase the separation between the independence model (indep) and the copula-based models (gumbel, joe) at high theta values.

When the data were generated under the MVN dependence structure, the conventional MVN model exhibited the smallest overall RMSE and the highest correlation with the true values (Table 6). Similarly, the conditional RMSE was the lowest for the MVN in both the lower tail and the mid-range of $\theta $ , although the differences were less substantial than in the Gumbel copula condition (Figure 6). These results confirm that the conventional MVN model provides the most accurate ability estimates when the underlying dependence structure is MVN.

Table 6

Accuracy of ability estimates under the MVN dependence condition

Table 6 Long description

A table with 15 columns and 8 data rows.

Headers:

- Condition: N, I, rho.

- Bias(theta): Gumbel, Joe, M V N, Indep.

- R M S E(theta): Gumbel, Joe, M V N, Indep.

- Cor(theta hat, theta): Gumbel, Joe, M V N, Indep.

Data Rows:

1. N 500, I 10, rho 0.45: Bias ranges from -0.004 to -0.001; R M S E ranges from 0.575 to 0.599; Cor ranges from 0.800 to 0.818.

2. N 1,000, I 10, rho 0.45: Bias ranges from -0.003 to 0.000; R M S E ranges from 0.572 to 0.595; Cor ranges from 0.803 to 0.820.

3. N 500, I 20, rho 0.45: Bias ranges from -0.004 to -0.001; R M S E ranges from 0.467 to 0.480; Cor ranges from 0.877 to 0.884.

4. N 1,000, I 20, rho 0.45: Bias ranges from -0.004 to 0.000; R M S E ranges from 0.464 to 0.477; Cor ranges from 0.879 to 0.886.

5. N 500, I 10, rho 0.71: Bias ranges from -0.001 to 0.013; R M S E ranges from 0.522 to 0.599; Cor ranges from 0.800 to 0.853.

6. N 1,000, I 10, rho 0.71: Bias ranges from 0.000 to 0.014; R M S E ranges from 0.519 to 0.595; Cor ranges from 0.803 to 0.854.

7. N 500, I 20, rho 0.71: Bias ranges from -0.002 to 0.002; R M S E ranges from 0.436 to 0.480; Cor ranges from 0.877 to 0.900.

8. N 1,000, I 20, rho 0.71: Bias ranges from -0.001 to 0.002; R M S E ranges from 0.432 to 0.477; Cor ranges from 0.879 to 0.901.

Note: N is sample size, I is test length, and Indep is the baseline independence 2 P L model.

Note: N represents the sample size. I represents the test length. The column headers (Gumbel, Joe, MVN, and Indep) indicate the estimating model, while Gumbel was the data-generating model. Indep represents the baseline independence 2PL model without response time.

Figure 6

Conditional RMSE of ability estimates under the MVN dependence condition.

Figure 6 Long description

A multi-panel figure containing four line graphs. All graphs share a common Y-axis labeled R M S E ranging from 0.0 to 2.0 and an X-axis labeled theta ranging from -4 to 4. A legend in the top-left panel identifies four models: gumbel (solid black line), m v n (dashed red line), indep (dotted green line), and joe (dash-dot blue line).

* Top-Left Panel: N = 500, I = 10, rho = 0.5. All four lines follow a U-shaped curve, overlapping closely with a minimum R M S E of approximately 0.5 near theta = 0. The curves rise more steeply as theta moves toward 4 than toward -4.

* Top-Right Panel: N = 1000, I = 10, rho = 0.5. The trends are nearly identical to the top-left panel, showing that increasing N from 500 to 1000 has minimal impact on the R M S E distribution for these models.

* Bottom-Left Panel: N = 500, I = 10, rho = 0.7. The U-shaped curves begin to diverge. The indep (green dotted) line shows a significantly higher R M S E as theta increases beyond 2, reaching above 2.0. The m v n (red dashed) line also rises higher than the gumbel and joe lines at the positive extreme of theta.

* Bottom-Right Panel: N = 1000, I = 10, rho = 0.7. Similar to the bottom-left, the indep model shows the highest error at high theta values. The gumbel and joe models remain the most stable with the lowest R M S E across the theta range.

We conducted additional simulations to evaluate model performance using data generated from $270^\circ $ -rotated Gumbel and Joe copulas, which reflect a negative association with upper-tail dependence. These simulations compared five models, including the two $270^\circ $ -rotated copulas, their unrotated counterparts, and the MVN model. Parameter recovery for the rotated models was comparable to the unrotated scenarios, showing minimal bias except for an underestimation of the $\delta $ parameter for the rotated Joe copula in the strong dependence condition. The marginal likelihood-based model selection effectively selected the true data-generating models, with selection accuracies ranging from 0.71 to 0.99 for the rotated Gumbel and 0.83 to 0.98 for the rotated Joe. Importantly, the misspecified unrotated models (with positive associations) were never selected. In terms of ability estimation, the rotated copula models again yielded the lowest conditional RMSE, particularly in the upper tail, outperforming the MVN model. The baseline independence model and misspecified unrotated copula models were the least accurate. (See Appendix E for full results.)

In summary, model parameters were recovered accurately, and the true data-generating model could be reliably identified using the marginal log-likelihood. Most importantly, the simulation results show that by correctly specifying the dependence structure, particularly tail dependencies, the copula models produced significantly more accurate ability estimates than the standard MVN and independence models. This advantage was most pronounced at the extremes of the ability scale, confirming that the model’s flexibility improved measurement accuracy.

5.1 Random-split replication analysis

In addition, we conducted a random-split replication analysis to evaluate the replicability of the model selection and to test whether asymmetric tail dependence is consistently observed across the PISA dataset. For both the math and science assessments, we aggregated country data within each cluster and randomly split the students into five folds. As a result, 30 datasets were fit to the same nine models, and the best model was selected based on the highest log-likelihood value. The analysis results suggest that the model selection is consistent within the PISA data, as the $90^\circ $ -rotated Joe copula model was selected for 28 of the 30 datasets. This provides strong evidence for the presence of asymmetric tail dependence in large-scale educational assessments. The log-likelihood values and selected models are summarized in Appendix F (Table F2).

5.2 Item response-level predictive accuracy

Using the InterModel Vigorish (IMV; Domingue, Kanopka, Kapoor, et al., Reference Domingue, Kanopka, Kapoor, Pohl, Chalmers, Rahal and Rhemtulla2024), the item response-level predictive accuracy of the $90^\circ $ -rotated Joe copula model and the MVN model was compared. IMV is an index that quantifies the predictive gain of a compared model $M1$ against a baseline model $M0$ for binary responses. Let $\hat p_{pi}$ denote the predicted item response for examinee p and item i from a model. For each prediction model, the uncertainty in binary prediction is quantified by finding a $w \in [0.5, 1]$ value such that

(29) $$ \begin{align} w \log (w) + (1 - w) \log (1- w) = \log \left( \prod_{p = 1}^N \prod_{i = 1}^I \hat p_{pi}^{y_{pi}} (1 - \hat p_{pi})^{1 - y_{pi}} \right)^{\frac{1}{N \times I}}. \end{align} $$

A w value closer to 0.5 indicates higher uncertainty in binary prediction. Then, the IMV is calculated as

(30) $$ \begin{align} \text{IMV} = \frac{w_1 - w_0}{w_0}. \end{align} $$

Here, $w_1$ and $w_0$ are the levels of uncertainty in $M1$ and $M0$ , respectively. IMV values greater than zero quantify the gain in replacing $M0$ with $M1$ for binary outcome prediction.

To evaluate the out-of-sample prediction accuracy of the models, the “missing response” paradigm was implemented for each of the 30 PISA datasets by randomly splitting the response-level data into five folds. Across five iterations, item responses and response times in each fold were treated as missing, and the models were fit using the remaining data. The predicted values for the missing responses were obtained using the item parameter estimates and the latent ability estimates from the 90J model, the MVN model, and the independence model. For each dataset, the IMV values were averaged across the five folds. Figure 11 shows the IMV values for the 30 PISA datasets, and the black diamonds indicate the average IMV across all 30 datasets. The two panels comparing each model against the independence model show that both hierarchical approaches outperformed the independence model. The average IMV of the 90J model against the MVN model was 0.00027. This suggests a weak improvement in predictive accuracy for the 90J relative to the MVN, a magnitude comparable to the predictive gain of a 3PL model over a 2PL model found in Domingue, Kanopka, Kapoor, et al. (Reference Domingue, Kanopka, Kapoor, Pohl, Chalmers, Rahal and Rhemtulla2024).

Figure 11

IMV values from the 2018 PISA Math and Science Data.

Figure 11 Long description

A horizontal dot plot titled I M V across 30 P I S A Datasets. The x-axis represents I M V values ranging from minus 0.002 to 0.004 with a dashed vertical line at 0.000. The y-axis lists three model comparisons.

* Top row: mvn:rot_90_joe in purple. The dots are clustered around the 0.000 line with a black diamond marker at 0.00027.

* Middle row: indep:rot_90_joe in red. The dots are shifted to the right of the 0.000 line with a black diamond marker at 0.00113.

* Bottom row: indep:mvn in blue. The dots are clustered slightly to the right of the 0.000 line with a black diamond marker at 0.00086.

A legend on the right identifies the colors: blue for indep:mvn, red for indep:rot_90_joe, and purple for mvn:rot_90_joe.

Figure 12 presents person-level predicted accuracy against the observed accuracy for the 90J, MVN, and independence models across the 30 PISA datasets. For each individual i, the predicted accuracy is the mean of the predicted response probabilities obtained from the five-fold out-of-sample prediction, while the observed accuracy is the proportion of correct responses across the items. Overall, the accuracy of low-performing students tends to be over-predicted, whereas the accuracy of high-performing students is under-predicted. However, the 90J and MVN models mitigate this bias in the lower tail by incorporating response time information. Specifically, the 90J model applies a stronger adjustment to fast-responding examinees within the low-performing group (observed accuracy < 0.3), drawing their predicted accuracy closer to the reference diagonal line. Similarly, Figure 13 plots the person-level predicted median response time against the observed median response time. For each individual, the predicted median response time was calculated as the median of the predicted median response time (i.e., $\text {Median}(\hat \beta _j - \hat \tau _i)$ ). The 90J model better predicts the median response time for slow-responding examinees (e.g., observed median RT > 150) compared to the MVN model. In terms of item-level accuracy and response time predictions, the two methods showed negligible differences (Appendix F, Figures F1 and F2).

Figure 12

Person-level predicted accuracy versus observed accuracy from the 2018 PISA Math and Science Data.

Figure 12 Long description

A multi-panel figure containing three scatter plots arranged horizontally. All three plots share identical axes. The x-axis is labeled Observed Accuracy and the y-axis is labeled Predicted Accuracy, both ranging from 0.0 to 1.0. A dashed diagonal identity line representing perfect prediction runs from the bottom-left to the top-right of each plot.

* The first plot on the left is titled 90 J and features black data points. The points are clustered in vertical bands, showing a strong positive linear correlation that closely follows the diagonal line.

* The middle plot is titled M V N and features red data points. Similar to the first plot, the data is organized in vertical clusters along the diagonal, though the density appears slightly more concentrated in the center.

* The third plot on the right is titled Indep and features light green data points. The distribution follows the same positive linear trend as the previous two, with data points clustered in vertical columns along the identity line.

In all three panels, the predicted accuracy values closely match the observed accuracy values, with the majority of data points falling between 0.2 and 0.9 on both axes.

Figure 13

Person-level predicted median RT versus observed median RT from the 2018 PISA Math and Science Data.

Figure 13 Long description

A two-panel scatter plot display. Both panels share identical axes: the x-axis represents Observed Median R T and the y-axis represents Predicted Median R T, with scales ranging from 0 to 300 in increments of 100. A dashed gray diagonal line representing a perfect one-to-one prediction (y equals x) runs from the origin to the top-right corner of each plot.

* The left panel, titled 90 J, features black semi-transparent data points. The points are densely clustered along the dashed identity line, showing a strong positive linear trend. The density is highest between 50 and 150 on both axes, with some scattering and increased variance as values exceed 200.

* The right panel, titled M V N, features red semi-transparent data points. Similar to the first panel, the data follows the dashed identity line closely, indicating high predictive accuracy. The spread of the red points appears slightly tighter around the identity line compared to the 90 J panel, particularly in the 100 to 200 range.