Bayesian Estimation of Latent Space Item Response Models with JAGSJAGSJAGS , StanStanStan , and NIMBLENIMBLENIMBLE in RRR

: The latent space item response model (LSIRM) is a newly-developed approach to analyzing and visualizing conditional dependencies in item response data, manifested as the interactions between respondents and items, between respondents, and between items. This paper provides a practical guide to the Bayesian estimation of LSIRM using three open-source software options, JAGS , Stan , and NIMBLE in R. By means of an empirical example, we illustrate LSIRM estimation, providing details on the model speciﬁcation and implementation, convergence diagnostics, model ﬁt evaluations


Introduction
Item Response Theory (IRT) is widely used for item response analysis in various fields, including education, psychology, medicine, and other social sciences, where assessments measure unobserved constructs, such as cognitive ability, personality, and mental well being.IRT provides a way to model the relationship between individuals' responses to test items and the underlying construct being measured by the test.Conventional IRT models are based on a set of strong assumptions, such as conditional independence (i.e., responses are independent of one another, given the latent trait).However, conditional independence may not hold in real-life data analysis when unobserved interactions between respondents and items exist and are conditional on the latent trait.
Recently, Jeon et al. [1] introduced a novel latent space item response modeling (LSIRM) framework that relaxes the conditional independence assumption.LSIRM captures conditional dependencies unaccounted for by conventional IRT models [2] in terms of respondentby-item interactions, where the interactions are modeled in the form of distances between respondents and items in a low-dimensional geometric space, called an interaction map.An estimated interaction shows conditional dependencies, or deviations from the Rasch model, for items and respondents, helping in the identification of intended or unintended response behaviors or in the evaluation of item grouping structures in comparison to the intended test design.In summary, the LSIRM approach has substantial technical advantages, due to its relaxing of the conditional independence assumption, and, by producing a geometrical representation of unobserved item-respondent interactions, it may supply important insights into the respondents and test items under investigation.
Jeon et al. [1] described a fully Bayesian approach with Markov chain Monte Carlo (MCMC) for estimating the original LSIRM.To facilitate the wider adoption and application of LSIRM in research and practice, we aim to provide a practical guide to LSIRM estimation using commonly used, free of charge, Bayesian software, such as JAGS [3][4][5], Stan [6,7], and NIMBLE [8,9] in R.
This paper begins with a brief description of LSIRM for binary response data and LSIRM parameter estimation with MCMC.We provide detailed R syntax for model specification and implementation in JAGS, Stan, and NIMBLE with a real-life data example.We also provide convergence and model fit diagnosis, run-time/sampling efficiency, and interaction map visualization comparison across the three packages.We conclude the paper with a summary and brief discussion.

Background 2.1. Latent Space Item Response Model
The original latent space item response model (LSIRM) [1] extended the Rasch model by introducing a respondent-by-item distance term from a p-dimensional space R p .LSIRM assumes that respondents and items have positions in the shared latent space R p .The distance between the respondent and item positions in the space indicates the interaction between the respondent and the item unexplained by the person and item parameters of the Rasch model.Specifically, LSIRM specifies the probability of giving a correct response by respondent j to item i as logit(P (Y ji = 1|θ j , β i , γ, z j , w i )) = θ j + β i − γ d(z j , w i ), where Y = {y ji } is a N by I item response matrix, and y ji is the response of respondent j to item i (j = 1, 2, . . ., N and i = 1, 2, . . ., I).When Y are all binary responses, we use y ji = 1 to denote an affirmative endorsement (e.g., correct, true, or yes) to the item, while y ji = 0 denotes no endorsement (e.g., incorrect, false, or no).θ j ∼ N(0, σ 2 ) is the person intercept or the main effect of respondent j that can be interpreted as respondent j's latent trait being measured by the test.The value β i is the item intercept or the main effect of item i that can be interpreted as the item i's easiness.The value γ > 0 is the weight of the distance term that captures the overall degree of conditional dependencies in the item response data; a larger γ indicates stronger evidence of conditional dependencies (i.e., violation of the conditional independence assumption) in the item response data under investigation, When γ = 0, the model reduces to the Rasch model.d(z j , w i ) is the distance function that measures pairwise distances between a respondent's latent position z j and an item's latent position w i in the p-dimensional Euclidean space R p with known dimension p.
To facilitate visualization and interpretation, in this paper we use a two-dimensional Euclidean space (p = 2) and Euclidean distance d(z j , w i ) = d(z j , w i ) = z j − w i , where a larger distance between z j and w i contributes to a lower probability of person j giving a correct response to item i, given the person's overall trait level and the item's overall easiness.As such, a distance between the respondent and item positions indicates a respondent-by-item interaction.When there is little interaction between respondents and items, all distances are near zero (i.e., all person and item positions are close to (0, 0) in a R 2 space), as shown in [1].In the sense that the latent space represents unobserved interactions, it is referred to as an interaction map [1].Note that the interaction map is not an ability space; i.e., the coordinates of the space do not represent ability dimensions.Instead, the latent space is a tool to represent respondent-by-item interactions unexplainable by the parameters of the Rasch model.Other additional details of the model specifications, assumptions, and interpretations can be found in the original paper [1].
The likelihood function of the observed item response data y can be written as where θ = {θ j }, β = {β i }, Z = {z j } and W = {w i }, j = 1, 2, . . ., N and i = 1, 2, . . ., I. The likelihood function assumes that the item responses are independent, given the person and item main effects and the distances (or interactions) between the respondent and item latent positions.This is a relaxation of the conditional independence assumption required by the Rasch model (that assumes no person-by-item interaction effects).

