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")

Selecting your First Item on a Computer Adaptive Test

Original Code

# Computer adaptive tests "adapt" to test taker ability by making assessments of the test taker's ability and providing questions that are meant to maximize the amount of information that can be inferred from the test.

# We often start by assuming that the ability score of an individual is at the average for the population to begin with.

# Once that assumption is made then the adapative test selects an item that maximizes the information at that point.


# Let's imagine that the true ability of the student is -1 and we would like to select items that get us from 0 to -1.

# Let's see how this works.

#######################################################
#  Method 1 - Item Information, Fisher Information Criteria

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

npool = 100

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)))

# and information function defined as:
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)

# We can use the previous post to find the information for any number theta values.
# For now we are only interested in the information for our itial "guess" at student ability:
theta.estimate=0

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

# and information function defined as:
PL3.info = function(theta.estimate, a, b, c) a^2 *(PL3(theta.estimate,a,b,c)-c)^2/(1-c)^2 * (1-PL3(theta.estimate,a,b,c))/PL3(theta.estimate,a,b,c)

# First I want to calculate the information at each theta for each item.

# In this first case theta only equals 0.
nthetas = length(theta.estimate)

# This matrix will contain the item pool information values
pool.info = matrix(NA, nrow=npool, ncol=nthetas)

colnames(pool.info) = paste("IT=",theta.estimate,sep="")

# These must be calculated for each item but can be for all thetas simultaneously.
for(i in 1:npool) pool.info[i,] = PL3.info(theta.estimate, pool[i,2]*1.7, pool[i,3], pool[i,4])

head(pool.info)

# Everything appears to be working well.  Let's find the max (note this code only works well for one theta).
cbind(pool,pool.info)[pool.info==max(pool.info),]

# Item 33 with am information of 1.2695 with a=1.35, b=-.14, and c=.007 is the best choice for the first item using the fisher information criteria.

#######################################################
#  Method 2 - Kullback–Leibler information divergence method

# Referencing Eggen of Cito (1999) the KL information divergence for a single item can be expressed as
# KL = p(true.theta)*log(p(true.theta)/p(estimated.theta)) + q(true.theta)*log(q(true.theta)/q(estimated.theta))
# where p is the probability of getting that item correct and q is the probability of getting that item wrong.

# KL can be thought of as an item selection criteria that is likely to give you the best item to distinguish between your current estimate and the true (if the true was known).

# In order to use information divergence we find the expected value of each item.
# Computationally this is approximated by inputing possible true.theta and their estimate and finding the average.


# Let's first let's define the KL function (drawing on PL3 defined previously)

KL = function(theta.true,theta.estimate, a, b, c) {
  # For the true value
  p.true = PL3(theta.true,a,b,c)
  q.true = 1-p.true
 
  # For the estimate
  p.estimate = PL3(theta.estimate,a,b,c)
  q.estimate = 1-p.estimate

  # The following line is the value to be returned to the KL function:
  p.true*log(p.true/p.estimate) + q.true*log(q.true/q.estimate)
  }
# The function is written to only take a single theta.estimate while multiple true theta's should not be a problem.

# Now let's specify a simplified discrete range that the true theta can take.

theta.true = c(-1,0,1)

nthetas = length(theta.true)

# This matrix will contain the item pool KL values
pool.KL = matrix(NA, nrow=npool, ncol=nthetas)

colnames(pool.KL) = paste("Theta=",theta.true,sep="")

# These must be calculated for each item but can be for all thetas simultaneously.
for(i in 1:npool) pool.KL[i,] = KL(theta.true, theta.estimate, pool[i,2]*1.7, pool[i,3], pool[i,4])
head(pool.KL)
# Note that for theta.true = 0 the entire column is equal to zero.
# This is because if the true is equal to the estimate than there is no item that maximize the ability to tell the difference between the estimate and the true because there is no difference.

# In order to select the best item for a true theta range of -1,0,1 we average across the three values (or sum them).

pool.avKL = apply(pool.KL, 1, mean)
head(pool.avKL)

cbind(pool,pool.info,pool.avKL)[pool.avKL==max(pool.avKL),]
# Interestingly, item 33 is once again the item which is selected.

