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

[Dimitris Rizopoulos]

On 11/2/2010 3:33 PM, Doran, Harold wrote:
>> 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.

well, the basic capabilities of ltm are described in the JSS paper 
(http://www.jstatsoft.org/v17/a5/paper), and more information about it 
including R script files with many example can be found in the R-wiki 
webpage: http://rwiki.sciviews.org/doku.php?id=packages:cran:ltm

Best,
Dimitris


>> 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
>

-- 
Dimitris Rizopoulos
Assistant Professor
Department of Biostatistics
Erasmus University Medical Center

Address: PO Box 2040, 3000 CA Rotterdam, the Netherlands
Tel: +31/(0)10/7043478
Fax: +31/(0)10/7043014
Web: http://www.erasmusmc.nl/biostatistiek/