Nothing
Mstep <- function(delta,item.parm,J,Kj,Lj,itr,Rg,Ng,model,ConstrType,linkfunc,
correction,lower.p,upper.p,algorithm,ConstrMatrix,solver, item.prior,
DesignMatrices,MstepMessage, est.bin,auglag_args,nloptr_args,solnp_args){
# print(solver)
if(!is.numeric(model))
model <- match(model, c("LOGGDINA","LOGITGDINA","UDF", "GDINA", "DINA", "DINO", "ACDM", "LLM", "RRUM", "MSDINA")) - 4
opts <- list()
for(j in which(est.bin)){ #for each item
optims <- GNLOptim_call(par = delta[[j]],solver=solver[j],modelj=model[j],correction=correction,
auglag_args=auglag_args,nloptr_args=nloptr_args,solnp_args = solnp_args,item.prior=item.prior,
Nj=Ng[j,1:Lj[j]],Rj=Rg[j,1:Lj[j]],designMj=DesignMatrices[[j]],
uPj=upper.p[j],lPj=lower.p[j],MstepMessage=MstepMessage,itr=itr,j=j,
ConstrMatrix=ConstrMatrix[[j]],linkfunc=linkfunc[j],eps=1e-16,ConstrType=ConstrType[j])
if(is.null(optims)) stop(paste("M-step fails for item",j,"at iteration",itr),call. = FALSE)
delta[[j]] <- c(optims$delta)
item.parm[j,seq_len(length(optims$phat))] <- optims$phat
opts[[j]] <- optims
}
return(list(item.parm=item.parm,delta=delta,optims=opts))
}
GNLOptim <- function(par,fn,gr=NULL,hin=NULL,hin.gr=NULL,optimizer=NULL,
modelj,correction,Nj,Rj,designMj,uPj,lPj,itr,j,auglag_args,nloptr_args,solnp_args,
ConstrMatrix,linkfunc,eps,ConstrType,MstepMessage,...){
ret <- list()
if(optimizer=="ClosedForm"){
if(linkfunc!=1) stop("Closed-form solutions are only available for the identity link GDINA, DINA and DINO models.",call. = FALSE)
if(modelj==0){
if (any(Nj<correction[2])){
Nj[which(Nj<correction[2])] <- Nj[which(Nj<correction[2])] + correction[2]
Rj[which(Nj<correction[2])] <- Rj[which(Nj<correction[2])] + correction[1]
}
phat <- c(Rj/Nj)
}else if(modelj==1 | modelj==2){
rNj <- c(rep(sum(Nj[designMj[,2]==0],na.rm = TRUE),sum(designMj[,2]==0)),
rep(sum(Nj[designMj[,2]==1],na.rm = TRUE),sum(designMj[,2]==1)))
rRj <- c(rep(sum(Rj[designMj[,2]==0],na.rm = TRUE),sum(designMj[,2]==0)),
rep(sum(Rj[designMj[,2]==1],na.rm = TRUE),sum(designMj[,2]==1)))
if (any(rNj<correction[2])){
rNj[which(rNj<correction[2])] <- rNj[which(rNj<correction[2])] + correction[2]
rRj[which(rNj<correction[2])] <- rRj[which(rNj<correction[2])] + correction[1]
}
phat <- c(rRj/rNj)
}else{
stop("Closed-form solutions are only available for the identity link GDINA, DINA and DINO models.",call. = FALSE)
}
phat[phat<lPj] <- lPj
phat[phat>uPj] <- uPj
ret$delta <- c(Calc_Dj(phat, designMj, linkfunc = linkfunc, boundary = 1))
ret$opt <- list(convergence=0,message="converged - closed-form solution.")
ret$phat <- phat
ret$convergence <- 0
}else if(optimizer=="BFGS"){
if(modelj>=3&&modelj<=5){
par <- initials_optim(par,modelj=modelj,Nj=Nj,Rj=Rj,designMj=designMj,uPj=uPj,ConstrMatrix=ConstrMatrix,
lPj=lPj,linkfunc=linkfunc,ConstrType=ConstrType)
# print(par)
}
Lj <- nrow(designMj)
if(ConstrType==1){
ui <- rbind(designMj,-1*designMj)
cii <- rep(c(lPj,-1*uPj),each=Lj)
ci <- list(gdina=cii,dina=cii, dino=cii, acdm=cii,
llm=rep(c(qlogis(lPj),-1*qlogis(uPj)),each=Lj),
rrum=rep(c(log(lPj),-1*log(uPj)),each=Lj))
}else if(ConstrType==3){
if(modelj>0){
ui=rbind(designMj,-1*designMj,diag(length(par)))
}else if(modelj==0){
ui=ConstrMatrix%*%designMj
}
cii <- c(rep(c(lPj,-1*uPj),each=Lj),c(lPj,rep(0,length(par)-1)))
ci <- list(gdina=rep(0,nrow(ConstrMatrix)),dina=cii,dino=cii,acdm=cii,
llm=c(rep(c(qlogis(lPj),-1*qlogis(uPj)),each=Lj),c(qlogis(lPj),rep(0,length(par)-1))),
rrum=c(rep(c(log(lPj),-1*log(uPj)),each=Lj),c(log(lPj),rep(0,length(par)-1))))
}
op <- tryCatch(stats::constrOptim(theta = par,f = fn,grad = gr,method = "BFGS",
ui=ui,ci=ci[[modelj+1]],control=list(trace=0),
Nj=Nj,Rj=Rj,designMj=designMj,uPj=uPj,ConstrMatrix=ConstrMatrix,
lPj=lPj,linkfunc=linkfunc,eps=eps,ConstrType=ConstrType),
error=function(e){
if(MstepMessage){
message(paste("M-step optimization failed for item",j,"at iteration",itr))
message("Error message from stats::constrOptim:")
message(e)
}
return("error")
})
if(length(op)==1 && op=="error"){
ret$convergence <- -1
}else{
ret$delta <- op$par
ret$opt <- op
ret$phat <- c(Calc_Pj(op$par,designMj,linkfunc,FALSE))
if(any(ret$phat<0)||any(ret$phat>1)){
warning(paste("Estimates are out of bounds for item",j,"at iteration",itr,"[stats::constrOptim]"),call. = FALSE)
ret$convergence <- -3
}else{
ret$convergence <- 0
}
}
# print(ret$convergence)
}else if(optimizer=="auglag"){
args.list <- c(list(par=par,fn=fn,gr=gr,hin = hin,hin.jac = hin.gr,
Nj=Nj,Rj=Rj,designMj=designMj,uPj=uPj,ConstrMatrix=ConstrMatrix,
lPj=lPj,linkfunc=linkfunc,eps=eps,ConstrType=ConstrType),auglag_args)
op <- tryCatch(do.call(alabama::auglag,args.list),
error=function(e){
if(MstepMessage){
message(paste("M-step optimization failed for item",j,"at iteration",itr))
message("Error message from alabama::auglag:")
message(e)
}
return("error")
})
if(length(op)==1&&op=="error"){
ret$convergence <- -1
}else{
ret$delta <- op$par
ret$opt <- op
ret$phat <- c(Calc_Pj(op$par,designMj,linkfunc,FALSE))
if(any(ret$phat<0)||any(ret$phat>1)){
warning(paste("Estimates are out of bounds for item",j,"at iteration",itr,"[alabama::auglag]"),call. = FALSE)
ret$convergence <- -3
}else{
ret$convergence <- 0
}
}
}else if(optimizer=="solnp"){
if(ConstrType==1){nineq <- 2*nrow(designMj)}else if(ConstrType==3){nineq <- 2*nrow(designMj)+nrow(ConstrMatrix)}
args.list <- list(pars=par,fun=fn,ineqfun = hin,ineqLB = rep(0,nineq),
ineqUB = rep(Inf,nineq),control=solnp_args,
Nj=Nj,Rj=Rj,designMj=designMj,uPj=uPj,ConstrMatrix=ConstrMatrix,
lPj=lPj,linkfunc=linkfunc,eps=eps,ConstrType=ConstrType)
# print(args.list)
op <- tryCatch(do.call(Rsolnp::solnp,args.list),
error=function(e){
if(MstepMessage){
message(paste("M-step optimization failed for item",j,"at iteration",itr,"\n"))
message("Error message from Rsolnp::solnp:\n")
message(e)
}
return("error")
})
if(length(op)==1&&op=="error"){
ret$convergence <- -1
}else{
print(op$pars)
ret$delta <- op$pars
ret$opt <- op
ret$phat <- c(Calc_Pj(c(op$pars),designMj,linkfunc,FALSE))
if(any(ret$phat<0)||any(ret$phat>1)){
message(paste("Estimates are out of bounds for item",j,"at iteration",itr,"[Rsolnp:solnp]"),call. = FALSE)
ret$convergence <- -3
}else{
ret$convergence <- 0
}
}
}else if(optimizer=="slsqp"){
op <- tryCatch(nloptr::nloptr(x0=par,eval_f=fn,eval_grad_f=gr,eval_g_ineq = hin,eval_jac_g_ineq = hin.gr,
Nj=Nj,Rj=Rj,designMj=designMj,uPj=uPj,ConstrMatrix=ConstrMatrix,greaterthan0 = FALSE,
lPj=lPj,linkfunc=linkfunc,eps=eps,ConstrType=ConstrType,
opts=nloptr_args,...),
error=function(e){
if(MstepMessage){
message(paste("M-step optimization failed for item",j,"at iteration",itr))
message("Error message from nloptr::slsqp:")
message(e)
}
return("error")
})
if(length(op)==1&&op=="error"){
ret$convergence <- -1
}else{
ret$delta <- op$solution
ret$opt <- op
ret$phat <- c(Calc_Pj(op$solution,designMj,linkfunc,FALSE))
if(any(ret$phat<0)||any(ret$phat>1)){
message(paste("Estimates are out of bounds for item",j,"at iteration",itr,"[nloptr:slsqp]"),call. = FALSE)
ret$convergence <- -3
}else{
ret$convergence <- 0
}
}
}else if(optimizer=="nloptr"){
op <- tryCatch(nloptr::nloptr(x0=par,eval_f=fn,eval_grad_f=gr,eval_g_ineq = hin,eval_jac_g_ineq = hin.gr,
Nj=Nj,Rj=Rj,designMj=designMj,uPj=uPj,ConstrMatrix=ConstrMatrix,greaterthan0 = FALSE,
lPj=lPj,linkfunc=linkfunc,eps=eps,ConstrType=ConstrType,
opts=nloptr_args,...),
error=function(e){
if(MstepMessage){
message(paste("M-step optimization failed for item",j,"at iteration",itr))
message("Error message from nloptr::nloptr:")
message(e)
}
return("error")
})
if(length(op)==1&&op=="error"){
ret$convergence <- -1
}else{
ret$delta <- op$solution
ret$opt <- op
ret$phat <- c(Calc_Pj(op$solution,designMj,linkfunc,FALSE))
if(any(ret$phat<0)||any(ret$phat>1)){
message(paste("Estimates are out of bounds for item",j,"at iteration",itr,"[nloptr:nloptr]"),call. = FALSE)
ret$convergence <- -3
}else{
ret$convergence <- 0
}
}
}
ret
}
initials_optim <- function(par,modelj,Nj,Rj,designMj,uPj,lPj, ConstrMatrix,linkfunc,eps=1e-5,ConstrType){
# print(par)
np <- length(par)
pj <- c(Calc_Pj(par,designMj,linkfunc,FALSE))
if(ConstrType==1){
if(any(pj<=lPj)||any(pj>=uPj)){
# print(pj)
pj[which(pj<=lPj)] <- lPj + eps
pj[which(pj>=uPj)] <- uPj - eps
tmpar <- Calc_Dj(pj[seq_len(np)],
designMj[seq_len(np),],
linkfunc, FALSE)
tmpj <- c(Calc_Pj(tmpar,designMj,linkfunc,FALSE))
if(any(tmpj<=lPj)||any(tmpj>=uPj)){
if(modelj==4){
pj <- qlogis(pj)
}else if(modelj==5){
pj <- log(pj)
}
tmpar <- c(min(pj),rep((max(pj)-min(pj))/(np-1),np-1))
}
par <- tmpar
# print(c(Calc_Pj(par,designMj,linkfunc,boundary = 0)))
}
}else if(ConstrType==3){
if(modelj==3){#ACDM
newlPj <- lPj
newuPj <- uPj
newpj <- pj
}else if (modelj==4){
newlPj <- qlogis(lPj)
newuPj <- qlogis(uPj)
newpj <- qlogis(pj)
}else if (modelj==5){
newlPj <- log(lPj)
newuPj <- log(uPj)
newpj <- log(pj)
}
if(par[1]<=newlPj || any(par[2:np]<=0)){
par[which(par<=newlPj)] <- newlPj + eps
par[which(par[2:np]<=0)] <- eps
}
if (sum(par)>=newuPj) {
par <- par - par * (sum(par) - newuPj + eps) / sum(par)
}
if (par[1]<=newlPj || any(par[2:np]<=0) || sum(par)>=newuPj) {
Lpar <- which.max(c(newlPj + eps,min(newpj,na.rm = TRUE)))
Upar <- which.min(c(newuPj - eps,max(newpj,na.rm = TRUE)))
par <- c(Lpar,rep((Upar-Lpar)/(np-1),np-1))
}
}
return(par)
}
GNLOptim_call <- function(par,solver,modelj,correction,auglag_args,nloptr_args,solnp_args,item.prior=item.prior,
Nj,Rj,designMj,uPj,lPj,MstepMessage,itr,j,ConstrMatrix,linkfunc,eps,ConstrType){
if(item.prior$on){ # adding prior to item parameters
ret <- list()
if(modelj==0){
phat <- c((Rj+item.prior$beta[1]-1)/(Nj+sum(item.prior$beta)-2))
phat[phat<lPj] <- lPj
phat[phat>uPj] <- uPj
ret$delta <- c(Calc_Dj(phat, designMj, linkfunc, TRUE))
ret$opt <- list(convergence=0,message="converged - closed-form solution.")
ret$phat <- phat
ret$convergence <- 0
return(ret)
}else if(modelj==1 | modelj==2){
rNj <- c(rep(sum(Nj[designMj[,2]==0]),sum(designMj[,2]==0)),
rep(sum(Nj[designMj[,2]==1]),sum(designMj[,2]==1)))
rRj <- c(rep(sum(Rj[designMj[,2]==0]),sum(designMj[,2]==0)),
rep(sum(Rj[designMj[,2]==1]),sum(designMj[,2]==1)))
phat <- c((rRj+item.prior$beta[1]-1)/(rNj+sum(item.prior$beta)-2))
phat[phat<lPj] <- lPj
phat[phat>uPj] <- uPj
ret$delta <- c(Calc_Dj(phat, designMj, linkfunc, TRUE))
ret$opt <- list(convergence=0,message="converged - closed-form solution.")
ret$phat <- phat
ret$convergence <- 0
return(ret)
}else if(linkfunc==2){ #logit link function
op <- tryCatch(nloptr::slsqp(x0=par,fn=Mstep_obj_fn_prior,
Nj=Nj,Rj=Rj,designMj=designMj,uPj=uPj,ConstrMatrix=ConstrMatrix,greaterthan0 = FALSE,
lPj=lPj,linkfunc=linkfunc,eps=eps,ConstrType=ConstrType,m=item.prior$normal[1],sd=item.prior$normal[2]),
error=function(e){
if(MstepMessage){
message(paste("M-step optimization failed for item",j,"at iteration",itr))
message("Error message from nloptr::slsqp:")
message(e)
}
return("error")
})
if(length(op)==1&&op=="error"){
ret$convergence <- -1
}else{
ret$delta <- op$par
ret$opt <- op
ret$phat <- c(Calc_Pj(op$par,designMj,linkfunc,FALSE))
if(any(ret$phat<0)||any(ret$phat>1)){
message(paste("Estimates are out of bounds for item",j,"at iteration",itr,"[nloptr:slsqp]"),call. = FALSE)
ret$convergence <- -3
}else{
ret$convergence <- 0
}
}
return(ret)
}else{
stop("Priors cannot be imposed.",call. = FALSE)
}
}else{ # No priors
if(solver=="auto"){
if(ConstrType==1){
if(modelj<3&modelj>=0){
solver <- "ClosedForm"
}else if(modelj>=3){
solver <- c("BFGS","slsqp","auglag","solnp")
}else if(modelj==-1||modelj==-2||modelj==-3){
solver <- c("slsqp","auglag","solnp")
}
}else{
if(modelj>0){
solver <- c("BFGS","auglag","solnp")
}else if(modelj<=0){
solver <- c("slsqp","auglag","solnp")
}
}
}
conv <- vector("numeric",length(solver))
for(s in solver){
if(s=="slsqp"){nloptr_args$algorithm = "NLOPT_LD_SLSQP"}
optims <- GNLOptim(par = par,fn=Mstep_obj_fn,gr=Mstep_obj_gr,hin=Mstep_ineq_fn,hin.gr=Mstep_ineq_jac,
optimizer=s,modelj=modelj,correction=correction,
auglag_args=auglag_args,nloptr_args=nloptr_args,solnp_args = solnp_args,
Nj=Nj,Rj=Rj,designMj=designMj,
uPj=uPj,lPj=lPj,MstepMessage=MstepMessage,itr=itr,j=j,
ConstrMatrix=ConstrMatrix,linkfunc=linkfunc,eps=1e-16,ConstrType=ConstrType)
conv[s] <- optims$convergence
if(optims$convergence==0) break
}
if(all(conv<0)){
return(NULL)
}else{
return(optims)
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.