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)
}

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)

Where we first define a new helper function check.sigma as check.sigma <- function( mcmc.conv , xylim) { par(mfrow=c(1,2)) traceplot(mcmc.conv[, 'sig2.theta'], smooth=TRUE, ylim=xylim, main='Trace plot sig2.theta') abline(h=sig2.theta , col='blue') densplot( mcmc.conv[, 'sig2.theta'], main='Density plot sig2.theta', xlab='sig2.theta', xlim=xylim,) abline(v=sig2.theta , col='blue') abline(v=var(theta.abl), col='purple') legend( "topright", c('true value','var(theta.abl)'), col=c('blue', 'purple'), lty='solid') } 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 , 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 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:

## Calculate the EAP estimates all.eap <- apply( run.D.mcmc.conv, MARGIN=2, mean ) ## Visualize check.sampler.graph( theta.abl, all.eap[ get.2pl.params(1:P.persons,NULL,NULL,NULL)], desc="Person Ability", ylab="EAP Estimates", col="blue" )

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.