Modified Bin and Union Method for Item Pool Design

# Reckase (2003) proposes a method for designing an item pool for a computer
# adaptive test that has been known as the bin and union method.  This method
# involves drawing a subject from a distribution of abilities. Then selecting
# the item that maximizes that subject's information from the possible set of 
# all items given a standard CAT proceedure. This is repeated until the test
# reaches the predifined stopping point.
 
# Then then next subject is drawn and a new set of items is drawn.  Items are
# divided into bins such that there is a kind of rounding.  Items which are
# sufficiently close to other items it terms of parameter fit are considered
# the same item and the two sets are unionized together into a larger pool.
 
# As more subjects are added more items are collected though at a decreasing
# rate as fewer new items become neccessary.
 
# In the original paper he uses a fixed length test though in a forthcoming 
# paper he and his student Wei He is also using a variable length test.
 
# I have modified his proceedure slightly in this simulation.  Rather than
# selecting optimal items for each subject based from the continuous pool
# of possible items I have the test look within the already constructed pool
# to see if any items are within bin length of the subject's estimated ability.
 
# If there is no item then I add an item that perfectly matches the subject's
# estimated ability. The reason I prefer this method is that I think it better
# represents the process that a CAT test typically must go through with items
# close to but rarely exactly at the level of the subjects. Thus the information
# for each subject will be slightly less as a result of this modified method
# relative to the original.
 
# As with the new paper this simulation uses a variable length test. My stopping
# rule is simple.  Once the test achieves a sufficiently high level of 
# information, then it stops.
 
# I have constructed this simulation as one with three nested loops.
# Over subjects within the item pool construction.
 
# It simulates the item pool construction a number of times to get the
# average number of items after each subject as well as a histogram
# of average number of items required at each difficulty level.
 
# I have also included a control for item exposure.  This control 
# dicatates that as the acceptable exposure rate is reduced, more items
# will be required since some are too frequently exposed.
 
# Overall this method is seems pretty great to me.  It allows for
# item selection criteria, stopping rules, and exposure controls
# to be easily modified to accomidate most any CAT design.
 
require("catR")
 
# Variable Length Test
 
# The number of times to repeat the simulation
nsim <- 10
 
# The number of subjects to simulate
npop <- 1000
 
# The maximum number of items
max.items <- 5000
 
# Maximum exposure rate of individual item
max.exposure <- .2
 
# Stop the test when information reaches this level
min.information <- 10
 
# How far away will the program reach for a new item (b-b_ideal)
bin.width <- .25
 
expect.a <- 1
 
p <- function(theta, b) exp(theta-b)/(1+exp(theta-b))
info <- function(theta, b, a=expect.a) p(theta,b)*(1-p(theta,b))*a^2
 
info(0,0)
 
# The choose.item funciton takes an input thetahat and searches
# available items to see if any already exist that can be used
# otherwise it finds a new item.
choose.item <- function(thetahat, item.b, items.unavailable, bin.width) {
  # Construct a vector of indexes of available items
  avail.n <- (1:length(item.b))
 
  # Remove any already make unusuable
  if (length(items.unavailable)>0) 
    avail.n <- (1:length(item.b))[-items.unavailable]
 
  # If there are no items available then generate the next item
  # equal to thetaest.
  if (length(avail.n)==0) 
    return(c(next.b=thetahat, next.n=length(item.b)+1))
 
  # Figure out how far each item is from thetahat
  avail.dist <- abs(item.b[avail.n]-thetahat)
 
  # Reorder the n's and dist in terms of proximity
  avail.n <- avail.n[order(avail.dist)]
  avail.dist <- sort(avail.dist)
 
  # If the closest item is within the bin width return it
  if (avail.dist[1]<bin.width) 
    return(c(next.b=item.b[avail.n[1]], next.n=avail.n[1]))
  # Otherwise generate a new item
  if (avail.dist[1]>=bin.width) 
    return(c(next.b=thetahat, next.n=length(item.b)+1))
}
 
# Define the simulation level vectors which will become matrices
Tnitems <- Ttest.length <- Titems.taken.N <- Titem.b<- NULL
 
