Bayesian Variable Selection in Dynamic Item Response Theory Models

Abstract

The recent surge in computerized testing brings challenges in the analysis of testing data with classic item response theory (IRT) models. To handle individually varying and irregularly spaced longitudinal dichotomous responses, we adopt a dynamic IRT model framework and then extend the model to link with individual characteristics at a hierarchical level. Further, we have developed an algorithm to select important characteristics of individuals that can capture the growth changes of one’s ability under this multi-level dynamic IRT model, where we can compute the Bayes factor of the proposed model including different covariates using a single Markov chain Monte Carlo output from the full model. In addition, we have shown the model selection consistency under the modified Zellner–Siow prior, and we have conducted simulations to illustrate the properties of the model selection consistency in finite samples. Finally, we have applied our proposed model and computational algorithms to a real data application, called EdSphere dataset, in educational testing.

Keywords

item response theory (IRT) model Bayes factor (BF)Markov chain Monte Carlo (MCMC)model consistency variable selection longitudinal testing data

1. Introduction

Item response theory (IRT) models play a critical role in the design and analysis of educational testing. The classic IRT models connect the correctness of the answer (a dichotomous random variable) for each item with the ability, the item difficulty, and potentially other item parameters via different types of link functions. The widely used Rasch model (Rasch, 1960) is a one-parameter IRT model with a logistic link. The two-parameter logistic IRT model (Birnbaum, 1968) adding a discrimination parameter to the Rasch model allows different scaling of items. Further, the two-parameter IRT model can be extended to a three-parameter IRT model (Lord, 1980), by involving an additional guessing parameter. More recent applications of the IRT models can be found in Van der Linden (2018) and Irwing et al. (2018) and references therein.

There are other extensions and variations of the IRT models. For example, multi-level IRT models (Bock & Mislevy, 1989; Fox, 2005) introduce extra layers for the ability parameter; the $β^{3}$ IRT model (Y. Chen et al., 2019) accommodates continuous responses and can generate a much-enriched family of item characteristic curves; and more recently, there is literature on parsimonious asymmetric IRT models (Shim et al., 2022) and count thresholds models (Tutz, 2022), where the latter one provides a flexible IRT model for count data without assuming a fixed distribution for the responses.

1.1 EdSphere Dataset and Dynamic IRT Model

A key assumption in the classic IRT models is the local independence assumption, which assumes the responses are independent given the item parameters and the individual latent ability. However, this assumption may not be true for our study, the EdSphere dataset. The EdSphere dataset is provided by Highroad Learning Company and it is obtained from a personalized literacy learning platform, which collects student performance and strategic behaviors each time when he/she reads an article. In this learning platform for a typical reading comprehensive test, a student can select a reading article from a pool of articles having the test difficulty in a range targeted to the current estimate of the student’s ability. Then, once the article is selected, the computer program will randomly select a sample of the eligible words to be clozed (i.e., removed and replaced by a blank) in the article and then displays the article to the student with these words clozed. The clozed words produced by this procedure are randomized and single time used items, which means if another student happens to select the same article to read, the cohort of clozed words presented to them are different. As a consequence, the replication of clozed words (i.e., items in IRT models) among students is highly improbable, so obtaining empirical estimates of item parameters is not feasible. Hence, in our proposed model, the item difficulty of EdSphere reading test is assumed to follow a measurement error model, where the ensemble mean is the test difficulty of the article estimated using the proprietary data and following method of moments techniques (c.f., Swartz et al., 2016 for more details).

Currently, the EdSphere dataset consists of 16,949 students from a school district in Mississippi who participated over 5 years (2007–2011) in the EdSphere learning platform. The students were in different grades and could enter and leave the program at different times. They were free to take tests on different days and had different time lapses between tests. An illustration of this kind of data structure is given in Table 1.

Table 1.

Illustration of the Unequally Spaced Response Data Structure

Date	Day 1	Day 2	Day 3	…	Day M	…	Day N	…
Examinee A	1 test	N/A	3 tests	…	5 tests	…	N/A	…
Examinee B	N/A	Enter	1 tests	…	N/A	…	4 tests	…
Examinee C	4 tests	2 tests	1 tests	…	N/A	…	Leave	…

This learning platform continuously yields longitudinal dichotomous observations located at individually varying and irregularly spaced time points, which suggests that we need to model the changes of latent traits using a dynamic structure. Further, according to the test design in EdSphere, there are several existing factors, such as overall comprehension, emotional status, and others (Wang et al., 2013), which would undermine the local independence assumption of classic IRT models. Therefore, we need some new IRT models that can accommodate the modern computerized testing (not merely EdSphere datasets) with the described features above, that is, randomized items, longitudinal observations, and local dependence.

To integrate these features and model the growth of latent traits in computerized testing, Wang et al. (2013) proposed a dynamic IRT model by combining the ideas of parametric functions of time as well as Markov chain models to capture the trajectory of latent traits. They embedded IRT models into a new class of state space models for analyzing longitudinal data located at individually varying and irregularly spaced time points. The trend of the ability growth has been governed by an individually varying parameter, $c_{i}$ , in the system equation of their dynamic IRT model, which can be viewed as a “growth factor.” To better explain this “growth factor” or advise the educators, we are very interested in finding out the important characteristics of each person, such as gender, grade, initial ability, and so on, that can help ability growth. Thus, based on the dynamic IRT model (Wang et al., 2013), we build up another hierarchical layer modeling the relationship between the “growth factor” and interesting covariates in the study. We then describe the Bayesian inference procedure to select the important covariates for this proposed dynamic IRT model (see Section 2.2).

1.2 Variable Selection With Bayes Factor

Finding variables associated with the “growth factor” in the proposed model is a natural and important question to ask, which can help better understand the features that encourage students’ ability growth. In educational testings, variable selection has been widely considered as it may influence the estimation and predictive accuracy of the model performance (Hamza et al., 2022; Muttakin et al., 2021; Rahman et al., 2017).

For fully Bayesian analysis of variable selection, there are some common techniques often used in literature, such as the deviance information criterion (DIC; Spiegelhalter et al., 2002), the Bayesian lasso (Park & Casella, 2008), the Bayesian spike-and-slab lasso (Ročková & George, 2018), the variational Bayes dynamic variable selection algorithm (Koop & Korobilis, 2018), the Bayes factor (BF) approach (Berger & Pericchi, 1996; Pericchi, 2005) and more recently, the variable selection using reversible jump based on simulating piecewise deterministic Markov processes samplers (Chevallier et al., 2022). Many of these criteria are initially discussed within a linear model framework, and they may be applied to hierarchical models such as the dynamic IRT model under consideration. In this paper, we will use the BF approach for variable selection in our proposed dynamic IRT model. It is not only because BF is a fundamental criterion for model comparison under Bayesian inference, which could be viewed as the Bayesian equivalent of the likelihood ratio test, but also because there is a steady increase in attention to the BF approach as a tool for hypothesis evaluation and model selection in psychological research (Heck et al., 2022; Pan & Yin, 2017; Tijmstra & Bolsinova, 2019).

In addition, the BF has been shown to enjoy good frequentist properties for variable selection in several different types of models. For example, Fernandez et al. (2001) and Liang et al. (2008) discussed some attractive model selection properties for BF under the linear model framework when $g$ -prior (Zellner, 1986) is assigned for the regression coefficients. In particular, Liang et al. (2008) showed that when the Zellner–Siow prior (Zellner & Siow, 1980), that is, the Inv-Gamma $(1 / 2, n / 2)$ prior is specified for the additional parameter $g$ in the $g$ -prior, the BF approach for variable selection has posterior model selection consistency. Wu et al. (2016) and Li and Clyde (2018) further discussed the model selection properties for the mixtures of $g$ -priors under the generalized linear model framework. Recently, Castellanos et al. (2021) developed a model selection with mixtures of $g$ -priors for censored data. These motivate us to consider utilizing the $g$ -priors for the variable selection problem in modeling the “growth factor” for the proposed dynamic IRT model. However, it is unclear whether or not some of the attractive model selection properties for the mixture of $g$ -priors especially model selection consistency will still hold under our proposed hierarchical models. Thus, we aim to be the first to establish the model selection consistency for the mixture of $g$ -priors in the hierarchical models.

However, due to the complex structure of our proposed dynamic IRT model and the high dimensionality of the parameter space, there are no analytical forms available for the BFs of the dynamic IRT models (while noticing there are analytical forms available for the linear regression model discussed in Liang et al., 2008 when $g$ is fixed). To solve this computational challenge, inspired by M.-H. Chen (2005), we develop an algorithm to compute the BFs under the dynamic IRT model using a single MCMC output, which is a computationally feasible approach if the number of covariates is moderate or even when the number of covariates is relatively large but still less than the sample size.

The rest of the paper is organized as follows. In Section 2, we introduce the proposed dynamic IRT model and Bayesian analysis techniques. In Section 3, we propose the method to compute BF with a single MCMC output. Section 4 establishes the model selection consistency for the proposed dynamic IRT model under the mixture of $g$ -priors. Section 5 contains the simulation study, and in Section 6, we provide the results for the real data analysis. Finally, we summarize the results in Section 7 and point out some future directions of the research.

