Rasch model

The above data may be analysed using the one-parameter Rasch model (see andersen:80, pp.253-254; boch:aitkin:81). The probability that student j responds correctly to item k is assumed to follow a logistic function parameterized by an `item difficulty' or threshold parameter and a latent variable representing the student's underlying ability. The ability parameters are assumed to have a Normal distribution in the population of students. That is:

The above model is equivalent to the following random effects logistic regression:

where corresponds to the scale parameter ( ) of the latent ability distribution. We assume a half-normal distribution with small precision for ; this represents vague prior information but constrains to be positive. Standard vague normal priors are assumed for the 's. Note that the location of the 's depend upon the mean of the prior distribution for which we have arbitrarily fixed to be zero. Alternatively, Boch and Aitkin ensure identifiability by imposing a sum-to-zero constraint on the 's. Hence we calculate to enable comparision of the BUGS posterior parameter estimates with the Boch and Aitkin marginal maximum likelihood estimates.

Boch and Aitkin compute the following likelihood ratio chi-square statistic (deviance) for testing the assumed model against a general multinomial alternative

on R - 1 - q degrees of freedom, where q is the number of paramters to be estimated, is the observed number of students scoring response pattern i = 1,...,R=32 and is the (unconditional) probability of a randomly selected student responding with pattern i. Conditional on ability level , this probabilty is

where and = 0 or 1 according to the value of the kth item in the ith pattern.

To obtain the unconditional probability , we integrate over as follows

where is the posterior ability distribution. Using Monte Carlo integration and a sample of simulated values for , we may approximate this integral by

where n indexes iteration. That is, we may estimate by the posterior mean of . To generate within BUGS, we just sample a random ability paramter from a standard normal at each iteration, calculate and substitute this into the equation for above. At the end of the BUGS run, we calculate the posterior mean of the 's and use this as our estimate of in the formula for .

The BUGS code for this model is given on the next page, and the graph is shown in Figure 14.

 
Figure 14:   Graphical model for LSAT example.

model LSAT;
const
   N = 1000,    # number students
   R = 32,      # number of possible test results
   T = 5;       # number of tests
var
   response[R,T], m[R], culm[R], alpha[T], a[T], theta[N], r[N,T], 
   p[N,T], beta, theta.new, p.theta[T], p.item[R,T], P.theta[R];
data response,m,culm in "lsat.dat";
inits in "lsat.in";
{
# Calculate individual (binary) responses to each test from multinomial data
   for (j in 1:culm[1]) {
     for (k in 1:T) {  r[j,k] <- response[1,k]; }
   }
   for (i in 2:R) {
     for (j in culm[i - 1] + 1:culm[i]) {
        for (k in 1:T) {  r[j,k] <- response[i,k]; }
     }
   }
# Rasch model
   for (j in 1:N) {
      for (k in 1:T) {
         logit(p[j,k]) <- beta*theta[j] - alpha[k];
         r[j,k] ~ dbern(p[j,k]);
      }
      theta[j] ~ dnorm(0,1);
   }
# Priors
   for (k in 1:T) {
      alpha[k] ~ dnorm(0,0.0001);   a[k] <- alpha[k] - mean(alpha[]);
   }
   beta ~ dnorm(0,0.0001) I(0,);

# Compute probability of response pattern i, for later use in computing G^2
  theta.new ~ dnorm(0,1);          # ability parameter for random student 
  for(k in 1:T) {
       logit(p.theta[k]) <- beta*theta.new - alpha[k];  
       for(i in 1:R) {
          p.item[i,k] <- pow(p.theta[k],response[i,k]) 
                          * pow((1-p.theta[k]),(1-response[i,k]));
       }
  } 
  for(i in 1:R) {    
     # P_i|theta = PROD_k p_k|theta   
     P.theta[i] <- p.item[i,1]*p.item[i,2]*p.item[i,3]*p.item[i,4]*p.item[i,5];
  }
}

Note that the data are read into BUGS in the original multinomial format to economize on space and effort. The 5 1000 individual binary responses for each item and student are then created within the BUGS code using the index variable culm (read in from the data file), where culm[i] = cumulative number of students recording response patterns 1, 2, ..., i; .

Analysis

A 2000 iteration BUGS run took 55 minutes after a 1000 iteration burn-in and produced the following results:

B&A =Boch & Aitkin (1981) Marginal maximum likelihood estimate