name <- 'old2PL'
par <- c(a = .5, b = -2)
est <- c(TRUE, TRUE)
P.old2PL <- function(par,Theta, ncat){
     a <- par[1]
     b <- par[2]
     P1 <- 1 / (1 + exp(-1*a*(Theta - b)))
     cbind(1-P1, P1)
}
x <- createItem(name, par=par, est=est, P=P.old2PL)
#So, let's estimate it!
dat <- expand.table(LSAT7)
sv <- mirt(dat, 1, c(rep('2PL',4), 'old2PL'), customItems=list(old2PL=x), pars = 'values')
tail(sv) #looks good
mod <- mirt(dat, 1, c(rep('2PL',4), 'old2PL'), customItems=list(old2PL=x))
coef(mod)
mod2 <- mirt(dat, 1, c(rep('2PL',4), 'old2PL'), customItems=list(old2PL=x), method = 'MHRM')
coef(mod2)
#several secondary functions supported
M2(mod, calcNull=FALSE)
itemfit(mod)
fscores(mod, full.scores=FALSE)
plot(mod)
# fit the same model, but specify gradient function explicitly (use of a brower() may be helpful)
gr <- function(x, Theta){
     # browser()
     a <- x@par[1]
     b <- x@par[2]
     P <- probtrace(x, Theta)
     PQ <- apply(P, 1, prod)
     r_P <- x@dat / P
     grad <- numeric(2)
     grad[2] <- sum(-a * PQ * (r_P[,2] - r_P[,1]))
     grad[1] <- sum((Theta - b) * PQ * (r_P[,2] - r_P[,1]))
     ## check with internal numerical form to be safe
     # numerical_deriv(x@par[x@est], mirt:::EML, obj=x, Theta=Theta, type='Richardson')
     grad
}
x <- createItem(name, par=par, est=est, P=P.old2PL, gr=gr)
mod <- mirt(dat, 1, c(rep('2PL',4), 'old2PL'), customItems=list(old2PL=x))
coef(mod, simplify=TRUE)
###non-linear
name <- 'nonlin'
par <- c(a1 = .5, a2 = .1, d = 0)
est <- c(TRUE, TRUE, TRUE)
P.nonlin <- function(par,Theta, ncat=2){
     a1 <- par[1]
     a2 <- par[2]
     d <- par[3]
     P1 <- 1 / (1 + exp(-1*(a1*Theta + a2*Theta^2 + d)))
     cbind(1-P1, P1)
}
x2 <- createItem(name, par=par, est=est, P=P.nonlin)
mod <- mirt(dat, 1, c(rep('2PL',4), 'nonlin'), customItems=list(nonlin=x2))
coef(mod)
###nominal response model (Bock 1972 version)
Tnom.dev <- function(ncat) {
   T <- matrix(1/ncat, ncat, ncat - 1)
   diag(T[-1, ]) <-  diag(T[-1, ]) - 1
   return(T)
}
name <- 'nom'
par <- c(alp=c(3,0,-3),gam=rep(.4,3))
est <- rep(TRUE, length(par))
P.nom <- function(par, Theta, ncat){
   alp <- par[1:(ncat-1)]
   gam <- par[ncat:length(par)]
   a <- Tnom.dev(ncat) %*% alp
   c <- Tnom.dev(ncat) %*% gam
   z <- matrix(0, nrow(Theta), ncat)
   for(i in 1:ncat)
       z[,i] <- a[i] * Theta + c[i]
   P <- exp(z) / rowSums(exp(z))
   P
}
nom1 <- createItem(name, par=par, est=est, P=P.nom, derivType = 'central')
nommod <- mirt(Science, 1, 'nom1', customItems=list(nom1=nom1))
coef(nommod)
Tnom.dev(4) %*% coef(nommod)[[1]][1:3] #a
Tnom.dev(4) %*% coef(nommod)[[1]][4:6] #dRun the code above in your browser using DataLab