An R Function to Calculate Response Probabilities in IRT

I’d like to write an R function on how to calculate item response probabilities for Item Response Theory (IRT) models. In this article I will consider only dichotomous items.

In the previous post I mentioned about three parameter model. In this post I will show a more general model: Four parameter model introduced by Barton and Lord (1981) to alleviate excessive punishment of high ability examinees when they make an error.

\displaystyle {{P}_{i}}\left( \theta \right)={{c}_{i}}+\left( d-{{c}_{i}} \right)\frac{{{e}^{D\cdot {{a}_{i}}\left( \theta -{{b}_{i}} \right)}}}{1+{{e}^{D\cdot {{a}_{i}}\left( \theta -{{b}_{i}} \right)}}}

In this model, θ is the ability of an examinee, P(θ) is the probability of correctly responding i’th item for this examinee, d is the upper asymptote, c is pseudo-guessing parameter, a is item discrimination parameter, b is item difficulty parameter for I‘th examinee. D is scaling constant, it’s value is 1.7.

If d parameter is set 1, we get three parameter model:

\displaystyle {{P}_{i}}\left( \theta \right)={{c}_{i}}+\left( 1-{{c}_{i}} \right)\frac{{{e}^{D\cdot {{a}_{i}}\left( \theta -{{b}_{i}} \right)}}}{1+{{e}^{D\cdot {{a}_{i}}\left( \theta -{{b}_{i}} \right)}}}

If pseudo-guessing parameter c is set to 0, we get two parameter model:

\displaystyle {{P}_{i}}\left( \theta \right)=\frac{{{e}^{D\cdot {{a}_{i}}\left( \theta -{{b}_{i}} \right)}}}{1+{{e}^{D\cdot {{a}_{i}}\left( \theta -{{b}_{i}} \right)}}}

If item discrimination value fixed for model and does not change from item to another item we get one parameter model:

\displaystyle {{P}_{i}}\left( \theta \right)=\frac{{{e}^{D\cdot a\left( \theta -{{b}_{i}} \right)}}}{1+{{e}^{D\cdot a\left( \theta -{{b}_{i}} \right)}}}

Finally, if we set item discrimination parameter to 1, we get Rasch model:

\displaystyle {{P}_{i}}\left( \theta \right)=\frac{{{e}^{D\left( \theta -{{b}_{i}} \right)}}}{1+{{e}^{D\left( \theta -{{b}_{i}} \right)}}}

Let’s write a basic R function to calculate probability of correct response to item with parameters (a,b,c,d) for an examinee with ability θ.

