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

[Doran, Harold]

Just noticed the bias for jml was wrong. It should be

> ### jml bias
> mean(coef(fm3) - trueVals)
[1] -0.173386111948639

> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-
> bounces at r-project.org] On Behalf Of Doran, Harold
> Sent: Tuesday, November 02, 2010 10:34 AM
> To: Paul Johnson; R-SIG-Mixed-Models at r-project.org
> Cc: Dimitris Rizopoulos
> Subject: Re: [R-sig-ME] IRT: discrimination parameters from LMER?
> 
> > 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))))
> }
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models