This post provides an introduction to fitting item response theory (IRT) models in R. From my experience, most scholars in the social sciences have heard about IRT as an alternative to classical test theory (and its methods such as EFA or CFA), but have never really worked with it. I believe that this is unfortunate as it offers a lot of advantages and insights into the validity and reliability of tests and items. In this tutorial, I will use the R packages mirt (multidimensional item response theory, version 1.36.1; Chalmers, 2022, CRAN page) and ggmirt (an extension I wrote to provide publication-ready plotting functions: https://github.com/masurp/ggmirt).

Introduction

IRT refers to a set of mathematical models which aim to explain the relationship between a latent ability, trait, or proficiency (denoted \(\theta\)) and its observable manifestations (e.g., multiple-choice questions, true-false items, items…). In contrast to classical test theory (CTT), IRT focuses on the pattern of responses and considers responses in probabilistic terms, rather than focusing on composite variables and linear regression theory. IRT thereby accounts for…

item discrimination: ability of an item to differentiate between respondents with different levels of proficiency.
item difficulty: the likelihood of a correct response, expressed as the proficiency level at which 50% of the participant sample is estimated to answer an item correctly.
Depending on the model some other parameters such as e.g., guessing probability…

Assessing item difficulties is useful in matching the test and trait levels of a target population and in ensuring that the entire range of the proficiency is covered. Therefore, there are several advantages to IRT models (compared to CTT models):

The stronger focus on item difficulty generally leads to scales that perform better at differentiating at the extremes of proficiency.
Item banks that are developed using IRT often provide a good basis for creating smaller scales or parallel tests that still exhibit a high validity and reliability.
It allows for computer adaptive testing (CAT). This approach can greatly approach the efficiency of testing. In CAT, the test items are selected to match each individual’s proficiency, so that the s/he will not be bored by easy items or frustrated by overly difficult items.
Specific IRT models (e.g., the Rasch model, see further below) has specific mathematical properties that are highly desirable, such as a) the number-correct score is a sufficient estimation of \(\theta\), b) specific objectivity, which means that item and person parameters will be similar even if only a subset of the item pool is used or a different population is studied.

The theoretical background of IRT models is to comprehensive to be covered here. Rather, this tutorial aims to introduce the main models and discuss their properties. For further study, we refer to the excellent book “Item Response Theory” by Christine DeMars. In this short introductory tutorial, we focus on three of the most used IRT models: the 3PL model, the 2PL model, and finally the 1PL or Rasch model. These models are named by the number of item parameter used in the function that models the relationship between \(\theta\) and the item response (0/1). Each model has unique properties but all of them are suited to estimate latent variables from binary items (e.g., knowledge tests), which we will deal with in this tutorial. There are also more complex IRT models, such as e.g., graded response models, which can be used for non-binary items (e.g., likert-type scales). Yet, these will be discussed in a more advanced tutorial.

Preparation and Data

Although we are going to focus on the the package mirt in this tutorial, there are actually several packages that can be used to estimated IRT models. The most common ones are ltm, eRm, or TAM. This article by Choi & Asilkalkan (2019) also provides an overview of all available packages and their specific advantages: https://doi-org.vu-nl.idm.oclc.org/10.1080/15366367.2019.1586404. For this tutorial, we are going to load the tidyverse (for data wrangling and visualization) and mirt (for the main IRT analyses). And as noted before, we further are going to load my package ggmirt, which represents an extension to mirt and provides functions to plot more publication-ready figures and helps to assess item, person, and model fit.

# Data wrangling
library(tidyverse)

# Very comprehensive package for IRT analyses
library(mirt)

# Extension for 'mirt' 
# devtools::install_github("masurp/ggmirt")
library(ggmirt)

For ggmirt package includes a convenient function to simulate data for IRT analyses. Let’s quickly create a data set with 500 observations and 10 items that should fit a 3PL, 2PL and perhaps even a 1PL model.

set.seed(42)
d <- sim_irt(500, 10, discrimination = .25, seed = 42)
head(d)

V2	V3	V4	V6	V7	V8	V9	V10
1	0	1	1	1	1	0	1
1	1	0	0	0	0	0	1
1	0	1	1	0	1	0	0
0	1	0	1	0	0	0	1
1	1	1	0	1	1	1	0
0	1	0	0	0	1	0	1

As you can see, each “participants” has answered 10 items that are binary. Imagine we would have administered a test (e.g., LSAT, PISA, knowledge test, exam,…) to 500 people. The score 1 means this person has answered a particular question/item correctly. The score 0 means, this person has answered this item falsely.

3PL model

The 3PL model takes item discrimination (first parameter: a), item difficulty (second parameter: b), and guessing probability (third parameter: c) into account. As such, the 2PL and 1PL model (discussed below, are special cases, or constrained versions of the 3PL model). Take a look at Fig. 1 below. It shows a typical item characteristic curve (ICC, but not to be mistaken for the intra-class correlation). The x-axis shows the latent ability (\(\theta\)) ranging from -4 to 4, with 0 being the average ability in the studied population. The y-axis shows the probability of solving the item. The curve thus represents the probability of answering this item given a certain level on the latent ability.

