IRT model with random items


Dear all,

I am trying to implement a 2PL IRT model with random effects for the discrimination (\alpha) and difficulty (\beta) parameters. I having trouble fitting the model to simulated data, and I am not sure I have coded the model correctly.

Let me start with my motivation behind the model. I have data on about 20 items for a hundred countries for the last 40 years. So I have country-years are nested on countries. I think the questions are not equally relevant across countries. I want to allow the discrimination and difficulty parameters to vary across countries. In IRT lingo, I want to generalize the usual 2PL setup to this: logit[Pr(y=1)]=\alpha_{item, country}\theta_{country,year}-\beta_{item,country}

I have been looking at the case studies and previous posts to get a sense of how to implement this model. I came up with the code below. I try to use a non-centered parameterization for \alpha to help with divergences, and some constraints on \beta to help with identification. In an effort to start simple, I’m currently just trying to get random effects for \alpha while keeping \beta constant across countries.

data { 
  int<lower=1> N;            // number of observations
  int<lower=1> I;            // number of items 
  int<lower=1> C;            // number of countries 
  int<lower=1> K;            // number of unique random effects
  int<lower=1> J;            // number of country-years 
  int<lower=1,upper=I> ii[N];  // item for response n 
  int<lower=0,upper=C> cc[N];  // country for response n  
  int<lower=0,upper=J> jj[N];  // country-year for response n 
  int<lower=1,upper=K> kk[N];  // random effect for response n
  int<lower=0,upper=1> y[N];   // response n

parameters {
  vector[I-1] beta_free; //freely estimated difficulty parameters
  vector[J] theta;    // Ability parameter for country-year              		
  vector[I] a_mean;   // Mean discriminations for each item 
  vector[K] a_offset; // Discrimination offsets for each country
  vector<lower=0> [I] a_sigma; //St dev. for discrimination offsets 
transformed parameters{
  vector [I] beta; // contains freely estimated and constrained betas
  vector[N] alpha; //Overal discriminations, including random effects 
  vector[N] eta;  // log odds of a correct probability
  beta[1:(I-1)] = beta_free;
  beta[I] = -1*sum(beta_free);  // Constrain  betas to help with identification 
  alpha = a_mean[ii]+ a_offset[kk] .* a_sigma[ii];  // non-centered
  eta = alpha .*theta[jj]-beta[ii];
model { 
  theta ~ normal(0,1);
  target += normal_lpdf(beta | 0, 3);
  a_mean ~ normal(0,1);
  a_offset ~ normal(0,1);
  a_sigma ~  inv_gamma(1,1);
  y ~ bernoulli_logit(eta); //likelihood  

I get lots of divergent iterations when I use a small simulated dataset to test the mode. Increasing adapt_delta to .99 helps a bit but I still get a lot of them. Bigger datasets don’t get divergences but the effective number of samples is really small and most \bar{R} remain above 1.1 despite using 5 chains with 3,000 iterations.

I had a lot of trouble with setting the random effects (line 6 under transformed parameters). I am not sure that I got it right. I tried using a non-centered parametrization but maybe I need to use Cholesky decomposition? Still trying to learn how to do this properly.

I am attaching my R script where I simulate the data. I think I did the simulation right but its my first time trying to simulate this kind of data.

I would appreciate any suggestions!
r script data simulation random2pl.R (2.7 KB)