[R] IRT Likelihood problem

[Doran, Harold]

I have the following item response theory (IRT) likelihood that I want
to maximize w.r.t. to theta (student ability).

L(\theta) = \prod(p(x))

Where p(x) is the 3-parameter logistic model when items are scored
dichotomously (x_{ij} = 0 or 1) and p(x) is Muraki's generalized partial
credit model when items are scored polytomously (x_{ij} = 0 \ldots J).

Now, I wrote the following two functions to maximize the likelihood. The
first one is for the 3PL and works when all items are scored
dichotmously. The second is for the GPCM and gives the MLE when all
items are scored polytomously.

In the code below, this requires that I first have in hand estimates of
the item parameters so that the maximization is only w.r.t. theta. 

# x = response pattern
# b = vector of location parameters
# a = vector of discrimination parameters 
# c = vector of guessing paramete
theta.3pl <- function(x, b, a, c){
 opt <- function(theta) -sum(dbinom(x, 1, c + ((1-c)/(1 +
exp(-1.7*a*(theta - b)))),  log = TRUE))
 start_val <- log(sum(x)/(length(x)/sum(x)))
 out <- optim(start_val , opt, method = "BFGS", hessian = TRUE)
 out$par
 }

# score = the category the student scored in for item i
# d     = the item parameters for the ith item 
# a     = the discrimination p
# Muraki's GPCM
pcm.max <- function(score, d, a){
    pcm <- function(theta, d, score, a)
exp(sum(a*(theta-d[1:score])))/sum(exp(cumsum(a*(theta-d))))
    opt <- function(theta) -sum(log(mapply(pcm, d, theta = theta, score=
score )))
    start_val <- log(sum(score-1)/(length(score-1)/sum(score-1)))
    out <- optim(start_val, opt, method = "BFGS", hessian = TRUE)
    round(out$par, 2)
  }

However, I have data for which there is a mixture of item types. Some
are dichotomous and others are polytomous. Therefore, I need to somehow
modify these functions to work in a conditional statement that first
evaluates whether the item is dichotomous or not and then uses the right
function to write out and then maximize the likelihood. However, I'm a
bit stumped on how I might code this. Can anyone suggest how that might
work?

For example, assume I have a test consisting of 3 items. The first two
are dichotomous and the last is polytomous. Assume the students score on
these three items is:

x <- c(1,0,3) # that is, 'right', 'wrong', 'scored in category 3'

And further assume the item parameters for these items are

Item 1 c = .11, b = 1.2, a = .58
Item 2 c = .20, b = .65, a = 1.2
Item 3 d = (0, -1.4, -.28, .95)

Now, my function pcm.max also reduces to Master's partial credit model
when a = 1 for all items. And, because Master's partial credit reduces
to the Rasch model when items are scored dichotmously, the only function
I need is pcm.max when working within the Rasch family of models.
However, Muraki's model does not reduce to the 3PL because of the
guessing parameter when items are scored dichotomously, so I think the
problem is slightly more complex in this scenario and requires the use
of both functions.

Last, and a bit tangentially, it is possible when maximizing the 3PL
that I will converge upon a local and not global maximum, a situation
that does not occur with Rasch or 2PL. I only know of two methods for
knowing whether I converged upon the right max. First, plot the
likelihood and look at it or second, use different starting points and
see if I converge to the same max. However, I am maximizing this
likelihood over 80,000 students (or so) and I don't think either of
these two methods are viable. If anyone has a suggestion on how I could
proceed thoughtfully on that issue I welcome that as well.

Many thanks,
Harold