2. The Proposed Dynamic IRT Model and Bayesian Inference

In this section, we follow the general framework of the dynamic IRT model proposed by Wang et al. (2013) to model the EdSphere dataset but with some changes, where the major change is to introduce an additional level for modeling the “growth factor.” There are several reasons to adopt the general framework of the dynamic IRT model proposed by Wang et al. (2013). The EdSphere dataset is generated from an online learning platform. Thus, the parametric model assumptions we make for the first and second levels come from how the real data was generated from the platform. To be specific, first, the EdSphere test design does not include any item discrimination parameter, and therefore, a one-parameter IRT model is more reasonable in the study of EdSphere datasets. Second, the EdSphere learning platform is in fact a Computer Adaptive Instruction and Testing (CAIT; c.f. Wang et al., 2013 and Swartz et al., 2016, 2011). With CAIT, a test pool of articles is selected for the student based on an estimate of his/her current ability. After the student selects an article from this test pool, the test questions (i.e., items) are then generated before the reading commences according to some prespecified computer protocols. In the current learning platform design, the computational algorithms for the CAIT to estimate his/her current ability are through an autoregressive (AR) model of order 1 to update his/her current ability. Thus, the assumption of an AR(1) model to capture the ability changes would be in alignment with the EdSphere platform design. Hence, based on these two reasons for the EdSphere learning platform design, we specify our proposed model with the one-parameter IRT model in the first level and the AR(1) model in the second level for capturing one’s ability growth, similar to the dynamic IRT models proposed by Wang et al. (2013).

The rest of this section is structured as follows. First, we will introduce our hierarchical model framework, and then discuss the prior choices for unknowns and the corresponding Bayesian computation schemes for our proposed model. Notice that although the discussion of the first and second levels of our proposed model in this section is specifically refers to one-parameter IRT models and AR(1) models, the computational algorithms of BF discussed in Section 3 and the property of the model selection consistency discussed in Section 4 are more general derivations and conclusions for any type of hierarchical models if they can satisfy the mild conditions mentioned in our theorems later. We have provided more details on this point in Section 7.

2.1 Model Settings and Prior Specifications

To analyze dichotomous responses collected at individually varying and irregularly spaced time points, we propose a three-level hierarchical IRT model to capture the changes in ability and associate the changes with certain characteristics of a person. In the first level of our model, the correctness of the items is linked to the latent ability of each individual, the difficulty of the items, and certain random effects. In the second level, it is assumed that the latent ability of the current day depends on the latent ability of the previous day and an increment term related to the time-lapse with uncertainty. These two levels are similar to the dynamic IRT models proposed by Wang et al. (2013). In the third level, we model the relationship of the individual “growth factor” $c_{i}$ s with the characteristics of test takers, such as socio-economic status and gender, in a log-normal linear regression. To summarize, the proposed dynamic IRT model is as below:

Level 1 : \Pr (X_{i, t, s, l} = 1 ∣ θ_{i, t}, η_{i, t, s}, d_{i, t, s, l}) = Φ (θ_{i, t} - d_{i, t, s, l} + η_{i, t, s}),

(1)

Level 2 : θ_{i, t} = θ_{i, t - 1} + c_{i} (1 - ρ θ_{i, t - 1}) Δ_{i, t}^{+} + ω_{i, t},

(2)

Level 3 : \log (c_{i}) = α + z_{i}^{T} β + ζ_{i},

(3)

where $X_{i, t, s, l}$ represents the correctness of the $l$ -th item in the $s$ -th test on the $t$ -th day taken by the $i$ th individual, for $i = 1, 2, \dots, n$ , $t = 1, 2, \dots, T_{i}$ , $s = 1, 2, \dots, S_{i, t}$ and $l = 1, 2, \dots, K_{i, t, s}$ ; $θ_{i, t}$ denotes the latent ability of individual $i$ on the test day $t$ (noting when $t = 1$ in Equation 2, $θ_{i, t - 1} = θ_{i, 0}$ on the right side, which is the initial ability of the $i$ th individual when she/he enters the study); $η_{i, t, s}$ is the test-wise random effect, and we assume $η_{i, t} = (η_{i, t, 1}, \dots, η_{i, t, S_{i, t}})' ~ N (0, τ_{i}^{- 1} I)$ for $i = 1, \dots, n$ , $t = 1, \dots, T_{i}$ (noting in the original dynamic IRT models proposed by Wang et al., 2013, they further decomposed the random effect term into daily effects and test effects, but here we follow the modeling of Observation Equation in Y. Liu and Wang, 2020 for EdSphere datasets and do not further divide the random resources since it is not our main purpose in this article); and since often there is no replication of the test items, we model the item difficulty using a measurement error model, that is, $d_{i, t, s, l} ~ N (a_{i, t, s}, σ^{2})$ , where the same assumption on $d_{i, t, s, l}$ was discussed in Wang et al. (2013); Swartz et al. (2016); and Y. Liu and Wang (2020). In this measurement error model, $a_{i, t, s}$ is called an ensemble of item difficulties associated with the article, which is given by its theoretical text complexity. The study of Stenner (2010) provides an empirical verification for the legitimate of using this theoretical value for the ensemble mean $a_{i, t, s}$ . Further, $σ$ is assumed to be known due to the test design of the EdSphere learning platform. In Level 1, similar to Y. Liu and Wang (2020), we use the Probit link $Φ^{- 1} (\cdot)$ , which represents the inverse cumulative distribution function of the standard normal distribution, to model the dichotomous data as it is more convenient in Bayesian computation.

In Level 2, $ρ$ is a parameter that reflects the “plateau effects” of the ability (i.e., the ability becomes much harder to grow when it approaches to mature), which can be used to slow down the linear growth of the ability when the ability becomes larger. Following the discussion of Swartz et al. (2011), we assume $ρ$ to be known, which is an empirical estimation from previous study and Wang et al. (2013) has used the same assumption for $ρ$ in their analysis. $Δ_{i, t}$ is the time lapse between the $t$ -th test day and the $(t - 1)$ -th test day for the individual $i$ , $t = 1, \dots, T_{i}$ . The term $Δ_{i, t}^{+} = \min (Δ_{T_{\max}}, Δ_{i, t})$ is the time lapse truncated by a pre-specified maximum time interval, and it reflects the fact the student’s ability may not be growing proportional to the time lapse when the gap is over certain number of days, that is, $Δ_{T_{\max}}$ here. In practice, it may occur when students take a vacation. We follow Wang et al. (2013) and set $Δ_{T_{\max}}$ = 14 days. Further, $c_{i}$ in Level 2 can be considered as the average growth rate for the $i$ -th person’s ability over time, which can be interpreted as the individual “growth factor” controlling the overall growth of one’s ability trajectory. Finally, $ω_{i, t}$ is used to capture any uncertainty that could not be fully explained by the first two terms of modeling $θ_{i, t}$ and we assume $ω_{i, t} ~ N (0, ϕ^{- 1} Δ_{i, t})$ . In fact, Level 2 can be rewritten as $λ_{i, t} = g_{i, t} λ_{i, t - 1} + ω_{i, t}$ with $λ_{i, t} = θ_{i, t} - ρ^{- 1}$ and $g_{i, t} = 1 - c_{i} {ρ Δ}_{i, t}^{+}$ , which is similar to the structure of an AR model with order 1 and $g_{i, t}$ s serve the role of parameters in this representation. In $g_{i, t}$ s, only $c_{i}$ s are unknown parameters and related to individuals. Thus, we consider modeling the relationship of these covariates with respect to $c_{i}$ s. Other parameters related to personal indexes, such as the precisions of random effects $η_{i, t}$ s, that is, $τ_{i}$ s, are less likely to be dependent on some features of covariates since $η_{i, t}$ are a vector of random effects to capture the randomness of test and daily effects of individuals that could not be explained by other factors. Thus, $c_{i}$ s are the most reasonable parameters in our study to investigate their relationship with other personal covariates.

In Level 3, we use $α$ to denote the intercept, $z_{i} = (z_{i 1}, \dots, z_{ip})'$ to indicate a $p$ -dimensional covariate vector for the $i$ -th individual and $β = (β_{1}, \dots, β_{p})'$ to be the $p$ -dimensional vector of coefficients for the covariates. Further, we define $Z = {(z_{1}, \dots, z_{n})}^{T}$ , which is a $n \times p$ design matrix of the covariates, and we assume each column has been centered, that is, $1^{T} Z = 0$ , where $1$ is a unit vector. Besides, $ζ_{i}$ is the random error, which is assumed to follow $ζ_{i} ~^{i . i . d .} N (0, ψ^{- 1})$ .

