# Longitudinal item response theory models in R

I'm trying to fit longitudinal item response theory (IRT) models in R. I have a test that was administered at multiple measurement occasions. I'd like to examine individuals' growth curves of factor scores (i.e., ability levels) from graded response models (GRMs). I have used the ltm package in R to fit cross-sectional GRM models in IRT, but it's not clear to me how (or whether it's even possible in ltm) to extend the models to handle repeated measures of the same items across time. How can I fit growth curves to longitudinal GRM factor scores to see changes in means/variances of ability levels over time? If this isn't possible with the ltm package, what packages/functions permit this? Specific code examples would be especially appreciated.

Here are empirical examples similar to what I'm looking to do (except the examples use a Rasch model for binary items where I'm using GRMs for polytomous data):

McArdle, J. J., Grimm, K. J., Hamagami, F., Bowles, R. P., & Meredith, W. (2009). Modeling life-span growth curves of cognition using longitudinal data with multiple samples and changing scales of measurement. Psychological Methods, 14, 126-149. doi: 10.1037/a0015857

McArdle, J. J., & Grimm, K. J. (2011). An empirical example of change analysis by linking longitudinal item response data from multiple tests. In A. A. Davier (Ed.), Statistical models for test equating, scaling, and linking (pp. 71-88): Springer: New York.

For example, I'd like to be able to estimate and plot individuals' growth curves to examine mean-level change in ability levels over time (from McArdle & Grimm, 2011):

And, I'd like to be able to estimate an average or prototypical growth curve for the sample (from McArdle & Grimm, 2011):

Here's a simulated data set with 20 polytomous items (1-3 response scale) at 3 different time points:

library(mirt)
library(mvtnorm)

set.seed(1)

numberItems <- 20
numberItemLevels <- 2
sampleSize <- 1000

a <- matrix(rlnorm(numberItems, .2, .2))
d <- matrix(rnorm(numberItems*numberItemLevels), numberItems)
d <- t(apply(d, 1, sort, decreasing=TRUE))

Theta <- mvtnorm::rmvnorm(n=sampleSize, 0, matrix(1))

t1 <- simdata(a, d, N=sampleSize, itemtype="graded", Theta=Theta)
t2 <- simdata(a, d, N=sampleSize, itemtype="graded", Theta=Theta+.5)
t3 <- simdata(a, d, N=sampleSize, itemtype="graded", Theta=Theta+1)

dat <- data.frame(t1, t2, t3)

• Could you perhaps give an example of your data, or something with the same structure? – Placidia Jun 24 '14 at 17:22
• Good suggestion, just added a simulated dataset to my question. – itpetersen Jun 24 '14 at 17:56
• I don't really understand the inspiration for this simulated data. If you are wondering how persons change over repeated test administrations, why not generate a Theta matrix for simdata() that actually changes over time and use that instead of just random values for each test? This just looks like 3 completely unrelated tests were administered. – philchalmers Jun 25 '14 at 23:47
• I was just giving the general format of the data (one should still be able to write the code to examine whether there is mean level change in factor scores across time regardless of whether there is actually mean level change). In any case, if it helps, I modified the simulated data set to show a mean level increase over time. I don't know how to use Theta to simulate the covariance structure, so I just calculated theta to have a mean level increase from T1 to T3. – itpetersen Jun 26 '14 at 13:26
• Are you comfortable with the following ideas: All items are treated on their own. You want to estimate discrete changes of the people on the latent scale over time for each item plus a general trend. Perhaps you even have groups? If you can also assert that all changes on items/persons are independent of each other and there is no unobserved heterogeneity, there may be a way with linear logistic models with relaxed assumptions, see here: statmath.wu-wien.ac.at/~hatz/psychometrics/11s/… – Momo Jun 26 '14 at 13:36

As a precursor, the IRT approach to this problem is very demanding computationally due to the higher dimensionality. It may be worthwhile to look into structural equation modeling (SEM) alternatives using the WLSMV estimator for ordinal data since I imagine less issues will exist. Plus, including external covariates is much easier within that framework. Both approaches I describe here are also possible in SEM.