Bayesian Estimation of LSIRM
For estimating LSIRM, Jeon et al. [1] applied a fully Bayesian approach with MCMC using Metropolis-Hasting and Gibbs sampling.The priors for the model parameters are set as follows: where 0 is a vector of zeros of length p and I p is the p-dimensional identity matrix.The priors of β i are set to be normal, and z j and w i are set to be multivariate normal with no covariances.The prior for the weight parameter γ > 0 is set to be log-normal.The prior for the variance parameter of the θ j is set to be inverse-Gamma.
The posterior of the model parameters is proportional to where f (•) denotes the prior and posterior probability density functions of the parameters and latent positions of respondents and items.For additional details on the MCMC sampling, we refer the reader to Jeon et al's paper [1].

Interaction Map
A unique advantage of the LSIRM approach is that it produces a visualization of item-by-person interactions in the form of distances in a low-dimensional geometric space, or an interaction map.Items and respondents are positioned on the interaction map, where the distances indicate the conditional dependencies unaccounted for by the Rasch model.Of note, one can also interpret distances between items and distances between respondents due to the triangular inequality property (even though the model does not explicitly model the item-item and respondent-respondent distances).Therefore, distances in an interaction map indicate unobserved similarities and dissimilarities between items, between respondents, and between respondents and items.
Two remarks can be made here in terms of producing and interpreting an interaction map.First, the distances are invariant to translations, reflections, and rotations of the latent positions [10,11].This means that the latent positions are not identifiable and there can be many latent space configurations that produce identical distances; thus, one should make interpretations of distances, but not particular positions of items or persons on the map.Second, because of the unidentifiability issue of latent positions, the posterior draws of the person and item positions (w and z) are not comparable across iterations.To match the latent positions, we post-process the MCMC output by applying Procrustes matching [12] to the posterior draws of the person and latent positions, by using the iteration that provides the best likelihood as the reference point [1].

JAGS
JAGS, short for Just Another Gibbs Sampler, is a software program developed by Martyn Plummer for conducting Bayesian inference using the Gibbs Sampling algorithm and BUGS language [5].This program is one of the most widely-used Bayesian inference software in the world, with over 200,000 downloads worldwide since its release [16].The jagsUI package is a set of wrappers around the rjags package to run Bayesian analyses with JAGS (specifically, via libjags).A benefit of this package is that it can automatically summarize posterior distributions and generate plots for assessing the model convergence (e.g., predictive check plots and trace plots) [3].Additionally, this package can output posterior means, standard deviations, and percentile intervals for each parameter, which is convenient for further output analysis and result interpretation.
Currently, JAGS is widely utilized for the estimation of latent variable models, such as IRT models [17], cognitive diagnosis models [18], latent class models [19], structural equation models [20,21], etc.In recent years, JAGS has also increasingly been applied for social network model estimation [22].However, JAGS has rarely been used for latent space models [10], let alone LSIRM, to the best of our knowledge.