# Now let's imagine that rather than picking just three true theta's to search for items over that we instead want to search across the standard normal standarized distribution of abilities theta.  This is very easy to do with the code we already have.  All we need do is take draw from the normal distribution.

theta.true = qnorm(seq(.01,.99,.01))

hist(theta.true)

nthetas = length(theta.true)

# This matrix will contain the item pool KL values
pool.KL = matrix(NA, nrow=npool, ncol=nthetas)

colnames(pool.KL) = paste("Theta=",theta.true,sep="")

# These must be calculated for each item but can be for all thetas simultaneously.
for(i in 1:npool) pool.KL[i,] = KL(theta.true, theta.estimate, pool[i,2]*1.7, pool[i,3], pool[i,4])
pool.avKL = apply(pool.KL, 1, mean)
head(pool.avKL)

cbind(pool,pool.info,pool.avKL)[pool.avKL==max(pool.avKL),]
# Surprisingly, once again item 33 is selected.  I am not sure why there is no difference in methods yet.  I suspect it is due to the symetric nature of the normal distribution.

# We can certainly play with the models at least.  I reran the code starting with theta.estimate=2 and the fisher information criteria did select a different item than the KL though the KL selected the same item with both choices of theta.  This is probably because the item pool is so small.  Given a larger item pool, I suspect that distributional choices become more imporant.

# Changing the item pool to 10,000 and theta.estimate=1 I ended up not finding any difference in item choice between only the three values -1,0,1 and the full normal distribution of values.

# Let's do one more thing.  Let us imagine that we think we know the true theta and we would like to use the KL to select an item that brings us closest to the true.

theta.true = -1
nthetas = length(theta.true)

pool.KL = matrix(NA, nrow=npool, ncol=nthetas)

colnames(pool.KL) = paste("Theta=",theta.true,sep="")

# These must be calculated for each item but can be for all thetas simultaneously.
for(i in 1:npool) pool.KL[i,] = KL(theta.true, theta.estimate, pool[i,2]*1.7, pool[i,3], pool[i,4])
pool.avKL = apply(pool.KL, 1, mean)
head(pool.avKL)

cbind(pool,pool.info,pool.avKL)[pool.avKL==max(pool.avKL),]

# In that case item 72 is selected with parameters a=1.5, b=-.43, c=.097.

# Ths KL difference criteria seems to be working well since it selected an item that is between the true -1 and the estimate 0.

Item Information - IRT

Original Code

# Interestingly, one of the founders of Item Response Theory (Fredrick Lord) developed his own concept of information.  His approach was both unique and interesting but ended up leaving us with the same Fisher Information Matrix.

# The Fisher Information equation is an equation that measures for any particular value of an estimate, what the particular amount of information you have on that estimate is.  It is a difficult subject that I struggle with.

# It might be useful to look at how it can be found.

# First keep in mind that for Fisher Information we are primarily concerned with Maximum Likelihood.

# For maximum likelihood we are primarily concerned in maximizing the ln(pdf(theta)) by choosing theta.

# Define L = ln(pdf(theta))

# The score is equal to the first derivative of the L with respect to thera.  Define Score = S = dL/dTheta = d ln(pdf(theta))/d theta

# We know that if the MLE maximization routine has worked properly then the score is equal to zero.

# Now the Fisher Information is equal to the expected value of the score squared (the second moment) given theta:

# I = E(S^2|theta)

# Now, the tricky part is thinking how to make sense of this thing we call information.

# First let's think how it is used.

# For one thing the standard error of estimation is se=1/(I)^.5
# In a similar vien, it can be shown that the variance of any unbiased estimator (f) of theta has a lower limit of 1/I: Var(f(theta)) >= 1/I

# Thus we can see that there is a direct inverse relationship between information and the variance of our estimator.

# How do we reconcile this?

# Item Response theory has some very nice applications of information that I believe shed light on other uses of information.

# For the Rasch Model the item information for one item is.  I(theta) = p(theta)(1-p(theta)

# Where p is the probability of getting the item correct.

# Let's map this out for a range of thetas assuming the difficulty parameter is one.

theta = seq(-4,6,.1)

# Now let's define the Rasch function:

rasch = function(theta,b) exp(theta-b)/(1+exp(theta-b))

# Let's also define the information value.

