Estimating IRT models by using nlme() function

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