In the proposed three-level dynamic IRT models, we aim to perform variable selection for the covariates at Level 3, that is, to find all non-zero $β_{j}$ s in $β$ . Liang et al. (2008) have showed some very attractive variable selection properties for the BF with mixtures of $g$ -priors in a linear model setting, thus it motivates us to apply the BF approach with mixtures of $g$ -priors for the variable selection in our hierarchical model and also to investigate whether the corresponding model selection properties for the mixtures of $g$ -priors will still hold for the hierarchical model settings in general.

Let us define a $p$ -dimensional vector of indicators $γ = {(γ_{1}, \dots, γ_{p})}^{T}$ with $γ_{j} = I (β_{j} \neq 0)$ , for $j = 1, \dots, p$ , and $I (\cdot)$ is an indicator function. Then, we can use $β_{γ}$ to present a $p_{γ}$ -dimensional vector of nonzero regression coefficients and $Z_{γ}$ to index the corresponding $n \times p_{γ}$ design matrix under the model $M_{γ}$ . For the precision parameter $ψ$ and the intercept $α$ in Level 3, we could not apply the non-informative priors used by Liang et al. (2008) anymore. It is because those non-informative priors could lead to an improper joint posterior distribution of unknown parameters for the proposed dynamic IRT model (the proof details are shown in the Supplemental File in the online version of the journal). Instead, we assign the following proper priors for $α$ , $β_{γ}$ , and $ψ$ of Level 3 for the proposed dynamic IRT model:

π (α ∣ ψ, M_{γ}) = {(\frac{2 π}{λ ψ})}^{- \frac{1}{2}} \exp {- \frac{1}{2} λ ψ α^{2}},

(4)

[β_{γ} ∣ ψ, M_{γ}] ~ N (0, \frac{g}{ψ} {(Z_{γ}^{T} Z_{γ})}^{- 1}),

(5)

π (ψ ∣ M_{γ}) = \frac{b^{a}}{Γ (a)} ψ^{a - 1} \exp {- b ψ} .

(6)

Thus, theorems such as model selection consistency proved by Liang et al. (2008) cannot be automatically extended in our model setting since the proper priors has to be used here.

Further, Liang et al. (2008) have pointed out the marginal likelihoods of all models have closed-form expressions when $g$ is fixed in their linear model settings, and the BF can also be computed via Laplace approximation when the Zellner–Siow prior is used, that is, $g ~ Inv - Gamma (1 / 2, n / 2)$ . However, under our hierarchical model, there are no analytical forms available for neither the marginal likelihoods nor BFs. In addition, it is not feasible to use Laplace approximation to compute the BF. Hence, a tractable computational algorithm is needed to obtain the BFs for our dynamic IRT models.

2.2 Posterior Distribution and Computation

To complete our Bayesian analysis, we further specify the priors for the remaining model parameters as follows: assign (a) $τ_{i} ~ Gamma (a_{τ}, b_{τ})$ with $a_{τ}$ and $b_{τ}$ being the shape and rate parameter, respectively, $i = 1, 2, \dots, n$ ; (b) $ϕ ~ Gamma (a_{ϕ}, b_{ϕ})$ with $a_{ϕ}$ and $b_{ϕ}$ as the shape and rate parameter, separately; and (c) $g ~ Inv - Gamma (1 / 2, n / 2)$ . Define $c = (c_{1}, \dots, c_{n})'$ , $θ = (θ_{1}, \dots, θ_{n})'$ with $θ_{i} = (θ_{i, 1}, \dots, θ_{i, t})'$ , $η = (η_{1}, \dots, η_{n})'$ with $η_{i} = (η_{i, 1}, \dots, η_{i, t})'$ and $τ = (τ_{1}, \dots, τ_{n})'$ . Also, we assign a prior for one’s initial ability $θ_{i, 0}$ , which is a normal distribution, with mean centered at $μ_{G_{j_{i}}}$ and variance being $V_{G_{j_{i}}}$ . Notice the index $j_{i}$ indicates the grade ( $j$ ) that the $i$ th individual belongs to. Then, the joint posterior distribution of $θ, c, τ, η, ϕ, β_{γ}, α, g$ , and $ψ$ given the data $X = {X_{i, t, s, l}}$ and covariates $Z_{γ}$ is

\begin{matrix} π (θ, c, τ, η, ϕ, β_{γ}, α, g, ψ ∣ X, Z_{γ}) \\ \propto Π_{i = 1}^{n} \frac{1}{\sqrt{2 π} V_{G_{j_{i}}}} \exp (- \frac{{(θ_{i, 0} - μ_{G_{j_{i}}})}^{2}}{2 V_{G_{j_{i}}}}) {Π_{i = 1}^{n} Π_{t = 1}^{T_{i}} Π_{s = 1}^{S_{i, t}} Π_{l = 1}^{K_{i, t, s}} \frac{1}{\sqrt{2 π σ}} \exp (- \frac{ε_{i, t, s, l}^{2}}{2 σ^{2}})} \\ \times {Π_{i = 1}^{n} Π_{t = 1}^{T_{i}} {(\frac{τ_{i}}{2 π})}^{(S_{i, t} - 1) / 2} \exp (- \frac{τ_{i} η_{i, t}^{T} η_{i, t}}{2})} \\ \times {Π_{i = 1}^{n} Π_{t = 1}^{T_{i}} Π_{s = 1}^{S_{i, t}} Π_{l = 1}^{K_{i, t, s}} Φ [{(- 1)}^{X_{i, t, s, l}} (θ_{i, t} - d_{i, t, s, l} + η_{i, t, s})]} \\ \times {Π_{i = 1}^{n} Π_{t = 1}^{T_{i}} \sqrt{\frac{ϕ}{2 π Δ_{i, t}}} \exp (- \frac{ϕ {θ_{i, t} - θ_{i, t - 1} - c_{i} (1 - ρ θ_{i, t - 1}) Δ_{i, t}^{+}}^{2}}{2 Δ_{i, t}})} \\ \times Π_{i = 1}^{n} \frac{\sqrt{ψ}}{c_{i}} \exp {- \frac{ψ {[(c_{i}) - (α + z_{i}^{T} β_{γ})]}^{2}}{2}} ψ^{- \frac{1}{2}} \exp {- \frac{1}{2} λ ψ α^{2}} ψ^{a - 1} \exp {- b ψ} \\ \times g^{- p / 2} \exp {- \frac{g^{- 1} ψ β_{γ}^{T} (Z_{γ}^{T} Z_{γ}) β_{γ}}{2}} g^{- 3 / 2} \exp (- \frac{n}{2 g}) Π_{i = 1}^{n} π (τ_{i}) π (ϕ), \end{matrix}

(7)

where $π (τ_{i})$ and $π (ϕ)$ are the gamma priors assigned at the beginning of this subsection. For the simulation and application studies in this article, we choose the hyperparameters in the priors of $τ_{i}$ and $ϕ$ with $a_{τ} = 0.01$ , $b_{τ} = 0.01$ and $a_{ϕ} = 1$ , $b_{ϕ} = 0.001$ , respectively.

Since the joint posterior (7) is not available in the closed form, we develop an MCMC sampling scheme to draw samples from the posterior (7). The sampling procedure of Level 1 and Level 2 are similar to Appendix A of Wang et al. (2013). Thus accordingly, we introduce a latent variable $Y_{i, t, s, l}$ using the data augmentation idea from Albert and Chib (1993), that is, letting $Y_{i, t, s, l} ~ N_{+} (θ_{i, t} - d_{i, t, s, l} + η_{i, t, s})$ if $X_{i, t, s, l} = 1$ and $Y_{i, t, s, l} ~ N_{-} (θ_{i, t} - d_{i, t, s, l} + η_{i, t, s})$ if $X_{i, t, s, l} = 0$ , where $N_{+} (\cdot, \cdot)$ indicates a normal distribution truncated by zero at the left side and $N_{-} (\cdot, \cdot)$ is a normal distribution truncated by zero on the right side. The sampling details of the proposed dynamic IRT model are in Supplemental Appendix A (available in the online version of this article).

3. Computing BFs With a Single MCMC Output

In our Bayesian analysis of the proposed IRT models, we can compare different choices of variables in Level 3 using BFs. But to compute BFs, we need to first calculate marginal likelihoods, which usually in practice, are not analytically available. Often, we can numerically calculate them using Monte Carlo methods, such as the Harmonic Mean method (Newton & Raftery, 1994), the Chib method (Chib, 1995), the Stepping Stone method (Xie et al., 2011), just naming a few. Recently, Y. Liu et al. (2019) have given a thorough comparison of the most commonly used Monte Carlo methods for computing the marginal likelihood of IRT models.

However, due to the complex structure and high dimensionality of the parameter space in the proposed dynamic IRT model, computing the marginal likelihood for each corresponding model in the variable selection with existing Monte Carlo methods is very challenging. In the paper of Verdinelli and Wasserman (1995), they presented the generalized version of Savage–Dickey (SD) density ratio to compute BF by multiplying a correction term with the SD density from Dickey and Lientz (1970). In M.-H. Chen (2005), they further developed an approach to obtain the BF of each model versus the full model using only a single MCMC output yielding from the full model based on the extension of the generalized SD ratio from Verdinelli and Wasserman (1995).

