Post 8: Sampling the variance of person ability with a Gibbs step

[This article was first published on Markov Chain Monte Carlo in Item Response Models, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 text{rest} stands in for the conditioning variables, f_star(dagger|dots) represents the density of the random variable named dagger (which has a star distribution), and pi_{pi} 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 sigma^2_theta only depends on the data (boldsymbol U) via the values of boldsymbol theta. This is common for parameters that are at higher levels of model hierarchy.

Since the complete conditional for sigma^2_theta 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 

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 

We can check burnin by visualizing the chain with coda:

run.D.mcmc.all
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 

and then make a parameter recovery plot for the variance:

run.D.sig2.param.recovery
Where we first define a new helper function check.sigma as

check.sigma 

Then we call check.sigma for this example to make the graph:

check.sigma( run.D.mcmc.conv, c(1.2,2))

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

  1. The true value is contained within the credible interval of the estimate.
  2. The MCMC estimator for sigma^2_theta is (roughly) an estimate of the sample variance of the theta_p 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 theta_p 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:

run.D.theta.param.recovery
## Calculate the EAP estimates
all.eap 

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.

To 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 about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)