R/02b-item_methods.R
In philchalmers/mirt: Multidimensional Item Response Theory
Defines functions
dP.ull
P.ull
d2P.cll
dP.cll
P.cll
d2P.5PL
dP.5PL
P.5PL
P.ggum
P.lca
P.ideal
P.nestlogit
P.comp
nominalParDeriv
P.nominal
P.rsm
P.gpcm
P.rating
P.poly
P.mirt
Use_R_Deriv
Use_R_ProbTrace
ordinal_itemtypes
Valid_iteminputs
Experimental_itemtypes
# valid itemtype inputs
# flag to indicate an experimental item type (requires an S4 initializer in the definitions below)
# note: cannot match Valid_iteminputs
Experimental_itemtypes <- function() c('experimental', 'grsmIRT', 'fivePL', 'cll', 'ull')
Valid_iteminputs <- function() c('Rasch', '2PL', '3PL', '3PLu', '4PL', '5PL', 'CLL', 'ULL',
'graded', 'grsm', 'gpcm', 'gpcmIRT',
'rsm', 'nominal', 'PC2PL','PC3PL', '2PLNRM', '3PLNRM', '3PLuNRM', '4PLNRM',
'ideal', 'lca', 'spline', 'monopoly', 'ggum', 'sequential', 'Tutz', Experimental_itemtypes())
ordinal_itemtypes <- function() c('dich', 'fivePL', 'graded', 'gpcm', 'sequential', 'cll', 'ull',
'ggum', 'rating', 'spline', 'monopoly',
'partcomp', 'rsm', 'ideal', 'gpcmIRT', 'grsmIRT')
# Indicate which functions should use the R function instead of those written in C++
Use_R_ProbTrace <- function() c('custom', 'spline', 'sequential', 'Tutz', Experimental_itemtypes())
Use_R_Deriv <- function() c('custom', 'rating', 'partcomp', 'nestlogit', 'spline', 'sequential', 'Tutz',
Experimental_itemtypes())
#--------------------------------------------------------------------
# Item model definitions
setClass("dich", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'dich'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'dich'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'dich'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'dich'),
definition = function(x){
par <- x@par
d <- par[x@nfact+1L]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'dich'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
rnorm(1L),
rnorm(1L, -2, .5),
rnorm(1L, 2, .5))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'dich'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.mirt <- function(par, Theta, asMatrix = FALSE, ot = 0)
{
return(.Call("traceLinePts", par, Theta, ot))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'dich', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
P <- P.mirt(x@par, Theta=Theta, ot=ot)
return(P)
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'dich', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call('dparsDich', x@par, Theta, estHess, x@dat, offterm)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'dich', Theta = 'matrix'),
definition = function(x, Theta){
N <- nrow(Theta)
nfact <- ncol(Theta)
parlength <- length(x@par)
u <- antilogit(x@par[parlength])
g <- antilogit(x@par[parlength - 1L])
d <- x@par[parlength - 2L]
a <- x@par[seq_len(nfact)]
Pstar <- P.mirt(c(a, d, -999, 999), Theta)[,2L]
grad <- hess <- vector('list', 2L)
grad[[1L]] <- grad[[2L]] <- hess[[1L]] <- hess[[2L]] <- matrix(0, N, nfact)
for(i in seq_len(nfact)){
grad[[2L]][ ,i] <- (u-g) * a[i] * (Pstar * (1 - Pstar))
grad[[1L]][ ,i] <- -1 * grad[[2L]][ ,i]
hess[[2L]][ ,i] <- 2 * (u - g) * a[i]^2 * ((1 - Pstar)^2 * Pstar) -
(u - g) * a[i]^2 * (Pstar * (1 - Pstar))
hess[[1L]][ ,i] <- -1 * hess[[2L]][ ,i]
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'dich', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta)
PQ <- apply(P, 1L, prod)
ret <- cbind(Theta * PQ, PQ, P)
ret
}
)
# ----------------------------------------------------------------
setClass("graded", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'graded'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'graded'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'graded'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'graded'),
definition = function(x){
par <- x@par
d <- par[-seq_len(x@nfact)]
d
}
)
setMethod(
f = "CheckIntercepts",
signature = signature(x = 'graded'),
definition = function(x){
z <- ExtractZetas(x)
return(all(sort(z, decreasing = TRUE) == z))
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'graded'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'graded'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.poly <- function(par, Theta, itemexp = FALSE, ot = 0)
{
return(.Call('gradedTraceLinePts', par, Theta, itemexp, ot, FALSE))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'graded', Theta = 'matrix'),
definition = function(x, Theta, itemexp = TRUE, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.poly(x@par, Theta=Theta, itemexp=itemexp, ot=ot))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'graded', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call("dparsPoly", x@par, Theta, offterm, x@dat,
length(x@par) - ncol(Theta), estHess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'graded', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
P <- ProbTrace(x, Theta, itemexp = FALSE)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(ncol(P)-1L)){
w1 <- P[,i] * (1-P[,i]) * a[j]
w2 <- P[,i+1L] * (1-P[,i+1L]) * a[j]
grad[[i]][ ,j] <- w1 - w2
hess[[i]][ ,j] <- a[j]^2 * (2 * P[ ,i] * (1 - P[,i])^2 -
P[ ,i] * (1 - P[,i]) -
2 * P[ ,i+1L] * (1 - P[,i+1L])^2 +
P[ ,i+1L] * (1 - P[,i+1L]))
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'graded', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta, itemexp=FALSE)
P <- P[,-c(1, ncol(P)), drop=FALSE]
PQ <- P * (1 - P)
ret <- cbind(Theta * rowSums(PQ), PQ)
ret
}
)
# ----------------------------------------------------------------
setClass("rating", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'rating'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'rating'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'rating'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'rating'),
definition = function(x){
par <- x@par
d <- par[-c(seq_len(x@nfact), length(par))]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'rating'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE),
rnorm(1L))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'rating'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.rating <- function(par, Theta, itemexp = FALSE, ot = 0)
{
return(.Call('gradedTraceLinePts', par, Theta, itemexp, ot, TRUE))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'rating', Theta = 'matrix'),
definition = function(x, Theta, itemexp = TRUE, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.rating(x@par, Theta=Theta, itemexp=itemexp, ot=ot))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'rating', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
hess <- matrix(0, length(x@par), length(x@par))
dat <- x@dat
nfact <- x@nfact
a <- x@par[seq_len(nfact)]
d <- ExtractZetas(x)
nzetas <- length(d)
shiftind <- length(x@par)
shift <- x@par[shiftind]
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
par <- c(a, d + shift)
P <- P.poly(par=par, Theta=Theta, ot=offterm)
ret <- .Call("dparsPoly", par, Theta, offterm, x@dat,
length(d), estHess)
grad <- ret$grad
hess <- ret$hess
hess <- cbind(hess, rep(0, nrow(hess)))
hess <- rbind(hess, rep(0, ncol(hess)))
dc <- numeric(1L)
Pfull <- P
PQfull <- Pfull * (1-Pfull)
P <- P.poly(c(a, d + shift), Theta, itemexp=TRUE, ot=offterm)
rs <- dat
for(i in seq_len(ncol(rs)))
dc <- dc + rs[,i]/P[,i] * (PQfull[,i] - PQfull[,i+1L])
dc <- sum(dc)
grad <- c(grad, dc)
if(estHess){
cind <- ncol(hess)
ddc <- numeric(nrow(P))
for(i in seq_len(ncol(rs)))
ddc <- ddc + rs[,i]/P[,i] * (Pfull[,i] - 3*Pfull[,i]^2 + 2*Pfull[,i]^3 -
Pfull[,i+1L] + 3*Pfull[,i+1L]^2 - 2*Pfull[,i+1L]^3) -
rs[,i]/P[,i]^2 * (PQfull[,i] - PQfull[,i+1L])^2
hess[cind, cind] <- sum(ddc)
for(i in seq_len(nzetas))
hess[cind, nfact + i] <- hess[nfact + i, cind] <-
sum((rs[,i]/P[,i] * (-Pfull[,i+1L] + 3*Pfull[,i+1L]^2 - 2*Pfull[,i+1L]^3) -
rs[,i]/P[,i]^2 * (PQfull[,i] - PQfull[,i+1L]) * (-PQfull[,i+1L]) +
rs[,i+1L]/P[,i+1L] * (Pfull[,i+1L] - 3*Pfull[,i+1L]^2 + 2*Pfull[,i+1L]^3) -
rs[,i+1L]/P[,i+1L]^2 * (PQfull[,i+1L] - PQfull[,i+2L]) * (PQfull[,i+1L])))
for(j in seq_len(nfact)){
tmp <- 0
for(i in seq_len(ncol(rs)))
tmp <- tmp + (rs[,i]/P[,i] * Theta[,j] *
(Pfull[,i] - 3*Pfull[,i]^2 + 2*Pfull[,i]^3 -
Pfull[,i+1L] + 3*Pfull[,i+1L]^2 - 2*Pfull[,i+1L]^3) -
rs[,i]/P[,i]^2 * (PQfull[,i] - PQfull[,i+1L]) * Theta[,j] *
(PQfull[,i] - PQfull[,i+1L]))
hess[cind, j] <- hess[j, cind] <- sum(tmp)
}
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'rating', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
P <- ProbTrace(x, Theta, itemexp = FALSE)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(ncol(P)-1L)){
w1 <- P[,i] * (1-P[,i]) * a[j]
w2 <- P[,i+1L] * (1-P[,i+1L]) * a[j]
grad[[i]][ ,j] <- w1 - w2
hess[[i]][ ,j] <- a[j]^2 * (2 * P[ ,i] * (1 - P[,i])^2 -
P[ ,i] * (1 - P[,i]) -
2 * P[ ,i+1L] * (1 - P[ ,i+1L])^2 +
P[ ,i+1L] * (1 - P[ ,i+1L]))
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'rating', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
setClass("gpcm", contains = 'AllItemsClass',
representation = representation(mat='logical'))
setMethod(
f = "print",
signature = signature(x = 'gpcm'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'gpcm'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'gpcm'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'gpcm'),
definition = function(x){
par <- x@par
d <- par[(length(par)-x@ncat+1L):length(par)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'gpcm'),
definition = function(x){
ns <- sum(grepl('ak', x@parnames))
par <- c(rlnorm(x@nfact, meanlog=-.2, sdlog=.5), rep(NA, ns), 0,
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'gpcm'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.gpcm <- function(par, Theta, ot = 0, mat = FALSE, returnNum = FALSE)
{
return(.Call("gpcmTraceLinePts", par, Theta, ot, FALSE, mat, returnNum))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'gpcm', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.gpcm(x@par, Theta=Theta, ot=ot, mat=x@mat))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'gpcm', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call("dparsNominal", x, Theta, offterm, FALSE, estHess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'gpcm', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
d <- ExtractZetas(x)
ak <- 0:(x@ncat - 1L)
P <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta)
Num <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta, returnNum = TRUE)
Den <- rowSums(Num)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(x@ncat)){
grad[[i]][ ,j] <- ak[i] * a[j] * P[ ,i] - P[ ,i] * (Num %*% (ak * a[j])) / Den
hess[[i]][ ,j] <- ak[i]^2 * a[j]^2 * P[ ,i] -
2 * ak[i] * a[j] * P[,i] * (Num %*% (ak * a[j])) / Den +
2 * P[,i] * ((Num %*% (ak * a[j])) / Den)^2 -
P[,i] * ((Num %*% (ak^2 * a[j]^2)) / Den)
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'gpcm', Theta = 'matrix'),
definition = function(x, Theta){
nfact <- ncol(Theta)
ncat <- x@ncat
num <- P.gpcm(x@par, Theta=Theta, returnNum = TRUE, mat = x@mat)
den <- rowSums(num)
P <- num/den
a <- x@par[seq_len(ncol(Theta))]
if(!x@mat){
ak <- matrix(x@par[(ncol(Theta)+1L):(ncol(Theta)+x@ncat)], nfact, ncat, byrow=TRUE)
} else {
ak <- matrix(x@par[(nfact+1L):(nfact+ncat*nfact)], nfact, ncat, byrow=TRUE)
}
dp <- matrix(0, nrow(Theta), length(x@par))
aknum <- eak <- vector('list', nfact)
e <- 0:(x@ncat-1)
for(i in seq_len(nfact)){
aknum[[i]] <- t(ak[i,] * t(num))
eak[[i]] <- e*ak[i,]
}
for(i in seq_len(nfact)){
for(j in 2L:ncat)
dp[,i] <- dp[,i] + eak[[i]][j]*Theta[,i]*P[,j] -
e[j]*P[,j]*rowSums(Theta[,i] * aknum[[i]])
}
for(j in seq_len(x@ncat)){
dp[,nfact + ncat + j] <- e[j] * P[,j] - e[j] * P[,j]^2 -
as.numeric(e[-j] %*% t(P[,-j]*P[,j]))
}
return(dp)
}
)
# ----------------------------------------------------------------
setClass("rsm", contains = 'AllItemsClass',
representation = representation(mat='logical'))
setMethod(
f = "print",
signature = signature(x = 'rsm'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'rsm'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'rsm'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'rsm'),
definition = function(x){
par <- x@par
d <- par[(length(par) - x@ncat):length(par)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'rsm'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5), 0:(x@ncat-1), 0,
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE),
rnorm(1L))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'rsm'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.rsm <- function(par, Theta, ot = 0, mat = FALSE, returnNum = FALSE)
{
return(.Call("gpcmTraceLinePts", par, Theta, ot, TRUE, mat, returnNum))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'rsm', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.rsm(x@par, Theta=Theta, ot=ot))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'rsm', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
dat <- x@dat
# nfact <- x@nfact
# nzetas <- ncol(dat)
a <- ExtractLambdas(x)
d <- ExtractZetas(x)
shift <- d[length(d)]
dshift <- d <- d[-length(d)]
dshift[-1L] <- d[-1L] + shift
ak <- 0:(length(d)-1L)
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
# P <- ProbTrace(x=x, Theta=Theta, useDesign = FALSE, ot=offterm)
tmp <- .Call("dparsNominal", x, Theta, offterm, TRUE, estHess)
# num <- P.nominal(c(a, ak, dshift), ncat=length(ak),
# Theta=Theta, returnNum=TRUE, ot=offterm)
grad <- tmp$grad
hess <- tmp$hess
#quick calcs for derivs
# nfact <- length(a)
# ncat <- length(d)
# akind <- nfact
# dind <- nfact + ncat*2 + 1L #go backwards
# ak2 <- ak^2
# P2 <- P^2
# P3 <- P^3
# aTheta <- as.vector(Theta %*% a)
# aTheta2 <- aTheta^2
# dat_num <- dat/num
# numsum <- rowSums(num)
# numD <- num %*% c(0, rep(1, ncol(num)-1L))
# numak <- matrix(num %*% ak, nrow(Theta), ncol(Theta))
# numakThetaD <- numak * Theta
# numD2 <- num %*% c(0, rep(1, ncol(num)-1L))
# numakThetaD2 <- numak * Theta
# ak0 <- ak
# ak0[1L] <- 0
# cind <- length(grad)
# if(estHess){
# tmp <- 0
# for(i in 1L:nzetas)
# tmp <- tmp + dat[,i]*numD^2 / numsum^2 - dat[,i]*numD2/numsum
# hess[cind, cind] <- sum(tmp)
# for(j in 1L:nzetas){
# tmp <- 0
# for(i in 1L:nzetas)
# tmp <- tmp + dat[,i]*P[,j]*numD/numsum - dat[,i]*P[,j]
# hess[cind, nfact+j] <- hess[nfact+j, cind] <- sum(tmp)
# }
# for(j in 1L:nfact){
# tmp <- 0
# for(i in 1L:nzetas)
# tmp <- tmp + dat[,i]*numD*numakThetaD[,j]/numsum^2 -
# dat[,i]* (num %*% ak0*Theta[,j])/numsum
# hess[cind, j] <- hess[j, cind] <- sum(tmp)
# }
# }
####
#TODO - can't seem to get the last value of the gradient quite right for some reason....
x2 <- x
x2@est <- c(rep(FALSE, length(x2@est)-1L), TRUE)
grad[x2@est] <- numerical_deriv(x@par[x2@est], EML, obj=x2, Theta=Theta)
if(estHess && any(x@est)){
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
}
####
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
ret
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'rsm', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
d <- ExtractZetas(x)
t <- d[length(d)]
d <- d[-length(d)]
d[-1L] <- d[-1L] + t
ak <- 0:(x@ncat - 1L)
P <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta)
Num <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta, returnNum = TRUE)
Den <- rowSums(Num)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(x@ncat)){
grad[[i]][ ,j] <- ak[i] * a[j] * P[ ,i] - P[ ,i] *
(Num %*% (ak * a[j])) / Den
hess[[i]][ ,j] <- ak[i]^2 * a[j]^2 * P[ ,i] -
2 * ak[i] * a[j] * P[,i] * (Num %*% (ak * a[j])) / Den +
2 * P[,i] * ((Num %*% (ak * a[j])) / Den)^2 -
P[,i] * ((Num %*% (ak^2 * a[j]^2)) / Den)
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'rsm', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
setClass("nominal", contains = 'AllItemsClass',
representation = representation(mat='logical'))
setMethod(
f = "print",
signature = signature(x = 'nominal'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'nominal'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'nominal'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'nominal'),
definition = function(x){
d <- x@par[(length(x@par) - x@ncat + 1L):length(x@par)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'nominal'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=-.2, sdlog=.5), 0,
abs(rnorm(x@ncat-1L, (x@ncat-1L) / 2, sd = 1)), 0,
rnorm(x@ncat-1L))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'nominal'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x@est[(x@nfact+1L):(x@nfact + x@ncat)] <- FALSE
x
}
)
P.nominal <- function(par, ncat, Theta, returnNum = FALSE, ot = 0){
return(.Call("nominalTraceLinePts", par, ncat, Theta, returnNum, ot))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'nominal', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.nominal(par=x@par, ncat=x@ncat, Theta=Theta, ot=ot))
}
)
nominalParDeriv <- function(a, ak, d, Theta, P, num, dat, estHess, gpcm = FALSE){
nfact <- length(a)
ncat <- length(d)
akind <- nfact
dind <- nfact + ncat
ak2 <- ak^2
P2 <- P^2
P3 <- P^3
aTheta <- as.vector(Theta %*% a)
aTheta2 <- aTheta^2
dat_num <- dat/num
numsum <- rowSums(num)
numakD <- num %*% ak
numak2D2 <- num %*% ak2
numakDTheta_numsum <- matrix(0, nrow(num), nfact)
for(i in seq_len(nfact))
numakDTheta_numsum[,i] <- (num %*% ak * Theta[, i])/ numsum
ret <- .Call('dparsNominal', a, ak, d, Theta, P, num, dat, nfact, ncat,
akind, dind, ak2, P2, P3, aTheta, aTheta2, dat_num, numsum, numakD,
numak2D2, numakDTheta_numsum, estHess)
ret
}
setMethod(
f = "Deriv",
signature = signature(x = 'nominal', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call("dparsNominal", x, Theta, offterm, FALSE, estHess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'nominal', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
d <- ExtractZetas(x)
ak <- x@par[(x@nfact+1):(x@nfact+x@ncat)]
P <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta)
Num <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta, returnNum = TRUE)
Den <- rowSums(Num)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(x@ncat)){
grad[[i]][ ,j] <- ak[i] * a[j] * P[ ,i] - P[ ,i] *
(Num %*% (ak * a[j])) / Den
hess[[i]][ ,j] <- ak[i]^2 * a[j]^2 * P[ ,i] -
2 * ak[i] * a[j] * P[,i] * (Num %*% (ak * a[j])) / Den +
2 * P[,i] * ((Num %*% (ak * a[j])) / Den)^2 -
P[,i] * ((Num %*% (ak^2 * a[j]^2)) / Den)
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'nominal', Theta = 'matrix'),
definition = function(x, Theta){
num <- P.nominal(x@par, ncat=x@ncat, Theta=Theta, returnNum=TRUE)
den <- rowSums(num)
P <- num/den
a <- x@par[seq_len(ncol(Theta))]
ak <- x@par[(ncol(Theta)+1L):(ncol(Theta)+x@ncat)]
dp <- matrix(0, nrow(Theta), length(x@par))
aknum <- t(ak * t(num))
aTheta <- as.numeric(a %*% t(Theta))
nfact <- ncol(Theta)
ncat <- x@ncat
e <- 0:(x@ncat-1)
eak <- e*ak
for(i in seq_len(nfact)){
for(j in 2L:ncat)
dp[,i] <- dp[,i] + eak[j]*Theta[,i]*P[,j] -
e[j]*P[,j]*rowSums(Theta[,i] * aknum)
}
for(j in seq_len(x@ncat)){
dp[,nfact + ncat + j] <- e[j] * P[,j] - e[j] * P[,j]^2 -
as.numeric(e[-j] %*% t(P[,-j]*P[,j]))
dp[,nfact + j] <- dp[,nfact + ncat + j] * aTheta
}
return(dp)
}
)
# ----------------------------------------------------------------
setClass("partcomp", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'partcomp'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'partcomp'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'partcomp'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'partcomp'),
definition = function(x){
d <- x@par[(x@nfact+1L):(length(x@par)-2L)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'partcomp'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
rlnorm(x@nfact, meanlog=.2, sdlog=.5),
rnorm(1L, -2, .5),
rnorm(1L, 2, .5))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'partcomp'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.comp <- function(par, Theta, ot = 0)
{
return(.Call('partcompTraceLinePts', par, Theta))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'partcomp', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
return(P.comp(x@par, Theta=Theta))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'partcomp', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
#local derivative from previous version with small mod
#u and g in logit form
dpars.comp <- function(lambda,zeta,g,r,f,Thetas,estHess)
{
nfact <- length(lambda)
pars <- c(zeta,lambda,g)
pgrad <- function(pars, r, thetas){
nfact <- ncol(thetas)
d <- pars[seq_len(nfact)]
a <- pars[(nfact+1L):(length(pars)-1L)]
c <- pars[length(pars)]
P <- P.comp(c(a,d,c,999), thetas)[,2L]
Pstar <- P.comp(c(a,d,-999,999),thetas)[,2L]
Qstar <- 1 - Pstar
Q <- 1 - P
c <- antilogit(c)
g_1g <- c * (1 - c)
const1 <- (r/P - (f-r)/Q)
dd <- da <- rep(0,nfact)
dc <- sum(r/P * (g_1g * (1 - Pstar)) + (f-r)/Q * (g_1g * (Pstar - 1)))
for(i in seq_len(nfact)){
Pk <- P.mirt(c(a[i],d[i],-999,999),matrix(thetas[,i]))[,2L]
Qk <- 1 - Pk
dd[i] <- sum((1-c)*Pstar*Qk*const1)
da[i] <- sum((1-c)*Pstar*Qk*thetas[,i]*const1)
}
return(c(dd,da,dc))
}
phess <- function(pars, r, thetas){
nfact <- ncol(thetas)
d <- pars[seq_len(nfact)]
a <- pars[(nfact+1L):(length(pars)-1L)]
c <- pars[length(pars)]
P <- P.comp(c(a,d,c,999), thetas)[,2L]
Pstar <- P.comp(c(a,d,-999,999),thetas)[,2L]
Qstar <- 1 - Pstar
Q <- 1 - P
g <- c <- antilogit(c)
g_1g <- c * (1 - c)
const1 <- (r/P - (f-r)/Q)
const2 <- (r/P^2 + (f-r)/Q^2)
hess <- matrix(0,nfact*2+1,nfact*2+1)
dNames <- paste("d",seq_len(nfact),sep='_')
aNames <- paste("a",seq_len(nfact),sep='_')
Names <- c(paste("d",seq_len(nfact),sep='_'),paste("a",1:nfact,sep='_'),'c_0')
for(i in seq_len(nfact*2+1L)){
for(j in seq_len(nfact*2+1L)){
if(i <= j){
d1 <- strsplit(Names[c(i,j)],"_")[[1L]]
d2 <- strsplit(Names[c(i,j)],"_")[[2L]]
k <- as.numeric(d1[2L])
m <- as.numeric(d2[2L])
Pk <- P.mirt(c(a[k],d[k],-999,999),matrix(thetas[,k]))[,2L]
Qk <- 1 - Pk
Pm <- P.mirt(c(a[m],d[m],-999,999),matrix(thetas[,m]))[,2L]
Qm <- 1 - Pm
if(i == j && d1[1L] == 'd'){
hess[i,i] <- sum(r/P * ((1-g)*(Pstar - 3*Pstar*Pk + 2*Pstar*Pk^2) ) -
r/P^2 * ((1-g) * (Pstar - Pstar*Pk))^2 +
(f-r)/Q * ((1-g)*(-Pstar + 3*Pstar*Pk - 2*Pstar*Pk^2)) -
(f-r)/Q^2 * ((1-g)*(-Pstar + Pstar*Pk))^2)
next
}
if(i == j && d1[1L] == 'a'){
hess[i,i] <- sum(r/P * ((1-g)*thetas[,k]^2*(Pstar - 3*Pstar*Pk + 2*Pstar*Pk^2) ) -
r/P^2 * ((1-g)*thetas[,k] * (Pstar - Pstar*Pk))^2 +
(f-r)/Q * ((1-g)*thetas[,k]^2*(-Pstar + 3*Pstar*Pk - 2*Pstar*Pk^2)) -
(f-r)/Q^2 * ((1-g)*thetas[,k] * (-Pstar + Pstar*Pk))^2)
next
}
if(i == j && d1[1L] == 'c'){
hess[i,i] <- sum(r/P * (g_1g * (2.0*(1-g) - 1.0 - 2.0*(1-g)*Pstar + Pstar)) -
r/P^2 * (g_1g * (1.0 - Pstar)) * (g_1g * (1.0 - Pstar)) +
(f-r)/Q * (g_1g * (-2.0*(1-g) + 1.0 + 2.0*(1-g)*Pstar - Pstar)) -
(f-r)/Q^2 * (g_1g * (-1.0 + Pstar)) * (g_1g * (-1.0 + Pstar)))
next
}
if(d1[1L] == 'a' && d2[1L] == 'a'){
hess[i,j] <- hess[j,i] <- sum((1-c)*thetas[,k]*thetas[,m]*
Qk*Pstar*Qm*(const1 - Pstar*(1-c)*const2))
next
}
if(d1[1L] == 'd' && d2[1L] == 'd'){
hess[i,j] <- hess[j,i] <- sum((1-c)*Qk*Pstar*Qm*(const1 - Pstar*(1-c)*const2))
next
}
if(d1[1L] == 'a' && d2[1L] == 'c'){
hess[i,j] <- hess[j,i] <- sum(r/P*(g_1g*thetas[,k]*(-Pstar + Pstar*Pk)) -
r/P^2*((1-g)*thetas[,k]*(Pstar - Pstar*Pk)*(g_1g*(1 - Pstar))) +
(f-r)/Q * (g_1g*thetas[,k]*(Pstar - Pstar*Pk))-
(f-r)/Q^2 * ((1-g)*thetas[,k]*(-Pstar + Pstar*Pk)*(g_1g*(-1 + Pstar))))
next
}
if(d1[1L] == 'd' && d2[1L] == 'c'){
hess[i,j] <- hess[j,i] <- sum(r/P*(g_1g*(-Pstar + Pstar*Pk)) -
r/P^2*((1-g)*(Pstar - Pstar*Pk)*(g_1g*(1 - Pstar))) +
(f-r)/Q * (g_1g*(Pstar - Pstar*Pk))-
(f-r)/Q^2 * ((1-g)*(-Pstar + Pstar*Pk)*(g_1g*(-1 + Pstar))))
next
}
if(d1[1L] == 'd' && d2[1L] == 'a' && d1[2] == d2[2]){
hess[i,j] <- hess[j,i] <- sum(
r/P * ((1-g)*thetas[,k]*(Pstar - 3*Pstar*Pk + 2*Pstar*Pk^2) ) -
r/P^2 * ((1-g)*thetas[,k] * (Pstar - Pstar*Pk) * (1-g)*(Pstar - Pstar*Pk)) +
(f-r)/Q * ((1-g)*thetas[,k]*(-Pstar + 3*Pstar*Pk - 2*Pstar*Pk^2)) -
(f-r)/Q^2 * ((1-g)*thetas[,k] * (-Pstar + Pstar*Pk) * (1-g)*(-Pstar + Pstar*Pk)))
next
}
if(d1[1L] == 'd' && d2[1L] == 'a' && d1[2] != d2[2]){
hess[i,j] <- hess[j,i] <- sum(r/P * ((1-g)*thetas[,m]*(Pstar - Pstar*Pk - Pstar*Pm + Pstar^2)) -
r/P^2 * ((1-g)*thetas[,m] * (Pstar - Pstar*Pm) * (1-g)*(Pstar - Pstar*Pk)) +
(f-r)/Q * ((1-g)*thetas[,m]*(-Pstar + Pstar*Pk + Pstar*Pm - Pstar^2)) -
(f-r)/Q^2 * ((1-g)*thetas[,m] * (-Pstar + Pstar*Pm) * (1-g)*(-Pstar + Pstar*Pk)))
next
}
}
}
}
return(hess)
}
#old pars in the form d, a, g
g <- pgrad(pars, r, Thetas)
if(estHess) h <- phess(pars, r, Thetas)
else h <- matrix(0, length(g), length(g))
#translate into current version
grad <- c(g[(nfact+1L):(nfact*2)], g[1L:nfact], g[length(g)], 0)
hess <- matrix(0, ncol(h) + 1L, ncol(h) + 1L)
if(estHess){
hess[seq_len(nfact), seq_len(nfact)] <- h[(nfact+1L):(nfact*2),(nfact+1L):(nfact*2)] #a block
hess[(nfact+1L):(nfact*2),(nfact+1L):(nfact*2)] <- h[seq_len(nfact), seq_len(nfact)] #d block
hess[nfact*2 + 1L, nfact*2 + 1L] <- h[nfact*2 + 1L, nfact*2 + 1L] #g
hess[nfact*2 + 1L, 1L:nfact] <- hess[1:nfact, nfact*2 + 1L] <-
h[nfact*2 + 1L, (nfact+1L):(nfact*2)] #ga
hess[nfact*2 + 1L, (nfact+1L):(nfact*2)] <- hess[(nfact+1L):(nfact*2), nfact*2 + 1L] <-
h[nfact*2 + 1L, seq_len(nfact)] #gd
hess[(nfact+1L):(nfact*2), seq_len(nfact)] <- t(h[(nfact+1L):(nfact*2), seq_len(nfact)])
hess[seq_len(nfact), (nfact+1L):(nfact*2)] <- t(h[seq_len(nfact), (nfact+1L):(nfact*2)]) #ads
}
return(list(grad=grad, hess=hess))
}
#####
f <- rowSums(x@dat)
r <- x@dat[ ,2L]
nfact <- x@nfact
a <- x@par[seq_len(nfact)]
d <- x@par[(nfact+1L):(nfact*2L)]
g <- x@par[length(x@par)-1L]
tmp <- dpars.comp(lambda=ExtractLambdas(x),zeta=ExtractZetas(x),g=x@par[nfact*2L + 1L],r=r, f=f,
Thetas=Theta, estHess=estHess)
ret <- list(grad=tmp$grad, hess=tmp$hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'partcomp', Theta = 'matrix'),
definition = function(x, Theta){
N <- nrow(Theta)
nfact <- ncol(Theta)
parlength <- length(x@par)
u <- x@par[parlength]
g <- x@par[parlength - 1L]
d <- ExtractZetas(x)
a <- ExtractLambdas(x)
Pstar <- P.comp(c(a, d, -999, 999), Theta)[,2L]
g <- antilogit(g)
u <- antilogit(u)
grad <- hess <- vector('list', 2L)
grad[[1L]] <- grad[[2L]] <- hess[[1L]] <- hess[[2L]] <- matrix(0, N, nfact)
for(j in seq_len(nfact)){
Pn <- P.mirt(c(a[j], d[j], -999, 999), Theta)[,2L]
grad[[2L]][ ,j] <- (u-g)*Pstar*a[j]*(1 - Pn)
grad[[1L]][ ,j] <- -1*grad[[2L]][ ,j]
hess[[2L]][ ,j] <- (u-g)*a[j]^2*Pstar*(1 - 3*Pn + 2*Pn^2)
hess[[1L]][ ,j] <- -1*hess[[2L]][ ,j]
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'partcomp', Theta = 'matrix'),
definition = function(x, Theta){
# message('partcomp derivatives not optimized') ##TODO
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
setClass("nestlogit", contains = 'AllItemsClass',
representation = representation(correctcat='integer'))
setMethod(
f = "print",
signature = signature(x = 'nestlogit'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'nestlogit'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'nestlogit'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'nestlogit'),
definition = function(x){
stop('not written')
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'nestlogit'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
rnorm(1L),
rnorm(1L, -2, .5),
rnorm(1L, 2, .5), 0,
abs(rnorm(x@ncat-2L, (x@ncat-2L) / 2, sd = 1)), x@ncat,
rnorm(x@ncat-2L))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'nestlogit'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x@est[(x@nfact+5L):(x@nfact + x@ncat + 1L)] <- FALSE
x
}
)
P.nestlogit <- function(par, Theta, correct, ncat)
{
return(.Call("nestlogitTraceLinePts", par, Theta, correct, ncat))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'nestlogit', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
return(P.nestlogit(x@par, Theta=Theta, correct=x@correctcat, ncat=x@ncat))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'nestlogit', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
#u and g in logit
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
dat <- x@dat
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta,
gradient = FALSE)
nfact <- x@nfact
a <- x@par[seq_len(x@nfact)]
d <- x@par[x@nfact+1L]
g <- x@par[x@nfact+2L]
u <- x@par[x@nfact+3L]
ak <- x@par[(x@nfact+4L):(x@nfact+4L+x@ncat-2L)]
dk <- x@par[(length(x@par)-length(ak)+1):length(x@par)]
correct <- x@correctcat
Pd <- P.mirt(c(a, d, g, u), Theta=Theta)[,2L]
g <- antilogit(g)
u <- antilogit(u)
Qd <- 1 - Pd
Pstar <- P.mirt(c(a, d, -999, 999), Theta=Theta)[,2L]
Qstar <- 1 - Pstar
num <- P.nominal(c(rep(1, nfact), ak, dk), ncat=length(dk),
Theta=Theta, returnNum=TRUE)
den <- rowSums(num)
Pn <- num/den
cdat <- dat[,correct]
idat <- dat[,-correct]
nd <- ncol(idat)
g_1g <- g * (1 - g)
u_1u <- u * (1 - u)
for(i in seq_len(nfact))
grad[i] <- sum( (u-g) * Theta[,i] * Qstar * Pstar * (
cdat / Pd - rowSums(idat/Qd)) )
grad[nfact+1L] <- sum( (u-g) * Qstar * Pstar * (
cdat / Pd - rowSums(idat/Qd)) )
grad[nfact+2L] <- sum( ((cdat * g_1g * (1-Pstar)/Pd) + rowSums(idat * g_1g * (Pstar - 1)/Qd)) )
grad[nfact+3L] <- sum( (cdat * u_1u * Pstar / Pd - rowSums(idat * u_1u * Pstar / Qd) ))
for(j in seq_len(nd)){
grad[nfact+3L+j] <- sum((
(idat[,j] * Qd * rowSums(Theta) * (Pn[,j] - Pn[,j]^2)) / (Qd * Pn[,j]) -
rowSums(idat[,-j, drop=FALSE]) * rowSums(Theta) * Pn[,j]))
grad[nfact+3L+nd+j] <- sum((
(idat[,j] * Qd * (Pn[,j] - Pn[,j]^2)) / (Qd * Pn[,j]) -
rowSums(idat[,-j, drop=FALSE]) * Pn[,j]))
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'nestlogit', Theta = 'matrix'),
definition = function(x, Theta){
a <- x@par[seq_len(x@nfact)]
d <- x@par[x@nfact+1L]
g <- x@par[x@nfact+2L]
u <- x@par[x@nfact+3L]
ak <- x@par[(x@nfact+4L):(x@nfact+4L+x@ncat-2L)]
dk <- x@par[(length(x@par)-length(ak)+1):length(x@par)]
Pn <- P.nominal(c(rep(1,ncol(Theta)), ak, dk), ncat=length(dk), Theta=Theta)
Num <- P.nominal(c(rep(1,ncol(Theta)), ak, dk), ncat=length(dk),
Theta=Theta, returnNum = TRUE)
Den <- rowSums(Num)
Pstar <- P.mirt(c(a, d, -999, 999), Theta)[,2L]
Q <- P.mirt(c(a, d, g, u), Theta)[,1]
g <- antilogit(g)
u <- antilogit(u)
Num2 <- P <- matrix(0, nrow(Theta), x@ncat)
P[,-x@correctcat] <- Pn
Num2[,-x@correctcat] <- Num
Num <- Num2
ak2 <- dk2 <- numeric(x@ncat)
ak2[-x@correctcat] <- ak
dk2[-x@correctcat] <- dk
ak <- ak2
dk <- dk2
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(x@ncat)){
if(i == x@correctcat){
grad[[i]][ ,j] <- (u-g) * a[j] * (Pstar * (1 - Pstar))
hess[[i]][ ,j] <- 2 * (u - g) * a[j]^2 * ((1 - Pstar)^2 * Pstar) -
(u - g) * a[j]^2 * (Pstar * (1 - Pstar))
} else {
grad[[i]][ ,j] <- -(u-g) * a[j] * (Pstar * (1 - Pstar)) * P[,i] +
Q * (ak[i] * P[ ,i] - P[ ,i] * (Num %*% (ak)) / Den)
hess[[i]][ ,j] <-
-2 * (u - g) * a[j]^2 * (1 - Pstar)^2 * Pstar * P[,i] +
(u - g) * a[j]^2 * Pstar * (1 - Pstar) * P[,i] -
2 * (u - g) * a[j] * ak[i] * (1 - Pstar) * Pstar * P[,i] +
2 * a[j] * (Pstar * (1 - Pstar)) * P[,i] * (Num %*% (ak)) / Den +
Q * ak[i]^2 * P[ ,i] -
2 * Q * ak[i] * P[,i] * (Num %*% (ak)) / Den +
2 * Q * P[,i] * ((Num %*% (ak)) / Den)^2 -
Q * P[,i] * ((Num %*% (ak^2)) / Den)
}
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'nestlogit', Theta = 'matrix'),
definition = function(x, Theta){
# message('nestlogit derivatives not optimized') ##TODO
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
setClass("ideal", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'ideal'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'ideal'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'ideal'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'ideal'),
definition = function(x){
d <- x@par[length(x@par)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'ideal'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
-rlnorm(1, meanlog=0, sdlog=.5))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'ideal'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.ideal <- function(par, Theta, ot = 0)
{
alpha <- par[length(par)]
beta <- par[-length(par)]
eta <- -0.5*(Theta %*% beta + alpha)^2
eta <- ifelse(eta < -20, -20, eta)
P <- exp(eta)
P <- cbind((1-P), P)
P <- ifelse(P < 1e-20, 1e-20, P)
P <- ifelse(P > (1 - 1e-20), (1 - 1e-20), P)
return(P)
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'ideal', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
return(P.ideal(x@par, Theta=Theta))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'ideal', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
d <- x@par[length(x@par)]
a <- x@par[-length(x@par)]
P <- ProbTrace(x, Theta=Theta)[,2L]
Q <- (1 - P)
Q <- ifelse(Q < 1e-7, 1e-7, Q)
int <- as.numeric(Theta %*% a + d)
for(i in seq_len(ncol(Theta)))
grad[i] <- -sum( x@dat[,1] * int * Theta[,i] * -P / Q +
x@dat[,2] * int * Theta[,i])
grad[i+1L] <- -sum(2 * x@dat[,1] * int * -P / Q +
2 * x@dat[,2] * int)/2
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
return(list(grad = grad, hess=hess))
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'ideal', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta)
# N <- nrow(Theta)
# nfact <- ncol(Theta)
# P <- ProbTrace(x, Theta=Theta)[,2L]
# d <- x@par[length(x@par)]
# a <- x@par[-length(x@par)]
# int <- as.numeric(t(a %*% t(Theta)) + d)
# grad <- hess <- vector('list', 2L)
# grad[[1L]] <- grad[[2L]] <- hess[[1L]] <- hess[[2L]] <- matrix(0, N, nfact)
# for(i in 1L:nfact){
# grad[[2L]][ ,i] <- 2 * a[i] * int * P
# grad[[1L]][ ,i] <- -1 * grad[[2L]][ ,i]
# hess[[2L]][ ,i] <- 2 * a[i]^2 * int * P + 4 * int^2 * a[i]^2 * P
# hess[[1L]][ ,i] <- -1 * hess[[2L]][ ,i]
# }
# return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'ideal', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta=Theta)[,2L]
dp <- matrix(0, nrow(Theta), length(x@par))
d <- x@par[length(x@par)]
a <- x@par[-length(x@par)]
int <- as.numeric(t(a %*% t(Theta)) + d)
for(i in seq_len(ncol(Theta)))
dp[,i] <- -2 * Theta[,i] * int * P
dp[,i+1L] <- -2 * int * P
dp
}
)
# ----------------------------------------------------------------
setClass("lca", contains = 'AllItemsClass',
representation = representation(score='numeric',
item.Q='matrix'))
setMethod(
f = "print",
signature = signature(x = 'lca'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'lca'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'lca'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'lca'),
definition = function(x){
stop('not written')
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'lca'),
definition = function(x){
par <- rnorm(length(x@par))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'lca'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.lca <- function(par, Theta, ncat, item.Q, returnNum = FALSE){
return(.Call('lcaTraceLinePts', par, Theta, item.Q, ncat, returnNum))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'lca', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
P <- P.lca(x@par, Theta=Theta, ncat=x@ncat, item.Q=x@item.Q)
return(P)
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'lca', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
# ret <- .Call('dparslca', x@par, Theta, x@item.Q,
# estHess, x@dat, offterm) #TODO use estHess in Rcpp
ret <- list(grad=numeric(length(x@par)),
hess=matrix(0, length(x@par), length(x@par)))
ret$grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta)
if(estHess && any(x@est)){
ret$hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
}
# if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'lca', Theta = 'matrix'),
definition = function(x, Theta){
stop('not written')
}
)
setMethod(
f = "dP",
signature = signature(x = 'lca', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta)
dp <- matrix(0, nrow(Theta), length(x@par))
ind <- 1L
for(j in 2L:x@ncat){
for(i in seq_len(ncol(Theta))){
dp[,ind] <- Theta[,i] * (P[,j] -
rowSums(P[,j,drop=FALSE] * P[,j,drop=FALSE]))
ind <- ind + 1L
}
}
dp
}
)
# ----------------------------------------------------------------
setClass("spline", contains = 'AllItemsClass',
representation = representation(stype='character',
Theta_prime='matrix',
item.Q='matrix',
sargs='list'))
setMethod(
f = "print",
signature = signature(x = 'spline'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'spline'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'spline'),
definition = function(x){
numeric(x@nfact)
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'spline'),
definition = function(x){
stop('not written')
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'spline'),
definition = function(x){
par <- rnorm(length(x@par), sd = 20)
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'spline'),
definition = function(x){
stop('spline null should not be run')
}
)
setMethod(
f = "ProbTrace",
signature = signature(x = 'spline', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(Theta) != nrow(x@Theta_prime))
x <- loadSplineParsItem(x, Theta)
P <- P.lca(x@par, item.Q=x@item.Q,
Theta=x@Theta_prime, ncat=x@ncat)
return(P)
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'spline', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
ret <- .Call('dparslca', x@par, x@Theta_prime, x@item.Q, FALSE, x@dat, offterm)
if(estHess && any(x@est)){
ret$hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=x@Theta_prime, gradient=FALSE)
}
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'spline', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta)
}
)
setMethod(
f = "dP",
signature = signature(x = 'spline', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta)
dp <- matrix(0, nrow(Theta), length(x@par))
ind <- 1L
for(j in 2L:x@ncat){
for(i in seq_len(ncol(Theta))){
dp[,ind] <- Theta[,i] * (P[,j] -
rowSums(P[,j,drop=FALSE] * P[,j,drop=FALSE]))
ind <- ind + 1L
}
}
dp
}
)
# ----------------------------------------------------------------
setClass("ggum", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'ggum'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'ggum'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'ggum'),
definition = function(x){
x@par[1L:x@nfact] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'ggum'),
definition = function(x){
x@par[(x@nfact+1L):length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'ggum'),
definition = function(x){
par <- c(rlnorm(x@nfact, .2, .2),
rnorm(x@nfact),
sort(rnorm(x@ncat-1L, 0, 1.5), decreasing = TRUE))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'ggum'),
definition = function(x){
x@par[1L:x@nfact] <- 0
x@est[1L:x@nfact] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'ggum', Theta = 'matrix'),
definition = function(x, Theta){
return(P.ggum(x@par, Theta=Theta, ncat=x@ncat))
}
)
# TODO write this with Rcpp
P.ggum <- function(par, Theta, ncat)
{
# C <- ncat - 1
# D <- (length(par)-C)/2
# M <- 2*C + 1
#
# sumtau <- 0
# sumdist<-0
# num <- matrix(nrow=nrow(Theta),ncol=(C+1))
# P <- matrix(nrow=nrow(Theta),ncol=(C+1))
#
# for (d in 1:D) {
# dist <- (par[d]**2)*((Theta[,d] - par[(D+d)])**2)
# sumdist <- dist+sumdist
# }
# dist<-sqrt(sumdist)
# z1 <- z2 <- num
#
# for (w in 0:C) {
#
# x1 <- w*dist
# x2 <- (M-w)*dist
#
# if (w>0) {
# for (d in 1:D) {
# tau <- (par[d]*par[(w+2*D)])
# sumtau <- tau + sumtau
# }
# }
#
# z1[,w + 1] <- x1 + sumtau
# z2[,w + 1] <- x2 + sumtau
#
# # num1 <- exp(x1 + sumtau)
# # num2 <- exp(x2 + sumtau)
# #
# # num[,(w+1)] <- num1 + num2
#
# }
#
# # numtot <- apply(num, 1, sum)
# #
# # for (k in 1:(C+1)) {
# # P[,k] <- num[,k]/numtot
# # }
#
# # less overflow this way
# maxz1 <- apply(z1, 1, max)
# maxz2 <- apply(z2, 1, max)
# maxz <- pmax(maxz1, maxz2)
# num2 <- exp(z1 - maxz) + exp(z2 - maxz)
# den <- rowSums(num2)
# P <- num2/den
# P <- ifelse(P < 1e-20, 1e-20, P)
# P <- ifelse(P > (1 - 1e-20), (1 - 1e-20), P)
# return(P)
return(.Call("ggumTraceLinePts", par, Theta, ncat))
}
# complete-data derivative used in parameter estimation (here it is done numerically)
setMethod(
f = "Deriv",
signature = signature(x = 'ggum', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
# NOT USED, but useful for numerical checking
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta,
type='Richardson')
if(estHess){
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta,
gradient=FALSE, type='Richardson')
}
return(list(grad=grad, hess=hess))
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'ggum', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'ggum', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass('custom', contains = 'AllItemsClass',
representation = representation(name='character',
P='function',
gr='function',
usegr='logical',
hss='function',
gen='function',
usehss='logical',
userdata='list',
useuserdata='logical',
derivType='character',
derivType.hss='character'))
setMethod(
f = "print",
signature = signature(x = 'custom'),
definition = function(x, ...){
cat('Custom item object named:', x@name)
}
)
setMethod(
f = "show",
signature = signature(object = 'custom'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'custom'),
definition = function(x){
a <- rep(.001, x@nfact)
a
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'custom'),
definition = function(x){
stop('not written')
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'custom'),
definition = function(x){
x@par <- x@gen(x)
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'custom'),
definition = function(x){
stop('calculating null models for custom item types is ambiguous.',
call. = FALSE)
}
)
setMethod(
f = "ProbTrace",
signature = signature(x = 'custom', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(x@useuserdata) return(x@P(x@par, Theta, x@ncat, x@userdata[[1L]]))
else return(x@P(x@par, Theta, x@ncat))
}
)
setMethod("initialize",
'custom',
function(.Object, name, par, est, parnames, lbound, ubound,
P, gr, hss, gen, userdata, derivType, derivType.hss, dps, dps2) {
dummyfun <- function(...) return(NULL)
names(est) <- names(par)
usegr <- usehss <- useuserdata <- TRUE
.Object@name <- name
.Object@par <- par
.Object@est <- est
.Object@parnames <- parnames
.Object@P <- P
.Object@derivType <- derivType
.Object@derivType.hss <- derivType.hss
.Object@itemclass <- 999L
.Object@dps <- dps
.Object@dps2 <- dps2
if(is.null(gr)){
.Object@gr <- dummyfun
usegr <- FALSE
} else .Object@gr <- gr
if(is.null(hss)){
.Object@hss <- dummyfun
usehss <- FALSE
} else .Object@hss <- hss
if(is.null(gen)){
.Object@gen <- function(object) object@par
} else .Object@gen <- gen
if(is.null(userdata)){
.Object@userdata <- list()
useuserdata <- FALSE
} else .Object@userdata <- userdata
.Object@usegr <- usegr
.Object@usehss <- usehss
.Object@useuserdata <- useuserdata
.Object@lbound <- if(!is.null(lbound)) lbound else rep(-Inf, length(par))
.Object@ubound <- if(!is.null(ubound)) ubound else rep(Inf, length(par))
.Object
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'custom', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(x@usegr) grad <- x@gr(x, Theta)
else grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, type=x@derivType)
if(estHess){
if(x@usehss) hess <- x@hss(x, Theta)
else hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, type=x@derivType.hss,
gradient=FALSE)
}
return(list(grad = grad, hess=hess))
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'custom', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta)
}
)
setMethod(
f = "dP",
signature = signature(x = 'custom', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
# itemtype='grsmIRT', 1-dim graded rating scale
# slope, thresholds and 1 location parameters
# ExtractZetas, GenRandomPar, ProbTrace, initialize MODIFIED
setClass("grsmIRT", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'grsmIRT'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'grsmIRT'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'grsmIRT'),
definition = function(x){
x@par[seq_len(x@nfact)] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'grsmIRT'),
definition = function(x){
#x@par[(x@nfact+1):(length(x@par))] #intercepts
x@par[(x@nfact+1):(x@ncat+1)] #intercepts
# if we set :(length(x@par)), when collapsed, may produce error?
# MAYBE DOESNT MATTER?
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'grsmIRT'),
definition = function(x){
par <- c(rlnorm(1,0,1),sort(rnorm(x@ncat-1, 0, 1), decreasing=TRUE), rnorm(1,0,0.5))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'grsmIRT'),
definition = function(x){
x@par[seq_len(x@nfact+x@ncat)] <- 0
x@est[seq_len(x@nfact+x@ncat)] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'grsmIRT', Theta = 'matrix'),
definition = function(x, Theta){
th1 = Theta[,1];
ncat = x@ncat
a1 = x@par[1]
d = x@par[2:ncat]
c = x@par[ncat+1]
nr = nrow(Theta); nc=ncat-1
if (all(d == sort(d, decreasing=TRUE))) {
D.star = matrix(c+d, nrow=nr, ncol=nc, byrow=TRUE)
TH1 = matrix(th1, nrow=nr, ncol=nc)
A = matrix(a1, nrow=nr, ncol=nc)
P.star = 1/(1+exp(-A*(D.star+TH1)))
P=cbind(1, P.star)-cbind(P.star, 0)
P <- ifelse(P < 1e-20, 1e-20, P)
P <- ifelse(P > (1 - 1e-20), (1 - 1e-20), P)
} else {
P <- matrix(1e-20, nrow=nr, ncol=ncat)}
return(P)
}
)
# complete-data derivative used in parameter estimation
# Analytic Derivation is provided...
# and it fails, numerical derivation is used.
setMethod(
f = "Deriv",
signature = signature(x = 'grsmIRT', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
dat <- x@dat;
th1 <- Theta[,1];
ncat <- x@ncat
a1 <- x@par[1]
d <- x@par[2:ncat]
c <- x@par[ncat+1]
nr <- length(th1);
nc <- ncat - 1
beta.mean <- mean(d)
beta.mat <- matrix(d, nrow=nr, ncol=nc, byrow=TRUE)
Sigma.mat <- (beta.mat - beta.mean)
th1.mat <- matrix(th1, nr, nc)
B.mat <- th1.mat + c + beta.mat
lam.star.mat <- exp(-a1*B.mat)
lam.mat <- cbind(lam.star.mat, 1)-cbind(0, lam.star.mat)
P.star <- 1/(1+lam.star.mat)
Q.star <- 1- P.star
P <- P.star[,-nc] - P.star[,-1]
inv.P <- 1/P
#when a1*B.mat >30, inv.P would be overflow...
PQ.star <- P.star*Q.star
P.star.a1 <- PQ.star*B.mat
P.star.beta <- list()
for (i in seq_len(ncat-1))
P.star.beta[[i]] <- matrix(0, nr, nc)
for (i in seq_len(ncat-1))
P.star.beta[[i]][,i] <- (PQ.star*a1)[,i]
P.star.c <- PQ.star*a1
P.a1 <- matrix(NA, nr, ncat);
P.c <- matrix(NA, nr, ncat);
P.beta <- list(); for (i in 1:(ncat-1)) P.beta[[i]] <- matrix(NA, nr, ncat);
P.a1[, 2:nc] <- inv.P*(P.star.a1[,-nc] - P.star.a1[, -1])
for (i in seq_len(ncat-1))
P.beta[[i]][, 2:nc] <- inv.P*(P.star.beta[[i]][, -nc] - P.star.beta[[i]][, -1])
P.c[, 2:nc] <- inv.P * (P.star.c[, -nc] - P.star.c[, -1])
# the first response category
P.a1[,1] <- -P.star[,1]*B.mat[,1]
for (i in 2:(ncat-1))
P.beta[[i]][,1] <- 0
P.beta[[1]][,1] <- -P.star[,1]*a1
P.c[,1] <- -P.star[,1]*a1
# the last response category
P.a1[,ncat]=Q.star[,ncat-1]*B.mat[,ncat-1]
for (i in seq_len(ncat-2))
P.beta[[i]][,ncat] <- 0
P.beta[[ncat-1]][,ncat] <- Q.star[,ncat-1]*a1
P.c[,ncat] <- Q.star[,ncat-1]*a1
grad[1] <- sum(P.a1 * dat)
for (i in seq_len(ncat-1))
grad[i+1] <- sum(P.beta[[i]] *dat)
grad[ncat+1] <- sum(P.c * dat)
grad[!x@est] <- 0
Nest <- is.nan(grad) | is.infinite(grad)
if (sum(Nest)>0) {
grad <- rep(0, length(x@par))
grad[x@est & Nest] <- numerical_deriv(x@par[x@est & Nest], EML, obj=x,
Theta=Theta)
}
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'grsmIRT', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'grsmIRT', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# defines how the item should be initiallized with the S4 function new()
# before the parameter estimates are updated
# (set starting values, lower/upper bounds, logicals indicating estimation, etc)
setMethod("initialize",
'grsmIRT',
function(.Object, nfact, ncat){
if(nfact != 1L)
stop('grsmIRT only possible for unidimensional models')
stopifnot(ncat >= 2L)
.Object@par <- c(rep(1, nfact), seq(1, -1, length.out=ncat-1), 0)
#.Object@par <- c(rep(1, nfact), seq(-3, 3, length.out=ncat-1), 0)
# -3 ~ 3 seems to be too far away
names(.Object@par) = c(paste("a", seq_len(nfact), sep=""),
paste("b", seq_len(ncat-1L), sep=""), "c")
.Object@est <- c(rep(TRUE, nfact), rep(TRUE,ncat-1L), TRUE)
.Object@parnames <- names(.Object@par)
.Object@lbound <- rep(-Inf, nfact+ncat)
.Object@ubound <- rep(Inf, nfact+ncat)
.Object
}
)
# ----------------------------------------------------------------
setClass("gpcmIRT", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'gpcmIRT'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'gpcmIRT'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'gpcmIRT'),
definition = function(x){
x@par[1L]
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'gpcmIRT'),
definition = function(x){
x@par[-c(1, length(x@par))]
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'gpcmIRT'),
definition = function(x){
par <- c(rlnorm(1,0,1),
sort(rnorm(x@ncat-1, 0, 1), decreasing=FALSE), 0)
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'gpcmIRT'),
definition = function(x){
x@par[1L] <- 0
x@est[1L] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'gpcmIRT', Theta = 'matrix'),
definition = function(x, Theta){
.Call('gpcmIRTTraceLinePts', x@par, Theta, FALSE, numeric(0))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'gpcmIRT', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call("dparsgpcmIRT", x@par, Theta, offterm, x@dat,
length(x@par) - ncol(Theta), FALSE)
hess <- matrix(0, length(x@par), length(x@par))
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
ret <- list(grad=ret$grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
# derivative of the model wrt to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'gpcmIRT', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta)
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'gpcmIRT', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("monopoly", contains = 'AllItemsClass',
representation = representation(k='integer'))
setMethod(
f = "print",
signature = signature(x = 'monopoly'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'monopoly'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'monopoly'),
definition = function(x){
x@par[1L]
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'monopoly'),
definition = function(x){
x@par[-c(1, length(x@par))]
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'monopoly'),
definition = function(x){
par <- c(rnorm(1),
sort(rnorm(x@ncat-1, 0, 2)),
rnorm(x@k*2))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'monopoly'),
definition = function(x){
x@par[1L] <- 0
x@est[1L] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'monopoly', Theta = 'matrix'),
definition = function(x, Theta){
.Call('monopolyTraceLinePts', x@par, Theta, x@ncat, x@k)
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'monopoly', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
grad <- numeric(length(x@par))
grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta)
hess <- matrix(0, length(x@par), length(x@par))
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'monopoly', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'monopoly', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("sequential", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'sequential'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'sequential'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'sequential'),
definition = function(x){
x@par[seq_len(x@nfact)] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'sequential'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'sequential'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'sequential'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'sequential', Theta = 'matrix'),
definition = function(x, Theta, itemexp = FALSE){
par <- x@par
a <- par[1L:ncol(Theta)]
d <- par[-c(1L:ncol(Theta))]
ncat <- length(d) + 1L
Z <- matrix(0, nrow(Theta), ncat - 1L)
for(i in seq_along(d))
Z[,i] <- Theta %*% a + d[i]
P <- matrix(0, nrow(Theta), ncat)
Psi <- matrix(0, nrow(Theta), ncat + 1L)
Psi[,1L] <- 1
for(i in seq_along(d))
Psi[,i+1L] <- 1 / (1 + exp(-Z[,i]))
Pkc <- t(apply(Psi, 1L, cumprod))
for(i in 1L:ncat)
P[,i] <- Pkc[,i] * (1 - Psi[,i+1L])
if(itemexp) return(list(Psi=Psi, P=P))
P
}
)
# complete-data derivative used in parameter estimation (here it is done numerically)
setMethod(
f = "Deriv",
signature = signature(x = 'sequential', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
dplist <- ProbTrace(x, Theta, itemexp=TRUE)
P <- dplist$P
Psi <- dplist$Psi
Psi2 <- Psi^2
ncat <- x@ncat
nfact <- ncol(Theta)
Pgrad <- vector("list", length(x@par))
storage <- matrix(0, nrow(Theta), ncat)
# as
ind <- 1L
for(fact in seq_len(nfact)){
for(cat in seq_len(ncat)){
collecta <- (cat - 1) * Theta[,fact] * P[,cat]
tmp <- Psi[,1L:cat, drop=FALSE]
prodtmp <- apply(tmp, 1L, prod)
if(cat > 1L){
for(tcat in seq_len(cat-1L)){
collecta <- collecta -
Theta[,fact] * (Psi[,tcat+1L] * P[,cat])
}
}
collecta <- collecta + prodtmp *
(Theta[,fact] * (Psi2[,cat + 1L] - Psi[,cat + 1L]))
storage[,cat] <- collecta
}
Pgrad[[ind]] <- storage
ind <- ind + 1L
}
#ds
for(cat in seq_len(ncat-1L)){
tmp <- Psi[,1L:cat, drop=FALSE]
prodtmp <- apply(tmp, 1L, prod)
for(tcat in 1L:ncat){
if(tcat < cat){
storage[,tcat] <- 0
} else if(tcat == cat){
storage[,tcat] <- prodtmp * (Psi2[,cat + 1L] - Psi[,cat + 1L])
} else {
storage[,tcat] <- P[,tcat] - Psi[,cat + 1L] * P[,tcat]
}
}
Pgrad[[ind]] <- storage
ind <- ind + 1L
}
Pgrad <- array(do.call(c, Pgrad), c(nrow(Theta), x@ncat, length(x@par)))
grad <- symbolicGrad_par(x, Theta, dp1=Pgrad, P=P)
hess <- matrix(0, length(x@par), length(x@par))
if(estHess && any(x@est)){
# grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta)
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
}
return(list(grad=grad, hess=hess)) #TODO replace with analytic derivatives
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'sequential', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'sequential', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("fivePL", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'fivePL'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'fivePL'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'fivePL'),
definition = function(x){
x@par[seq_len(x@nfact)] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'fivePL'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'fivePL'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
rnorm(1L),
rnorm(1L, -2, .5),
rnorm(1L, 2, .5),
1)
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'fivePL'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.5PL <- function(x, Theta, ncat = 2L){
a <- x@par[1L]
d <- x@par[2L]
g <- 1 / (1 + exp(-x@par[1L + 2L]))
u <- 1 / (1 + exp(-x@par[1L + 3L]))
S <- exp(x@par[5L])
P <- g + (u-g)/ (1 + exp(-1*(a * Theta + d)) )^S
cbind(1-P, P)
}
dP.5PL <- function(x, Theta, ncat = 2L)
{
.e1 <- x@par
.e3 <- exp(-(Theta * .e1[1L] + .e1[2L]))
.e4 <- exp(.e1[5L])
.e5 <- 1 + .e3
.e8 <- .e5^.e4
.e13 <- 1/(1 + exp(-.e1[4L])) - 1/(1 + exp(-.e1[3L]))
.e14 <- .e5^(1 + .e4)
.e15 <- 1/.e8
.e18 <- .e13 * .e3 * .e4/.e14
.e20 <- .e13 * log1p(.e3)/.e8
.e21 <- 1 - .e15
.e25 <- Theta * .e13 * .e3 * .e4/.e14
c(a = cbind(-.e25, .e25),
d = cbind(-.e18, .e18), g = cbind(-.e21, .e21),
u = cbind(-.e15, .e15), S = cbind(.e20, -.e20))
}
d2P.5PL <- function(x, Theta, ncat = 2L)
{
.e1 <- x@par
.e3 <- exp(-(Theta * .e1[1L] + .e1[2L]))
.e4 <- exp(.e1[5L])
.e5 <- 1 + .e4
.e6 <- 1 + .e3
.e7 <- .e6^.e5
.e14 <- 1/(1 + exp(-.e1[4L])) - 1/(1 + exp(-.e1[3L]))
.e15 <- log1p(.e3)
.e16 <- 1/.e7
.e17 <- 2 * .e5
.e18 <- .e6^.e4
.e22 <- .e6^(.e4 - .e17) * .e5 * .e3 - .e16
.e24 <- .e3 * .e4/.e7
.e25 <- .e15/.e18
.e28 <- Theta * .e3 * .e4/.e7
.e29 <- -.e24
.e30 <- -.e25
.e31 <- -.e28
.e35 <- .e16 - .e6^(.e5 - .e17) * .e4 * .e15
.e38 <- .e4 * .e15/.e7 - .e16
.e42 <- Theta * .e22 * .e14 * .e3 * .e4
.e46 <- .e22 * .e14 * .e3 * .e4
.e48 <- .e35 * .e14 * .e3
.e50 <- .e14 * .e3 * .e38
.e52 <- .e14 * .e15^2/.e18
.e53 <- cbind(.e29, .e24)
.e54 <- cbind(.e30, .e25)
.e55 <- cbind(-.e42, .e42)
.e56 <- cbind(.e31, .e28)
.e57 <- cbind(.e24, .e29)
.e58 <- cbind(.e25, .e30)
.e59 <- cbind(.e28, .e31)
.e62 <- Theta * .e35 * .e14 * .e3
.e65 <- Theta * .e14 * .e3 * .e38
.e70 <- Theta^2 * .e22 * .e14 * .e3 * .e4
c(a = c(a = cbind(-.e70, .e70), d = .e55, g = .e59, u = .e56,
S = cbind(.e65, -.e65)),
d = c(a = .e55, d = cbind(-.e46, .e46), g = .e57, u = .e53,
S = cbind(.e50, -.e50)),
g = c(a = .e59, d = .e57, g = c(0, 0), u = c(0, 0), S = .e54),
u = c(a = .e56, d = .e53, g = c(0, 0), u = c(0, 0), S = .e58),
S = c(a = cbind(-.e62, .e62), d = cbind(-.e48, .e48), g = .e54, u = .e58,
S = cbind(-.e52, .e52)))
}
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'fivePL', Theta = 'matrix'),
definition = function(x, Theta, itemexp = FALSE){
ret <- P.5PL(x=x, Theta=Theta)
ret
}
)
# complete-data derivative used in parameter estimation
setMethod(
f = "Deriv",
signature = signature(x = 'fivePL', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(any(x@est)){
dp1 <- array(dP.5PL(x, Theta), c(nrow(Theta),x@ncat,length(x@par)))
grad <- symbolicGrad_par(x, Theta, dp1=dp1)
if(estHess){
dp2 <- array(d2P.5PL(x, Theta), c(nrow(Theta),x@ncat,length(x@par), length(x@par)))
hess <- symbolicHessian_par(x, Theta, dp1=dp1, dp2=dp2)
}
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
ret
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'fivePL', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'fivePL', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("cll", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'cll'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'cll'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'cll'),
definition = function(x){
NA
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'cll'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'cll'),
definition = function(x){
par <- rnorm(1L)
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'cll'),
definition = function(x){
x
}
)
P.cll <- function(x, Theta, ncat = 2L){
b <- x@par[1L]
P <- 1 - exp(-1 * exp(Theta - b))
P <- ifelse(P > 1 - 1e-16, 1 - 1e-16, P)
cbind(1-P, P)
}
dP.cll <- function(x, Theta, ncat = 2L)
{
.e1 <- exp(Theta - x@par[1L])
.e4 <- exp(-.e1) * .e1
cbind(.e4, -.e4)
}
d2P.cll <- function(x, Theta, ncat = 2L)
{
.e1 <- exp(Theta - x@par[1L])
.e5 <- exp(-.e1) * .e1 * (.e1 - 1)
cbind(.e5, -.e5)
}
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'cll', Theta = 'matrix'),
definition = function(x, Theta, itemexp = FALSE){
ret <- P.cll(x=x, Theta=Theta)
ret
}
)
# complete-data derivative used in parameter estimation
setMethod(
f = "Deriv",
signature = signature(x = 'cll', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(any(x@est)){
dp1 <- array(dP.cll(x, Theta), c(nrow(Theta),x@ncat,length(x@par)))
grad <- symbolicGrad_par(x, Theta, dp1=dp1)
if(estHess){
dp2 <- array(d2P.cll(x, Theta), c(nrow(Theta),x@ncat,length(x@par), length(x@par)))
hess <- symbolicHessian_par(x, Theta, dp1=dp1, dp2=dp2)
}
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
ret
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'cll', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'cll', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("ull", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'ull'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'ull'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'ull'),
definition = function(x){
NA
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'ull'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'ull'),
definition = function(x){
par <- c(rlnorm(1L, .2, .3), sort(rnorm(x@ncat-1L), decreasing = TRUE))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'ull'),
definition = function(x){
x@par[1L] <- 0
x@est[1L] <- FALSE
x
}
)
P.ull <- function(x, Theta, ncat){
eta <- x@par[1L]
lambdas <- exp(x@par[2L:length(x@par)])
TheEta <- Theta^eta
ret <- matrix(0, nrow(Theta), ncat)
PS <- matrix(0, nrow(Theta), ncat+1L)
PS[,1L] <- 1
for(i in 1L:(ncat-1L))
PS[,i+1L] <- lambdas[i] * TheEta / (1 + lambdas[i] * TheEta)
for(i in 1L:ncat)
ret[,i] <- PS[,i] - PS[,i+1L]
ret
}
dP.ull <- function(x, Theta, ncat)
{
N <- nrow(Theta)
eta <- x@par[1L]
lambdas <- exp(x@par[2L:length(x@par)])
TheEta <- Theta^eta
Theta2Eta <- Theta^(2 * eta)
dlambdas <- vector('list', ncat-1L)
for(j in 2L:ncat - 1L){
mat <- matrix(0, N, ncat)
tmp <- (TheEta - Theta2Eta * lambdas[j] /
(1 + TheEta * lambdas[j])) / (1 + TheEta * lambdas[j])
mat[,j] <- -tmp
mat[,j+1L] <- tmp
dlambdas[[j]] <- mat
}
dlambdas <- do.call(cbind, dlambdas)
.e9 <- log(Theta)
.e10_e14 <- 1 + as.vector(TheEta) * matrix(lambdas, N, ncat-1L, byrow=TRUE)
mat <- matrix(0, N, ncat)
mat[,1L] <- -(.e9 * lambdas[1L] *
(TheEta - Theta2Eta * lambdas[1L] / .e10_e14[,1L])/.e10_e14[,1L])
.e25 <- TheEta - Theta2Eta * lambdas[1L] / .e10_e14[,1L]
.e26 <- TheEta - Theta2Eta * lambdas[ncat-1L] / .e10_e14[,ncat-1L]
mat[,ncat] <- lambdas[ncat-1L] * .e9 * .e26 / .e10_e14[,ncat-1L]
if(ncat > 2L){
for(j in 2:(ncat - 2L)){
tmp <- lambdas[j] * (TheEta - Theta2Eta * lambdas[j] /
.e10_e14[,j])/.e10_e14[,j]
if(j == 2L)
mat[,2L] <- (lambdas[1L] * .e25 / .e10_e14[,1L] - tmp) * .e9
if(ncat == 3L) break
tmp2 <- lambdas[j+1L] * (TheEta - Theta2Eta * lambdas[j+1L] /
.e10_e14[,j+1L])/.e10_e14[,j+1L]
mat[,j+1L] <- (tmp - tmp2) * .e9
}
}
c(cbind(mat, dlambdas))
}
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'ull', Theta = 'matrix'),
definition = function(x, Theta, itemexp = FALSE){
ret <- P.ull(x=x, Theta=Theta, ncat=x@ncat)
ret
}
)
# complete-data derivative used in parameter estimation
setMethod(
f = "Deriv",
signature = signature(x = 'ull', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(any(x@est)){
dp1 <- array(dP.ull(x, Theta, ncat=x@ncat), c(nrow(Theta),x@ncat,length(x@par)))
grad <- symbolicGrad_par(x, Theta, dp1=dp1)
if(estHess){
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,Theta=Theta,
gradient=FALSE)
}
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
ret
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'ull', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'ull', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
# experimental itemtype (used as a template to create custom IRT models). Edit this if you want
# to experiment with your own customized IRT models
setClass("experimental", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'experimental'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'experimental'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'experimental'),
definition = function(x){
x@par[seq_len(x@nfact)] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'experimental'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'experimental'),
definition = function(x){
par <- c(rlnorm(1, .2, .2), rnorm(1))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'experimental'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'experimental', Theta = 'matrix'),
definition = function(x, Theta){
a <- x@par[1L]
b <- x@par[2L]
p <- exp(a * (Theta - b)) / (1 + exp(a * (Theta - b)))
p <- ifelse(p < 1e-10, 1e-10, p) #numerical constraints to avoid log() problems
p <- ifelse(p > 1 - 1e-10, 1 - 1e-10, p)
P <- cbind(1-p, p)
return(P)
}
)
# complete-data derivative used in parameter estimation (here it is done numerically)
setMethod(
f = "Deriv",
signature = signature(x = 'experimental', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(any(x@est)){
grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta)
if(estHess){
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
}
}
return(list(grad=grad, hess=hess)) #TODO replace with analytic derivatives
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'experimental', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'experimental', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# defines how the item should be initiallized with the S4 function new()
# before the parameter estimates are updated
# (set starting values, lower/upper bounds, logicals indicating estimation, etc)
setMethod("initialize",
'experimental',
function(.Object, nfact, ncat){
stopifnot(nfact == 1L)
stopifnot(ncat == 2L)
.Object@par <- c(a=1, b=0)
.Object@est <- c(TRUE, TRUE)
.Object@parnames <- names(.Object@par)
.Object@lbound <- rep(-Inf, 2L)
.Object@ubound <- rep(Inf, 2L)
.Object
}
)
philchalmers/mirt documentation built on April 14, 2024, 6:41 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions dP.ull P.ull d2P.cll dP.cll P.cll d2P.5PL dP.5PL P.5PL P.ggum P.lca P.ideal P.nestlogit P.comp nominalParDeriv P.nominal P.rsm P.gpcm P.rating P.poly P.mirt Use_R_Deriv Use_R_ProbTrace ordinal_itemtypes Valid_iteminputs Experimental_itemtypes
# valid itemtype inputs
# flag to indicate an experimental item type (requires an S4 initializer in the definitions below)
# note: cannot match Valid_iteminputs
Experimental_itemtypes <- function() c('experimental', 'grsmIRT', 'fivePL', 'cll', 'ull')
Valid_iteminputs <- function() c('Rasch', '2PL', '3PL', '3PLu', '4PL', '5PL', 'CLL', 'ULL',
'graded', 'grsm', 'gpcm', 'gpcmIRT',
'rsm', 'nominal', 'PC2PL','PC3PL', '2PLNRM', '3PLNRM', '3PLuNRM', '4PLNRM',
'ideal', 'lca', 'spline', 'monopoly', 'ggum', 'sequential', 'Tutz', Experimental_itemtypes())
ordinal_itemtypes <- function() c('dich', 'fivePL', 'graded', 'gpcm', 'sequential', 'cll', 'ull',
'ggum', 'rating', 'spline', 'monopoly',
'partcomp', 'rsm', 'ideal', 'gpcmIRT', 'grsmIRT')
# Indicate which functions should use the R function instead of those written in C++
Use_R_ProbTrace <- function() c('custom', 'spline', 'sequential', 'Tutz', Experimental_itemtypes())
Use_R_Deriv <- function() c('custom', 'rating', 'partcomp', 'nestlogit', 'spline', 'sequential', 'Tutz',
Experimental_itemtypes())
#--------------------------------------------------------------------
# Item model definitions
setClass("dich", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'dich'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'dich'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'dich'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'dich'),
definition = function(x){
par <- x@par
d <- par[x@nfact+1L]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'dich'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
rnorm(1L),
rnorm(1L, -2, .5),
rnorm(1L, 2, .5))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'dich'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.mirt <- function(par, Theta, asMatrix = FALSE, ot = 0)
{
return(.Call("traceLinePts", par, Theta, ot))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'dich', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
P <- P.mirt(x@par, Theta=Theta, ot=ot)
return(P)
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'dich', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call('dparsDich', x@par, Theta, estHess, x@dat, offterm)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'dich', Theta = 'matrix'),
definition = function(x, Theta){
N <- nrow(Theta)
nfact <- ncol(Theta)
parlength <- length(x@par)
u <- antilogit(x@par[parlength])
g <- antilogit(x@par[parlength - 1L])
d <- x@par[parlength - 2L]
a <- x@par[seq_len(nfact)]
Pstar <- P.mirt(c(a, d, -999, 999), Theta)[,2L]
grad <- hess <- vector('list', 2L)
grad[[1L]] <- grad[[2L]] <- hess[[1L]] <- hess[[2L]] <- matrix(0, N, nfact)
for(i in seq_len(nfact)){
grad[[2L]][ ,i] <- (u-g) * a[i] * (Pstar * (1 - Pstar))
grad[[1L]][ ,i] <- -1 * grad[[2L]][ ,i]
hess[[2L]][ ,i] <- 2 * (u - g) * a[i]^2 * ((1 - Pstar)^2 * Pstar) -
(u - g) * a[i]^2 * (Pstar * (1 - Pstar))
hess[[1L]][ ,i] <- -1 * hess[[2L]][ ,i]
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'dich', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta)
PQ <- apply(P, 1L, prod)
ret <- cbind(Theta * PQ, PQ, P)
ret
}
)
# ----------------------------------------------------------------
setClass("graded", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'graded'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'graded'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'graded'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'graded'),
definition = function(x){
par <- x@par
d <- par[-seq_len(x@nfact)]
d
}
)
setMethod(
f = "CheckIntercepts",
signature = signature(x = 'graded'),
definition = function(x){
z <- ExtractZetas(x)
return(all(sort(z, decreasing = TRUE) == z))
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'graded'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'graded'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.poly <- function(par, Theta, itemexp = FALSE, ot = 0)
{
return(.Call('gradedTraceLinePts', par, Theta, itemexp, ot, FALSE))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'graded', Theta = 'matrix'),
definition = function(x, Theta, itemexp = TRUE, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.poly(x@par, Theta=Theta, itemexp=itemexp, ot=ot))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'graded', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call("dparsPoly", x@par, Theta, offterm, x@dat,
length(x@par) - ncol(Theta), estHess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'graded', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
P <- ProbTrace(x, Theta, itemexp = FALSE)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(ncol(P)-1L)){
w1 <- P[,i] * (1-P[,i]) * a[j]
w2 <- P[,i+1L] * (1-P[,i+1L]) * a[j]
grad[[i]][ ,j] <- w1 - w2
hess[[i]][ ,j] <- a[j]^2 * (2 * P[ ,i] * (1 - P[,i])^2 -
P[ ,i] * (1 - P[,i]) -
2 * P[ ,i+1L] * (1 - P[,i+1L])^2 +
P[ ,i+1L] * (1 - P[,i+1L]))
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'graded', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta, itemexp=FALSE)
P <- P[,-c(1, ncol(P)), drop=FALSE]
PQ <- P * (1 - P)
ret <- cbind(Theta * rowSums(PQ), PQ)
ret
}
)
# ----------------------------------------------------------------
setClass("rating", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'rating'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'rating'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'rating'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'rating'),
definition = function(x){
par <- x@par
d <- par[-c(seq_len(x@nfact), length(par))]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'rating'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE),
rnorm(1L))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'rating'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.rating <- function(par, Theta, itemexp = FALSE, ot = 0)
{
return(.Call('gradedTraceLinePts', par, Theta, itemexp, ot, TRUE))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'rating', Theta = 'matrix'),
definition = function(x, Theta, itemexp = TRUE, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.rating(x@par, Theta=Theta, itemexp=itemexp, ot=ot))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'rating', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
hess <- matrix(0, length(x@par), length(x@par))
dat <- x@dat
nfact <- x@nfact
a <- x@par[seq_len(nfact)]
d <- ExtractZetas(x)
nzetas <- length(d)
shiftind <- length(x@par)
shift <- x@par[shiftind]
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
par <- c(a, d + shift)
P <- P.poly(par=par, Theta=Theta, ot=offterm)
ret <- .Call("dparsPoly", par, Theta, offterm, x@dat,
length(d), estHess)
grad <- ret$grad
hess <- ret$hess
hess <- cbind(hess, rep(0, nrow(hess)))
hess <- rbind(hess, rep(0, ncol(hess)))
dc <- numeric(1L)
Pfull <- P
PQfull <- Pfull * (1-Pfull)
P <- P.poly(c(a, d + shift), Theta, itemexp=TRUE, ot=offterm)
rs <- dat
for(i in seq_len(ncol(rs)))
dc <- dc + rs[,i]/P[,i] * (PQfull[,i] - PQfull[,i+1L])
dc <- sum(dc)
grad <- c(grad, dc)
if(estHess){
cind <- ncol(hess)
ddc <- numeric(nrow(P))
for(i in seq_len(ncol(rs)))
ddc <- ddc + rs[,i]/P[,i] * (Pfull[,i] - 3*Pfull[,i]^2 + 2*Pfull[,i]^3 -
Pfull[,i+1L] + 3*Pfull[,i+1L]^2 - 2*Pfull[,i+1L]^3) -
rs[,i]/P[,i]^2 * (PQfull[,i] - PQfull[,i+1L])^2
hess[cind, cind] <- sum(ddc)
for(i in seq_len(nzetas))
hess[cind, nfact + i] <- hess[nfact + i, cind] <-
sum((rs[,i]/P[,i] * (-Pfull[,i+1L] + 3*Pfull[,i+1L]^2 - 2*Pfull[,i+1L]^3) -
rs[,i]/P[,i]^2 * (PQfull[,i] - PQfull[,i+1L]) * (-PQfull[,i+1L]) +
rs[,i+1L]/P[,i+1L] * (Pfull[,i+1L] - 3*Pfull[,i+1L]^2 + 2*Pfull[,i+1L]^3) -
rs[,i+1L]/P[,i+1L]^2 * (PQfull[,i+1L] - PQfull[,i+2L]) * (PQfull[,i+1L])))
for(j in seq_len(nfact)){
tmp <- 0
for(i in seq_len(ncol(rs)))
tmp <- tmp + (rs[,i]/P[,i] * Theta[,j] *
(Pfull[,i] - 3*Pfull[,i]^2 + 2*Pfull[,i]^3 -
Pfull[,i+1L] + 3*Pfull[,i+1L]^2 - 2*Pfull[,i+1L]^3) -
rs[,i]/P[,i]^2 * (PQfull[,i] - PQfull[,i+1L]) * Theta[,j] *
(PQfull[,i] - PQfull[,i+1L]))
hess[cind, j] <- hess[j, cind] <- sum(tmp)
}
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'rating', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
P <- ProbTrace(x, Theta, itemexp = FALSE)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(ncol(P)-1L)){
w1 <- P[,i] * (1-P[,i]) * a[j]
w2 <- P[,i+1L] * (1-P[,i+1L]) * a[j]
grad[[i]][ ,j] <- w1 - w2
hess[[i]][ ,j] <- a[j]^2 * (2 * P[ ,i] * (1 - P[,i])^2 -
P[ ,i] * (1 - P[,i]) -
2 * P[ ,i+1L] * (1 - P[ ,i+1L])^2 +
P[ ,i+1L] * (1 - P[ ,i+1L]))
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'rating', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
setClass("gpcm", contains = 'AllItemsClass',
representation = representation(mat='logical'))
setMethod(
f = "print",
signature = signature(x = 'gpcm'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'gpcm'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'gpcm'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'gpcm'),
definition = function(x){
par <- x@par
d <- par[(length(par)-x@ncat+1L):length(par)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'gpcm'),
definition = function(x){
ns <- sum(grepl('ak', x@parnames))
par <- c(rlnorm(x@nfact, meanlog=-.2, sdlog=.5), rep(NA, ns), 0,
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'gpcm'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.gpcm <- function(par, Theta, ot = 0, mat = FALSE, returnNum = FALSE)
{
return(.Call("gpcmTraceLinePts", par, Theta, ot, FALSE, mat, returnNum))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'gpcm', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.gpcm(x@par, Theta=Theta, ot=ot, mat=x@mat))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'gpcm', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call("dparsNominal", x, Theta, offterm, FALSE, estHess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'gpcm', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
d <- ExtractZetas(x)
ak <- 0:(x@ncat - 1L)
P <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta)
Num <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta, returnNum = TRUE)
Den <- rowSums(Num)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(x@ncat)){
grad[[i]][ ,j] <- ak[i] * a[j] * P[ ,i] - P[ ,i] * (Num %*% (ak * a[j])) / Den
hess[[i]][ ,j] <- ak[i]^2 * a[j]^2 * P[ ,i] -
2 * ak[i] * a[j] * P[,i] * (Num %*% (ak * a[j])) / Den +
2 * P[,i] * ((Num %*% (ak * a[j])) / Den)^2 -
P[,i] * ((Num %*% (ak^2 * a[j]^2)) / Den)
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'gpcm', Theta = 'matrix'),
definition = function(x, Theta){
nfact <- ncol(Theta)
ncat <- x@ncat
num <- P.gpcm(x@par, Theta=Theta, returnNum = TRUE, mat = x@mat)
den <- rowSums(num)
P <- num/den
a <- x@par[seq_len(ncol(Theta))]
if(!x@mat){
ak <- matrix(x@par[(ncol(Theta)+1L):(ncol(Theta)+x@ncat)], nfact, ncat, byrow=TRUE)
} else {
ak <- matrix(x@par[(nfact+1L):(nfact+ncat*nfact)], nfact, ncat, byrow=TRUE)
}
dp <- matrix(0, nrow(Theta), length(x@par))
aknum <- eak <- vector('list', nfact)
e <- 0:(x@ncat-1)
for(i in seq_len(nfact)){
aknum[[i]] <- t(ak[i,] * t(num))
eak[[i]] <- e*ak[i,]
}
for(i in seq_len(nfact)){
for(j in 2L:ncat)
dp[,i] <- dp[,i] + eak[[i]][j]*Theta[,i]*P[,j] -
e[j]*P[,j]*rowSums(Theta[,i] * aknum[[i]])
}
for(j in seq_len(x@ncat)){
dp[,nfact + ncat + j] <- e[j] * P[,j] - e[j] * P[,j]^2 -
as.numeric(e[-j] %*% t(P[,-j]*P[,j]))
}
return(dp)
}
)
# ----------------------------------------------------------------
setClass("rsm", contains = 'AllItemsClass',
representation = representation(mat='logical'))
setMethod(
f = "print",
signature = signature(x = 'rsm'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'rsm'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'rsm'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'rsm'),
definition = function(x){
par <- x@par
d <- par[(length(par) - x@ncat):length(par)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'rsm'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5), 0:(x@ncat-1), 0,
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE),
rnorm(1L))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'rsm'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.rsm <- function(par, Theta, ot = 0, mat = FALSE, returnNum = FALSE)
{
return(.Call("gpcmTraceLinePts", par, Theta, ot, TRUE, mat, returnNum))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'rsm', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.rsm(x@par, Theta=Theta, ot=ot))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'rsm', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
dat <- x@dat
# nfact <- x@nfact
# nzetas <- ncol(dat)
a <- ExtractLambdas(x)
d <- ExtractZetas(x)
shift <- d[length(d)]
dshift <- d <- d[-length(d)]
dshift[-1L] <- d[-1L] + shift
ak <- 0:(length(d)-1L)
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
# P <- ProbTrace(x=x, Theta=Theta, useDesign = FALSE, ot=offterm)
tmp <- .Call("dparsNominal", x, Theta, offterm, TRUE, estHess)
# num <- P.nominal(c(a, ak, dshift), ncat=length(ak),
# Theta=Theta, returnNum=TRUE, ot=offterm)
grad <- tmp$grad
hess <- tmp$hess
#quick calcs for derivs
# nfact <- length(a)
# ncat <- length(d)
# akind <- nfact
# dind <- nfact + ncat*2 + 1L #go backwards
# ak2 <- ak^2
# P2 <- P^2
# P3 <- P^3
# aTheta <- as.vector(Theta %*% a)
# aTheta2 <- aTheta^2
# dat_num <- dat/num
# numsum <- rowSums(num)
# numD <- num %*% c(0, rep(1, ncol(num)-1L))
# numak <- matrix(num %*% ak, nrow(Theta), ncol(Theta))
# numakThetaD <- numak * Theta
# numD2 <- num %*% c(0, rep(1, ncol(num)-1L))
# numakThetaD2 <- numak * Theta
# ak0 <- ak
# ak0[1L] <- 0
# cind <- length(grad)
# if(estHess){
# tmp <- 0
# for(i in 1L:nzetas)
# tmp <- tmp + dat[,i]*numD^2 / numsum^2 - dat[,i]*numD2/numsum
# hess[cind, cind] <- sum(tmp)
# for(j in 1L:nzetas){
# tmp <- 0
# for(i in 1L:nzetas)
# tmp <- tmp + dat[,i]*P[,j]*numD/numsum - dat[,i]*P[,j]
# hess[cind, nfact+j] <- hess[nfact+j, cind] <- sum(tmp)
# }
# for(j in 1L:nfact){
# tmp <- 0
# for(i in 1L:nzetas)
# tmp <- tmp + dat[,i]*numD*numakThetaD[,j]/numsum^2 -
# dat[,i]* (num %*% ak0*Theta[,j])/numsum
# hess[cind, j] <- hess[j, cind] <- sum(tmp)
# }
# }
####
#TODO - can't seem to get the last value of the gradient quite right for some reason....
x2 <- x
x2@est <- c(rep(FALSE, length(x2@est)-1L), TRUE)
grad[x2@est] <- numerical_deriv(x@par[x2@est], EML, obj=x2, Theta=Theta)
if(estHess && any(x@est)){
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
}
####
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
ret
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'rsm', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
d <- ExtractZetas(x)
t <- d[length(d)]
d <- d[-length(d)]
d[-1L] <- d[-1L] + t
ak <- 0:(x@ncat - 1L)
P <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta)
Num <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta, returnNum = TRUE)
Den <- rowSums(Num)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(x@ncat)){
grad[[i]][ ,j] <- ak[i] * a[j] * P[ ,i] - P[ ,i] *
(Num %*% (ak * a[j])) / Den
hess[[i]][ ,j] <- ak[i]^2 * a[j]^2 * P[ ,i] -
2 * ak[i] * a[j] * P[,i] * (Num %*% (ak * a[j])) / Den +
2 * P[,i] * ((Num %*% (ak * a[j])) / Den)^2 -
P[,i] * ((Num %*% (ak^2 * a[j]^2)) / Den)
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'rsm', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
setClass("nominal", contains = 'AllItemsClass',
representation = representation(mat='logical'))
setMethod(
f = "print",
signature = signature(x = 'nominal'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'nominal'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'nominal'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'nominal'),
definition = function(x){
d <- x@par[(length(x@par) - x@ncat + 1L):length(x@par)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'nominal'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=-.2, sdlog=.5), 0,
abs(rnorm(x@ncat-1L, (x@ncat-1L) / 2, sd = 1)), 0,
rnorm(x@ncat-1L))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'nominal'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x@est[(x@nfact+1L):(x@nfact + x@ncat)] <- FALSE
x
}
)
P.nominal <- function(par, ncat, Theta, returnNum = FALSE, ot = 0){
return(.Call("nominalTraceLinePts", par, ncat, Theta, returnNum, ot))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'nominal', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(x@fixed.design) > 1L && useDesign)
Theta <- cbind(x@fixed.design, Theta)
return(P.nominal(par=x@par, ncat=x@ncat, Theta=Theta, ot=ot))
}
)
nominalParDeriv <- function(a, ak, d, Theta, P, num, dat, estHess, gpcm = FALSE){
nfact <- length(a)
ncat <- length(d)
akind <- nfact
dind <- nfact + ncat
ak2 <- ak^2
P2 <- P^2
P3 <- P^3
aTheta <- as.vector(Theta %*% a)
aTheta2 <- aTheta^2
dat_num <- dat/num
numsum <- rowSums(num)
numakD <- num %*% ak
numak2D2 <- num %*% ak2
numakDTheta_numsum <- matrix(0, nrow(num), nfact)
for(i in seq_len(nfact))
numakDTheta_numsum[,i] <- (num %*% ak * Theta[, i])/ numsum
ret <- .Call('dparsNominal', a, ak, d, Theta, P, num, dat, nfact, ncat,
akind, dind, ak2, P2, P3, aTheta, aTheta2, dat_num, numsum, numakD,
numak2D2, numakDTheta_numsum, estHess)
ret
}
setMethod(
f = "Deriv",
signature = signature(x = 'nominal', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call("dparsNominal", x, Theta, offterm, FALSE, estHess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'nominal', Theta = 'matrix'),
definition = function(x, Theta){
a <- ExtractLambdas(x)
d <- ExtractZetas(x)
ak <- x@par[(x@nfact+1):(x@nfact+x@ncat)]
P <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta)
Num <- P.nominal(c(a, ak, d), ncat=length(d), Theta=Theta, returnNum = TRUE)
Den <- rowSums(Num)
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(x@ncat)){
grad[[i]][ ,j] <- ak[i] * a[j] * P[ ,i] - P[ ,i] *
(Num %*% (ak * a[j])) / Den
hess[[i]][ ,j] <- ak[i]^2 * a[j]^2 * P[ ,i] -
2 * ak[i] * a[j] * P[,i] * (Num %*% (ak * a[j])) / Den +
2 * P[,i] * ((Num %*% (ak * a[j])) / Den)^2 -
P[,i] * ((Num %*% (ak^2 * a[j]^2)) / Den)
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'nominal', Theta = 'matrix'),
definition = function(x, Theta){
num <- P.nominal(x@par, ncat=x@ncat, Theta=Theta, returnNum=TRUE)
den <- rowSums(num)
P <- num/den
a <- x@par[seq_len(ncol(Theta))]
ak <- x@par[(ncol(Theta)+1L):(ncol(Theta)+x@ncat)]
dp <- matrix(0, nrow(Theta), length(x@par))
aknum <- t(ak * t(num))
aTheta <- as.numeric(a %*% t(Theta))
nfact <- ncol(Theta)
ncat <- x@ncat
e <- 0:(x@ncat-1)
eak <- e*ak
for(i in seq_len(nfact)){
for(j in 2L:ncat)
dp[,i] <- dp[,i] + eak[j]*Theta[,i]*P[,j] -
e[j]*P[,j]*rowSums(Theta[,i] * aknum)
}
for(j in seq_len(x@ncat)){
dp[,nfact + ncat + j] <- e[j] * P[,j] - e[j] * P[,j]^2 -
as.numeric(e[-j] %*% t(P[,-j]*P[,j]))
dp[,nfact + j] <- dp[,nfact + ncat + j] * aTheta
}
return(dp)
}
)
# ----------------------------------------------------------------
setClass("partcomp", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'partcomp'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'partcomp'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'partcomp'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'partcomp'),
definition = function(x){
d <- x@par[(x@nfact+1L):(length(x@par)-2L)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'partcomp'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
rlnorm(x@nfact, meanlog=.2, sdlog=.5),
rnorm(1L, -2, .5),
rnorm(1L, 2, .5))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'partcomp'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.comp <- function(par, Theta, ot = 0)
{
return(.Call('partcompTraceLinePts', par, Theta))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'partcomp', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
return(P.comp(x@par, Theta=Theta))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'partcomp', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
#local derivative from previous version with small mod
#u and g in logit form
dpars.comp <- function(lambda,zeta,g,r,f,Thetas,estHess)
{
nfact <- length(lambda)
pars <- c(zeta,lambda,g)
pgrad <- function(pars, r, thetas){
nfact <- ncol(thetas)
d <- pars[seq_len(nfact)]
a <- pars[(nfact+1L):(length(pars)-1L)]
c <- pars[length(pars)]
P <- P.comp(c(a,d,c,999), thetas)[,2L]
Pstar <- P.comp(c(a,d,-999,999),thetas)[,2L]
Qstar <- 1 - Pstar
Q <- 1 - P
c <- antilogit(c)
g_1g <- c * (1 - c)
const1 <- (r/P - (f-r)/Q)
dd <- da <- rep(0,nfact)
dc <- sum(r/P * (g_1g * (1 - Pstar)) + (f-r)/Q * (g_1g * (Pstar - 1)))
for(i in seq_len(nfact)){
Pk <- P.mirt(c(a[i],d[i],-999,999),matrix(thetas[,i]))[,2L]
Qk <- 1 - Pk
dd[i] <- sum((1-c)*Pstar*Qk*const1)
da[i] <- sum((1-c)*Pstar*Qk*thetas[,i]*const1)
}
return(c(dd,da,dc))
}
phess <- function(pars, r, thetas){
nfact <- ncol(thetas)
d <- pars[seq_len(nfact)]
a <- pars[(nfact+1L):(length(pars)-1L)]
c <- pars[length(pars)]
P <- P.comp(c(a,d,c,999), thetas)[,2L]
Pstar <- P.comp(c(a,d,-999,999),thetas)[,2L]
Qstar <- 1 - Pstar
Q <- 1 - P
g <- c <- antilogit(c)
g_1g <- c * (1 - c)
const1 <- (r/P - (f-r)/Q)
const2 <- (r/P^2 + (f-r)/Q^2)
hess <- matrix(0,nfact*2+1,nfact*2+1)
dNames <- paste("d",seq_len(nfact),sep='_')
aNames <- paste("a",seq_len(nfact),sep='_')
Names <- c(paste("d",seq_len(nfact),sep='_'),paste("a",1:nfact,sep='_'),'c_0')
for(i in seq_len(nfact*2+1L)){
for(j in seq_len(nfact*2+1L)){
if(i <= j){
d1 <- strsplit(Names[c(i,j)],"_")[[1L]]
d2 <- strsplit(Names[c(i,j)],"_")[[2L]]
k <- as.numeric(d1[2L])
m <- as.numeric(d2[2L])
Pk <- P.mirt(c(a[k],d[k],-999,999),matrix(thetas[,k]))[,2L]
Qk <- 1 - Pk
Pm <- P.mirt(c(a[m],d[m],-999,999),matrix(thetas[,m]))[,2L]
Qm <- 1 - Pm
if(i == j && d1[1L] == 'd'){
hess[i,i] <- sum(r/P * ((1-g)*(Pstar - 3*Pstar*Pk + 2*Pstar*Pk^2) ) -
r/P^2 * ((1-g) * (Pstar - Pstar*Pk))^2 +
(f-r)/Q * ((1-g)*(-Pstar + 3*Pstar*Pk - 2*Pstar*Pk^2)) -
(f-r)/Q^2 * ((1-g)*(-Pstar + Pstar*Pk))^2)
next
}
if(i == j && d1[1L] == 'a'){
hess[i,i] <- sum(r/P * ((1-g)*thetas[,k]^2*(Pstar - 3*Pstar*Pk + 2*Pstar*Pk^2) ) -
r/P^2 * ((1-g)*thetas[,k] * (Pstar - Pstar*Pk))^2 +
(f-r)/Q * ((1-g)*thetas[,k]^2*(-Pstar + 3*Pstar*Pk - 2*Pstar*Pk^2)) -
(f-r)/Q^2 * ((1-g)*thetas[,k] * (-Pstar + Pstar*Pk))^2)
next
}
if(i == j && d1[1L] == 'c'){
hess[i,i] <- sum(r/P * (g_1g * (2.0*(1-g) - 1.0 - 2.0*(1-g)*Pstar + Pstar)) -
r/P^2 * (g_1g * (1.0 - Pstar)) * (g_1g * (1.0 - Pstar)) +
(f-r)/Q * (g_1g * (-2.0*(1-g) + 1.0 + 2.0*(1-g)*Pstar - Pstar)) -
(f-r)/Q^2 * (g_1g * (-1.0 + Pstar)) * (g_1g * (-1.0 + Pstar)))
next
}
if(d1[1L] == 'a' && d2[1L] == 'a'){
hess[i,j] <- hess[j,i] <- sum((1-c)*thetas[,k]*thetas[,m]*
Qk*Pstar*Qm*(const1 - Pstar*(1-c)*const2))
next
}
if(d1[1L] == 'd' && d2[1L] == 'd'){
hess[i,j] <- hess[j,i] <- sum((1-c)*Qk*Pstar*Qm*(const1 - Pstar*(1-c)*const2))
next
}
if(d1[1L] == 'a' && d2[1L] == 'c'){
hess[i,j] <- hess[j,i] <- sum(r/P*(g_1g*thetas[,k]*(-Pstar + Pstar*Pk)) -
r/P^2*((1-g)*thetas[,k]*(Pstar - Pstar*Pk)*(g_1g*(1 - Pstar))) +
(f-r)/Q * (g_1g*thetas[,k]*(Pstar - Pstar*Pk))-
(f-r)/Q^2 * ((1-g)*thetas[,k]*(-Pstar + Pstar*Pk)*(g_1g*(-1 + Pstar))))
next
}
if(d1[1L] == 'd' && d2[1L] == 'c'){
hess[i,j] <- hess[j,i] <- sum(r/P*(g_1g*(-Pstar + Pstar*Pk)) -
r/P^2*((1-g)*(Pstar - Pstar*Pk)*(g_1g*(1 - Pstar))) +
(f-r)/Q * (g_1g*(Pstar - Pstar*Pk))-
(f-r)/Q^2 * ((1-g)*(-Pstar + Pstar*Pk)*(g_1g*(-1 + Pstar))))
next
}
if(d1[1L] == 'd' && d2[1L] == 'a' && d1[2] == d2[2]){
hess[i,j] <- hess[j,i] <- sum(
r/P * ((1-g)*thetas[,k]*(Pstar - 3*Pstar*Pk + 2*Pstar*Pk^2) ) -
r/P^2 * ((1-g)*thetas[,k] * (Pstar - Pstar*Pk) * (1-g)*(Pstar - Pstar*Pk)) +
(f-r)/Q * ((1-g)*thetas[,k]*(-Pstar + 3*Pstar*Pk - 2*Pstar*Pk^2)) -
(f-r)/Q^2 * ((1-g)*thetas[,k] * (-Pstar + Pstar*Pk) * (1-g)*(-Pstar + Pstar*Pk)))
next
}
if(d1[1L] == 'd' && d2[1L] == 'a' && d1[2] != d2[2]){
hess[i,j] <- hess[j,i] <- sum(r/P * ((1-g)*thetas[,m]*(Pstar - Pstar*Pk - Pstar*Pm + Pstar^2)) -
r/P^2 * ((1-g)*thetas[,m] * (Pstar - Pstar*Pm) * (1-g)*(Pstar - Pstar*Pk)) +
(f-r)/Q * ((1-g)*thetas[,m]*(-Pstar + Pstar*Pk + Pstar*Pm - Pstar^2)) -
(f-r)/Q^2 * ((1-g)*thetas[,m] * (-Pstar + Pstar*Pm) * (1-g)*(-Pstar + Pstar*Pk)))
next
}
}
}
}
return(hess)
}
#old pars in the form d, a, g
g <- pgrad(pars, r, Thetas)
if(estHess) h <- phess(pars, r, Thetas)
else h <- matrix(0, length(g), length(g))
#translate into current version
grad <- c(g[(nfact+1L):(nfact*2)], g[1L:nfact], g[length(g)], 0)
hess <- matrix(0, ncol(h) + 1L, ncol(h) + 1L)
if(estHess){
hess[seq_len(nfact), seq_len(nfact)] <- h[(nfact+1L):(nfact*2),(nfact+1L):(nfact*2)] #a block
hess[(nfact+1L):(nfact*2),(nfact+1L):(nfact*2)] <- h[seq_len(nfact), seq_len(nfact)] #d block
hess[nfact*2 + 1L, nfact*2 + 1L] <- h[nfact*2 + 1L, nfact*2 + 1L] #g
hess[nfact*2 + 1L, 1L:nfact] <- hess[1:nfact, nfact*2 + 1L] <-
h[nfact*2 + 1L, (nfact+1L):(nfact*2)] #ga
hess[nfact*2 + 1L, (nfact+1L):(nfact*2)] <- hess[(nfact+1L):(nfact*2), nfact*2 + 1L] <-
h[nfact*2 + 1L, seq_len(nfact)] #gd
hess[(nfact+1L):(nfact*2), seq_len(nfact)] <- t(h[(nfact+1L):(nfact*2), seq_len(nfact)])
hess[seq_len(nfact), (nfact+1L):(nfact*2)] <- t(h[seq_len(nfact), (nfact+1L):(nfact*2)]) #ads
}
return(list(grad=grad, hess=hess))
}
#####
f <- rowSums(x@dat)
r <- x@dat[ ,2L]
nfact <- x@nfact
a <- x@par[seq_len(nfact)]
d <- x@par[(nfact+1L):(nfact*2L)]
g <- x@par[length(x@par)-1L]
tmp <- dpars.comp(lambda=ExtractLambdas(x),zeta=ExtractZetas(x),g=x@par[nfact*2L + 1L],r=r, f=f,
Thetas=Theta, estHess=estHess)
ret <- list(grad=tmp$grad, hess=tmp$hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'partcomp', Theta = 'matrix'),
definition = function(x, Theta){
N <- nrow(Theta)
nfact <- ncol(Theta)
parlength <- length(x@par)
u <- x@par[parlength]
g <- x@par[parlength - 1L]
d <- ExtractZetas(x)
a <- ExtractLambdas(x)
Pstar <- P.comp(c(a, d, -999, 999), Theta)[,2L]
g <- antilogit(g)
u <- antilogit(u)
grad <- hess <- vector('list', 2L)
grad[[1L]] <- grad[[2L]] <- hess[[1L]] <- hess[[2L]] <- matrix(0, N, nfact)
for(j in seq_len(nfact)){
Pn <- P.mirt(c(a[j], d[j], -999, 999), Theta)[,2L]
grad[[2L]][ ,j] <- (u-g)*Pstar*a[j]*(1 - Pn)
grad[[1L]][ ,j] <- -1*grad[[2L]][ ,j]
hess[[2L]][ ,j] <- (u-g)*a[j]^2*Pstar*(1 - 3*Pn + 2*Pn^2)
hess[[1L]][ ,j] <- -1*hess[[2L]][ ,j]
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'partcomp', Theta = 'matrix'),
definition = function(x, Theta){
# message('partcomp derivatives not optimized') ##TODO
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
setClass("nestlogit", contains = 'AllItemsClass',
representation = representation(correctcat='integer'))
setMethod(
f = "print",
signature = signature(x = 'nestlogit'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'nestlogit'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'nestlogit'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'nestlogit'),
definition = function(x){
stop('not written')
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'nestlogit'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
rnorm(1L),
rnorm(1L, -2, .5),
rnorm(1L, 2, .5), 0,
abs(rnorm(x@ncat-2L, (x@ncat-2L) / 2, sd = 1)), x@ncat,
rnorm(x@ncat-2L))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'nestlogit'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x@est[(x@nfact+5L):(x@nfact + x@ncat + 1L)] <- FALSE
x
}
)
P.nestlogit <- function(par, Theta, correct, ncat)
{
return(.Call("nestlogitTraceLinePts", par, Theta, correct, ncat))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'nestlogit', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
return(P.nestlogit(x@par, Theta=Theta, correct=x@correctcat, ncat=x@ncat))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'nestlogit', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
#u and g in logit
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
dat <- x@dat
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta,
gradient = FALSE)
nfact <- x@nfact
a <- x@par[seq_len(x@nfact)]
d <- x@par[x@nfact+1L]
g <- x@par[x@nfact+2L]
u <- x@par[x@nfact+3L]
ak <- x@par[(x@nfact+4L):(x@nfact+4L+x@ncat-2L)]
dk <- x@par[(length(x@par)-length(ak)+1):length(x@par)]
correct <- x@correctcat
Pd <- P.mirt(c(a, d, g, u), Theta=Theta)[,2L]
g <- antilogit(g)
u <- antilogit(u)
Qd <- 1 - Pd
Pstar <- P.mirt(c(a, d, -999, 999), Theta=Theta)[,2L]
Qstar <- 1 - Pstar
num <- P.nominal(c(rep(1, nfact), ak, dk), ncat=length(dk),
Theta=Theta, returnNum=TRUE)
den <- rowSums(num)
Pn <- num/den
cdat <- dat[,correct]
idat <- dat[,-correct]
nd <- ncol(idat)
g_1g <- g * (1 - g)
u_1u <- u * (1 - u)
for(i in seq_len(nfact))
grad[i] <- sum( (u-g) * Theta[,i] * Qstar * Pstar * (
cdat / Pd - rowSums(idat/Qd)) )
grad[nfact+1L] <- sum( (u-g) * Qstar * Pstar * (
cdat / Pd - rowSums(idat/Qd)) )
grad[nfact+2L] <- sum( ((cdat * g_1g * (1-Pstar)/Pd) + rowSums(idat * g_1g * (Pstar - 1)/Qd)) )
grad[nfact+3L] <- sum( (cdat * u_1u * Pstar / Pd - rowSums(idat * u_1u * Pstar / Qd) ))
for(j in seq_len(nd)){
grad[nfact+3L+j] <- sum((
(idat[,j] * Qd * rowSums(Theta) * (Pn[,j] - Pn[,j]^2)) / (Qd * Pn[,j]) -
rowSums(idat[,-j, drop=FALSE]) * rowSums(Theta) * Pn[,j]))
grad[nfact+3L+nd+j] <- sum((
(idat[,j] * Qd * (Pn[,j] - Pn[,j]^2)) / (Qd * Pn[,j]) -
rowSums(idat[,-j, drop=FALSE]) * Pn[,j]))
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'nestlogit', Theta = 'matrix'),
definition = function(x, Theta){
a <- x@par[seq_len(x@nfact)]
d <- x@par[x@nfact+1L]
g <- x@par[x@nfact+2L]
u <- x@par[x@nfact+3L]
ak <- x@par[(x@nfact+4L):(x@nfact+4L+x@ncat-2L)]
dk <- x@par[(length(x@par)-length(ak)+1):length(x@par)]
Pn <- P.nominal(c(rep(1,ncol(Theta)), ak, dk), ncat=length(dk), Theta=Theta)
Num <- P.nominal(c(rep(1,ncol(Theta)), ak, dk), ncat=length(dk),
Theta=Theta, returnNum = TRUE)
Den <- rowSums(Num)
Pstar <- P.mirt(c(a, d, -999, 999), Theta)[,2L]
Q <- P.mirt(c(a, d, g, u), Theta)[,1]
g <- antilogit(g)
u <- antilogit(u)
Num2 <- P <- matrix(0, nrow(Theta), x@ncat)
P[,-x@correctcat] <- Pn
Num2[,-x@correctcat] <- Num
Num <- Num2
ak2 <- dk2 <- numeric(x@ncat)
ak2[-x@correctcat] <- ak
dk2[-x@correctcat] <- dk
ak <- ak2
dk <- dk2
grad <- hess <- vector('list', x@ncat)
for(i in seq_len(x@ncat))
grad[[i]] <- hess[[i]] <- matrix(0, nrow(Theta), x@nfact)
for(j in seq_len(x@nfact)){
for(i in seq_len(x@ncat)){
if(i == x@correctcat){
grad[[i]][ ,j] <- (u-g) * a[j] * (Pstar * (1 - Pstar))
hess[[i]][ ,j] <- 2 * (u - g) * a[j]^2 * ((1 - Pstar)^2 * Pstar) -
(u - g) * a[j]^2 * (Pstar * (1 - Pstar))
} else {
grad[[i]][ ,j] <- -(u-g) * a[j] * (Pstar * (1 - Pstar)) * P[,i] +
Q * (ak[i] * P[ ,i] - P[ ,i] * (Num %*% (ak)) / Den)
hess[[i]][ ,j] <-
-2 * (u - g) * a[j]^2 * (1 - Pstar)^2 * Pstar * P[,i] +
(u - g) * a[j]^2 * Pstar * (1 - Pstar) * P[,i] -
2 * (u - g) * a[j] * ak[i] * (1 - Pstar) * Pstar * P[,i] +
2 * a[j] * (Pstar * (1 - Pstar)) * P[,i] * (Num %*% (ak)) / Den +
Q * ak[i]^2 * P[ ,i] -
2 * Q * ak[i] * P[,i] * (Num %*% (ak)) / Den +
2 * Q * P[,i] * ((Num %*% (ak)) / Den)^2 -
Q * P[,i] * ((Num %*% (ak^2)) / Den)
}
}
}
return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'nestlogit', Theta = 'matrix'),
definition = function(x, Theta){
# message('nestlogit derivatives not optimized') ##TODO
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
setClass("ideal", contains = 'AllItemsClass')
setMethod(
f = "print",
signature = signature(x = 'ideal'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'ideal'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'ideal'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'ideal'),
definition = function(x){
d <- x@par[length(x@par)]
d
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'ideal'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
-rlnorm(1, meanlog=0, sdlog=.5))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'ideal'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.ideal <- function(par, Theta, ot = 0)
{
alpha <- par[length(par)]
beta <- par[-length(par)]
eta <- -0.5*(Theta %*% beta + alpha)^2
eta <- ifelse(eta < -20, -20, eta)
P <- exp(eta)
P <- cbind((1-P), P)
P <- ifelse(P < 1e-20, 1e-20, P)
P <- ifelse(P > (1 - 1e-20), (1 - 1e-20), P)
return(P)
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'ideal', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
return(P.ideal(x@par, Theta=Theta))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'ideal', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
d <- x@par[length(x@par)]
a <- x@par[-length(x@par)]
P <- ProbTrace(x, Theta=Theta)[,2L]
Q <- (1 - P)
Q <- ifelse(Q < 1e-7, 1e-7, Q)
int <- as.numeric(Theta %*% a + d)
for(i in seq_len(ncol(Theta)))
grad[i] <- -sum( x@dat[,1] * int * Theta[,i] * -P / Q +
x@dat[,2] * int * Theta[,i])
grad[i+1L] <- -sum(2 * x@dat[,1] * int * -P / Q +
2 * x@dat[,2] * int)/2
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
return(list(grad = grad, hess=hess))
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'ideal', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta)
# N <- nrow(Theta)
# nfact <- ncol(Theta)
# P <- ProbTrace(x, Theta=Theta)[,2L]
# d <- x@par[length(x@par)]
# a <- x@par[-length(x@par)]
# int <- as.numeric(t(a %*% t(Theta)) + d)
# grad <- hess <- vector('list', 2L)
# grad[[1L]] <- grad[[2L]] <- hess[[1L]] <- hess[[2L]] <- matrix(0, N, nfact)
# for(i in 1L:nfact){
# grad[[2L]][ ,i] <- 2 * a[i] * int * P
# grad[[1L]][ ,i] <- -1 * grad[[2L]][ ,i]
# hess[[2L]][ ,i] <- 2 * a[i]^2 * int * P + 4 * int^2 * a[i]^2 * P
# hess[[1L]][ ,i] <- -1 * hess[[2L]][ ,i]
# }
# return(list(grad=grad, hess=hess))
}
)
setMethod(
f = "dP",
signature = signature(x = 'ideal', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta=Theta)[,2L]
dp <- matrix(0, nrow(Theta), length(x@par))
d <- x@par[length(x@par)]
a <- x@par[-length(x@par)]
int <- as.numeric(t(a %*% t(Theta)) + d)
for(i in seq_len(ncol(Theta)))
dp[,i] <- -2 * Theta[,i] * int * P
dp[,i+1L] <- -2 * int * P
dp
}
)
# ----------------------------------------------------------------
setClass("lca", contains = 'AllItemsClass',
representation = representation(score='numeric',
item.Q='matrix'))
setMethod(
f = "print",
signature = signature(x = 'lca'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'lca'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'lca'),
definition = function(x){
x@par[seq_len(x@nfact)]
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'lca'),
definition = function(x){
stop('not written')
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'lca'),
definition = function(x){
par <- rnorm(length(x@par))
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'lca'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.lca <- function(par, Theta, ncat, item.Q, returnNum = FALSE){
return(.Call('lcaTraceLinePts', par, Theta, item.Q, ncat, returnNum))
}
setMethod(
f = "ProbTrace",
signature = signature(x = 'lca', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
P <- P.lca(x@par, Theta=Theta, ncat=x@ncat, item.Q=x@item.Q)
return(P)
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'lca', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
# ret <- .Call('dparslca', x@par, Theta, x@item.Q,
# estHess, x@dat, offterm) #TODO use estHess in Rcpp
ret <- list(grad=numeric(length(x@par)),
hess=matrix(0, length(x@par), length(x@par)))
ret$grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta)
if(estHess && any(x@est)){
ret$hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
}
# if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'lca', Theta = 'matrix'),
definition = function(x, Theta){
stop('not written')
}
)
setMethod(
f = "dP",
signature = signature(x = 'lca', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta)
dp <- matrix(0, nrow(Theta), length(x@par))
ind <- 1L
for(j in 2L:x@ncat){
for(i in seq_len(ncol(Theta))){
dp[,ind] <- Theta[,i] * (P[,j] -
rowSums(P[,j,drop=FALSE] * P[,j,drop=FALSE]))
ind <- ind + 1L
}
}
dp
}
)
# ----------------------------------------------------------------
setClass("spline", contains = 'AllItemsClass',
representation = representation(stype='character',
Theta_prime='matrix',
item.Q='matrix',
sargs='list'))
setMethod(
f = "print",
signature = signature(x = 'spline'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'spline'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'spline'),
definition = function(x){
numeric(x@nfact)
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'spline'),
definition = function(x){
stop('not written')
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'spline'),
definition = function(x){
par <- rnorm(length(x@par), sd = 20)
x@par[x@est] <- par[x@est]
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'spline'),
definition = function(x){
stop('spline null should not be run')
}
)
setMethod(
f = "ProbTrace",
signature = signature(x = 'spline', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(nrow(Theta) != nrow(x@Theta_prime))
x <- loadSplineParsItem(x, Theta)
P <- P.lca(x@par, item.Q=x@item.Q,
Theta=x@Theta_prime, ncat=x@ncat)
return(P)
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'spline', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
ret <- .Call('dparslca', x@par, x@Theta_prime, x@item.Q, FALSE, x@dat, offterm)
if(estHess && any(x@est)){
ret$hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=x@Theta_prime, gradient=FALSE)
}
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'spline', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta)
}
)
setMethod(
f = "dP",
signature = signature(x = 'spline', Theta = 'matrix'),
definition = function(x, Theta){
P <- ProbTrace(x, Theta)
dp <- matrix(0, nrow(Theta), length(x@par))
ind <- 1L
for(j in 2L:x@ncat){
for(i in seq_len(ncol(Theta))){
dp[,ind] <- Theta[,i] * (P[,j] -
rowSums(P[,j,drop=FALSE] * P[,j,drop=FALSE]))
ind <- ind + 1L
}
}
dp
}
)
# ----------------------------------------------------------------
setClass("ggum", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'ggum'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'ggum'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'ggum'),
definition = function(x){
x@par[1L:x@nfact] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'ggum'),
definition = function(x){
x@par[(x@nfact+1L):length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'ggum'),
definition = function(x){
par <- c(rlnorm(x@nfact, .2, .2),
rnorm(x@nfact),
sort(rnorm(x@ncat-1L, 0, 1.5), decreasing = TRUE))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'ggum'),
definition = function(x){
x@par[1L:x@nfact] <- 0
x@est[1L:x@nfact] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'ggum', Theta = 'matrix'),
definition = function(x, Theta){
return(P.ggum(x@par, Theta=Theta, ncat=x@ncat))
}
)
# TODO write this with Rcpp
P.ggum <- function(par, Theta, ncat)
{
# C <- ncat - 1
# D <- (length(par)-C)/2
# M <- 2*C + 1
#
# sumtau <- 0
# sumdist<-0
# num <- matrix(nrow=nrow(Theta),ncol=(C+1))
# P <- matrix(nrow=nrow(Theta),ncol=(C+1))
#
# for (d in 1:D) {
# dist <- (par[d]**2)*((Theta[,d] - par[(D+d)])**2)
# sumdist <- dist+sumdist
# }
# dist<-sqrt(sumdist)
# z1 <- z2 <- num
#
# for (w in 0:C) {
#
# x1 <- w*dist
# x2 <- (M-w)*dist
#
# if (w>0) {
# for (d in 1:D) {
# tau <- (par[d]*par[(w+2*D)])
# sumtau <- tau + sumtau
# }
# }
#
# z1[,w + 1] <- x1 + sumtau
# z2[,w + 1] <- x2 + sumtau
#
# # num1 <- exp(x1 + sumtau)
# # num2 <- exp(x2 + sumtau)
# #
# # num[,(w+1)] <- num1 + num2
#
# }
#
# # numtot <- apply(num, 1, sum)
# #
# # for (k in 1:(C+1)) {
# # P[,k] <- num[,k]/numtot
# # }
#
# # less overflow this way
# maxz1 <- apply(z1, 1, max)
# maxz2 <- apply(z2, 1, max)
# maxz <- pmax(maxz1, maxz2)
# num2 <- exp(z1 - maxz) + exp(z2 - maxz)
# den <- rowSums(num2)
# P <- num2/den
# P <- ifelse(P < 1e-20, 1e-20, P)
# P <- ifelse(P > (1 - 1e-20), (1 - 1e-20), P)
# return(P)
return(.Call("ggumTraceLinePts", par, Theta, ncat))
}
# complete-data derivative used in parameter estimation (here it is done numerically)
setMethod(
f = "Deriv",
signature = signature(x = 'ggum', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
# NOT USED, but useful for numerical checking
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta,
type='Richardson')
if(estHess){
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta,
gradient=FALSE, type='Richardson')
}
return(list(grad=grad, hess=hess))
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'ggum', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'ggum', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass('custom', contains = 'AllItemsClass',
representation = representation(name='character',
P='function',
gr='function',
usegr='logical',
hss='function',
gen='function',
usehss='logical',
userdata='list',
useuserdata='logical',
derivType='character',
derivType.hss='character'))
setMethod(
f = "print",
signature = signature(x = 'custom'),
definition = function(x, ...){
cat('Custom item object named:', x@name)
}
)
setMethod(
f = "show",
signature = signature(object = 'custom'),
definition = function(object){
print(object)
}
)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'custom'),
definition = function(x){
a <- rep(.001, x@nfact)
a
}
)
setMethod(
f = "ExtractZetas",
signature = signature(x = 'custom'),
definition = function(x){
stop('not written')
}
)
setMethod(
f = "GenRandomPars",
signature = signature(x = 'custom'),
definition = function(x){
x@par <- x@gen(x)
x
}
)
setMethod(
f = "set_null_model",
signature = signature(x = 'custom'),
definition = function(x){
stop('calculating null models for custom item types is ambiguous.',
call. = FALSE)
}
)
setMethod(
f = "ProbTrace",
signature = signature(x = 'custom', Theta = 'matrix'),
definition = function(x, Theta, useDesign = TRUE, ot=0){
if(x@useuserdata) return(x@P(x@par, Theta, x@ncat, x@userdata[[1L]]))
else return(x@P(x@par, Theta, x@ncat))
}
)
setMethod("initialize",
'custom',
function(.Object, name, par, est, parnames, lbound, ubound,
P, gr, hss, gen, userdata, derivType, derivType.hss, dps, dps2) {
dummyfun <- function(...) return(NULL)
names(est) <- names(par)
usegr <- usehss <- useuserdata <- TRUE
.Object@name <- name
.Object@par <- par
.Object@est <- est
.Object@parnames <- parnames
.Object@P <- P
.Object@derivType <- derivType
.Object@derivType.hss <- derivType.hss
.Object@itemclass <- 999L
.Object@dps <- dps
.Object@dps2 <- dps2
if(is.null(gr)){
.Object@gr <- dummyfun
usegr <- FALSE
} else .Object@gr <- gr
if(is.null(hss)){
.Object@hss <- dummyfun
usehss <- FALSE
} else .Object@hss <- hss
if(is.null(gen)){
.Object@gen <- function(object) object@par
} else .Object@gen <- gen
if(is.null(userdata)){
.Object@userdata <- list()
useuserdata <- FALSE
} else .Object@userdata <- userdata
.Object@usegr <- usegr
.Object@usehss <- usehss
.Object@useuserdata <- useuserdata
.Object@lbound <- if(!is.null(lbound)) lbound else rep(-Inf, length(par))
.Object@ubound <- if(!is.null(ubound)) ubound else rep(Inf, length(par))
.Object
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'custom', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(x@usegr) grad <- x@gr(x, Theta)
else grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, type=x@derivType)
if(estHess){
if(x@usehss) hess <- x@hss(x, Theta)
else hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, type=x@derivType.hss,
gradient=FALSE)
}
return(list(grad = grad, hess=hess))
}
)
setMethod(
f = "DerivTheta",
signature = signature(x = 'custom', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta)
}
)
setMethod(
f = "dP",
signature = signature(x = 'custom', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta)
}
)
# ----------------------------------------------------------------
# itemtype='grsmIRT', 1-dim graded rating scale
# slope, thresholds and 1 location parameters
# ExtractZetas, GenRandomPar, ProbTrace, initialize MODIFIED
setClass("grsmIRT", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'grsmIRT'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'grsmIRT'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'grsmIRT'),
definition = function(x){
x@par[seq_len(x@nfact)] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'grsmIRT'),
definition = function(x){
#x@par[(x@nfact+1):(length(x@par))] #intercepts
x@par[(x@nfact+1):(x@ncat+1)] #intercepts
# if we set :(length(x@par)), when collapsed, may produce error?
# MAYBE DOESNT MATTER?
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'grsmIRT'),
definition = function(x){
par <- c(rlnorm(1,0,1),sort(rnorm(x@ncat-1, 0, 1), decreasing=TRUE), rnorm(1,0,0.5))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'grsmIRT'),
definition = function(x){
x@par[seq_len(x@nfact+x@ncat)] <- 0
x@est[seq_len(x@nfact+x@ncat)] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'grsmIRT', Theta = 'matrix'),
definition = function(x, Theta){
th1 = Theta[,1];
ncat = x@ncat
a1 = x@par[1]
d = x@par[2:ncat]
c = x@par[ncat+1]
nr = nrow(Theta); nc=ncat-1
if (all(d == sort(d, decreasing=TRUE))) {
D.star = matrix(c+d, nrow=nr, ncol=nc, byrow=TRUE)
TH1 = matrix(th1, nrow=nr, ncol=nc)
A = matrix(a1, nrow=nr, ncol=nc)
P.star = 1/(1+exp(-A*(D.star+TH1)))
P=cbind(1, P.star)-cbind(P.star, 0)
P <- ifelse(P < 1e-20, 1e-20, P)
P <- ifelse(P > (1 - 1e-20), (1 - 1e-20), P)
} else {
P <- matrix(1e-20, nrow=nr, ncol=ncat)}
return(P)
}
)
# complete-data derivative used in parameter estimation
# Analytic Derivation is provided...
# and it fails, numerical derivation is used.
setMethod(
f = "Deriv",
signature = signature(x = 'grsmIRT', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
dat <- x@dat;
th1 <- Theta[,1];
ncat <- x@ncat
a1 <- x@par[1]
d <- x@par[2:ncat]
c <- x@par[ncat+1]
nr <- length(th1);
nc <- ncat - 1
beta.mean <- mean(d)
beta.mat <- matrix(d, nrow=nr, ncol=nc, byrow=TRUE)
Sigma.mat <- (beta.mat - beta.mean)
th1.mat <- matrix(th1, nr, nc)
B.mat <- th1.mat + c + beta.mat
lam.star.mat <- exp(-a1*B.mat)
lam.mat <- cbind(lam.star.mat, 1)-cbind(0, lam.star.mat)
P.star <- 1/(1+lam.star.mat)
Q.star <- 1- P.star
P <- P.star[,-nc] - P.star[,-1]
inv.P <- 1/P
#when a1*B.mat >30, inv.P would be overflow...
PQ.star <- P.star*Q.star
P.star.a1 <- PQ.star*B.mat
P.star.beta <- list()
for (i in seq_len(ncat-1))
P.star.beta[[i]] <- matrix(0, nr, nc)
for (i in seq_len(ncat-1))
P.star.beta[[i]][,i] <- (PQ.star*a1)[,i]
P.star.c <- PQ.star*a1
P.a1 <- matrix(NA, nr, ncat);
P.c <- matrix(NA, nr, ncat);
P.beta <- list(); for (i in 1:(ncat-1)) P.beta[[i]] <- matrix(NA, nr, ncat);
P.a1[, 2:nc] <- inv.P*(P.star.a1[,-nc] - P.star.a1[, -1])
for (i in seq_len(ncat-1))
P.beta[[i]][, 2:nc] <- inv.P*(P.star.beta[[i]][, -nc] - P.star.beta[[i]][, -1])
P.c[, 2:nc] <- inv.P * (P.star.c[, -nc] - P.star.c[, -1])
# the first response category
P.a1[,1] <- -P.star[,1]*B.mat[,1]
for (i in 2:(ncat-1))
P.beta[[i]][,1] <- 0
P.beta[[1]][,1] <- -P.star[,1]*a1
P.c[,1] <- -P.star[,1]*a1
# the last response category
P.a1[,ncat]=Q.star[,ncat-1]*B.mat[,ncat-1]
for (i in seq_len(ncat-2))
P.beta[[i]][,ncat] <- 0
P.beta[[ncat-1]][,ncat] <- Q.star[,ncat-1]*a1
P.c[,ncat] <- Q.star[,ncat-1]*a1
grad[1] <- sum(P.a1 * dat)
for (i in seq_len(ncat-1))
grad[i+1] <- sum(P.beta[[i]] *dat)
grad[ncat+1] <- sum(P.c * dat)
grad[!x@est] <- 0
Nest <- is.nan(grad) | is.infinite(grad)
if (sum(Nest)>0) {
grad <- rep(0, length(x@par))
grad[x@est & Nest] <- numerical_deriv(x@par[x@est & Nest], EML, obj=x,
Theta=Theta)
}
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'grsmIRT', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'grsmIRT', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# defines how the item should be initiallized with the S4 function new()
# before the parameter estimates are updated
# (set starting values, lower/upper bounds, logicals indicating estimation, etc)
setMethod("initialize",
'grsmIRT',
function(.Object, nfact, ncat){
if(nfact != 1L)
stop('grsmIRT only possible for unidimensional models')
stopifnot(ncat >= 2L)
.Object@par <- c(rep(1, nfact), seq(1, -1, length.out=ncat-1), 0)
#.Object@par <- c(rep(1, nfact), seq(-3, 3, length.out=ncat-1), 0)
# -3 ~ 3 seems to be too far away
names(.Object@par) = c(paste("a", seq_len(nfact), sep=""),
paste("b", seq_len(ncat-1L), sep=""), "c")
.Object@est <- c(rep(TRUE, nfact), rep(TRUE,ncat-1L), TRUE)
.Object@parnames <- names(.Object@par)
.Object@lbound <- rep(-Inf, nfact+ncat)
.Object@ubound <- rep(Inf, nfact+ncat)
.Object
}
)
# ----------------------------------------------------------------
setClass("gpcmIRT", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'gpcmIRT'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'gpcmIRT'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'gpcmIRT'),
definition = function(x){
x@par[1L]
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'gpcmIRT'),
definition = function(x){
x@par[-c(1, length(x@par))]
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'gpcmIRT'),
definition = function(x){
par <- c(rlnorm(1,0,1),
sort(rnorm(x@ncat-1, 0, 1), decreasing=FALSE), 0)
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'gpcmIRT'),
definition = function(x){
x@par[1L] <- 0
x@est[1L] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'gpcmIRT', Theta = 'matrix'),
definition = function(x, Theta){
.Call('gpcmIRTTraceLinePts', x@par, Theta, FALSE, numeric(0))
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'gpcmIRT', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
ret <- .Call("dparsgpcmIRT", x@par, Theta, offterm, x@dat,
length(x@par) - ncol(Theta), FALSE)
hess <- matrix(0, length(x@par), length(x@par))
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
ret <- list(grad=ret$grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
# derivative of the model wrt to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'gpcmIRT', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta)
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'gpcmIRT', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("monopoly", contains = 'AllItemsClass',
representation = representation(k='integer'))
setMethod(
f = "print",
signature = signature(x = 'monopoly'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'monopoly'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'monopoly'),
definition = function(x){
x@par[1L]
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'monopoly'),
definition = function(x){
x@par[-c(1, length(x@par))]
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'monopoly'),
definition = function(x){
par <- c(rnorm(1),
sort(rnorm(x@ncat-1, 0, 2)),
rnorm(x@k*2))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'monopoly'),
definition = function(x){
x@par[1L] <- 0
x@est[1L] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'monopoly', Theta = 'matrix'),
definition = function(x, Theta){
.Call('monopolyTraceLinePts', x@par, Theta, x@ncat, x@k)
}
)
setMethod(
f = "Deriv",
signature = signature(x = 'monopoly', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
if(nrow(x@fixed.design) > 1L && ncol(x@fixed.design) > 0L)
Theta <- cbind(x@fixed.design, Theta)
grad <- numeric(length(x@par))
grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta)
hess <- matrix(0, length(x@par), length(x@par))
if(estHess && any(x@est))
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
return(ret)
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'monopoly', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'monopoly', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("sequential", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'sequential'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'sequential'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'sequential'),
definition = function(x){
x@par[seq_len(x@nfact)] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'sequential'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'sequential'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
sort(rnorm(x@ncat-1L, sd = 2), decreasing=TRUE))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'sequential'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'sequential', Theta = 'matrix'),
definition = function(x, Theta, itemexp = FALSE){
par <- x@par
a <- par[1L:ncol(Theta)]
d <- par[-c(1L:ncol(Theta))]
ncat <- length(d) + 1L
Z <- matrix(0, nrow(Theta), ncat - 1L)
for(i in seq_along(d))
Z[,i] <- Theta %*% a + d[i]
P <- matrix(0, nrow(Theta), ncat)
Psi <- matrix(0, nrow(Theta), ncat + 1L)
Psi[,1L] <- 1
for(i in seq_along(d))
Psi[,i+1L] <- 1 / (1 + exp(-Z[,i]))
Pkc <- t(apply(Psi, 1L, cumprod))
for(i in 1L:ncat)
P[,i] <- Pkc[,i] * (1 - Psi[,i+1L])
if(itemexp) return(list(Psi=Psi, P=P))
P
}
)
# complete-data derivative used in parameter estimation (here it is done numerically)
setMethod(
f = "Deriv",
signature = signature(x = 'sequential', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
dplist <- ProbTrace(x, Theta, itemexp=TRUE)
P <- dplist$P
Psi <- dplist$Psi
Psi2 <- Psi^2
ncat <- x@ncat
nfact <- ncol(Theta)
Pgrad <- vector("list", length(x@par))
storage <- matrix(0, nrow(Theta), ncat)
# as
ind <- 1L
for(fact in seq_len(nfact)){
for(cat in seq_len(ncat)){
collecta <- (cat - 1) * Theta[,fact] * P[,cat]
tmp <- Psi[,1L:cat, drop=FALSE]
prodtmp <- apply(tmp, 1L, prod)
if(cat > 1L){
for(tcat in seq_len(cat-1L)){
collecta <- collecta -
Theta[,fact] * (Psi[,tcat+1L] * P[,cat])
}
}
collecta <- collecta + prodtmp *
(Theta[,fact] * (Psi2[,cat + 1L] - Psi[,cat + 1L]))
storage[,cat] <- collecta
}
Pgrad[[ind]] <- storage
ind <- ind + 1L
}
#ds
for(cat in seq_len(ncat-1L)){
tmp <- Psi[,1L:cat, drop=FALSE]
prodtmp <- apply(tmp, 1L, prod)
for(tcat in 1L:ncat){
if(tcat < cat){
storage[,tcat] <- 0
} else if(tcat == cat){
storage[,tcat] <- prodtmp * (Psi2[,cat + 1L] - Psi[,cat + 1L])
} else {
storage[,tcat] <- P[,tcat] - Psi[,cat + 1L] * P[,tcat]
}
}
Pgrad[[ind]] <- storage
ind <- ind + 1L
}
Pgrad <- array(do.call(c, Pgrad), c(nrow(Theta), x@ncat, length(x@par)))
grad <- symbolicGrad_par(x, Theta, dp1=Pgrad, P=P)
hess <- matrix(0, length(x@par), length(x@par))
if(estHess && any(x@est)){
# grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta)
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
}
return(list(grad=grad, hess=hess)) #TODO replace with analytic derivatives
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'sequential', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'sequential', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("fivePL", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'fivePL'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'fivePL'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'fivePL'),
definition = function(x){
x@par[seq_len(x@nfact)] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'fivePL'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'fivePL'),
definition = function(x){
par <- c(rlnorm(x@nfact, meanlog=0, sdlog=.5),
rnorm(1L),
rnorm(1L, -2, .5),
rnorm(1L, 2, .5),
1)
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'fivePL'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
P.5PL <- function(x, Theta, ncat = 2L){
a <- x@par[1L]
d <- x@par[2L]
g <- 1 / (1 + exp(-x@par[1L + 2L]))
u <- 1 / (1 + exp(-x@par[1L + 3L]))
S <- exp(x@par[5L])
P <- g + (u-g)/ (1 + exp(-1*(a * Theta + d)) )^S
cbind(1-P, P)
}
dP.5PL <- function(x, Theta, ncat = 2L)
{
.e1 <- x@par
.e3 <- exp(-(Theta * .e1[1L] + .e1[2L]))
.e4 <- exp(.e1[5L])
.e5 <- 1 + .e3
.e8 <- .e5^.e4
.e13 <- 1/(1 + exp(-.e1[4L])) - 1/(1 + exp(-.e1[3L]))
.e14 <- .e5^(1 + .e4)
.e15 <- 1/.e8
.e18 <- .e13 * .e3 * .e4/.e14
.e20 <- .e13 * log1p(.e3)/.e8
.e21 <- 1 - .e15
.e25 <- Theta * .e13 * .e3 * .e4/.e14
c(a = cbind(-.e25, .e25),
d = cbind(-.e18, .e18), g = cbind(-.e21, .e21),
u = cbind(-.e15, .e15), S = cbind(.e20, -.e20))
}
d2P.5PL <- function(x, Theta, ncat = 2L)
{
.e1 <- x@par
.e3 <- exp(-(Theta * .e1[1L] + .e1[2L]))
.e4 <- exp(.e1[5L])
.e5 <- 1 + .e4
.e6 <- 1 + .e3
.e7 <- .e6^.e5
.e14 <- 1/(1 + exp(-.e1[4L])) - 1/(1 + exp(-.e1[3L]))
.e15 <- log1p(.e3)
.e16 <- 1/.e7
.e17 <- 2 * .e5
.e18 <- .e6^.e4
.e22 <- .e6^(.e4 - .e17) * .e5 * .e3 - .e16
.e24 <- .e3 * .e4/.e7
.e25 <- .e15/.e18
.e28 <- Theta * .e3 * .e4/.e7
.e29 <- -.e24
.e30 <- -.e25
.e31 <- -.e28
.e35 <- .e16 - .e6^(.e5 - .e17) * .e4 * .e15
.e38 <- .e4 * .e15/.e7 - .e16
.e42 <- Theta * .e22 * .e14 * .e3 * .e4
.e46 <- .e22 * .e14 * .e3 * .e4
.e48 <- .e35 * .e14 * .e3
.e50 <- .e14 * .e3 * .e38
.e52 <- .e14 * .e15^2/.e18
.e53 <- cbind(.e29, .e24)
.e54 <- cbind(.e30, .e25)
.e55 <- cbind(-.e42, .e42)
.e56 <- cbind(.e31, .e28)
.e57 <- cbind(.e24, .e29)
.e58 <- cbind(.e25, .e30)
.e59 <- cbind(.e28, .e31)
.e62 <- Theta * .e35 * .e14 * .e3
.e65 <- Theta * .e14 * .e3 * .e38
.e70 <- Theta^2 * .e22 * .e14 * .e3 * .e4
c(a = c(a = cbind(-.e70, .e70), d = .e55, g = .e59, u = .e56,
S = cbind(.e65, -.e65)),
d = c(a = .e55, d = cbind(-.e46, .e46), g = .e57, u = .e53,
S = cbind(.e50, -.e50)),
g = c(a = .e59, d = .e57, g = c(0, 0), u = c(0, 0), S = .e54),
u = c(a = .e56, d = .e53, g = c(0, 0), u = c(0, 0), S = .e58),
S = c(a = cbind(-.e62, .e62), d = cbind(-.e48, .e48), g = .e54, u = .e58,
S = cbind(-.e52, .e52)))
}
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'fivePL', Theta = 'matrix'),
definition = function(x, Theta, itemexp = FALSE){
ret <- P.5PL(x=x, Theta=Theta)
ret
}
)
# complete-data derivative used in parameter estimation
setMethod(
f = "Deriv",
signature = signature(x = 'fivePL', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(any(x@est)){
dp1 <- array(dP.5PL(x, Theta), c(nrow(Theta),x@ncat,length(x@par)))
grad <- symbolicGrad_par(x, Theta, dp1=dp1)
if(estHess){
dp2 <- array(d2P.5PL(x, Theta), c(nrow(Theta),x@ncat,length(x@par), length(x@par)))
hess <- symbolicHessian_par(x, Theta, dp1=dp1, dp2=dp2)
}
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
ret
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'fivePL', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'fivePL', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("cll", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'cll'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'cll'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'cll'),
definition = function(x){
NA
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'cll'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'cll'),
definition = function(x){
par <- rnorm(1L)
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'cll'),
definition = function(x){
x
}
)
P.cll <- function(x, Theta, ncat = 2L){
b <- x@par[1L]
P <- 1 - exp(-1 * exp(Theta - b))
P <- ifelse(P > 1 - 1e-16, 1 - 1e-16, P)
cbind(1-P, P)
}
dP.cll <- function(x, Theta, ncat = 2L)
{
.e1 <- exp(Theta - x@par[1L])
.e4 <- exp(-.e1) * .e1
cbind(.e4, -.e4)
}
d2P.cll <- function(x, Theta, ncat = 2L)
{
.e1 <- exp(Theta - x@par[1L])
.e5 <- exp(-.e1) * .e1 * (.e1 - 1)
cbind(.e5, -.e5)
}
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'cll', Theta = 'matrix'),
definition = function(x, Theta, itemexp = FALSE){
ret <- P.cll(x=x, Theta=Theta)
ret
}
)
# complete-data derivative used in parameter estimation
setMethod(
f = "Deriv",
signature = signature(x = 'cll', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(any(x@est)){
dp1 <- array(dP.cll(x, Theta), c(nrow(Theta),x@ncat,length(x@par)))
grad <- symbolicGrad_par(x, Theta, dp1=dp1)
if(estHess){
dp2 <- array(d2P.cll(x, Theta), c(nrow(Theta),x@ncat,length(x@par), length(x@par)))
hess <- symbolicHessian_par(x, Theta, dp1=dp1, dp2=dp2)
}
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
ret
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'cll', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'cll', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
setClass("ull", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'ull'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'ull'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'ull'),
definition = function(x){
NA
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'ull'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'ull'),
definition = function(x){
par <- c(rlnorm(1L, .2, .3), sort(rnorm(x@ncat-1L), decreasing = TRUE))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'ull'),
definition = function(x){
x@par[1L] <- 0
x@est[1L] <- FALSE
x
}
)
P.ull <- function(x, Theta, ncat){
eta <- x@par[1L]
lambdas <- exp(x@par[2L:length(x@par)])
TheEta <- Theta^eta
ret <- matrix(0, nrow(Theta), ncat)
PS <- matrix(0, nrow(Theta), ncat+1L)
PS[,1L] <- 1
for(i in 1L:(ncat-1L))
PS[,i+1L] <- lambdas[i] * TheEta / (1 + lambdas[i] * TheEta)
for(i in 1L:ncat)
ret[,i] <- PS[,i] - PS[,i+1L]
ret
}
dP.ull <- function(x, Theta, ncat)
{
N <- nrow(Theta)
eta <- x@par[1L]
lambdas <- exp(x@par[2L:length(x@par)])
TheEta <- Theta^eta
Theta2Eta <- Theta^(2 * eta)
dlambdas <- vector('list', ncat-1L)
for(j in 2L:ncat - 1L){
mat <- matrix(0, N, ncat)
tmp <- (TheEta - Theta2Eta * lambdas[j] /
(1 + TheEta * lambdas[j])) / (1 + TheEta * lambdas[j])
mat[,j] <- -tmp
mat[,j+1L] <- tmp
dlambdas[[j]] <- mat
}
dlambdas <- do.call(cbind, dlambdas)
.e9 <- log(Theta)
.e10_e14 <- 1 + as.vector(TheEta) * matrix(lambdas, N, ncat-1L, byrow=TRUE)
mat <- matrix(0, N, ncat)
mat[,1L] <- -(.e9 * lambdas[1L] *
(TheEta - Theta2Eta * lambdas[1L] / .e10_e14[,1L])/.e10_e14[,1L])
.e25 <- TheEta - Theta2Eta * lambdas[1L] / .e10_e14[,1L]
.e26 <- TheEta - Theta2Eta * lambdas[ncat-1L] / .e10_e14[,ncat-1L]
mat[,ncat] <- lambdas[ncat-1L] * .e9 * .e26 / .e10_e14[,ncat-1L]
if(ncat > 2L){
for(j in 2:(ncat - 2L)){
tmp <- lambdas[j] * (TheEta - Theta2Eta * lambdas[j] /
.e10_e14[,j])/.e10_e14[,j]
if(j == 2L)
mat[,2L] <- (lambdas[1L] * .e25 / .e10_e14[,1L] - tmp) * .e9
if(ncat == 3L) break
tmp2 <- lambdas[j+1L] * (TheEta - Theta2Eta * lambdas[j+1L] /
.e10_e14[,j+1L])/.e10_e14[,j+1L]
mat[,j+1L] <- (tmp - tmp2) * .e9
}
}
c(cbind(mat, dlambdas))
}
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'ull', Theta = 'matrix'),
definition = function(x, Theta, itemexp = FALSE){
ret <- P.ull(x=x, Theta=Theta, ncat=x@ncat)
ret
}
)
# complete-data derivative used in parameter estimation
setMethod(
f = "Deriv",
signature = signature(x = 'ull', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(any(x@est)){
dp1 <- array(dP.ull(x, Theta, ncat=x@ncat), c(nrow(Theta),x@ncat,length(x@par)))
grad <- symbolicGrad_par(x, Theta, dp1=dp1)
if(estHess){
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,Theta=Theta,
gradient=FALSE)
}
}
ret <- list(grad=grad, hess=hess)
if(x@any.prior) ret <- DerivativePriors(x=x, grad=ret$grad, hess=ret$hess)
ret
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'ull', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'ull', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# ----------------------------------------------------------------
# experimental itemtype (used as a template to create custom IRT models). Edit this if you want
# to experiment with your own customized IRT models
setClass("experimental", contains = 'AllItemsClass',
representation = representation())
setMethod(
f = "print",
signature = signature(x = 'experimental'),
definition = function(x, ...){
cat('Item object of class:', class(x))
}
)
setMethod(
f = "show",
signature = signature(object = 'experimental'),
definition = function(object){
print(object)
}
)
#extract the slopes (should be a vector of length nfact)
setMethod(
f = "ExtractLambdas",
signature = signature(x = 'experimental'),
definition = function(x){
x@par[seq_len(x@nfact)] #slopes
}
)
#extract the intercepts
setMethod(
f = "ExtractZetas",
signature = signature(x = 'experimental'),
definition = function(x){
x@par[length(x@par)] #intercepts
}
)
# generating random starting values (only called when, e.g., mirt(..., GenRandomPars = TRUE))
setMethod(
f = "GenRandomPars",
signature = signature(x = 'experimental'),
definition = function(x){
par <- c(rlnorm(1, .2, .2), rnorm(1))
x@par[x@est] <- par[x@est]
x
}
)
# how to set the null model to compute statistics like CFI and TLI (usually just fixing slopes to 0)
setMethod(
f = "set_null_model",
signature = signature(x = 'experimental'),
definition = function(x){
x@par[seq_len(x@nfact)] <- 0
x@est[seq_len(x@nfact)] <- FALSE
x
}
)
# probability trace line function. Must return a matrix with a trace line for each category
setMethod(
f = "ProbTrace",
signature = signature(x = 'experimental', Theta = 'matrix'),
definition = function(x, Theta){
a <- x@par[1L]
b <- x@par[2L]
p <- exp(a * (Theta - b)) / (1 + exp(a * (Theta - b)))
p <- ifelse(p < 1e-10, 1e-10, p) #numerical constraints to avoid log() problems
p <- ifelse(p > 1 - 1e-10, 1 - 1e-10, p)
P <- cbind(1-p, p)
return(P)
}
)
# complete-data derivative used in parameter estimation (here it is done numerically)
setMethod(
f = "Deriv",
signature = signature(x = 'experimental', Theta = 'matrix'),
definition = function(x, Theta, estHess = FALSE, offterm = numeric(1L)){
grad <- rep(0, length(x@par))
hess <- matrix(0, length(x@par), length(x@par))
if(any(x@est)){
grad[x@est] <- numerical_deriv(x@par[x@est], EML, obj=x, Theta=Theta)
if(estHess){
hess[x@est, x@est] <- numerical_deriv(x@par[x@est], EML, obj=x,
Theta=Theta, gradient=FALSE)
}
}
return(list(grad=grad, hess=hess)) #TODO replace with analytic derivatives
}
)
# derivative of the model wft to the Theta values (done numerically here)
setMethod(
f = "DerivTheta",
signature = signature(x = 'experimental', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_DerivTheta(x, Theta) #replace with analytical derivatives
}
)
# derivative of the probability trace line function wrt Theta (done numerically here)
setMethod(
f = "dP",
signature = signature(x = 'experimental', Theta = 'matrix'),
definition = function(x, Theta){
numDeriv_dP(x, Theta) #replace with analytical derivatives
}
)
# defines how the item should be initiallized with the S4 function new()
# before the parameter estimates are updated
# (set starting values, lower/upper bounds, logicals indicating estimation, etc)
setMethod("initialize",
'experimental',
function(.Object, nfact, ncat){
stopifnot(nfact == 1L)
stopifnot(ncat == 2L)
.Object@par <- c(a=1, b=0)
.Object@est <- c(TRUE, TRUE)
.Object@parnames <- names(.Object@par)
.Object@lbound <- rep(-Inf, 2L)
.Object@ubound <- rep(Inf, 2L)
.Object
}
)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.