Stan
Stan, a software program for performing Bayesian inference, is named after Stanislaw Ulam, a mathematician who helped develop the Monte Carlo method in the 1940s [23].Developed in 2012, Stan promises computational and algorithmic efficiency, which is mainly thanks to the No-U-Turn Sampler (NUTS; [24]), an adaptive variant of Hamiltonian Monte Carlo (HMC; [25]) used within the program.HMC is a generalization of the Metropolis algorithm that allows for more efficient movement through the posterior distribution by performing multiple steps per iteration.Stan is based on C++ and can be run from the command line, R, Python, Matlab, or Julia.Although Stan performs similarly to other Bayesian programs, such as JAGS, for simple models, it reportedly outperforms other software as model complexity grows [26] and scales better for large data sets [6], although JAGS seems to earn more MCMC efficiency from conjugate priors than Stan [27].Stan has been used for latent variable model estimation, but not often for latent space models.Recently, ref. [28] used Stan for latent space model estimation, reporting that HMC retains as much of the posterior dependency structure as possible compared to the variational Bayes approach [29].Taken together, Stan's efficiency and flexibility in handling complex models make it a useful tool for Bayesian inference for LSIRM.

NIMBLE
NIMBLE [8], short for Numerical Inference of Statistical Models for Bayesian and Likelihood Estimation, is another software program for Bayesian inference.NIMBLE allows a combination of high-level processing in R and low-level processing in compiled C++, which helps write fast numerical functions with or without the involvement of BUGS models.Additionally, NIMBLE adopts and extends the BUGS language for specifying models, allowing users to easily transition from BUGS to NIMBLE.The efficiency of NIMBLE depends on the specific model and priors being used.In some cases, NIMBLE may be faster and produce better quality chains than JAGS and Stan, while in some other cases, packages such as Stan may be more efficient [30].NIMBLE has been used for semi-parametric and non-parametric Bayesian IRT models [31,32] and other latent variable models [33].NIMBLE has not been used for estimating LSIRM or latent space models so far.

Estimating LSIRM with JAGS JAGS JAGS, Stan Stan Stan, and NIMBLE NIMBLE NIMBLE
Here we illustrate the Bayesian estimation of LSIRM with JAGS, Stan, and NIMBLE.We applied the priors established in Section 2.2 in all data analysis, with the following values: Users may consider a different set of priors based on prior knowledge and specific assumptions about the data under investigation.Prior choice is a critical task in Bayesian inference, and one may want to consider sensitivity analysis to validate the appropriateness of the chosen priors [34][35][36][37] .We used 1000 thinned MCMC iterations in all analyses.For Stan, we used a burn-in period of 10,000 and a thinning interval of 5, resulting in a total number of iterations of 15,000.For comparable results in effective sample sizes of model parameters, we used a burn-in period of 40,000 and a thinning period of 50 for NIMBLE and JAGS.Note that the optimal number of iterations, burn-in period, and thinning intervals vary for different data sets.All of the analysis was done within the R environment on the Hoffman2 Linux Compute Cluster hosted at UCLA and each run was conducted within a node with two cores and 16GB RAM.
In Listing 1, we provide the code for loading and the verbal data for analysis.The input data is a N by I response matrix, where N is the number of respondents, and I is the number of items.

LSIRM Estimation in JAGS Model Specification
Listing 2 presents JAGS code that can be used to fit LSIRM through a .bugfile using the BUGS language.The code is a realization of Equation ( 1) with the priors specified in Equation (3).
• Line 10 and line 15 specify bivariate normal priors for the person and item latent positions, respectively, assuming a two-dimensional latent space.The specification can be adjusted if the model specifies a higher-dimensional latent space.• Lines 19 and 20 set the variance parameter of respondents' latent trait θ j to follow an Inverse Gamma distribution with the function invsigma ˜dgamma(a_sigma, b_sigma) .
We then use the reciprocal of invsigma to obtain the desired Inverse Gamma distribution for sigma2 .
• Lines 34 to 36 save the log-likelihood.

