2.1 First level: The observation equations in DIR-RT models

Equations (1.1) and (1.2) are based on a one-time exam for each test taker. To consider a much complex data structure in the computerized test (as in the EdSphere testbed), we first expand the labels of notations. Let $X_{i,t,s,l}$ be the item response to indicate the correctness for the answer of the lth item in the sth test on the tth day given by the ith person, where $i=1,\ldots ,n$ (number of subjects); $t=1,\ldots ,T_i$ (number of test dates); $s=1,\ldots ,S_{i,t}$ (number of tests in a day); and $l=1,\ldots ,K_{i,t,s}$ (number of items in a test). Likewise, denote the difficulty of the lth item as $d_{i,t,s,l}$ . It is ideal to record the time for each examinee spending on a single item, however, in practice, more often the time spent on the entire exam is merely stored for each individual (a case for reading comprehension tests in EdSphere datasets). Thus, in our proposed models, the response time is defined at the test level, i.e., $R_{i,t,s}$ , implying the time spent on the sth test for the ith individual on the tth day; whereas, our models can be easily revised to cope with the response time stored for each item whenever such data are available.

The label extension illustrates two major features of computerized (adaptive) testing: 1) there are few replications of items among different time, tests, and test takers and 2) the observations are recorded at individually-varying and irregularly-spaced time points. In our study, only $X_{i,t,s,l}$ ’s and $R_{i,t,s}$ ’s are observed. Usually, the response time is naturally bounded away from zero, and a logarithmic transformation of $R_{i,t,s}$ is taken to remove its skewness (Fox & Marianti, Reference Fox and Marianti2016; Liu et al., Reference Liu, Wang, Hancock and Chen2022; Van der Linden et al., Reference Van der Linden, Klein Entink and Fox2010).

2.1.1 The observation equations of item responses

Often in a design of computerized tests, item difficulty, i.e., $d_{i,t,s,l}$ , is a randomized parameter, assumed to be randomly selected from a bank of items with a certain ensemble mean. Thus, we model $d_{i,t,s,l}=a_{i,t,s}+\epsilon _{i,t,s,l}$ as a measurement error model with $a_{i,t,s}$ being an ensemble mean difficulty of items in the sth test, and $\epsilon _{i,t,s,l}\sim \mathcal {N}(0,\sigma ^2)$ with $\sigma ^2$ known according to the test design. Here, $\mathcal {N} (\cdot , \cdot )$ denotes a normal distribution. Similarly, as Wang et al. (Reference Wang, Berger and Burdick2013) did, we extend classic one-parameter IRT models as

(2.1) $$ \begin{align} &\mathrm{Pr}(X_{i,t,s,l}=1\mid \theta_{i,t}, \varphi_{i,t}, \eta_{i,t,s}, a_{i,t,s}) =\mathrm{F}(\theta_{i,t}-a_{i,t,s}+\varphi_{i,t}+\eta_{i,t,s}+\epsilon_{i,t,s,l}), \end{align} $$

where $\theta _{i,t}$ represents the ith person’s ability on day t, with the assumption that one’s ability is constant over a given day; $\varphi _{i,t}$ and $\eta _{i,t,s}$ take account of daily and test random effects, respectively, to explain the possible local dependence of item responses. Further, assume $\varphi _{i,t}\sim \mathcal {N}(0,\delta _{i}^{-1})$ with its precision unknown and different for each person. To make $\eta _{i,t,s}$ well separate from $\varphi _{i,t}$ in inference, we denote ${\eta }_{i,t}=(\eta _{i,t,1},\ldots ,\eta _{i,t,S_{i,t}})'$ as the vector of test random effects on day t for the individual i and let $\mathbf {I}_{S_{i,t}}$ be an $S_{i,t}\times S_{i,t}$ identity matrix, then assume ${\eta }_{i,t} \sim \mathcal {N}_{S_{i,t}}(0,\tau _i^{-1}\mathbf { I}_{S_{i,t}} \mid \sum _{s=1}^{S_{i,t}}\eta _{i,t,s}=0)$ . The reason for letting $\eta _{i,t}$ be a singular multivariate normal (by setting the sum of test random effects to be zero on a day t) is to remove any possibility of unidentifiability issues between daily and test random effects. Notice that we use precision parameters instead of variance parameters for normal distributions as it is more convenient to draw precision parameters in the MCMC. Here, we choose $\mathrm {F}(x)$ to be a logistic link in the analysis as per the convention in the EdSphere platform.

2.1.2 The observation equations of response times

When working on a task, a subject has the choice to work faster or slower. However, when we work faster, we often incline to make more unnoticeable mistakes. On the other hand, when one’s ability is not comparable to the test taken, it always takes a longer time for him/her to finish the test. Thus, we model the response time as below:

(2.2) $$ \begin{align} &\text{log} (R_{i,t,s})=\mu_{i}- \nu_{i,t}+\beta L(\theta_{i,t}-a_{i,t,s})+\zeta_{i,t,s}. \end{align} $$

Equation (2.2) is an extension of Thissen’s model (Thissen, Reference Thissen and Weiss1983) and its variations (Ferrando & Lorenzo-Seva, Reference Ferrando and Lorenzo-Seva2007; Ranger & Kuhn, Reference Ranger and Kuhn2012). In Equation (2.2), $\mu _{i}$ reflects the average response time for the ith respondent in general. $\nu _{i,t}$ implies the variation for the speed of the respondent i on the tth day, with the negative sign indicating that the slower the speed is, the more time one needs to spend on the exam. We further assume that the speed for an examinee will not change much during one day, thus the index of the speed only varies according to individuals and days and $\nu _{i,t}\sim \mathcal {N}(0, \kappa _i^{-1})$ , which has an individual-specific precision parameter $\kappa _i$ and zero mean (to ensure identifiability in the presence of $\mu _i$ ). In the third term, $L(\theta _{i,t}-a_{i,t,s})$ is mapping how the distance of $\theta _{i,t}-a_{i,t,s}$ in relationship to the response time; and $\beta $ is a regression coefficient to adjust this relationship. Further, we assume the residual term $\zeta _{i,t,s},\sim \mathcal {N}(0,\varrho ^{-1})$ with $\varrho $ as a common precision parameter to borrow the strength of the information across different tests and individuals. Although $\varrho $ varying across individuals may be an alternative, such an assumption might cause unidentifiability issue with precision parameter $\kappa _i$ when we encounter the situation that an examinee only takes one test per day.