irt.3pl.basic <- function(t,a=1,b,c=0,d=1,D=1.7) {
   return( c + (d-c) / (1 + exp(-a * D * (t - b))) )

Let’s run this function with two examples. In the first example I’ll calculate the probability of correct response for an examinee with ability 0 to a question with item parameters a=1,b=0,c=0,d=1.

theta <- 0;
a <- 1;
b <- 0;
c <- 0;
d <- 1;

[1] 0.5

Since examinee’s ability and items difficulty corresponds there is a 50% chance that this examinee will correctly answer this question. In the next example let’s consider three examinees with theta’s -1,0 and 1.

theta <- c(-1,0,1)
a <- 1
b <- 0
c <- 0
d <- 1;

[1] 0.1544653 0.5000000 0.8455347

Results are as expected. Examinee with ability -1 has a low probability of correctly answering this question. Reverse is true for examinee with high ability.

In the next step, I’ll graph item characteristic curve of an item with parameters a=0.7; b=0.5; c=0.2; d=.98:

theta <- seq(-4,4,by=.1)
Probability <- irt.3pl.basic(t=theta,a=0.7, b=0.5, c=0.2, d=.98)

In the next plot I will extend this function to a multi item and multi ability case.


Barton, M.A., & Lord, F.M. (1981). An Upper Asymptote for the Three-Parameter Logistic Item-Response Model. Arlington, Va.: Office of Naval Research, Personnel and Training Research Programs Office.

A Basic Demonstration of Three Parameter Logistic Model in R

Item response theory’s three parameter logistic model is well known. I leave the explanation of this model to major sources (such as Lord (1980), page 12). The equation is as follows:

P\left(\theta\right)=c+\frac{1-c}{1+{{e}^{-1.7\cdot a\cdot\left(\theta-b\right)}}}

where θ denotes ability, a denotes item discrimination parameter, b denotes item difficulty parameter and c denotes guessing parameter. For an examinee who has ability θ, P(θ) is the probability of correctly answering a question with parameters a, b, c. D=1.7 is put to make this logistic curve approach to normal ogive curve. (more on this in later posts)

Let’s write a basic code in R for calculating probability of correctly answering a question with (a,b,c) parameters. (I’m assuming R is already installed: )

If we enter the following code to R, it will give the probability of correctly answering an item with parameters (a=1,b=0.5,c=0.2). (You can enter just the parts after “>”, red fonts are outputs, green fonts are comments)

> # First enter parameters
> a = 1
> b = 0.5
> c = 0.2
> theta = 1
> # P is the probability of correctly answering the item with (a,b,c) parameters
> P = c + (1-c) / (1+exp(-1.7*a*(theta-b)))

> P

[1] 0.7604537

As can be seen, P is 0.76. An examinee with ability 1 (θ=1), would answer this item with probability 0.76.

How about examinees with abilities ranging from -3 to 3? For these examinees, what are the probabilities of correctly answering this question?

First, let’s give new values to theta:

> theta = -3:3
> theta
[1] -3 -2 -1 0 1 2 3

Then, let’s calculate the probabilities corresponding to each theta and show it. (I used “cbind” function of R for better presentation)

> P = c + (1-c) / (1+exp(-1.7*a*(theta-b)))
> cbind(theta,P)
theta P
[1,] -3 0.2020793
[2,] -2 0.2112509
[3,] -1 0.2579412
[4,] 0 0.4395463

[5,] 1 0.7604537

[6,] 2 0.9420588

[7,] 3 0.9887491

This seems good. Lets’ graph it with following code:

> plot(theta,P)

Here is the graph:

Graph is fine but it can be more detailed. So, I define more dense theta’s and calculate the probabilities each. Then graph it. Here is the code for doing this:

> # First define new thetas
> # “seq(-3,3,0.1)” means: sequence of numbers from -3 to 3 with interval 0.1.
> theta = seq(-3,3,0.1)
> # Then calculate probabilities
> P = c + (1-c) / (1+exp(-1.7*a*(theta-b)))
> # Then graph thetas on x axis and probabilities on y axis:
> plot(theta,P)

Here is the graph:

All of these are for the item with parameters (a=1,b=0.5,c=0.2). Next, let’s change the item parameters to (a=.8,b=0,c=0.1), and plot the probabilities of correctly answering this item for various thetas.

> # Enter new parameters
> a = 0.8
> b = 0
> c = 0.1
> theta = seq(-3,3,0.1)
> # Calculate probability of correctly answering this new item for all of these thetas
> P = c + (1-c) / (1+exp(-1.7*a*(theta-b)))

This time, I will plot probabilities as a line instead of points. So, I add ‘ type=”l” ‘ to the plot command:

> plot(theta,P,type=”l”)

Here is the plot:

I’d like to see both items on same graph. So, let’s define first item’s parameters as (a1,b1,c1) and probabilities as P1. Similarly, define second item’s parameters as (a2,b2,c2) and probabilities as P2:

> # First items’ parameters
> a1 = 1
> b1 = 0.5
> c1 = 0.2
> P1 = c1 + (1-c1) / (1+exp(-1.7*a1*(theta-b1)))
> # Second items’ parameters
> a2 = 0.8

> b2 = 0

> c2 = 0.1

> P2 = c2 + (1-c2) / (1+exp(-1.7*a2*(theta-b2)))

Next, plot the two graphs. Here, I plot the first item first, then second item on top of it. First item is shown by blue line, second item by red line. Also, I define the limits of y axis as 0 to 1 (by: “ylim=c(0,1)”) to show guessing parameter better. Code for plot:

> plot(theta,P1,col=”red”,type=”l”,ylim=c(0,1))
> lines(theta,P2,col=”blue”)

Here is the plot:


In the next post, I’m planning to create an R function to calculate probabilities and create plots.


Lord, F.M. (1980). Applications of item response theory to practical testing problems: L. Erlbaum Associates Hillsdale, NJ.