Motivated by the idea from M.-H. Chen (2005) to compute BF efficiently in a single MCMC output, we formally define notations of the models under comparison as follows. Let $M_{N}$ be the null model where all $β_{j}$ coefficients are $0$ , that is, $β = 0$ ; $M_{F}$ denotes the full model where all $β_{j}$ coefficients are considered to be non-zero; let $M_{γ}$ denote the model with $β_{- γ} = 0$ , where ${(β_{γ}, β_{- γ})}^{T} = β$ . Let $p_{γ}$ be the length of $β_{γ}$ , thus we have $p_{γ} + p_{- γ} = p$ for any $γ$ . Notice for the variable selection problem in a hierarchical model, the setups for Level 1 and Level 2 are the same for every model in the comparison.

We will use priors specified in Section 2 for all unknown parameters and apply a uniform prior on the model space following Liang et al. (2008). Based on the current linear structure at Level 3, we have the identity

π (c ∣ α, β_{γ}, ψ, M_{γ}) = π (c ∣ α, β_{γ}, β_{- γ} = 0, ψ, M_{F})

(8)

holds, which leads to the proposition below.

Proposition 1. For $M_{γ} \neq M_{N}$ , the BF of model $M_{γ}$ to the full model $M_{F}$ under the proposed dynamic IRT model framework is

B F_{H} [M_{γ} : M_{F}] = E [\frac{π (β_{γ} ∣ ψ, g, M_{γ})}{π (β_{γ}, β_{- γ} = 0 ∣ ψ, g, M_{F})} π (β_{- γ} = 0 ∣ c, α, β_{γ}, ψ, g, X, M_{F})],

(9)

and the BF of the null model $M_{N}$ to the full model $M_{F}$ is

B F_{H} [M_{N} : M_{F}] = E [\frac{π (β = 0 ∣ c, α, ψ, g, X, M_{F})}{π (β = 0 ∣ ψ, g, M_{F}) π (g)}],

(10)

where both expectations (9) and (10) are taken w.r.t. $π (c, α, β_{γ}, ψ, g ∣ X, M_{F})$ . The $π (c, α, β_{γ}, ψ, g ∣ X, M_{F})$ is the joint posterior distribution of $c, α, β_{γ}, ψ, g$ under the full model $M_{F}$ .

