Simulating 3 parameter IRT data

# Some Easy Item Response Theory Data Generating Commands

# The following code will draw random test items from several common item response theory forms.

# The command mat.binom will be useful in drawing a matrix of binomial results. I will use a command I programmed in a previous post. (http://www.econometricsbysimulation.com/2012/09/item-response-theory-estimation.html)
mat.binom <- function(n,p) {
  bin.mat <- p*0
  for (i in 1:nrow(p)) {
    for (ii in 1:ncol(p)) {
    # This will draw a random binomial for each of the probabilities.
    bin.mat[i,ii] <- rbinom(1,n,p[i,ii])
  }
  }
 return(bin.mat)
}


# Define a function for the 3 parameter item response model. D is automatically set to 1.7 to create equivalent responses to that of using the normal distribution.
rirt3 <- function(theta=0,a=1,b=0,c=0, D=1.7) {
  nitems = max(length(a),length(b),length(c))
  # Note theta, a, b, or c can be vectors.  The length of the a,b,c vector(s) tells R how many items to draw.
   
  # This creates a matrix of theta values
  m.theta <- matrix(theta, nrow=length(theta), ncol=nitems)
 
  # This creates a matrices of item parameter values for parameters a,b,c called m.a,m.b, and m.c
  for (v in c("a","b","c")) assign(paste("m.",v,sep=""),
         t(matrix(get(v), ncol=length(theta), nrow=nitems)))

  # This code does the same thing as the following code.
  # m.a <- t(matrix(a, ncol=length(theta), nrow=nitems))
  # m.b <- t(matrix(b, ncol=length(theta), nrow=nitems))
  # m.c <- t(matrix(c, ncol=length(theta), nrow=nitems))

  # Now we calculate probabilities of getting the problem right:
  m.p = m.c+(1-m.c)/(1+exp(-D*m.a*(theta-m.b)))
 

 
  # We just need to pass to the probability matrix into the mat.binom function.
  Y = mat.binom(1, m.p)
 
  # Now let's return the items that we have calculated.
  return(Y)
}

# Let's see this command in action
theta=seq(-8,8,1)
a=1
b=seq(-10,10,5)
c=seq(.1,.5,.1)
# b will be matched with c such that -10 is with .1, -5 with .2, 0 with .3, etc.
# We can see how these will be matched by observing
cbind(a,b,c)

rirt3(theta,a,b,c)

# We can see that though the last item is very difficult there are a few who have guessed it.  That is because the guessing parameter is very high.

# Of course the one parameter and the two IRT models are easily gotten from the three parameter if a=1 and c=0.  The one parameter IRT model is also known as the Rasch model.

# Let's now define the Generalized Graded Response model.

rgrm <- function(theta=0,a=cbind(rep(1,5),rep(1,5)),b=cbind(rep(0,5),rep(1,5)), D=1.7) {
  # b is now an input matrix with each row representing a different item and each column representing a different grade for that item.  The number of columns should be equal to the maximum number of grades for all of the items.  Items may have less than the full number of grades indicated by NA values at upper levels.

  # We will use the grm Cululative Grade Function defined in a previous post as what we will use to generate our probabilities.
  # http://www.econometricsbysimulation.com/2012/09/generalize-graded-response-model.html
  grm <- function(theta=1, a=1, b = c(0,1,2), cplot=T, pplot=T, stackpoly=F ) {
    ngrade = max(length(b),length(a))
    CGF = matrix(NA, ncol=ngrade  , nrow=length(theta))
    for (i in 1:length(theta)) CGF[i,] = exp(a*(theta[i]-b))/(1+exp(a*(theta[i]-b)))
    return(CGF)
  }

  # This will be the number of items to generate.
  if (sum(dim(a)!=dim(b))>0) warning()
  nitems = nrow(b)
 
  # This matrix will hold the results of the items.
  Y <- matrix(NA, nrow=length(theta), ncol=nitems)

  print(nitems)
  for (i in 1:nitems) {
    a.sub = !is.na(a[i,])
    b.sub = !is.na(b[i,])
    print(grm(theta=theta, a=a[i,a.sub], b=b[i,b.sub]))
  }
}

rgrm()
rgrm(theta=c(1,2,9,2,3,4))
rgrm(theta=c(1,2,9,2,3,4), a=cbind(rep(1,6),rep(1,6),rep(c(NA,1),3)),b=cbind(rep(1,6),rep(2,6),rep(c(NA,3),3)))