**Markov Chain Monte Carlo in Item Response Models**, and kindly contributed to R-bloggers)

The Bayesian 2PL IRT model we defined in Post 1 set a hyper-prior on the variance of the person ability parameters. This post implements the sampler for this parameter as a Gibbs step. We will check that Gibbs step is working by running it on the fake data from Post 2, visualizing the results, and checking that the true value is recovered.

# Implementing the variance sampler

Recall from Post 1 that the relevant complete conditional densities are:

begin{aligned}

f(theta_p|text{rest}) & propto

prod_{i=1}^I

pi_{pi}^{u_{pi}} (1 – pi_{pi})^{1-u_{pi}}

f_text{Normal}(theta_p| 0,sigma_theta^2) ~, \

f(sigma^2_theta|text{rest}) & propto

f_text{Inverse-Gamma}left(sigma^2_theta left|

alpha_theta + frac{P}{2},

beta_theta + frac{1}{2} sum_{p=1}^P theta^2_p right.

right) ,

end{aligned}

where stands in for the conditioning variables, represents the density of the random variable named (which has a distribution), and is defined as:

begin{equation}

ln{frac{pi_{pi}}{1+pi_{pi}}} = a_i ( theta_p – b_i) label{eq:pipi} quad.

end{equation}

Note that the complete conditional for only depends on the data () via the values of . This is common for parameters that are at higher levels of model hierarchy.

Since the complete conditional for is distributed as an inverse-gamma, we can sample from the complete conditional directly without the need for the Metropolis-Hastings algorithm. This is referred to as a Gibbs step.

We implement the sampler in R by noting that an Inverse-Gamma is defined to be the reciprocal of a Gamma random variable, and that R implements a draw from a Gamma with the `rgamma`

function:

sample.s2 <- function(U.data, old) { ## Grab the appropriate values alpha.th <- old$hyperpar$alpha.th beta.th <- old$hyperpar$beta.th P.persons <- length(old$th) ## Update the chain cur <- old cur$s2 <- 1/rgamma(1, shape = alpha.th + P.persons/2, rate = beta.th + sum((old$th)^2)/2 ) return(cur) }

Note that we do not need to calculate an acceptance rate because it is always 1.

# Testing that the sampler works

Since the complete conditional for the variance depends only on the person ability parameters, we simplify our test of the variance sampler code by setting the item parameter samplers to the dummy sampler from Post 3.

## Load the necessary code from this and previous posts source("http://mcmcinirt.stat.cmu.edu/setup/post-8.R") ## Reset the item sampler to dummy samplers sample.a <- dummy.cc.sampler sample.b <- dummy.cc.sampler set.seed( 314159 ) run.D <- run.chain.2pl( M.burnin = 0, M.keep = 1000, M.thin = 1, U.data = U, hyperpars = hyperpars, ## Generate starting values uniformly from -5 to 5 th.init = runif( length(theta.abl), min=-5, max=5 ), ## Keep item parameters at their true values a.init = a.disc, b.init = b.diff, ## Generate starting values uniformly from 0 to 5 s2.init = runif( 1, min=0, max=5), ## MH.th=1, MH.a=NA, MH.b=NA, verbose = TRUE ) ## Average acceptance rates: ## theta.abl: 0.482 ## a.disc: NaN ## b.diff: NaN ## Convert to coda run.D.mcmc <- mcmc( t(run.D) )

We can check burnin by visualizing the chain with `coda`

:

plot( run.D.mcmc[, get.2pl.params(1,1,1,1)], density=FALSE, smooth=TRUE ) |

It looks like the sampler burns in almost immediately. We discard the first 100 iterations to be conservative

run.D.mcmc.conv <- window( run.D.mcmc, start=100)

and then make a parameter recovery plot for the variance:

The MCMC estimate is off from the true value of , which is plotted in blue. This is not a concern for two reasons:

- The true value is contained within the credible interval of the estimate.
- The MCMC estimator for is (roughly) an estimate of the
*sample variance*of the values (see Post 1 for details). Note that our MCMC estimate of the posterior mode is quite close to the sample variance of the generated values, which is plotted in purple. Indeed, if we generate*different data*, the posterior mode will find the*sample variance*of the person ability parameters each time (not shown).

As an additional check that the Gibbs step is working, we can make sure that the person ability parameters are still recovered well:

The person ability parameters are, indeed, still recovered well.

# Conclusion

The sampler for the variance of the person parameters is working well.

In the next post we combine all of the samplers and tune them.

**leave a comment**for the author, please follow the link and comment on their blog:

**Markov Chain Monte Carlo in Item Response Models**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...