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.

Relating Classics Test Theory Parameters to IRT Parameters

* Relating Classics Test Theory Parameters to IRT Parameters

* This post follows some of the discussion in Multidimensional Item Response Theory by Mark Reckase's chapter 2.

* First let's relate Classical Test theory ideas of difficulty to that of IRT parameter b - difficulty.

* We will use for our underlying true DGP the three parameter logistic model.

* (2.13) P(U=u) = c + (1-c) * exp(a(t - b))/(1 + exp(a(t - b)))

* That is, the probability of getting the problem right, is a function of c (the guessing parameter of the item), a (the discriminatory parameter of the item), b (the difficulty parameter of the item), and t (the ability of the test taker).

* We want to equate this to the classical test theory idea of difficult.  In classical test theory the probability of getting a item correct is the difficulty of the item.

* Let us imagine a heterogenous group of 1000 students

clear
set obs 1000

* Create a student ID
gen stud_id = _n

gen t = rnormal()

* Now let's see how we can compare classical difficulties (CD) to IRT difficulty parameter b.

* Let's imagine that all of our students test the same test with 100 items that range in b value which is independent of the choices of parameters a and c.

* Let's have all of our data listed vertically.

* Create 200 test items for each student
expand 200

* Give each item a different ID
bysort stud_id: gen item_id = _n

* There are many ways to make sure all of the items have the same parameters.  I will use a for loop though generating a seperate data set for all of the items and merging it in would be another good way or drawing all of the parameters from distributions then taking the average accross all of the items of the same ID would probably be the most efficient code wise but would make it difficult to specify exact distributional parameters.
gen a = .
gen b = .
gen c = .

qui forv i = 1/200 {
  * This will only draw one random variable for each local macro
  local a = runiform()/2+.4
  local b = rnormal()*2
  local c = runiform()/4

  * This will assign that draw to the item `i'
  replace a = `a' if item_id==`i'
  replace b = `b' if item_id==`i'
  replace c = `c' if item_id==`i'
}

* Now let's generate the probability of getting that problem correct given the parameter values and the student t scores.

gen P = c + (1-c) * exp(a*(t - b))/(1 + exp(a*(t - b)))

* If we try to do a direct scatterplot then we are overwhelmed.

* Instead we want to know the probability of a correct answer for each item (given the population being tested).

* Let's us first preserve the current state of our data.
preserve

* So we collapse the data set to item level.

* The default of the collapse command is to take the mean.
collapse a b c P t, by(item_id)
* I included the t value just as a debugging test.  t should be constant accross all items.

label var b "IRT b"
label var P "Difficulty (Probability of Correctly Answer)"

scatter P b

* The reason things start fanning out as b gets large is due to the guessing parameter c.

* Even when b is so large that the probability of getting the answer correct based on knowledge is close to zero there is still the chance of guessing the correct answer.

restore

* In order to compare the discrimination parameter to classical test theory we will look at the Point Biserial correlation.  Which is the correlation between a the responses to an item on the test and the total test score.

* First we need to draw actual item reponses

gen u = rbinomial(1, P)

* Now let's generate total test scores

bysort stud_id: egen total_score = sum(u)

* Now let's generate the point-biserial values for our items

gen pbiserial = .

qui forv i = 1/200 {
  corr u total_score if item_id == `i'
  replace pbiserial = r(rho) if item_id == `i'
}

preserve
collapse a pbiserial, by(item_id)

label var a "IRT discriminatory parameter (a)"
label var pbiserial "Classical test theory point-biserial correlation"
twoway (lfitci pbiserial a)  (scatter pbiserial a)

restore


* In order to approximate c using just classical test scores we will look at the lowest 10% of students in terms of total scores and see how they perform on each question on average.

xtile score_pct = total_score, nquantiles(10)
* This should have created 5 groups ranked accounting to total_score

bysort item_id: egen c_hat = mean(u) if score_pct == 1

preserve
collapse c c_hat, by(item_id)

label var c "IRT guessing parameter (c)"
label var c_hat "A guess at the guessing parameter"
twoway (lfitci c_hat c)  (scatter c_hat c)