# Loop through the number of simulations
for (j in 1:nsim) {
 
 
  # Seems to be working well
  choose.item(3, c(0,4,2,2,3.3), NULL, .5)
 
  # This is the initial item pool
  item.b <- 0
 
  # This is the initial number of items taken
  items.taken.N <- rep(0,max.items)
 
  # A vector to record the individual test lengths 
  test.length <- NULL
 
  # Number of total items after each individual
  nitems <- NULL
 
  # Draw theta from a population distribution
  theta.pop <- rnorm(npop)
 
  # Start the individual test
  for (i in 1:npop) {
    # The this person has a theta of:
    theta0 <- theta.pop[i]
 
    # Our initial guess at theta = 0
    thetahat <- 0
 
    print(paste("Subject:", i,"- Item Pool:", length(item.b)))
    response <- items.taken <- NULL
 
    # Remove any items that would have been overexposed
    items.unavailable <- (1:length(item.b))[!(items.taken.N < max.exposure*npop)]
 
    # The initial imformation on each subject is zero
    infosum <- 0
 
    # Loop through each subject
    while(infosum < min.information) {
 
      chooser <- choose.item(thetahat, item.b, items.unavailable, bin.width)
 
      nextitem <- chooser[2]
      nextb <- chooser[1]
        names(nextitem) <- names(nextb) <- NULL
 
      items.unavailable <- c(items.unavailable,nextitem)
      item.b[nextitem] <- nextb
 
      response <- c(response, runif(1)<p(theta0, nextb))
 
      items.taken <- c(items.taken, nextitem)
 
      it <- cbind(1, item.b[items.taken], 0,1)
 
      thetahat <- thetaEst(it, response)
 
      infosum <- infosum+info(theta0, nextb)
    }
 
    # Save individual values
    nitems <- c(nitems, length(item.b))
 
    test.length <- c(test.length, length(response))
 
    items.taken.N[items.taken] <- items.taken.N[items.taken]+1
  }
 
  # Save into matrices the results of each simulation
  Titem.b <- c(Titem.b, sort(item.b))
  Tnitems <- cbind(Tnitems, nitems)
  Ttest.length <- cbind(Ttest.length, test.length)
  Titems.taken.N <- cbind(Titems.taken.N, items.taken.N)
 
}
 
plot(apply(Tnitems, 1, max), type="n",
     xlab = "N subjects", ylab = "N items",
     main = paste(nsim, "Different Simulations"))
for (i in 1:nsim) lines(Tnitems[,i], col=grey(.3+.6*i/nsim))
 

# We can see that the number of items is a function of the number of

# subjects taking the exam.  This relationship becomes relaxed

# when the number of subjects becomes large and the exposure controls

# are removed.

hist(Titem.b, breaks=30)

 
hist(Ttest.length, breaks=20)

Created by Pretty R at inside-R.org

Estimating Person Characteristics from IRT Data - 3PL Model

Original Code

# One of the basic tasks in item response theory is estimating the ability of a test taker from the responses to a series of items (questions).

# Let's draw the same pool of items that we have used on several previous posts:

# First let's imagine we have a pool of 100 3PL items.
set.seed(101)

npool = 500

pool = as.data.frame(cbind(item=1:npool, a=abs(rnorm(npool)*.25+1.1), b=rnorm(npool), c=abs(rnorm(npool)/7.5+.1)))
summary(pool)

# Drawing on code from a previous post we can calculate several useful functions:

# (http://www.econometricsbysimulation.com/2012/11/matching-item-information.html)

# Each item has a item characteristic curve (ICC) of:
PL3 = function(theta,a, b, c) c+(1-c)*exp(a*(theta-b))/(1+exp(a*(theta-b)))

# Let's imagine that we have a single test taker and that test taker has taken the first 15 items from the pool.

items.count = 15
items.taken = pool[1:items.count,]

# And that the person has a latent theta ability of 1
theta = 1.3

# Let's calculate the cut points for each of the items.
items.cut = PL3(theta, items.taken$a, items.taken$b, items.taken$c)