In this example, this 3PL model provides a smooth curve for this item because it is based on item discrimination (steepness of the slope), difficulty (point of inflexion at which the probability of answering the item correctly is 50%) and the guessing probability (slightly raised lower asymptote). The mathematical form of the 3PL model is:

\[P(\theta_j) = c_i + (1-c_i)\frac{e^{1.7a_i(\theta_j-b_i)}}{1+e^{1.71a_i(\theta_j-b_i)}}\]

Here, the probability of a correct response \(P(\theta)\) given a person’s ability \(\theta\) is expressed as a function of all three item parameters and a person’s ability, where \(i\) indicates the item and \(j\) indicates an individual person. In the figure, the steeper the slope, the better the item is at differentiating people clearly in close proximity to its difficulty. The lower the asymptote, the lower the likelihood of selecting the right response by chance. As the item difficulty (point of inflexion is slightly above the average: 0.71), the item can be solved by a majority of potential participants.

So how do we fit such a model?

Fitting the model

To fit any IRT model, we simply use the function mirt(). In a first step, we specify the model itself. In this case, we simply want to estimate a uni-dimensional model, so we use the following string command "Factor = item_1 - item_n". The syntax for specifying models is comparatively simple, but for more complex (e.g., multidimensional models), see the some examples using ?mirt. It makes sense to save this string in a separate object.

Next, we provide the function itself with

the data set, which only contains the relevant variables (i.e., potential id variables have to be excluded)
the previously specifed model (string)
the type of model we want to estiamte (in this case “3PL”)

I am adding verbose = FALSE to not print information about the iterations. But you can remove this as well.

unimodel <- 'F1 = 1-10'

fit3PL <- mirt(data = d, 
               model = unimodel,  # alternatively, we could also just specify model = 1 in this case
               itemtype = "3PL", 
               verbose = FALSE)
fit3PL

## 
## Call:
## mirt(data = d, model = unimodel, itemtype = "3PL", verbose = FALSE)
## 
## Full-information item factor analysis with 1 factor(s).
## Converged within 1e-04 tolerance after 24 EM iterations.
## mirt version: 1.33.2 
## M-step optimizer: BFGS 
## EM acceleration: Ramsay 
## Number of rectangular quadrature: 61
## Latent density type: Gaussian 
## 
## Log-likelihood = -2738.689
## Estimated parameters: 30 
## AIC = 5537.377; AICc = 5541.343
## BIC = 5663.816; SABIC = 5568.594
## G2 (993) = 497.87, p = 1
## RMSEA = 0, CFI = NaN, TLI = NaN

The created object is of class “SingleGroupClass” and contains all necessary information and data to assess the model. By running the object itself, we get some information about the type of estimation as well as some model fit indices (including AIC and BIC), which can be used to compare models with one another.

Understanding IRT parameter

To a certain degree, an IRT analysis is similar to a factor analysis in CTT. If we use the summary() function, we get the so-called factor solution including factor loadings (F1) and the communalities (h2), which are squared factor loadings and are interpreted as the variance accounted for in an item by the latent trait. Almost all of the items in this case have a substantive relationship (loadings > .50) with the latent trait.

# Factor solution
summary(fit3PL)

##        F1    h2
## V1  0.707 0.500
## V2  0.513 0.264
## V3  0.535 0.286
## V4  0.727 0.529
## V5  0.574 0.329
## V6  0.566 0.320
## V7  0.603 0.364
## V8  0.443 0.196
## V9  0.718 0.516
## V10 0.480 0.230
## 
## SS loadings:  3.534 
## Proportion Var:  0.353 
## 
## Factor correlations: 
## 
##    F1
## F1  1

In IRT, however, we are usually more interested in the actual IRT parameters as discussed above (discrimination, difficulty and guessing probability). They can be extracted as follows:

params3PL <- coef(fit3PL, IRTpars = TRUE, simplify = TRUE)
round(params3PL$items, 2) # g = c = guessing parameter

	a	b	g	u
V1	1.70	1.17	0.00	1
V2	1.02	-0.58	0.00	1
V3	1.08	0.35	0.00	1
V4	1.80	0.81	0.09	1
V5	1.19	0.52	0.00	1
V6	1.17	-0.14	0.00	1
V7	1.29	1.87	0.00	1
V8	0.84	0.03	0.00	1
V9	1.76	2.11	0.00	1
V10	0.93	-0.05	0.01	1

The values of the slope (a-parameters) parameters ranged from 0.84 to 1.80. This parameter is a measure of how well an item differentiates individuals with different theta levels. Larger values, or steeper slopes, are better at differentiating people. A slope also can be interpreted as an indicator of the strength of a relationship between and item and latent trait, with higher slope values corresponding to stronger relationships.