Run MCMC in JAGS
The following code in Listing 3 can be used to run the JAGS model specification in R.
Here we list some comments to clarify the code.Note that in Stan the codes are defined in blocks, and, unlike BUGS, the declarations and statements are executed in order.
• Lines 1 to 13 display the data block, where data and hyperparameters are defined.Note that it is mandatory to specify the data type, such as: int as an integer, real as a real number, cov_matrix as a covariance matrix.This follows the convention of programming languages like C++. • Lines 14 to 21 contain the parameters block, where the parameters to estimate are defined with the type and range of values.

•
The range of values that each object can take can be defined inside the <> symbol.
For instance, we know that the number of respondents and items, tau_gamma and tau_beta can only be positive, so the lower bound is set to 0, while, for responses, the lower bound is set to 0 and the upper bound to 1.
• Lines 22 to 39 define the model.Note that the normal distribution in Stan is parameterized with mean and standard deviation, while the inverse gamma distribution is specified with shape and scale.• Lines 41 to 47 present the generated quantity block that computes the log-likelihood.

•
The code in Listing 5 can be used in different ways, depending on the choice of MCMC sampling function.It can be written within brackets in the R script, saved in a mymodel object, for Stan model compilation with the stan_model() function.
Users can also save it as mymodel.stanfile, and directly call the model from stan() (wrapper) function.We demonstrate the two procedures in Listings 6 and 7.

Run MCMC in Stan
The following code in Listing 6 can be used to run the Stan model specification in R.
Listing 6: Run MCMC in Stan.

•
In Line 4, inside the " " , the code for the model specification that we presented in Listing 5 should be copied as text.Stan also provides a wrapper that allows users to compile and estimate the model in one step, as shown in Listing 7. In this function, it is possible to specify the number of cores ( cores=2 ).This option sets the number of processors for executing the MCMC chains in parallel.It is useful to set as many processors as the hardware and RAM allow (up to the number of chains).Note that the option can be used in the sampling() function shown in Listing 6, but this option is not used in this paper to ensure coherence across the packages.

Run MCMC in NIMBLE
The following code, shown in Listing 9, can be used to run the NIMBLE model specification in R.

•
Lines 6 to 14 define the constants used in the model specification and pos_mu and pos_covmat define the two-dimensional identity matrix.• Lines 18 to 24 provide the initial values for the MCMC estimation.In NIMBLE, log(gamma) is treated as an independent parameter log_gamma .Therefore, an initial value is provided for log_gamma .
• Line 26 defines an array of parameters of interest to be saved in the posterior samples.The parallelization of the MCMC chain sampling can be easily accomplished with a base package parallel in R [13].The code for parallelizing NIMBLE code can be found at https://anonymous.4open.science/r/LSIRM_Estimation-14EE(accessed on 8 May 2023).

MCMC Convergence Diagnostics
There are many options for MCMC diagnostics after fitting a Bayesian model through JAGS, Stan, and NIMBLE within the R environment.For example, JAGS users can directly use jagsUI [3] or rjags [4] packages, while Stan users can either work directly within the rstan package [7] or load the bayesplot package [41] that allows visualizing unique diagnostics permitted by HMC and the NUTS algorithm [42,43].Of note, the recently developed shinystan package [44] is a useful option for analyzing posterior samples.It provides a graphical user interface for interactive plots and tables powered by the Shiny web application framework by R Studio [45].It works with the output of MCMC programs written in any programming language, and it has extended functionality for the rstan package and the No-U-Turn sampler.
Here, we provide the steps for convergence checking through the ggmcmc [46] package, using γ as an example, as shown in Listing 11.The ggmcmc [46] package can be used with any Bayesian software output and provides flexible visualization, bringing the design and implementation of ggplot2 [47] to MCMC diagnostics.
We evaluated posterior samples of the model parameters with some visual diagnostics to assess the convergence of each chain.Figure 1 presents the density plots, traceplots, and autocorrelation plots for γ using the MCMC samples from JAGS, Stan and NIMBLE.The same sets of diagnostics were run for β and σ 2 , and all diagnostics agreed that the model showed good convergence.

Density plots
Trace plots Autocorrelation plots