We use the $B F_{H}$ symbol above to emphasize the BF is defined under the hierarchical model framework. The detailed proof for Proposition 1 is postponed to Supplemental Appendix B (available in the online version of this article). Based on Proposition 1, we can compute the BF of any model $M_{γ}$ to the full model $M_{F}$ using a single MCMC output sampling from the joint posterior in Equation 7. Then, the BF of an arbitrary pair of models ( $M_{γ'}$ , $M_{γ}$ ) is obtained by

B F_{H} [M_{γ'} : M_{γ}] = \frac{B F_{H} [M_{γ'} : M_{F}]}{B F_{H} [M_{γ} : M_{F}]} .

(11)

Notice that Equations (9) and (11) in Proposition 1 are in a similar structure as what was presented in Theorem 2 of M.-H. Chen (2005), but we have shown that it is applicable in the hierarchical model framework. Further, unlike the generalized version of the SD density ratio in Verdinelli and Wasserman (1995), both expectations (9) and (10) are taken w.r.t. the posterior of $M_{F}$ , which facilitates us to use a single MCMC output from the joint posterior under the full model $M_{F}$ to compute the ratios of BF. This is an essential step to make the computation of BF feasible for complex hierarchical models as the one presented in the paper.

Next, using the assigned prior distribution for $β$ in Equation (5), we can simplify the quantity inside the expectation of Equation (9) as below

\begin{matrix} \frac{π (β_{γ} ∣ ψ, g, M_{γ})}{π (β_{γ}, β_{- γ} = 0 ∣ ψ, g, M_{F})} π (β_{- γ} = 0 ∣ c, α, β_{γ}, ψ, g, X, M_{F}) \\ = \frac{{(2 π)}^{- p_{γ}} \det {g ψ^{- 1} {(Z_{γ}^{T} Z_{γ})}^{- 1}}^{- 1 / 2} \exp {- \frac{1}{2} g^{- 1} {ψ β}_{γ}^{T} (Z_{γ}^{T} Z_{γ}) β_{γ}}}{{(2 π)}^{- p} \det {g ψ^{- 1} {(Z^{T} Z)}^{- 1}}^{- 1 / 2} \exp {- \frac{1}{2} g^{- 1} {ψ β}_{γ}^{T} (Z_{γ}^{T} Z_{γ}) β_{γ}}} \\ \times π (β_{- γ} = 0 ∣ c, α, β_{γ}, ψ, g, X, M_{F}) \\ = {(2 π)}^{p - p_{γ}} {(g ψ^{- 1})}^{(p - p_{γ}) / 2} \det {(Z^{T} Z / Z_{γ}^{T} Z_{γ})}^{- 1 / 2} \\ \times π (β_{- γ} = 0 ∣ c, α, β_{γ}, ψ, g, X, M_{F}) . \end{matrix}

(12)

Here, $π (β_{- γ} ∣ c, α, β_{γ}, ψ, g, X, M_{F})$ is a $p_{- γ}$ dimensional multivariate normal distribution, that is,

\begin{matrix} [β_{- γ} ∣ c, α, β_{γ}, ψ, g, X, M_{F}] \\ ~ N_{p_{- γ}} (h_{- γ} + W_{21} W_{11}^{- 1} (β_{γ} - h_{γ}), W_{22} - W_{21} W_{11}^{- 1} W_{12}), \end{matrix}

with $h_{γ} = g / (g + 1) {(Z_{γ}^{T} Z_{γ})}^{- 1} Z_{γ}^{T} \log (c)$ , $h_{- γ} = g / (g + 1) {(Z_{- γ}^{T} Z_{- γ})}^{- 1} Z_{- γ}^{T} \times \log (c)$ , $W_{11}$ is a $p_{γ} \times p_{γ}$ matrix, $W_{22}$ is a $p_{- γ} \times p_{- γ}$ matrix and $W_{12} = W_{21}^{T}$ is a $p_{γ} \times p_{- γ}$ matrix, while $W_{21}$ is a transpose of $W_{12}$ , which is a $p_{- γ} \times p_{γ}$ matrix. Notice $W_{11}$ , $W_{12}$ , $W_{21}$ , and $W_{22}$ are defined through the partition of the matrix $W = ψ^{- 1} g / (g + 1) {(Z^{T} Z)}^{- 1}$ , that is,

W = (\begin{matrix} W_{11} & W_{12} \\ W_{21} & W_{22} \end{matrix}) .

(13)

Similarly, we can derive,

[β ∣ c, α, ψ, g, X, M_{F}] ~ N_{p} (g / (g + 1) {(Z^{T} Z)}^{- 1} Z^{T} \log (c), W) .

(14)

Then, after some math simplification, we can write out the quantity inside the expectation of Equation (10) as the expression below:

\begin{matrix} \frac{π (β = 0 ∣ c, α, ψ, g, X, M_{F})}{π (β = 0 ∣ ψ, g, M_{F}) π (g)} = {(g + 1)}^{p / 2} \exp {- \frac{ψ}{2} \frac{g}{g + 1} {(c)}^{T} Z {(Z^{T} Z)}^{- 1} Z^{T} (c)} \\ \frac{Γ (1 / 2)}{{(n / 2)}^{1 / 2}} g^{3 / 2} e^{n / (2 g)}, \end{matrix}

(15)

by noticing that $π (β = 0 ∣ c, α, ψ, g, X, M_{F})$ is the density of the multivariate normal distribution for $β$ evaluated at $0$ in Equation (14) and $π (β = 0 ∣ ψ, g, M_{F})$ is the prior distribution (5) of $β$ evaluated at $β = 0$ . Thus, according to Equations (12) and (15), we are able to conveniently compute the BFs for all model pairs using a single MCMC output from the joint posterior of the full model. For how to draw the MCMC samples from the joint posterior of $M_{F}$ , we have described in Section 2.2.

4. Model Selection Consistency

A major motivation to use the BF as the variable selection tool is its attractive model selection properties discussed in Liang et al. (2008). In Liang et al. (2008), they pointed out under the linear model setup, $pli m_{n} \Pr (M_{γ} ∣ D) = 1$ , where $M_{γ}$ denotes the true underlying model, $D$ presents the data, “plim” means the convergence in probability with respect to the sampling distribution under the true model $M_{γ}$ , and $\Pr (M_{γ} ∣ D)$ , the posterior probability of a model $M_{γ}$ can be rewritten as a representation only involving BFs and priors of models in the computation, that is,

\Pr (M_{γ} ∣ D) = \frac{\Pr (D ∣ M_{γ}) π (M_{γ})}{\sum_{γ'} \Pr (D ∣ M_{γ'}) π (M_{γ'})} = \frac{π (M_{γ})}{\sum_{γ'} π (M_{γ'}) BF [M_{γ'} : M_{γ}]} .

Further, in Theorem 3 of Liang et al. (2008), they showed that when the Zellner–Siow prior is used and for any model $M_{γ'}$ that does not contain $M_{γ}$ with the assumption below holds, that is,

li m_{n \to \infty} \frac{β_{γ}^{T} Z_{γ}^{T} (I - P_{γ'}) Z_{γ} β_{γ}}{n} = b_{γ'} \in (0, \infty),

(16)

where $P_{γ'}$ is the projection matrix onto the column space of $Z_{γ'}$ , then we have

pli m_{n} BF [M_{γ'} : M_{γ}] = 0, \forall M_{γ'} \neq M_{γ} .

(17)

In the proposed dynamic IRT model, modeling $\log (c_{i})$ at Level 3 exactly mimics the role of the response variable in a linear regression model, though $c_{i}$ is an unknown parameter rather than the observed response. Thus, for our great interest, we want to investigate whether the posterior model selection consistency still holds under this hierarchical model framework.

However, the proof for the model selection consistency in Liang et al. (2008) could not be directly applied here. To note, we could not employ the non-informative priors assigned in Liang et al. (2008) for the precision parameter $ψ$ and the intercept $α$ in our hierarchical model, which would lead to an improper posterior under the proposed model (c.f., the Supplemental File for details). Hence, we have to change the priors for $ψ$ and $α$ to be proper priors as specified in (4) and (6). Then, it is necessary to first show (17) still holds for the linear model under the proper priors, which Liang et al. (2008) did not prove in their paper. After that, we can apply the result under the proper priors to show the model selection consistency in the framework of our proposed dynamic IRT model.

4.1 Model Consistency of the Proposed Dynamic IRT Model

Let us first focus on the general linear model used by Liang et al. (2008) and add the superscript in $M_{γ}^{L}$ to distinguish the model associated with the general linear model instead of our proposed model, that is,

M_{γ}^{L} : μ = α 1_{n} + X_{γ} β_{γ},

(18)

where $μ$ is the mean vector for the responses in $M_{γ}^{L}$ , $α$ is an intercept, which is common to all possible linear models, $X_{γ}$ is a $n \times p_{γ}$ design matrix under the model $M_{γ}^{L}$ and $β_{γ}$ is the $p_{γ}$ -dimensional vector of nonzero regression coefficients.

Without loss of generality, we assume the columns of $X_{γ}$ have been centered, that is, $1^{T} X_{γ} = 0$ . Let us consider $M_{N}^{L} : μ = α 1_{n}$ to be the null model and to compare each model $M_{γ}^{L}$ with the null model $M_{N}^{L}$ through the hypotheses $H_{0} : β_{γ} = 0$ and $H_{1} : β_{γ} \in R^{p_{γ}}$ . Then, following the similar ideas of Liang et al. (2008), we can show the result stated in Liang et al. (2008) still holds and we summarize it in the lemma below. To save the space, we postpone the proof of Lemma 1 in Supplemental Appendix C (available in the online version of this article).

Lemma 1. Assume the condition (16) holds, the posterior probability under Zellner–Siow prior is consistent for model selection if we use proper priors (4) to (6) for the parameters $α$ , $ψ$ , and $β_{γ}$ in the general linear model (18).

Now, using Lemma 1, we can further prove the model selection consistency of the dynamic IRT model. The following Theorem 2 states the model selection consistency of the dynamic IRT model, and we provide the details of the proof in Supplemental Appendix D (available in the online version of this article).

Theorem 2. Assume marginal likelihoods of all possible models defined under our proposed hierarchical framework of dynamic IRT models are finite for any sample size $n$ . Denote $M_{γ}$ as the true model and for any model $M_{γ'}$ that does not contain the true model, we assume the condition (16) holds.

Further, follow the definition of $B F_{H} [M_{γ'} : M_{γ}]$ in (11), we have

pli m_{n} B F_{H} [M_{γ'} : M_{γ}] = 0,

(19)

for all $M_{γ'} \neq M_{γ}$ .

5. Simulation Study

Though the proposed BF approach is supported by Theorem 2 for the posterior model selection consistency, it is desirable to examine the performance of this approach in finite samples. Thus, we design and carry out several simulations in this section.

Let us consider four choices for the number of sample size of individuals, that is, $n = 30, 60, 90, 180$ with three settings of the study period $T_{i} = 50, 100, 150$ . Notice that we make $T_{i}$ be the same for each individual $i = 1, \dots, n$ in the simulation. According to the proposed dynamic IRT model (1) to (3), we further need to specify the number of items per test $K_{i, t, s}$ , the number of tests per day $S_{i, t}$ , and the time lapse $Δ_{i, t}$ between two distinctive test dates before we are able to generate the data from the proposed model. Following the settings of Wang et al. (2013), we set $K_{i, t, s} = 10$ , $S_{i, t} = 4$ , and $Δ_{i, t} = t + 10$ for $t \leq T_{i} / 2$ , and $Δ_{i, t} = t - 10$ for $T_{i} / 2 + 1 \leq t \leq T_{i}$ for any $i = 1, \dots, n$ , $t = 1, \dots, T_{i}$ and $s = 1, \dots, S_{i, t}$ . Also, because of the design of the EdSphere learning platform, following Wang et al. (2013), we use $ρ = 0.118$ and $σ^{2} = 0 . 7333^{2}$ . For the unknown parameters in the model, we randomly choose $τ_{i}^{- 1 / 2} ~^{i . i . d .} U (0.2, 0.4)$ , set $α = - 5$ , $ψ = 100$ , and $ϕ^{- 1 / 2} = 0.0055$ . We choose $p = 10$ with the three different setups of coefficients $β$ , as illustrated in Sections 5.1 to 5.3.

In addition, the first five covariates $z_{i 1}, \dots, z_{i 5}$ , for $i = 1, \dots, n$ are drawn from a multivariate normal distribution $N_{5} (μ_{2}, Σ_{2})$ , where $μ_{2} = (0, 0, 0, 0, 0)'$ and the variance-covariance matrix $Σ_{2}$ is

Σ_{2} = (\begin{matrix} 1 & 0.1 & 0.3 & 0.2 & 0.13 \\ 0.1 & 1 & 0.2 & 0.15 & 0.22 \\ 0.3 & 0.2 & 1 & 0.12 & 0.1 \\ 0.2 & 0.15 & 0.12 & 1 & 0.3 \\ 0.13 & 0.22 & 0.1 & 0.3 & 1 \end{matrix});

while the last five covariates $z_{i 6}, \dots, z_{i 10}$ for $i = 1, \dots, n$ are independently generated from a Bernoulli distribution with a success probability equal to .5. Once the covariates matrix is produced, we center each column of this matrix to be zero and scale each column to have unit variance in order to build the design matrix $Z$ . Then, using the values of $α$ , $β$ , and the design matrix $Z$ , we can generate $c_{i}$ ’s in Equation 3, after that in each case discussed below, the values of $θ_{i, t}$ ’s are created according to Equation (2), and finally, we generate the data $X_{i, t, s, l}$ ’s in our model from Equation (1).

Using the data generation scheme described above, we carry out the simulations with 100 different random seeds for each case of coefficients $β$ specified from Sections 5.1 to 5.3. In each generated dataset, we run 50,000 MCMC iterations to sample the posterior distribution of $M_{F}$ , and we discard the first 10,000 samples as a burn-in period.

5.1 Case 1: $β = (0.3, 0.5, 0, - 0.4, - 0.18, 0, 0, 0.3, 0, 0.15)'$

In this simulation study, from Table 2, we see that the percentage of using BF approach to select the right model among 100 simulations is steadily increasing when the sample size increases for all three settings of $T_{i}$ in our simulations. When $n = 180$ , we find that the percentages in all three $T_{i}$ choices are close to 1, which is consistent with our theoretical results in Theorem 2.

Table 2.

Percentage of the True Model Ranked as the “Best” Model by the Bayes Factor Approach Among 100 Simulations

n	T_i
n	50	100	150
30	0.65	0.66	0.69
60	0.76	0.78	0.77
90	0.81	0.85	0.83
180	0.93	0.92	0.95

For $β$ , we have collected the quantities for assessing the accuracy of its point estimation and interval estimation, which include the bias of the mean, root mean square error (RMSE) and coverage probability (CP). The CP computes whether the true value contains in the 95% credible intervals or not. Since we have 10 $β_{i}$ ’s in each setting, to save the space for creating a table, for example, in the column of Bias of Mean as shown in Table 3, we average the bias of the mean values of all $β_{i}$ ’s, and then in the column of Bias of Min and Bias of Max, we provide the minimum values and maximum values among all $β_{i}$ ’s, respectively, among the 100 simulations. Similar interpretations will apply to the rest columns in Table 3. From Table 3, we can see the biases of the estimation of all $β_{i}$ ’s are pretty small and the mean of the CPs are very close to the normal level of 0.95 or higher, which demonstrates our MCMC algorithms listed in the Supplemental Appendix A (available in the online version of this article) works well to recover the true values of $β$ .

Table 3.

Quantities for Assessing Estimation of $β$ , Case 1

	Bias			RMSE			CP
( $n$ , $T_{i}$ )	Mean	Min	Max	Mean	Min	Max	Mean	Min	Max
(30, 50)	0.0097	−0.0078	0.0107	0.1311	0.0251	0.3420	0.987	0.970	1.000
(30, 100)	0.0049	−0.0327	0.0316	0.1953	0.0564	0.3638	0.970	0.940	1.000
(30, 150)	0.0011	−0.0992	0.0944	0.6029	0.3393	0.7914	0.984	0.970	1.000
(60, 50)	−0.0011	−0.0610	0.0522	0.1790	0.0167	0.2170	0.967	0.940	0.990
(60, 100)	−0.0168	−0.0701	0.0119	0.2330	0.0598	0.3409	0.976	0.960	1.000
(60, 150)	0.0079	−0.0821	0.0722	0.4571	0.0913	0.5633	0.954	0.950	0.990
(90, 50)	−0.0076	−0.0320	0.0105	0.1463	0.0185	0.3132	0.968	0.940	1.000
(90, 100)	−0.0039	−0.0431	0.0295	0.1843	0.0671	0.3166	0.974	0.950	1.000
(90, 150)	0.0068	−0.0832	0.0613	0.3977	0.1923	0.6091	0.956	0.940	0.980
(180, 50)	0.0008	−0.0022	0.0026	0.0707	0.0087	0.1125	0.947	0.920	0.960
(180, 100)	−0.0017	−0.0063	0.0092	0.0962	0.0079	0.1022	0.965	0.95	1.000
(180, 150)	0.0009	−0.0052	0.0075	0.0913	0.0103	0.1201	0.982	0.970	1.000

Note. CP = coverage probability; RMSE = root mean square error.

5.2 Case 2: $β = (0.3, 0.5, 0.6, - 0.4, - 0.18, - 0.5, 0.2, 0.3, 0.33, 0.15)'$

In this simulation study, we assume the full model be the true model. Similar to the output in Case 1, from Table 4, we have seen that when $n$ is increasing, the probability of selecting the true model is also increasing, which further confirms the empirical and theoretical results agree with each other. In analogy to Case 1, we put the quantities for assessing $β$ in Table 5, where we can see the biases of the estimation of all $β_{i}$ ’s are also small and the mean of the CPs are very close to the normal level of 0.95 or higher, which suggests that our MCMC algorithms in Supplemental Appendix A (available in the online version of this article) work well in recovering the true values of $β$ when the full model is the truth.

Table 4.

Percentage of the True Model Ranked as the “Best” Model by the Bayes Factor Approach Among 100 Simulations

n	T _i
n	50	100	150
30	0.71	0.72	0.72
60	0.80	0.79	0.81
90	0.88	0.86	0.89
180	0.97	1.00	0.98

Table 5.

Quantities for Assessing Estimation of $β$ , Case 2

	Bias			RMSE			CP
( $n$ , $T_{i}$ )	Mean	Min	Max	Mean	Min	Max	Mean	Min	Max
(30, 50)	−0.0082	−0.0175	0.0253	0.2257	0.0972	0.3643	0.981	0.970	1.000
(30, 100)	−0.0055	−0.0592	0.0478	0.2311	0.0755	0.3533	0.975	0.950	1.000
(30, 150)	0.0017	−0.0855	0.0902	0.7712	0.2766	0.8253	0.962	0.940	0.990
(60, 50)	−0.0042	−0.0533	0.0682	0.1878	0.0822	0.3511	0.952	0.940	0.980
(60, 100)	0.0023	−0.0311	0.0204	0.2062	0.0374	0.5358	0.971	0.950	1.000
(60, 150)	0.0083	−0.0472	0.0329	0.4564	0.2312	0.7052	0.986	0.960	1.000
(90, 50)	−0.0069	−0.0635	0.0694	0.1712	0.0208	0.3515	0.961	0.940	0.990
(90, 100)	−0.0032	−0.0815	0.0624	0.2178	0.0377	0.5180	0.959	0.940	0.990
(90, 150)	0.0067	−0.0559	0.0610	0.4608	0.1762	0.6410	0.987	0.970	1.000
(180, 50)	0.0021	−0.0038	0.0109	0.1023	0.0180	0.2187	0.946	0.920	0.970
(180, 100)	0.0009	−0.0043	0.0086	0.0955	0.0104	0.1532	0.991	0.980	1.000
(180, 150)	−0.0016	−0.0036	0.0053	0.1276	0.0096	0.2239	0.983	0.960	1.000

Note. CP = coverage probability; RMSE = root mean square error.

5.3 Case 3: $β = (0, 0.5, 0, 0, 0, 0, 0, 0.3, 0, 0)'$

In this simulation, we consider a case when there are a fewer number of non-zero coefficients. In Table 6, we have also seen that there is an increasing trend for the probability of selecting the true models, which conforms to the theoretical trend indicated in Section 4. In Table 7, we have represented all quantities related to how well the estimation of $β$ is, as accompanied with the previous two cases, it further supports our MCMC algorithms can recover the truth of $β$ very well.

Table 6.

Percentage of the True Model Ranked as the “Best” Model by the Bayes Factor Approach Among 100 Simulations

n	T_i
n	50	100	150
30	0.62	0.63	0.67
60	0.73	0.75	0.77
90	0.79	0.85	0.83
180	0.91	0.93	0.93

Table 7.

Quantities for Assessing Estimation of $β$ , Case 3

	Bias			RMSE			CP
( $n$ , $T_{i}$ )	Mean	Min	Max	Mean	Min	Max	Mean	Min	Max
(30, 50)	−0.0050	−0.0483	0.0364	0.1607	0.0186	0.2250	0.961	0.950	1.000
(30, 100)	0.0079	−0.0225	0.0275	0.2674	0.0812	0.4471	0.967	0.940	1.000
(30, 150)	0.0023	−0.0846	0.0664	0.4861	0.3179	0.6939	0.982	0.970	1.000
(60, 50)	−0.0049	−0.0565	0.0369	0.2650	0.0223	0.4644	0.986	0.970	1.000
(60, 100)	0.0047	−0.0623	0.0409	0.3005	0.0270	0.3858	0.971	0.950	0.990
(60, 150)	0.0085	−0.0934	0.0886	0.5100	0.4425	0.7355	0.983	0.970	1.000
(90, 50)	−0.0021	−0.0519	0.0245	0.2298	0.0583	0.4057	0.992	0.980	1.000
(90, 100)	0.0079	−0.0390	0.0231	0.1645	0.0422	0.3196	0.995	0.980	1.000
(90, 150)	−0.0037	−0.0789	0.0619	0.3947	0.2282	0.5386	0.962	0.940	0.990
(180, 50)	−0.0009	−0.0033	0.0045	0.0806	0.0072	0.1282	0.996	0.990	1.000
(180, 100)	0.0013	−0.0061	0.0101	0.0913	0.0063	0.1103	0.987	0.960	1.000
(180, 150)	0.0007	−0.0082	0.0077	0.0792	0.0077	0.1176	0.997	0.990	1.000

Note. CP = coverage probability; RMSE = root mean square error.

6. Real Data Application

The simulation study in Section 5 shows the finite sample performance of our proposed method is consistent with the theoretical results we derive. Thus, we are confident to apply the proposed method to the EdSphere dataset. In this analysis, 54 students are randomly selected in the data, and 4 candidate covariates are used: gender, economic status, race, and preLexile score, where we denote them by $Z_{1}$ , $Z_{2}$ , $Z_{3}$ , and $Z_{4}$ , respectively. Among these covariates, the gender $Z_{1}$ is binary, where $Z_{1} = 1$ for male and $Z_{1} = 0$ for female; the economic status $Z_{2}$ has two categories, that is, disadvantaged and non-disadvantaged, based on the lunch status, homelessness state, and other district-level indicators of economic situation of a person, and thus we denote $Z_{2} = 1$ for disadvantaged and 0 otherwise; for the race, $Z_{3} = 1$ is used for White and $Z_{3} = 0$ for Non-white; $Z_{4}$ is a continuous covariate, which is a preLexile score used for quantifying individual’s initial ability before they enter the EdSphere learning platform. Among these 54 students, their distinctive test days range from 45 to 122, the number of tests per day range from 3 to 18, and the number of questions per test range from 2 to 78.

Then, we run the MCMC procedure with 120,000 iterations for the joint posterior of the unknown parameters in the full model and discard the first 20,000 MCMC samples as a burn-in period. Then, we employ the proposed BF approach in Section 3 to compute the BFs for different combinations of covariates used in Equation 3, and then the results of the BF are summarized in Table 8. A total of 120,000 iterations of the MCMC algorithm will take about 12 hr for a Core i5—8GB RAM laptop but the comparison between BFs only takes less than 2 min. Generally speaking, the computation speed for the model comparison part is very fast. We have employed the Geweke Convergence test (Geweke, 1992) and the Heidelberger & Welch stationary test (Heidelberger & Welch, 1983) to test the convergence of MCMC chains for all parameters and we have used the R package boa to perform these tests (Smith, 2007). For the Geweke Convergence test, the $p$ -values of the MCMC chains for all parameters are greater than 0.05; while regarding the Heidelberger & Welch stationary test, the test results of all parameters indicate “passing.” Thus, the MCMC chains of our parameters have successfully passed both the Geweke Convergence test and the Heidelberger & Welch stationary test.

Table 8.

Estimated the Logarithm of BF ( $\log (BF)$ ) of Each Model $M_{γ}$ Versus the Full Model $M_{F}$

Model	$\log (BF)$
{ $Z_{1}$ , $Z_{2}$ , $Z_{3}$ , $Z_{4}$ }	0
{ $Z_{1}$ , $Z_{2}$ , $Z_{3}$ }	−28.42
{ $Z_{2}$ , $Z_{3}$ , $Z_{4}$ }	−41.38
{ $Z_{1}$ , $Z_{2}$ , $Z_{4}$ }	−66.16
{ $Z_{2}$ , $Z_{3}$ }	−67.79
{ $Z_{2}$ , $Z_{4}$ }	−79.42
{ $Z_{1}$ , $Z_{3}$ , $Z_{4}$ }	−85.07
{ $Z_{2}$ }	−106.08
{ $Z_{1}$ , $Z_{4}$ }	−126.37
{ $Z_{3}$ , $Z_{4}$ }	−138.31
{ $Z_{1}$ , $Z_{2}$ }	−157.08
{ $Z_{4}$ }	−152.34
{ $Z_{1}$ , $Z_{3}$ }	−161.67
{ NULL}	−172.69
{ $Z_{1}$ }	−188.83
{ $Z_{3}$ }	−192.86

From Table 8, we see that the model with the largest BF comparing to the full model is the full model itself, that is, including the gender, race, economic status and preLexile in terms of explaining the growth factor $c_{i}$ . The BF of the full model $M_{F}$ to the second “preferred” model ${M_{1, 2, 3} : Z_{1}, Z_{2}, Z_{3}}$ is equal to $e^{28.42}$ , which is considered as having decisive evidence to favor the full model from Kass and Raftery (1995). We summarize the posterior mean of regression coefficients and their corresponding 95% credible interval (CI) in Table 9. Unlike in the simulations, we don’t know the ground truth of our selection in the application. But it is easy to see that all 95% CI of $β$ coefficients in Table 9 do not include zeroes, which supports from another perspective that the choice of the full model as the best seems correct.

Table 9.

Posterior Mean of $β$ and 95% Credible Intervals

$β$	Posterior mean	95% CI
$β_{1}$	−0.2215	−0.8422 to −0.1977
$β_{2}$	0.8721	0.4443 to 1.5412
$β_{3}$	0.3395	0.0955 to 0.9843
$β_{4}$	−0.2572	−0.7064 to −0.1062

7. Discussion

From the simulation results, we can see that the BF approach seems to be a reliable tool for variable selection in the proposed dynamic IRT model, especially when the sample size is relatively moderate to large. In general, computing the marginal likelihoods of hierarchical models with multi-levels is difficult if no analytical forms are available, and the usage of the BF may be limited. However, if one wants to perform the variable selection under a hierarchical model, at a single, “stand-alone” layer that consists of the linear structure similar to Level 3 as in our proposed hierarchical model, Proposition 1 may still be applied in a similar fashion to compute the BFs conveniently. The proof for posterior model selection consistency may also hold if (a) similar assumptions as those stated in Theorem 2 hold and (b) the layer we consider to do variable selections links with other layers in certain ways, for example, the parameters in the considered layer are not included in any other layers.

There are also other model selection criteria for hierarchical models, such as the widely applicable information criterion (WAIC; Watanabe, 2010) and DIC (Spiegelhalter et al., 2002) in the literature. When applying the DIC for model selection, often the “focus” of the parameters of interest may affect the results. One may integrate out all “intermediate” parameters to obtain the likelihood with respect to the parameters of interest, which is usually challenging in hierarchical models such as the complex IRT models discussed in the paper, or one may treat the likelihood conditioning on the “intermediate” parameters. But the latter way has several criticisms (F. Liu et al., 2022 and the reference therein). In our future work, it would be great interesting to compare the model selection performances of these criteria in comparison to the BF approach we used here.

The current variable selection algorithm requires comparisons of BFs for all possible submodels. Under a low-dimensional setting of predictors, the time taken to explore the model space is not a issue. For large amount of predictors, to go through the whole model space within a finite time interval might be impossible. Thus, how to explore the model space in order to find the best model using BFs will be an interesting topic to study.

Supplemental Material

sj-pdf-1-jeb-10.3102_10769986251314527 – Supplemental material for Bayesian Variable Selection in Dynamic Item Response Theory Models

Supplemental material, sj-pdf-1-jeb-10.3102_10769986251314527 for Bayesian Variable Selection in Dynamic Item Response Theory Models by Jingyu Sun, Yang Liu, Xiaojing Wang and Ming-Hui Chen in Journal of Educational and Behavioral Statistics

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The research of Dr. J. Sun and Dr. X. Wang was supported by the National Science Foundation under Grant No. 1848451.

ORCID iDs

Xiaojing Wang

Ming-Hui Chen

Authors

JINGYU SUN is a PhD student who graduated from the Department of Statistics at the University of Connecticut in 2023. He is currently working as a data scientist in Revolution Medicines, located at 700 Saginaw Drive, Redwood City, CA 94063; e-mail: jingyu.sun@uconn.edu. His research mainly focuses on Bayesian analysis, item response theory model, and cluster analysis.

YANG LIU is a staff research scientist at Upstart Network, 2950 S Delaware St, St #410, San Mateo, CA 94403, e-mail: yang.5.liu@uconn.edu. He graduated from the Department of Statistics at the University of Connecticut in 2019. His research mainly focuses on machine Learning, Bayesian statistics, item response theory modeling, and subgroup analysis.

XIAOJING WANG is an associate professor at the Department of Statistics at the University of Connecticut, 215 Glenbrook Road, Unit 4120, Storrs, CT 06269-4120, e-mail: xiaojing.wang@uconn.edu. Her research mainly focuses on the latent variable models, nonparametric and semiparametric methods, subgroup analysis, Bayesian modeling and methodology, and personalized learning.

MING-HUI CHEN is a Board of Trustees Distinguished Professor and Head of the Department of Statistics at the University of Connecticut, 215 Glenbrook Road, Unit 4120, Storrs, CT 06269-4120, e-mail: ming-hui.chen@uconn.edu. His research mainly focuses on Bayesian Statistical Methodology, Bayesian Computation, Big Data, Categorical Data Analysis, Design of Bayesian Clinical Trials, DNA Microarray Data Analysis, Missing Data Analysis (EM, MCEM, and Bayesian), Monte Carlo Methodology, Prior Elicitation, Statistical Modeling, Survival Data Analysis, and Variable Selection.

References

Albert

J. H.

Chib

(1993). Bayesian analysis of binary and polychotomous response data. Journal of the American statistical Association, 88(422), 669–679.

Berger

J. O.

Pericchi

L. R.

(1996). The intrinsic Bayes factor for model selection and prediction. Journal of the American Statistical Association, 91(433), 109–122.

Birnbaum

A. L.

(1968). Some latent trait models and their use in inferring an examinee’s ability. In Lord

F. M.

Novick

M. R.

(Eds.), Statistical theories of mental test scores (pp. 397–479). Addison-Wesley, Reading.

Bock

R. D.

Mislevy

R. J.

(1989). A hierarchical item response model for educational testing. In Bock

R. Darrell

(Ed.), Multilevel analysis of educational data (pp. 57–74). Academic Press, San Diego.

Castellanos

M. E.

Garcia-Donato

Cabras

(2021). A model selection approach for variable selection with censored data. Bayesian Analysis, 16(1), 271–300.

Chen

M.-H.

(2005). Computing marginal likelihoods from a single MCMC output. Statistica Neerlandica, 59(1), 16–29.

Chen

Silva Filho

Prudencio

R. B.

Diethe

Flach

(2019).

β^{3}

—irt: A new item response model and its applications. In The 22nd international conference on artificial intelligence and statistics, PMLR: Volume 89, Naha, Okinawa, Japan (pp. 1013–1021).

Chevallier

Fearnhead

Sutton

(2022). Reversible jump PDMP samplers for variable selection. Journal of the American Statistical Association, 18, 1–31.

Chib

(1995). Marginal likelihood from the Gibbs output. Journal of the American Statistical Association, 90(432), 1313–1321.

10.

Dickey

J. M.

Lientz

(1970). The weighted likelihood ratio, sharp hypotheses about chances, the order of a Markov chain. The Annals of Mathematical Statistics, 41, 214–226.

11.

Fernandez

Ley

Steel

M. F.

(2001). Benchmark priors for Bayesian model averaging. Journal of Econometrics, 100(2), 381–427.

12.

Fox

J.-P.

(2005). Multilevel IRT using dichotomous and polytomous response data. British Journal of Mathematical and Statistical Psychology, 58(1), 145–172.

13.

Geweke

(1992). Evaluating the accuracy of sampling-based approaches to the calculations of posterior moments. Bayesian Statistics, 4, 641–649.

14.

Hamza

M. A.

Hassine

S. H.

Abunadi

Al-Wesabi

F. N.

Alsolai

Hilal

A. M.

Motwakel

(2022). Feature selection with optimal stacked sparse autoencoder for data mining. Computers, Materials & Continua, 72(2), 2581–2596.

15.

Heck

D. W.

Boehm

Böing-Messing

Bürkner

P.-C.

Derks

Dienes

Karimova

Kiers

Klugkist

Kuiper

R. M.

Lee

M. D.

Leenders

Leplaa

H. J.

Linde

Meijerink

Moerbeek

…Hoijtink

(2022). A review of applications of the Bayes factor in psychological research. Psychological Methods, 28(3), 558–579.

16.

Heidelberger

Welch

P. D.

(1983). Simulation run length control in the presence of an initial transient. Operations Research, 31(6), 1109–1144.

17.

Irwing

Booth

Hughes

D. J.

(2018). The Wiley handbook of psychometric testing: A multidisciplinary reference on survey, scale and test development. John Wiley & Sons.

18.

Kass

R. E.

Raftery

A. E.

(1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773–795.

19.

Koop

Korobilis

(2018). Bayesian dynamic variable selection in high dimensions. arXiv preprint arXiv:1809.03031.

20.

Clyde

M. A.

(2018). Mixtures of g-priors in generalized linear models. Journal of the American Statistical Association, 113(524), 1828–1845.

21.

Liang

Paulo

Molina

Clyde

M. A.

Berger

J. O.

(2008). Mixtures of g priors for Bayesian variable selection. Journal of the American Statistical Association, 103(481), 410–423.

22.

Liu

Wang

Hancock

Chen

M.-H.

(2022). Bayesian model assessment for jointly modeling multidimensional response data with application to computerized testing. Psychometrika, 87(4), 1290–1317.

23.

Liu

Cao

Wang

Chen

M.-H.

(2019). A comparison of Monte Carlo methods for computing marginal likelihoods of item response theory models. Journal of the Korean Statistical Society, 48(4), 503–512.

24.

Liu

Wang

(2020). Bayesian nonparametric monotone regression of dynamic latent traits in item response theory models. Journal of Educational and Behavioral Statistics, 45(3), 274–296.

25.

Lord

F. M.

(1980). Applications of item response theory to practical testing problems. Lawrence Erlbaum Associates.

26.

Muttakin

Wang

J.-T.

Mulyanto

Leu

J.-S.

(2021). Evaluation of feature selection methods on psychosocial education data using additive ratio assessment. Electronics, 11(1), 114.

27.

Newton

M. A.

Raftery

A. E.

(1994). Approximate Bayesian inference by the weighted likelihood Bootstrap. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 56, 3–48.

28.

Park

Casella

(2008). The Bayesian lasso. Journal of the American Statistical Association, 103(482), 681–686.

29.

Pan

Yin

(2017). Using the Bayes factors to evaluate person fit in the item response theory. Applied Measurement in Education, 30(3), 213–227.

30.

Pericchi

L. R.

(2005). Model selection and hypothesis testing based on objective probabilities and bayes factors. Handbook of Statistics, 25, 115–149.

31.

Rahman

Setiawan

N. A.

Permanasari

A. E.

(2017). Feature selection methods in improving accuracy of classifying students’ academic performance. In 2017 2nd international conferences on information technology, information systems and electrical engineering (ICITISEE), IEEE, Yogyakarta, Indonesia, pp. 267–271.

32.

Rasch

(1960). Studies in mathematical psychology: I. probabilistic models for some intelligence and attainment tests. Danish Institute for Educational Research.

33.

Ročková

George

E. I.

(2018). The spike-and-slab lasso. Journal of the American Statistical Association, 113(521), 431–444.

34.

Shim

Bonifay

Wiedermann

(2022). Parsimonious asymmetric item response theory modeling with the complementary log-log link. Behavior Research Methods, 55, 200–219.

35.

Smith

B. J.

(2007). boa: An r package for MCMC output convergence assessment and posterior inference. Journal of Statistical Software, 21(11), 1–37.

36.

Spiegelhalter

D. J.

Best

N. G.

Carlin

B. P.

Van Der Linde

(2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64(4), 583–639.

37.

Stenner

(2010). Using technology to merge assessment and instruction. In 2nd international conference for teaching and learning with technology. Available from MetaMetrics, Inc.

38.

Swartz

C. W.

Hanlon

S. T.

Childress

E. L.

Stenner

A. J.

(2016). An approach to design-based implementation research to inform development of edsphere®: A brief history about the evolution of one personalized learning platform. In Rosen

Ferrara

Mosharraf

(Eds.) Handbook of research on technology tools for real-world skill development (pp. 284–318). IGI Global Scientific Publishing, Pennsylvania, PA.

39.

Swartz

C. W.

Hanlon

S. T.

Stenner

A. J.

Burdick

D. S.

Emerson

(2011). Edsphere™: Using technology to enhance literacy through deliberate practice. Metametrics Research Brief.

40.

Tijmstra

Bolsinova

(2019). Bayes factors for evaluating latent monotonicity in polytomous item response theory models. Psychometrika, 84(3), 846–869.

41.

Tutz

(2022). Flexible item response models for count data: The count thresholds model. Applied Psychological Measurement, 46, 643–661.

42.

Van der Linden

W. J.

(2018). Handbook of item response theory: Three volume set. CRC Press.

43.

Verdinelli

Wasserman

(1995). Computing Bayes factors using a generalization of the savage-dickey density ratio. Journal of the American Statistical Association, 90(430), 614–618.

44.

Wang

Berger

J. O.

Burdick

D. S.

(2013). Bayesian analysis of dynamic item response models in educational testing. The Annals of Applied Statistics, 7(1), 126–153.

45.

Watanabe

(2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11, 3571–3594.

46.

H.-H.

Ferreira

M. A.

Gompper

M. E.

(2016). Consistency of hyper-g-prior-based bayesian variable selection for generalized linear models. Brazilian Journal of Probability and Statistics, 30(4), 691–709.

47.

Xie

Lewis

P. O.

Fan

Kuo

Chen

M.-H.

(2011). Improving marginal likelihood estimation for Bayesian phylogenetic model selection. Systematic Biology, 60(2), 150–160.

48.

Zellner

(1986). On assessing prior distributions and Bayesian regression analysis with g-prior distributions. In Goel

Zellner A

(Eds.), Bayesian inference and decision techniques (pp. 233-243). Elsevier Science Publishers, Inc., New York, NY.

49.

Zellner

Siow

(1980). Posterior odds ratios for selected regression hypotheses. Trabajos de estadística y de investigación operativa, 31(1), 585–603.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.30 MB

Bayesian Variable Selection in Dynamic Item Response Theory Models

Abstract

Keywords

1. Introduction

1.1 EdSphere Dataset and Dynamic IRT Model

1.2 Variable Selection With Bayes Factor

2. The Proposed Dynamic IRT Model and Bayesian Inference

2.1 Model Settings and Prior Specifications

2.2 Posterior Distribution and Computation

3. Computing BFs With a Single MCMC Output

4. Model Selection Consistency

4.1 Model Consistency of the Proposed Dynamic IRT Model

5. Simulation Study

5.1 Case 1: β = ( 0 . 3 , 0 . 5 , 0 , − 0 . 4 , − 0 . 18 , 0 , 0 , 0 . 3 , 0 , 0 . 15 ) ′

5.2 Case 2: β = ( 0 . 3 , 0 . 5 , 0 . 6 , − 0 . 4 , − 0 . 18 , − 0 . 5 , 0 . 2 , 0 . 3 , 0 . 33 , 0 . 15 ) ′

5.3 Case 3: β = ( 0 , 0 . 5 , 0 , 0 , 0 , 0 , 0 , 0 . 3 , 0 , 0 ) ′

6. Real Data Application

7. Discussion

Supplemental Material

sj-pdf-1-jeb-10.3102_10769986251314527 – Supplemental material for Bayesian Variable Selection in Dynamic Item Response Theory Models

Footnotes

Declaration of Conflicting Interests

Funding

ORCID iDs

Authors

References

Supplementary Material

5.1 Case 1: $β = (0.3, 0.5, 0, - 0.4, - 0.18, 0, 0, 0.3, 0, 0.15)'$

5.2 Case 2: $β = (0.3, 0.5, 0.6, - 0.4, - 0.18, - 0.5, 0.2, 0.3, 0.33, 0.15)'$

5.3 Case 3: $β = (0, 0.5, 0, 0, 0, 0, 0, 0.3, 0, 0)'$