In the IRT literature, it is well accepted that when prior information is not readily available, large data sizes are needed to estimate IRT parameters (e.g., Swaminathan and Gifford, 1983), as Bayesian estimation with a flat prior is equivalent to a maximum likelihood estimation, and that the accuracy of item (person) parameter estimates is related to the number of subjects (items) (e.g., Natesan et al., 2016; Sheng, 2010). Hence, in this simulation study, we evaluated the effect of sample size, test length, actual variability of item slope and intercept parameters on the accuracy with which model parameters are estimated considering a half-normal, a non-informative uniform or an inverse gamma prior.

Item responses for k items (k = 10, 20, and 40) and n individuals (n = 100, 300, 500, and 1000) were generated according to the 3PNO model, as defined in (1). Ability parameters were generated as samples from a standard normal distribution, pseudo-chance-level parameters were generated from a uniform distribution, γ_j ~ U(0.05, 0.4), and item slope and intercept parameters were generated as samples from uniform distributions so that

Sim1: α_j ~ U(0, 2), β_j ~ U(−1, 1);

When implementing the MCMC procedure, a diffuse prior were assumed for γ_j, μ_α and μ_β so that γ_j ~ Beta(1, 1), p(μ_α)∝1 and p(μ_β)∝1. In addition, three ways of setting the prior distributions for α_j and β_j were considered such that both had

non-informative prior distributions that are uniform on σ, i.e., p(σα2)∝1/σα and p(σβ2)∝1/σβ ;