There are two ways that I know of which you can estimate unidimensional longitudinal IRT models that are not Rasch in nature. The the first approach requires a unique latent factor for each time block and a specific residual variation term for each item. A different approach, similar to what one would find in the SEM literature, is via a latent growth curve model whereby only a fixed number of factors are estimated (three if the relationship over time is believed to be linear). Fixed loadings are used in this approach, so computationally it may be much more stable due to the reduced number of estimated parameters, so I would tend to prefer the growth curve model for both the smaller dimensionality and fewer estimated parameters.

The idea for both approaches is to set up latent time factors indicating how person level $\theta$ values change over each test administration, and constrain the influence of their loadings across time as well so that their hyper parameters can be estimated (i.e., the latent mean and covariances). Item constraints must also be imposed across time to remain invariable so that the person differences are only captured in the hyper parameters. Since this approach can require a huge number of integration dimensions, so you'll need to use something like the dimensional reduction algorithm which is available in mirt under the bfactor() function.

Instead of going through a worked example here, which would take a lot of code, I'll simply point to a worked versions of these analyses. A word of warning though, these are very computationally demanding and may take more than an hour to converge on your computer since you have 4 dimensions of integration in the first case and 3 dimensions in the second. Or, if you don't have much RAM you could experience issues when increase the number of quadpts.

Data simulation script: https://github.com/philchalmers/mirt/blob/gh-pages/data-scripts/Longitudinal-IRT.R

Analysis output: http://philchalmers.github.io/mirt/html/Longitudinal-IRT.html

In the first example, if you save the factor scores by using fscores() you'll obtain estimates for each time point regarding how individual $\theta$ values are changing. In the second example, using the linear growth curve approach, the first column of the factor scores will represent the initial $\theta$ estimates while the second column will indicate the slope/change occurring on average over time. In the example, I set up a constant mean change of .5, so the values in fscores() should all be around 0.5 for each individual. Both analyses give roughly the same conclusions but are somewhat different approaches to the problem. However, if you are familiar with longitudinal models in SEM then these should be fairly natural to interpret.

• This is really helpful, Phil. Any suggestions where to look for SEM alternatives? I'm fairly partial to the latent growth curve model framework. – itpetersen Jun 26 '14 at 20:48
• I'm sure the Mplus manual has a few from what I remember. If not, the categorical variables only really add the problem of what to do with intercepts to make them invariant across time so as to not affect the hyper parameters. So if you know LGM fairly well, that's probably all you really need to remember when working in the SEM framework. – philchalmers Jun 26 '14 at 21:10
• It might also be worthwhile to ask the lavaan Google group this question since that package supports the WLSMV estimator, perhaps pointing to the example and background information asked here. They are a nice bunch of folk to talk to, and I'm sure some of them have worked with categorical latent growth curve models before. groups.google.com/forum/#!forum/lavaan – philchalmers Jun 27 '14 at 15:53
• Thanks again, Phil. For people browsing this thread in the future, I found a relevant example (6.15) in the Mplus User's Guide v7: statmodel.com/download/usersguide/… – itpetersen Jun 29 '14 at 22:28
• thank you for providing the mirt example. Would you recommend any literature regarding the advanced use of the mirt package. I was trying to follow your code, but I had difficulties following the syntax regarding the constraints, etc. I have read the mirt package manual, but I was looking for a more verbose source. – user151310 Nov 3 '14 at 15:46

In the IRT literature for complicated IRT models (multiple groups, longitudinal/repeated measures, multidimensional) the recommended framework is Bayesian, because of relative easiness of estimation. I have had good experience using R package "rstan", which implements a flavor of Hamiltonian Monte Carlo. I had a data set with almost 3000 subjects measured at 6 time points, had 5 dimensions and 3 groups, and rstan worked really well. I was fitting a GRM model to ordinal test with 30 items. Check out the stan user group for example code.