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

October 10, 2013
By

(This article was first published on 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 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 <- 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:

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 <- window( run.D.mcmc, start=100)

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 <- 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 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 <- 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.

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 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...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

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)