rasch.info = function(theta,b) rasch(theta,b)*(1-rasch(theta,b))

plot(theta, rasch.info(theta,1), ylab = "Information", main="The most information is gained at theta=1 (1PL)", type="l")

# We can see that information peaks at theta=1.  What does this mean for practical purposes?  If you want to know the most about a student's ability give them a question in which the difficulty of the item is at their ability level.  If you though give them a question that is far too easy or far too hard then even if they do well on the question or poor on the question you have not learned that much about their ability level on the theta scale.

# Why?  Because the likelihood of someone else doing equally well on that question who is close on the theta scale is equivalent.

# We can see this by looking at the item characteristic curve:

plot(theta, rasch(theta,1), main="Item Characteristic Curve", ylab="Probability of Correct Response", type="l")

# We can see that the change in the probability of getting the item correct is largest at theta=1.  However, as theta gets very large or very small there is very little change in the probabilities (within each prospective range) of getting the item correct as a result of a small change in theta.

# Another way of thinking about this is: if we had two students in the room and both of them were about the same ability and we wanted to know which was stronger, we would want to give them the question which was most closely aligned at the difficulty related to their ability.

# Let's look at a few more different IRT models.

# The 2 parameter logistic model is very similar to the single parameter:

# I(theta) = a^2*p(theta)*(1-p(theta))

# The only difference exists in that the two parameter model allows for their to be more or less information generated from each item as a result of the "discriminatory" power of the item.  Though it is worth noting that a high a parameter model does not strictly dominate a low a parameter model in terms of information for all values of theta.  This is because the p*(1-p) is also a function of theta.  Let's see this in action:

PL2 = function(theta,a, b) exp(a*(theta-b))/(1+exp(a*(theta-b)))

PL2.info = function(theta, a, b) a^2*PL2(theta,a,b)*(1-PL2(theta,a,b))

plot(theta, PL2.info(theta,2,1), ylab = "Information", main="Larger a creates a larger peak but shallower tails", type="l")

lines(theta, PL2.info(theta,.75,1))

# We can see that we greatly prefer a=2 for theta values between -1 and 3.

# However, the item with less discriminatory power may have more information in the tails.  This is somewhat difficult to understand.  On a test, a good example may be imagine two different questions.  One question, is a arithmatic question for which students must demonstrate knowledge of long division.  The alternative, is a question in which students answer a word problem by piecing together components and understanding the concepts behind the math.  The first question may be a better question for identifying if a student either knows or does not know long division.  However, the second question, may be less good at identifying specific mastery at a particular skill level, but rather demonstrates the ability of the student to pull together various math concepts into a coherent answer.  Thus we might not be able to infer much about arithmatic ability if a student answers this question correctly.  But, any student answering this question correctly whatever, their ability level tells us something about their ability.

# We will also look briefly at the 3 parameter Logistic Model:

# I(theta) = a^2 *(p-c)^2/(1-c)^2 * (1-p)/p

# It can be shown that as c increases, the information function monotonically decreases.  This makes sense in that the guessing parameter c is the opposite of information.  As c gets larger, the likelihood of the student getting the question right by pure chance also gets larger.

PL3 = function(theta,a, b, c) c+(1-c)*exp(a*(theta-b))/(1+exp(a*(theta-b)))

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)

plot(theta, PL3.info(theta,1.2,1,.4), ylab = "Information", main="3PL information a=1.2, b=1, c=.4", type="l")

# We can see that the 3 parameter logistic model has an asymetric form.  This is the result of the guessing parameter.  The problem with respondents having the ability to guess correctly is that we end up having less information on low ability respondents than we would like.

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.

Item Response Theory Estimation

# Item response theory (IRT) is very tricky in a sense because it requires both estimation of the regressors (student ability) as well as estimation of the parameters of the items.

# IRT, I have heard that this somehwhat complicated task can be accomplished though a series of maximum likelihood estimations.  First by specifying specific student abilities then by estimating the item parameters.  Then estimating the student abilities by the item parameters, ect.

# Let's start by simulating some data.

nitem = 40

nstud = 100

# To begin with let's start with a single parameter IRT model

# Probability of getting the item i right for person j is is R(x_i=1)=exp(theta_j-b_i)/(1+exp(theta_j-b_i))

# We will range our item difficulties from -4 to 4

