[R] Estimating IRT models by using nlme() function

[Dimitris Rizopoulos]

Note that there are several IRT dedicated packages in R, such as the ltm 
and eRm packages. For more info you may check the Psychometrics Task 
View at: http://cran.r-project.org/web/views/Psychometrics.html


I hope it helps.

Best,
Dimitris



On 11/11/2011 9:38 PM, Cengiz Zopluoğlu wrote:
> Hi,
>
> I have a question about estimating IRT models by using nlme, not just rasch
> model, but also other models.
>
> Behavior Research Methods
> <http://www.springerlink.com/content/1554-351x/>  Volume
> 37, Number 2<http://www.springerlink.com/content/1554-351x/37/2/>, 202-218,
> DOI: 10.3758/BF03192688
> Using SAS PROC NLMIXED to fit item response theory models (2005). Ching-Fan
> Sheu<https://mail.google.com/mail/u/0/html/compose/static_files/goog_851944013>,
> Cheng-Te Chen<https://mail.google.com/mail/u/0/html/compose/static_files/goog_851944013>,
> Ya-Hui Su and Wen-Chung Wang
>
> An example of this is provided in the paper above. They use SAS PROC
> NLMIXED, and estimate all possible binary and Polytomous IRT models. I want
> to replicate what they did using R and nlme package. I am able to fit rasch
> model, but have some error messages for 2PL IRT model. So, I could not go
> beyond that.
>
> This is really important for me, because any nonlinear model can be fitted
> by using this approach. I would like to see how well nlme() package
> recovers true parameters, not just regular IRT models but also some other
> less used IRT models.
>
> Here is the example. Let's say I have the following dataset. 1000 people
> responds five items.
>
>
>             ID  Item Response d1 d2 d3 d4 d5
>
>       1.1  1    1        0  1  0  0  0  0
>
>       1.2  1    2        0  0  1  0  0  0
>
>       1.3  1    3        0  0  0  1  0  0
>
>       1.4  1    4        1  0  0  0  1  0
>
>       1.5  1    5        0  0  0  0  0  1
>
>       2.1  2    1        1  1  0  0  0  0
>
>       2.2  2    2        0  0  1  0  0  0
>
>       2.3  2    3        0  0  0  1  0  0
>
>       2.4  2    4        1  0  0  0  1  0
>
>       2.5  2    5        0  0  0  0  0  1
>
>       3.1  3    1        0  1  0  0  0  0
>
>       3.2  3    2        0  0  1  0  0  0
>
>       3.3  3    3        1  0  0  1  0  0
>
>       3.4  3    4        1  0  0  0  1  0
>
>       3.5  3    5        0  0  0  0  0  1
>
>       ..............................
>
>       ..............................
>
> 1000.1 1000    1        1  1  0  0  0  0
>
> 1000.2 1000    2        0  0  1  0  0  0
>
> 1000.3 1000    3        1  0  0  1  0  0
>
> 1000.4 1000    4        1  0  0  0  1  0
>
> 1000.5 1000    5        0  0  0  0  0  1
>
> The items are nested within students. Response is the actual dependent
> variable. d1, d2, d3, d4, and d5 are the corresponding dummy codes for each
> item. We treat persons as random effects, and item parameters as fixed
> effects. To fit Rasch model, I run the following code:
>
>
>    d<- read.table("data.txt", header=TRUE)
> d1<- d$d1
> d2<- d$d2
> d3<- d$d3
> d4<- d$d4
> d5<- d$d5
>   ###########################################################
> onePL<- function(b1,b2,b3,b4,b5,theta) { #nonlinear model to fit
>      b= b1*d1+b2*d2+b3*d3+b4*d4+b5*d5
>      exp((theta-b))/(1+exp((theta-b)))
> }
> #####################################################################
> nlme(model=Response ~ onePL(b1,b2,b3,b4,b5,theta),
> data = d,
> fixed = b1+b2+b3+b4+b5 ~ 1,
> random = theta ~ 1 | ID,
> start=c(b1=0,b2=0,b3=0,b4=0,b5=0),
> )
>
>    *OUTPUT *
>
>   Nonlinear mixed-effects model fit by maximum likelihood
>    Model: Response ~ onePL(b1, b2, b3, b4, b5, theta)
>    Data: d
>    Log-likelihood: -2597.344
>    Fixed: b1 + b2 + b3 + b4 + b5 ~ 1
>          b1         b2         b3         b4         b5
>   -1.1240371  1.5931634  0.2574785 -2.0993236  0.8881437
>   Random effects:
>   Formula: theta ~ 1 | ID
>              theta  Residual
>   StdDev: 0.9765999 0.3780802
>   Number of Observations: 5000
>   Number of Groups: 1000
>
>
> This make sense to me. Actually, the data was simulated data, and the b
> parameters were close to true values used to generate data. Also the
> standard deviation of random effects (theta or latent ability level in this
> case) was close to 1 that was used to generate IRT person ability
> parameters.
>
> Then, I run the following code to estimate 2 PL IRT model. It was all same,
> just add an additional "a" parameter for each item.
>
>
>    twoPL<- function(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,theta) {
>      a= a1*d1+a2*d2+a3*d3+a4*d4+a5*d5
>      b= b1*d1+b2*d2+b3*d3+b4*d4+b5*d5
>      exp(a*(theta-b))/(1+exp(a*(theta-b)))
> }
>   nlme(model=Response ~ twoPL(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,theta),
> data = d,
> fixed = a1+a2+a3+a4+a5+b1+b2+b3+b4+b5 ~ 1,
> random = theta ~ 1 | ID,
> start=c(a1=1,a2=1,a3=1,a4=1,a5=1,b1=0,b2=0,b3=0,b4=0,b5=0),
> )
>
>
> It gives the following error:
>
> *Error in MEEM(object, conLin, control$niterEM) :
>    Singularity in backsolve at level 0, block 1*
>
> Is there anyone who have an idea what I am doing wrong? Is there any error
> in the syntax? I never used nlme package before, so I might be missing some
> components in the code. Or, is this just estimation problem and nlme() can
> not fit this model for some reason?
>
> I would appreciate any help.
>
> Thanks
>
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

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