It's not really clear what is going on in this paper as the authors seem to have breezed over some important information. From what it appears, they seem to want to estimate all the parameters in IRT models through a single call to R's optim
function, but unless further constrains are supplied to the optimization process this simply will not work (a similar approach could be used for linear factor analysis, and will suffer the same problem).
In the typical approach to using logistic regression in IRT, the latent trait values ($\theta$) are realized by simply using the (potentially rescaled) total/sum score as a proxy for each participants ability, $\theta_i$. Admittedly this is a bit of a hack, but for bias detection methods in unidimensional tests where a grouping variable is included it has seen some success.
Here the authors are trying at a joint maximum-likelihood approach but do not seem to be imposing any additional constraints on the distribution of the ability parameters. Hence, the optimizer really won't really know where to go and essentially just returns the starting values (i.e., the mean of the abilities could arbitrarily increase to explain successful item endorsement, or alternatively the difficulty of the items could decrease to indicate that the items were just really easy. Both are valid because the model is not sufficiently identified).
Here's some R code to try and reproduce what the authors are doing. Note that the cost function only seems to change the value of the objective function.
library(mirt)
dat <- key2binary(SAT12,
key = c(1,4,5,2,3,1,2,1,3,1,2,4,2,1,5,3,4,4,1,4,3,3,4,1,3,5,1,3,1,5,4,5))
ProbFun <- function(X, theta, tol = 1e-10){
Z <- X %*% t(theta)
ret <- t(1 / (1 + exp(-Z)))
ret <- ifelse(ret < tol, tol, ret)
ret <- ifelse(ret > 1 - tol, 1 - tol, ret)
ret
}
L <- function(p, dat, lambda = 0){
X <- matrix(p[1:(ncol(dat))], ncol(dat))
theta <- cbind(1, matrix(p[-c(1:ncol(dat))], nrow(dat)))
if(ncol(thetas) > ncol(X)) X <- cbind(X, 1)
PP <- ProbFun(X, theta)
-sum(log(PP) * dat + log(1 - PP) * (1 - dat)) + sum(lambda * p^2)
}
ability <- scale(rowSums(dat))
intercepts <- qlogis(colMeans(dat))
p <- c(intercepts, ability)
optim(p, L, dat=dat, lambda=0)
optim(p, L, dat=dat, lambda=.1)
optim(p, L, dat=dat, lambda=.5)
optim(p, L, dat=dat, lambda=1)
# marginal ML run
mod <- mirt(dat, 1, 'Rasch')
coef(mod, simplify=TRUE)