vec.b = seq(-4,4,length.out=nitem)

# We will also range our student ability from -4 to 4

vec.theta = seq(-4,4,length.out=nstud)

# Each student will have a probability p of getting each question right.

# Let's first make matrices from theta and the b parameter sets.
b = t(matrix(vec.b,nrow=nitem,ncol=nstud))
theta = matrix(vec.theta,nrow=nstud,ncol=nitem)
  # Inputing vectors into the matrix command will cause them to be duplicated to fill out the values of the matrices.  The only trick is making sure that they are read in the right direction.  The command writes the input vector row by row.  Thus the b matrix becuase the rows are constant must be transposed while the theta matrix does not.

# Each column is an item and each row is a different student.

# Now let's calculate the probability of any person answering any of the questions correctly.
p = exp(theta-b)/(1+exp(theta-b))

# The next graph requires the plotrix package.
require(plotrix)
stackpoly(t(p),col=gray(seq(0,1,length=nstud)),
    border=gray(seq(.3,1,length=nstud)-.2),
    xlab="b - Difficulty Parameter",
    ylab= "Probability of Correct Response",
    main="Mapping of Probability of Correct Responses for Homogenous Items")

# Each line represents the response probability by a singe student.  All students have lower probability of correct response as the items get more difficult.

# Create a function that will draw a matrix of binomial responses
mat.binom <- function(n,p) {
  bin.mat <- p*0
  for (i in 1:nrow(p)) {
    for (ii in 1:ncol(p)) {
      # This
      bin.mat[i,ii] <- rbinom(1,n,p[i,ii])
    }
  }
  return(bin.mat)
}

# Now let's generate the actual responses to the probabilities
y = mat.binom(1,p)

total.score <- apply(y,1,sum)

plot(vec.theta, total.score, xlab=~ theta, ylab="Total Score",
                    main="Total Score as a Function of Ability")

bhat <- rep(NA,nitem)
for (i in 1:nitem) {
  bhat[i] <- (glm(y[,i] ~ 1 , family=binomial("logit")))[[1]]
}

plot(vec.b,-bhat, xlab= "b - Item Difficulty", ylab= "Estimated Difficulty")
# Actually our graph is looking pretty good.  The scale is off but that is really unimportant since latent trait scales are arbitrary anyways.

# To scale the b I first make it so that the min is zero and the range is 1.  Then I multiply by the desired range (8) and add 4 to get the mean to 0.
b.scaled <- (bhat-min(bhat))/(max(bhat)-min(bhat))*-8+4

plot(vec.b,b.scaled, xlab= "b - Item Difficulty", ylab= "Estimated Difficulty", main="Estimates of Item Difficulty")

fit.line = lm(b.scaled~vec.b)

abline(fit.line, col="red")

# Interestingly, it is actually unneeded to do anything more to estimate student ability since the one parameter Rasch model is fully identified as a function of just the number of items each student got correct.

# Let's see it in action.

theta.hat <- rep(NA,nstud)
for (i in 1:nstud) {
  theta.hat[i] <- print(glm(y[i,] ~ 1 , family=binomial("logit"))[[1]])
  print(paste(i))
}


plot(vec.theta, theta.hat , xlab= ~ theta, ylab= "Estimated Student Ability", main="Estimated Student Ability Unscaled")
# Actually our graph is looking pretty good.  The scale is off but that is really unimportant since latent trait scales are arbitrary anyways.

theta.scaled <- (theta.hat-min(theta.hat))/(max(theta.hat)-min(theta.hat))*8+4

plot(vec.theta, theta.scaled, xlab= ~ theta, ylab= "Estimated Student Ability", main="Estimated Student Ability Scaled")

fit.line = lm(theta.scaled ~vec.theta)

abline(fit.line, col="red")

# This shows us that even a simple estimation technique, (two series of logit regressions) can we quite effective when we are estimating only one parameter for both students and one for the items.

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.

Playing around with IRT Graphs

# It is often useful to plot item response theory (IRT) graphs

# I will use the three parameter IRT model:

# ie. P(X=1) = c + (1-c)*(exp(a*(theta-b))/(1+(exp(a*(theta-b)))

# Let's first declare it as a function

