## Bayesian binary item response theory models using bayesmh

This post was written jointly with Yulia Marchenko, Executive Director of Statistics, StataCorp.

Table of Contents

__Overview__

1PL model

2PL model

3PL model

4PL model

5PL model

Conclusion

Overview

Item response theory (IRT) is used for modeling the relationship between the latent abilities of a group of subjects and the examination items used for measuring their abilities. Stata 14 introduced a suite of commands for fitting IRT models using maximum likelihood; see, for example, the blog post Spotlight on irt by Rafal Raciborski and the [IRT] Item Response Theory manual for more details. In this post, we demonstrate how to fit Bayesian binary IRT models by using the **redefine()** option introduced for the bayesmh command in Stata 14.1. We also use the likelihood option **dbernoulli()** available as of the update on 03 Mar 2016 for fitting Bernoulli distribution. If you are not familiar with the concepts and jargon of Bayesian statistics, you may want to watch the introductory videos on the Stata Youtube channel before proceeding.

Introduction to Bayesian analysis, part 1 : The basic concepts

Introduction to Bayesian analysis, part 2: MCMC and the Metropolis-Hastings algorithm

We use the abridged version of the mathematics and science data from DeBoeck and Wilson (2004), **masc1**. The dataset includes 800 student responses to 9 test questions intended to measure mathematical ability.

The irt suite fits IRT models using data in the wide form – one observation per subject with items recorded in separate variables. To fit IRT models using bayesmh, we need data in the long form, where items are recorded as multiple observations per subject. We thus reshape the dataset in a long form: we have a single binary response variable, **y**, and two index variables, **item** and **id**, which identify the items and subjects, respectively. This allows us to formulate our IRT models as multilevel models. The following commands load and prepare the dataset.

. webuse masc1 (Data from De Boeck & Wilson (2004)) . generate id = _n . quietly reshape long q, i(id) j(item) . rename q y

To ensure that we include all levels of **item** and **id** in our models, we use **fvset base none** to keep the base categories.

. fvset base none id item

In what follows, we present eight Bayesian binary IRT models increasing in complexity and explanatory power. We perform Bayesian model comparison to gain insight into what would be the more appropriate model for the data at hand.

For high-dimensional models such as IRT models, you may see differences in the estimation results between different platforms or different flavors of Stata because of the nature of the Markov chain Monte Carlo (MCMC) sampling and finite numerical precision. These differences are not a source of concern; they will be within the range of the MCMC variability and will lead to similar inferential conclusions. The differences will diminish as the MCMC sample size increases. The results in this post are obtained from Stata/SE on the 64-bit Linux platform using the default 10,000 MCMC sample size.

Let the items be indexed by \(i=1,\dots,9\) and the subjects by \(j=1,\dots,800\). Let \(\theta_j\) be the latent mathematical ability of subject \(j\), and let \(Y_{ij}\) be the response of subject \(j\) to item \(i\).

1PL model

In the one-parameter logistic (1PL) model, the probability of getting a correct response is modeled as an inverse-logit function of location parameters \(b_i\), also called item difficulties, and a common slope parameter \(a\), also called item discrimination:

\[

P(Y_{ij}=1) = {\rm InvLogit}\{a(\theta_j-b_i)\} =

\frac{\exp\{a(\theta_j-b_i)\}}{1+\exp\{a(\theta_j-b_i)\}}

\]

Typically, the abilities are assumed to be normally distributed:

\[

\theta_j \sim {\rm N}(0,1)

\]

In a multilevel framework, the \(\theta_j\)’s represent random effects. In a Bayesian framework, we use the term “random effects” to refer to the parameters corresponding to levels of grouping variables identifying the hierarchy of the data.

A Bayesian formulation of the 1PL model also requires prior specification for the model parameters \(a\) and \(b_i\). The discrimination parameter \(a\) is assumed to be positive and is often modeled in the log scale. Because we have no prior knowledge about the discrimination and difficulty parameters, we assume that the prior distributions of \(\ln(a)\) and \(b_i\) have support on the whole real line, are symmetric, and are centered at 0. A normal prior distribution is thus a natural choice. We furthermore assume that \(\ln(a)\) and \(b_i\) are close to 0 and have prior variance of 1, which is an entirely subjective decision. We thus assign \(\ln(a)\) and \(b_i\) standard normal prior distributions:

\[\ln(a) \sim {\rm N}(0, 1)\] \[b_i \sim {\rm N}(0, 1) \]

To specify the likelihood function of the 1PL model in **bayesmh**, we use a nonlinear equation specification for the response variable **y**. The direct nonlinear specification for this model is

bayesmh y = ({discrim}*({subj:i.id}-{diff:i.item})), likelihood(logit) ...

where **{discrim}** is the discrimination parameter \(a\), **{subj:i.id}** are latent abilities \(\theta_j\), and **{diff:i.item}** are item difficulties \(b_i\). The logit model is used for the probability of a success, \(P(Y_{ij}=1)\). Specification **{subj:i.id}** in the above nonlinear expression is viewed as a substitutable expression for linear combinations of indicators associated with the **id** variable and parameters \(\theta_j\). This specification may be computationally prohibitive with a large number of subjects. A more efficient solution is to use the **redefine()** option to include subject random effects \(\theta_j\) in the model. The same argument may apply to the **{diff:i.item}** specification when there are many items. Thus, it may be computationally convenient to treat the \(b_i\) parameters as “random effects” in the specification and use the **redefine()** option to include them in the model.

A more efficient specification is thus

bayesmh y = ({discrim}*({subj:}-{diff:})), likelihood(logit) /// redefine(subj:i.id) redefine(diff:i.item) ... ///

where **{subj:}** and **{diff:}** in the nonlinear specification now represent the \(\theta_j\) and \(b_i\) parameters, respectively, without using expansions into linear combinations of indicator variables.

Below, we show the full **bayesmh** specification of the 1PL model and the output summary. In our examples, we treat the abilities **{subj:i.id}** as nuisance parameters and exclude them from the final results. The discrimination model parameter **{discrim}** must be positive and is thus initialized with 1. A longer burn-in period, **burnin(5000)**, allows for longer adaptation of the MCMC sampler, which is needed given the large number of parameters in the model. Finally, the estimation results are stored for later model comparison.

