[R-sig-ME] IRT: discrimination parameters from LMER?

[Doran, Harold]

> I'm glad to learn I'm not grossly misunderstanding IRT/lmer. What do
> you think about my observation about the fact that the difficulty
> estimate from ltm seems to be systematically different/better than the
> lmer estimate?

I don't think this is true necessarily. I've done a small experiment all reproducible code is at bottom of email. I simulate some rasch data and then estimate item parameters using ltm, lmer, jml (in MiscPsycho), and mml (a function I wrote that is below using the traditional mml methods). To make a long story short, I compute the bias under each model.

> ### ltm bias
> mean((coef(fm1)[,1] - mean(coef(fm1)[,1])) - trueVals)
[1] -0.173386111948639
> 
> ### lmer bias
> mean((fixef(fm2) - mean(fixef(fm2))) - trueVals)
[1] -0.173386111948639
> 
> ### jml bias
> mean(coef(fm1) - trueVals)
[1] 0.39007242168395
> 
> ### mml bias
> mean((coef(fm4)[1:20] - mean(coef(fm4)[1:20])) - trueVals)
[1] -0.173386111948639

It seems to me that jml (as expected) performs the worst, but all other functions are comparable and seem to recover parameter estimates rather well. If you run the mml function below, it is slow. But, notice that the standard deviation in mml is exactly the same as what lmer gives. That is because a multilevel rasch model with the parameterization I have used (items as fixed respondents as random) is exactly the same as mml. But, lmer's capacity can go much, much further. Had I not constrained the ltm function to have a discrimination of 1, the common discrimination parameter would also be the sd given by mml and lmer.
 
> Do you care to weigh in on the question of whether lmer could ever be
> made to fit a 2PL or 3PL model?  I mean, could I write a family for
> lmer that would do that, or is the calculation required too grossly
> different?

Yes, I'm sure it can. The issue is who will write the code. Doug Bates is extremely generous with his help and his openness to others making contributions to lme4. But, (as he might do so on this thread), he will forewarn you that the underlying code in lmer is very complex and doesn't resemble methods commonly used in other packages. So, because Doug is so amazingly smart, and I am just a mortal, I have decided to steer away from contributing to lmer other than being a "toy tester" for Doug.

I am not a fan of the 3PL for many reasons. First, it is a model that cannot be identified without putting a very slim prior on the guessing parameter. In the end, because the prior is so slim, guess what your parameter estimates are close to? Yes, the prior. So, I just have never been made comfortable with that.

> If 2PL or 3PL is really a whole different structure, I don't mind
> using "ltm."  But defeats the ideal that you proposed in that 2007 JSS
> article (which I really like!), that we might be better off fitting
> specialized models with general purpose software.

I think ltm is a really great package. It may have some capacity to do what you want, I don't fully know its capacities. Perhaps Dimitris can offer some thought. 
 
> Here's a similar example.  Sometimes I have wanted to fit ordinal
> dependent variables, and have learned that those models do not fall
> into the family of GLM and a tool different from lmer is necessary to
> work on them.   I wish it weren't so, but have accepted reality.
> 
> pj

This is very true. I have been doing some work on polytomous models, but am working very slowly on functions to estimate them. There are packages in R that can currently do this. I wonder if MCMCglmm can do it; I don't know.


library(ltm)
library(lmer)
library(MiscPsycho)

Nitems <- 20
Nperson <- 500

set.seed(12345)
dat <- simRasch(Nperson, Nitems) 
itemDat <- dat$data

### Fit data using rasch in LTM
fm1 <- rasch(itemDat, constraint = cbind(ncol(itemDat) + 1, 1))

### Fit using lmer
itemDat$id <- 1:nrow(itemDat)
### reshape into long format
testScores <- reshape(itemDat, idvar='id', varying=list(names(itemDat)[1:Nitems]), v.names=c('answer'), timevar='question', direction='long')

# Treat items as fixed but students as random
fm2 <- lmer(answer ~ factor(question) - 1 + (1|id), testScores, family = binomial(link='logit'))


### Fit using JML
fm3 <- jml(dat$data, bias=TRUE)

### Fit using mml function below
fm4 <- mml(dat$data, Q =10, startVal = coef(fm1)[,1])

### Compute bias
trueVals <- dat$gen

### ltm bias
mean((coef(fm1)[,1] - mean(coef(fm1)[,1])) - trueVals)

### lmer bias
mean((fixef(fm2) - mean(fixef(fm2))) - trueVals)

### jml bias
mean(coef(fm1) - trueVals)

### mml bias
mean((coef(fm4)[1:20] - mean(coef(fm4)[1:20])) - trueVals)



mml <- function(data, Q, startVal = NULL, ...){
	if(!is.null(startVal) && length(startVal) != ncol(data) ){
		stop("Length of argument startVal not equal to the number of parameters estimated")
	}	
	if(!is.null(startVal)){
		startVal <- startVal
		} else {
		p <- colMeans(data)
		startVal <- as.vector(log((1 - p)/p))
	}
	qq <- gauss.quad.prob(Q, dist = 'normal', mu = 0, sigma=1)
	rr1 <- matrix(0, nrow = Q, ncol = nrow(data))
	data <- as.matrix(data)
	L <- nrow(data)
	C <- ncol(data)
	u <- 0
	fn <- function(b){
	s <- b[C+1]
	b <- b[1:C]	
		for(j in 1:Q){
			for(i in 1:L){
				rr1[j,i] <- exp(sum(dbinom(data[i,], 1, plogis(qq$nodes[j]-b), log = TRUE))) * 
				((1/(s*sqrt(2*pi)))  * exp(-((qq$nodes[j]-u)^2/(2*s^2))))/dnorm(qq$nodes[j]) * qq$weights[j]
			}
		}
	-sum(log(colSums(rr1)))
	}	
	opt <- optim(c(startVal,1), fn, method = "BFGS", hessian = TRUE)
	list("coefficients" = opt$par, "LogLik" = -opt$value, "Std.Error" = sqrt(diag(solve(opt$hessian))))
}