restore

* We can see there is some relationship between our guess at the guessing parameter using the item responses for the lowest 10% of students and the true guessing parameter.

* The problem with this graph is that items have different difficulties



* Finally let's look at out estimates of student ability t relative to that of their total test score

preserve
collapse total_score t, by(stud_id)

label var total_score  "Total test score"
label var t "IRT univeratiate ability (t)"
twoway (lfitci total_score t)  (scatter total_score t)

restore

* It seems that total test score provides a reasonable linear approximation for student ability even when ability is drawn using a IRT data generating process.

* The largest advantage of IRT relative to that of classical test theory total test performance estimators is the external applicability of the estimates.  IRT is supposed to predict future performance on different tests.  While, classical test theory only predicts performance on the same or similar tests (if I understand this properly).

Why Item Response Theory is very Cool!

Currently I am reading Multidimensional Item Response Theory by one of my mentors Mark Reckase and I realize how very cool item reponse theory (IRT) is.

The underlying purpose of item response theory is measure two things simultaneously.

1. Measure the item's parameters.  This means measure the characteristics of the test item such as difficulty of the item, likelihood that someone will correctly guess the item, and typically the discriminatory power of the item, ie. how much power the item has to identify the difference in ability between one student and another student.

2. While attempting to measure the items' parameters on a test, IRT methods must also estimate the ability of the students taking the test simultaneously.  That is to say that each student has a different level of ability when entering the test and generally that ability level is unknown a-priori.  Of course after taking the test the ability level is still not "known" but at least some reasonable approximation of the ability level of the student should be known at that time.

So why is this cool?  Well, imagine having a data set with K items (variables) and N students (observations). Now, only knowing that either the students got the answers correct or they failed do an estimation which is largely similar to a logit model except that your exogenous variables (Xs) also need to be estimated.

At this point many of you will probably say, well, this is clearly possible if difficult you are willing to make some identifying assumptions.

So far, the only identifying assumptions that I see needed are 1. The probability that a student answers a question correctly is greater the greater the ability of the student, dP(Y=1|X,T,D)/dX>0. And 2. The probability that the question is answered incorrectly is greater the greater the difficulty of the  item dP(Y=1|X,T,D)/dD<0 .="." and="and" assumption="assumption" form="form" know="know" make="make" model.="model." nbsp="nbsp" of="of" p="p" structural="structural" that="that" the="the" underlying="underlying" we="we" well="well">
I will be posting more on IRT as I learn.

Calculate an Empirical CDF for a random variable

# First lets define the function. ecdf (Empirical Cumulative Distribution Function)
ecdf <- function(x, bins=100) {

  # We will sort the imput variable from smallest to largest
  x_sort = sort(x)

  # Now we will create a percentile variable that will rank x from slightly above zero to slightly below 1.
  pct =  (1:length(x)/(length(x)+1))

  # We will now ceate a number of groups equal to length(x)/bins
  grps = ceiling(pct*bins)

  # First let's create an empty vector of zeros with length = length(pct) = length(x)
  avpct = 0*pct
  # Next, let's replace avpct with the mean of pct for each group.
  for (i in 1:bins) avpct[grps==i] = mean(pct[grps==i])

  # Create a bunch of empty value to hold the CDF values
  cdf_bins = cdf_pcts = cdf_values = 1:bins

  # Now loop through each bin value
  for (i in 1:bins) {
    cdf_pcts[i] = mean(pct[grps==i])
    cdf_values[i] = mean(x_sort[grps==i])
  }

  # Because there is no
  data.frame(bins=cdf_bins, pcts=cdf_pcts, xvar=cdf_values)

}

# We will first generate some data, y is 4000 normal draws divided by a uniform(.1,1.1) draw.
y <- rnorm(4000)/(runif(4000)+.1)

cdf1 = ecdf(y, 100) # Calculate an Empirical CDF for variable y with 100 bins
cdf1

# We can see the ECDF is pretty steep around the mean 0.  This is because most of the observations fall between -10 and 10, yet a very few are either large (around 30) or small (around -30), forcing the x dimension to be large.
plot(cdf1$xvar, cdf1$pcts, type="l")