. set seed 14 . bayesmh y = ({discrim}*({subj:}-{diff:})), likelihood(logit) /// > redefine(diff:i.item) redefine(subj:i.id) /// > prior({subj:i.id}, normal(0, 1)) /// > prior({discrim}, lognormal(0, 1)) /// > prior({diff:i.item}, normal(0, 1)) /// > init({discrim} 1) exclude({subj:i.id}) /// > burnin(5000) saving(sim1pl, replace) Burn-in ... Simulation ... Model summary ------------------------------------------------------------------------------ Likelihood: y ~ logit({discrim}*(xb_subj-xb_diff)) Priors: {diff:i.item} ~ normal(0,1) (1) {subj:i.id} ~ normal(0,1) (2) {discrim} ~ lognormal(0,1) ------------------------------------------------------------------------------ (1) Parameters are elements of the linear form xb_diff. (2) Parameters are elements of the linear form xb_subj. Bayesian logistic regression MCMC iterations = 15,000 Random-walk Metropolis-Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 7,200 Acceptance rate = .3074 Efficiency: min = .02691 avg = .06168 Log marginal likelihood = . max = .09527 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval] -------------+---------------------------------------------------------------- diff | item | 1 | -.6934123 .0998543 .003576 -.6934789 -.8909473 -.4917364 2 | -.1234553 .0917187 .002972 -.1241642 -.3030341 .0597863 3 | -1.782762 .1323252 .00566 -1.781142 -2.05219 -1.534451 4 | .3152835 .0951978 .003289 .3154714 .1279147 .4981263 5 | 1.622545 .127213 .005561 1.619388 1.377123 1.883083 6 | .6815517 .0978777 .003712 .6788345 .4911366 .881128 7 | 1.303482 .1173994 .005021 1.302328 1.084295 1.544913 8 | -2.353975 .1620307 .008062 -2.351207 -2.672983 -2.053112 9 | -1.168668 .1120243 .004526 -1.163922 -1.392936 -.9549209 -------------+---------------------------------------------------------------- discrim | .8644787 .0439804 .002681 .8644331 .7818035 .9494433 ------------------------------------------------------------------------------ file sim1pl.dta saved . estimates store est1pl

The sampling efficiency is acceptable, about 6% on average, with no indication of convergence problems. Although detailed convergence inspection of all parameters is outside the scope of this post, we recommend that you do so by using, for example, the bayesgraph diagnostics command.

Though we used informative priors for the model parameters, the estimation results from our Bayesian model are not that different from the maximum likelihood estimates obtained using the **irt 1pl** command (see example 1 in [IRT] irt 1pl). For example, the posterior mean estimate for **{discrim}** is 0.86 with an MCMC standard error of 0.003, whereas **irt 1pl** reports 0.85 with a standard error of 0.05.

The log-marginal likelihood is reported missing because we have excluded the **{subj:i.id}** parameters from the simulation results and the Laplace-Metropolis estimator of the log-marginal likelihood is not available in such cases. This estimator requires simulation results for all model parameters to compute the log-marginal likelihood.

2PL model

The two-parameter logistic (2PL) model extends the 1PL model by allowing for item-specific discrimination. The probability of correct response is now modeled as a function of item-specific slope parameters \(a_i\):

\[

P(Y_{ij}=1) = {\rm InvLogit}\{a_i(\theta_j-b_i)\} =

\frac{\exp\{a_i(\theta_j-b_i)\}}{1+\exp\{a_i(\theta_j-b_i)\}}

\]

The prior specification for \(\theta_j\) stays the same as in the 1PL model. We will, however, apply more elaborate prior specifications for the \(a_i\)’s and \(b_i\)’s. It is a good practice to use proper prior specifications without overwhelming the evidence from the data. The impact of the priors can be controlled by introducing additional hyperparameters. For example, Kim and Bolt (2007) proposed the use of a normal prior for the difficulty parameters with unknown mean and variance. Extending this approach to the discrimination parameters as well, we apply a hierarchical Bayesian model in which the \(\ln(a_i)\) and \(b_i\) parameters have the following prior specifications:

\[ \ln(a_i) \sim {\rm N}(\mu_a, \sigma_a^2) \] \[ b_i \sim {\rm N}(\mu_b, \sigma_b^2) \]

The mean hyperparameters, \(\mu_a\) and \(\mu_b\), and variance hyperparameters, \(\sigma_a^2\) and \(\sigma_b^2\), require informative prior specifications. We assume that the means are centered at 0 with a variation of 0.1:

\[

\mu_a, \mu_b \sim {\rm N}(0, 0.1)

\]

To lower the variability of the \(\ln(a_i)\) and \(b_i\) parameters, we apply an inverse-gamma prior with shape 10 and scale 1 for the variance parameters:

\[

\sigma_a^2, \sigma_b^2 \sim {\rm InvGamma}(10, 1)

\]

Thus, the prior mean of \(\sigma_a^2\) and \(\sigma_b^2\) is about 0.1.

In the **bayesmh** specification, the hyperparameters \(\mu_a\), \(\mu_b\), \(\sigma_a^2\), and \(\sigma_a^2\) are denoted as **{mu_a}**, **{mu_b}**, **{var_a}**, and **{var_b}**, respectively. We use the **redefine(discrim:i.item)** option to include in the model the discrimination parameters \(a_i\), referred to as **{discrim:}** in the likelihood specification.

Regarding the MCMC simulation, we change some of the default options. The hyperparameters **{mu_a}**, **{mu_b}**, **{var_a}**, and **{var_b}** are placed in separate blocks to improve the simulation efficiency. The discrimination parameters **{discrim:i.item}** must be positive and are thus initialized with 1s.