The location or difficulty parameters (b-parameter) is also listed for each item. Location parameters are interpreted as the value of theta that corresponds to a .50 probability of responding correctly at or above that location on an item. The location parameters show that the items cover a wide range of the latent trait.

Model fit, person fit, and item fit evaluation

Similar to factor analytical approaches, we can assess how well the model fits the data. Rather than using a a \(\chi^2\) statisic, we use a specific index, M2, which is specifically designed to assess the fit of item response models.

M2(fit3PL)

	M2	df	p	RMSEA	RMSEA_5	RMSEA_95	SRMSR	TLI	CFI
stats	17.77133	25	0.8519424	0	0	0.0203603	0.033713	1.018438	1

As we can see, the M2 statistic is comparatively low and non-significant. So there are no concerning differences between the model and the data. This is further supported by a very low RMSEA and a CFA and TLI of 1.

In IRT, however, we are usually more interested in item– and person-fit indices. IRT allows us to assess how well each items fits the model and whether the indiviual response patterns align with the model.

Let us start with item fit: Similar to many other areas, different indices have been proposed to assess item fit. We can use the function itemfit() to get a variety of them. By default, we receive the S_X2 by Orlando and Thissen (2000) and the corresponding dfs, RMSEA and p-values. This test should be non-significant to indicate good item fit. As we can see here, only item V9 shows a lower fit with the model. Proponents of the “Rasch Model” (see further below) often rather report infit and outfit statistics. We can get those by adding the argument fit_stats = "infit"). We get both mean-squared and standardized versions of these measures (Linacre provides some guidelines to interpret these: https://www.rasch.org/rmt/rmt162f.htm). Roughly speaking the non-standardized values should be between .5 and 1.5 to not be degrading. In the package ggmirt, we can also use the function ´itemfitPlot()` to inspect this visually.

itemfit(fit3PL)

item	S_X2	df.S_X2	RMSEA.S_X2	p.S_X2
V1	4.708662	4	0.0188425	0.3185173
V2	2.503355	5	0.0000000	0.7759897
V3	6.761645	5	0.0265720	0.2389791
V4	3.700650	5	0.0000000	0.5932672
V5	4.093206	5	0.0000000	0.5360762
V6	2.824089	5	0.0000000	0.7270838
V7	8.648072	5	0.0382381	0.1239519
V8	2.901227	5	0.0000000	0.7152104
V9	11.610311	4	0.0617477	0.0204970
V10	2.693987	5	0.0000000	0.7470379

itemfit(fit3PL, fit_stats = "infit") # typical for Rasch modeling

item	outfit	z.outfit	infit	z.infit
V1	0.6386014	-2.365002	0.8743268	-1.7362580
V2	0.8481622	-3.058789	0.8705611	-3.7125518
V3	0.8320215	-3.792355	0.8794736	-3.4444109
V4	0.8362322	-2.532776	0.8729352	-2.5664209
V5	0.8035407	-3.634037	0.8727943	-3.1862993
V6	0.8104688	-4.158982	0.8432826	-4.6042826
V7	0.7308953	-1.594022	0.9935245	-0.0331537
V8	0.8929857	-3.419761	0.9110881	-3.2959228
V9	0.5264292	-1.359471	1.0527710	0.3834486
V10	0.8713525	-3.689040	0.8945440	-3.6438640

itemfitPlot(fit3PL)

We again see that item V9 has a lower fit (outfit value close to .5), but according to Linacre’s guidelines, this should not be concerning.

Assessing person fit

We can technically produce the exact same measures for each person to assess how well each person’s response patterns aligns with this model. Think about it this way: If a person with a high theta (that is high latent ability) answers does not answer an easy item correctly, this person does not fit the model. Conversely, if a person with a low ability answer a very difficult question correctly, it likewise doesn’t fit the model. In practice, there will most likely be a few people who do not fit the model well. But as long as the number of non-fitting respondents is low, we are good. We mostly look again at infit and outfit statistics. If less than 5% of the respondents have higher or lower infit and outfit values than 1.96 and -1.96, we are good.

head(personfit(fit3PL))

outfit	z.outfit	infit	z.infit	Zh
1.0161479	0.1749541	1.0183924	0.1584039	-0.0732308
0.6331922	-0.1327510	0.8411215	-0.4976128	0.5615934
0.7456650	-0.1905736	0.8912849	-0.3611853	0.4390657
0.7209540	-0.0392887	0.9571507	-0.0637470	0.2601261
2.1363674	2.3695393	1.7713727	2.1832656	-2.6948633
0.7538809	0.0610415	1.0221023	0.1724740	0.0843561

personfit(fit3PL) %>%
  summarize(infit.outside = prop.table(table(z.infit > 1.96 | z.infit < -1.96)),
            outfit.outside = prop.table(table(z.outfit > 1.96 | z.outfit < -1.96))) # lower row = non-fitting people

infit.outside	outfit.outside
0.958	0.98
0.042	0.02

personfitPlot(fit3PL)

How to run IRT analyses in R