# In a future post I hope to use the empirical CDF information to calculate an inverse CDF that can be used to draw correlated non-normally distributed random variables.

LIE: Law of Iterative Expectations, Word Problem

The law of iterative expectations is extremely useful and I have been trying to think of ways of explaining it.

It basically state that E(Y)=E(E(Y|X))=E(E(Y|X,Z)) where X and Z are sets of covariates.

See my previous post for more information: Law of Iterative Expectations/Law of Total Expectations

I have been trying to explain it in a manner that is intuitively appealing.  Often times I find notation to be a little difficult to dissect.

Imagine Bob and Susie are brother and sister, they go to the same school and always buy lunch.  Their parents randomly decide how much to give them for lunch each day.  But their parents like Susie more than Bob so however much money they give Bob, Susie expects to get a larger amount of money.  There are a few things that you might want to find out.

a. What is the expected amount of money that Bob E(X) will get?
b. What is the expected amount of money that Susie E(Y) will get?
c. What is the distribution of lunch money to Bob (pdf(X))?
d. How much money does Susie expect to get given that we know how much Bob got E(Y|X)?

We know already from the setup that E(Y)>E(X).  {This is not important for the example}

If we know a and b, can we infer c or d? No because in general E(Y)=f(E(X)) and unless we make some assumption about f, then cannot identify the relationship between a and b.

What about if we knew a,b,c? No.  Imagine, if we had a simple distribution p=1/2, x=1: p=1/2, x = 5, thus E(X) = 3, because 1/2*1+1/2*5 = 3.  Imagine that we know E(Y)=6.  This still does not infer the relationship d.  Why, because E(Y|X)=2*X or E(Y|X)=X+3 or any other of a set of infinite functional forms.

* So what can we infer?

* 1. If we know the distribution of X (c), then we know the expected value of X (a).  This is because the expected value function is defined as only the integration of x across the pdf of x.

* 2. If we know the distribution of X (how frequently Bob gets each amount of money) and we know the expected relationship between how much money Bob gets and how much Susie gets (d) then we can figure out how much money Susie gets.  How does this work?  Well imagine that Bob gets $2 with probability 1/3, $4 with probability 1/3, and $6 with probability 1/3.  Imagine also that Susie expects to gets the squared amount of whatever Bob got E(Y|X)=X^2.  Thus we can figure out how much Susie expects to get on any day without knowing how much Bob has gotten yet E(Y).  E(Y)=E(E(Y|X)) = 1/3*2^2 + 1/3*4^2 + 1/3*6^2 = (56)/3 ~ 18.6.

Matrix Operations in Mata

* This post demonstrates a few methods for how to input matrices into Mata and how to do some basic matrix operations.  Mata is a matrix programming language so basic matrix operations are extremely easy.

mata
// Let's first build some matrices in Mata
// They can be built directly
A = ( 2 , -1, 5 \ 3, 0 , -1 \ 3, 3 , -1)
A

B1 = (0 , 3 ,-1)
B2 = (3 , 0 , 0)
B3 = (0 , 2 , 0)

// Or through a combination of vectors
B = (B1 \ B2 \ B3)

B

// We can also start with an empty matrix and fill in values

C = J(3,2,4)
// The J command creates a matrix with 3 rows, 2 colums, and with default values of 4

C
// We can replace individual elements once the matrix is created
C[1,2] = 3
C[3,1] = 7

// Or entire submatrices
C[2,] = (2,6)
C

// Now let's see how various matrix operations perform in Mata:

// a. AB
A*B

// b. BA
B*A

// c. A+B
A+B

// d. A'B'
A'*B'

// e. B'A'
B'*A'

// f. A'B
A'*B

// g. (AB)'
(A*B)'

// h. ABC
A*B*C

// i. C'AC
C'*A*C

// j. CAC'
C*A*C'
// This does not exist because C*A cannot be multiplied as C is 3x2 and A is 3x3