. set seed 14 . bayesmh y = ({discrim:}*({subj:}-{diff:})), likelihood(logit) /// > redefine(discrim:i.item) redefine(diff:i.item) /// > redefine(subj:i.id) /// > prior({subj:i.id}, normal(0, 1)) /// > prior({discrim:i.item}, lognormal({mu_a}, {var_a})) /// > prior({diff:i.item}, normal({mu_b}, {var_b})) /// > prior({mu_a} {mu_b}, normal(0, 0.1)) /// > prior({var_a} {var_b}, igamma(10, 1)) /// > block({mu_a mu_b var_a var_b}, split) /// > init({discrim:i.item} 1) /// > exclude({subj:i.id}) burnin(5000) saving(sim2pl, replace) Burn-in ... Simulation ... Model summary ------------------------------------------------------------------------------ Likelihood: y ~ logit(xb_discrim*(xb_subj-xb_diff)) Priors: {discrim:i.item} ~ lognormal({mu_a},{var_a}) (1) {diff:i.item} ~ normal({mu_b},{var_b}) (2) {subj:i.id} ~ normal(0,1) (3) Hyperpriors: {mu_a mu_b} ~ normal(0,0.1) {var_a var_b} ~ igamma(10,1) ------------------------------------------------------------------------------ (1) Parameters are elements of the linear form xb_discrim. (2) Parameters are elements of the linear form xb_diff. (3) Parameters are elements of the linear form xb_subj. Bayesian logistic regression MCMC iterations = 15,000 Random-walk Metropolis-Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 7,200 Acceptance rate = .3711 Efficiency: min = .01617 avg = .04923 Log marginal likelihood = . max = .1698 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval] -------------+---------------------------------------------------------------- discrim | item | 1 | 1.430976 .1986011 .010953 1.413063 1.089405 1.850241 2 | .6954823 .1081209 .004677 .6897267 .4985004 .9276975 3 | .9838528 .1343908 .009079 .9780275 .7506566 1.259427 4 | .8167792 .1169157 .005601 .8136229 .5992495 1.067578 5 | .9402715 .1351977 .010584 .9370298 .6691103 1.214885 6 | .9666747 .1420065 .008099 .9616285 .7038868 1.245007 7 | .5651287 .0864522 .006201 .5617302 .3956216 .7431265 8 | 1.354053 .2048404 .015547 1.344227 .9791096 1.761437 9 | .7065096 .1060773 .006573 .6999745 .5102749 .9271799 -------------+---------------------------------------------------------------- diff | item | 1 | -.5070314 .0784172 .003565 -.507922 -.671257 -.3596057 2 | -.1467198 .117422 .003143 -.1456633 -.3895978 .0716841 3 | -1.630259 .1900103 .013494 -1.612534 -2.033169 -1.304171 4 | .3273735 .1073891 .003565 .3231703 .1248782 .5492114 5 | 1.529584 .1969554 .01549 1.507982 1.202271 1.993196 6 | .6325194 .115724 .005613 .6243691 .4272131 .8851649 7 | 1.827013 .2884057 .019582 1.79828 1.349654 2.490633 8 | -1.753744 .1939559 .014743 -1.738199 -2.211475 -1.438146 9 | -1.384486 .2059005 .012105 -1.361195 -1.838918 -1.059687 -------------+---------------------------------------------------------------- mu_a | -.1032615 .1148176 .003874 -.102376 -.3347816 .1277031 var_a | .1129835 .0356735 .001269 .1056105 .063403 .1981331 mu_b | -.0696525 .2039387 .004949 -.072602 -.4641566 .3298393 var_b | .6216005 .2023137 .008293 .5843444 .3388551 1.101153 ------------------------------------------------------------------------------ file sim2pl.dta saved . estimates store est2pl

The average simulation efficiency is about 5%, but some of the parameters converge slower than the others, such as **{diff:7.item}**, which has the largest MCMC standard error (0.02) among the difficulty parameters. If this was a rigorous study, to lower the MCMC standard errors, we would recommend longer simulations with MCMC sample sizes of at least 50,000.

We can compare the 1PL and 2PL models by using the deviance information criterion (DIC) available with the bayesstats ic command.

. bayesstats ic est1pl est2pl, diconly Deviance information criterion ------------------------ | DIC -------------+---------- est1pl | 8122.428 est2pl | 8055.005 ------------------------

DIC is often used in Bayesian model selection as an alternative to AIC and BIC criteria and can be easily obtained from an MCMC sample. Larger MCMC samples produce more reliable DIC estimates. Because different MCMC samples produce different sample DIC values and the sample approximation error in calculating DIC is not known, one should not rely solely on DIC when choosing a model.

Lower DIC values indicate better fit. The DIC of the 2PL model (8,055) is markedly lower than the DIC of the 1PL model (8,122), implying better fit of the 2PL model.

3PL model

The three-parameter logistic (3PL) model introduces lower asymptote parameters \(c_i\), also called guessing parameters. The probability of giving a correct response is given by

\[

P(Y_{ij}=1) = c_i + (1-c_i){\rm InvLogit}\{a_i(\theta_j-b_i)\} ,\ c_i > 0

\]

The guessing parameters may be difficult to estimate using maximum likelihood. Indeed, the **irt 3pl** command with the **sepguessing** option fails to converge, as you can verify by typing

. irt 3pl q1-q9, sepguessing

on the original dataset.

It is thus important to specify an informative prior for \(c_i\). We assume that the prior mean of the guessing parameters is about 0.1 and thus apply

\[

c_i \sim {\rm InvGamma}(10, 1)

\]

Similarly to the discrimination and difficulty parameters, the \(c_i\)’s are introduced as random-effects parameters in the **bayesmh** specification and are referred to as **{gues:}** in the likelihood specification.

Unlike 1PL and 2PL models, we cannot use the **likelihood(logit)** option to model the probability of success because the probability of correct response is no longer an inverse-logit transformation of the parameters. Instead, we use **likelihood(dbernoulli())** to model the probability of success of a Bernoulli outcome directly.

To have a valid initialization of the MCMC sampler, we assign the \(c_i\)’s positive starting values, 0.1.