Model Fit Indices: WAIC and LPPD
In the current paper, we worked with a single model specification; hence model comparisons may be unnecessary.For illustration purposes, here we show how model fit indices can be produced, such as Watanabe-Akaike or Widely Applicable Information Criterion (WAIC) [48] and Log Pointwise Predictive Density (LPPD) [49] (Listing 12).WAIC is reportedly preferred over conventional statistics, such as BIC [50] and DIC [51], because of its desirable property of averaging over the posterior distribution, rather than conditioning on a point [49].An additional advantage is that WAIC produces a stable and positive effective number of parameters ( pwaic ) to evaluate the complexity of the model [52].LPPD estimates the expected log predictive density for each data point directly, providing a comprehensive evaluation of the model fit to the data regardless of the model's complexity [49].These indices are available from rstan and NIMBLE.One can also calculate them using R packages such as loo [53], as used in this paper.Table 1 shows WAIC, pwaic and LPPD from all three packages, which suggests that all these packages produce similar model fit information for the specified LSIRM.
Listing 12: Model fit indices from loo.
# transform the fitted MCMC object to matrix fitted _ samples <-as .matrix ( fitted _ LSIRM ) # extract columns for log likelihood llCols <-grep ( " log _ ll " , colnames ( fitted _ samples ) ) log _ ll .samples <-fitted _ samples [ , llCols ] # calculate the model fit indices library ( loo ) waic _ result <-loo :: waic ( log _ ll .samples ) lppd _ result <-loo :: elpd ( log _ ll .samples ) We evaluated the effective sample size and run time, and, hence, sampling efficiency, of the three packages for LSIRM estimation.All three packages were run on the UCLA hoffman2 cluster computing node that was run on the Intel Xeon Gold 6240 CPU processor, which had a base clock speed of 2.60 GHz and could reach a maximum turbo frequency of 3.90 GHz.The RAM requested for the computing task was fixed at 16 GB.
Table 2 displays the run time for compiling and sampling stages under the chosen iteration/burn-in/thinning conditions of the three packages.All three packages showed satisfactory sampling quality, with Stan performing well with smaller thinning intervals and fewer iterations, while NIMBLE was faster in the sampling stage than the other packages.Note that the compilation time of JAGS included both the model compilation and 1,000 iterations of model adaptation time.As can be seen, parallelization significantly reduced the sampling time for all three packages, compared to the case where a single core was used for sampling the two chains.Note: 'one core' and 'parallel on two cores' indicate the cases where the two chains were run separately and simultaneously (using parallel computing), respectively.The time was based on 15,000 iterations and 10,000 burn-in periods, with a thinning interval of 5 for Stan, and 90,000 iterations and 40,000 burn-in periods, with a thinning interval of 50 for JAGS and NIMBLE The effective sample size (ESS) was a measure of the number of independent draws from the posterior distribution that were equivalent to the actual MCMC draws [54,55].ESS can be computed with the coda package [15] using the posterior samples of parameters as an input, as shown in Listing 13.
Listing 13: Effective sample size from coda.
# transform the fitted MCMC object to matrix fitted _ samples <-as .matrix ( fitted _ LSIRM ) # extract columns for beta betaCols <-grep ( " beta " , colnames ( fitted _ samples ) ) beta .samples <-fitted _ samples [ , betaCols ] # calculate effective sample size for beta library ( coda ) ess _ b <-coda :: effectiveSize ( beta .samples ) ess _ mean _ b <-mean ( ess _ b ) Table 3 lists the ESS of the LSIRM model parameters from the three packages.While we observed comparable ESS, posterior samples from Stan showed the highest overall ESS, particularly for β and γ parameters.We further assessed the efficiency of the MCMC estimation in the three packages by calculating the ratio of the ESS to the run time for sampling (ESS/time in seconds) [26,56].As Table 3 shows, NIMBLE and Stan showed pretty high sampling efficiency.Note that the sample size, model complexity, and choice of priors can significantly impact the run time and sampling efficiency.Users can use the provided tools to evaluate the sampling efficiency for the model of interest with the software of their choice.