It is noted that although the inverse-gamma (0.001, 0.001) prior density was not suggested by Gelman (2006), it was considered in this study because of its popularity (see e.g., Xu et al., 2009; O'Brien et al., 2015). With each of the prior specifications considered, the Gibbs sampling procedure was implemented where 10,000–50,000 iterations were obtained with the first half set as burn-in.

In the first simulation study, only two levels of variability (i.e., σα,β2=0.083 and σα,β2=0.333), and three different hyperpriors (one uniform, one inverse-gamma and one half-normal) were considered. Given that the range for item intercept parameters (β_j) is generally wider than [−1, 1] in practice, it would be interesting to see how the weakly informative half-t family performs when the variance for β goes beyond 1. Hence, a second simulation study was conducted where two factors were manipulated, namely, actual variability of the item intercept parameters and specifications of prior distributions for the item scale hyperparameters.

Item responses for 20 items and 1000 individuals were generated according to the 3PNO model, as defined in (1). Ability parameters were generated as samples from a standard normal distribution, item slope and pseudo-chance-level parameters were generated from uniform distributions, α_j ~ U(0, 2) and γ_j ~ U(0.05, 0.4), and item intercept parameters were generated as samples from uniform distributions so that

Sim1: β ~ U(−0.5, 0.5);

It is noted that in the four simulations, the uniform distributions from which the intercept parameters were sampled from have an increasingly large σ_β ranging from 0.29−2.31 (with the corresponding variance ranging from 0.083−5.333).

In addition, six ways of setting the prior distributions for σα2 and σβ2 were considered such that both had

non-informative prior distributions that are uniform on σ, i.e., p(σα2)∝1/σα and p(σβ2)∝1/σβ;

The Gibbs sampling procedure was implemented where 20,000 iterations were obtained with the first 10,000 as burn-in.

For both simulation studies, convergence was evaluated using the Gelman-Rubin R (Gelman and Rubin, 1992) statistic. The usual practice is using multiple Markov chains from different starting points. Alternatively, a single chain can be divided into sub-chains so that convergence is assessed by comparing the between and within sub-chain variance (Fox, 2007). Since a single chain is less wasteful in the number of iterations needed, the latter approach was adopted. For each Markov chain, the initial values were set to be α_j = 1, β_j = 0, and γ_j = 0.2 for all items j and θ_i = 0 for all persons i. After discarding the burn-in samples, the chain was then separated into five sub-chains of equal length and the R statistic was calculated following the procedure by Gelman and Rubin (1992). Convergence can also be monitored visually using time series graphs of the simulated sequence, such as the trace plot, the running mean plot, and the autocorrelation plot shown in Figure 1 for one item. It is observed that for this item, the autocorrelations between successive parameter draws became negligible at lags > 600, suggesting burn-in for a single chain should not take longer than that (Geyer, 1992). Indeed, the trace plot and the running mean plot both suggest possible convergence within 10,000 iterations. Inspection of such plots has, however, been criticized for being unreliable and unwieldy in the presence of a large number of model parameters (Gelman et al., 2003; Nylander et al., 2008). The R statistic obtained from using a single chain was hence the major approach for assessing convergence in this study.

For each simulated scenario, 25 replications, as recommended by Harwell et al. (1996), were conducted to avoid erroneous results in estimation due to sampling error. The accuracy of item/person parameter estimates was evaluated using the root mean square error (RMSE) and bias. Let τ denote the true value of a parameter (e.g., α_j, β_j, γ_j, or θ_i) and t_r its estimate in the rth replication (r = 1, …, R). The RMSE is defined as

and the bias is defined as

These quantities were averaged over items/persons to provide summary indices.

In simulation study 2 where six model specifications were compared, the adequacy of the fit of the hierarchical 3PNO model with a given prior density on the simulated data was evaluated using Bayesian deviance. It should be noted that this measure provides a model comparison criterion. Hence, it evaluates the fit of a model in a relative, not absolute, sense. The Bayesian deviance information criterion (DIC) was introduced by Spiegelhalter et al. (2002) who generalized the classical information criteria (e.g., AIC, BIC) to one that is based on the posterior distribution of the deviance. This criterion is defined as DIC=D̄+pD, where D̄=E[-2logLϑ|y(y|ϑ)] is the posterior expectation of the deviance (with L being the likelihood function), and pD=Eϑ|y(D)-D[Eϑ|y(ϑ)]=D̄-D(ϑ̄) is the effective number of parameters (Carlin and Louis, 2000). In addition, let D(ϑ̄)=-2log[Lϑ|y(y|ϑ̄)], where ϑ̄ is the posterior mean. To compute Bayesian DIC, MCMC samples of the parameters, ϑ⁽¹⁾, …, ϑ^(G), can be drawn from the Gibbs sampler, then D̄ is approximated as D̄=1G[-2log∏g = 1GL(y|ϑ(g))]. Generally more complicated models tend to provide better fit. Hence, penalizing for the number of parameters makes DIC a more reasonable measure to use. In this study, the Bayesian deviance estimate was obtained after each implementation and averaged across the 25 replications for each model specification.

Each implementation of the Gibbs sampler gave rise to Gelman-Rubin R statistics close to 1, indicating possible convergence of the Markov chains within the simulated number of iterations. Hence, the posterior estimates were obtained as the posterior expectations of the Gibbs samples and the average RMSE and bias values for α_j, β_j, γ_j, and θ_i are summarized in Tables 1–4. A close examination of these values leads to the following observations:

In Sim1 where the actual variances for α and β were both 0.333, the half-normal distribution performed relatively less well in recovering item and person parameters than the uniform or inverse-gamma distributions when data sizes (n × k) were < 10,000, although it has to be noted that this distribution consistently resulted in smaller bias in estimating α regardless of sample size or test length. The uniform prior density performed similarly as the inverse-gamma prior density, with a slight advantage for data with small sample sizes (i.e., n = 100) especially where k = 10.

From these observations, it is hence noted that the half-normal prior outperforms the uniform or inverse-gamma family of prior distributions in estimating the corresponding parameters in the hierarchical 3PNO model when the variability for item slope or intercept parameters is close to zero. Further, the inverse-gamma (0.001, 0.001) prior tends to perform better in parameter recovery than the uniform prior when the respective item parameters have a small variance. This can be explained by the fact that the non-informative uniform prior density is flat and hence does not restrict σ away from large values, whereas inverse-gamma (0.001, 0.001) and half-normal distributions have a peak around zero and do perform such shrinkage for variances near zero. It is further noted that the estimation error and bias in estimating item parameters (α, β, or γ) reduce with the increase of sample sizes, and that the error and bias in estimating person parameters (θ) reduce with the increase of test lengths. This is consistent with findings from other studies (e.g., Swaminathan and Gifford, 1983; Sheng, 2010; Natesan et al., 2016).

Each implementation of the Gibbs sampler gave rise to Gelman-Rubin R statistics close to 1, indicating possible convergence within 20,000 iterations. Hence, the posterior estimates were obtained as the posterior expectations of the Gibbs samples and the average RMSE and bias values for α_j, β_j, and γ_j are summarized in Table 5. A close examination of these values leads to the following observations:

In all four simulations where the actual variability for the intercept parameters (σ_β) ranged between 0.29−2.31, the three weakly informative half-t prior densities, namely, the half-Cauchy, the half-normal and the half-t with 4 degrees of freedom had consistently smaller RMSE if not bias in recovering item and person parameters, and in particular in recovering β_j, compared with the other three prior distributions.

From these observations, it is noted that the weakly informative half-t family works well in a wider range of situations than the non-informative or the informative prior densities in recovering the 3PNO model item parameters. Specifically, when the actual scale hyperparameter is close to 0, the non-informative prior density does not work well compared with informative or weakly informative prior densities, as prior information helps to restrict σ_β away from large values. Even with prior information, the informative inverse-gamma (3,2) prior density does not outperform the weakly informative half-t distribution because the latter has a better behavior near 0. Figure 2 graphically illustrates this point. The non-informative inverse-gamma prior distribution, in general, is not recommended for large σ values because it is not non-informative for these values (see Figure 2). On the other hand, prior information has to be correctly specified, as misspecification leads to large bias and hence large estimation error, e.g., the actual σ_β in Sim4 was clearly not in the range of the inverse-gamma (3, 2) distribution (see Figure 2). When comparing among the three half-t distributions, the half-Cauchy is more likely to allow for occasionally large values than the half-t density with 4 degrees of freedom, or the half-normal, which has the smallest tail (see Figure 2). Hence, if it is known a priori that the scale hyperparameter might be far away from 0, a half-Cauchy or a half-t distribution with a small degrees of freedom is suggested. Otherwise, a half-normal may be considered.

In addition to parameter recovery, the model comparison results in each simulation were averaged over the 25 replications and are summarized in Table 6, which shows the averaged estimates for the posterior expectation of the deviance (D̄), the deviance of the posterior expectation [D(ϑ̄)] values, the effective number of parameters (p_D), and the Bayesian DIC, respectively. The model with the half-Cauchy or half-t prior density shows consistently smaller D̄, D(ϑ̄), and DIC than those with other prior specifications. Since small deviance values indicate better model fit, models with such prior distributions for the scale hyperparameters are shown to provide a better description of the simulated data when the actual σ_β was especially larger than 0.29, even after penalizing for model complexities, i.e., the effective number of parameters.

After a close examination and comparison of the values shown in the table, a few remarks can be drawn from these results:

In all four simulations where σ_β ranged between 0.29−2.31, Bayesian deviances consistently preferred the models with half-t prior densities, whose deviance values were increasingly smaller compared with other model specifications when σ_β became larger.

In summary, the results in parameter recovery and model comparisons using Bayesian deviances suggest that the hierarchical 3PNO model with weakly informative half-t densities for the item scale hyperparameters does show advantages compared with those with non-informative or informative prior distributions in various test situations where the actual variability ranges from small to large values. When the actual variability of item parameters is not close to 0, the half-Cauchy or the half-t with small degrees of freedom is preferred over the half-normal distribution because of its flexibility in allowing for occasional large variability.

The overall CBASE exam contains 25 multiple-choice items on English reading/literature. The data used in this study were from college students who took the LP form of CBASE in years 2001 and 2002. After removing those who attempted the exam multiple times and removing missing responses, a sample of 1,200 examinees was randomly selected. To assess the model goodness-of-fit, hierarchical 3PNO models with the three half-t prior distributions were compared with each other and further compared with those with uniform and inverse-gamma prior densities described in Section 4.

Each of the six model specifications was implemented to the CBASE English data using the Gibbs sampling procedure, where 20,000 iterations were obtained with the first 10,000 set as burn-in. The Gelman-Rubin R statistics were used to assess convergence and they were found to be around or close to 1, suggesting that stationarity had been reached within the simulated Markov chains for the model. The Bayesian deviance estimates were subsequently obtained for each model specification and the results are summarized in Table 7. Among the six model specifications considered, the ones with half-t and half-Cauchy prior densities had the smallest DIC and/or expected posterior deviance (D̄) values. Therefore, the 3PNO model with a half-t or half-Cauchy hyperprior provided the best description of the data compared with other model specifications, even after penalizing for a large effective number of parameters (e.g., p_D = 987.24 for the half-t, and p_D = 993.80 for the half-Cauchy). On the other hand, with slightly larger deviances, the model with an informative inverse-gamma (3, 2) hyperprior described the data less adequately than those with other prior specifications. The relatively small differences in deviances suggested the actual variability might be close to 0. Furthermore, the model with a half-normal or half-t hyperprior had a larger p_D than that with a non-informative prior density. Given the findings from the second simulation study in Section 4, this indicated that the actual standard deviation for the item slope and/or intercept parameters with the data was not larger than 1.

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.