For $L(\cdot )$ , there are two popular choices in the literature. One is $L(\cdot )= \cdot $ , the monotone relationship and the other is $L(\cdot )=|\cdot |$ , the inverted U-shaped relationship. These two relationships represent two different psychological behaviors of test takers during an exam as discussed in Section 1.2. As far as we know, there is no formal statistical comparison of these two relationships in real examples of educational tests. We aim to fill in this gap by developing statistical technique for such comparison and applying to EdSphere datasets.

2.2 Second level: System equations in the DIR-RT models

Following the idea of Wang et al. (Reference Wang, Berger and Burdick2013), we combine both parametric growth models and Markov chain models to capture one’s ability growth over time, that is,

(2.3) $$ \begin{align} &\theta_{i,t}=\theta_{i,t-1}+c_i(1-\rho\theta_{i,t-1})\Delta_{i,t}^{+}+w_{i,t}. \end{align} $$

The first term in Equation (2.3) denotes the ability at the previous time point, $\theta _{i,t-1}$ . The second term represents a parametric growth model with $c_i$ as the average growth rate of the ith person’s ability over time; $\Delta _{i, t}^{+}=\min \{\Delta _{i,t}, T_{\max }\}$ is the time lapse between two test dates for the ith individual (i.e., $\Delta _{i, t}$ ) but truncated by a pre-specified maximum time interval $T_{\max }$ ( $T_{\max }=14$ is used in the application, reflecting typical holiday time for students in school); $\rho $ is the parameter to control the rate of one’s growth, which reduces the growth rate when the ability becomes mature. Note that $\rho $ is known from empirical experiments in EdSphere datasets (Swartz et al., Reference Swartz, Hanlon, Childress and Stenner2016; Wang et al., Reference Wang, Berger and Burdick2013). In principle, $\rho $ should be individual-specific, but it is distinguishable from $c_i$ only when one’s ability level is reaching maturation; our investigation of ability growth in the EdSphere data focuses on early age students, so only the $c_i$ are made individual-specific. Lastly, $w_{i,t}\sim \mathcal {N}(0,\phi ^{-1}\Delta _{i, t})$ represents the uncertainty that cannot be explained by the first two terms in Equation (2.3), where $\phi $ is a common unknown parameter to help borrow information and avoid substantial risks of confounding in the likelihood between $\delta _i$ ’s and $\phi ^{-1}\Delta _{i,t}$ when the time lapses between tests are equally spaced for a student. The variance of $w_{i,t}$ presumes that one’s ability changes become more uncertain, if he/she is absent for a long period. Moreover, we can write Equation (2.3) as a first-order Markov process (see Step 2 of Appendix A.1), which is beneficial for conducting MCMC later.

2.3 A summary of DIR-RT models

To summarize, the proposed one-parameter DIR-RT model has two levels:

$$ \begin{align*} \textit{2nd level:} ~~~& \ \theta_{i,t}=\theta_{i,t-1}+c_i(1-\rho\theta_{i,t-1})\Delta_{i,t}^{+}+w_{i,t}, \\ &\hspace{-32pt} \textit{1st level:} ~~~ \mathrm{Pr}(X_{i,t,s,l}=1\mid \theta_{i,t}, \varphi_{i,t}, \eta_{i,t,s}, a_{i,t,s}, \epsilon_{i,t,s,l}) \\ &~~~~~~~~~~~~~~~~~~~=\frac{\exp(\theta_{i,t}-a_{i,t,s}+\varphi_{i,t}+\eta_{i,t,s}+\epsilon_{i,t,s,l})}{1+\exp(\theta_{i,t}-a_{i,t,s}+\varphi_{i,t}+\eta_{i,t,s}+\epsilon_{i,t,s,l})}, \\ & \ \log(R_{i,t,s})=\mu_{i}-\nu_{i,t}+\beta L(\theta_{i,t}-a_{i,t,s})+\zeta_{i,t,s}, \end{align*} $$

where $R_{i,t,s}$ and $X_{i,t,s,l}$ are observed; $a_{i,t,s}$ ’s, $\rho $ , $\Delta _{i,t}^{+}$ ’s, and $\Delta _{i,t}$ ’s are known and $\epsilon _{i,t,s,l}\sim \mathcal {N}(0,\sigma ^2)$ with known $\sigma ^2$ . Also, we have the assumptions $\varphi _{i,t}\sim \mathcal {N}(0,\delta _{i}^{-1})$ , ${\eta }_{i,t} \sim \mathcal {N}_{S_{i,t}}(0,\tau _i^{-1}\mathbf {I}_{S_{i,t}} \mid \sum _{s=1}^{S_{i,t}}\eta _{i,t,s}=0)$ , $w_{i,t}\sim \mathcal {N}(0, \phi ^{-1}\Delta _{i,t})$ , $\zeta _{i,t,s}\sim \mathcal {N} (0,\varrho ^{-1}),$ and $\nu _{i,t}\sim \mathcal {N} (0, \kappa _i^{-1})$ . Here, $ L(\theta _{i,t}-a_{i,t,s})$ can either monotone relationship $(\theta _{i,t}-a_{i,t,s})$ or inverted U-shaped relationship $|\theta _{i,t}-a_{i,t,s}|$ .

3.1 Prior distributions for the unknown parameters

