R-package LNIRT for joint modeling of response accuracy and times

Level 1: measurement models

Level 2 and 3: structural person and item parameter models

Parameter estimation and model identification

A summarized description of the estimation method and model identification rules is presented. A more detailed description of the parameter estimation can be found in Fox, Klein Entink & van der Linden (2007), Klein Entink, Fox & van der Linden (2009) and Fox (2010). IRT models are usually identified by fixing the mean and variance of the latent variable to zero and one, respectively. Typically, this can be done directly by restricting the prior mean ${\mu }_{\theta }$ and variance ${\sigma }_{\theta }^2$, or by putting restrictions on the item parameters. The joint model can be identified in the same way. However, for the identification rules in the package LNIRT restricting the variance of a person parameter has been avoided. When restricting the variance of a (random) person parameter, the covariance matrix in Eq. (6) is also restricted, and the inverse-Wishart distribution does not apply to a restricted covariance matrix. For this scenario, Klein Entink, Fox & van der Linden (2009) redefined the prior for the person parameters, where $\rho$ becomes a regression parameter in the regression of working speed on ability with the working speed variance included in the error variance. However, this identification procedure would also complicate other modeling features (e.g., model-fit tools, variable working speed).

The LNIRT package provides two options to identify the model. For both options, the variance of the latent scales are identified by restricting the product of discriminations and of the time discriminations to one, ${\prod }_k{a}_k=1$ and ${\prod }_k{\varphi }_k=1$, respectively. For option one (referred to as (R-code) ident=1), the mean of the scales are identified by restricting the sum of the difficulty and of the time-intensity parameters to zero, ${\sum }_k{b}_k=0$ and ${\sum }_k{\lambda }_k=0$, respectively. For option two (referred to as (R-code) ident=2), the mean of the scales are identified by fixing the mean of the ability parameter to zero, ${\mu }_{\theta }=0$, and of the speed parameter to zero, ${\mu }_{\zeta }=0$.

Model parameters are estimated through Gibbs sampling from their joint posterior distribution. The procedure involves the division of all unknown parameters into blocks, with iterative sampling of the conditional posterior distributions of the parameters in each block given the preceding draws for the parameters in all other blocks (Fox, 2010). In various simulation studies and real data analysis, the Gibbs samplers in the LNIRT package showed good convergence properties and generally returned efficient MCMC samples (e.g., Fox & Marianti, 2016, 2017; Klein Entink, Fox & van der Linden, 2009).

Logistic vs probit model item parameter estimates

The LNIRT package uses the normal-ogive IRT model (Probit model) to model the probability of a correct/positive response. The item parameter estimates under the normal ogive (Probit) model can be transformed to a Logistic scale and vice versa. Therefore, the logistic scale factor is used (1.7) to transform parameters on a Logistic scale to those on a Probit scale. Let ${\mu }_L$, ${\sigma }_L$ and ${\mu }_P$, ${\sigma }_P$ denote the mean and variance of the latent scale under the Logistic and Probit model, respectively. Then, the item parameters $a_k$ and $b_k$ are invariant under both scales, when applying the logistic scale factor

(10) $$\underset{ProbitModel}{\underbrace{a_k\frac{\left({\theta }_i-{\mu }_P\right)}{{\sigma }_P}-b_k}}=\underset{LogisticModel}{\underbrace{1.7\left(a_k\frac{\left({\theta }_i-{\mu }_L\right)}{{\sigma }_L}-b_k\right)}}.$$

Assume that data were generated under the Logistic IRT model, and item parameter estimates were obtained under the Probit model. The Probit model estimates can be transformed to the Logistic model estimates;

$$\frac{{\hat{a}}_k}{{\sigma }_P}=\frac{1.7a_k}{{\sigma }_L}$$and, it follows that

$$\frac{\left(\frac{{\hat{a}}_k}{{\sigma }_P}\right)}{1.7}=\frac{a_k}{{\sigma }_L}.$$

In the same way, the difficulty parameters under the Probit model can be transformed. When assuming that ${\mu }_P$ and ${\mu }_L$ are zero, the transformation is

$$\begin{array}{rl}{\hat{b}}_k& =b_k\\ 1.7{\hat{b}}_k& =b_k^{\ast },\end{array}$$since $b_k^{\ast }=1.7b_k$ is the difficulty parameter under the Logistic model.

Person fit

The general idea is (1) to define a person-fit statistic for an RA and RT pattern, (2) to define classification variables, as a function of the statistics, to represent a non-aberrant and an aberrant state, and (3) to compute the posterior probability of the aberrant state. A person-fit test for RA and RT patterns is discussed. For each test, a dichotomous classification variable is introduced, which states whether a pattern is considered extreme. Then, the person-fit test for the joint model is constructed from the dichotomous classification variables for RA and RT patterns.

Item fit

Without reproducing the equations for item-fit statistics, the person-fit tests for RA and RT patterns can be modified to examine the fit of RA and RT item patterns. Therefore, the log-likelihood of RA and of RTs of an item is considered to define an item-fit statistic for RA and RTs, respectively. The estimated posterior probability of each item-fit statistic represents the probability that the statistic is extreme under the model, which means that the pattern of responses to the item is extreme.

