## File Name: tam_mml_2pl_mstep_slope.R
## File Version: 9.579
tam_mml_2pl_mstep_slope <- function (B_orig, B, B_obs, B.fixed, max.increment, nitems, A,
AXsi, xsi, theta, nnodes, maxK, itemwt, Msteps, ndim, convM,
irtmodel, progress, est.slopegroups, E, basispar, se.B,
equal.categ, B_acceleration, trim_increment="half", iter, eps=1E-10,
maxcat=NULL, use_rcpp=FALSE, use_rcpp_calc_prob=TRUE )
{
old_increment <- array( max.increment, dim=c(nitems, maxK, ndim) )
xbar2 <- xxf <- xbar <- array(0, dim=c( nitems, maxK ) )
xtemp <- matrix(0, nrow=1, ncol=1)
converge <- FALSE
Biter <- 1
mK <- 1:maxK
items.temp <- 1:nitems
items.conv <- NULL
items.conv1 <- c(-1)
increments_msteps <- rep(NA,Msteps)
B_old <- B
eps <- 1E-10
NI <- nitems
TP <- nrow(theta)
if (progress){
cat("\nM Step Slopes |")
utils::flush.console()
}
if (irtmodel %in% c("GPCM","GPCM.groups") ){
old_increment.temp <- matrix( .3, nrow=NI, ncol=ndim )
}
if (irtmodel=="GPCM.design" ){
Nlambda <- ncol(E) # number of lambda parameters
old_increment.temp <- matrix( .3, nrow=Nlambda, ncol=ndim )
}
if (irtmodel=="2PL.groups"){
ES <- length( unique( est.slopegroups) )
old_increment.temp <- array( .3, dim=c(ES, maxK,ndim ) )
use_rcpp <- FALSE
}
#------------------------------------------------
# begin algorithm M-steps
while (!converge & ( Biter <=Msteps ) ) {
#compute expectation
res <- tam_mml_calc_prob( iIndex=items.temp, A=A, AXsi=AXsi, B=B, xsi=xsi,
theta=theta, nnodes=nnodes, maxK=maxK, maxcat=maxcat,
use_rcpp=use_rcpp_calc_prob )
rprobs <- res$rprobs
LIT <- length(items.temp)
###### D I M E N S I O N S ######
for (dd in 1:ndim){
if ( irtmodel %in% c("GPCM","GPCM.design","GPCM.groups") ){
xtemp <- matrix(0, nrow=LIT, ncol=TP )
}
if ( irtmodel=="2PL.groups"){
xtemp <- array(0, dim=c( LIT, TP, maxK) )
}
##### C A T E G O R I E S #########
if ( ! equal.categ ){
rprobs[ is.na(rprobs) ] <- 0
}
if (use_rcpp){
res <- tam_rcpp_mml_2pl_mstep_item_slopes_suffstat(
rprobs=as.vector(rprobs), items_temp=items.temp, theta=theta, dd=dd-1, LIT=LIT,
TP=TP, nitems=nitems, maxcat=maxcat, maxK=maxK, itemwt=itemwt, xxf_=xxf,
xbar_=xbar, xbar2_=xbar2, irtmodel=irtmodel, xtemp_=xtemp,
items_conv=items.conv1 )
} else {
res <- tam_mml_2pl_mstep_item_slopes_suffstat_R( rprobs=rprobs, maxK=maxK,
LIT=LIT, TP=TP, itemwt=itemwt, theta=theta, dd=dd, items.temp=items.temp,
items.conv=items.conv, xbar=xbar, xbar2=xbar2, xxf=xxf, xtemp=xtemp,
irtmodel=irtmodel )
}
xbar <- res$xbar
xbar2 <- res$xbar2
xxf <- res$xxf
xtemp <- res$xtemp
#----------------
if ( irtmodel %in% c("GPCM","GPCM.design", "GPCM.groups")){ # begin GPCM / GPCM.design
B_obs.temp <- B_obs
B_obs.temp[ items.temp,,dd] <- tam_matrix2( mK-1, LIT, maxK ) * B_obs[ items.temp,,dd]
xbar[ is.na(xbar) ] <- 0
xbar.temp <- tam_matrix2( mK-1, LIT, maxK ) *xbar
diff.temp <- rowSums(B_obs.temp[,,dd]) - rowSums(xbar.temp)
xbar2.temp <- diag( xtemp^2 %*% itemwt[, items.temp ] )
xxf.temp <- rowSums( xxf )
if ( ! is.null(items.conv) ){
diff.temp[ items.conv, ] <- 0
}
deriv.temp <- xbar2.temp - xxf.temp
} # loop GPCM and GPCM.design
#----------------
if (irtmodel %in% c("GPCM","GPCM.groups")){ # begin GPCM
deriv.temp[ is.na(deriv.temp)] <- eps
if (irtmodel=="GPCM.groups"){
diff.temp <- tam_aggregate_derivative_information(deriv=diff.temp,
groups=est.slopegroups)
deriv.temp <- tam_aggregate_derivative_information(deriv=deriv.temp,
groups=est.slopegroups)
}
increment.temp <- diff.temp*abs(1/( deriv.temp + eps ) )
increment.temp <- tam_trim_increment(increment=increment.temp,
max.increment=abs(old_increment.temp[,dd]),
trim_increment=trim_increment)
old_increment.temp[,dd] <- increment.temp
increment <- tam_outer( increment.temp, mK - 1)
if (Biter==1){
se.B[,,dd] <- tam_outer( sqrt( 1 / abs( deriv.temp )), mK-1 )
}
} # end GPCM
#----------------
if ( irtmodel=="GPCM.design"){
diff.temp <- ( t(E) %*% diff.temp )[,1]
deriv.temp <- ( t(E^2) %*% deriv.temp )[,1]
increment.temp <- diff.temp*abs(1/( deriv.temp + eps ) )
increment.temp <- tam_trim_increment(increment=increment.temp,
max.increment=abs(old_increment.temp[,dd]),
trim_increment=trim_increment)
old_increment.temp[,dd] <- increment.temp
basispar[,dd] <- basispar[,dd] + increment.temp
increment.temp <- ( E %*% increment.temp )[,1]
increment <- tam_outer( increment.temp, mK-1)
d1 <- tam_outer( 1 / abs( deriv.temp ), mK-1 )
LL <- ncol(d1)
for (ll in 1:LL){
m1 <- sqrt( diag( E %*% d1[,ll] %*% t( d1[,ll] ) %*% t(E) ) )
if (Biter==1){
se.B[,ll,dd] <- m1
}
}
nB <- dim(B)
B_ind <- 1 * ( B_orig !=0 )
for (dd in 1:nB[3]){
EB <- E %*% basispar[,dd]
for (cc in 1:(nB[2]-1)){
B[,cc+1,dd] <- cc * EB * B_ind[,cc+1,dd]
}
}
B00 <- B
} # end GPCM.design
#----------------
if (irtmodel=="2PL.groups"){ # begin 2PL slopegroups
a1 <- stats::aggregate( B_obs[,,dd] - xbar, list(est.slopegroups), sum )
a2 <- stats::aggregate( xbar2 - xxf, list(est.slopegroups), sum )
deriv.temp <- as.matrix(a2[,-1])
diff.temp <- as.matrix(a1[,-1])
increment.temp <- diff.temp*abs(1/( deriv.temp + eps ) )
increment.temp <- tam_trim_increment(increment=increment.temp,
max.increment=abs(old_increment.temp[,,dd]),
trim_increment=trim_increment)
old_increment.temp[,,dd] <- increment.temp
ind.match <- match( est.slopegroups, a1[,1] )
increment <- increment.temp[ ind.match, ]
if (Biter==1){ se.B[,,dd] <- sqrt( 1 / abs( deriv.temp[ind.match, ] )) }
} # end 2PL slope groups
#----------------
if (irtmodel=="2PL"){ # begin 2PL
diff <- B_obs[,,dd] - xbar
if ( ! is.null( items.conv) ){
diff[ items.conv, ] <- 0
}
deriv <- xbar2 - xxf
increment <- diff*abs(1/( deriv + eps ) )
increment <- tam_trim_increment(increment=increment,
max.increment=abs(old_increment[,,dd]),
trim_increment=trim_increment)
if (Biter==1){
se.B[,,dd] <- sqrt( 1 / abs( deriv ))
}
} # end 2PL
#----------------
increment[B_orig[,,dd]==0] <- 0 # B[i,k,] could be zero for some dimensions
old_increment[,,dd] <- increment
B[,,dd] <- B[,,dd] + increment # update B parameter
if (irtmodel=="GPCM.design"){
B <- B00
}
} # end dimensions dd
#----------------------------
if ( irtmodel=="2PL" ){
items.temp <- which( apply( old_increment, 1,
FUN=function(ll){ ! ( max(abs(ll)) < convM ) } ) )
items.conv <- setdiff( 1:nitems, items.temp )
if ( length(items.conv) > 0 ){
items.conv1 <- items.conv
}
}
max_change <- max(abs(old_increment))
increments_msteps[Biter] <- max_change
if ( max_change < convM ) {
converge <- TRUE
}
Biter <- Biter + 1
if (progress){
cat( "-" ) ; utils::flush.console()
}
} # end while loop
#**** end algorithm M-steps
#---------------------------------------------
se.B[ B_orig==0 ] <- 0
# acceleration
if ( B_acceleration$acceleration !="none" ){
B_acceleration <- tam_accelerate_parameters( xsi_acceleration=B_acceleration,
xsi=as.vector(B), iter=iter, itermin=3)
B <- array( B_acceleration$parm, dim(B) )
}
# change in B parameters
B_change <- max( abs( B- B_old) )
#---- OUTPUT
res <- list( "B"=B, "basispar"=basispar, "se.B"=se.B, Biter=Biter,
B_acceleration=B_acceleration, B_change=B_change,
increments_msteps=increments_msteps)
return(res)
}
#----------------------------------------------------------------------
Mstep_slope.v2 <- tam_mml_2pl_mstep_slope
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.