# We can see how the cut point works by graphing
plot(0,0, type="n", xlim=c(-3,3),ylim=c(0,1), xlab=~theta,
           ylab="Probability of Correct Response", yaxs = "i", xaxs = "i" , main="Item Characteristics Curves and Ability Level")

for(i in 1:items.count) {
  lines(seq(-3,3,.1),PL3(seq(-3,3,.1), items.taken$a[i], items.taken$b[i], items.taken$c[i]), lwd=2)
  abline(h=items.cut[i], col="blue")
}
abline(v=theta,col="red", lwd=3)

# Now let's draw a uniform draw that we will use to calculate whether each item as passed.

rdraw = runif(items.count)

# Finally, we will calculate item responses
item.responses = 0
item.responses = items.cut > rdraw

###############################################
# Done with Simulation - Time for Estimation

# We want to use the information we know about the items (the item parameters and the responses) in order to estimate a best guess at the true ability of the test taker.

# First we must check if the person got all of the items either correct or incorrect.

sum(item.responses)
# If this is either a 0 or a number equal to the number of items then we cannot esimate an interior maximum without additional assumptions.

# We will attempt to recover our theta value using the r command optim

# First we need to specify the function to optimize over.

MLE = function(theta) sum(log((item.responses==T)*PL3(theta, items.taken$a, items.taken$b, items.taken$c) +
                              (item.responses==F)*(1-PL3(theta, items.taken$a, items.taken$b, items.taken$c))))
# The optimization function takes as its argument the choice variables to be optimized (theta).
# The way the above optimization works is that you specify the probility of each response piecewise.
# If the response is correct, then you count the CDF of theta up to that point as contributing to the probability of observing a correct outcome.  If the response is negative, then you count it as contributing the the probability of a incorrect outcome.  You then choose the theta that produces the greatest total probability.

MLEval = 0
theta.range = seq(-3,3,.1)
for(i in 1:length(theta.range)) MLEval[i] =MLE(theta.range[i])

plot(theta.range, MLEval, type="l", main="Maximim Likelihood function", xlab= ~theta, ylab="Sum of Log Likelihood")
abline(v=theta, col="blue")
# We can visually see that the maximum of the slope will not be at the true value though it will be close.

optim(0,MLE, method="Brent", lower=-6, upper=6, control=list(fnscale = -1))
abline(v=optim(0,MLE, method="Brent", lower=-6, upper=6, control=list(fnscale = -1))$par, col="red")

# We can see that we can estimate theta reasonably well with just 15 items from a paper test (red line estimate, blue line true).  However, looking at the graph of the ICCs, we can see that for most of the items, the steepest point (where they have the most discriminating power) is at an ability set lower than the test taker's ability.  Thus, this test provides the most information about a person who has a lower ability than the person with a theta=1.2.

# We can use R's optim function to find the ideal theta that would maximize the information from this test.
# Item information is:
PL3.info = function(theta, a, b, c) a^2 *(PL3(theta,a,b,c)-c)^2/(1-c)^2 * (1-PL3(theta,a,b,c))/PL3(theta,a,b,c)

# Notice, this is not the best way of defining the test information function since the items are not arguments.
test.info = function(theta) sum(PL3.info(theta, items.taken$a, items.taken$b, items.taken$c))

# Construct a vector to hold the test information
info = 0
for(i in 1:length(theta.range)) info[i]=test.info(theta.range[i])

plot(theta.range, info, type="l", main="Information Peaks Slightly Above 0", xlab= ~theta, ylab="Information")
abline(v=theta, col="blue")
# But we want to know about the test taker at theta

optim(0,test.info, method="Brent", lower=-6, upper=6, control=list(fnscale = -1))
# The person this test would be best suited to evaluate would have an ability rating of .19
abline(v=optim(0,test.info, method="Brent", lower=-6, upper=6, control=list(fnscale = -1))$par, col="red")

Generating Tests of Particular Properties

# This post will demonstrate how to generate total test score distributions from a set of 50 items.