// k. trace(C'AC)
trace(C'*A*C)

// l. trace(CAC')
trace(C*A*C')
// Likewise the trace does not exist

end

* Matrices can also be input into Mata from data sets.

clear
set obs 5
gen y = 3
replace y = 6  if _n == 2
replace y = 10 if _n == 3
replace y = 8  if _n == 4
replace y = 2  if _n == 5

mata
// The command st_data retrieves data from stata variables to be used in Mata
y = st_data(. , "y")
y

// Once input, matrices can easily be manipulated
y*y'

// The square of the norm of y
norm(y)^2

// This happens to be identical to:
sum(y:^2)

// The norm is the equclidean norm which is the square square root of the sum of all of the squares of a vector.
// Thus the square of it is just the sum of the squares.

// As should be clear, manipulating matrices in Mata is extremely easy.
// Thus Stata is able to pack a powerful Matrix programming language inside an effective high level user language.

// The largest frustration that I have had with Mata is the relative quality of the documentation.
// I find the documentation of Mata much harder to use than that of Stata (at least in version 11, perhaps version 12 has better documentation).
end

Mata speed gains over Stata

* The inclusion of Mata as an available alternative programming language for Stata users was a great move by Stata.

* Mata in general runs much quicker than programming on the surface level in Stata.

* In Stata each loop that runs is compiled (interpretted into machine code) as it runs creating a lot of work for the machine.

* In Mata on the other hand, the entire loop is compiled prior to running.

* Let's see how this works.

* Let's say we want to add up the square of the numbers 1 through 100000

* Method 1: Surface loop

timer clear 1
timer on 1
local x2 = 0

forv i = 1/1000000 {
  local x2 = `x2'+`i'^2
}

di `x2'

timer off 1
timer list 1

* On my laptop, this takes about 13.5 seconds

* Method 2: Mata loop
timer clear 1
timer on 1
mata
  x2=0
  // This command can be read as start i at 1,
  // keep looping so long as i is less than 1000000,
  // the third argument looks a little fishy but it is syntax
  // that has been around for a while (at least since C).
  // It would be identical to writing i=i+1, in other words, add 1 to i.
  // Following the for loop we can immediately place a since line command.
  for (i = 1; i <= 1000000; i++) x2=x2+i^2
  // If there is nothing done with the value x2 then mata displays this value.
  // R handles this identically
  x2
end

timer off 1
timer list 1

* In contrast, my computer completed the loop using mata in .27 seconds, many magnitudes of speed faster.

* However this does not mean you need to learn to use mata (since it has its own limitations and syntax) in order to speed up your commands.

* Method 3: Use Stata's data structure to accomplish vector tasks
timer clear 1
timer on 1

clear
set obs 1000000
gen x2 = _n^2

* The sum command will calculate the mean of x2 which is the same as the sum of x2 divided by it's number of observations.
sum x2
* We can reverse that operation easily.
di r(N)*r(mean)

timer off 1
timer list 1
* Using a little knowledge of how Stata stores post command information this method does the same trick in .2 seconds

* Method 4: The speed gains in 3 was as a result of using the vector structure of data columns.  Mata can do very similar things even easier.

timer clear 1
timer on 1
// This command looks a little fishy, but it is easy to understand.
// Order of operations must be taken into account.
// First the 10^6 is evaluated which equals 1000000
// Then the vector 1..10^6 is made which looks like 1 2 3 ... 1000000
// The .. tells mata to make a count vector.
// If I had written :: then mata would have made a column vector instead.
// Once the vector is made then the command :^2 tells stata to do a piece wise squaring of each term in the vector.
// Finally the sum command adds all of the elements of the vector together to generate the result we were looking for.
mata: sum((1..10^6):^2)
timer off 1
timer list 1
* The result is that this command only took .04 seconds to run through efficient coding in Mata.

# As a matter of comparison, this command
system.time(sum((1:10^6)^2))
# took .04 seconds in R

# And the loop:
x=0
system.time(for(i in 1:10^6) x=x+i^2)
# 1.3 seconds

# Thus Mata in this example is significantly faster than Stata and about the same speed as R.