irt = function(theta,a,b,c) c + (1-c)*(exp(a*(theta-b))/(1+(exp(a*(theta-b)))))

# Thus we can calculate the probability that a person with a theta = 1 will get a problem right that has parameters a=.3,b=2,c=.2.

irt(theta=1,a=.3,b=2,c=.2)

# R is convient in that functions can take either vectors or scalars.  Thus if irt is given a scalar and a vector it will reuse the scalar multiple times to be of length equal to the vector.

# So to map out an ability range that we are interested in, let's first define a vector for the relevant range of theta.
theta_map = seq(-4,4,.1)

# Let's see how the IRT function will return a probability response vector of length equal to theta:
irt(theta=theta_map,a=.3,b=2,c=.2)

# Let's see it plotted:
plot(theta_map,irt(theta=theta_map,a=.3,b=2,c=.2), type="l", ylim=c(0,1))

# This would probably not be a very useful item becuase there is not much variation in response probability based on ability level of respondents.

##########################    VARYING a      ############################
# Let's map a few competing items. First let's start with an empty plot:
plot(theta_map, 0*theta_map, type="n", ylim=c(0,1),
        xaxs = "i", yaxs = "i",
        ylab = "Probability of Correct Response",
        xlab = ~ theta,
        main="The Effect of Varying a")

# In many IRT models the a parameter is the "discrimination parameter"
# Let's see what happens when we vary it.
for (i in seq(0,2,.1)) {
  lines(theta_map,irt(theta=theta_map,a=i,b=2,c=.2), type="l", ylim=c(0,1))
}

# We can see that varying a from 0 to 2 leads the item to have more ability to "discriminate" between respondents with ability less than b and those with ability more than b.

# The lowest a can go and still make sense is 0 while there is no upper limit.  As a approaches infinity, the function becomes a step function at the point b.
##########################    VARYING b      ############################
plot(theta_map, 0*theta_map, type="n", ylim=c(0,1),
        xaxs = "i", yaxs = "i",
        ylab = "Probability of Correct Response",
        xlab = ~ theta,
        main="The Effect of Varying b")

# In many IRT models the b parameter is the "difficulty parameter"
# Any student with ability above that parameter should find the chance of answering it greater than  c+(1-c)/2 = (1+c)/2.  That is the probability of guessing at the right answer plus half the remaining probability needed to be spanned.
for (i in seq(-3,3,.3)) {
  lines(theta_map,irt(theta=theta_map,a=1,b=i,c=.2), type="l", ylim=c(0,1))
}

# By varying b we map out the ranges of abilities that students posses.  Thus varying difficulty of problems allows us to guess at were on the theta spectrum any student may lie.

##########################    VARYING c      ############################
plot(theta_map, 0*theta_map, type="n", ylim=c(0,1),
        xaxs = "i", yaxs = "i",
        ylab = "Probability of Correct Response",
        xlab = ~ theta,
        main="The Effect of Varying c")

# In many IRT models the b parameter is the "difficulty parameter"
# Any student with ability above that parameter should find the chance of answering it greater than  c+(1-c)/2 = (1+c)/2.  That is the probability of guessing at the right answer plus half the remaining probability needed to be spanned.
for (i in seq(0,1,.05)) {
  lines(theta_map,irt(theta=theta_map,a=1,b=0,c=i), type="l", ylim=c(0,1))
}

# As c gets large the probability span of the IRT function goes from 1 to 0.  This means that the difference in probabily as a result of ability is decreasing as the guessing parameter gets large.  Thus, a test item that has a high guessing parameter in general has less ability to discriminate between low ability respondents and high ability respondents.  Thus, true or false questions may be uniformative indicators of student ability in many cases.


##########################    Overview      ############################

# What do these graphs tell us about the ideal parameter set of our items?

# In general higher a is better so as to help us discriminate between ability levels.

# There is no obvious ideal level of b.  If we have a specific ability level that we would like to make sure all of our students reach then having all of our items share a similar b might be ideal.  However, if on the other hand we would like to be able to asses the overall ability of students then having items with bs that vary along our entire range of interest might be preferable.

# As for the guessing parameter, I am pretty sure we always want that to be smaller.  That is, if it is harder for students to guess the right answer then when students get the right answer, we have more confidence that that reflects their true ability rather than a good guess.