Model Parameters
After convergence and model fit checking, posterior samples can be extracted and further analyzed using various options for post-processing packages in R.Here we demonstrate a universal procedure with coda to summarize the parameters of interest.Listing 14 presents the code for this step.
Listing 14: Extracting information from MCMC chains.

Interaction Map Visualization
For interaction map visualization, the posterior samples of the person and item latent positions (z and w) needed to be post-processed to resolve the unidentifiability and match the samples across MCMC iterations, as discussed in Section 2.3.We applied Procrustes matching, as shown in Listing 15.One can use the matched samples to draw the interaction map.Note that zz.samples is a list object that has the dimensions of the number of chains ×, the length of the posterior samples ×, the number of coordinates of z.The number of coordinates equals the product of the sample size (N) and the size of the latent space dimension (p = 2).JAGS and NIMBLE sort the last dimension of the zz.samples as z 11 , . . ., z N1 , z 12 , . . ., z N2 , which can be directly processed with the code provided in Listing 15.However, Stan sorts the estimated coordinates as z 11 , z 12 , . . ., z N1 , z N2 ; thus, this should be reordered before proceeding with the provided function.The same logic applies to ww.samples .The functions used in Listing 15 can be directly retrieved using the source function with the link provided in Line 13, and can also be obtained from https://anonymous.4open.science/r/LSIRM_Estimation-14EE(accessed on 8 May 2023).
Listing 15: Procrustes matching and visualization of interaction map and strength (inverse distance) plots.Figure 2 displays the estimated interaction maps from the three packages.In all plots, solid black dots represent respondents, and colored numbers represent items, where Cursing items are in red, Scolding in green, and Shouting in blue.The exact positions of items and respondents were not identical across the three packages, due to the latent positions' unidenfiability discussed earlier.That said, it can be observed that the overall configuration of the item positions was similar across the three plots.For example, in all plots, the blue item groups were distant from the red and green item groups.There was some separation between the red and green items, although, overall, the two groups (red and green) were close to each other.This indicated strong dependencies (or similarities) among the items within each item group.There were some similarities between Cursing and Scolding items (red and green), which were stronger than similarities between Cursing and Shouting items (red and blue).As explained earlier, a large distance to an item indicates that the respondent has a lower success probability for the item given her/his overall trait level and the item's overall difficulty.In other words, a respondent's distance to an item can be understood as her/his weakness (or lower likelihood) to the item, given her/his overall trait level and the item's difficulty level.Thus, it can be informative to quantify a person's distances to the test items to evaluate which items the respondent shows particular strengths or weaknesses (or likelihood or unlikelihood) towards.Figure 3 illustrates Respondent 1's inverse distances, and, therefore, indicates her/his likelihood of endorsing the 24 verbal aggression items, given the respondent's overall latent trait level and the item's difficulty.The items are ordered from the largest to the shortest inverse distances.The three item groups are colored the same way as in Figure 2. Overall, the results were similar across the three packages.To be specific, the three packages gave the same order for the top 5 items.While there were some minor differences, the rank correlations were nearly 1.0 between any package pair, confirming strong similarity in the person-item distances across the three packages.

