Generalized Graded Response Model

# Graded Response Model

# The graded response model (grm) by Fumiko Samejima is an extension of the dichotomous item response model developed by Fredrick Lord.

# The graded response model allows for answers to items (test questions) to have values that are between the full value for the answer and no value.  An example of this would be an algebra question that gave points for steps successfully completed even if the final answer was not correct.

# This post will define a function that generates the probability each score value or lower as well as the probability of getting given a specific ability level (this is very closely related to CDFs and PDFs).

grm <- function(b = c(0,1,2), a=1, theta=1, cplot=T, pplot=T, stackpoly=F ) {
  # a is the discrimination parameter for the levels of the item.  If a is a constant then the item is "homogenous".  Otherwise it is heterogenous and each grade level has its own specified a.

  # Let's first check if the b are arranged from lowest to highest.
  if (sum(rank(b)!=1:length(b))>0) warning ("b must be ranked from lowest to highest")
    # The rank function will rank b from 1 to the number of observations in b
    # If that ranking does not align with the vector from 1 to length(b) then we have a problem.
 
  # c=T will cause the cumulative grade graph to be plotted.
  # p=T will cause the probability graph graph to be plotted.

  # Count the number of grades.
  ngrade = max(length(b),length(a))

  # Expand a to or b to be the same length
  if ((length(a)==1)&(length(b)>1)) a<-rep(a,length(b))
  if ((length(b)==1)&(length(a)>1)) b<-rep(b,length(a))

  # First let us define a matrix that will return results
  CGF = matrix(NA, ncol=ngrade  , nrow=length(theta))
  PG = matrix(NA, ncol=ngrade+1, nrow=length(theta))
    # Each column will be for a different grade while each row will be for a different theta (if input)

 for (i in 1:length(theta)) {

  # The inverse cumulative grade function (this function returns the probability that grade of a particular value or higher will be returned.
  CGF[i,] = exp(a*(theta[i]-b))/(1+exp(a*(theta[i]-b)))

  # The probability of getting a 0 on an item is the probability of not getting any points.
  PG[i,1] <- 1 - CGF[i,1]
  # The probability of getting the highest grade is the same as the probability of getting the highest grade or more.
  PG[i,ngrade+1] <- CGF[i,ngrade]

   # For the grades in between max and min the values are more dynamically generated.
    for (ii in 1:(ngrade-1)) {
     PG[i,ii+1] <- CGF[i,ii] - CGF[i,ii+1]
    }
  }

  # Plot the graphs of interest.

  # First let's set the number of plot frames to 1 as default
  par(mfrow=c(1,1))
  if ((cplot==T)&(pplot==T)) par(mfrow=c(2,1))

  plottitle <- paste("G1", " a=",a[1]," b=",b[1], sep="")
  for (i in 2:ngrade) {
            plottitle = paste(plottitle, "; G", i, sep="")
            if (sd(a)!=0) plottitle = paste(plottitle, " a=",a[i], sep="")
            if (sd(b)!=0) plottitle = paste(plottitle, " b=",b[i], sep="")
  }


  # If cumulative grade plot is to be graphed but stackpoly is not (the shaded graph) then do the following
  if ((cplot==T) & (stackpoly==F)) {
    # Plot an empty graph
    plot(c(min(theta),max(theta)),c(0,1), type="n",
        ylim=c(0,1), ylab = "Probability", xlab= ~ theta,
        main=plottitle)
    # Plot the individual graded probabilities
    for (ii in 1:ngrade) {
      lines (theta,CGF[,ii])
    }
  }

  # If stack poly is set to on we will use the stack poly setting to make a graph with shades
  if ((cplot==T) & (stackpoly==T)) {
   require(plotrix)
   stackpoly(matrix(theta,nrow=length(theta),ncol=ngrade),CGF[,ncol(CGF):1],
            col=gray(seq(0.1,0.9,length=ngrade)), border="black",
            ylab="Probability", xlab= ~ theta,
            main=plottitle)
  }

  # If the plot densities is set to on then this will plot the densities in their own graph.
  if (pplot==T) {
    plot(c(min(theta),max(theta)),c(0,1), type="n",
        ylim=c(0,1), ylab = "Probability", xlab= ~ theta,
        main=plottitle )
    for (ii in 1:(ngrade+1)) {
      lines (theta,PG[,ii])
    }
  }

  # This will return both the probability of each grade value for each theta value
  # as well as the probability of that grade or higher in the C graph.
  return(list(PG=PG,CGF=CGF))
}


test1 = grm(b=c(-2,-1,0,1,2), a=2, theta=seq(-4,4,.1))
# A homogenous item with a somewhat low discrimination parameter leading to relatively small probabilities of any partial outcome if the item taking population has theta parameters uniformly distributed between -4 and 4

test2 = grm(b=c(-2,-1,0,1,2), a=2, theta=seq(-4,4,.1), stackpoly=T )
# Setting stackpoly on will generate better looking graphs.
# However this requires that the package plotrix is installed

test3 = grm(b=c(-2,-1,0,1,2), a=4, theta=seq(-4,4,.1), stackpoly=T )
# Increasing the discrimination parameter will increase the steepness of the IRT curves and thus the peak size of the graded responses.

test4 = grm(b=c(-2,-1,1,2), a=4, theta=seq(-4,4,.1), stackpoly=T )
# Removing the middle scoring will cause the probability of the new middle outcome to double.

test5 = grm(b=c(-3,-1,0,1,3), a=c(4,2,3,4,5), theta=seq(-4,4,.1), stackpoly=T )
# Creating a heterogenious item will cause the curves to no longer be symetric and the peaks to vary in size.