Prior choice is crucial in any Bayesian analysis. In the absence of expert’s knowledge or historical information, objective priors are used for the unknown parameters to avoid the large impacts of priors on the inference but to have some good frequentist properties (Berger, Reference Berger2006). Whenever scientific knowledge is available, we instead want to incorporate the information into the prior specifications.

Following these rules, a natural choice for the prior of one’s initial latent ability is $ \theta _{i,0}\sim \mathcal {N}(\mu _{G_{j_i}}, V_{G_{j_i}}), $ where $\mu _{G_{j_i}}$ and $V_{G_{j_i}}$ are the mean and the variance of the subpopulation (j) to which an individual i belongs. Since $c_i$ ’s in Equation (2.3) are the average growth rates (or the average learning rate in an educational context), it is often assumed to be positive, thus we specify $c_i$ as $ \pi (c_i)\propto {\mathcal {I}}(c_i\geq 0), \ \text {for all} \ i, $ where ${\mathcal {I}}(\cdot )$ is an indicator function. The speed $\kappa _i$ ’s and the random error $\varrho $ are the precision parameters in the log-normal model for response times, we follow the suggestion of Sun et al. (Reference Sun, Tsutakawa and He2001), using $ \pi (\varrho )\propto 1/\varrho ^{3/2}, $ for $\varrho $ and for all i, $\pi (\kappa _i)\propto 1/\kappa _i^{3/2}$ for $\kappa _i$ . In addition, for the prior choice of scale parameters $\delta _i$ ’s, $\tau _i$ ’s, and $\phi $ , we follow the discussion of Wang et al. (Reference Wang, Berger and Burdick2013) and use $ \pi (\phi ) \propto 1/\phi ^{3/2}, \ \pi (\delta _i)\propto 1/\delta _i^{3/2}, \ \text {and} \ \pi (\tau _i)\propto 1/\tau _i^{3/2} \ \text {for all}\ i. $ A natural objective prior for $\mu _i$ , the average response time of each individual, is a constant prior, $ \pi (\mu _i)\propto 1, \ \text {for all}\ i. $ Similarly, assume the prior of $\beta $ is $ \pi (\beta ) \propto 1. $ Intuitively, the regression coefficient $\beta $ in Equation (2.2) with a negative sign makes more sense though, we prefer to let the data determine the value and the sign of $\beta $ .

3.2 Posterior distribution and data augmentation scheme