Residual analysis

Bayesian residual computation has been considered by Albert & Chib (1993), Johnson & Albert (2006), and Fox (2010, Chapter 5) to evaluate the fit of an IRT model. This approach is extended to the joint model, and latent residuals $e_{ik}$ are defined, which represent the difference between a latent continuous RA and the mean component. The latent continuous RA is defined through a data augmentation method to facilitate a Gibbs sampling algorithm (Fox, 2010, Chapter 3 and 4). Conditional expectation of a latent residual is derived by integrating out the augmented latent variable to obtain a Rao-Blackwell estimator for the latent residuals. The estimated latent residuals are used to quantify the total percentage of outliers per item and per test taker.

For the RA, the conditional expected latent residual is given by

(18) $$E\left(e_{ik}\mid a_k,b_k,{\theta }_i\right)=\{\begin{array}{cc}\frac{-\phi \left(b_k-a_k{\theta }_i\right)}{\mathrm{\Phi }\left(b_k-a_k{\theta }_i\right)}\hfill & y_{ik}=0\hfill \\ \frac{\phi \left(b_k-a_k{\theta }_i\right)}{\mathrm{\Phi }\left(b_k-a_k{\theta }_i\right)}\hfill & y_{ik}=1\hfill \end{array}$$where $\phi$ and $\mathrm{\Phi }$ are the normal density function and the cumulative distribution function, respectively. Subsequently, the posterior probability that a latent residual is greater than a threshold C is given by

(19) $$P\left(\left|e_{ik}\right|>C\mid a_k,b_k,{\theta }_i\right)=\{\begin{array}{cc}\frac{\mathrm{\Phi }\left(C\right)}{1-\mathrm{\Phi }\left(a_k{\theta }_i-b_k\right)}\hfill & y_{ik}=0\hfill \\ \frac{\mathrm{\Phi }\left(-C\right)}{\mathrm{\Phi }\left(a_k{\theta }_i-b_k\right)}\hfill & y_{ik}=1\hfill \end{array}$$

The residuals for the RTs, referred to as ${ϵ}_{ik}$, can be estimated directly as the difference between the $R{T}_{ik}$ and the mean, ${ϵ}_{ik}=(r{t}_{ik}-({\lambda }_k-{\varphi }_k{\zeta }_i))$. The extremeness of an RT residual is expressed as the posterior probability that the residual is greater than a threshold C. This posterior probability can be expressed as

(20) $$P\left(\left|{ϵ}_{ik}\right|>C\mid {\zeta }_i,{\lambda }_k,{\varphi }_k,r{t}_{ik}\right)=\mathrm{\Phi }\left(-C-\frac{{ϵ}_{ik}}{{\sigma }_k}\right)+1-\mathrm{\Phi }\left(C-\frac{{ϵ}_{ik}}{{\sigma }_k}\right).$$

Differential Working Speed

The generalization of the joint model to allow differential working speed has been introduced by Fox & Marianti (2016). A summarized version of their procedure is given, which shows how the joint model can be modified to include differential working speed.

Differential working speed is defined as a change in working speed of a test taker during the test. This relaxes a basic assumption of the joint model to work with a constant speed throughout the test. The change in working speed is modeled using a latent growth model. A random intercept, a linear trend component, and a quadratic time component are considered to model the speed trajectory of each test taker. The random intercept, ${\zeta }_{i0}$, represents the initial value of working speed. The linear trend component, ${\zeta }_{i1}$, is used to model a linear change in speed. Test takers can start slowly (fast) to increase (decrease) their speed later on. The quadratic time component, ${\zeta }_{i2}$, is used to decelerate or accelerate the linear trend. For instance, a positive linear trend in working speed can be decelerated by a negative quadratic component.

The log-normal differential speed model with a random trend and quadratic time variable is represented by

(23) $$\begin{array}{rl}R{T}_{ik}& ={\lambda }_k-{\varphi }_k\left({\zeta }_{i0}+{\zeta }_{1i}{X}_{ik}+{\zeta }_{2i}{X}_{ik}^2\right)+{ϵ}_{ik}\\ {\zeta }_{i0}& \sim N\left(0,{\sigma }_{{\zeta }_0}^2\right)\\ {\zeta }_{i1}& \sim N\left(0,{\sigma }_{{\zeta }_1}^2\right)\\ {\zeta }_{i2}& \sim N\left(0,{\sigma }_{{\zeta }_2}^2\right),\end{array}$$where the time variable $X_{ik}$ represents the order in which the test items are solved. The beginning of the test is represented by the time-point $X_{i1}=0$ for each test taker $i$. This does not mean that each test taker should start the test at the same time. It is simply a reference point for the speed process of each test taker. The time scale is defined on an equidistant scale to reduce the MCMC computations. Let ${\mathbf{X}}_{(i)}=X_{(i1)},X_{(i2)},\dots ,X_{(iK)}$ represent the order in which items are solved by test taker $i$. Then, a convenient time scale is defined by $X_{ik}=(X_{(ik)}-1)/K$—the times are defined on a scale from 0 to 1, with 1, the upper bound, representing an infinite number of items. The time scale addresses the order in which the items are made and assumes an equidistant property of the speed measurements.

The average of the time intensities defines the average time to complete the test, since the random effects have a population mean of zero. The random intercept represents the speed measured at the first item. The population-average speed trajectory is constant, and shows no change in speed. Test takers can work faster (slower) than this population-average level, which corresponds to a positive (negative) initial speed level. A negative (positive) growth rate shows a decrease (increase) in speed, which can be decelerated (accelerated) by the quadratic time component. The random component variances represent the variance in growth parameters of the working speed trajectories. The covariance between random speed effects is modeled through the time-discrimination parameters. When the time-discrimination parameters are all fixed to one, the full covariance matrix of the random speed components is freely estimated.

Input

The following arguments are mandatory for the LNIRT function:

• RT: A matrix RT containing the log-normally transformed response times in wide format, (N) persons in rows and (K) items in columns.

• Y: A matrix Y containing the (binary) RA data in wide format, (N) persons in rows and (K) items in columns. The binary RA data is coded as zero (incorrect) and one (correct) (missing values as NA).

• data (optional if RT and Y are given): A list or a simLNIRT object (output object of the function simLNIRT) containing the RT and RA matrices and optionally predictors for the item and person parameters. If a simLNIRT object is provided, in the summary the simulated item and time parameters are shown alongside of the estimates. If the required variables cannot be found in the list, or if no data object is given, then the variables are taken from the environment from which LNIRT is called.

• XG: The integer number of MCMC iterations (the default is 1,000), this includes the burn-in period.

The remaining arguments are optional for the LNIRT function:

• burnin: The percentage of the total number of MCMC iterations (XG) which will serve as the burn-in period of the chains. The default is 10%.

• ident: Identification rule, (ident=1) restrict sum of difficulties and sum of time intensities to zero and (ident=2) restrict mean person parameters to zero. The default is (ident=2). The product of (time) discriminations is restricted to one.

• guess: A logical variable with the default FALSE, where TRUE (FALSE) represents (not) a guessing parameter in the IRT model.

• par1: A logical variable with the default FALSE, where TRUE represents the bracket notation for the mean term of the IRT component as in Eq. (2), and FALSE the non-bracket notation as in Eq. (1). In general, the MCMC performance is better for the non-bracket parameterization.

• XGresid: the number of MCMC iterations (default is 1,000) to be done before starting the residual computation.

• residual: A logical variable with the default FALSE. When TRUE, a complete residual analysis is done together with the estimation of the joint model parameters, as discussed in “Model Fit”. The residual computations are started after XGresid MCMC draws. Therefore, XG should be greater than XGresid, and a sufficient number of MCMC iterations should be made after XGresid MCMC iterations to obtain accurate residual estimates. Preferably, at least 5,000 MCMC iterations are made, when doing a residual analysis.

• td: A logical variable with the default TRUE. When TRUE, the time-discrimination parameter is estimated. When FALSE, the time discrimination is restricted to one.

• WL: A logical variable with the default FALSE. When TRUE, the time-discrimination parameter represents the inverse of the measurement error variance parameter according to the parameterization of van der Linden (2007).

• alpha: An optional vector of length K of pre-defined item discrimination parameters.

• beta: An optional vector of length K of pre-defined item difficulty parameters.

• phi: An optional vector of length K of pre-defined time-discrimination parameters.

• lambda: An optional vector of length K of pre-defined time-intensity parameters.

• XPA: An optional matrix of predictor variables for the ability parameters, where the columns represent the predictor variables. Categorical predictor variables need to be dummy coded.

• XPT: An optional matrix of predictor variables for the speed parameters, where the columns represent the predictor variables. Categorical predictor variables need to be dummy coded.

• XIA: An optional matrix of predictor variables for the item difficulty parameters, where the columns represent the predictor variables. Categorical predictor variables need to be dummy coded.

• XIT: An optional matrix of predictor variables for the time-intensity parameters, where the columns represent the predictor variables. Categorical predictor variables need to be dummy coded.

• MBDY: An optional missing-by-design indicator matrix—of the same size as Y—for missing values (coded NA) in the Y matrix due to the test design (0 is missing by design, 1 is not missing by design). Multiple imputations are simulated for missing values (missing at random) in the Y matrix which are not assigned in MBDY as missing-by-design.

• MBDT: An optional missing-by-design indicator matrix—of the same size as RT—for missing values (coded NA) in the RT matrix due to the test design (0 is missing by design, 1 is not missing by design). Multiple imputations are simulated for missing values (missing at random) in the RT matrix which are not assigned in MBDT as missing-by-design.

Output

The LNIRT function creates an object of class LNIRT, which stores the MCMC output, posterior draws and posterior mean estimates. The output of the function LNIRT is described in an itemized way, without including a residual computation.

• data: If available a data object from the function simLNIRT.

• burnin: Percentage from XG representing the burn-in MCMC iterations.

• ident: Same as the input variable ident.

• guess: Same as the input variable guess.

• MAB: The MCMC sampled values for the item parameters (discrimination, difficulty, time discrimination, time intensity), object is an array of dimension XG (number of MCMC iterations) by K (number of items) by 4 (number of item parameters).

• MCMC.Samples: This object contains the sampled MCMC values from LNIRT, where the following objects are stored:

• – Cov.Person.Ability.Speed: Samples of covariance of ability and speed.

• – CovMat.Item: Array of dimension XG (MCMC iterations) by K (items) by 2 (time-discrimination and time-intensity parameters), it contains the sampled variance of time discrimination and time intensity as well as their sampled covariance.

• – Item.Dificulty: Samples of difficulty parameters (XG by K).

• – Item.Discrimination: Samples of discrimination parameters (XG by K).

• – Item.Guessing: Samples of guessing parameters (XG by K).

• – Mu.Item.Difficulty: Sampled values of mean item difficulty parameter.

• – Mu.Item.Discrimination: Sampled values of mean item discrimination parameter.

• – Mu.Person.Ability: Sampled values of mean ability parameter.

• – Mu.Person.Speed: Sampled values of mean speed parameter.

• – Mu.Time.Discrimination: Sampled values of mean time-discrimination parameter.

• – Mu.Time.Intensity: Sampled values of mean time-intensity parameter.

• – Person.Ability: Sampled values of person’s ability parameter (XG by N).

• – Person.Speed: Sampled values of person’s speed parameter (XG by N).

• – Sigma2: Sampled values of measurement error variance parameters (XG by K).

• – Time.Discrimination: Sampled time-discrimination parameters (XG by K).

• – Time.Intensity: Sampled time-intensity parameters (XG by K).

• – Var.Person.Ability: Sampled variance ability parameters.

• – Var.Person.Speed: Sampled variance speed parameters.

• Mguess: Sampled values for the guessing parameter (under the default Beta prior B(20,80)).

• MmuI: Sampled values of the mean discrimination, difficulty, time discrimination, and mean time intensity (in that order).

• MmuP: Sampled values of the mean ability and mean speed parameters.

• MSI: Sampled values of the covariance matrix of item parameters. Array of dimension XG by 4 by 4 (discrimination, difficulty, time discrimination, time intensity), see also CovMat.Item.

• Msigma2: Sampled values of measurement error variance parameters.

• MSP: Sampled values of the covariance matrix of person parameters (ability and speed), in an array of dimension XG by 2 by 2.

• Mtheta: Posterior mean estimates of ability and speed of dimension N by 2.

• MTSD: Posterior standard deviation of ability and speed of dimension N by 2.

• par1: Same as the input variable par1.

• Post.Means: Posterior mean estimates of the following parameters:

• – Cov.Person.Ability.Speed: Covariance ability and speed.

• – CovMat.Item: Covariance matrix item parameters.

• – Item.Difficulty: Item difficulty estimates.

• – Item.Discrimination: Item discrimination estimates.

• – Mu.Item.Difficulty: Mean item difficulty estimate.

• – Mu.Item.Discrimination: Mean item discrimination estimate.

• – Mu.Person.Ability: Mean ability estimate.

• – Mu.Person.Speed: Mean speed parameter.

• – Mu.Time.Discrimination: Mean time discrimination.

• – Mu.Time.Intensity: Mean time intensity.

• – Person.Ability: Ability parameters.

• – Person.Speed: Speed parameters.

• – Sigma2: Measurement error variance.

• – Time.Discrimination: Time-discrimination parameters.

• – Time.Intensity: Time-intensity parameters.

• – Var.Person.Ability: Variance ability parameter.

• – Var.Person.Speed: Variance speed parameter.

• RT: Logarithm of the RT data.

• td: Same as the input variable td.

• WL: Same as the input variable WL.

• XIA: Same as the input variable XIA.

• XIT: Same as the input variable XIT.

• XPA: Same as the input variable XPA.

• XPT: Same as the input variable XPT.

• Y: Same as the input variable Y.

When the residual computation is included (residual=TRUE), then the LNIRT object includes the following additional output variables:

• EAPCP1: Posterior probability that the RT pattern is flagged as aberrant according to Eq. (16), using the posterior probability that the person-fit statistic $l^t$ is extreme as defined in Eq. (15) with a significance level of 0.05.

• EAPCP2: Posterior probability that the RA pattern is flagged as aberrant according to Eq. (13), using the posterior probability that the person-fit statistic, $l_s^y$, is significant (with a significance level of 0.05).

• EAPCP3: Posterior probability that both patterns (RA and RT) are flagged as aberrant according to Eq. (17).

• EAPKS: Posterior probability of an extreme KS-test result according to Eq. (22), which indicates that the RT residuals are not normally distributed.

• EAPKSA: A significance probability of an extreme KS-test that the latent residuals of RA items are not normally distributed. This significance test has no power.

• EAPresid: Posterior probability of an extreme (standardized) residual (RT data), which is greater than plus or minus two in absolute value.

• EAPresidA: Posterior probability of an extreme (standardized) latent residual (RA data), which is greater than plus or minus two in absolute value.

• IFl: The (negative) log-likelihood statistic for item fit under the IRT model (high values represent misfit).

• IFlp: The posterior significance probability of an extreme item fit under the IRT model.

• lZI: Item-fit statistic representing the posterior probability that the squared sum of item residuals are extreme under the model.

• EAPl0: The log-likelihood contribution of each RA observation, where a low value represents a misfit.

• PFl: The standardized (negative) log-likelihood contribution of each RA pattern, where a high value represents a misfit.

• PFlp: Posterior (significance) probability of observing a more extreme person-fit statistic for the RA pattern than the observed one.

• lZPA: Posterior significance probability of an extreme person-fit test for RA pattern based on latent residuals. This significance test has no power (under construction).

• lZPT: The (unstandardized) estimated person-fit statistic according to Eq. (14).

• lZP: Posterior (significance) probability of observing a more extreme person-fit statistic for the RT pattern than the observed one, according to Eq. (15).

MCMC output analysis

The posterior draws in the output of LNIRT are obtained through a Gibbs sampler. MCMC convergence tests are necessary to make sure that the MCMC chains converged before making statistical inferences. The coda package (Plummer et al., 2006) can be used to do an MCMC convergence analysis. The LNIRT function returns a single MCMC chain for each model parameter, which can be directly translated to an MCMC object with the function as.mcmc. The single-chain convergence diagnostics, for instance Geweke, and Heidelberg-Welch, can be used to examine if the chains did not converge. The multiple-chain convergence diagnostics, for instance Gelman and Rubin’s convergence statistic, requires multiple calls to LNIRT to create multiple MCMC chains with different starting values—LNIRT uses random starting values when parameter values are not pre-specified.

With the summary function of the package LNIRT, posterior means and standard deviations are computed and reported. Although the LNIRT Gibbs samplers produce efficient MCMC samples with low autocorrelation, it is possible that for a specific dataset some of the chains have a medium to high autocorrelation. Then, the reported standard deviation estimates may underestimate the standard deviation, since autocorrelation in the chains are ignored. The coda and the mcmcse package can be used to compute the autocorrelation and the Monte Carlo standard error. Note that the autocorrelation has no effect on the posterior mean estimate. The effective sample size—the sample size of independently distributed values with the same variance as the autocorrelated MCMC sample—can also be computed to examine if the run was long enough to make accurate and reliable inferences. A reasonable rule of thumb is to have an effective sample size of 400. Then, the Monte Carlo standard error is less than 5% of the overall uncertainty of the posterior mean.

The default burn-in period is 10% of the total number of MCMC iterations. This burn-in period is used in the computation of the posterior estimates, which are also reported with the summary function. Extensive simulation studies showed that the burn-in period is usually below 100 MCMC iterations, but this could be higher for a specific data set. Furthermore, the MCMC properties are better for the RA model in Eq. (1) than for the RA model with the additional brackets given in Eq. (2).

The credential (Form 1) data

The credential data set of Cizek & Wollack (2016) is analyzed to illustrate the functionality of the package LNIRT. The credential data set concerns 1,636 test takers who applied for licensure. Form 1 of the test was administered, which consisted of 170 licensure exam items, and 30 pretest items. A total of 10 background variables of the test takers was stored—this includes the country where the candidate received his/her educational training, the state in which the test taker applied for licensure, and the center where the candidate took the exam. The RA and RT data were also stored. The collected data followed from a year of testing using a computer-based program that tests continuously.

Each candidate completed one of the three pretest forms, each consisting of 10 items. In Table 1, the test design is given, which shows the incomplete test design of in total 200 items for three groups.

The RA and RT data were extracted from the data to be used in the LNIRT() function. The RA data consisted of 1,636 test takers (in rows) and 170 exam items (in columns).

 R> Y<- as.matrix(data[c(which(colnames(data)=="iraw.1")

 +:which(colnames(Cdata)=="iraw.170"))])

 R> head(Y[,1:5],3)

      iraw.1 iraw.2 iraw.3 iraw.4 iraw.5

 [1,]     1    1    1    0    1        

 [2,]     1    0    1    0    0

 [3,]     0    0    0    0    1    

The RT data also consisted of 1,636 test-takers and 170 exam items. The RTs were transformed to a logarithmic scale. A total of 105 RTs were equal to zero. Possibly the item was skipped and a default incorrect response was recorded, since the corresponding RA were all incorrect. The zero RTs were converted into NAs, since it was unknown why a zero RT was recorded. The minimum RT was 2 s, and the 100 highest RTs ranged from 6 to 12 min.

 R> RT<-as.matrix(data[c(which(colnames(data)=="idur.1"):

 + which(colnames(data)=="idur.170"))])

 R> RT[RT==0]<-NA

 R> RT<-log(RT)

 R> head(RT[,1:5],3)

     idur.1  idur.2  idur.3  idur.4  idur.5

 [1,] 4.094345 3.555348 3.555348 3.465736 3.295837

 [2,] 4.007333 4.060443 4.343805 3.850148 4.110874

 [3,] 4.234107 3.761200 4.007333 3.178054 4.127134

In the first analysis, the (default) joint model was fitted to the data to explore the item parameter estimates and the item and person covariance matrix. The model was identified by restricting the population means of ability and speed to zero and by restricting the product of time discriminations and discriminations to one. A total of 5,000 MCMC iterations were computed, and the burn-in period was 500 (i.e., 10% of the total number of MCMC iterations).

 R> library(LNIRT)

 R> out0 <- LNIRT(RT=RT, Y=Y, XG=5000,ident=2,burnin=10)

 R> summary(out0)

Output analysis

The general summary function of LNIRT shows the item number (which corresponds to the column number of Y and of RT), and the posterior mean (labeled EAP) and standard deviation (labeled SD) estimates of the item and of the prior parameters of the items and persons. In Fig. 1, the estimation results for the first five item parameters are given.

In Fig. 2, the item and person prior parameter estimates are given. The mean of the items, ${\mu }_I$, is labeled in the output as mu_a (mean discrimination), mu_b (mean difficulty), mu_phi (mean time discrimination) and mu_lam (mean time intensity). The posterior mean estimate of the covariance matrix of the items, ${\mathbf{\text{Σ}}}_I$, is given under the label Sigma_I. The estimated mean item discrimination is around 1.19 with a variance of around 0.32. Around 10 items have an estimated item discrimination below 0.40. The item difficulty estimates ranges from −1.84 to.75, with a mean of −0.70 and a variance of 0.27. The estimated average time-discrimination and time-intensity is around 1.03 and 3.96, with a variance of 0.05 and 0.11, respectively. The time discrimination ranged from 0.47 to 1.42, and it appears that the items show better discriminating performance in speed than in ability. However, the estimated working-speed variance is very small and around 0.03. The variance in RTs is mostly explained by the time-intensity parameters and hardly by differences in working speed across test takers. Thus, the time discriminations are somewhat higher than the item discriminations, but the time discriminations have a minor contribution in explaining variance in RTs.

The covariance matrix of the item parameters shows that the less difficult items are more discriminating (correlation of −0.43), and the less time-intensive items are also more time-discriminating (correlation of −0.40). Differences in ability and speed are better measured with less difficult and less time-intensive items. The more difficult items are also more time-intensive (correlation of 0.46). The time-discriminating items are positively correlated with the item discriminations (correlation of 0.49).

The covariance matrix of the person parameters (ability and speed), ${\mathbf{\text{Σ}}}_P$, is given under the label Sigma_P. The variance in ability and speed across test takers are both small. There is a positive correlation between ability and speed of around 0.40, which states that more able test takers also work faster than the less able ones.

Model fit

The joint model-fit tools are illustrated by re-running the analysis with explanatory variables for the 170 items and activating the residual computation (R-code residual=TRUE).

 R> out1 <- LNIRT(RT=RT,Y=Y,XG=5000,XPA=XA,XPT=XT,residual=TRUE)

 R> summary(out1)

The summary report contains an overview of the extreme residuals under the header Residual Analysis. A total of 0.07% of RTs residuals were considered extreme with at least 95% posterior probability (Eq. (20)). This concerns RTs that were small and around 2–6 s or much higher than the item-average RTs. For the RA residuals, around 0.02% was estimated to be extreme with 95% posterior probability (Eq. (19)). This concerns incorrect responses from test takers with an above-average ability. For 62 items (36.5%) the assumption of log-normally distributed residuals was violated for which a significant probability of the KS test was computed (Eq. (22)). The variance in working speed and time intensities are small, and the estimated residual variance is around 0.26. Therefore, RT outliers more easily affect the fit of the log-normal distribution. The item-fit statistics (RA and RT data) did not identify a significant misfit of an item. Despite the outliers, log-likelihoods of item patterns were not significant under the joint model.

The EAPCP1 is reported and around 18.34%, representing the percentage of RT patterns which are considered extreme with 95% posterior probability. The percentage of significant extreme RT patterns is around 19.5%, when using a significance level of 0.05. The reported EAPCP2 is around 1.59% and the significant RA patterns is around 1.65%, which shows that there are only a few RA patterns identified as extreme. Finally, around 0.31% of the joint patterns (RA and RT) are extreme (object EAPCP3).

The heterogeneity in person-fit statistics for RA an RT patterns is examined with a linear regression of the statistics on the number of test attempts, country, and pretest group. The country variable was (dummy) recoded (US (Cgroup1=1),Philippines (Cgroup2=1),India (Cgroup3=1),Others (intercept)). The pretest groups were already represented by two dummy coded variables.

  R> summary(lm(out1$PFl ~ as.factor(XFT$Attempt)+(XFT$Cgroup)+(XFT$Pgroup)))

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -0.2042450.045153-4.5236.53e-06 ***

as.factor(XFT$Attempt)2-0.0043230.060106 -0.0720.94267

as.factor(XFT$Attempt)3-0.1024010.079525 -1.2880.19805

as.factor(XFT$Attempt)4-0.2534210.114568 -2.2120.02711 *

as.factor(XFT$Attempt)4+0.0176700.080984 0.2180.82731

XFT$Cgroup10.0518250.046087 1.1240.26097

XFT$Cgroup2-0.2159070.054834 -3.9378.58e-05 ***

XFT$Cgroup3-0.1094580.053455 -2.0480.04075 *

XFT$Pgroup10.0117850.034349 0.343 0.73157

XFT$Pgroup20.0953900.029290 3.257 0.00115 **

It follows that the RA patterns of test takers with more test attempts are less extreme than those with a fewer attempts. Test takers from the Philippines and India are less likely to have an aberrant RA pattern. Furthermore, test takers from the third pretest group (variable Pgroup2) have higher person-fit statistic scores. They correspond to RA patterns that are less likely to be observed under the model than RA patterns with lower person-fit statistics. For the person-fit test for RA patterns, for a significance level of 0.05 the critical value is 1.645. The average-statistic scores (number of attempts, country, pretest group) are much lower than the critical value, and test takers with an aberrant RA pattern are also outliers in their groups.

The heterogeneity in person-fit statistic scores for the RT patterns is also explored through a linear regression with the same explanatory variables.

 R> summary(lm(out1$lZPT ~ as.factor(XFT$Attempt)+(XFT$Cgroup)+(XFT$Pgroup)))

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 169.92313.085555.071< 2e-16 ***

as.factor(XFT$Attempt)21.15064.10730.2800.7794

as.factor(XFT$Attempt)310.24105.43431.8850.0597.

as.factor(XFT$Attempt)42.69487.82900.3440.7307

as.factor(XFT$Attempt)42.32465.53410.4200.6745

XFT$Cgroup1-7.75633.1494-2.4630.0139 *

XFT$Cgroup23.01383.74710.8040.4213

XFT$Cgroup315.40763.65284.2182.6e-05 ***

XFT$Pgroup12.51522.34731.0720.2841

XFT$Pgroup20.26812.00160.1340.8935

For the $l_i^t$ statistic, the critical value is 201.4 when the significance level is 0.05 ( $l_i^t$ is chi-square distributed with 170 degrees of freedom). The intercept corresponds to test takers from pretest group 1 from citizens outside the US, India, and the Philippines, who took the test for the first time. For test takers with more attempts the person-fit statistic is on average higher and closer to the critical value. Those test takers are more likely to give a very fast or much slower response. Test takers from the US have on average much lower person-fit scores, and test takers from India have a much higher person-fit score.

In Fig. 3, the person-fit scores for RA patterns are plotted against those of RT patterns to provide a more comprehensive overview of the aberrant and non-aberrant patterns per country. For this example, the R-code is given below

set1 <- which(XFT$Country=="USA")

set2 <- which(XFT$Country=="India")

set3 <- which(XFT$Country=="Philippines")

plot(out1$PFl,out1$lZPT,xlab="Person-fit Statistic RA",

ylab="Person-fit Statistic RT",col="black",cex=.5,bty="l",xlim=c(-3,3),

ylim=c(0,500),cex.main=.8,cex.axis=.7,cex.lab=.8,pch=15)

points(out1$PFl[set1],out1$lZPT[set1],col="blue",pch=10,cex=.5)

points(out1$PFl[set2],out1$lZPT[set2],col="red",pch=13,cex=.5)

points(out1$PFl[set3],out1$lZPT[set3],col="green",pch=16,cex=.5)

abline(h = qchisq(.95, df= 170),lty = 2,col="red")

abline(v = qnorm(.95),lty = 2,col="red")

legend(-3,500,c("India","US","Philippines","Other"), col=c("red","blue",

"green","black"),pch = c(13,10,16,15), bg = "gray95",cex=.7)

For both person-fit statistics, the threshold value of the significant area is marked with a dotted red line. With respect to aberrant RT patterns, a serious number of test takers are marked as aberrant, since their value is above the threshold of 201.4. Test takers from India have relatively more aberrant persons with respect to RT patterns, and the US test takers relatively less, although many of the aberrant test takers are US citizens. A few test takers are marked as aberrant with respect to their RA pattern, since their statistic value is above 1.645. Only five persons have been marked as aberrant in terms of the RA and RT patterns. None of the four quadrants seems to be dominated by one specific country.

The amsterdam chess data

The Amsterdam Chess data was collected by Van der Maas & Wagenmakers (2005) and is used to demonstrate a differential working speed analysis with the package LNIRT. The data study of Fox & Marianti (2016) is followed who also used this data set to illustrate their differential working speed model. In the Amsterdam Chess study, RA and RT patterns are collected of 259 test takers who responded to 40 chess tasks. The chess items have three subdimensions; tactical skill (20 items), positional skill (10 items), and end-game skill (10 items). Each item represents a chess board situation, and the problem-solving task was to select the best possible move. Both RAs and RTs were recorded, where RA was coded as 1 (correct) and 0 (incorrect). In this analysis, the differential working speed model in Eq. (23) is fitted to explore heterogeneity in speed trajectories and to examine the between-person relationship between ability and random speed components.

The LNIRTQ function is used to run the joint model with a differential working speed model with a random trend and quadratic time component. The RA and RT data were extracted from the Amsterdam Chess data to be used in the LNIRTQ function. The RA data consisted of 259 test takers (in rows) and 40 chess items (in columns).

 R> data(AmsterdamChess)

 R> N <- nrow(AmsterdamChess)

 R> Y <- as.matrix(AmsterdamChess[c(which(colnames(AmsterdamChess)=="Y1"):

         which(colnames(AmsterdamChess)=="Y40"))])

 R> K <- ncol(Y)

 R> Y[Y==9] <- NA

 R> RT <- as.matrix(AmsterdamChess[c(which(colnames(AmsterdamChess)=="RT1"):

         which(colnames(AmsterdamChess)=="RT40"))])

 R> RT[RT==10000.000] <- NA

 R> RT <- log(RT) #logarithm of RTs

The missing RA and RT values coded as 9 and 10,000 are replaced by NAs. There are three records with all missing values (rows 147,201, and 209). In LNIRTQ, imputations are generated under the model for the missing data even if complete records are missing. A time scale is defined by starting at zero and taking steps of 1/K to end in (K−1)/K. The LNIRTQ model is identified by restricting the difficulty and discrimination parameters to zero and one, respectively. Furthermore, the mean of each random speed parameter (intercept, trend, quadratic) is restricted to zero, and the product of time discriminations is restricted to one. Finally, the covariance of the speed components is restricted to zero, since the covariance among the speed components is modeled by the time discriminations. The LNIRTQ function is ran for 10,000 MCMC iterations, with a default burnin-period of 10%.

 R> X <- 1:K

 R> X <- (X - 1)/K

 R> outchess <- LNIRTQ(Y=Y,RT=RT,X=X,XG=10000)

 R> summary(outchess)

MCMC chains can be checked in the same way for convergence and run length by computing the Geweke and Heidelberger statistics, the effective sample sizes, the MCMC standard errors. The output of LNIRTQ contains MCMC chains of item parameters and hyper prior parameters for items and persons. The examined chains did not have any non-convergence problems. The effective sample size of the chains ranged from 466 (time discrimination item 1) to 9,000 (time intensity item 1). The naive posterior standard deviation were often close to the time-series standard error. The 10,000 MCMC iterations were sufficiently long and the chains showed convergence after 1,000 iterations.

The summary function of LNIRTQ provides a similar output as the LNIRT function. The posterior mean estimates and (naive) standard deviations of the item parameters are given, in addition to the covariance matrix of the item and person parameters. For reasons of brevity, the correlation matrix of the items and the covariance matrix of the persons are discussed.

The correlation matrix of the items is again given in the order of $a$, $b$, $\varphi$, and $\lambda$. The correlation between discrimination and difficulty is around 0.37, and between time discrimination and intensity around −0.53. The difficult items tend to discriminate better in ability than the less difficult items. The high time-intensive items do not discriminate well between speed levels. RTs were hardly affected by increasing working speed for the high time-intensive items. For low time-intensive items this effect on the RTs is much higher. There is a strong correlation between item difficulty and time intensity of 0.79, which states that the difficult items also took more time to be completed.

  --- Covariance matrix Items ---

   Item Matrix Correlation

[,1] [,2] [,3] [,4]

 [1,] 1.000  0.371  0.050  0.406

 [2,] 0.371  1.000-0.602  0.785

 [3,] 0.050-0.605  1.000-0.534

 [4,] 0.406  0.785-0.534  1.000

    --- Covariance matrix Person ---

Estimated Value

  Covariance Matrix

Theta Intercept Slope1 Slope2

 Theta0.4230.114-0.004-0.014

 Intercept0.1140.0600.0000.000

 Slope1-0.0040.0000.1140.000

 Slope2-0.0140.0000.0000.063

The covariance matrix of the random person parameters are given under the label Theta (ability), Intercept (speed intercept), Slope1 (speed trend), and Slope2 (speed quadratic component). The variance in working speed intercepts (i.e., the variance in starting speed of persons) is small (0.06). The speed trajectories of the persons show a variance of 0.11 in trends. The variance in deceleration/acceleration in working speed is around 0.06. The positive covariance between ability and the speed intercept shows that high-ability persons are more likely to start with a higher speed. The negative covariance between ability and the speed trend shows that the speed trajectories of high-ability persons has a negative trend (decrease in speed), in comparison to low-ability persons. The covariance between ability and the quadratic speed component is around −0.014, which means that high-ability persons are more likely to show an acceleration in the negative trend in speed.

For the records with missing values, the imputed data leads to population-average estimates of the random person effects. The ability estimate is around the population average of −0.162, and the speed components are around zero. For instance, for the record of 147 the estimated random person effects are (ability, intercept, trend, quadratic):

 R> outchess$Mtheta[147,]

 [1] -0.228 -0.018 -0.004 -0.004