lsirm12pl: An R Package for the Latent Space Item Response Model
Dongyoung Go
Department of Statistics and Data ScienceGwanghee Kim
Department of Statistics and Data ScienceJina Park
Department of Statistics and Data ScienceJunyong Park
Samsung Electronics, Formerly at Yonsei UniversityMinjeong Jeon
School of Education and Information StudiesIck Hoon Jin
Department of Statistics and Data Science2025-10-21
Abstract
The item response model in latent space (LSIRM; Jeon et al. (2021)) uncovers unobserved interactions between respondents and items in the item response data by embedding both in a shared latent metric space. The R package lsirm12pl implements Bayesian estimation of the LSIRM and its extensions for various response types, base model specifications, and missing data handling. Furthermore, the lsirm12pl package provides methods to improve model utilization and interpretation, such as clustering item positions on an estimated interaction map. The package also offers convenient summary and plotting options to evaluate and process the estimated results. In this paper, we provide an overview of the LSIRM’s methodological foundation and describe several extensions included in the package. We then demonstrate the use of the package with real data examples contained within it.
Introduction
Item response theory (IRT) models are a widely-used statistical approach to analyze assessment data in various fields, e.g., medical, educational, psychological, health, and marketing research (An and Yung 2014; Zanon et al. 2016). IRT models are designed to establish a relationship between observed item response data and unobserved person characteristics, commonly referred to as latent traits, e.g., competencies, attitudes, or personality (Ayala 2009; Brzezińska 2018). IRT models can predict the probability of a correct (or positive) response as a function of the respondents’ latent traits and item features, such as item difficulty and discrimination. Additional technical details of IRT models are provided in the subsequent section.
Several R packages are available for estimating IRT models. The ltm package (Rizopoulos 2006) is available to analyze dichotomous and polytomous item response data, including a one-parameter logistic (1PL) model (or the Rasch model), a two-parameter logistic (2PL) model, a three-parameter model (Rasch 1960; Birnbaum 1968) and a graded response model (Samejima 1968). The eRm package (Mair and Hatzinger 2007) estimates various extensions of the Rasch model, such as the rating scale model (RSM) (Andrich 1978), partial credit model (PCM) (Masters 1982), linear logistic test model (LLTM) (Scheiblechner 1972), the linear rating scale model (LRSM) (G. H. Fischer and Parzer 1991), and the linear partial credit model (LPCM) (Glas and Verhelst 1989; Gerhard H. Fischer and Ponocny 1994). The mirt package (Chalmers 2012) can estimate a wide range of IRT models, including exploratory and confirmatory multidimensional item response models. The pcIRT package (Hohensinn 2018) provides functions for estimating IRT models for polytomous (nominal) and continuous data, including the multidimensional polytomous Rasch model (Andersen 1973) and the continuous rating scale model (Müller 1987), using conditional maximum likelihood (CML) estimation (Baker and Kim 2023).
Conventional IRT models are typically based on two assumptions: conditional independence and homogeneity. Conditional independence assumes that item responses are independent of each other conditional on respondents’ latent traits. The homogeneity assumption is, e.g., that respondents with the same ability have the same probability of answering a question correctly and that respondents have the same probability of answering a question with the same item difficulty. However, these assumptions are often violated in practice due to unobserved interactions between respondents and items, e.g., when particular items are more similar to each other (e.g., testlets) or when particular respondents show different probabilities of giving correct responses to certain items compared to other respondents with similar ability (e.g., differential item functioning). Violations of these assumptions can lead to biased parameter estimates and inferences (Chen and Wang 2007; Braeken 2010; Myszkowski and Storme 2024). While some existing methods can address known violations prior to data analysis, there is currently no approach that enables the detection or management of unknown sources of violations in item response analysis, to the best of our knowledge.
Jeon et al. (2021) proposed a latent space item response model (LSIRM) that addresses such limitations of conventional IRT models. The LSIRM aims to estimate inherent interactions between respondents and items, alongside the latent traits of both. A key feature is its visual representation of these interactions in a low-dimensional latent space, called an interaction map in the form of distances between them. Interaction maps offer a clear and intuitive interpretation of complex respondent-item relationships, allowing users to identify patterns and clusters based on spatial proximity on the map. Additional details of the LSIRM are provided in a later section.
This paper presents the lsirm12pl package in R that offers Bayesian estimation of the LSIRM and its extensions. Jeon et al. (2021) focused on the response data for binary items and the Rasch model as the base model, and currently, no package is available to estimate the LSIRMs. To broaden the applicability of latent space item response modeling, lsirm12pl enables: (1) modeling continuous item responses; (2) missing data handling under different missing mechanisms assumptions (Rubin 1976); and (3) an extended base model specification using the 2PL model both for binary and continuous item response data. The package also offers options to cluster latent positions of items in the estimated interaction map using spectral clustering and the Neyman-Scott process model. The lsirm12pl supplies convenient summary and plotting options for interaction maps, model assessment, diagnosis, and result process and interpretation, aiming to improve the utilization of the LSIRM in practice.
The subsequent sections of the paper are structured as follows. To begin, we provide a concise overview of the 1PL and 2PL IRT models. We then delve into the LSIRM for binary response data and demonstrate how to fit the model with the package functions using a real dataset. Next, we extend the LSIRM to accommodate continuous data and provide guidance on fitting this extended model using the lsirm12pl package. In the end, we conclude the paper with final remarks and a discussion of future developments.
Item response theory models
IRT models are essential for analyzing item response data, which comprises respondents’ answers to items on tests, surveys, or questionnaires. This section briefly discusses two widely utilized IRT models for dichotomous item response data.
1PL IRT Model
The 1PL model, also known as the Rasch model (Rasch 1961), is a classic IRT model for analyzing dichotomous item response data. Suppose \(\boldsymbol{Y} = \left\{ y_{k,i} \right\} \in \left\{ 0, 1 \right\}^{N \times P}\) is the \(N \times P\) binary item response matrix under analysis, where \(y_{k,i}=1\) indicates a correct (or positive) response of the respondent \(k\) to the item \(i\). In the Rasch model, the probability of the correct response for item \(i\) given by respondent \(k\) is given as follows: \[\text{logit}(\mathbb{P}(y_{k,i}=1|\theta_{k},\beta_{i}))=\theta_{k} +\beta_{i}, \; \; \theta_{k} \sim N(0, \sigma^2), \label{rasch_formula} (\#eq:rasch-formula)\] where \(\theta_k \in \mathbb{R}\) is the respondent intercept parameter of respondent \(k\), \(\beta_i \in \mathbb{R}\) is the item intercept parameter of item \(i\). The person intercept parameter represents the latent trait of interest, such as cognitive ability. The item intercept parameter represents the item easiness (or minus difficulty). As the respondent’s ability increases, her/his likelihood of giving correct responses increases. On the other hand, the likelihood of correct responses decreases as the difficulty of the item increases. Note that based on the model, the probability of a respondent’s giving a correct response to an item is a function of the two parameters – the person’s ability and the item’s difficulty. This means that respondents with the same level of ability are assumed to have the same probability of giving a correct response to an item. Similarly, items with the same level of difficulty are assumed to have the same probability of being correctly responded to. In other words, interactions between persons and items are not allowed in this model. However, in reality, respondents with the same level of ability and items with the same level of difficulty may show a different success probability due to unobserved characteristics they may have, violating the assumptions.
2PL IRT Model
The 2PL model (Birnbaum 1968) extends the 1PL model by incorporating a discrimination parameter for each item. This parameter reflects the ability of the item to differentiate among respondents with similar abilities. The item discrimination parameters are also considered as item slopes, indicating how quickly the probability of correct responses increases as the respondent’s abilities increase (An and Yung 2014). A larger discrimination parameter indicates that the item is better at distinguishing between respondents with similar ability levels. In the 2PL model, the probability of person \(k\) giving a correct response to item \(i\) is given as follows: \[\text{logit}(\mathbb{P}(y_{k,i}=1|\theta_{k}, \alpha_{i}, \beta_{i}))=\alpha_{i}\theta_{k} + \beta_{i}, \; \; \theta_{k} \sim N(0, \sigma^2), \label{birnbaum2pl_formula} (\#eq:birnbaum2pl-formula)\] where \(\alpha_i \in \mathbb{R}\) is the discrimination parameter for the item \(i\). For model identifiability, one of the item slope parameters is fixed at 1, e.g., \(\alpha_{1} =1\). With the term \(\alpha_{i}\theta_{k}\), the 2PL model allows for an interaction between respondents and items to some degree. However, it is not likely to capture all interactions between respondents and items that might not depend on respondents’ ability (Jeon et al. 2021).
Standard lsirm: 1pl lsirm for dichotomous data
The LSIRM (Jeon et al. 2021) has been proposed as an extension of the conventional IRT models. The key idea of the LSIRM is to place respondents and items in a low-dimensional, shared metric space, called an interaction map, so that their unobserved interactions can be captured in the form of distances between them. Below we provide a brief overview of the LSIRM in its original form, presented for dichotomized data.
Statistical framework
To capture unobserved interactions between respondents and items, the original LSIRM (Jeon et al. 2021) embeds items and respondents in an interaction map, i.e., a \(D\)-dimensional latent Euclidean space. Note that the interaction map is used as a tool to represent pairwise interactions between respondents and items; thus, the dimensions of an interaction map do not represent any specific quantity or have a substantive meaning. The original LSIRM is built on the 1PL IRT model; thus, we refer to the original LSIRM as the 1PL LSIRM.
The 1PL LSIRM assumes that the probability of giving a correct answer
to item \(i\) by the respondent \(k\) is determined by a linear combination
of the main effect of the respondent \(k\), the main effect of the item \(i\), and the pairwise distance between the
latent position of item \(i\) and the
latent position of respondent \(k\).
Then, the model is given by: \[\text{logit}(\mathbb{P}(y_{k,i}=1|\theta_{k},\beta_{i},\gamma,\boldsymbol{z_{k},w_{i}}))=\theta_{k}+\beta_{i}-\gamma
d (\boldsymbol{z_{k}},\boldsymbol{w_{i}}),
\label{1pl_formula} (\#eq:1pl-formula)\] where \(\theta_k \in \mathbb{R}\) and \(\beta_i \in \mathbb{R}\) represent the
respondent’s latent trait and the item’s easiness, respectively, same as
in the conventional 1PL model. These two terms can also be seen as the
main effects of respondents and items. The third term, \(-\gamma d
(\boldsymbol{z_{k}},\boldsymbol{w_{i}})\), captures the
interactions between respondents and items, where \(\boldsymbol{z_{k}}\in\mathbb{R}^{D}\) and
\(\boldsymbol{w_{i}}\in\mathbb{R}^{D}\)
are the latent position of the respondent \(k\), and the latent position of the item
\(i\), respectively, where \(\gamma \geq 0\) is the weight of the
distance term \(d(z_{k},w_{i})\). For a
distance function \(d:\) \(\mathbb{R}^{D} \times \mathbb{R}^{D} \mapsto
[0,\infty)\), we use a Euclidean norm (that is, \(d(z_{k},w_{i})=||z_k - w_i||\)) as a
distance function in lsirm1pl to maintain simplicity and
enhance the interpretability of the model (Hoff,
Raftery, and Handcock 2002). The weight of the distance term
\(\gamma \geq 0\) indicates the amount
of interactions between respondents and items in the data, such that
large \(\gamma\) implies stronger
evidence for respondent by item interactions in the data. In contrast,
near zero \(\gamma\) implies little
evidence of respondent-item interactions in the data, thus suggesting
that one may go with the conventional IRT model without the interaction
term.
The likelihood function of the observed data with the 1PL LSIRM is \[\mathbb{L}\left(\boldsymbol{Y=y | \theta, \beta}, \gamma, \boldsymbol{Z, W} \right) = \prod_{K=1}^{N} \prod_{i=1}^{P} \mathbb{P}(Y_{k,i} = y_{k,i} | \theta_{k}, \beta_{i}, \gamma, \boldsymbol{z_{k}, w_{i}}), \label{likelihood} (\#eq:likelihood)\] where \(\boldsymbol{\theta} = \left(\theta_{1}, \ldots, \theta_{N}\right)\), \(\boldsymbol{\beta}=\left(\beta_{1}, \ldots, \beta_{P}\right)\), \(\boldsymbol{Z} = \left(\boldsymbol{z_{1}}, \ldots, \boldsymbol{z_{N}}\right)\) and \(\boldsymbol{W}=\left(\boldsymbol{w_{1}}, \ldots, \boldsymbol{w_{P}}\right)\). Here, item responses are assumed to be independent conditional on the positions of respondents and items in an interaction map, as well as the main effects of respondent and item. This means that the traditional conditional independence assumption is alleviated with the model by accounting for respondent-item interactions. The detailed parameter estimations are described in the following Section.
Parameter Estimation
The R package lsirm12pl applies a fully Bayesian approach using the Markov chain Monte Carlo (MCMC) for estimation of the LSIRM. Following is the posterior distribution of the 1PL LSIRM. \[\begin{aligned} \pi\left(\boldsymbol{\theta, \beta}, \gamma, \boldsymbol{Z, W} | \boldsymbol{Y=y} \right) & \propto \prod_{k=1}^{N} \prod_{i=1}^{P} \mathbb{P}(Y_{k,i} = y_{k,i} | \theta_{k}, \beta_{i}, \gamma, \boldsymbol{z_{k}, w_{i}}) \\ &\times \prod_{k=1}^N \pi(\theta_k) \times \prod_{i=1}^P \pi(\beta_i) \times \pi(\gamma) \times \prod_{k=1}^N \pi(\boldsymbol{z_{k}}) \prod_{i=1}^P \pi(\boldsymbol{w_{i}}) \label{posterior} \end{aligned} (\#eq:posterior)\] The priors for the model parameters are specified as follows: \[\label{1pl-prior} \begin{split} \theta_{k} \mid \sigma^{2} & \sim N(0,\sigma^{2}),\ \sigma^{2}>0\\ \beta_{i} \mid \tau_{\beta}^{2} & \sim N(0,\tau_{\beta}^{2}),\ \tau_{\beta}^{2}>0\\ \log \gamma \mid \mu_{\gamma},\tau_{\gamma}^2 & \sim N(\mu_{\gamma},\tau_{\gamma}^2),\ \mu_{\gamma}\in\mathbb{R},\ \tau_{\gamma}^2>0\\ \sigma^{2} \mid a_{\sigma},b_{\sigma} & \sim\text{Inv-Gamma}(a_{\sigma},b_{\sigma}),a_{\sigma}>0,b_{\sigma}>0\\ \boldsymbol{z_{k}} & \sim\text{MVN}_{D}(\boldsymbol{0,I_{D}})\\ \boldsymbol{w_{i}} & \sim\text{MVN}_{D}(\boldsymbol{0,I_{D}}). \end{split} (\#eq:1pl-prior)\] where \(\boldsymbol{0}\) is a \(D\)-vector of zeros and \(\boldsymbol{I_{D}}\) is \(D \times D\) identity matrix. The argument names and default values for the prior specifications in the lsirm12pl are described in our Github site 1.
To generate posterior samples for \(\boldsymbol{\theta,\ \beta,\ \gamma,\ Z}\) and \(\boldsymbol{W}\), we use the Metropolis-Hastings-within-Gibbs sampler (Chib and Greenberg 1995). The conditional posterior distribution for each parameter is given in Equations in the Github site. We list the arguments, the default values for the jumping rules, and the standard deviations of the Gaussian proposal distributions in the Github site.
The log-odds of the probability of giving correct responses depend on the latent positions through distances, as discussed in Jeon et al. (2021). Because distances are invariant to translations, reflections, and rotations of the positions of respondents and items, the likelihood function is invariant under these transformations. To resolve the identifiability issue of latent positions, the lsirm12pl package applies Procrustes transformation (Gower 1975) as a post-processing of posterior samples, which is a standard procedure in the latent space modeling literature (Hoff, Raftery, and Handcock 2002; Sewell and Chen 2015; Jeon et al. 2021).
An Illustrated Example
Here, we demonstrate how to apply the 1PL LSIRM to real datasets using the lsirm12pl package. To this end, we use the Inductive Reasoning Developmental Test (TDRI) dataset (Golino 2016) from the package, which contains item responses from 1,803 Brazilians (52.5% female) of ages ranging from 5 to 85 years (M = 15.75; SD = 12.21). TDRI is a pencil-and-paper test consisting of 56 items that are designed to assess developmentally sequenced and hierarchically organized inductive reasoning.
The lsirm12pl
package provides several functions for fitting the LSIRM, which requires
setting hyperparameter values for prior distributions and/or tuning
parameters for MCMC chains. By default, the lsirm function
runs with the default settings unless otherwise specified by the user.
The default MCMC run setting includes 15,000 iterations, 2,500 burn-ins,
and 5 thinning. The base function for running the LSIRM is
lsirm(A\(\sim\)<term 1>(<term 2>, <term 3>, ...))
where A is an item response matrix to be analyzed,
<term1> is for the model option – either ‘lsirm1pl’
or ‘lsirm2pl’, while <term 2> and
<term 3> are other specific options for the chosen
model, which are detailed in the documentation of the lsirm12pl
package. The following is an example of how to fit the 1PL LSIRM to the
TDRI dataset (with no missing) with a default estimation setting with 4
MCMC chains using 2 multi-core processors. Additional details of other
default values can be found in the documentation of the lsirm12pl
package.
R > library("lsirm12pl")
R > data <- lsirm12pl::TDRI
R > data <- data[complete.cases(data),]
R > head(data)
i1 i2 i3 i4 ... i54 i55 i56
1 1 1 1 1 ... 0 0 0
2 1 1 1 1 ... 0 0 0
3 1 1 1 1 ... 0 0 0
4 1 1 1 1 ... 0 0 0
5 1 1 1 1 ... 0 0 0
R > lsirm_result <- lsirm(data ~ lsirm1pl(chains = 4, multicore = 2, seed = 2025))The estimation results of the model parameters \(\boldsymbol{\theta,\ \beta},\ \gamma,\ \boldsymbol {Z}\), and \(\boldsymbol{W}\) are stored in individual lists per chain, while all results across the chains are stored in a single list.
The summary() function generates a summary of the first
chain by default, but users can obtain summaries of other chains by
setting the chain.idx option. The function provides
posterior estimates of the “covariate coefficients” (model parameters
such as \(\beta\)), allowing the users
to select the "mean", "median", or
"mode" through the estimate option. It also
provides the highest posterior density interval (HPD) for the model
parameters with the CI option that allows the users to set
different significance levels. Additionally, the function supplies the
Bayesian information criterion (BIC) and the maximum log posterior
value. When the column names are available in the input data, the
summary() function uses these names in the summary of the
results.
R > summary(lsirm_result, chain.idx = 1, estimate = 'mean', CI = 0.95)
==========================
Summary of model
==========================
Call: lsirm.formula(formula = data ~ lsirm1pl(chains = 4, multicore = 2, seed = 2025))
Model: lpl LSIRM
Data type: binary
Variable Selection: FALSE
Missing: NA
MCMC sample of size 15000, after burnin of 2500 iteration
Covariate coefficients posterior means of chain 1 :
Estimate 2.5% 97.5%
i1 6.42354 5.82468 7.0642
...
i56 -1.01604 -2.58822 0.8133
---------------------------
Overall BIC (Smaller is better) : 43661.52
Maximum Log-posterior Iteration:
value iter
[1,] -11487 137The diagnostic() function checks the convergence of MCMC
for each parameter using various diagnostic tools, such as trace plots,
posterior density distributions, autocorrelation functions (ACF), and
Gelman-Rubin-Brooks plots. The diagnostic() function has
options: draw.item, and gelman.diag. The
draw.item option in the diagnostic() function
specifies the names and indexes of the parameters to diagnose. The
draw.item option is set to a list where a key represents
each parameter such as
“beta”, “theta”, “gamma”, “alpha”, “sigma”, and
“sigma_sd”, and the values indicate the indices of these
parameters. The indexes can be expressed as vectors. In the following
example code, the draw.item option is set as
list(“beta” = c(1)) to check the diagnostic of the first
index of the beta parameter, i.e. \(\beta_1\). With
gelman.diag = TRUE, the Gelman-Rubin convergence
diagnostic, known as potential scale reduction factors (PSRF), is
obtained.
R > diagnostic(lsirm_result,
draw.item = list("beta" = c(1)),
gelman.diag = TRUE)
Potential scale reduction factors:
Point est. Upper C.I.
beta [i1] 1.01 1.04
diagnostic() function on the results of the 1PL LSIRM
fitted to the TDRI dataset. The top left is the trace plot, the top
right is the posterior density plot, the bottom left is the
autocorrelation plot, and the bottom right is the Gelman-Rubin-Brooks
plot. Different colors represent different MCMC chains.Figure 1 displays the output of
diagnostic() for \(\beta_1\): trace plot (top left), posterior
density plot (top right), autocorrelation plot (bottom left), and
Gelman-Rubin-Brooks plot (bottom right). Different colors indicate
different MCMC chains. Trace plots visualize the mixing of the MCMC
chains. In a well-converged model, the trace plots should show that
chains fluctuate consistently around a constant value, which indicates
the parameter space is thoroughly explored. Density plots display the
distribution of the sampled parameter values. In a converged model, the
density plots should exhibit smooth, overlapping curves across different
chains, which demonstrates that the samples are drawn from the same
posterior distribution. Autocorrelation plots measure the correlation
between samples at different lags. For a converged model, the
autocorrelation should decrease rapidly as the lag increases, which
suggests the samples are independent and the parameter space has been
effectively explored. Gelman-Rubin-Brooks plots show a shrink factor,
known as the potential scale reduction, which compares the variance
within each chain to the variance between chains. A shrink factor value
close to 1 indicates the within-chain variance is similar to the
between-chain variance, providing evidence of good convergence. By
examining these diagnostics collectively, the convergence of the model
parameters from the LSIRM can be ensured.
The lsirm12pl
package supplies the gof() function that assesses the
goodness of fit of the LSIRM to the item response data being analyzed. A
main diagnostic tool is a boxplot that displays the average posterior
‘predicted’ response value for each item (item-wise means) with the
average ‘observed’ response value for each item indicated by red dots.
When the model fits well, the red dot is close to the midline of the
boxplot. The gof() function includes a
chain.idx option to choose a specific chain for assessing
the goodness of fit. By default, the first MCMC chain is selected
(i.e. chain.idx=1). When the diagnostic()
function shows good convergence for all parameters, it does not matter
which chain is chosen because the parameter estimates must be consistent
across all chains. For illustration, we use the first MCMC chain by
setting chain.idx=1.
For binary data, the gof() function additionally offers
a receiver operating characteristic (ROC) curve. The ROC curve is a
graphical representation that assesses the performance of a binary
classifier. The area under the curve (AUC) quantifies the overall
ability of the model to distinguish between classes, with values ranging
from 0.5 to 1. A higher AUC and a curve close to the top-left corner
indicate better performance.
R > gof(lsirm_result, chain.idx = 1)
gof() function. The box plots of
the average posterior predicted response values for items (item-wise
means) with the average observed response value for each item indicated
by red dots (left) and the receiver operating characteristic (ROC) curve
(right).Figure 2 presents the output of the gof()
function for the TDRI dataset, showing the results for the first MCMC
chain, including both the boxplot (left) and the ROC curve (right). In
this example, all of the red dots are closely located at the midline of
the boxplot and the area under the curve (AUC) of the ROC curve
approaches 1, suggesting that the fit of the LSIRM to the data is
satisfactory.
The main outcome of the LSIRM fitting is the estimates of the
respondent and item main effects and the interaction map. For a visual
summary of these outputs, the user can use the plot()
function with the option argument.
R > plot(lsirm_result, option = "beta", chain.idx = 1)
R > plot(lsirm_result, option = "theta", chain.idx = 1)
plot() function based on the 1PL LSIRM fitted to the TDRI
dataset. (a) Boxplots of the posterior samples for βi; outliers are
suppressed in the boxplots for the sake of simplicity. (b) Boxplots of
the point estimate for θk as a function
of the total sum scores of positive responses.
Figures 5 display the plot()
results for summarizing \(\beta_i\) and
\(\theta_k\) by setting the
option argument as the "beta" and
"theta", respectively. option="beta" generates
the boxplots of the posterior samples for \(\beta_i\), while
option="theta" generates the boxplots of the point
estimates for \(\theta_k\), plotted
against the total sum scores of positive responses (ranging from 0 to
\(P\), where \(P\) is the number of items). In Figure 3, as the item number increases (from
left to right on the x axis), the \(\beta_k\) estimates decrease, indicating
that earlier items are easier. In Figure 4, individuals
who correctly answer more items (i.e., higher sum scores) have higher
\(\theta_i\) estimates. It is sensible
that those who answer more items correctly have higher \(\theta_i\) values, as \(\theta_i\) represents their ability
levels.
The plot() also returns the interaction map with
option="interaction". The interaction map is created based
on the estimated latent position of the item and respondent, \(\boldsymbol{w}_i\) and \(\boldsymbol{z}_k\), respectively. Note that
the primary and unique advantage of the LSIRM lies in deriving intuitive
information from the interaction map, based on the distances between
respondents and items, between respondents, and between items. In the
interaction map, a shorter distance between the latent position of the
item \(i\) (\(\boldsymbol{w_i}\)) and the latent position
of the respondent \(k\) (\(\boldsymbol{z_k}\)) indicates a stronger
dependence (or interactions), which implies that the respondent \(k\) is more likely to respond correctly (or
positively) to the item \(i\), given
the person’s ability. In the lsirm12pl
package, the interaction map is set to a two-dimensional space, as in
(Jeon et al. 2021), for parsimony and
interpretability.
R > plot(lsirm_result, option = "interaction", chain.idx = 1)
Figure 6 (top left) displays the interaction map based on the 1PL LSIRM for the TDRI data. The latent positions of items are represented as red numbers and respondents as black dots on the map. Note that the two coordinates (dimensions) of the map do not represent specific quantities or substantive meaning, as the interaction map is simply a tool to represent respondent-by-item interactions (Jeon et al. 2021). In Figure 6, we observe that respondents are widely spread around the center of the interaction map, while the items are separated by some clusters. For example, items located in the far north, such as items 50, 52, 53, 54, and 55 (red numbers), are far from most respondents (black dots) on the map. This suggests that most respondents are less likely to respond correctly to those specific items, regardless of their ability levels. In other words, those items located in the south are difficult items, which is also confirmed in Figure 5(a).
R > plot(lsirm_result, option = "interaction", rotation = TRUE, chain.idx = 1)If desired, one can rotate the interaction map to improve the
interpretability of the coordinates of the map with
rotation=TRUE.
The lsirm12pl package offers the oblimin rotation (Jennrich 2002) using the GPArotation package in R (Bernaards and Jennrich 2005). Figure 7 (top left) shows the rotated version of the original map shown in Figure 6 (top right). In this example, the original and rotated maps appear pretty similar, although on the right there seems to be slight clockwise rotation compared to the one on the left; after rotation, the x- and y- coordinates may be interpreted based on the items that are closely placed to the respective coordinates.
Further, the lsirm12pl
package offers an option to cluster latent positions of items. Two types
of clustering methods are available: spectral clustering (Ng, Jordan, and Weiss 2001; Luxburg 2007) and
the Neyman-Scott process modeling approach (Thomas 1949; Neyman and Scott 1952; Yi et al.
2024), with cluster option as spectral
and neyman, respectively. Spectral clustering is a
technique that clusters points based on the eigenvalues of a similarity
matrix, using the spectral properties of the data to identify clusters.
The implementation of spectral clustering is based on the kernlab
package (Karatzoglou et al. 2004) in R
where the number of clusters is determined using the average silhouette
width (Batool and Hennig 2021). The
Neyman-Scott point process modeling approach is a method to cluster
points in time or space. Parent points (cluster center) are generated
using a Poisson process, with offspring points clustered around each
parent based on a defined probability distribution. The lsirm12pl
package implements this method directly; it applies the MCMC algorithm a
number of times independently to determine the distribution of the
cluster number and cluster centers. Then, the mode of the cluster number
distribution is chosen as the optimal cluster number. With the optimal
cluster number, the cluster center that minimizes the Bayesian
information criterion is selected, and items are assigned to the nearest
cluster based on Euclidean distances.
R > plot(lsirm_result, cluster = "spectral", chain.idx = 1)
Clustering result (Spectral Clustering):
group item
A 49, 50, 51, 52, 53, 54, 55, 56
B 33, 34, 35, 36, 37, 38, 39, 40
C 41, 42, 43, 44, 45, 46, 47, 48
D 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18,
19, 20, 21, 22, 23, 24
E 25, 26, 27, 28, 29, 30, 31, 32
R > plot(lsirm_result, cluster = "neyman", chain.idx = 1)
|==================================================| 100%
Clustering result (Neyman-Scott process):
group item
A 33, 34, 35, 36, 37, 38, 39, 40
B 49, 52, 54, 55
C 25, 26, 27, 28, 29, 30, 31, 32
D 17, 18, 19, 20, 21, 22, 23, 24
E 50, 51, 53, 56
F 9, 10, 11, 12, 13, 14, 15, 16
G 1, 2, 3, 4, 5, 6, 7, 8
H 41, 42, 43, 44, 45, 46, 47, 48 Figures 11 and 12 are item clustering results based on spectral clustering and the Neyman-Scott process model, respectively. In both plots, the gray dots indicate respondents, where numbers in colors indicate items with different cluster memberships. The Neyman-Scott process modeling approach additionally displays the center of the item cluster (alphabets A to H) and a contour for each item cluster. In this example, the spectral clustering method identifies five clusters, while the Neyman-Scott process modeling approach identifies eight clusters. Overall, the Neyman-Scott process modeling approach further split the items near the center (red items on the left) and in the north (light green items on the left) compared to the spectral clustering method.
Flexible Modeling Options
The lsirm12pl package provides a range of flexible modeling options for the LSIRM.
First, with fixed_gamma option, one can fix the distance
weight \(\gamma\) to a constant value
of 1. By doing so, the scale of the latent space can be standardized,
easing comparison between different interaction maps. The default value
of fixed_gamma option is FALSE (i.e., \(\gamma\) is freely estimated).
R > lsirm_result <- lsirm(data ~ lsirm1pl(fixed_gamma = TRUE))Second, with spikenslab option, one can apply a model
selection with the spike-and-slab prior for the distance weight \(\gamma\) (Ishwaran
and Rao 2005). The spike-and-slab prior is a mixture of two
log-normal priors: one is densely concentrated near zero, while another
is more broadly spread across positive values. With this option, one can
determine whether \(\gamma\) is zero or
not, which in turn determines whether the distance term is needed for
the data being analyzed. If \(\gamma\)
is not zero, it is evidence for item-by-respondent interactions in the
data being analyzed, and thus it would be useful to investigate the
interactions between respondents and items with the LSIRM approach. The
default value of spikenslab option is
FALSE.
The posterior probability of \(\gamma\) being non-zero is indicated by the
return value of pi_estimate. A pi_estimate
value greater than 0.5 suggests that \(\gamma\) is likely non zero. Note that the
two options fixed_gamma and spikenslab cannot
be used simultaneously.
R > lsirm_result <- lsirm(data ~ lsirm1pl(spikenslab = TRUE))
R > lsirm_result$pi_estimate
[1] 0.9984Third, the package offers options for handling missing data. In the
context of missing data, three key assumptions are considered (Rubin 1976): missing completely at random
(MCAR), missing at random (MAR), and missing not at random (MNAR). MCAR
occurs when the probability of missing data on a variable is unrelated
to any other measured variables or the variable itself, making the
missingness entirely random. MAR happens when the probability of missing
data on a variable is related to some of the observed data but not the
missing data itself, meaning the missingness can be explained by other
observed variables. Unlike MCAR and MAR, MNAR assumes that the missing
data mechanism depends on the unobserved data itself, making it
challenging to estimate the model without strong assumptions or
additional information about the missing data. Therefore, we focus on
two types of missing data, MCAR and MAR, in the lsirm12pl
package with the missing_data option.
With missing_data = "mcar", the missing data are assumed
to follow MCAR and the parameters are estimated solely based on the
observed elements of the dataset being analyzed. On the other hand, with
missing_data="mar", the missing data are assumed to be MAR,
and the data augmentation algorithm (Tanner and
Wong 1987) is applied to impute missing values. With
missing_data="mar", the function returns the posterior
samples of the imputed responses with imp and the
probability of a correct response with imp_estimate.
Imputed values are listed in the order of respondents. Note that all
missing values should be recoded via missing.val. The
percentage of missing data in the TDRI dataset is 30%, and missing
values are replaced with 99 which is the default coding of missing data.
The missing_data option can be used in combination with
other options such as
(spikenslab = TRUE, missing_data = "mar").
R > data <- lsirm12pl::TDRI
R > data[is.na(data)] <- 99
R > lsirm_result <- lsirm(data ~ lsirm1pl(missing_data = "mcar"))
R > lsirm_result <- lsirm(data ~ lsirm1pl(missing_data = "mar"))
R > lsirm_result$imp_estimate
[1] 0.9900 0.9868 0.9768 0.9884 0.9716 0.9716
...
[997] 0.9860 0.9788 0.8240 0.9924
[ reached getOption("max.print") -- omitted 32327 entries ]
R > plot(lsirm_result, option = "interaction")Figures 10(c) and 10(d) show the resulting interaction maps with MAR and MCAR assumptions, respectively. Note that a reflection is applied to the interaction map with MCAR assumptions (c) across the y-axis to ease comparisons with the other interaction plots. Reflection or rotation does not change the interpretations of the interaction map as interpretations are based on the distances, not the positions themselves. If the interaction maps of two missing assumptions are considerably different, further investigation is necessary to determine which assumption would be more suitable for the analyzed data. In this example, the interaction maps with the missing data options are similar to the original map presented in Figure 6. Note that the original interaction map (a) is based on the complete item response data with no missing values; that is, the data includes only respondents who answered all items. In contrast, the interaction maps with the missing data option are based on MAR and MCAR assumptions. The similarity between these maps suggests that the observed only or imputed item response data provides a reasonable representation of the original item response data.
2pl lsirm for dichotomous data
The two-parameter LSIRM (2PL LSIRM) extends the 1PL LSIRM with item discrimination parameters.
Statistical framework
The 2PL LSIRMs the probability of a correct response by respondent \(k\) to item \(i\) as follows: \[\text{logit}\Big(\mathbb{P}(Y_{k, i}=1 | \theta_{k}, \alpha_{i}, \beta_{i}, \gamma, \boldsymbol{z_{k}, w_{i}})\Big) = \alpha_{i} \theta_{k}+ \beta_{i} - \gamma d (\boldsymbol{z_{k}}, \boldsymbol{w_{i}}), \; \; \theta_{k} \sim N(0, \sigma^2),\] where \(\theta_k\), \(\beta_i\), \({\bf z}_k\), \({\bf w}_i\) and \(\gamma\) have similar interpretations to Equation @ref(eq:1pl-formula), while \(\alpha_i\) represents the item discrimination parameters and one of the \(\alpha_i\) parameters is fixed at 1, e.g., \(\alpha_{1} =1\) to ensure identifiability.
The observed data likelihood function under the 2PL LSIRM is given as \[\begin{aligned} \mathbb{L}\left(\boldsymbol{Y=y|\theta,\alpha,\beta,\gamma,Z,W}\right)=\prod_{k=1}^{N}\prod_{i=1}^{P}\mathbb{P}(Y_{k,i}=y_{k,i}|\theta_{k},\alpha_{i}, \beta_{i},\gamma,\boldsymbol{z_{k},w_{i}}). \label{likelihood_2pl} \end{aligned} (\#eq:likelihood-2pl)\] -1cm
Parameter Estimation
Following is the posterior distribution of the 2PL LSIRM. \[\begin{aligned}
\pi\left(\boldsymbol{\theta, \beta}, \gamma, \boldsymbol{Z, W} |
\boldsymbol{Y=y} \right) & \propto
\prod_{k=1}^{N}\prod_{i=1}^{P}\mathbb{P}(Y_{k,i}=y_{k,i}|\theta_{k},\alpha_{i},
\beta_{i},\gamma,\boldsymbol{z_{k},w_{i}}) \\
&\times \prod_{k=1}^N \pi(\theta_k) \times \prod_{i=1}^P
\pi(\alpha_i) \times \prod_{i=1}^P \pi(\beta_i) \times \pi(\gamma)
\times \prod_{k=1}^N \pi(\boldsymbol{z_{k}}) \prod_{i=1}^P
\pi(\boldsymbol{w_{i}})
\end{aligned}\] The prior distributions for \(\theta_k\), \(\beta_i\), \({\bf
z}_k\), \({\bf w}_i\), and \(\gamma\) are identical for the 1PL LSIRM.
For item discrimination parameters \(\alpha\), a log-normal distribution is used
with a mean of \(\mu_{\alpha}\) and a
variance of \(\tau_{\alpha}^2\), as
\(\alpha\) is typically assumed to be
positive. The argument names and default values for the prior
specification are shown in the Github site. For the item discrimination
parameters \(\alpha\), the arguments
and default values are
pr_mean_alpha = 0.5, pr_sd_alpha = 1.
The conditional posterior distribution for each parameter follows the
same form as the Equation in the Github site. The jumping rule for each
parameter is given in the Github site. The default jumping rule for
\(\alpha\) is
jump_alpha = 1.
An Illustrated Example
We apply the 2PL LSIRM to the TDRI data. The default settings of the 2PL LSIRM are the same as the 1PL LSIRM except for \(\alpha\). The base function for the 2PL LSIRM is
lsirm(A\(\sim\)lsirm2pl())
The 2PL LSRIRM for dichotomous data was fitted with 4 MCMC chains
using 2 multi-core processors, as demonstrated in the following example
code. Similarly to the lsirm1pl results, the estimation
results for the model parameters \(\boldsymbol{\theta,\ \beta,\ \alpha}\),
\(\gamma,\ \boldsymbol{Z}\), and \(\boldsymbol{W}\) for each chain are
provided in individual lists.
R > library("lsirm12pl")
R > data <- lsirm12pl::TDRI
R > data <- data[complete.cases(data),]
R > lsirm_result <- lsirm(data ~ lsirm2pl(chains = 4, multicore = 2, seed = 2025))To diagnose the results of the 2PL LSIRM analysis, we can use the
diagnostic() function. Similarly to the diagnostic of the
1PL LSIRM, the convergence of MCMC for each parameter can be checked
using various diagnostic tools, such as trace plots, posterior density
distributions, autocorrelation functions (ACF), and Gelman-Rubin-Brooks
plots. By setting the draw.item option to
list(’beta’= c(1)), we perform diagnostics for the \(\beta_i\) parameter of the first item
(\(i=1\)).
R > diagnostic(lsirm_result,
draw.item = list(beta = c(1)),
gelman.diag = T)
Potential scale reduction factors:
Point est. Upper C.I.
beta [i1] 1 1
diagnostic() function on the results of the 2PL LSIRM
fitted to the TDRI dataset. The top left is the trace plot, the top
right is the posterior density plot, the bottom left is the
autocorrelation plot, and the bottom right is the Gelman-Rubin-Brooks
plot. Different colors represent different MCMC chains.Figure 14 displays the diagnostic results for
\(\beta_1\) of the results of the 2PL
LSIRM fitted to the TDRI dataset. These plots help us assess the
convergence of MCMC for \(\beta_1\).
The interpretation of each plot from the diagnostic()
function is the same as mentioned for the 1PL LSIRM in the previous
Section. In Figure 14, the convergence of \(\beta_1\) was confirmed by various
diagnostic tools.
The gof() function assesses the goodness-of-fit of the
LSIRM. Figure 15 visualizes the box plots of the posterior
predicted response values (item-wise means) compared with the observed
item-wise means and the ROC curve to check the performance of the 2PL
LSIRM fitted to the TDRI dataset. In the figure, most of the red dots
are located close to the midline of the boxplots and the AUC is 0.97,
indicating the fit of the 2PL LSIRM to the TDRI dataset is
satisfactory.
R > gof(lsirm_result, chain.idx = 1)
gof() function. The box plots of
the average posterior predicted response values for items against the
average observed response value for each item indicated by red dots
(left) and the receiver operating characteristic (ROC) curve
(right).Similarly to the 1PL LSIRM, the plot() function can
generate different graphs for the 2PL LSIRM results with the
option argument:
R > plot(lsirm_result, option = "beta", chain.idx = 1)
R > plot(lsirm_result, option = "theta", chain.idx = 1)
R > plot(lsirm_result, option = "alpha", chain.idx = 1) 
plot() function on the results of the 2PL LSIRM fitted to
the TDRI dataset. (a) Boxplots of the posterior samples for βi. (b) Boxplots
of the point estimate for θk as a function
of the total sum scores of positive responses. (c) Box plots of the
posterior samples for αi.
Figure 19 illustrates the results of the
plot() function with the "beta",
"theta", and "alpha" options for the 2PL
LSIRM, respectively. Similarly to the 1PL LSIRM, \(\beta_i\) decreases for later items,
indicating that later items are more difficult than earlier items. In
addition, respondents with more correctly answered items (that is,
higher sum scores) have higher \(\theta_k\) estimates. Lastly, the
distribution of the posterior samples for \(\alpha_i\) is similar across all items.
The plot() function can be used to create a
visualization of the interaction map based on the 2PL LSIRM by setting
option argument as "interaction". Figures 20 and 21 show the original and rotated
interaction maps, respectively. In the interaction maps based on the 1PL
and 2PL LSIRM, we notice that although the clustering of items is
similar, the distances between item groups appear smaller in the
interaction map of the 2PL LSIRM. This difference arises because the 2PL
LSRIM takes item discrimination (\(\alpha\)) into account, which explains some
degree of item-by-person interactions. The interaction map with the
oblimin rotation (b) shows a slight clockwise rotation compared to the
original interaction map (a), similarly to the 1PL LSIRM’s case.
R > plot(lsirm_result, option = "interaction", chain.idx = 1)
R > plot(lsirm_result, option = "interaction", rotation = TRUE, chain.idx = 1)
The plot() function visualizes the cluster of the latent
positions of items by using spectral clustering and the Neyman-Scott
process modeling approach, achieved by setting the cluster
option to spectral and neyman,
respectively.
R > plot(lsirm_result, cluster = "spectral", chain.idx = 1)
Clustering result (Spectral Clustering):
group item
A 49, 50, 51, 52, 53, 54, 55, 56
B 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40
C 41, 42, 43, 44, 45, 46, 47, 48
D 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16
E 17, 18, 19, 20, 21, 22, 23, 24
R > plot(lsirm_result, cluster = "neyman", chain.idx = 1)
|==================================================| 100%
Clustering result (Neyman-Scott process):
group item
A 33, 34, 35, 36, 37, 38, 39, 40
B 49, 50, 51, 52, 53, 54, 55, 56
C 41, 42, 43, 44, 45, 46, 47, 48
D 17, 18, 19, 20, 21, 22, 23, 24
E 25, 26, 27, 28, 29, 30, 31, 32
F 1, 2, 3, 4, 5, 6, 7, 8
G 9, 10, 11, 12, 13, 14, 15, 16 
Figure 23 and 24 display the result of spectral clustering and the Neyman-Scott process approach, respectively. The interpretation of gray dots, numbers, alphabets, and contours is the same as the 1PL LSIRM case. Spectral clustering resulted in 5 clusters, whereas the Neyman-Scott process approach resulted in 7 clusters. That is, the Neyman-Scott process approach identified more item groups than the spectral clustering. The overall clustering results are similar to the 1PL LSIRM case.
The flexible modeling options discussed with the 1PL LSIRM can be
applied to the 2PL LSIRM. That is, fixed_gamma fixes the
distance weight \(\gamma\) to 1 and
spikenslab assigns a spike and slab prior to \(\gamma\). Additionally, the two missing
data options, missing_data = "mcar" and
missing_data = "mar", are available for the 2PL LSIRM.
R > lsirm_result <- lsirm(data ~ lsirm2pl(fixed_gamma = TRUE))
R > lsirm_result <- lsirm(data ~ lsirm2pl(spikenslab = TRUE))
R > lsirm_result <- lsirm(data ~ lsirm2pl(missing_data = "mcar"))
R > lsirm_result <- lsirm(data ~ lsirm2pl(missing_data = "mar"))Lsirm for continuous item responses data
We consider an extension of the LSIRM for continuous item response data (LSIRM-continuous), which is done by using an appropriate link function following the generalized linear model framework (McCullagh and Nelder 2019).
Statistical framework
Consider the continuous item response data consisting of the \(N \times P\) matrix \(\boldsymbol{Y} = \left\{ y_{k,i} \right\} \in \mathbb{R}^{N\times P}\). By choosing the identity link function for continuous item responses, the 1PL LSIRM-continuous is given as follows: \[\begin{aligned} y_{k,i} \mid \boldsymbol{\theta}, \boldsymbol{\beta}, \boldsymbol{\gamma},\boldsymbol{Z}, \boldsymbol{W},\sigma_{\epsilon}^2 = \theta_{k}+\beta_{i}-\gamma d (\boldsymbol{z_{k}},\boldsymbol{w_{i}}) + \epsilon_{k,i}, \; \; \theta_{k} \sim N(0, \sigma^2), \; \; \epsilon_{k,i} \sim N (0, \sigma_{\epsilon}^2). \end{aligned}\] where \(\theta_k\) and \(\beta_i\) are the main effects of the respondent \(k\) and the item \(i\), respectively. \(\boldsymbol{z_k}\) and \(\boldsymbol{w_i}\) are the latent positions of the respondent \(k\) and the item \(i\), respectively. An additional error term \(\epsilon_{k,i} \sim N(0, \sigma_{\epsilon}^2)\) explains residuals unexplained by the 1PL LSIRM-continuous. A shorter distance between \(\boldsymbol{w_i}\) and \(\boldsymbol{z_k}\) indicates that the respondent \(k\) is likely to give a higher response value to item \(i\) given the main effects of the respondent and the item.
It is straightforward to extend the 1PL LSIRM-continuous to the 2PL version by adding item discrimination (or slope) parameters \(\alpha_i\). The 2PL LSIRM-continuous for \(y_{k,i}\) is given as follows: \[\begin{aligned} y_{k,i} \mid \boldsymbol{\theta},\boldsymbol{\alpha}, \boldsymbol{\beta}, \boldsymbol{\gamma},\boldsymbol{Z}, \boldsymbol{W},\sigma_{\epsilon}^2 = \alpha_{i}\theta_{k}+\beta_{i}-\gamma d (\boldsymbol{z_{k}},\boldsymbol{w_{i}}) + \epsilon_{k,i}, \; \; \theta_{k} \sim N(0, \sigma^2), \; \; \epsilon_{k,i} \sim N (0, \sigma_{\epsilon}^2). \end{aligned}\] The interpretations of the model parameters in the 2PL LSIRM-continuous are similar to the case of the 1PL LSIRM-continuous and the 2PL LSIRM. For model identifiability, one of the item slopes is fixed to 1, e.g., \(\alpha_{1} =1\).
The likelihood function of the 1PL LSIRM-continuous is given as \[\begin{aligned} \mathbf{L}\left(\boldsymbol{Y=y|\theta,\beta},\gamma,\boldsymbol{Z,W},\sigma_{\epsilon}^2\right) & =\prod_{k=1}^{N}\prod_{i=1}^{P}\mathbf{L}(Y_{k, i}=y_{k, i}|\theta_{k}, \beta_{i},\gamma,\boldsymbol{z_{k},w_{i}},\sigma_{\epsilon}^2), \\ & = \prod_{k=1}^{N}\prod_{i=1}^{P} N(y_{k, i,};\theta_k +\beta_i - \gamma || {\boldsymbol{z_k}} - {\boldsymbol{w_i}} ||, \sigma_{\epsilon}^2), \end{aligned}\] and the likelihood function of the 2PL LSIRM-continuous is given as \[\begin{aligned} \mathbf{L}\left(\boldsymbol{Y=y|\theta,\alpha,\beta},\gamma,\boldsymbol{Z,W},\sigma_{\epsilon}^2\right) & = \prod_{k=1}^{N}\prod_{i=1}^{P} N(y_{j, i};\alpha_i\theta_k +\beta_i - \gamma|| {\boldsymbol{z_k}} - {\boldsymbol{w_i}} ||, \sigma_{\epsilon}^2). \end{aligned}\]
-0.75cm
Parameter Estimation
Following is the posterior distribution of LSRIM-continuous. \[\begin{aligned} \pi\left(\boldsymbol{\theta, \beta}, \gamma, \boldsymbol{Z, W}, \sigma_{\epsilon}^2 | \boldsymbol{Y=y} \right) &\propto \prod_{k=1}^{N}\prod_{i=1}^{P}N(y_{k, i,};\theta_k +\beta_i - \gamma || {\boldsymbol{z_k}} - {\boldsymbol{w_i}} ||, \sigma_{\epsilon}^2) \\ &\times \prod_{k=1}^N \pi(\theta_k) \times \prod_{i=1}^P \pi(\beta_i) \times \pi(\gamma) \times \prod_{k=1}^N \pi(\boldsymbol{z_{k}}) \times \prod_{i=1}^P \pi(\boldsymbol{w_{i}}) \times \pi(\sigma_{\epsilon}^2) \end{aligned}\] The priors of the model parameters in the LSIRM-continuous are set the same as the priors for the LSIRM. For the error term \(\sigma_{\epsilon}^2\), we set the prior as follows: \[\begin{aligned} \sigma_{\epsilon}^{2}|a_{\sigma_{\epsilon}},b_{\sigma_{\epsilon}} & \sim\text{Inv-Gamma}(a_{\sigma_{\epsilon}},b_{\sigma_{\epsilon}}),\ a_{\sigma_{\epsilon}}>0,\ b_{\sigma_{\epsilon}}>0. \end{aligned}\] The conditional posterior distributions of the LSIRM-continuous are similar to the LSIRM, which are given in the Github site. The jumping rule defaults remain the same as in the binary case (shown in the Github site).
An Illustrated Example
We demonstrate the application of the LSIRM-continuous using the Big Five Personality Test (FPT) dataset (Goldberg 1992), which is included in the package. Data were collected from 1,015,342 respondents using an interactive online personality test from 2016 to 2018. The “Big-Five Factor Markers” from the international personality item pool were used in the test, which consists of 50 questions with response categories 1 = disagree, 3 = neutral, and 5 = agree based on a five-point Likert scale. Negatively worded items were reverse-coded. For illustration purposes, we randomly selected a sample of 3,000 respondents from the original data. We treated the ordinal item responses as continuous data, which is a common practice in applied research and is generally considered acceptable for other models, such as factor analysis.
The main fitting function for running the 1PL and 2PL LSIRM-continuous is identical to the 1PL and 2PL LSIRM for binary data. The function automatically identifies the data type and applies the appropriate models. Following is the code for data pre-processing and 1PL LSIRM-continuous fitting.
R > data <- lsirm12pl::BFPT
R > data[(data==0)|(data==6)] = NA
R > reverse <- c(2, 4, 6, 8, 10, 11, 13, 15, 16, 17,
18, 19, 20, 21, 23, 25, 27, 32, 34,
36, 42, 44, 46)
R > data[, reverse] <- 6 - data[, reverse]
R > data <- data[complete.cases(data),]
R > head(data)
EXT1 EXT2 EXT3 ... OPN8 OPN9 OPN10
1 2 3 2 ... 2 4 4
2 1 3 3 ... 4 3 4
3 4 3 3 ... 2 4 4
5 1 4 3 ... 3 5 3
6 1 3 2 ... 3 5 3
8 3 5 4 ... 4 3 4
R > lsirm_result <- lsirm(data ~ lsirm1pl(niter = 25000, nburn = 5000, nthin = 10,
jump_beta = 0.08, jump_theta = 0.3,
jump_gamma = 1.0,
chains = 4, multicore = 2, seed = 2025))To ensure convergence of the parameters, a different number of
iterations and jumping rules of parameters are applied in this example.
The estimated results for each chain are returned in the list containing
estimated information on \(\boldsymbol{\theta,\ \beta},\ \gamma,\
\boldsymbol{Z,\ W},\ {\sigma}\), and \(\sigma_\epsilon\). The functions
summary(), diagnostics(), and
gof() described for the 1PL LSIRM can be applied similarly
to obtain a summary, diagnosis, and goodness-of-fit results,
respectively.
R > diagnostic(lsirm_result, draw.item = list(beta = c("AGR1")))
diagnostic() function on the results of the 1PL
LSIRM-continuous fitted to the FPT dataset. The top left is the trace
plot, the top right is the posterior density plot, the bottom left is
the autocorrelation plot, and the bottom right is the
Gelman-Rubin-Brooks plot. Different colors represent different MCMC
chains.Figure 26 displays the diagnostic results
for \(\beta_{AGE_1}\) obtained using
the diagnostic() function. The results suggest the
convergence of \(\beta_{AGE_1}\) is
achieved for this model.
Figure 27 shows the result of the goodness of fit
assessment for the continuous example data. Unlike binary data, ROC is
not available for continuous data, so the results of the
gof() function include only the boxplots of the average
predicted response values (item-wise means) against the observed
item-wise means (marked with red dots). In this example, the red dots
are located close to the midlines of the boxplots, implying a
satisfactory model fit.
R > gof(lsirm_result, chain.idx = 1)
gof()
function. The box plots of the average posterior predicted response
value for items against the average observed response value for each
item indicated by red dots.Same as before, the plot() function can be used to
summarize the parameter estimate of \(\beta_i\) and \(\theta_k\).
R > plot(lsirm_result, option = "beta", chain.idx = 1)
R > plot(lsirm_result, option = "theta", chain.idx = 1) Figure 28 shows the boxplots of the posterior samples for \(\beta_i\). The parameter estimate of \(\beta_{35}\) is the smallest , while the estimates for \(\beta_{41}\) to \(\beta_{50}\) are relatively high. Figure 29 shows the boxplots of the point estimates of \(\theta_k\) as a function of the sum scores of the observed responses binned into 10 groups to prevent overlap and improve readability on the x-axis. As expected, higher sum scores are aligned with higher \(\theta_k\) values (indicating socially desirable characteristics).
plot() function on the 1PL LSIRM-continuous fitted to the
FPT dataset. (a) Boxplots of the posterior samples for βi. (b) Boxplots
of the point estimates of θk as a function
of the total sum scores of the responses binned into 10 groups
Figure 31 illustrates the interaction map derived from the 1PL LSIRM-continuous. In the interaction map, the latent positions of the respondents are positioned around the center, while the items are scattered in a triangular pattern showing roughly three clusters (in the north, west, and east) in addition to the cluster around the center. The original interaction map is slightly rotated counterclockwise in the rotated interaction map, so that the items are more closely placed to the two coordinates.
R > plot(lsirm_result, chain.idx = 1)
R > plot(lsirm_result, rotation = TRUE, chain.idx = 1)
Figure 34 and Figure 35 depict the result of spectral clustering and the Neyman-Scott process modeling approach for the BFPT example dataset, respectively. Gray dots, numbers with colors, alphabets, and contours have the same interpretations as the clustering results presented earlier with the binary LSIRM models. In Figure 35, the items in cluster C, highlighted in purple and located near the center of the interaction map, are closer to the latent positions of most people compared to other items. This implies that these items are more likely to receive higher response values. Conversely, items in clusters that are positioned farther from the center, such as clusters A, B, and E, are distant from many (but different groups of) respondents, indicating they are likely to receive lower response values by those who are far away from the corresponding item clusters.
R > plot(lsirm_result, cluster = "spectral", chain.idx = 1)
Clustering result (Spectral Clustering):
group item
A 11, 13, 16, 17, 18, 20
B 38, 41, 42, 43, 44, 45, 46, 48, 49, 50
C 1, 2, 3, 4, 5, 7, 8, 9, 10
D 6, 21, 22, 24, 25, 26, 27, 28, 29, 30
E 12, 14, 15, 19, 23, 31, 32, 33, 34, 35,
36, 37, 39, 40, 47
R > plot(lsirm_result, cluster = "neyman", chain.idx = 1)
|==================================================| 100%
Clustering result (Neyman-Scott process):
group item
A 38, 41, 42, 43, 44, 45, 46, 48, 49, 50
B 1, 2, 3, 4, 5, 7, 8, 9, 10
C 6, 21, 22, 24, 25, 26, 27, 28, 29, 30
D 23, 31, 32, 33, 34, 35, 36, 37, 39, 40, 47
E 11, 12, 13, 14, 15, 16, 17, 18, 19, 20 
The flexible modeling options discussed with the binary LSIRM can be applied to the functions for the LSIRM-continuous.
R > lsirm_result <- lsirm(data ~ lsirm1pl(fixed_gamma = TRUE))
R > lsirm_result <- lsirm(data ~ lsirm1pl(spikenslab = TRUE))
R > lsirm_result <- lsirm(data ~ lsirm1pl(missing_data = "mcar"))
R > lsirm_result <- lsirm(data ~ lsirm1pl(missing_data = "mar"))Conclusion
In this paper, we introduced the R package lsirm12pl for estimating the LSIRM (Jeon et al. 2021) and its extensions. The LSIRM is a powerful extension of conventional IRT models that allow for the estimation and visualization of potential interactions between respondents and items in an interaction map, a low dimensional Euclidean space. The original LSIRM framework was proposed for binary item response data based on the 1PL IRT base model. To broaden its applicability, we extended the original LSIRM to cover different response types and model specifications. Further, we added several useful options, e.g., for handling missing data, item clustering, and model assessment.
It is worth noting that one may observe some variability in clustering results across multiple chains. Such variability stems from at least three factors: First, clustering is an unsupervised method, meaning that results are not uniquely determined and may naturally vary across Markov chains. Second, we apply the Neyman-Scott Process clustering to the posterior distributions of the LSIRM model, and the clustering outcomes are therefore directly influenced by the estimated LSIRM results. Third, Bayesian models inherently incorporate uncertainty, reflected in the posterior distributions and subsequently in the clustering outcomes. Having said that, we noted in our empirical analysis that the overall clustering structure was pretty stable across multiple chains. For critical applications, we recommend implementing a consensus clustering approach by aggregating results across multiple chains to identify the most robust and stable groupings.
A fully Bayesian approach was used for model estimation with a Metropolis-Hastings-within-Gibbs sampler. The lsirm12pl package offers default estimation settings for priors, jumping rules, number of iterations, burn-ins, and thinning. The default estimation setting works reasonably well in a broad range of situations, but users can manually revise the estimation settings if desired. In addition, the package lsirm12pl offers convenient supplemental functions to evaluate, summarize, visualize, diagnose, and interpret the estimated results. We provided detailed illustrations of the package with real data examples available in the package. We hope that these illustrations guide researchers in using the lsirm12pl package for the analysis of their own datasets.
The R package lsirm12pl
is our first step in making the LSIRM approach more applicable and
usable in practice. The code is written using Rcpp
(Eddelbuettel and François 2011; Eddelbuettel
2013; Eddelbuettel and Balamuta 2018) and RcppArmadillo
(Eddelbuettel and Sanderson 2014) in R for
efficient computation. For example, in our first data example with 726
respondents and 56 test items, the computation time for one chain using
a single core lsirm1pl was 1.93 minutes on an Apple M1
laptop. We provide additional details of our package implementation in
the package GitHub site (https://github.com/jiniuslab/lsirm12pl).
We will continue to update the package by incorporating additional
modeling, data analysis, and visualization options to make the lsirm12pl
package more useful in a wider range of situations. For example, we are
currently engaged in the development of other extensions, such as for
ordinal and longitudinal data. As advancements are made, we will
systematically include them in the software package. This ongoing
development will further enhance the utility and applicability of the
LSIRM in various practical scenarios.
Acknowledgement
We thank the Editor, Associate Editor, and reviewers for their constructive comments. We also thank Dr. Won Chang for constructive comments on our work. This work was partially supported by the National Research Foundation of Korea [grant number NRF-2020R1A2C1A01009881, RS-2023-00217705, RS-2024-00333701; Basic Science Research Program awarded to IHJ], [grant number NRF-2025S1A5C3A0200632411; awarded to IHJ and MJ], and the ICAN (ICT Challenge and Advanced Network of HRD) support program [grant number RS-2023-00259934], supervised by the IITP (Institute of Information & Communications Technology Planning & Evaluation). Correspondence should be addressed to Ick Hoon Jin, Department of Applied Statistics, Department of Statistics and Data Science, Yonsei University, Seoul, Republic of Korea. Go, Kim, and Park are co-first authors and listed in alphabetical order.