# It will build on a previous post (http://www.econometricsbysimulation.com/2012/09/simulating-3-parameter-irt-data.html) that demonstrated how to easily draw a matrix of 3 parameter univariate item responses given a vector of thetas.

# In order to run this code, first run the previous code in order for the command rirt3 to be defined.

# rirt takes a vector of thetas and vectors of the a,b, and c parameters and constructs a response matrix with rows equal to the number of test takers and columns equal to the number of items.

Item.Scores <- rirt3(theta=seq(-2,2,.1) ,a=c(1,.5,1,1,1),b=c(1,2,3,4,5),c=c(.15,.1,.2,.025,.1))
Total.Score <- apply(Item.Scores,1,sum)
hist(Total.Score, breaks=0:max(Total.Score), main="Total test score - 5 item test")

# Thus will create three rows for forty individuals with two items taken per individual.

# Let's look at a practical example.  Imagine we have 10000 test takers with ability drawn from a normal distribution.

theta = rnorm(10000)

# And we get to choose 50 items (or rather their parameters) to build a test from.

# We know on average about 95% of our sample falls within two standard deviations of the mean zero.

# First let's assume that we are restricted to items that have equal discriminatory power a=1 and equal guessing probability = .1.

# 1. Assignment one.  Build a test that has total scores that looks normally distributed.

# Since the underlying parameter distribution is normally distributed we should only need offer items that span the relevant range.

b=seq(-2.5,2.5, length.out=50)

Item.Scores <- rirt3(theta=theta ,a=2,b=b,c=.1)
Total.Score <- apply(Item.Scores,1,sum)
hist(Total.Score, breaks=0:max(Total.Score), main="Total test score - 50 item test")

plot(theta, Total.Score, main="Theta is an excellent predictor of total score")
# From the steady angle of this plot we can see at all levels this test does generally equally well at discriminating between different levels of theta.

# However, we may be interested less in how well average students perform and more in how well we can discriminate between top students who we would like to give scholarships to as well as bottom students which we would like to refer to remedial studies as necessary.

# 2. Assignment two.  Build a test that has more power to discriminate at the lower end and at the upper end of the ability distribution.

# In order to do this let us divide the students into three groups.  Low group and high group each get 20 questions while middle group instead gets 10 questions.  Low group starts as -3 and goest to -1.1, middle -1 to 1, and high 1.1 to 3.

b=c(seq(-3,-1.1,length.out=22), seq(-1,1,length.out=6), seq(1.1,3, length.out=22))

Item.Scores <- rirt3(theta=theta ,a=2,b=b,c=.1)
Total.Score <- apply(Item.Scores,1,sum)
hist(Total.Score, breaks=0:max(Total.Score), main="Total test score - 50 item test")

# It is not obvious that this is doing what we want it to be doing from the histogram.  This is because the middle most bins are now more crowded making all of the other bins look small.

plot(theta, Total.Score, main="Around the middle it is hard to tell between ability level")

# We can see from the plot that we have less movement in total scores around the mean ability level while high and low ability levels tend to be well discriminated.

# 3. Assignment three.  Construct a test in which most total scores greater than 10 and less than 40 are equally likely to occure.

# In order to do this we want to make sure parts of the theta distribution in which there are many students stacked also have a large number of questions causing spread within that group.

# We should be able to do this by using the inverse of the normal CDF.

b=qnorm(seq(.01,.99,length.out=50))

Item.Scores <- rirt3(theta=theta ,a=2,b=b,c=.1)
Total.Score <- apply(Item.Scores,1,sum)
hist(Total.Score, breaks=0:max(Total.Score), main="Total test score - 50 item test")

plot(theta, Total.Score, main="Detecting differences between students is equally distributed")

# It is not obvious from this plot because there are more thetas grouped in the middle where the plot is steepest.  If we rank the thetas on the other hand it becomes very clear.

plot(rank(theta), Total.Score, main="Detecting differences between students is equally distributed")

# Try ranking the thetas on the other plots.  You will find that the other plots will fatten out even more near the center of the plot.  That is because without stacking extra items near the center of the distribution of abilities it is hard to tell the difference between densely stacked students.