Equation (2.1) introduces a logit link for modeling the dichotomous item responses $X_{i,t,s,l}$ , and according to Andrews & Mallows (Reference Andrews and Mallows1974), the density of a standard logistic distribution can be represented as a scale mixture of normals (Andrews & Mallows, Reference Andrews and Mallows1974; Wang et al., Reference Wang, Berger and Burdick2013). Then, applying the data augmentation idea, we can introduce a latent variable $Y_{i,t,s,l}$ for each response variable $X_{i,t,s,l}$ . Notice that in the data augmentation, $Y_{i,t,s,l}$ follows a normal distribution with mean $\theta _{i,t}-a_{i,t,s}+\varphi _{i,t}+\eta _{i,t,s}+\epsilon _{i,t,s,l}$ and variance $4\gamma _{i,t,s,l}^2$ , where $\gamma _{i,t,s,l}$ is a scale parameter assumed to have a Kolmogorov–Smirnov (K-S) distribution (i.e., $\pi (\gamma )=8\sum _{\alpha =1}^{\infty }(-1)^{(\alpha +1)}\alpha ^2\gamma \exp (-2\alpha ^2\gamma ^2)$ , $\gamma \geq 0$ . Here, we use the fact that the K-S mixture of normals will be the logistic distribution and we can verify that $ \mathrm { Pr}(X_{i,t,s,l}=1\mid \theta _{i,t}, a_{i,t,s}, \varphi _{i,t},\eta _{i,t,s},\epsilon _{i,t,s,l})=\mathrm {Pr}(Y_{i,t,s,l}>0 \mid \theta _{i,t}, a_{i,t,s}, \varphi _{i,t},\eta _{i,t,s},\epsilon _{i,t,s,l})$ . Though introducing more unknowns $Y_{i,t,s,l}$ ’s, it facilitates the MCMC computation. Then, we can rewrite the one-parameter DIR-RT models (2.1)–(2.3) as

$$ \begin{align*} \theta_{i,t}&=\theta_{i,t-1}+c_i(1-\rho\theta_{i,t-1})\Delta_{i,t}^{+}+w_{i,t}, \\ \log (R_{i,t,s})&=\mu_{i}- \nu_{i,t}+\beta L(\theta_{i,t}-a_{i,t,s})+\zeta_{i,t,s},\\ Y_{i,t,s,l}&=\theta_{i,t}-a_{i,t,s}+\varphi_{i,t}+\eta_{i,t,s}+\xi_{i,t,s,l}, \end{align*} $$

where $\xi _{i,t,s,l}\sim \mathcal {N}(0, \psi _{i,t,s,l}^{-1})$ with $\psi _{i,t,s,l}^{-1}=4\gamma _{i,t,s,l}^2+\sigma ^2$ and $\gamma _{i,t,s,l} \sim \text {K-S distribution}$ , $w_{i,t}\sim \mathcal {N}(0,\phi ^{-1}\Delta _{i, t})$ , $\varphi _{i,t}\sim \mathcal {N}(0,\delta _i^{-1})$ , ${\eta }_{i,t} \sim \mathcal {N}_{S_{i,t}}(0,\tau _i^{-1}\mathbf {I}_{S_{i,t}} \mid \sum _{s=1}^{S_{i,t}}\eta _{i,t,s}=0)$ , $\nu _{i,t}\sim \mathcal {N}(0,\kappa _{i}^{-1})$ , and $\zeta _{i,t,s}\sim \mathcal {N}(0, \varrho ^{-1})$ .

Further, we define $\theta =(\theta _1,\ldots , \theta _n)'$ with $\theta _i=(\theta _{i,0},\theta _{i,1},\ldots ,\theta _{i,T_i})'$ , $c=(c_1,\ldots , c_n)'$ , $\tau =(\tau _1,\ldots ,\tau _n)'$ , $\delta =(\delta _1, \ldots , \delta _n)'$ , $\mu =(\mu _1, \ldots , \mu _n)'$ , and $\kappa =(\kappa _1, \ldots , \kappa _n)'$ . For $l=1,\ldots ,K_{i,t,s}$ , $s=1,\ldots , S_{i,t}$ , $t=1,\ldots , T_i$ , and $i=1,\ldots , n$ , we have $Y=\{Y_{i,t,s,l}\}$ , $\gamma =\{\gamma _{i,t,s,l}\},$ and $X=\{X_{i,t,s,l}\}$ ; $\varphi =\{\varphi _{i,t}\}$ and $\nu =\{\nu _{i,t}\}$ ; $\log R= \{\log R_{i,t,s}\}, \eta =\{\eta _{i,t,s}\}$ and particularly, we use $\eta _{i,t}^*=(\eta _{i,t,1},\ldots ,\eta _{i,t, S_{i,t}-1})'$ . Then, given the data $(X, \log R)$ , the joint posterior of $(\theta , Y, c, \tau , \eta ,\varphi ,\delta , \phi , \beta , \nu , \kappa , \mu , \varrho , \gamma )$ for the proposed DIR-RT models is

$$ \begin{align*} &\pi( \theta, Y, c, \tau, \eta, \varphi, \delta, \phi,\beta, \nu, \kappa, \mu, \varrho, \gamma \mid X, \log R) \\ \propto & \left\{\prod_{i=1}^n\prod_{t=1}^{T_i}\prod_{s=1}^{S_{i,t}}\prod_{l=1}^{K_{i,t,s}} \pi(\gamma_{i,t,s,l})\right\} \left\{\prod_{i=1}^n \pi(\theta_{i,0})\pi(c_i)\pi(\delta_i)\pi(\tau_i)\pi(\kappa_i)\pi(\mu_i)\right\} \pi(\beta)\pi(\phi)\pi(\varrho) \end{align*} $$

(3.1) $$ \begin{align} \times &\left\{\prod_{i=1}^n\prod_{t=1}^{T_i}\prod_{s=1}^{S_{i,t}}\prod_{l=1}^{K_{i,t,s}} \left[{\mathcal{I}}(Y_{i,t,s,l}>0){\mathcal{I}}(X_{i,t,s,l} =1){\mathcal{I}}(Y_{i,t,s,l}\leq0){\mathcal{I}}(X_{i,t,s,l}=0)\right]\right. \nonumber \\ \times & \left.\sqrt{\frac{\psi_{i,t,s,l}}{2\pi}}\exp\left(-\frac{\psi_{i,t,s,l}(Y_{i,t,s,l}-\theta_{i,t}+a_{i,t,s}-\varphi_{i,t}-\eta_{i,t,s})^2}{2}\right){\mathcal{I}}\left(\eta_{i,t,S_{i,t}}=-\sum_{s=1}^{S_{i,t}-1}\eta_{i,t,s}\right)\right\} \nonumber \\ \times & \left\{\prod_{i=1}^n \prod_{t=1}^{T_i}\sqrt{\frac{\delta_i}{2\pi}}\exp\left(-\frac{\delta_i\varphi_{i,t}^2}{2}\right)\right\}\left\{\prod_{i=1}^n\prod_{t=1}^{T_i}\left(\frac{\tau_i}{2\pi}\right)^{(S_{i,t}-1)/2}\exp\left(-\frac{\tau_i\eta_{i,t}^{*'}\Sigma_{i,t}^{-1}\eta_{i,t}^*}{2}\right)\right\}\nonumber \\ \times& \left\{\prod_{i=1}^n \prod_{t=1}^{T_i}\sqrt{\frac{\phi}{2\pi\Delta_{i,t}}}\exp\left(-\frac{\phi[\theta_{i,t}-\theta_{i,t-1}-c_i(1-\rho\theta_{i,t-1})\Delta_{i,t}^+]^2 }{2\Delta_{i,t}}\right)\right\}\nonumber \\ \times & \prod_{i=1}^n \prod_{t=1}^{T_i} \prod_{s=1}^{S_{i,t}}\sqrt{\frac{\varrho}{2\pi}}\exp\left(-\frac{\varrho(\log(R_{i,t,s})-\mu_i+\nu_{i,t}-\beta L(\theta_{i,t}-a_{i,t,s}))^2}{2}\right)\nonumber \\ \times & \prod_{i=1}^n \prod_{t=1}^{T_i}\sqrt{\frac{\kappa_i}{2\pi}}\exp\left(-\frac{\kappa_i\nu_{i,t}^2}{2}\right), \end{align} $$

where $\Sigma _{i,t}^{-1}=\mathbf {1}_{S_{i,t}-1}\mathbf {1}_{S_{i,t}-1}'+\mathbf {I}_{S_{i,t}-1}, $ with $\mathbf {1}_{S_{i,t}-1}$ being a $(S_{i,t}-1)$ -dimensional unit vector. Note that $\pi (\theta _{i,0})$ , $\pi (c_i)$ , $\pi (\delta _i)$ , $\pi (\tau _i)$ , $\pi (\kappa _i)$ , $\pi (\mu _i)$ , $\pi (\beta )$ , $\pi (\phi )$ , and $\pi (\varrho )$ are the priors specified in Section 3.1 and $\pi (\gamma _{i,t,s,l})$ is the K-S density. To verify the posterior propriety of DIR-RT models under the objective priors, we can follow the similar steps as Wang et al. (Reference Wang, Berger and Burdick2013) did in the DIR models. The major difference in the proof is that DIR-RT models contain an additional part, which is to conjointly model with response times. Since the logarithm of response times is modeled as a normal linear mixed regression model, then the well-established facts for the posterior propriety of normal linear mixed models in Bayesian literature (cf. Sun et al., Reference Sun, Tsutakawa and He2001) can be coupled with results from Appendix C in Wang et al. (Reference Wang, Berger and Burdick2013) to show the posterior propriety of DIR-RT models.

3.3 MCMC computation of DIR-RT models

The computation is carried out by the MCMC scheme, where we sample the posterior (3.1) via block Gibbs sampling schemes. The difficulty of the sampling scheme arises in sampling the posterior distribution of latent ability $\theta _i=(\theta _{i,0}, \ldots , \theta _{i,T_i})'$ for each individual i, for $i=1, \ldots , n$ , where coordinates of $\theta _i$ are typically high dimensional and strongly correlated. When $L(x)=x$ , using the data augmentation idea, the proposed model is transformed so that $\theta _i$ (the vector) could be block sampled (see Step 2 of Appendix A.1)—within a single Gibbs sampling step conditional on the other parameters (excluding $\theta _{i,t}$ ’s)—by the highly efficient forward filtering and backward sampling algorithm (West & Harrison, Reference West and Harrison1997). However, if $L(x)=|x|$ , the computation becomes more challenging as $\theta _i$ cannot be drawn as a block. Instead, we utilize the fact that the full conditional distribution of each component of $\theta _i$ , i.e., $\theta _{i, t}$ , follows a mixture of truncated Gaussians, then $\theta _{i,t}$ ’s are drawn one at a time (cf. Step 2 of Appendix A.1).

The details of MCMC steps are given in Appendix A.1. The Gibbs sampling starts at Step 1 in Appendix A.1, with initial values $\theta ^{(0)}=\vec {0}$ , $c^{(0)}=\vec {0}$ , $\phi ^{(0)}=1$ , $\varphi ^{(0)}=\vec {0}$ , $\eta ^{(0)}=\vec {0}$ , $\delta ^{(0)}=\vec {1}$ , $\tau ^{(0)}=\vec {1}$ , $\gamma ^{(0)}=\vec {1}$ , $\mu ^{(0)}=\vec {1}$ , $\nu ^{(0)}=\vec {0}$ , and $\beta ^{(0)}=0$ in the applications and simulations, then loops through Step 15 in Appendix A.1, until the MCMC has converged. The convergence is evaluated informally by looking at trace plots.

The statistical inferences are made straightforward from the MCMC samples. For example, an estimate and $95\%$ credible interval (CI) for the latent trajectory of one’s ability $\theta _{i,t}$ can be plotted from the median, $2.5\%$ , and $97.5\%$ empirical quantiles of the corresponding MCMC realizations. In examples, ability will be graphed as a function of t, so that the dynamic changes of an examinee are apparent.

4.1 DIR-RT models simulation

Following the simulation study of DIR models in Wang et al. (Reference Wang, Berger and Burdick2013), we consider multiple individuals taking a series of tests scheduled at individually-varying and irregularly-spaced time points. Assume there are 10 individuals, each of them taking four tests on 50 different test dates, where each test contains 10 items. This specification means $K_{i,t,s}=10$ , for $s=1, \ldots , S_{i,t}$ , $t= 1,\ldots , T_i$ , $i=1, \ldots , n$ with $S_{i,t}=4$ , $T_i =50,$ and $n=10$ . The time lapse between two consecutive test dates $\Delta _{i,t} = t+10 $ if $t \le T_i / 2 $ and $ \Delta _{i,t} = t-10 $ otherwise, creating an irregularly spaced gap between two test dates.

To compare the results of DIR-RT models with DIR models, we assign same values of the parameters, $\phi $ , $\delta _i$ , $\tau _i, c_i$ as used in Wang et al. (Reference Wang, Berger and Burdick2013), where $\phi =1/0.0218^2$ , leading to standard deviation of $w_{i,t}$ in Equation (2.3) being $0.0218\sqrt {\Delta _{i,t}}$ and the values of $\delta _i$ , $\tau _i$ , and $c_i$ are specified in Table 1. For the modeling part of response times, the parameter values of $\kappa _i$ and $\mu _i$ are listed in Table 1, $\beta =-0.17$ and $\varrho =1.25$ . The parameters of DIR-RT models are chosen in such a way that they mimic the magnitude of EdSphere datasets.

Table 1 Values of unknowns used in the simulation

We consider the inverted U-shape linkage, i.e., $L(x)=|x|$ . Then, we simulate random effects and latent variables in DIR-RT models using the assigned parameter values above. Once we generated the values of $\theta _{i,t}$ using Equation (2.3), we set the test difficulties, $a_{i,t,s}$ ’s in Equation (2.1) to be $\theta _{i,t}+\zeta ^*$ , where $\zeta ^*$ is a random variable with uniform distribution on $(-0.1 , 0.1)$ . Next, the values of $\epsilon _{i,t,s,l}$ are drawn from $\mathcal {N}(0,\sigma ^2)$ with $\sigma =0.7333$ and we choose $\rho =0.1180$ . Notice that the values of $\sigma $ and $\rho $ are the same values as used in the EdSphere platform. Finally, the dichotomous data of item responses and continuous data of response times generated from the joint model are our observations, and we use the Bayesian machinery from Section 3 to estimate the model parameters of DIR-RT models. To analyze the simulation data, we utilize the priors for unknowns as specified in Section 3.1 except for $\beta $ , where we use $\pi (\beta )\propto \mathbf {1}(\beta <0)$ to make the MCMC computation converge faster.

The unknown parameters are estimated through the posterior median calculated from their corresponding MCMC samples. Each MCMC was run for 50,000 iterations with a 25,000 burn-in period. In nested or hierarchical models involving many precision parameters ( $\delta _i$ ’s, $\tau _i$ ’s, $\psi _{i,t,s,l}$ ’s, $\phi ,$ and others), a large number of burn-in iterations is often required to stabilize MCMC mixing. To run the MCMC with 50,000 iterations, it takes about 15 hours and 22 minutes by using a 3.40-GHz processor with 16-GB RAM. The MCMC convergence is assessed through a graphical examination of trace plots and through Geweke’s Diagnostics test statistics (Geweke, Reference Geweke1991). The absolute values of all parameters’ Geweke’s test statistics (in both simulation and application) are below 2, suggesting our MCMC chain has been converged. Please refer to Sections 1and 3 of the Supplementary Material for more details and additional information of auto-correlation of MCMC samples and the effective samples size are also provided there. Figures 1a–d give posterior median estimates (red squares) along with 95% CIs (red bars) of $c_i$ ’s, $\tau _i^{-1/2}$ ’s, $\delta _i^{-1/2}$ ’s, $\kappa _i^{-1/2}$ ’s, and $\mu _i$ ’s, respectively, and illustrate their true values (black dot). Clearly seen from Figure 1, the true values of parameters are contained within their corresponding $95\%$ CIs. For the posterior median estimates of parameters $\phi ^{-1/2}$ , $\varrho ^{-1/2}$ , and $\beta $ are $0.0190$ , $0.9075$ , and $-0.1815$ , respectively, with their corresponding 95% CIs being $[0.0159, 0.0229]$ , $[0.8753, 0.9427]$ , and $[-0.7028, -0.0071]$ , all of which contain their truth.

Figure 1 Posterior summary of $c_i$ ’s, $\tau _i^{-1/2}$ , $\delta _i^{-1/2}$ ’s, $\kappa _i^{-1/2}$ ’s, and $\mu _i$ ’s.

Note: The black dots represent truth, red squares are posterior median estimates, and red bars indicate $95\%$ CIs.

Next, we discuss our primary interest of estimating latent ability trajectories. Figures 2a–d illustrate four types of growth curves in our simulation, where (a) represents an individual with steady growth; (b) indicates an individual with increasing growth but nearly flat region at the end; (c) shows an individual with interrupted growth (with true ability drops in certain period); and (d) displays monotonic growth with decreasing growth rate in the middle. In Figure 2, the true ability curves (black dots) have been plotted along with our posterior median estimates of ability (blue circles) and their corresponding $95\%$ credible band (CB) (starred lines). Notice that in each subfigure, a very small proportion of true values (less than $5\%$ ) are outside of $95\%$ CBs.

Figure 2 The latent trajectory of one’s ability growth, where black dots, blue circles, and starred lines represent true ability, the posterior median estimates, and the 95% CBs, respectively.

To better assess how well the Bayesian model actually captures the truth, we can calculate the frequency of true values falling in the corresponding CIs of MCMC runs over different random trials, i.e., the frequentist coverage probability (CP) (Wang et al., Reference Wang, Berger and Burdick2013). Thus, to evaluate CPs, we conduct the simulation with the same set-up values of parameters specified earlier but generate datasets with 100 different random seeds. The CPs for $\varrho $ , $\beta $ , and $\phi $ are $96\%$ , $97\%,$ and $87\%$ , respectively, and the average CPs over all individuals for $\kappa _i$ ’s, $c_i$ ’s, $\tau _i$ ’s, $\delta _i$ ’s, and $\mu _i$ ’s are 95.4%, 92.1%, 95.4%, 93.8%, and 95%, respectively (refer to the CPs of these individual parameters in additional Table S.6 in the Supplementary Material). In addition, the average CPs for one’s ability over the study period, i.e., for $\theta _1, \ldots , \theta _{10}$ , are 94.88%, 94.90%, 95.08%, 95.62%, 95.40%, 95.54%, 95.50%, 96.22%, 95.78%, and 95.22%, respectively. Thus, while the inferential method is Bayesian, it seems to yield sets that have good frequentist coverage.

4.2 Benefits of joint modeling response times with item responses

To the best of our knowledge, the two-level DIR-RT model is the first attempt to jointly model response times and accuracy in the analysis of longitudinal data for latent traits. We are, thus, interested in knowing the benefits of introducing response times as an extra source of information. An interesting investigation is to compare the estimates of one’s ability trajectory relative to the truth with and without using response times.

Figure 3 displays the growth curve of two selected individuals, where the statistical inference is based on the simulated example in Section 4.1. For other individuals, results are similar and thus, their plots are omitted due to space limitations. In Figures 3a,b, $95\%$ CIs of DIR models (dashed red lines) encompass 95% CIs of DIR-RT models (starred blue lines), both $95\%$ CIs contain the true values (black dots). The average length of the 95% CB of ability estimates for DIR-RT models is much shorter than that of DIR models ( $0.6454$ versus $1.0370$ for the 2nd individual; and $0.6772$ versus $1.1401$ for the 6th individual, respectively). In addition, in Figure 3, both graphs of DIR-RT for ability estimates (blue circles) adhere more closely to true ability (black dots) in relative to DIR ability estimates (red dots). The average mean-squared distance between the truth and the posterior median ability estimates over time for DIR-RT models are $0.0240$ for $\theta _2$ and $0.0187$ for $\theta _6$ , in comparison to that of DIR models are $0.0711$ for $\theta _2$ , and $0.0653$ for $\theta _6$ , both of which are at least three times larger than DIR-RT models. To conclude, the results of DIR-RT models illustrate that we can largely improve the precision and remarkably reduce the bias of the estimates of one’s ability by incorporating response times.

Figure 3 The comparison of ability estimates between DIR-RT and DIR models, where black dots, blue circles, and red dots represent true ability, DIR-RT ability estimates, and DIR ability estimates, respectively; starred lines (blue) and dashed lines (red) represent 95% CBs for DIR-RT and DIR models, respectively.

5.1 DIR-RT models with inverted U-shape vs. monotone linkage

To the best of our knowledge, there has not been any empirical work to check the goodness of fit of the two linkages (i.e., monotone or inverted U-shape) for which one indeed fits the data better in conjointly modeling response times and item responses, especially when the testing data are collected at irregular and individual varying time points for a series of computerized (adaptive) testing. The lack of this research motivates us to propose some statistical measures that can be used to conduct an empirical study of EdSphere datasets to identify a better linkage and in this study, we still use the prior assigned for unknowns in Section 3.1.

Figure 4 illustrates the ability trajectories using monotone (left) and inverted U-shape (right) linkages side by side, where red dots present the posterior median trajectory of ability estimation, red dashed lines correspond to their CBs and black plus points are the “raw score,” which is a rough estimate of one’s ability obtained by solving the equation that the expectation of expected score for a person’s ability is equivalent to the observed score (Stenner, Reference Stenner2022; Swartz et al., Reference Swartz, Hanlon, Childress and Stenner2016). As a side for reference, we also include the EdSphere estimates (blue dots), the estimates employed in the design of the EdSphere learning platform, in Figure 4, which can be viewed as a preliminary version of DIR models without considering local dependence. However, EdSphere estimates use the data only available to that date. Thus, EdSphere estimates (blue dots) are much varying than ours in Figure 4.

Figure 4 The posterior summary of ability growth for the 10th individual in two linkages, where red circles, black plus, and blue dots are posterior median estimates of the ability, raw score, and EdSphere estimates, respectively, and red dashed lines are 95% CBs of our estimates.

In comparison of (a) and (b) in Figure 4, the estimates of ability growth are comparatively more robust (in an underlying increasing trend) for the inverted U-shape than that for a monotone linkage. This phenomenon is shown, for example, in the estimation of ability trajectory in Figure 4, during the period of 350–500 days and 700–800 days, where the ability trajectory estimated by an inverted U-shape linkage in the joint model is more reluctant to change its increasing trend unless there is strong information from the data (noting raw scores (black plus) in Figure 4 can be regarded as the raw data since they have the same scale as $\theta _{i,t}$ ). In the following, we have proposed two statistical methods to rigorously compare the two linkages in fitting the joint model.

5.1.1 Comparison of two linkages using Lindley’s method

The regression slope $\beta $ plays a key role in controlling the influence of the ability–difficulty distance function $L(\cdot )$ to the response time. When $\beta =0$ , it implies the distance between ability and difficulty does not affect the time an individual spends on a test and the corresponding linkage function $L(\cdot )$ can be ignored. Thus, we are interested in testing $H_0: \beta =0$ versus $H_1:\beta \neq 0$ . Lindley’s method (see Lindley Reference Lindley1965, Section 5.6), advocated by authors including Zellner (Reference Zellner1971), is an ad hoc way to test this hypothesis. According to Lindley’s method, we reject the hypothesis of $\beta =0$ at the $\alpha $ level of significance if the $100(1-\alpha )\%$ highest posterior density interval does not include $0$ . The posterior density of $\beta $ is in bell shapes (clearly seen from histograms in Figure 5). Thus, $100(1-\alpha )\%$ highest posterior density interval of $\beta $ is the same as its $100(1-\alpha )\%$ CI, but the latter one is much easier to obtain from the MCMC samples. From the table in Figure 5, it suggests $\beta =0$ cannot be rejected at $\alpha =1\%$ for monotone linkage since $99\%$ CI of $\beta $ under monotone linkage includes $0$ , while $\beta =0$ is rejected at both $\alpha =1\%$ and $\alpha =5\%$ for inverted U-shape. This indicates that the monotone linkage has a weaker correlation with response times than the inverted U-shape linkage for EdSphere datasets.

Figure 5 Posterior histogram (left) and posterior summary (right) for $\beta $ under monotone linkage and inverted U-shape, where “PM” in the table is short for “posterior median.”

5.1.2 Comparison of two linkages using deviance information criterion (DIC)

The hypothesis of $\beta =0$ can also be viewed as a Bayesian model selection question and then, we can employ DIC, one of the most popular model selection criteria in the Bayesian literature, for the comparison. Let $\boldsymbol {y}$ denote the data and $\boldsymbol {\Theta }$ indicate parameters in the model. Define ${\pi (\boldsymbol {\Theta }\mid {\boldsymbol {y}})}$ as the posterior distribution of $\boldsymbol {\Theta }$ and $\overline {\boldsymbol {\Theta }} = \mathrm {E}_{\pi (\boldsymbol {\Theta }\mid \mathbf {y})} [\boldsymbol {\Theta } \mid \boldsymbol {y}]$ represents the posterior mean of $\boldsymbol {\Theta }$ . Then, $\mathrm {Dev}(\Theta )=-2\log f({\boldsymbol {y}} \mid \boldsymbol {\Theta })$ is the deviance function for the likelihood $f({\boldsymbol {y}} \mid \boldsymbol {\Theta })$ . Following Spiegelhalter et al. (Reference Spiegelhalter, Best, Carlin and Van Der Linde2002), we have

(5.1) $$ \begin{align} & \mathrm{DIC}=\mathrm{Dev}(\overline{\boldsymbol{\Theta}})+2p_{D}, \end{align} $$

where $p_D=\mathrm {E}_{\pi (\boldsymbol {\Theta }\mid {\boldsymbol {y}})}[\mathrm {Dev}(\boldsymbol {\Theta }) \mid {\boldsymbol {y}}]-\mathrm {Dev} [\bar {\boldsymbol {\Theta }}]$ represents the effective number of model parameters. DIC is a measure to trade off between model adequacy ( $\mathrm {Dev}(\bar {\boldsymbol {\Theta }})$ part) and model complexity ( $p_D$ part), and thus the DIC can be used to compare the linkage choices in DIR-RT models. The smaller the DIC value is, the better the model does fit the data.

There are many variants of DIC, such as conditional DIC and integrated DIC. Nevertheless, recent studies have cautioned against the use of the conditional DIC, whose computation is based on the likelihood conditioning on the latent variables and is sensitive to transformations of latent variables and distributions (Millar, Reference Millar2009). Zhang et al. (Reference Zhang, Tao, Wang and Shi2019) also found that the conditional DIC often selects a model that is more complex than the true model. Instead, the integrated DIC by marginalizing out latent variables in the likelihood often performs better. Without doubt, the more the latent variables can be integrated out, the more efficiency of DIC can be achieved. Several papers showed the integrated DIC often had much smaller Monte Carlo errors compared to the conditional DICs (Celeux et al., Reference Celeux, Forbes, Robert and Titterington2006; Chan & Grant, Reference Chan and Grant2016; Merkle et al., Reference Merkle, Furr and Rabe-Hesketh2019). But usually, the analytical expressions for the integrated likelihood are hard to obtain and rather often, one resorts back to numerical integration, which is typically time-consuming.

In the proposed DIR-RT models, when $L(\cdot )$ is a monotone linkage, DIR-RT models can be written as a framework of dynamic linear models (see Step 2 of Appendix A.2) by introducing K-S random variables ( $\gamma _{i,t,s,l}$ ’s) in the data augmentation step. Then, following an analogy of Chan & Grant (Reference Chan and Grant2016) for derivations of dynamic linear models, we can show that the integrated likelihood of DIR-RT models under monotone linkage can be approximated by taking one-dimensional Monte Carlo integration over $\gamma _{i,t,s,l}$ ’s drawn from the K-S prior distribution. However, when $L(\cdot )$ is an inverted U-shape linkage, the simplification of the integrated likelihood of DIR-RT models becomes much harder since given $\gamma _{i,t,s,l}$ ’s, the DIR-RT models cannot be written as dynamic linear models. Thus, to evaluate the integrated likelihood under an inverted U-shape becomes much more computationally intensive.

However, conditioning on ability $\theta _{i,t}$ ’s, the modeling of response times part and item response part in DIR-RT models are independent. Thus, no matter what the linkage is chosen, Equations (2.1) and (2.3) of the DIR-RT models are the same under different linkages. Then, motivated by the partial DIC idea of Yao et al. (Reference Yao, Kim, Chen, Ibrahim, Shah and Lin2015) in meta-analysis, an alternative to integrated DIC for our joint DIR-RT models is to consider the integrated DIC based on the part of response times only. Based on this partial DIC derived (see details in Appendix B), the inverted U-shape linkage turns out to fit the EdSphere data better, where the partial DIC for the inverted U-shape is 5661.4 in comparison to that of the monotone linkage, which is 5775.3.

5.2 Retrospective estimation of ability growth under inverted U-shaped linkage

Both Lindley’s method and partial DIC criterion support the choice of an inverted U-shape linkage for the analysis of EdSphere data using DIR-RT models. Therefore, to investigate the first two goals mentioned in the beginning of Section 5, we are going to use an inverted U-shape throughout the rest of the article. After we run MCMC for the EdSphere datasets using the inverted U-shape, we have conducted a posterior predictive check (Gelman et al., Reference Gelman, Meng and Stern1996) using a mean test statistic for the response time and its p-value is $0.4978$ , which shows the inverted U-shape indeed fit the response time well.

Figure 6 presents a retrospective analysis of the reading ability for 3rd, 12th, 18th, and 23rd individuals (using all data recorded for each individual during the study period). In Figure 6, the red circles are the posterior median estimates of one’s ability, the red dashed lines correspond to the 2.5% and 97.5% quantiles of the posterior distributions of the abilities, and the black plus points are raw scores. Similar to Wang et al. (Reference Wang, Berger and Burdick2013), we find all these growth trajectories have an overall increasing trend but such kind of growth can be interrupted. In particular, when there is a large time gap between subsequent tests, the ability appears to drop for some individuals, which is clearly seen from Figure 6. A natural explanation might be that during vacations, students do not read and could actually lose ability or they become less familiar with computerized tests after a long break.

Figure 6 The posterior summary of the ability growth for $\theta _3$ , $\theta _{12}$ , $\theta _{18,}$ and $\theta _{23}$ , where red circles, black plus, and blue dots are posterior median estimates of the ability, raw score, and EdSphere estimates, respectively, and red dashed lines represent 95% CBs of our estimates.

In Figure 7, we summarize the posterior median (red square) and 95% CI (red bars at two ends) of the average growth rates $c_i$ ’s, the standard deviations of test random effects $\tau _i^{-1/2}$ ’s, the standard deviations of the daily random effects $\delta _i^{-1/2}$ ’s, the standard deviations of speediness, $\kappa _i^{-1/2}$ ’s, and the average response time for each individual, $\mu _i$ , for $i = 1, \ldots , 25$ . Moreover, the estimated posterior median of $\phi ^{-1/2}$ is 0.0708 and its 95% CI is $[0.0608, 0.0831]$ and the result of $\beta $ is summarized in the table of Figure 5.

Figure 7 The posterior summary of c, $\tau _i^{-1/2}$ ’s, $\delta _i^{-1/2}$ ’s, $\kappa _i^{-1/2}$ ’s, and $\mu _i$ ’s.

Figures 7b,c show that the standard deviations of test and daily random effects (i.e., $\eta _{i,t,s}$ ’s and $\varphi _{i,t}$ ’s) are almost all quite large with 95% CIs that are well separated from zero except 22nd and 23rd individuals. Recall that these random effects were included in the model to account for a possible lack of the local independence; the evidence is thus strong that the local independence is, indeed, not tenable for this data and that both types of random effects are presented for most individuals. Similarly, Figure 7d illustrates that speediness of individuals is different on the daily basis except that of the 22nd individual (whose speed is almost steady during the studying period). Moreover, there are clearly some patterns in the variation of speediness across individuals; as for some of them, the difference of their speediness on a daily basis is more crucial than that of the others. The differentiation of the average response time in Figure 7e suggests that some individuals incline to take a longer time to finish a test than others no matter what the difficulty level of the test is. As well, it is not surprising the average growth rates are quite different among individuals, as shown in Figure 7a.