. set seed 14 . bayesmh y, likelihood(dbernoulli({gues:}+(1-{gues:})* /// > invlogit({discrim:}*({subj:}-{diff:})))) /// > redefine(discrim:i.item) redefine(diff:i.item) /// > redefine(gues:i.item) redefine(subj:i.id) /// > prior({subj:i.id}, normal(0, 1)) /// > prior({discrim:i.item}, lognormal({mu_a}, {var_a})) /// > prior({diff:i.item}, normal({mu_b}, {var_b})) /// > prior({gues:i.item}, igamma(10, 1)) /// > prior({mu_a} {mu_b}, normal(0, 0.1)) /// > prior({var_a} {var_b}, igamma(10, 1)) /// > block({mu_a mu_b var_a var_b}, split) /// > init({discrim:i.item} 1 {gues:i.item} 0.1) /// > exclude({subj:i.id}) burnin(5000) saving(sim3pls, replace) Burn-in ... Simulation ... Model summary ------------------------------------------------------------------------------ Likelihood: y ~ binomial(xb_gues+(1-xb_gues)*invlogit(xb_discrim*(xb_subj-xb_diff)),1) Priors: {discrim:i.item} ~ lognormal({mu_a},{var_a}) (1) {diff:i.item} ~ normal({mu_b},{var_b}) (2) {gues:i.item} ~ igamma(10,1) (3) {subj:i.id} ~ normal(0,1) (4) Hyperpriors: {mu_a mu_b} ~ normal(0,0.1) {var_a var_b} ~ igamma(10,1) ------------------------------------------------------------------------------ (1) Parameters are elements of the linear form xb_discrim. (2) Parameters are elements of the linear form xb_diff. (3) Parameters are elements of the linear form xb_gues. (4) Parameters are elements of the linear form xb_subj. Bayesian Bernoulli model MCMC iterations = 15,000 Random-walk Metropolis-Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 7,200 Acceptance rate = .3496 Efficiency: min = .0148 avg = .03748 Log marginal likelihood = . max = .2044 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval] -------------+---------------------------------------------------------------- discrim | item | 1 | 1.712831 .2839419 .018436 1.681216 1.232644 2.351383 2 | .8540871 .1499645 .008265 .8414399 .6058463 1.165732 3 | 1.094723 .1637954 .01126 1.081756 .817031 1.454845 4 | 1.090891 .2149095 .013977 1.064651 .7488589 1.588164 5 | 1.363236 .2525573 .014858 1.338075 .9348136 1.954695 6 | 1.388325 .3027436 .024245 1.336303 .9466695 2.068181 7 | .9288217 .2678741 .021626 .8750048 .5690308 1.603375 8 | 1.457763 .2201065 .01809 1.438027 1.068937 1.940431 9 | .7873631 .127779 .007447 .7796568 .563821 1.06523 -------------+---------------------------------------------------------------- diff | item | 1 | -.2933734 .0976177 .006339 -.2940499 -.4879558 -.0946848 2 | .2140365 .157158 .008333 .2037788 -.0553537 .5550411 3 | -1.326351 .1981196 .013101 -1.326817 -1.706671 -.9307443 4 | .6367877 .1486799 .007895 .6277349 .3791045 .9509913 5 | 1.616056 .1799378 .00966 1.606213 1.303614 2.006817 6 | .8354059 .124184 .00656 .8191839 .614221 1.097801 7 | 2.066205 .3010858 .018377 2.034757 1.554484 2.709601 8 | -1.555583 .1671435 .012265 -1.54984 -1.89487 -1.267001 9 | -.9775626 .2477279 .016722 -.9936727 -1.431964 -.4093629 -------------+---------------------------------------------------------------- gues | item | 1 | .1078598 .0337844 .0019 .1020673 .0581353 .1929404 2 | .1128113 .0372217 .002162 .1065996 .0596554 .2082417 3 | .123031 .0480042 .002579 .1127147 .0605462 .2516237 4 | .1190103 .0390721 .002369 .1123544 .0617698 .2095427 5 | .0829503 .0185785 .001275 .0807116 .0514752 .1232547 6 | .1059315 .0289175 .001708 .1022741 .0584959 .1709483 7 | .1235553 .0382661 .002964 .1186648 .0626495 .2067556 8 | .1142118 .0408348 .001733 .1062507 .0592389 .2134006 9 | .1270767 .0557821 .003939 .113562 .0621876 .2825752 -------------+---------------------------------------------------------------- mu_a | .109161 .1218499 .005504 .1126253 -.135329 .3501061 var_a | .108864 .0331522 .001053 .1030106 .0604834 .1860996 mu_b | .0782094 .1974657 .004367 .0755023 -.3067717 .4638104 var_b | .5829738 .1803167 .006263 .5562159 .3260449 1.034225 ------------------------------------------------------------------------------ file sim3pls.dta saved . estimates store est3pls

The estimated posterior means of the \(c_i\)’s range between 0.08 and 0.13. Clearly, the introduction of guessing parameters has an impact on the item discrimination and difficulty parameters. For example, the estimated posterior means of \(\mu_a\) and \(\mu_b\) shift from -0.10 and -0.07, respectively, for the 2PL model to 0.11 and 0.08, respectively, for the 3PL model.

Because the estimated guessing parameters are not that different, one may ask whether item-specific guessing parameters are really necessary. To answer this question, we fit a model with a common guessing parameter, **{gues}**, and compare it with the previous model.

. set seed 14 . bayesmh y, likelihood(dbernoulli({gues}+(1-{gues})* /// > invlogit({discrim:}*({subj:}-{diff:})))) /// > redefine(discrim:i.item) redefine(diff:i.item) /// > redefine(subj:i.id) /// > prior({subj:i.id}, normal(0, 1)) /// > prior({discrim:i.item}, lognormal({mu_a}, {var_a})) /// > prior({diff:i.item}, normal({mu_b}, {var_b})) /// > prior({gues}, igamma(10, 1)) /// > prior({mu_a} {mu_b}, normal(0, 0.1)) /// > prior({var_a} {var_b}, igamma(10, 1)) /// > block({mu_a mu_b var_a var_b gues}, split) /// > init({discrim:i.item} 1 {gues} 0.1) /// > exclude({subj:i.id}) burnin(5000) saving(sim3pl, replace) Burn-in ... Simulation ... Model summary ------------------------------------------------------------------------------ Likelihood: y ~ binomial({gues}+(1-{gues})*invlogit(xb_discrim*(xb_subj-xb_diff)),1) Priors: {discrim:i.item} ~ lognormal({mu_a},{var_a}) (1) {diff:i.item} ~ normal({mu_b},{var_b}) (2) {subj:i.id} ~ normal(0,1) (3) {gues} ~ igamma(10,1) Hyperpriors: {mu_a mu_b} ~ normal(0,0.1) {var_a var_b} ~ igamma(10,1) ------------------------------------------------------------------------------ (1) Parameters are elements of the linear form xb_discrim. (2) Parameters are elements of the linear form xb_diff. (3) Parameters are elements of the linear form xb_subj. Bayesian Bernoulli model MCMC iterations = 15,000 Random-walk Metropolis-Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 7,200 Acceptance rate = .3753 Efficiency: min = .01295 avg = .03714 Log marginal likelihood = . max = .1874 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval] -------------+---------------------------------------------------------------- discrim | item | 1 | 1.692894 .2748163 .021944 1.664569 1.232347 2.299125 2 | .8313512 .1355267 .00606 .8218212 .5928602 1.125729 3 | 1.058833 .1611742 .014163 1.054126 .7676045 1.393611 4 | 1.041808 .1718472 .008782 1.029867 .7398569 1.397073 5 | 1.534997 .3208687 .023965 1.497019 1.019998 2.266078 6 | 1.38296 .2581948 .019265 1.355706 .9559487 1.979358 7 | .8310222 .1698206 .012896 .8107371 .5736484 1.248736 8 | 1.442949 .2266268 .017562 1.431204 1.066646 1.930829 9 | .77944 .1159669 .007266 .7750891 .5657258 1.014941 -------------+---------------------------------------------------------------- diff | item | 1 | -.3043161 .0859905 .005373 -.2968324 -.4870583 -.1407109 2 | .1814508 .1289251 .006543 .1832146 -.0723988 .4313265 3 | -1.391216 .1924384 .014986 -1.373093 -1.809343 -1.050919 4 | .5928491 .1262631 .006721 .5829347 .356614 .857743 5 | 1.617348 .1929263 .011604 1.601534 1.293032 2.061096 6 | .817635 .1172884 .006125 .812838 .5990503 1.064322 7 | 2.006949 .2743517 .01785 1.981052 1.556682 2.594236 8 | -1.576235 .1747855 .013455 -1.559435 -1.952676 -1.272108 9 | -1.039362 .1840773 .01138 -1.02785 -1.432058 -.7160181 -------------+---------------------------------------------------------------- gues | .1027336 .0214544 .001753 .1022211 .0627299 .1466367 mu_a | .1009741 .123915 .006567 .0965353 -.1343028 .3510697 var_a | .1121003 .0344401 .001154 .1059563 .0628117 .1970842 mu_b | .0632173 .1979426 .004572 .0666684 -.3292497 .4482957 var_b | .5861236 .1818885 .006991 .5574743 .3239369 1.053172 ------------------------------------------------------------------------------ file sim3pl.dta saved . estimates store est3pl

We can again compare the two 3PL models by using the **bayesstats ic** command:

. bayesstats ic est3pls est3pl, diconly Deviance information criterion ------------------------ | DIC -------------+---------- est3pls | 8049.425 est3pl | 8049.426 ------------------------

Although the estimated DICs of the two 3PL models are essentially the same, we decide for demonstration purposes to proceed with the model with item-specific guessing parameters.

4PL model

The four-parameter logistic (4PL) model extends the 3PL model by adding item-specific upper asymptote parameters \(d_i\):

\[

P(Y_{ij}=1) = c_i + (d_i-c_i){\rm InvLogit}\{a_i(\theta_j-b_i)\}

,\ c_i < d_i < 1
\]
The \(d_i\) parameter can be viewed as an upper limit on the probability of correct response to the \(i\)th item. The probability of giving correct answers by subjects with very high ability can thus be no greater than \(d_i\).
We restrict the \(d_i\)'s to the (0.8,1) range and assign them a \({\rm Uniform}(0.8,1)\) prior. For other parameters, we use the same priors as in the 3PL model.
In the **bayesmh** specification of the model, the condition \(c_i < d_i\) is incorporated in the likelihood, and the condition \(d_i < 1\) is implied by the specified prior for the \(d_i\)'s. We initialize the \(d_i\)'s to 0.9. We use the **notable** option to suppress the long table output.

. set seed 14 . bayesmh y, likelihood(dbernoulli(({gues:}+({d:}-{gues:})* /// > invlogit({discrim:}*({subj:}-{diff:})))* /// > cond({gues:}<{d:},1,.))) /// > redefine(discrim:i.item) redefine(diff:i.item) /// > redefine(gues:i.item) redefine(d:i.item) redefine(subj:i.id) /// > prior({subj:i.id}, normal(0, 1)) /// > prior({discrim:i.item}, lognormal({mu_a}, {var_a})) /// > prior({diff:i.item}, normal({mu_b}, {var_b})) /// > prior({gues:i.item}, igamma(10, 1)) /// > prior({d:i.item}, uniform(0.8, 1)) /// > prior({mu_a} {mu_b}, normal(0, 0.1)) /// > prior({var_a} {var_b}, igamma(10, 1)) /// > block({mu_a mu_b var_a var_b}, split) /// > init({discrim:i.item} 1 {gues:i.item} 0.1 {d:i.item} 0.9) /// > exclude({subj:i.id}) burnin(5000) saving(sim4pls, replace) notable Burn-in ... Simulation ... Model summary ------------------------------------------------------------------------------ Likelihood: y ~ binomial(<expr1>,1) Priors: {discrim:i.item} ~ lognormal({mu_a},{var_a}) (1) {diff:i.item} ~ normal({mu_b},{var_b}) (2) {gues:i.item} ~ igamma(10,1) (3) {d:i.item} ~ uniform(0.8,1) (4) {subj:i.id} ~ normal(0,1) (5) Hyperpriors: {mu_a mu_b} ~ normal(0,0.1) {var_a var_b} ~ igamma(10,1) Expression: expr1 : (xb_gues+(xb_d-xb_gues)*invlogit(xb_discrim*(xb_subj-xb_diff)))* con d(xb_gues<xb_d,1,.) ------------------------------------------------------------------------------ (1) Parameters are elements of the linear form xb_discrim. (2) Parameters are elements of the linear form xb_diff. (3) Parameters are elements of the linear form xb_gues. (4) Parameters are elements of the linear form xb_d. (5) Parameters are elements of the linear form xb_subj. Bayesian Bernoulli model MCMC iterations = 15,000 Random-walk Metropolis-Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 7,200 Acceptance rate = .3639 Efficiency: min = .004825 avg = .0219 Log marginal likelihood = . max = .1203 Note: There is a high autocorrelation after 500 lags. file sim4pls.dta saved . estimates store est4pls

We use **bayesstats summary** to display results of selected model parameters.

. bayesstats summary {d:i.item} {mu_a var_a mu_b var_b} Posterior summary statistics MCMC sample size = 10,000 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval] -------------+---------------------------------------------------------------- d | item | 1 | .9598183 .0255321 .001948 .9621874 .9044441 .9981723 2 | .9024564 .0565702 .007407 .9019505 .8066354 .9944216 3 | .9525519 .0281878 .002845 .9551054 .8972454 .9971564 4 | .8887963 .0561697 .005793 .8859503 .8036236 .9916784 5 | .8815547 .0588907 .007215 .8708021 .8031737 .9926549 6 | .8891188 .0586482 .006891 .881882 .8024593 .9935512 7 | .874271 .0561718 .008087 .8635082 .8018176 .9880433 8 | .9663644 .0147606 .001121 .9667563 .9370666 .9950912 9 | .889164 .0486038 .005524 .8834207 .8084921 .9857415 -------------+---------------------------------------------------------------- mu_a | .3336887 .1436216 .009742 .334092 .0562924 .6164115 var_a | .1221547 .0406908 .002376 .1144729 .0642768 .2229326 mu_b | -.0407488 .1958039 .005645 -.0398847 -.4220523 .3323791 var_b | .4991736 .1612246 .00629 .4660071 .2802531 .9023824 ------------------------------------------------------------------------------

The **bayesmh** command issued a note indicating high autocorrelation for some of the model parameters. This may be related to slower MCMC convergence or more substantial problems in the model specification. It is thus worthwhile to inspect the individual autocorrelation of the parameters. We can do so by using the **bayesstats ess** command. The parameters with lower estimated sample size (ESS) have higher autocorrelation and vice versa.

. bayesstats ess {d:i.item} {mu_a var_a mu_b var_b} Efficiency summaries MCMC sample size = 10,000 ---------------------------------------------------- | ESS Corr. time Efficiency -------------+-------------------------------------- d | item | 1 | 171.82 58.20 0.0172 2 | 58.33 171.43 0.0058 3 | 98.17 101.87 0.0098 4 | 94.02 106.36 0.0094 5 | 66.62 150.11 0.0067 6 | 72.44 138.05 0.0072 7 | 48.25 207.26 0.0048 8 | 173.30 57.70 0.0173 9 | 77.41 129.19 0.0077 -------------+-------------------------------------- mu_a | 217.35 46.01 0.0217 var_a | 293.34 34.09 0.0293 mu_b | 1203.20 8.31 0.1203 var_b | 656.92 15.22 0.0657 ----------------------------------------------------

We observe that the parameters with ESS lower than 200 are among the asymptote parameter’s \(d_i\)’s. This may be caused, for example, by overparameterization of the likelihood model and subsequent nonidentifiability, which is not resolved by the specified priors.

We can also fit a model with a common upper asymptote parameter, \(d\), and compare it with the model with the item-specific upper asymptote.

. set seed 14 . bayesmh y, likelihood(dbernoulli(({gues:}+({d}-{gues:})* /// > invlogit({discrim:}*({subj:}-{diff:})))* /// > cond({gues:}<{d},1,.))) /// > redefine(discrim:i.item) redefine(diff:i.item) /// > redefine(gues:i.item) redefine(subj:i.id) /// > prior({subj:i.id}, normal(0, 1)) /// > prior({discrim:i.item}, lognormal({mu_a}, {var_a})) /// > prior({diff:i.item}, normal({mu_b}, {var_b})) /// > prior({gues:i.item}, igamma(10, 1)) /// > prior({d}, uniform(0.8, 1)) /// > prior({mu_a} {mu_b}, normal(0, 0.1)) /// > prior({var_a} {var_b}, igamma(10, 1)) /// > block({mu_a mu_b var_a var_b d}, split) /// > init({discrim:i.item} 1 {gues:i.item} 0.1 {d} 0.9) /// > exclude({subj:i.id}) burnin(5000) saving(sim4pl, replace) notable Burn-in ... Simulation ... Model summary ------------------------------------------------------------------------------ Likelihood: y ~ binomial(<expr>>,1) Priors: {discrim:i.item} ~ lognormal({mu_a},{var_a}) (1) {diff:i.item} ~ normal({mu_b},{var_b}) (2) {gues:i.item} ~ igamma(10,1) (3) {subj:i.id} ~ normal(0,1) (4) {d} ~ uniform(0.8,1) Hyperpriors: {mu_a mu_b} ~ normal(0,0.1) {var_a var_b} ~ igamma(10,1) Expression: expr1 : (xb_gues+({d}-xb_gues)*invlogit(xb_discrim*(xb_subj-xb_diff)))* cond (xb_gues<{d},1,.) ------------------------------------------------------------------------------ (1) Parameters are elements of the linear form xb_discrim. (2) Parameters are elements of the linear form xb_diff. (3) Parameters are elements of the linear form xb_gues. (4) Parameters are elements of the linear form xb_subj. Bayesian Bernoulli model MCMC iterations = 15,000 Random-walk Metropolis-Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 7,200 Acceptance rate = .3877 Efficiency: min = .0107 avg = .03047 Log marginal likelihood = . max = .1626 file sim4pl.dta saved . estimates store est4pl . bayesstats summary {d mu_a var_a mu_b var_b} Posterior summary statistics MCMC sample size = 10,000 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval] -------------+---------------------------------------------------------------- d | .9664578 .0144952 .001293 .9668207 .9371181 .9924572 mu_a | .2206696 .1387873 .01113 .2208302 -.0483587 .4952625 var_a | .1245785 .0391551 .001806 .1188779 .0658243 .2187058 mu_b | .0371722 .2020157 .00501 .0331742 -.3481366 .4336587 var_b | .5603447 .1761812 .006817 .5279243 .3157048 .9805077 ------------------------------------------------------------------------------

We now compare the two 4PL models by using the **bayesstats ic** command:

. bayesstats ic est4pls est4pl, diconly Deviance information criterion ------------------------ | DIC -------------+---------- est4pls | 8050.805 est4pl | 8037.075 ------------------------

The DIC of the more complex 4PL model (8,051) is substantially higher than the DIC of the simpler model (8,037). This and the potential nonidentifiability of the more complex **est4pls** model, indicated by high autocorrelation in the simulated MCMC sample, compel us to proceed with the model with a common upper asymptote, **est4pl**.

The posterior distribution of \(d\) has an estimated 95% equal-tailed credible interval of (0.93, 0.99) and is concentrated about 0.97. The \({\rm Uniform}(0.8,1)\) prior on \(d\) does not seem to be too restrictive. The estimated DIC of the **est4pl** model (8,037) is lower than the DIC of the **est3pls** 3PL model from the previous section (8,049), implying that the introduction of the upper asymptote parameter \(d\) does improve the model fit.

5PL model

The five-parameter logistic (5PL) model extends the 4PL model by adding item-specific asymmetry parameters \(e_i\):

\[

P(Y_{ij}=1) = c_i + (d_i-c_i){\rm InvLogit}\big[{\{a_i(\theta_j-b_i)\}}^{e_i}\big]

,\ c_i < d_i < 1, \ 0 < e_i < 1
\]
In the previous section, we found the 4PL model with common upper asymptote \(d\), **est4pl**, to be the best one so far. We thus consider here a 5PL model with common upper asymptote \(d\).

Typically, we expect the \(e_i\) parameters to be close to 1. Similarly to the upper asymptote parameter \(d\), the \(e_i\) parameters are assumed to be in the (0.8,1) range and are assigned \({\rm Uniform}(0.8,1)\) prior. We initialize the \(e_i\)s to 0.9. We again use the **notable** option to suppress the long table output, and we display a subset of results by using **bayesstats summary**. (We could have used **bayesmh**'s **noshow()** option instead to achieve the same result.)

. set seed 14 . bayesmh y, likelihood(dbernoulli(({gues:}+({d}-{gues:})* /// > (invlogit({discrim:}*({subj:}-{diff:})))^{e:})* /// > cond({gues:}<{d},1,.))) /// > redefine(discrim:i.item) redefine(diff:i.item) /// > redefine(gues:i.item) redefine(e:i.item) redefine(subj:i.id) /// > prior({subj:i.id}, normal(0, 1)) /// > prior({discrim:i.item}, lognormal({mu_a}, {var_a})) /// > prior({diff:i.item}, normal({mu_b}, {var_b})) /// > prior({gues:i.item}, igamma(10, 1)) /// > prior({d}, uniform(0.8, 1)) /// > prior({e:i.item}, uniform(0.8, 1)) /// > prior({mu_a} {mu_b}, normal(0, 0.1)) /// > prior({var_a} {var_b}, igamma(10, 1)) /// > block({mu_a mu_b var_a var_b d}, split) /// > init({discrim:i.item} 1 {gues:i.item} 0.1 {d} {e:i.item} 0.9) /// > exclude({subj:i.id}) burnin(5000) saving(sim5pls, replace) notable Burn-in ... Simulation ... Model summary ------------------------------------------------------------------------------ Likelihood: y ~ binomial(<expr1>,1) Priors: {discrim:i.item} ~ lognormal({mu_a},{var_a}) (1) {diff:i.item} ~ normal({mu_b},{var_b}) (2) {gues:i.item} ~ igamma(10,1) (3) {e:i.item} ~ uniform(0.8,1) (4) {subj:i.id} ~ normal(0,1) (5) {d} ~ uniform(0.8,1) Hyperpriors: {mu_a mu_b} ~ normal(0,0.1) {var_a var_b} ~ igamma(10,1) Expression: expr1 : (xb_gues+({d}-xb_gues)*(invlogit(xb_discrim*(xb_subj-xb_diff)))^xb_e )* cond(xb_gues<{d},1,.) ------------------------------------------------------------------------------ (1) Parameters are elements of the linear form xb_discrim. (2) Parameters are elements of the linear form xb_diff. (3) Parameters are elements of the linear form xb_gues. (4) Parameters are elements of the linear form xb_e. (5) Parameters are elements of the linear form xb_subj. Bayesian Bernoulli model MCMC iterations = 15,000 Random-walk Metropolis-Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 7,200 Acceptance rate = .3708 Efficiency: min = .007341 avg = .02526 Log marginal likelihood = . max = .1517 file sim5pls.dta saved . estimates store est5pls . bayesstats summary {e:i.item} {d mu_a var_a mu_b var_b} Posterior summary statistics MCMC sample size = 10,000 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval] -------------+---------------------------------------------------------------- e | item | 1 | .897859 .0578428 .006083 .8939272 .8050315 .9957951 2 | .9042669 .0585023 .005822 .90525 .8053789 .9956565 3 | .88993 .0562398 .005013 .887011 .803389 .9930454 4 | .9010241 .0574186 .006492 .9042044 .8030981 .9925598 5 | .9126369 .0545625 .00521 .9178927 .8098596 .9964487 6 | .9037269 .0583833 .006814 .9086704 .8054932 .9961268 7 | .9136308 .0558911 .005373 .9203899 .8112029 .996217 8 | .889775 .0568656 .005119 .8849938 .803912 .9938777 9 | .8808435 .056257 .004743 .8727194 .8030522 .9904972 -------------+---------------------------------------------------------------- d | .9671374 .0144004 .001165 .9670598 .9382404 .9933374 mu_a | .2770211 .1353777 .00832 .2782552 .0141125 .5418087 var_a | .122635 .0404159 .002148 .1160322 .0666951 .2208711 mu_b | .1211885 .1929743 .004955 .1199136 -.2515431 .503733 var_b | .5407642 .1747674 .006353 .5088269 .3016315 .9590086 ------------------------------------------------------------------------------

We also want to compare the above model with a simpler one using a common asymmetry parameter \(e\).

. set seed 14 . bayesmh y, likelihood(dbernoulli(({gues:}+({d}-{gues:})* /// > (invlogit({discrim:}*({subj:}-{diff:})))^{e})* /// > cond({gues:}<{d},1,.))) /// > redefine(discrim:i.item) redefine(diff:i.item) /// > redefine(gues:i.item) redefine(subj:i.id) /// > prior({subj:i.id}, normal(0, 1)) /// > prior({discrim:i.item}, lognormal({mu_a}, {var_a})) /// > prior({diff:i.item}, normal({mu_b}, {var_b})) /// > prior({gues:i.item}, igamma(10, 1)) /// > prior({d} {e}, uniform(0.8, 1)) /// > prior({mu_a} {mu_b}, normal(0, 0.1)) /// > prior({var_a} {var_b}, igamma(10, 1)) /// > block({mu_a mu_b var_a var_b d e}, split) /// > init({discrim:i.item} 1 {gues:i.item} 0.1 {d e} 0.9) /// > exclude({subj:i.id}) burnin(5000) saving(sim5pl, replace) notable Burn-in ... Simulation ... Model summary ------------------------------------------------------------------------------ Likelihood: y ~ binomial(<expr1>,1) Priors: {discrim:i.item} ~ lognormal({mu_a},{var_a}) (1) {diff:i.item} ~ normal({mu_b},{var_b}) (2) {gues:i.item} ~ igamma(10,1) (3) {subj:i.id} ~ normal(0,1) (4) {d e} ~ uniform(0.8,1) Hyperpriors: {mu_a mu_b} ~ normal(0,0.1) {var_a var_b} ~ igamma(10,1) Expression: expr1 : (xb_gues+({d}-xb_gues)*(invlogit(xb_discrim*(xb_subj-xb_diff)))^{e}) * cond(xb_gues<{d},1,.) ------------------------------------------------------------------------------ (1) Parameters are elements of the linear form xb_discrim. (2) Parameters are elements of the linear form xb_diff. (3) Parameters are elements of the linear form xb_gues. (4) Parameters are elements of the linear form xb_subj. Bayesian Bernoulli model MCMC iterations = 15,000 Random-walk Metropolis-Hastings sampling Burn-in = 5,000 MCMC sample size = 10,000 Number of obs = 7,200 Acceptance rate = .3805 Efficiency: min = .008179 avg = .02768 Log marginal likelihood = . max = .08904 file sim5pl.dta saved . estimates store est5pl . bayesstats summary {e d mu_a var_a mu_b var_b} Posterior summary statistics MCMC sample size = 10,000 ------------------------------------------------------------------------------ | Equal-tailed | Mean Std. Dev. MCSE Median [95% Cred. Interval] -------------+---------------------------------------------------------------- e | .9118363 .0558178 .004194 .9175841 .8063153 .9960286 d | .9655166 .0147373 .001495 .9659029 .9354708 .9924492 mu_a | .2674271 .1368926 .008485 .270597 .0102798 .5443345 var_a | .1250759 .0428095 .002635 .1173619 .0654135 .2340525 mu_b | .1015121 .2048178 .006864 .103268 -.3052377 .4934158 var_b | .5677309 .1824591 .006981 .5331636 .3079868 1.016762 ------------------------------------------------------------------------------

We use **bayesstats ic** to compare the DIC values of the two 5PL models:

. bayesstats ic est5pls est5pl, diconly Deviance information criterion ------------------------ | DIC -------------+---------- est5pls | 8030.894 est5pl | 8034.517 ------------------------

The estimated DIC of the more complex **est5pls** model (8,031) is lower than the DIC of the simpler model (8,035), suggesting a better fit.

Conclusion

Finally, we compare all eight fitted models.

. bayesstats ic est1pl est2pl est3pl est3pls est4pl est4pls est5pl est5pls, /// > diconly Deviance information criterion ------------------------ | DIC -------------+---------- est1pl | 8122.428 est2pl | 8055.005 est3pl | 8049.426 est3pls | 8049.425 est4pl | 8037.075 est4pls | 8050.805 est5pl | 8034.517 est5pls | 8030.894 ------------------------

The **est5pls** model has the lowest overall DIC. To confirm this result, we run another set of simulations with a larger MCMC sample size of 50,000. (We simply added the **mcmcsize(50000)** option to the **bayesmh** specification of the above eight models.) The following DIC values, based on the larger MCMC sample size, are more reliably estimated.

. bayesstats ic est1pl est2pl est3pl est3pls est4pl est4pls est5pl est5pls, /// > diconly Deviance information criterion ------------------------ | DIC -------------+---------- est1pl | 8124.015 est2pl | 8052.068 est3pl | 8047.067 est3pls | 8047.738 est4pl | 8032.417 est4pls | 8049.712 est5pl | 8031.375 est5pls | 8031.905 ------------------------

Again, the 5PL models have the lowest DIC values and seem to provide the best fit. However, the DIC differences between models **est4pl**, **est5pl**, and **est5pls** are minimal and may very well be within the estimation error. Regardless, these three models appear to be better than the simpler 1PL, 2PL, and 3PL models.

Additional model checking may be needed to assess the models' fit, and we should not rely solely on the DIC values to make our final model selection. A practitioner may still prefer the simpler **est4pl** 4PL model to the 5PL models even though it has a slightly higher DIC. In fact, given that the posterior mean estimate of the upper asymptote parameter \(d\) is 0.96 with a 95% equal-tailed credible interval of (0.94, 0.99), some practitioners may prefer the even simpler **est3pl** 3PL model.

**References**

De Boeck, P., and M. Wilson, ed. 2004. *Explanatory Item Response Models: A Generalized Linear and Nonlinear Approach*. New York: Springer.

Kim, J.-S., and D. M. Bolt. 2007. *Estimating item response theory models using Markov chain Monte Carlo methods*. Educational Measurement: Issues and Practice 26: 38-51.

Pingback: The Stata Blog » Fitting distributions using bayesmh()