Discussion
In this paper, we presented a practical guide to Bayesian estimation of latent space item response models (LSIRMs), a recent development in IRT, using three popular, free of charge, Bayesian packages, namely, JAGS, Stan, and NIMBLE, run from R through specific packages.With an empirical example, we provided detailed instructions for specifying the model and running the code with each package.We also evaluated model convergence, model fit, sampling quality, and run time.Additionally, we presented the visualization of interaction maps with the three packages.
In regard to the three packages, JAGS is relatively easy to learn and implement, making it suitable for researchers who may be new to Bayesian modeling.With a relatively long history, JAGS boasts an extensive user community that can provide a wide range of support as open-source software.JAGS also supports a wide range of probability distributions, making it a convenient choice for diverse modeling scenarios.Stan, on the other hand, offers an advanced and flexible modeling framework, allowing for efficient estimation of highly complex models involving high-dimensional model parameter space.One key advantage of Stan is its built-in support for parallel computation.NIMBLE strikes a balance between the ease of use of JAGS and the flexibility of Stan.Using a syntax similar to BUGS, NIMBLE allows users to write their own algorithms with the nimbleFunction system, which can be helpful for research and experimentation.In the context of LSIRM estimation, we found that the results were consistent across the three packages with the empirical data under investigation.NIMBLE was somewhat faster than JAGS and Stan, but NIMBLE needed longer chains with larger thinning intervals than Stan.Thus, one may choose a package depending on specific needs and expertise.It is worth noting that, in our illustrative examples, the number of iterations, burning periods, and the length of thinning intervals were chosen to yield posterior samples that had comparable qualities between the three chosen packages.Longer chains were needed for JAGS and NIMBLE to match with Stan's results in our example.However, different set ups, potentially with a shorter chain, could be considered when comparability between the packages is not of concern.
In closing, we demonstrated, in the current paper, the estimation of LSIRM in an original specification for binary item responses.As LSIRM was introduced as an extension of the Rasch model for conditional dependencies, one may be interested in model selection between LSIRM and the Rasch model.The original paper, [1], addressed the model selection question by applying a spike-and-slab prior to the weight parameter γ, which can also be implemented in the discussed Bayesian packages.Lastly, to further broaden the applicability of LSIRM, additional research would be needed to show the LSIRM estimation in various other specifications, such as for polytomous data, which is possible due to the flexibility of the Bayesian packages discussed in this paper.

• From Lines 6
to 15 data and hyperparameters are set • In Lines 19 to 24 the two chains of MCMC are initialized with the same set of values.• In Lines 30 to 38 the MCMC algorithm is run to sample from the posterior distribution of the model parameters, given the data.The samples are stored in an object of class stanfit .

Listing 10 :
A Wrapper function for NIMBLE.fitted _ LSIRM _ w <-nimble :: nimbleMCMC ( code = mymodel , # code specifies the model data = mydata , # list of data constants = myconst , # list of constants used in the model inits = myinits , # list of initial values monitors = myparams , # an array of parameters of interest to be saved niter = 90000 , # number of total iterations nburnin = 40000 , # number of burn -in periods thin = 50 , # size of thinning interval nchains = 2 , # number of Markov chains progressBar = TRUE , # show the progress while sampling s a m p l e s A s C o d a M C M C = T ) # save the samples as an mcmc .list object

Figure 1 .
Figure 1.Visual diagnostics of γ from the three packages.The first, second, and third rows show density plots, trace plots, and autocorrelation plots, respectively.The three columns indicate JAGS, Stan, and NIMBLE results, respectively.Note that the first MCMC chain is displayed in red and the second MCMC chain in light blue.

Figure 2 .
Figure 2.Estimated interaction maps of the matched latent positions from the three packages.In all plots, solid black dots represent respondents, and colored numbers represent items, where Cursing items are in red, Scolding in green, and Shouting in blue.

Figure 3 .
Figure 3. Respondent 1's inverse distances to the 24 verbal aggression items estimated from the three packages.The height of the bars indicates respondent 1's likelihood of endorsing the item, given her/his overall trait level and the item's difficulty.The items are ordered from the highest to lowest endorsement likelihood for Respondent 1.The coloring of the bars is based on item types, where Cursing items are in red, Scolding in green and Shouting in blue.

Listing 4 :
A Wrapper function for JAGS.
Listing 5 presents the LSIRM model specification written in Stan language.Listing 5: LSIRM model specification in Stan.
Lines 38 to 44 run MCMC, and save the posterior samples of the parameters specified in myparams vector.NIMBLE also provides a wrapper that allows users to combine Lines 27 to 44 in one function as shown in Listing 10.
• In Lines 28 to 31, the nimbleModel() function processes the model code, constants, data, and initial values and returns a NIMBLE model.NIMBLE checks the model specification, model sizes, and dimensions at this step.• Line 32 compileNimble() compiles the specified model.• Lines 33 to 35 create and return an uncompiled executable MCMC function.• Line 36 compiles the executable MCMC function regarding the complied specified model.•

Table 1 .
WAIC and LPPD for LSIRM model from the three packages.

Table 2 .
The compilation time and sampling time of the three packages.

Table 3 .
The ESS and sampling efficiency, ESS/time (s), of the three packages.

Table 4 .
The posterior means (Post.m)and 95% credible intervals (CI) of the model parameters of the three packages.