#targetPath<-"\\\\libfs1\\users$\\home\\zjiang17\\My Documents\\GitHub\\StanDCM\\R\\";setwd(targetPath) #AllFiles<-list.files(targetPath) #AllFiles<-AllFiles[AllFiles!="StanDCM.R"] #for(i in 1:length(AllFiles)){ # line <- readLines(AllFiles[i]) # write(line, "StanDCM.txt",append = T) # write("\n#SEPERATION#", "StanDCM.txt",append = T) #} #SEPERATION# #' @title Generate data list #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param response.matrix save the .stan file to somewhere; the default path is getwd() #' @param GroupID name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\data.RData") ;Qmatrix<-cbind(Qmatrix,rep(1,9));Qmatrix[1,1]<-0 Generate.datalist<-function(Qmatrix,response.matrix,GroupID=NA,SubjectID=NA,TimestampID=NA,quad.structure=F){ if(quad.structure==T){Nz=3}else{Nz=2} if(!is.na(TimestampID)[1]&!is.na(SubjectID)[1]){ if(sum(is.na(response.matrix))!=0){print('Stop!The response dataset contains missing value(s)')}else{ Generate.dataList<-list(Y=response.matrix, Nr=nrow(response.matrix), Na=ncol(Qmatrix), Np = length(unique(SubjectID)), Ni = ncol(response.matrix), Nz = Nz, No = length(unique(TimestampID)), SubjectID = SubjectID, OccasionID = rank(unique(TimestampID))[match(TimestampID,unique(TimestampID))], TimeVec = unique(TimestampID)[order(unique(TimestampID))] ) } }else{ if(sum(is.na(response.matrix))!=0){print('Stop!The response dataset contains missing value(s)')}else{ if( length((unique(unlist(unique(c(response.matrix))))))>2 ){ Generate.dataList<-list(Y=response.matrix, Na=ncol(Qmatrix), Np = nrow(response.matrix), Ni = ncol(response.matrix),Nc=2^(ncol(Qmatrix)), Ns = length((unique(unlist(unique(c(response.matrix))))))) }else{ Generate.dataList<-list(Y=response.matrix, Na=ncol(Qmatrix), Np = nrow(response.matrix), Ni = ncol(response.matrix),Nc=2^(ncol(Qmatrix))) } } if(!is.na(GroupID)[1]){Generate.dataList$GroupID=GroupID} } if(quad.structure==T){Generate.dataList$TimeQuadVec= (unique(TimestampID)[order(unique(TimestampID))])^2} Generate.dataList } #SEPERATION# #' @title Install and load packages. #' #' @description #' \code{Install.package} allows other functions to install and load packages. #' #' @param needed_packages the packages to be installed and loaded. #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @examples #' \dontrun{ #' Install.package("ggplot2") #' } #' Install.package = function(needed_packages){ for (i in 1:length(needed_packages)){ haspackage = require(needed_packages[i], character.only = TRUE) if (haspackage==FALSE){ install.packages(needed_packages[i]) require(needed_packages[i], character.only = TRUE) } } } StanPosterior.item<-function(stan.model,lambda.parm=TRUE){ x<-stan.model summaryStan<-summary(x)$'summary' temp<-!grepl("PImat|log_lik|contributionsI|posteriorIC|Vc|posteriorPC|lp",rownames(summaryStan)) summaryStan<-summaryStan[temp,] if(lambda.parm==TRUE){ summaryStan[grepl("l",rownames(summaryStan)),] }else{ summaryStan[!grepl("l",rownames(summaryStan)),] } } #SEPERATION# #' @title Generate Stan code and Run the estimation for LCDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param savepath save the .stan file to somewhere; the default path is getwd() #' @param savename name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\data.RData") ;Qmatrix<-cbind(Qmatrix,rep(1,9));Qmatrix[1,1]<-0 Parm.name<-function(Qmatrix){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j])} } } Classp.exp1<-colnames(OUTPUT[[1]]) #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) trueParmName<-itemParmName out.list<-list() out.list$parm.name<-trueParmName out.list$intercept.name<-OUTPUT[[5]] out.list$main.name<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1] out.list$interaction.name<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])!=1] out.list$class.expression<-Classp.exp1 out.list } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanCRUM.run<-function(Qmatrix,response.matrix,script.path=NA,save.path=getwd(),save.name="CRUM_uninf",iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ #break stop() } Cdm.init<-F if(init.list=='cdm'){ Cdm.init<-T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =respMatrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list( paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') inilist1<-eval(parse(n =2000000 ,text=inilist1)) for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) StanCRUM.script(Qmatrix,save.path=save.path,save.name=save.name) script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path) }else{ compiled_model<-stan_model(script.path) } if(Cdm.init==T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanCRUM.script<-function(Qmatrix,save.path=getwd(),save.name="CRUM_uninf"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j])} } } #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',intercept,';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp[loopi,loopc]),');\n',sep='') } } ####################053119update#################### for(loopi in 1:nr){ for(loopc in 1:Nc){ Reparm[loopi,loopc]<- str_replace_all(Reparm[loopi,loopc],paste((name.inter),collapse="|"),"0") } } Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',c(itemParmName[1:numMainEffect],intercept),'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) ####################053119update END ############## #Likelihood Stan code Likelihood<-' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); } contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ' #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; } ' #Parameter Specification parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') #Reparameter Specification transparm.spec<-paste(c(' transformed parameters{ matrix[Ni, Nc] PImat;\n', paste0(unlist(Reparm)),'}\n'),collapse='') #Model Specification model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,Parmprior,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] #Generated Quantities Specification generatedQuantities.spec<-' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ' if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file= filename,append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanCRUM_mG.run<-function(Qmatrix, response.matrix, GroupID, fixeditem.vector=NA, class.equal=T, script.path=NA,save.path=getwd(),save.name="CRUM_uninf_multiG", iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ group.num<-length(unique(GroupID)) rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init<-F if(init.list=='cdm'){ Cdm.init<-T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =response.matrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list( paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') IniList1<-NULL if(!class.equal){ temp.inilist1<-eval(parse(n =2000000 ,text=inilist1)) eval(parse(n =2000000 ,text=paste(paste('IniList1$Vc_g',1:group.num,sep=''),"<-temp.inilist1$Vc",sep='') )) inilist1<-IniList1 }else{ inilist1<-eval(parse(n =2000000 ,text=inilist1))} for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix,GroupID) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) #Need to update script StanCRUM_mG.script(Qmatrix=Qmatrix, group.num=group.num, fixeditem.vector=fixeditem.vector, class.equal=class.equal, save.path=save.path,save.name=save.name) script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path)} else{ compiled_model<-stan_model(script.path) } if(Cdm.init==T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanCRUM_mG.script<-function(Qmatrix, group.num, fixeditem.vector=NA, class.equal=T, save.path=getwd(),save.name="CRUM_uninf_multiG"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate CRUM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] Kernel.exp.detect<-OUTPUT[[1]] #052719updates Kernel.exp.CRUM<-OUTPUT[[1]] #052719updates for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j]) Kernel.exp.detect[i,j]<-NA} #052719updates } } for (i in 1:nrow(OUTPUT[[1]])){ #052719updates theClosestEffect<-which(is.na(Kernel.exp.detect[i,]))[1] #052719updates useToReplaceLonger<-Kernel.exp[i,theClosestEffect] #052719updates Kernel.exp.CRUM[i,is.na(Kernel.exp.detect[i,])]<-useToReplaceLonger #052719updates } #052719updates #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) ############################################################# #######060319update: Multiple Group########################## ############################################################# if(!is.na(fixeditem.vector)[1]){ fixedItem.vector<-c(1:nr)[-fixeditem.vector] }else{fixedItem.vector<-c(1:nr)} ############################################################# #######060719update: Multiple Group########################## ############################################################# Kernel.exp.CRUM<-Kernel.exp for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp[loopi,loopc]),');\n',sep='') } } for(loopi in 1:nr){ for(loopc in 1:Nc){ Reparm[loopi,loopc]<- str_replace_all(Reparm[loopi,loopc],paste((name.inter),collapse="|"),"0") } } Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',c(itemParmName[1:numMainEffect],intercept),'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) ############################################################# #######060719update: Multiple Group End###################### ############################################################# Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) ############################################################# Kernel.exp.CRUM.groupName<-paste("Kernel.exp.CRUM_g",c(1:group.num),sep='') for(i in 1:group.num){ tempfill.Kernel.exp.CRUM<-Kernel.exp.CRUM temp.Kernel.exp.CRUM<-Kernel.exp.CRUM[fixedItem.vector,] for(j in 1:nrow(temp.Kernel.exp.CRUM)){ for(z in 1:ncol(temp.Kernel.exp.CRUM)){ temp.Kernel.exp.CRUM[j,z]<-paste(temp.Kernel.exp.CRUM[j,z],'_g',i,sep='') } } for(j in 1:nrow(temp.Kernel.exp.CRUM)){ for(z in 1:ncol(temp.Kernel.exp.CRUM)){ temp.Kernel.exp.CRUM[j,z]<-str_replace_all(temp.Kernel.exp.CRUM[j,z],"\\+",paste("_g",i,"+",sep='')) } } tempfill.Kernel.exp.CRUM[fixedItem.vector,]<-temp.Kernel.exp.CRUM assign(Kernel.exp.CRUM.groupName[i],tempfill.Kernel.exp.CRUM) } Kernel.exp.CRUM.list<-list() for(i in 1:group.num){Kernel.exp.CRUM.list[[i]]<-eval(parse(text=paste("Kernel.exp.CRUM_g",i,sep='')))} ############################################################# ##########060719update: Multiple Group End################### ############################################################# ############################################################# ##########060719update: Multiple Group####################### ############################################################# PImat.groupName<-paste("PImat_g",c(1:group.num),sep='') Reparm.multigroup<-array(0,dim = c(nr,nclass,group.num)) #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ for (loopg in 1:group.num){ Reparm.multigroup[loopi,loopc,loopg]<-paste(' PImat_g',loopg,'[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp.CRUM.list[[loopg]][loopi,loopc]),');\n',sep='') } } } ############################################################# ##########060719update: Multiple Group ###################### ############################################################# intercept.multigroup<-array(0,dim = c(nr,1,group.num)) mainEff.multigroup<-array(0,dim = c(numMainEffect,1,group.num)) #Group Invariant Parameter Name fixedParmName<-NULL if(!is.na(fixeditem.vector)[1]){ for(i in fixeditem.vector){ fixedParmName<-c(fixedParmName,out[[3]][[i]]) } fixedParmName<-c(fixedParmName,out[[5]][fixeditem.vector]) } #Group Variant Parameter Name freeParmName<-c(out[[5]],itemParmName)[!c(out[[5]],itemParmName)%in%fixedParmName] for(i in 1:group.num){ tempfill.intercept<-out[[5]] tempfill.mainEff<-itemParmName[1:numMainEffect] temp.intercept<-tempfill.intercept[fixedItem.vector] temp.mainEff<-tempfill.mainEff[!tempfill.mainEff%in%fixedParmName] for(j in 1:length(temp.intercept)){ temp.intercept[j]<-paste(temp.intercept[j],"_g",i,sep='') } for(j in 1:length(temp.mainEff)){ temp.mainEff[j]<-paste(temp.mainEff[j],"_g",i,sep='') } tempfill.intercept[fixedItem.vector]<-temp.intercept tempfill.mainEff[!tempfill.mainEff%in%fixedParmName]<-temp.mainEff intercept.multigroup[,,i]<-tempfill.intercept mainEff.multigroup[,,i]<-tempfill.mainEff } Constrain.List<-paste(' real',unique(mainEff.multigroup),';\n ') Unconstrain.List<-paste(' real',unique(c(intercept.multigroup)),';\n ') ############################################################# ##########060319update: Multiple Group End################### ############################################################# ############################################################# #####060719update:Change from update.Parmprior&fix ###### ############################################################# update.Parmprior.multiGroup<-paste(paste(' ',unique(c(intercept.multigroup)),'~normal(0,5)',';\n',sep=''), paste(' ',unique(c(mainEff.multigroup)),'~normal(0,5)',';\n',sep='')) update.Parmprior.multiGroup<-unique(update.Parmprior.multiGroup) update.Parmprior.multiGroup<-c(" //Prior\n ",paste(' ',fixedParmName,'~normal(0,5)',';\n',sep=''),update.Parmprior.multiGroup ) if(class.equal){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup,paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep='') ) }else{ for(i in 1:group.num){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup, paste(' Vc_g',i,'~dirichlet(rep_vector(2.0, Nc));\n',sep='') ) } } ##therefore we can use: fix.Parmprior,update.Parmprior ############################################################# ############################################################# ############################################################# #####060419update:Likelihood Add PImat_g##################### ############################################################# PImat.likelihood1<-NULL PImat.likelihood0<-NULL Vc.likelihood<-NULL for(loopg in 1:group.num){ temp.PImat.likelihood1<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(1|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood1<-paste(PImat.likelihood1,temp.PImat.likelihood1,sep='\n') temp.PImat.likelihood0<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(0|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood0<-paste(PImat.likelihood0,temp.PImat.likelihood0,sep='\n') } if(!class.equal){ for(loopg in 1:group.num){ temp.Vc.likelihood<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsC[iterc]=log(Vc_g',loopg,'[iterc])+sum(contributionsI);',sep=''), sep='\n') Vc.likelihood<-paste(Vc.likelihood,temp.Vc.likelihood,sep='\n') } } #Likelihood Stan code if(class.equal){ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '} contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ',sep='')}else{ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '}\n', Vc.likelihood ,' } target+=log_sum_exp(contributionsC); } ',sep='') } ############################################################# #####060419update:Likelihood Add PImat_g end################ ############################################################# #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; vector[Np] GroupID; } ' ############################################################# #####060419update: Stan script############################## ############################################################# #Parameter Specification if(class.equal){parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='')}else{ parm.spec<-paste(c(' parameters{\n ', paste(paste(' simplex[Nc] Vc_g',1:group.num, ";",sep=''),"\n"), paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') } #Reparameter Specification transparm.spec<-paste(c(' transformed parameters{ ', paste(' matrix[Ni, Nc] PImat_g',1:group.num,';\n',sep=''), paste0(c(Reparm.multigroup)),'}\n'),collapse='') #Model Specification update052619 model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,update.Parmprior.multiGroup,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] #Generated Quantities Specification IC.generatedquantities<-NULL if(!class.equal){ for(loopg in 1:group.num){ temp.IC.generatedquantities<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsIC[iteri,iterc]=log(Vc_g',loopg,'[iterc])+contributionsI[iteri];',sep=''), sep='\n') IC.generatedquantities<-paste(IC.generatedquantities,temp.IC.generatedquantities,sep='\n') } } if(class.equal){ generatedQuantities.spec<-paste(' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, ' contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ',sep='')}else{ generatedQuantities.spec <- paste( ' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)', PImat.likelihood1,'\n', ' else', PImat.likelihood0,'\n', " ",IC.generatedquantities,'\n', ' posteriorIC[iteri,iterc]=contributionsI[iteri];','\n }\n', ' log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ', sep = '' ) } if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file=filename, append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# #' @title A function to generate leave-one-out cross-validation for Stan Model #' #' @description #' The StanLCDM.loofit Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param savepath save the .stan file to somewhere; the default path is getwd() #' @param savename name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export StanDCM.fit <- function(x, pars ="log_lik", cores = 2, save_psis =TRUE) { loo1 <- loo(x = x, pars = pars, save_psis = save_psis, cores = cores) loo1 } #SEPERATION# #' @title A function to calculate pecentage of extreme p-values for sumscore distribution #' #' @description #' The StanDCM.ppmc Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param stan.model A rStan object #' @param response.matrix the response matrix used by RStan Object #' @param n.sim number of simulations for Posterior Predictive Model Checking #' @param n.burnin number of burn-ins #' @param plot.option logical. whether to provide a plot for ppmc #' @return p-values tables #' #' @author {Jihong Zhang, University of Iowa, \email{jihong-zhang@uiowa.edu}} #' #' @export #' @examples #' load("data.RData") #' Qmatrix<-cbind(Qmatrix,rep(1,9)); Qmatrix[1,1]<-0 #' dim(respMatrix) #' misspecifiedQmatrix <- Qmatrix #' misspecifiedQmatrix[1:6,] <- 1-Qmatrix[1:6,] #' misspecifiedQmatrix[1,3] = 0 #' mod2 <- StanDINA.run(misspecifiedQmatrix,response.matrix = respMatrix,iter=100,init.list='cdm', chain.num = 3, warmup = 20) #' StanDCM.ppmc(stan.model = mod2, response.matrix = respMatrix, n.sim = 1000, n.burnin = 1, plot.option = FALSE) #' end - start StanDCM.ppmc <- function(stan.model, response.matrix, n.sim = 6000, n.burnin = 20, plot.option=FALSE) { Install.package(c("Rlab", "MCMCpack", "tidyr", "dplyr", "pbapply")) if(plot.option == TRUE) Install.package("ggplot2") mod1 <- stan.model posterior.matrix <- as.matrix(mod1) Ni <- ncol(response.matrix) Np <- nrow(response.matrix) PImat.name <- grep(names(mod1), pattern = "PImat", value = T) Nc <- length(PImat.name) / Ni Vc.name <- grep(names(mod1), pattern = "Vc", value = T) n.sim = n.sim n.iter = mod1@stan_args[[1]]$iter if (n.burnin >= n.iter) { stop('Number of iteration less than number of burn-ins') } Vc.posterior.matrix <- posterior.matrix[(n.burnin +1): n.iter, Vc.name] PImat.posterior.matrix <- posterior.matrix[(n.burnin +1): n.iter, PImat.name] # simulate response matrix based on PImat and Vc pseudo.sumscore.dist.matrix = NULL if (is.na(n.sim)) { n.sim = Np } time = 0 # create progress bar # pb <- txtProgressBar(min = 0, max = n.sim, style = 3) pseudo.sumscore.extract <- function() { time <<- time + 1 iter <- sample(1:(n.iter-n.burnin), 1, replace = TRUE) # select Predicted Probability Matrix PImat.select <- matrix(PImat.posterior.matrix[iter,], Ni, Nc) # draw from dirichlet distribution and generate class for each person Vc.draw <- rdirichlet(1, rep(1, Nc)) pseudo.class.vector<-apply(rmultinom(Np, 1, Vc.draw), 2, function(x){which(x==1)}) pseudo.response.matrix <- t(PImat.select[, pseudo.class.vector]) pseudo.response.matrix <- t(apply(pseudo.response.matrix, 1, function(x) rbinom(Ni, 1, x))) pseudo.sumscore.vector <- apply(pseudo.response.matrix, 1, sum) pseudo.sumscore.dist.vector <- data.frame(table(pseudo.sumscore.vector)) pseudo.sumscore.dist.vector$time = time pseudo.sumscore.dist.vector } # (2) Set up the style for progree bar pboptions(type = "txt", style = 3, char = "=") #（1）extract the pseudo sumscore for n.sim times pseudo.sumscore.dist.matrix = pbreplicate(n.sim, pseudo.sumscore.extract(), simplify = FALSE) pseudo.sumscore.dist.matrix.long <- do.call(rbind,pseudo.sumscore.dist.matrix) pseudo.sumscore.dist.matrix.wide <- spread(key = pseudo.sumscore.vector, value = Freq , pseudo.sumscore.dist.matrix.long) pseudo.sumscore.dist.df <- pseudo.sumscore.dist.matrix.wide[,-1] #(2) Compute the observed raw score distribution rep.quantile <- apply(pseudo.sumscore.dist.df, 2, quantile,na.rm = TRUE, probs = c(0.05, 0.5, 0.95)) rep.quantile <- data.frame(t(rep.quantile)) colnames(rep.quantile) <- c("Q5", "Q50", "Q95") rep.quantile$score <- as.numeric(rownames(rep.quantile)) response.sumscore.obs <- apply(response.matrix, 1, sum) observed <- data.frame(table(response.sumscore.obs)) colnames(observed) <- c("score", "obs") observed$score <- as.numeric(as.character(observed$score)) # (3) combine them together df <- left_join(observed, rep.quantile, by = "score") df$'WithinPosterior(95%)' <- ifelse(df$obs > df$Q95 | df$obs < df$Q5, FALSE, TRUE) if(plot.option == TRUE){ p <- ggplot(df, aes(x = score + 1, y = obs)) + theme_bw() + labs(x = "Raw score", y = "Number of examinees", title = "The observed and replicated raw score distributions") # (2) Add the entire replicated raw score distributions using violin plot # and jittered data points p <- p + geom_violin(data = pseudo.sumscore.dist.matrix.long, aes(x = pseudo.sumscore.vector, y = Freq)) p <- p + geom_jitter(data = pseudo.sumscore.dist.matrix.long, aes(x = pseudo.sumscore.vector, y = Freq, color = pseudo.sumscore.vector), position = position_jitter(width = 0.12), alpha = 0.15) + theme(legend.position = "none") # (3) Overlay 5%, 50% (median), and 95% quantiles of each replicated raw # score distribution p <- p + geom_point(aes(y = Q50), size = 4, shape = 2) + geom_line(aes(y = Q50),size = 1, linetype = "dashed") p <- p + geom_line(aes(y = Q5), size = 1, linetype = "dotted") p <- p + geom_line(aes(y = Q95), size = 1, linetype = "dotted") # (4) Overlay the observed raw score distribution p <- p + geom_point(size = 4) + geom_line(size = 1) print(p) } df } #SEPERATION# #' @title DCM calibration under the DINA model via Stan #' #' #' @description #' \code{StanDINA} uses Stan program to calibrate the deterministic inputs, noisy and gate model for dichotomous responses, and its #' extension #' #' In addition, users are allowed to specify design matrix and link function for each item, and distinct models may be used #' in a single test for different items. The attributions can be either dichomous or polytomous. #' @usage #' StanDINA.run<-function(Qmatrix,response.matrix,script.path=NA,save.path=getwd(),save.name="DINA_uninf",iter=1000,warmup = floor(iter/2), #' chain.num=3,init.list='random',control.list=NA) #' #' @param Qmatrix A required matrix #' @param response.matrix save the .stan file to somewhere; the default path is getwd() #' @param script.path save the .stan file to somewhere; the default path is getwd() #' @param save.name the name of saved stan file #' @param iter number of iteration of MCMC estimation. defalts to 1000. #' @param chain.num number of MCMC chain.num. #' @param init.list the initial values. 'random' or 'CDM' #' @param control.list the controlled parameters #' #' @return StanDINA returens an object of class StanDINA. Methods for StanDINA objects include #' \code{\link{extract}} for extract for extracting various components, \code{\link{coef}} for #' extracting strctural parameters. #' #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu} \cr Jihong Zhang, University of Iowa, \email{jihong-zhang@uiowa.edu}} #' #' @export #' @examples #' \dontrun{ #' #----------- DINA model-----------# #' mod1<-StanDINA.run(Qmatrix, respMatrix, iter=20, init.list='cdm', chain.num = 3) #' summary(mod1) #' } StanDINA.run<-function(Qmatrix,response.matrix, script.path=NA,save.path=getwd(),save.name="DINA_uninf", iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init<-F if(init.list=='cdm'){ Cdm.init<-T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =respMatrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list( paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') inilist1<-eval(parse(n =2000000 ,text=inilist1)) for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) StanDINA.script(Qmatrix,save.path=save.path,save.name=save.name) script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path) }else{ compiled_model<-stan_model(script.path) } if(Cdm.init==T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title A function to generate leave-one-out cross-validation for Stan Model #' #' @description #' The StanLCDM.loofit Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages StanDINA.script<-function(Qmatrix,save.path=getwd(),save.name="DINA_uninf"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j])} } } #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) ############################################################# ###########052619update:The highest-interactionn############ hi.interaction<-rep(1,nr) zero.list<-constraints.list fixparm.vec<-NULL for(i in 1:nr){ hi.interaction[i]<-constraints.list[[i]][length(constraints.list[[i]])] if(length(zero.list[[i]])==1){zero.list[[i]]=NA}else{ zero.list[[i]]<-constraints.list[[i]][ 1:(length(constraints.list[[i]])-1) ] fixparm.vec<-c(fixparm.vec,zero.list[[i]]) } } #intercept, hi.interaction,zero.list/fixparm.vec are what we need Constrain.List<-paste(' real',hi.interaction,';\n ') Unconstrain.List<-paste(' real',intercept,';\n ') ############################################################# ############################################################# #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp[loopi,loopc]),');\n',sep='') } } Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) ############################################################# ###########052619update:The highest-interactionn############ update.Parmprior<-Parmprior fix.Parmprior<-NULL for(i in 1:length(Parmprior)){ if(grepl(paste(fixparm.vec, collapse = "|"), Parmprior[i])){ update.Parmprior[i]<-"" } } fix.Parmprior<-c(paste(' real',fixparm.vec,';\n '), paste(' ',fixparm.vec,"=0",';\n ') ) #########052719update:create g and s parameters gParm<-rep(0,nr) sParm<-rep(0,nr) for(loopi in 1:nr){ gParm[loopi]<-paste(' gParm[',loopi,']=inv_logit(',paste(Kernel.exp[loopi,1]),');\n',sep='') sParm[loopi]<-paste(' sParm[',loopi,']=1-inv_logit(',paste(Kernel.exp[loopi,nclass]),');\n',sep='') } ##therefore we can use: fix.Parmprior,update.Parmprior ############################################################# ############################################################# #Likelihood Stan code Likelihood<-' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); } contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ' #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; } ' #Parameter Specification parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') #Reparameter Specification transparm.spec<-paste(c(' transformed parameters{ matrix[Ni, Nc] PImat; vector[Ni] gParm; vector[Ni] sParm;\n', fix.Parmprior, gParm, #052719update sParm, #052719update paste0(unlist(Reparm)),'}\n'),collapse='') #Model Specification update052619 model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,update.Parmprior,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] #Generated Quantities Specification generatedQuantities.spec<-' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ' if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file=filename,append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanDINA_mG.run<-function(Qmatrix, response.matrix, GroupID, fixeditem.vector=NA, class.equal=T, script.path=NA,save.path=getwd(),save.name="DINA_uninf_multiG", iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ group.num<-length(unique(GroupID)) rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init<-F if(init.list=='cdm'){ Cdm.init<-T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =response.matrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list( paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') IniList1<-NULL if(!class.equal){ temp.inilist1<-eval(parse(n =2000000 ,text=inilist1)) eval(parse(n =2000000 ,text=paste(paste('IniList1$Vc_g',1:group.num,sep=''),"<-temp.inilist1$Vc",sep='') )) inilist1<-IniList1 }else{ inilist1<-eval(parse(n =2000000 ,text=inilist1))} for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix,GroupID) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) #Need to update script StanDINA_mG.script(Qmatrix=Qmatrix, group.num=group.num, fixeditem.vector=fixeditem.vector, class.equal=class.equal, save.path=save.path,save.name=save.name) script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path)} else{ compiled_model<-stan_model(script.path) } if(Cdm.init==T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanDINA_mG.script<-function(Qmatrix, group.num, fixeditem.vector=NA, class.equal=T, save.path=getwd(),save.name="DINA_uninf_multiG"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] Kernel.exp.detect<-OUTPUT[[1]] #052719updates Kernel.exp.dina<-OUTPUT[[1]] #052719updates for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j]) Kernel.exp.detect[i,j]<-NA} #052719updates } } for (i in 1:nrow(OUTPUT[[1]])){ #052719updates theClosestEffect<-which(is.na(Kernel.exp.detect[i,]))[1] #052719updates useToReplaceLonger<-Kernel.exp[i,theClosestEffect] #052719updates Kernel.exp.dina[i,is.na(Kernel.exp.detect[i,])]<-useToReplaceLonger #052719updates } #052719updates #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) ############################################################# ###########052619update:The highest-interactionn############ hi.interaction<-rep(1,nr) zero.list<-constraints.list fixparm.vec<-NULL for(i in 1:nr){ hi.interaction[i]<-constraints.list[[i]][length(constraints.list[[i]])] if(length(zero.list[[i]])==1){zero.list[[i]]=NA}else{ zero.list[[i]]<-constraints.list[[i]][ 1:(length(constraints.list[[i]])-1) ] fixparm.vec<-c(fixparm.vec,zero.list[[i]]) } } #intercept, hi.interaction,zero.list/fixparm.vec are what we need Constrain.List<-paste(' real',hi.interaction,';\n ') Unconstrain.List<-paste(' real',intercept,';\n ') ############################################################# #######060319update: Multiple Group########################## ############################################################# if(!is.na(fixeditem.vector)[1]){ fixedItem.vector<-c(1:nr)[-fixeditem.vector] }else{fixedItem.vector<-c(1:nr)} ############################################################# #######060719update: Multiple Group########################## ############################################################# Kernel.exp.dina<-Kernel.exp for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp[loopi,loopc]),');\n',sep='') } } for(loopi in 1:nr){ for( loopc in 1:nclass){ Kernel.exp.dina[loopi,loopc]<-str_replace_all(Kernel.exp.dina[loopi,loopc],paste(fixparm.vec, collapse = "|"),"0") Kernel.exp.dina[loopi,loopc]<-str_replace_all(Kernel.exp.dina[loopi,loopc],"\\+0","") Reparm[loopi,loopc]<-str_replace_all(Reparm[loopi,loopc],paste(fixparm.vec, collapse = "|"),"0") Reparm[loopi,loopc]<-str_replace_all(Reparm[loopi,loopc],"\\+0","") } } ############################################################# #######060719update: Multiple Group End###################### ############################################################# Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) ############################################################# Kernel.exp.dina.groupName<-paste("Kernel.exp.dina_g",c(1:group.num),sep='') for(i in 1:group.num){ tempfill.Kernel.exp.dina<-Kernel.exp.dina temp.Kernel.exp.dina<-Kernel.exp.dina[fixedItem.vector,] for(j in 1:nrow(temp.Kernel.exp.dina)){ for(z in 1:ncol(temp.Kernel.exp.dina)){ temp.Kernel.exp.dina[j,z]<-paste(temp.Kernel.exp.dina[j,z],'_g',i,sep='') } } for(j in 1:nrow(temp.Kernel.exp.dina)){ for(z in 1:ncol(temp.Kernel.exp.dina)){ temp.Kernel.exp.dina[j,z]<-str_replace_all(temp.Kernel.exp.dina[j,z],"\\+",paste("_g",i,"+",sep='')) } } tempfill.Kernel.exp.dina[fixedItem.vector,]<-temp.Kernel.exp.dina assign(Kernel.exp.dina.groupName[i],tempfill.Kernel.exp.dina) } Kernel.exp.dina.list<-list() for(i in 1:group.num){Kernel.exp.dina.list[[i]]<-eval(parse(text=paste("Kernel.exp.dina_g",i,sep='')))} ############################################################# ##########060719update: Multiple Group End################### ############################################################# ############################################################# ##########060719update: Multiple Group####################### ############################################################# PImat.groupName<-paste("PImat_g",c(1:group.num),sep='') Reparm.multigroup<-array(0,dim = c(nr,nclass,group.num)) #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ for (loopg in 1:group.num){ Reparm.multigroup[loopi,loopc,loopg]<-paste(' PImat_g',loopg,'[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp.dina.list[[loopg]][loopi,loopc]),');\n',sep='') } } } ############################################################# ##########060319update: Multiple Group End################### ############################################################# #########052719update:create g and s parameters #########052719update:create g and s parameters gParm<-rep(0,nr) sParm<-rep(0,nr) for(loopi in 1:nr){ gParm[loopi]<-paste(' gParm[',loopi,']=inv_logit(',paste(Kernel.exp.dina[loopi,1]),');\n',sep='') sParm[loopi]<-paste(' sParm[',loopi,']=1-inv_logit(',paste(Kernel.exp.dina[loopi,nclass]),');\n',sep='') } ############################################################# ##########060319update: Multiple Group ###################### ############################################################# gParm.multigroup<-array(0,dim = c(nr,1,group.num)) sParm.multigroup<-array(0,dim = c(nr,1,group.num)) for(i in 1:group.num){ tempfill.gParm<-gParm tempfill.sParm<-sParm tempfill.gParm<-str_replace_all(tempfill.gParm,"gParm",paste("gParm_g",i,sep='')) tempfill.sParm<-str_replace_all(tempfill.sParm,"sParm",paste("sParm_g",i,sep='')) temp.gParm<-gParm[fixedItem.vector] temp.sParm<-sParm[fixedItem.vector] for(j in 1:length(temp.gParm)){ temp.gParm[j]<-str_replace_all(temp.gParm[j],"\\)",paste("_g",i,")",sep='')) temp.sParm[j]<-str_replace_all(temp.sParm[j],"\\)",paste("_g",i,")",sep='')) } for(j in 1:length(temp.gParm)){ temp.gParm[j]<-str_replace_all(temp.gParm[j],"\\+",paste("_g",i,")",sep='')) temp.sParm[j]<-str_replace_all(temp.sParm[j],"\\+",paste("_g",i,"+",sep='')) } temp.gParm<-str_replace_all(temp.gParm,"gParm",paste("gParm_g",i,sep='')) temp.sParm<-str_replace_all(temp.sParm,"sParm",paste("sParm_g",i,sep='')) tempfill.gParm[fixedItem.vector]<-temp.gParm tempfill.sParm[fixedItem.vector]<-temp.sParm gParm.multigroup[,,i]<-tempfill.gParm sParm.multigroup[,,i]<-tempfill.sParm } ############################################################# ##########060319update: Multiple Group ###################### ############################################################# keep.oneMainEffect.multigroup<-NULL intercept.multigroup<-array(0,dim = c(nr,1,group.num)) mainEff.multigroup<-array(0,dim = c(numMainEffect,1,group.num)) interaction.multigroup<-array(0,dim = c(length(name.inter),1,group.num)) #Group Invariant Parameter Name fixedParmName<-NULL if(!is.na(fixeditem.vector)[1]){ for(i in fixeditem.vector){ fixedParmName<-c(fixedParmName,out[[3]][[i]]) } fixedParmName<-c(fixedParmName,out[[5]][fixeditem.vector]) } #Group Variant Parameter Name freeParmName<-c(out[[5]],itemParmName)[!c(out[[5]],itemParmName)%in%fixedParmName] for(i in 1:group.num){ tempfill.intercept<-out[[5]] tempfill.mainEff<-itemParmName[1:numMainEffect] tempfill.interaction<-name.inter tempfill.keep.oneMainEffect<-hi.interaction temp.intercept<-tempfill.intercept[fixedItem.vector] temp.mainEff<-tempfill.mainEff[!tempfill.mainEff%in%fixedParmName] temp.interaction<-tempfill.interaction[!tempfill.interaction%in%fixedParmName] temp.keep.oneMainEffect<-tempfill.keep.oneMainEffect[!tempfill.keep.oneMainEffect%in%fixedParmName] for(j in 1:length(temp.intercept)){ temp.intercept[j]<-paste(temp.intercept[j],"_g",i,sep='') } for(j in 1:length(temp.mainEff)){ temp.mainEff[j]<-paste(temp.mainEff[j],"_g",i,sep='') } for(j in 1:length(temp.interaction)){ temp.interaction[j]<-paste(temp.interaction[j],"_g",i,sep='') } for(j in 1:length(temp.keep.oneMainEffect)){ temp.keep.oneMainEffect[j]<-paste(temp.keep.oneMainEffect[j],"_g",i,sep='') } tempfill.intercept[fixedItem.vector]<-temp.intercept tempfill.mainEff[!tempfill.mainEff%in%fixedParmName]<-temp.mainEff tempfill.interaction[!tempfill.interaction%in%fixedParmName]<-temp.interaction tempfill.keep.oneMainEffect[!tempfill.keep.oneMainEffect%in%fixedParmName]<-temp.keep.oneMainEffect intercept.multigroup[,,i]<-tempfill.intercept mainEff.multigroup[,,i]<-tempfill.mainEff interaction.multigroup[,,i]<-tempfill.interaction keep.oneMainEffect.multigroup<-c(keep.oneMainEffect.multigroup,tempfill.keep.oneMainEffect) } Constrain.List<-paste(' real',unique(keep.oneMainEffect.multigroup),';\n ') Unconstrain.List<-paste(' real',unique(c(intercept.multigroup)),';\n ') ############################################################# ##########060319update: Multiple Group End################### ############################################################# Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) update.Parmprior<-Parmprior fix.Parmprior<-NULL for(i in 1:length(Parmprior)){ if(grepl(paste(fixparm.vec, collapse = "|"), Parmprior[i])){ update.Parmprior[i]<-"" } } update.Parmprior<-update.Parmprior[update.Parmprior!=''] fix.Parmprior<-c(paste(' real',fixparm.vec,';\n '), paste(' ',fixparm.vec,"=0",';\n ') ) ############################################################# #####060319update:Change from update.Parmprior&fix ###### ############################################################# update.Parmprior.multiGroup<-NULL for(i in 1:length(update.Parmprior)){ for (j in 1:length(freeParmName)){ for (z in 1:group.num){ if(sum(grepl(freeParmName[j],update.Parmprior[i]))>=1){ temp.update.Parmprior<-str_replace_all(update.Parmprior[i],freeParmName[j],paste(freeParmName[j],"_g",z,sep='')) update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup, temp.update.Parmprior) } } } } update.Parmprior.multiGroup<-unique(update.Parmprior.multiGroup) update.Parmprior.multiGroup<-c(" //Prior\n",paste(' ',fixedParmName,'~normal(0,5)',';\n',sep=''),update.Parmprior.multiGroup ) if(class.equal){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup,paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep='') ) }else{ for(i in 1:group.num){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup, paste(' Vc_g',i,'~dirichlet(rep_vector(2.0, Nc));\n',sep='') ) } } ##therefore we can use: fix.Parmprior,update.Parmprior ############################################################# ############################################################# ############################################################# #####060419update:Likelihood Add PImat_g##################### ############################################################# PImat.likelihood1<-NULL PImat.likelihood0<-NULL Vc.likelihood<-NULL for(loopg in 1:group.num){ temp.PImat.likelihood1<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(1|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood1<-paste(PImat.likelihood1,temp.PImat.likelihood1,sep='\n') temp.PImat.likelihood0<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(0|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood0<-paste(PImat.likelihood0,temp.PImat.likelihood0,sep='\n') } if(!class.equal){ for(loopg in 1:group.num){ temp.Vc.likelihood<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsC[iterc]=log(Vc_g',loopg,'[iterc])+sum(contributionsI);',sep=''), sep='\n') Vc.likelihood<-paste(Vc.likelihood,temp.Vc.likelihood,sep='\n') } } #Likelihood Stan code if(class.equal){ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '} contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ',sep='')}else{ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '}\n', Vc.likelihood ,' } target+=log_sum_exp(contributionsC); } ',sep='') } ############################################################# #####060419update:Likelihood Add PImat_g end################ ############################################################# #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; vector[Np] GroupID; } ' ############################################################# #####060419update: Stan script############################## ############################################################# #Parameter Specification if(class.equal){parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='')}else{ parm.spec<-paste(c(' parameters{\n ', paste(paste(' simplex[Nc] Vc_g',1:group.num, ";",sep=''),"\n"), paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') } #Reparameter Specification transparm.spec<-paste(c(' transformed parameters{ ', paste(' matrix[Ni, Nc] PImat_g',1:group.num,';\n',sep=''), paste(' vector[Ni] gParm_g',1:group.num,';\n',sep=''), paste(' vector[Ni] sParm_g',1:group.num,';\n',sep=''), c(gParm.multigroup), #060419update c(sParm.multigroup), #060419update paste0(c(Reparm.multigroup)),'}\n'),collapse='') #Model Specification update052619 model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,update.Parmprior.multiGroup,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] #Generated Quantities Specification IC.generatedquantities<-NULL if(!class.equal){ for(loopg in 1:group.num){ temp.IC.generatedquantities<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsIC[iteri,iterc]=log(Vc_g',loopg,'[iterc])+contributionsI[iteri];',sep=''), sep='\n') IC.generatedquantities<-paste(IC.generatedquantities,temp.IC.generatedquantities,sep='\n') } } if(class.equal){ generatedQuantities.spec<-paste(' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, ' contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ',sep='')}else{ generatedQuantities.spec <- paste( ' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)', PImat.likelihood1,'\n', ' else', PImat.likelihood0,'\n', " ",IC.generatedquantities,'\n', ' posteriorIC[iteri,iterc]=contributionsI[iteri];','\n }\n', ' log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ', sep = '' ) } if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file=filename, append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# #' @title DCM calibration under the DINA model via Stan #' #' #' @description #' \code{StanDINA} uses Stan program to calibrate the deterministic inputs, noisy and gate model for dichotomous responses, and its #' extension #' #' In addition, users are allowed to specify design matrix and link function for each item, and distinct models may be used #' in a single test for different items. The attributions can be either dichomous or polytomous. #' @usage #' StanDINA.run<-function(Qmatrix,response.matrix,script.path=NA,save.path=getwd(),save.name="DINA_uninf",iter=1000,warmup = floor(iter/2), #' chain.num=3,init.list='random',control.list=NA) #' #' @param Qmatrix A required matrix #' @param response.matrix save the .stan file to somewhere; the default path is getwd() #' @param script.path save the .stan file to somewhere; the default path is getwd() #' @param save.name the name of saved stan file #' @param iter number of iteration of MCMC estimation. defalts to 1000. #' @param chain.num number of MCMC chain.num. #' @param init.list the initial values. 'random' or 'CDM' #' @param control.list the controlled parameters #' #' @return StanDINA returens an object of class StanDINA. Methods for StanDINA objects include #' \code{\link{extract}} for extract for extracting various components, \code{\link{coef}} for #' extracting strctural parameters. #' #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu} \cr Jihong Zhang, University of Iowa, \email{jihong-zhang@uiowa.edu}} #' #' @export #' @examples #' \dontrun{ #' #----------- DINO model-----------# #' mod1<-StanDINO.run(Qmatrix, respMatrix, iter=20, init.list='cdm', chain.num = 3) #' summary(mod1) #' } StanDINO.run<-function(Qmatrix,response.matrix,script.path=NA,save.path=getwd(),save.name="DINO_uninf",iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init<-F if(init.list=='cdm'){ Cdm.init<-T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =respMatrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list( paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') inilist1<-eval(parse(n =2000000 ,text=inilist1)) for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) StanDINO.script(Qmatrix,save.path=save.path,save.name=save.name) script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path) }else{ compiled_model<-stan_model(script.path) } if(Cdm.init==T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for DINO model #' #' #' @description #' The \code{StanDINO.run} Function allows to automate Stan code geneartion for DINO model with binary resposnes #' #' @usage #' StanDINA.run<-function(Qmatrix,response.matrix,script.path=NA,save.path=getwd(), #' save.name="DINA_uninf",iter=1000,warmup = floor(iter/2), #' chain.num=3,init.list='random',control.list=NA) #' #' @param Qmatrix A required matrix #' @param response.matrix save the .stan file to somewhere; the default path is getwd() #' @param script.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @param iter name the .stan #' @param chain.num name the .stan #' @param init.list name the .stan #' @param control.list name the .stan #' #' @return StanDINA returens an object of class StanDINA. Methods for StanDINA objects include #' \code{\link{extract}} for extract for extracting various components, \code{\link{coef}} for #' extracting strctural parameters. #' #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu} \cr Jihong Zhang, University of Iowa, \email{jihong-zhang@uiowa.edu}} #' #' @export #' @examples #' \dontrun{ #' #----------- DINO model-----------# #' mod1<-StanDINO.run(Qmatrix, respMatrix, iter=20, init.list='cdm', chain.num = 3) #' summary(mod1) #' } StanDINO.script<-function(Qmatrix,save.path=getwd(),save.name="DINO_uninf"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] Kernel.exp.detect<-OUTPUT[[1]] #052719updates Kernel.exp.dino<-OUTPUT[[1]] #052719updates for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j]) Kernel.exp.detect[i,j]<-NA} #052719updates } } for (i in 1:nrow(OUTPUT[[1]])){ #052719updates theClosestEffect<-which(is.na(Kernel.exp.detect[i,]))[1] #052719updates useToReplaceLonger<-Kernel.exp[i,theClosestEffect] #052719updates Kernel.exp.dino[i,is.na(Kernel.exp.detect[i,])]<-useToReplaceLonger #052719updates } #052719updates #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) ############################################################# ###########052719update: Only one main effect is needed ##### keep.oneMainEffect<-rep(1,nr) #052719update: zero.list<-constraints.list fixparm.vec<-NULL #052719update: later they will be fixed to zero for(i in 1:nr){ keep.oneMainEffect[i]<-OUTPUT[[3]][[i]][max(which(OUTPUT[[4]][[i]]==1))]#052719update: if(length(zero.list[[i]])==1){zero.list[[i]]=NA}else{#052719update: zero.list[[i]]<-constraints.list[[i]][-1]#052719update: fixparm.vec<-c(fixparm.vec,zero.list[[i]])#052719update: } } #intercept, hi.interaction,zero.list/fixparm.vec are what we need fixparm.vec<-unlist(constraints.list)[!unlist(constraints.list)%in%keep.oneMainEffect] Constrain.List<-paste(' real',keep.oneMainEffect,';\n ')#052719update: Unconstrain.List<-paste(' real',intercept,';\n ') ############################################################# ############################################################# #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp.dino[loopi,loopc]),');\n',sep='') } } #########052719update:create g and s parameters gParm<-rep(0,nr) sParm<-rep(0,nr) for(loopi in 1:nr){ gParm[loopi]<-paste(' gParm[',loopi,']=inv_logit(',paste(Kernel.exp.dino[loopi,1]),');\n',sep='') sParm[loopi]<-paste(' sParm[',loopi,']=1-inv_logit(',paste(Kernel.exp.dino[loopi,nclass]),');\n',sep='') } Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) ############################################################# ###########052719update: Only one main effect is needed ##### update.Parmprior<-Parmprior fix.Parmprior<-NULL for(i in 1:length(Parmprior)){ if(grepl(paste(fixparm.vec, collapse = "|"), Parmprior[i])){ update.Parmprior[i]<-"" } } update.Parmprior<-update.Parmprior[update.Parmprior!=''] fix.Parmprior<-c(paste(' real',fixparm.vec,';\n '), paste(' ',fixparm.vec,"=0",';\n ') ) ##therefore we can use: fix.Parmprior,update.Parmprior ############################################################# ############################################################# #Likelihood Stan code Likelihood<-' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); } contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ' #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; } ' #Parameter Specification parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') #Reparameter Specification transparm.spec<-paste(c(' transformed parameters{ matrix[Ni, Nc] PImat; vector[Ni] gParm; vector[Ni] sParm;\n', gParm, #052719update sParm, #052719update paste0(unlist(Reparm)),'}\n'),collapse='') #Model Specification update052619 model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,update.Parmprior,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] #Generated Quantities Specification generatedQuantities.spec<-' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ' if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file=filename, append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanDINO_mG.run<-function(Qmatrix, response.matrix, GroupID, fixeditem.vector=NA, class.equal=T, script.path=NA,save.path=getwd(),save.name="DINO_uninf_multiG", iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ group.num<-length(unique(GroupID)) rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init<-F if(init.list=='cdm'){ Cdm.init<-T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =respMatrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list( paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') IniList1<-NULL if(!class.equal){ temp.inilist1<-eval(parse(n =2000000 ,text=inilist1)) eval(parse(n =2000000 ,text=paste(paste('IniList1$Vc_g',1:group.num,sep=''),"<-temp.inilist1$Vc",sep='') )) inilist1<-IniList1 }else{ inilist1<-eval(parse(n =2000000 ,text=inilist1))} for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix,GroupID) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) #Need to update script StanDINO_mG.script(Qmatrix=Qmatrix, group.num=group.num, fixeditem.vector=fixeditem.vector, class.equal=class.equal, save.path=save.path,save.name=save.name) script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path)} else{ compiled_model<-stan_model(script.path) } if(Cdm.init==T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for StanLCDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanDINO_mG.script<-function(Qmatrix, group.num, fixeditem.vector=NA, class.equal=T, save.path=getwd(),save.name="DINO_uninf_multiG"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] Kernel.exp.detect<-OUTPUT[[1]] #052719updates Kernel.exp.dino<-OUTPUT[[1]] #052719updates for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j]) Kernel.exp.detect[i,j]<-NA} #052719updates } } for (i in 1:nrow(OUTPUT[[1]])){ #052719updates theClosestEffect<-which(is.na(Kernel.exp.detect[i,]))[1] #052719updates useToReplaceLonger<-Kernel.exp[i,theClosestEffect] #052719updates Kernel.exp.dino[i,is.na(Kernel.exp.detect[i,])]<-useToReplaceLonger #052719updates } #052719updates #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) ############################################################# ###########052719update: Only one main effect is needed ##### keep.oneMainEffect<-rep(1,nr) #052719update: zero.list<-constraints.list fixparm.vec<-NULL #052719update: later they will be fixed to zero for(i in 1:nr){ keep.oneMainEffect[i]<-OUTPUT[[3]][[i]][max(which(OUTPUT[[4]][[i]]==1))]#052719update: if(length(zero.list[[i]])==1){zero.list[[i]]=NA}else{#052719update: zero.list[[i]]<-constraints.list[[i]][-1]#052719update: fixparm.vec<-c(fixparm.vec,zero.list[[i]])#052719update: } } #intercept, hi.interaction,zero.list/fixparm.vec are what we need fixparm.vec<-unlist(constraints.list)[!unlist(constraints.list)%in%keep.oneMainEffect] Constrain.List<-paste(' real',keep.oneMainEffect,';\n ')#052719update: Unconstrain.List<-paste(' real',intercept,';\n ') ############################################################# #######060319update: Multiple Group########################## ############################################################# if(!is.na(fixeditem.vector)[1]){ fixedItem.vector<-c(1:nr)[-fixeditem.vector] }else{fixedItem.vector<-c(1:nr)} ############################################################# #######060319update: Multiple Group########################## ############################################################# Kernel.exp.dino.groupName<-paste("Kernel.exp.dino_g",c(1:group.num),sep='') for(i in 1:group.num){ tempfill.Kernel.exp.dino<-Kernel.exp.dino temp.Kernel.exp.dino<-Kernel.exp.dino[fixedItem.vector,] for(j in 1:nrow(temp.Kernel.exp.dino)){ for(z in 1:ncol(temp.Kernel.exp.dino)){ temp.Kernel.exp.dino[j,z]<-paste(temp.Kernel.exp.dino[j,z],'_g',i,sep='') } } for(j in 1:nrow(temp.Kernel.exp.dino)){ for(z in 1:ncol(temp.Kernel.exp.dino)){ temp.Kernel.exp.dino[j,z]<-str_replace_all(temp.Kernel.exp.dino[j,z],"\\+",paste("_g",i,"+",sep='')) } } tempfill.Kernel.exp.dino[fixedItem.vector,]<-temp.Kernel.exp.dino assign(Kernel.exp.dino.groupName[i],tempfill.Kernel.exp.dino) } Kernel.exp.dino.list<-list() for(i in 1:group.num){Kernel.exp.dino.list[[i]]<-eval(parse(text=paste("Kernel.exp.dino_g",i,sep='')))} ############################################################# ##########060319update: Multiple Group End################### ############################################################# ############################################################# ##########060319update: Multiple Group####################### ############################################################# PImat.groupName<-paste("PImat_g",c(1:group.num),sep='') Reparm.multigroup<-array(0,dim = c(nr,nclass,group.num)) #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp.dino[loopi,loopc]),');\n',sep='') } } for(loopi in 1:nr){ for( loopc in 1:nclass){ for (loopg in 1:group.num){ Reparm.multigroup[loopi,loopc,loopg]<-paste(' PImat_g',loopg,'[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp.dino.list[[loopg]][loopi,loopc]),');\n',sep='') } } } ############################################################# ##########060319update: Multiple Group End################### ############################################################# #########052719update:create g and s parameters gParm<-rep(0,nr) sParm<-rep(0,nr) for(loopi in 1:nr){ gParm[loopi]<-paste(' gParm[',loopi,']=inv_logit(',paste(Kernel.exp.dino[loopi,1]),');\n',sep='') sParm[loopi]<-paste(' sParm[',loopi,']=1-inv_logit(',paste(Kernel.exp.dino[loopi,nclass]),');\n',sep='') } ############################################################# ##########060319update: Multiple Group ###################### ############################################################# gParm.multigroup<-array(0,dim = c(nr,1,group.num)) sParm.multigroup<-array(0,dim = c(nr,1,group.num)) for(i in 1:group.num){ tempfill.gParm<-gParm tempfill.sParm<-sParm tempfill.gParm<-str_replace_all(tempfill.gParm,"gParm",paste("gParm_g",i,sep='')) tempfill.sParm<-str_replace_all(tempfill.sParm,"sParm",paste("sParm_g",i,sep='')) temp.gParm<-gParm[fixedItem.vector] temp.sParm<-sParm[fixedItem.vector] for(j in 1:length(temp.gParm)){ temp.gParm[j]<-str_replace_all(temp.gParm[j],"\\)",paste("_g",i,")",sep='')) temp.sParm[j]<-str_replace_all(temp.sParm[j],"\\)",paste("_g",i,")",sep='')) } for(j in 1:length(temp.gParm)){ temp.gParm[j]<-str_replace_all(temp.gParm[j],"\\+",paste("_g",i,")",sep='')) temp.sParm[j]<-str_replace_all(temp.sParm[j],"\\+",paste("_g",i,"+",sep='')) } temp.gParm<-str_replace_all(temp.gParm,"gParm",paste("gParm_g",i,sep='')) temp.sParm<-str_replace_all(temp.sParm,"sParm",paste("sParm_g",i,sep='')) tempfill.gParm[fixedItem.vector]<-temp.gParm tempfill.sParm[fixedItem.vector]<-temp.sParm gParm.multigroup[,,i]<-tempfill.gParm sParm.multigroup[,,i]<-tempfill.sParm } ############################################################# ##########060319update: Multiple Group ###################### ############################################################# keep.oneMainEffect.multigroup<-NULL intercept.multigroup<-array(0,dim = c(nr,1,group.num)) mainEff.multigroup<-array(0,dim = c(numMainEffect,1,group.num)) interaction.multigroup<-array(0,dim = c(length(name.inter),1,group.num)) #Group Invariant Parameter Name fixedParmName<-NULL if(!is.na(fixeditem.vector)[1]){ for(i in fixeditem.vector){ fixedParmName<-c(fixedParmName,out[[3]][[i]]) } fixedParmName<-c(fixedParmName,out[[5]][fixeditem.vector]) } #Group Variant Parameter Name freeParmName<-c(out[[5]],itemParmName)[!c(out[[5]],itemParmName)%in%fixedParmName] for(i in 1:group.num){ tempfill.intercept<-out[[5]] tempfill.mainEff<-itemParmName[1:numMainEffect] tempfill.interaction<-name.inter tempfill.keep.oneMainEffect<-keep.oneMainEffect temp.intercept<-tempfill.intercept[fixedItem.vector] temp.mainEff<-tempfill.mainEff[!tempfill.mainEff%in%fixedParmName] temp.interaction<-tempfill.interaction[!tempfill.interaction%in%fixedParmName] temp.keep.oneMainEffect<-tempfill.keep.oneMainEffect[!tempfill.keep.oneMainEffect%in%fixedParmName] for(j in 1:length(temp.intercept)){ temp.intercept[j]<-paste(temp.intercept[j],"_g",i,sep='') } for(j in 1:length(temp.mainEff)){ temp.mainEff[j]<-paste(temp.mainEff[j],"_g",i,sep='') } for(j in 1:length(temp.interaction)){ temp.interaction[j]<-paste(temp.interaction[j],"_g",i,sep='') } for(j in 1:length(temp.keep.oneMainEffect)){ temp.keep.oneMainEffect[j]<-paste(temp.keep.oneMainEffect[j],"_g",i,sep='') } tempfill.intercept[fixedItem.vector]<-temp.intercept tempfill.mainEff[!tempfill.mainEff%in%fixedParmName]<-temp.mainEff tempfill.interaction[!tempfill.interaction%in%fixedParmName]<-temp.interaction tempfill.keep.oneMainEffect[!tempfill.keep.oneMainEffect%in%fixedParmName]<-temp.keep.oneMainEffect intercept.multigroup[,,i]<-tempfill.intercept mainEff.multigroup[,,i]<-tempfill.mainEff interaction.multigroup[,,i]<-tempfill.interaction keep.oneMainEffect.multigroup<-c(keep.oneMainEffect.multigroup,tempfill.keep.oneMainEffect) } Constrain.List<-paste(' real',unique(keep.oneMainEffect.multigroup),';\n ') Unconstrain.List<-paste(' real',unique(c(intercept.multigroup)),';\n ') ############################################################# ##########060319update: Multiple Group End################### ############################################################# Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) update.Parmprior<-Parmprior fix.Parmprior<-NULL for(i in 1:length(Parmprior)){ if(grepl(paste(fixparm.vec, collapse = "|"), Parmprior[i])){ update.Parmprior[i]<-"" } } update.Parmprior<-update.Parmprior[update.Parmprior!=''] fix.Parmprior<-c(paste(' real',fixparm.vec,';\n '), paste(' ',fixparm.vec,"=0",';\n ') ) ############################################################# #####060319update:Change from update.Parmprior&fix ###### ############################################################# update.Parmprior.multiGroup<-NULL for(i in 1:length(update.Parmprior)){ for (j in 1:length(freeParmName)){ for (z in 1:group.num){ if(sum(grepl(freeParmName[j],update.Parmprior[i]))>=1){ temp.update.Parmprior<-str_replace_all(update.Parmprior[i],freeParmName[j],paste(freeParmName[j],"_g",z,sep='')) update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup, temp.update.Parmprior) } } } } update.Parmprior.multiGroup<-unique(update.Parmprior.multiGroup) update.Parmprior.multiGroup<-c(" //Prior\n",paste(' ',fixedParmName,'~normal(0,5)',';\n',sep=''),update.Parmprior.multiGroup ) if(class.equal){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup,paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep='') ) }else{ for(i in 1:group.num){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup, paste(' Vc_g',i,'~dirichlet(rep_vector(2.0, Nc));\n',sep='') ) } } ##therefore we can use: fix.Parmprior,update.Parmprior ############################################################# ############################################################# ############################################################# #####060419update:Likelihood Add PImat_g##################### ############################################################# PImat.likelihood1<-NULL PImat.likelihood0<-NULL Vc.likelihood<-NULL for(loopg in 1:group.num){ temp.PImat.likelihood1<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(1|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood1<-paste(PImat.likelihood1,temp.PImat.likelihood1,sep='\n') temp.PImat.likelihood0<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(0|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood0<-paste(PImat.likelihood0,temp.PImat.likelihood0,sep='\n') } if(!class.equal){ for(loopg in 1:group.num){ temp.Vc.likelihood<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsC[iterc]=log(Vc_g',loopg,'[iterc])+sum(contributionsI);',sep=''), sep='\n') Vc.likelihood<-paste(Vc.likelihood,temp.Vc.likelihood,sep='\n') } } #Likelihood Stan code if(class.equal){ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '} contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ',sep='')}else{ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '}\n', Vc.likelihood ,' } target+=log_sum_exp(contributionsC); } ',sep='') } ############################################################# #####060419update:Likelihood Add PImat_g end################ ############################################################# #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; vector[Np] GroupID; } ' ############################################################# #####060419update: Stan script############################## ############################################################# #Parameter Specification if(class.equal){parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='')}else{ parm.spec<-paste(c(' parameters{\n ', paste(paste(' simplex[Nc] Vc_g',1:group.num, ";",sep=''),"\n"), paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') } #Reparameter Specification transparm.spec<-paste(c(' transformed parameters{ ', paste(' matrix[Ni, Nc] PImat_g',1:group.num,';\n',sep=''), paste(' vector[Ni] gParm_g',1:group.num,';\n',sep=''), paste(' vector[Ni] sParm_g',1:group.num,';\n',sep=''), c(gParm.multigroup), #060419update c(sParm.multigroup), #060419update paste0(c(Reparm.multigroup)),'}\n'),collapse='') #Model Specification update052619 model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,update.Parmprior.multiGroup,Likelihood),sep=''),'\n}',sep='')) #Generated Quantities Specification IC.generatedquantities<-NULL if(!class.equal){ for(loopg in 1:group.num){ temp.IC.generatedquantities<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsIC[iteri,iterc]=log(Vc_g',loopg,'[iterc])+contributionsI[iteri];',sep=''), sep='\n') IC.generatedquantities<-paste(IC.generatedquantities,temp.IC.generatedquantities,sep='\n') } } if(class.equal){ generatedQuantities.spec<-paste(' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, ' contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ',sep='')}else{ generatedQuantities.spec <- paste( ' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)', PImat.likelihood1,'\n', ' else', PImat.likelihood0,'\n', " ",IC.generatedquantities,'\n', ' posteriorIC[iteri,iterc]=contributionsI[iteri];','\n }\n', ' log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ', sep = '' ) } if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file=filename, append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# #' @title Generate Stan code and Run the estimation for LCDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\data.RData") StanLCDM.run<-function(Qmatrix,response.matrix,script.path=NA,save.path=getwd(),save.name="LCDM_uninf",iter=1000,warmup = 0, chain.num=3, init.list='random',control.list=NA){ rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init = F if(init.list=='cdm'){ Cdm.init = T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =response.matrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est # remove nagetive values in initial values CDM.parm.est[CDM.parm.est < 0] <- 0.01 parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) parm.ini[abs(parm.ini)>10] <- 10*sign(parm.ini[abs(parm.ini)>10]) parm.ini[parm.ini==0] <- 0.05 CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list(paste(trueParmName,'=',round(parm.ini,2)), paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') inilist1<-eval(parse(n =2000000 ,text=inilist1)) for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) StanLCDM.script(Qmatrix,save.path=save.path,save.name=save.name) if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } script.path<-filename options(warn=0) compiled_model<-stan_model(script.path) }else{ compiled_model<-stan_model(script.path) } if(Cdm.init == T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains= chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains= chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for LCDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanLCDM.script<-function(Qmatrix,save.path=getwd(),save.name="LCDM_uninf"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j])} } } #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp[loopi,loopc]),');\n',sep='') } } Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) #Likelihood Stan code Likelihood<-' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); } contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ' #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; } ' #Parameter Specification parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') #Reparameter Specification transparm.spec<-paste(c(' transformed parameters{ matrix[Ni, Nc] PImat;\n', paste0(unlist(Reparm)),'}\n'),collapse='') #Model Specification model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,Parmprior,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] #Generated Quantities Specification generatedQuantities.spec<-' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ' if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file= filename,append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# StanLCDM_gC.run<-function(Qmatrix,response.matrix,SubjectID=NA,TimestampID=NA,quad.structure=F,variance.equal=F, script.path=NA,save.path=getwd(),save.name="LCDM_gC_uninf",iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ data.list<-Generate.datalist(Qmatrix,respMatrix,SubjectID=SubjectID,TimestampID=TimestampID) time.vector<-data.list$TimeVec if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) StanLCDM_gC.script(Qmatrix,time.vector=time.vector,quad.structure=quad.structure,variance.equal=variance.equal, save.path=save.path,save.name="LCDM_gC_uninf") script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path) }else{ compiled_model<-stan_model(script.path) } estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) estimated_model } #SEPERATION# StanLCDM_gC.script<-function(Qmatrix,time.vector,quad.structure=F,variance.equal=F, save.path=getwd(),save.name="LCDM_gC_uninf"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix);Na<-nc nr<-nrow(Qmatrix); #R commend attrSamplingStan<-NULL for(j in 1:nc){ attrSamplingStan<-c(attrSamplingStan, paste("A",j,"[iterp,itero] ~ bernoulli_logit(Beta0[",j,"]+Beta1[",j,"]*Theta[iterp,itero])",sep='')) } attrSamplingStan<-paste(attrSamplingStan ,"; \n") attrPosteriorSamplingStan<-NULL for(j in 1:nc){ attrPosteriorSamplingStan<-c(attrPosteriorSamplingStan, paste("A",j,"[iterp,itero] <- bernoulli_rng(inv_logit(Beta0[",j,"]+Beta1[",j,"]*Theta[iterp,itero]))",sep='')) } attrPosteriorSamplingStan<-paste(attrPosteriorSamplingStan ,"; \n") defineA<-paste(paste(paste('matrix[Np,No] A',c(1:Na),sep=''),";",sep='')) defineAlpha<-paste(paste(paste('matrix[Np,No] Alpha',c(1:Na),sep=''),";",sep='')) defineAThreshold<-paste(paste("real AThres",1:Na,sep=""),"[No];",sep='') AThresholdPrior<-paste(paste(" AThres",1:Na,sep=''),"~beta(2,2);",sep="") ####0617update: New parameters for GC######### occasion.num<-length(time.vector) group.num<-occasion.num if(quad.structure){ thridOrder.num<-3 timeQuad.vector<-(time.vector)^2 }else{(thridOrder.num<-2) timeQuad.vector<-NA } No<-occasion.num temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] Kernel.exp.detect<-OUTPUT[[1]] #052719updates Kernel.exp.LCDM<-OUTPUT[[1]] #052719updates for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j]) Kernel.exp.detect[i,j]<-NA} #052719updates } } Kernel.exp.LCDM<-Kernel.exp #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real ',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) ############################################################# #######061719update: Attribute Match Table################### ############################################################# substrRight <- function(x, n){ substr(x, nchar(x)-n+1, nchar(x)) } #for main effect Attr4Main<-cbind(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1], paste("A",substrRight(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],1),sep='')) Attr4Main<-cbind(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1], apply(Attr4Main,1,function(x){paste(x,collapse = "*")})) #for interaction effect Attr4Interaction<-NULL if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] Attr4Interaction<-rbind(Attr4Interaction,c(temp.nm, paste(temp.nm,paste(paste("A",temp.subnm,sep=""),collapse = '*'),sep= '*'))) }}else{ Attr4Interaction<-NULL } replaceWithAttribute<-rbind(Attr4Main,Attr4Interaction) replaceWithAttributeLoop<-replaceWithAttribute replaceWithAttributeLoop[,1]<-paste(replaceWithAttribute[,1],"[itero]",sep='') replaceAttributeLoop<-cbind(paste("A",1:nc,sep=''), paste(paste("A",1:nc,sep=''),"[iterp,itero]",sep='')) for(j in 1:nrow(replaceWithAttributeLoop)){ replaceWithAttributeLoop[j,2]<-str_replace_all(replaceWithAttributeLoop[j,2],replaceWithAttribute[j,1],replaceWithAttributeLoop[j,1]) } for(j in 1:nrow(replaceWithAttributeLoop)){ for(z in 1:nc){ replaceWithAttributeLoop[j,2]<-str_replace_all(replaceWithAttributeLoop[j,2],replaceAttributeLoop[z,1],replaceAttributeLoop[z,2]) } } befAttributeKernel<-Kernel.exp.LCDM[,Nc] befAttributeKernel<-paste(befAttributeKernel,"[itero]",sep='') befAttributeKernel<-str_replace_all(befAttributeKernel,"\\+","[itero]+") ############################################################ ############################################################ ####0617update: move after Nt####### ############################################################# #######060719update: Multiple Group########################## ############################################################# aftAttributeKernel<-befAttributeKernel for(j in 1:length(befAttributeKernel)){ for(z in 1:nrow(replaceWithAttribute)){ aftAttributeKernel[j]<-str_replace_all(aftAttributeKernel[j],fixed(replaceWithAttributeLoop[z,1]),replaceWithAttributeLoop[z,2]) } } for(j in 1:nr){ aftAttributeKernel[j]<-paste('Eta[',j,",itero][iterp]=",aftAttributeKernel[j],";",sep='') } aftAttributeKernel<-str_replace_all(aftAttributeKernel,"A","Alpha") ############################################################ #061819update:Parameter definition ############################################################ itemParmStan<-c(paste('real ', Attr4Main[,1],'[No];',sep=''), paste('real ', intercept,'[No];',sep=''), paste('real ', Attr4Interaction[,1],'[No];',sep=''), 'real Beta0[Na];', 'real Beta1[Na];') #cat(paste(itemParmStan,"\n")) Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',c(itemParmName,"Beta0" ,"Beta1" ),'~normal(0,5)',';\n',sep=''))) data.block<- 'data{ int Np; int Nr; int Nz; int No; int Na; int Ni; int TimeVec[No]; int OccasionID[Nr]; int SubjectID[Nr]; int Y[Nr, Ni]; }' parameter.block<-paste("\nparameters{\n matrix[Np,Nz] Zeta; matrix[Np,No] Theta;", '\n', paste(defineAThreshold,collapse='\n'), '\n', paste(itemParmStan,collapse='\n'), '\n', '\ncov_matrix[Nz] Sigma;', '\nvector[Nz] Mu;', #R commend if(variance.equal){ '\nreal ResVec;' }else{ '\nreal ResMat[No];' },"}") transformedParm.block<-' transformed parameters{ vector[Nz] zeros; matrix[Nz,Nz] identity_Sigma; matrix[Nz,Nz] identity_Mu; zeros = rep_vector(0, Nz); identity_Sigma =diag_matrix(rep_vector(1.0,Nz)); identity_Mu =diag_matrix(rep_vector(10,Nz)); } ' ###################################################################################### ###################################################################################### if(variance.equal){ ResPrior<-'\nResVec~exponential(.1);' }else{ ResPrior<-'\nResMat~exponential(.1);' } AnAlpha_model<-c(paste(paste("A",1:Na,sep=""), "[iterp,itero]= inv_logit(Beta0[1]+Beta1[1]*Theta[iterp,itero]) ;\n",sep='') ,paste(paste("Alpha",1:Na,sep=""), "[iterp,itero]=ceil(fdim(A1[iterp,itero],AThres1[itero]));\n",sep='')) model.block<-paste( "model{\n", paste(defineAlpha,collapse='\n'), "\n", paste(defineA,collapse='\n'), "\n", "vector[Np] Eta[Ni,No];", "\n", paste(Parmprior,collapse =''), "Sigma~inv_wishart(Nz,identity_Sigma); Mu~multi_normal(zeros,identity_Mu);", '\n', paste(AThresholdPrior,collapse ='\n'), ResPrior, '\n', "for (iterp in 1:Np){ Zeta[iterp,] ~ multi_normal(Mu, Sigma); }", "\n", "for (iterp in 1:Np){ for(itero in 1:No){\n ", paste(" ",AnAlpha_model,collapse =' '), " }\n }\n ", "for (iterp in 1:Np){\n for( itero in 1:No){\n ", paste(aftAttributeKernel,collapse ='\n'), " \n }\n }\n ", if(!quad.structure){ if(variance.equal){ 'for (iterp in 1:Np){ for(itero in 1:No){ Theta[iterp,itero]~normal(Zeta[iterp,1]+Zeta[iterp,2]*TimeVec[itero],ResVec); } }' }else{ 'for (iterp in 1:Np){ for(itero in 1:No){ Theta[iterp,itero]~normal(Zeta[iterp,1]+Zeta[iterp,2]*TimeVec[itero],ResMat[itero]); } }' } }else{ if(variance.equal){ 'for (iterp in 1:Np){ for(itero in 1:No){ Theta[iterp,itero]~normal(Zeta[iterp,1]+Zeta[iterp,2]*TimeVec[itero]+Zeta[iterp,3]*TimeQuadVec[itero],ResVec); } }' }else{ 'for (iterp in 1:Np){ for(itero in 1:No){ Theta[iterp,itero]~normal(Zeta[iterp,1]+Zeta[iterp,2]*TimeVec[itero]+Zeta[iterp,3]*TimeQuadVec[itero],ResMat[itero]); } }' } }, "\n", " for (iterr in 1:Nr){ for (iteri in 1:Ni){ Y[iterr,iteri] ~ bernoulli_logit(Eta[iteri,OccasionID[iterr]][SubjectID[iterr]]); } } } ",sep='') ########################################################################### ####################For posterior############################################# if(variance.equal){ ResPosterior<-'\nreal ResVec;' ResPosteriorSampling<-'\nResVec<-exponential_rng(.1);' }else{ ResPosterior<-'\nreal ResMat[No];' ResPosteriorSampling<-'\nResMat<-exponential_rng(.1);' } if(!quad.structure){ if(variance.equal){ ThetaPosteriorSampling<-'for (iterp in 1:Np){ for(itero in 1:No){ Theta[iterp,itero]<-normal_rng(Zeta[iterp,1]+Zeta[iterp,2]*TimeVec[itero],ResVec); } }' }else{ ThetaPosteriorSampling<-'for (iterp in 1:Np){ for(itero in 1:No){ Theta[iterp,itero]<-normal_rng(Zeta[iterp,1]+Zeta[iterp,2]*TimeVec[itero],ResMat[itero]); } }' } }else{ if(variance.equal){ ThetaPosteriorSampling<-'for (iterp in 1:Np){ for(itero in 1:No){ Theta[iterp,itero]<-normal_rng(Zeta[iterp,1]+Zeta[iterp,2]*TimeVec[itero]+Zeta[iterp,3]*TimeQuadVec[itero],ResVec); } }' }else{ ThetaPosteriorSampling<-'for (iterp in 1:Np){ for(itero in 1:No){ Theta[iterp,itero]<-normal_rng(Zeta[iterp,1]+Zeta[iterp,2]*TimeVec[itero]+Zeta[iterp,3]*TimeQuadVec[itero],ResMat[itero]); } }' }} #####################End: For posterior################################################### ########################################################################################## generatedQuant.block <- paste("\n generated quantities {", "\n", paste(defineAlpha,collapse='\n'), '\n', paste(defineA,collapse='\n'), "\n", "vector[Np] Eta[Ni,No]; real log_lik[Nr,Ni]; ", "\n", "for (iterp in 1:Np){ for(itero in 1:No){\n ", paste(" ",AnAlpha_model,collapse =' '), " }\n }\n ", '\n', "for (iterp in 1:Np){\n for( itero in 1:No){\n ", paste(aftAttributeKernel,collapse ='\n'), " \n }\n }\n ", "for (iterr in 1:Nr){ for (iteri in 1:Ni){ log_lik[iterr,iteri] = bernoulli_log(Y[iterr,iteri],inv_logit(Eta[iteri,OccasionID[iterr]][SubjectID[iterr]])); } } }",sep="") if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file= filename,append=FALSE) cat( paste(c(' ', data.block,parameter.block,transformedParm.block,model.block,generatedQuant.block) )) sink(NULL) } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanLCDM_mG.run<-function(Qmatrix, response.matrix, GroupID, fixeditem.vector=NA, class.equal=T, script.path=NA,save.path=getwd(),save.name="LCDM_uninf_multiG", iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ group.num<-length(unique(GroupID)) rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init<-F if(init.list=='cdm'){ Cdm.init<-T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gLCDM( data =response.matrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list( paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') IniList1<-NULL if(!class.equal){ temp.inilist1<-eval(parse(n =2000000 ,text=inilist1)) eval(parse(n =2000000 ,text=paste(paste('IniList1$Vc_g',1:group.num,sep=''),"<-temp.inilist1$Vc",sep='') )) inilist1<-IniList1 }else{ inilist1<-eval(parse(n =2000000 ,text=inilist1))} for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix,GroupID) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) #Need to update script StanLCDM_mG.script(Qmatrix=Qmatrix, group.num=group.num, fixeditem.vector=fixeditem.vector, class.equal=class.equal, save.path=save.path,save.name=save.name) script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path)} else{ compiled_model<-stan_model(script.path) } if(Cdm.init==T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanLCDM_mG.script<-function(Qmatrix, group.num, fixeditem.vector=NA, class.equal=T, save.path=getwd(),save.name="LCDM_uninf_multiG"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] Kernel.exp.detect<-OUTPUT[[1]] #052719updates Kernel.exp.LCDM<-OUTPUT[[1]] #052719updates for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j]) Kernel.exp.detect[i,j]<-NA} #052719updates } } for (i in 1:nrow(OUTPUT[[1]])){ #052719updates theClosestEffect<-which(is.na(Kernel.exp.detect[i,]))[1] #052719updates useToReplaceLonger<-Kernel.exp[i,theClosestEffect] #052719updates Kernel.exp.LCDM[i,is.na(Kernel.exp.detect[i,])]<-useToReplaceLonger #052719updates } #052719updates #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) ############################################################# #######060319update: Multiple Group########################## ############################################################# if(!is.na(fixeditem.vector)[1]){ fixedItem.vector<-c(1:nr)[-fixeditem.vector] }else{fixedItem.vector<-c(1:nr)} ############################################################# #######060719update: Multiple Group########################## ############################################################# Kernel.exp.LCDM<-Kernel.exp for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp[loopi,loopc]),');\n',sep='') } } ############################################################# #######060719update: Multiple Group End###################### ############################################################# Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) ############################################################# Kernel.exp.LCDM.groupName<-paste("Kernel.exp.LCDM_g",c(1:group.num),sep='') for(i in 1:group.num){ tempfill.Kernel.exp.LCDM<-Kernel.exp.LCDM temp.Kernel.exp.LCDM<-Kernel.exp.LCDM[fixedItem.vector,] for(j in 1:nrow(temp.Kernel.exp.LCDM)){ for(z in 1:ncol(temp.Kernel.exp.LCDM)){ temp.Kernel.exp.LCDM[j,z]<-paste(temp.Kernel.exp.LCDM[j,z],'_g',i,sep='') } } for(j in 1:nrow(temp.Kernel.exp.LCDM)){ for(z in 1:ncol(temp.Kernel.exp.LCDM)){ temp.Kernel.exp.LCDM[j,z]<-str_replace_all(temp.Kernel.exp.LCDM[j,z],"\\+",paste("_g",i,"+",sep='')) } } tempfill.Kernel.exp.LCDM[fixedItem.vector,]<-temp.Kernel.exp.LCDM assign(Kernel.exp.LCDM.groupName[i],tempfill.Kernel.exp.LCDM) } Kernel.exp.LCDM.list<-list() for(i in 1:group.num){Kernel.exp.LCDM.list[[i]]<-eval(parse(text=paste("Kernel.exp.LCDM_g",i,sep='')))} ############################################################# ##########060719update: Multiple Group End################### ############################################################# ############################################################# ##########060719update: Multiple Group####################### ############################################################# PImat.groupName<-paste("PImat_g",c(1:group.num),sep='') Reparm.multigroup<-array(0,dim = c(nr,nclass,group.num)) #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ for (loopg in 1:group.num){ Reparm.multigroup[loopi,loopc,loopg]<-paste(' PImat_g',loopg,'[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp.LCDM.list[[loopg]][loopi,loopc]),');\n',sep='') } } } ############################################################# ##########060719update: Multiple Group ###################### ############################################################# intercept.multigroup<-array(0,dim = c(nr,1,group.num)) mainEff.multigroup<-array(0,dim = c(numMainEffect,1,group.num)) interaction.multigroup<-array(0,dim = c(length(name.inter),1,group.num)) #Group Invariant Parameter Name fixedParmName<-NULL if(!is.na(fixeditem.vector)[1]){ for(i in fixeditem.vector){ fixedParmName<-c(fixedParmName,out[[3]][[i]]) } fixedParmName<-c(fixedParmName,out[[5]][fixeditem.vector]) } #Group Variant Parameter Name freeParmName<-c(out[[5]],itemParmName)[!c(out[[5]],itemParmName)%in%fixedParmName] for(i in 1:group.num){ tempfill.intercept<-out[[5]] tempfill.mainEff<-itemParmName[1:numMainEffect] tempfill.interaction<-name.inter temp.intercept<-tempfill.intercept[fixedItem.vector] temp.mainEff<-tempfill.mainEff[!tempfill.mainEff%in%fixedParmName] temp.interaction<-tempfill.interaction[!tempfill.interaction%in%fixedParmName] for(j in 1:length(temp.intercept)){ temp.intercept[j]<-paste(temp.intercept[j],"_g",i,sep='') } for(j in 1:length(temp.mainEff)){ temp.mainEff[j]<-paste(temp.mainEff[j],"_g",i,sep='') } for(j in 1:length(temp.interaction)){ temp.interaction[j]<-paste(temp.interaction[j],"_g",i,sep='') } tempfill.intercept[fixedItem.vector]<-temp.intercept tempfill.mainEff[!tempfill.mainEff%in%fixedParmName]<-temp.mainEff tempfill.interaction[!tempfill.interaction%in%fixedParmName]<-temp.interaction intercept.multigroup[,,i]<-tempfill.intercept mainEff.multigroup[,,i]<-tempfill.mainEff interaction.multigroup[,,i]<-tempfill.interaction } Constrain.List<-paste(' real',unique(mainEff.multigroup),';\n ') Unconstrain.List<-paste(' real',unique(c(intercept.multigroup,interaction.multigroup)),';\n ') ############################################################# ##########060319update: Multiple Group End################### ############################################################# Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) update.Parmprior<-Parmprior fix.Parmprior<-NULL update.Parmprior<-update.Parmprior[update.Parmprior!=''] ############################################################# #####060719update:Change from update.Parmprior&fix ###### ############################################################# update.Parmprior.multiGroup<-NULL for(i in 1:length(update.Parmprior)){ for (j in 1:length(freeParmName)){ for (z in 1:group.num){ if(sum(grepl(freeParmName[j],update.Parmprior[i]))>=1){ temp.update.Parmprior<-str_replace_all(update.Parmprior[i],freeParmName[j],paste(freeParmName[j],"_g",z,sep='')) update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup, temp.update.Parmprior) } } } } update.Parmprior.multiGroup<-unique(update.Parmprior.multiGroup) update.Parmprior.multiGroup<-c(" //Prior\n",paste(' ',fixedParmName,'~normal(0,5)',';\n',sep=''),update.Parmprior.multiGroup ) if(class.equal){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup,paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep='') ) }else{ for(i in 1:group.num){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup, paste(' Vc_g',i,'~dirichlet(rep_vector(2.0, Nc));\n',sep='') ) } } ##therefore we can use: fix.Parmprior,update.Parmprior ############################################################# ############################################################# ############################################################# #####060419update:Likelihood Add PImat_g##################### ############################################################# PImat.likelihood1<-NULL PImat.likelihood0<-NULL Vc.likelihood<-NULL for(loopg in 1:group.num){ temp.PImat.likelihood1<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(1|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood1<-paste(PImat.likelihood1,temp.PImat.likelihood1,sep='\n') temp.PImat.likelihood0<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(0|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood0<-paste(PImat.likelihood0,temp.PImat.likelihood0,sep='\n') } if(!class.equal){ for(loopg in 1:group.num){ temp.Vc.likelihood<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsC[iterc]=log(Vc_g',loopg,'[iterc])+sum(contributionsI);',sep=''), sep='\n') Vc.likelihood<-paste(Vc.likelihood,temp.Vc.likelihood,sep='\n') } } #Likelihood Stan code if(class.equal){ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '} contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ',sep='')}else{ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '}\n', Vc.likelihood ,' } target+=log_sum_exp(contributionsC); } ',sep='') } ############################################################# #####060419update:Likelihood Add PImat_g end################ ############################################################# #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; vector[Np] GroupID; } ' Constrain.List<-unique(Constrain.List);Unconstrain.List<-unique(Unconstrain.List) ############################################################# #####060419update: Stan script############################## ############################################################# #Parameter Specification if(class.equal){parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='')}else{ parm.spec<-paste(c(' parameters{\n ', paste(paste(' simplex[Nc] Vc_g',1:group.num, ";",sep=''),"\n"), paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') } #Reparameter Specification transparm.spec<-paste(c(' transformed parameters{ ', paste(' matrix[Ni, Nc] PImat_g',1:group.num,';\n',sep=''), paste0(c(unique(Reparm.multigroup) )),'}\n'),collapse='') update.Parmprior.multiGroup<-update.Parmprior.multiGroup[!startsWith(str_remove_all(update.Parmprior.multiGroup," "),"~")] #Model Specification update052619 model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,update.Parmprior.multiGroup,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] #Generated Quantities Specification IC.generatedquantities<-NULL if(!class.equal){ for(loopg in 1:group.num){ temp.IC.generatedquantities<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsIC[iteri,iterc]=log(Vc_g',loopg,'[iterc])+contributionsI[iteri];',sep=''), sep='\n') IC.generatedquantities<-paste(IC.generatedquantities,temp.IC.generatedquantities,sep='\n') } } if(class.equal){ generatedQuantities.spec<-paste(' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, ' contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ',sep='')}else{ generatedQuantities.spec <- paste( ' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)', PImat.likelihood1,'\n', ' else', PImat.likelihood0,'\n', " ",IC.generatedquantities,'\n', ' posteriorIC[iteri,iterc]=contributionsI[iteri];','\n }\n', ' log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ', sep = '' ) } if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file=filename, append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# #' @title Generate Stan code and Run the estimation for NCRUM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanNCRUM.run<-function(Qmatrix,response.matrix,script.path=NA,save.path=getwd(),save.name="NCRUM_uninf",iter=1000,warmup = 0, chain.num=3, init.list='random',control.list=NA){ rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init = F if(init.list=='cdm'){ Cdm.init = T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =response.matrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est # remove nagetive values in initial values CDM.parm.est[CDM.parm.est < 0] <- 0.01 parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) parm.ini[abs(parm.ini)>10] <- 10*sign(parm.ini[abs(parm.ini)>10]) parm.ini[parm.ini==0] <- 0.05 CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list(paste(trueParmName,'=',round(parm.ini,2)), paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') inilist1<-eval(parse(n =2000000 ,text=inilist1)) for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) StanNCRUM.script(Qmatrix,save.path=save.path,save.name=save.name) if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } script.path<-filename options(warn=0) compiled_model<-stan_model(script.path) }else{ compiled_model<-stan_model(script.path) } if(Cdm.init == T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains= chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains= chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanNCRUM.script<-function(Qmatrix,save.path=getwd(),save.name="NCRUM_uninf"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j])} } } #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',c(itemParmName[1:numMainEffect],name.inter),';\n ') Unconstrain.List<-paste(' real',intercept,';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp[loopi,loopc]),');\n',sep='') } } #########053019update:create pistar and rstar parameters########################### pistarParm<-rep(0,nr) for(loopi in 1:nr){ pistarParm[loopi]<-paste(' piParm[',loopi,']=inv_logit(',paste(Kernel.exp[loopi,nclass]),');\n',sep='') } piParm<-str_replace_all(intercept,"l","pi");piParm<-str_replace_all(piParm,"_0","") rParm<-c(itemParmName[1:numMainEffect],name.inter) rParm<-str_replace_all(rParm,"l","r") Kernel.exp.NCRUM<-Kernel.exp for(i in 1:Nc){Kernel.exp.NCRUM[,i]<-piParm} #for(i in 1:nr){ # for (j in 1:Nc){ # Kernel.exp.NCRUM.inv[i,j]<-str_replace_all(Kernel.exp[i,j], # intercept[i], # piParm[i]) #} #} allrParm<-NULL for(i in 1:nc){ allrParm<-cbind(allrParm,paste(paste("r",1:nr,"_1",sep=''),i,sep='') ) } allrParm[!t(apply(allrParm ,1,function(x){x%in%rParm}))]<-1 noneOneallrParm<-allrParm[allrParm!='1'] noneOneAllrParm<-noneOneallrParm noneOneAllrParm<-paste(' real',noneOneAllrParm, ';\n',sep=' ') for(i in 1:Nc){ for(j in 1:nr){ Kernel.exp.NCRUM[j,i]<-paste(paste(allrParm[j,which(strsplit(colnames(Kernel.exp)[i],'')[[1]]=='0')],"*",collapse ='',sep=''), Kernel.exp.NCRUM[j,i],sep='') } } Kernel.exp.NCRUM[,Nc]<-str_replace_all(c(Kernel.exp.NCRUM[,Nc]), "[*]", '') piTransf<-Kernel.exp.NCRUM[,Nc] piTransfwithoutName<-piTransf for(loopi in 1:nr){ piTransf[loopi]<-paste(' piParm[',loopi,']=inv_logit(',paste(Kernel.exp[loopi,Nc]),');\n',sep='') piTransfwithoutName[loopi]<-paste('inv_logit(',paste(Kernel.exp[loopi,Nc]),')',sep='') } rTransf<-NULL for(loopi in 1:nr){ for(loopc in 1:Nc){ if(str_count(Kernel.exp.NCRUM[loopi,loopc],"\\*")==1){ if(strsplit(Kernel.exp.NCRUM[loopi,loopc],"\\*")[[1]][1]=='1'){next} temp_r<-noneOneallrParm[(grep(strsplit(Kernel.exp.NCRUM[loopi,loopc],"\\*")[[1]][1],noneOneallrParm))] if(temp_r%in%noneOneallrParm){ rTransf<-c(rTransf,paste(" ",temp_r,"=",paste( paste('inv_logit(',Kernel.exp[loopi,loopc],")",sep=''), "/", piTransfwithoutName[loopi]) )) } } } } rTransf<-paste(rTransf,');\n',sep='') #########053019update:create pistar and rstar parameters END######### Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) #Likelihood Stan code Likelihood<-' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); } contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ' #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; } ' #Parameter Specification parm.spec<-paste(c('parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') #Reparameter Specification transparm.spec<-paste(c(' transformed parameters{ matrix[Ni, Nc] PImat; vector[Ni] piParm;\n', noneOneAllrParm, piTransf,#053019update rTransf,#053019update paste0(unlist(Reparm)),'}\n'),collapse='') #Model Specification update052619 model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,Parmprior,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] #Generated Quantities Specification generatedQuantities.spec<-' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1) contributionsI[iteri]=bernoulli_lpmf(1|PImat[iteri,iterc]); else contributionsI[iteri]=bernoulli_lpmf(0|PImat[iteri,iterc]); contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ' if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file=filename, append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanNCRUM_mG.run<-function(Qmatrix, response.matrix, GroupID, fixeditem.vector=NA, class.equal=T, script.path=NA,save.path=getwd(),save.name="NCRUM_uninf_multiG", iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ group.num<-length(unique(GroupID)) rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init<-F if(init.list=='cdm'){ Cdm.init<-T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =response.matrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list( paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') IniList1<-NULL if(!class.equal){ temp.inilist1<-eval(parse(n =2000000 ,text=inilist1)) eval(parse(n =2000000 ,text=paste(paste('IniList1$Vc_g',1:group.num,sep=''),"<-temp.inilist1$Vc",sep='') )) inilist1<-IniList1 }else{ inilist1<-eval(parse(n =2000000 ,text=inilist1))} for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix,GroupID) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) #Need to update script StanNCRUM_mG.script(Qmatrix=Qmatrix, group.num=group.num, fixeditem.vector=fixeditem.vector, class.equal=class.equal, save.path=save.path,save.name=save.name) script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path)} else{ compiled_model<-stan_model(script.path) } if(Cdm.init==T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanNCRUM_mG.script<-function(Qmatrix, group.num, fixeditem.vector=NA, class.equal=T, save.path=getwd(),save.name="NCRUM_uninf_multiG"){ #Load packages Install.package("plyr") Install.package('stringr') nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] Kernel.exp.detect<-OUTPUT[[1]] #052719updates Kernel.exp.NCRUM<-OUTPUT[[1]] #052719updates for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j]) Kernel.exp.detect[i,j]<-NA} #052719updates } } for (i in 1:nrow(OUTPUT[[1]])){ #052719updates theClosestEffect<-which(is.na(Kernel.exp.detect[i,]))[1] #052719updates useToReplaceLonger<-Kernel.exp[i,theClosestEffect] #052719updates Kernel.exp.NCRUM[i,is.na(Kernel.exp.detect[i,])]<-useToReplaceLonger #052719updates } #052719updates #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==3], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==4], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==5], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==6], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==7], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==8], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==9], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==10], unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==11],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-as.data.frame(matrix(0,nr,nclass)) #########053019update:create pistar and rstar parameters########################### pistarParm<-rep(0,nr) for(loopi in 1:nr){ pistarParm[loopi]<-paste(' piParm[',loopi,']=inv_logit(',paste(Kernel.exp[loopi,nclass]),');\n',sep='') } piParm<-str_replace_all(intercept,"l","pi");piParm<-str_replace_all(piParm,"_0","") rParm<-c(itemParmName[1:numMainEffect],name.inter) rParm<-str_replace_all(rParm,"l","r") Kernel.exp.NCRUM<-Kernel.exp for(i in 1:Nc){Kernel.exp.NCRUM[,i]<-piParm} #for(i in 1:nr){ # for (j in 1:Nc){ # Kernel.exp.NCRUM.inv[i,j]<-str_replace_all(Kernel.exp[i,j], # intercept[i], # piParm[i]) #} #} ############################################################# ###########060719update:The rParm update##################### allrParm<-NULL for(i in 1:nc){ allrParm<-cbind(allrParm,paste(paste("r",1:nr,"_1",sep=''),i,sep='') ) } allrParm[!t(apply(allrParm ,1,function(x){x%in%rParm}))]<-1 noneOneallrParm<-allrParm[allrParm!='1'] noneOneAllrParm<-noneOneallrParm noneOneAllrParm<-paste(' real',noneOneAllrParm, ';\n',sep=' ') for(i in 1:Nc){ for(j in 1:nr){ Kernel.exp.NCRUM[j,i]<-paste(paste(allrParm[j,which(strsplit(colnames(Kernel.exp)[i],'')[[1]]=='0')],"*",collapse ='',sep=''), Kernel.exp.NCRUM[j,i],sep='') } } Kernel.exp.NCRUM[,Nc]<-str_replace_all(c(Kernel.exp.NCRUM[,Nc]), "[*]", '') piTransf<-Kernel.exp.NCRUM[,Nc] piTransfwithoutName<-piTransf for(loopi in 1:nr){ piTransf[loopi]<-paste(' piParm[',loopi,']=inv_logit(',paste(Kernel.exp[loopi,Nc]),');\n',sep='') piTransfwithoutName[loopi]<-paste('inv_logit(',paste(Kernel.exp[loopi,Nc]),')',sep='') } rTransf<-NULL for(loopi in 1:nr){ for(loopc in 1:Nc){ if(str_count(Kernel.exp.NCRUM[loopi,loopc],"\\*")==1){ if(strsplit(Kernel.exp.NCRUM[loopi,loopc],"\\*")[[1]][1]=='1'){next} temp_r<-noneOneallrParm[(grep(strsplit(Kernel.exp.NCRUM[loopi,loopc],"\\*")[[1]][1],noneOneallrParm))] if(temp_r%in%noneOneallrParm){ rTransf<-c(rTransf,paste(" ",temp_r,"=",paste( paste('inv_logit(',Kernel.exp[loopi,loopc],")",sep=''), "/", piTransfwithoutName[loopi]) )) } } } } rTransf<-paste(rTransf,');\n',sep='') rTransf<-str_replace_all(rTransf,"\\)\\);","\\);") ############################################################# #######060319update: Multiple Group########################## ############################################################# if(!is.na(fixeditem.vector)[1]){ fixedItem.vector<-c(1:nr)[-fixeditem.vector] }else{fixedItem.vector<-c(1:nr)} ############################################################# #######060719update: Multiple Group########################## ############################################################# Kernel.exp.NCRUM<-Kernel.exp for(loopi in 1:nr){ for( loopc in 1:nclass){ Reparm[loopi,loopc]<-paste(' PImat[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp[loopi,loopc]),');\n',sep='') } } ############################################################# #######060719update: Multiple Group End###################### ############################################################# Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) ############################################################# Kernel.exp.NCRUM.groupName<-paste("Kernel.exp.NCRUM_g",c(1:group.num),sep='') for(i in 1:group.num){ tempfill.Kernel.exp.NCRUM<-Kernel.exp.NCRUM temp.Kernel.exp.NCRUM<-Kernel.exp.NCRUM[fixedItem.vector,] for(j in 1:nrow(temp.Kernel.exp.NCRUM)){ for(z in 1:ncol(temp.Kernel.exp.NCRUM)){ temp.Kernel.exp.NCRUM[j,z]<-paste(temp.Kernel.exp.NCRUM[j,z],'_g',i,sep='') } } for(j in 1:nrow(temp.Kernel.exp.NCRUM)){ for(z in 1:ncol(temp.Kernel.exp.NCRUM)){ temp.Kernel.exp.NCRUM[j,z]<-str_replace_all(temp.Kernel.exp.NCRUM[j,z],"\\+",paste("_g",i,"+",sep='')) } } tempfill.Kernel.exp.NCRUM[fixedItem.vector,]<-temp.Kernel.exp.NCRUM assign(Kernel.exp.NCRUM.groupName[i],tempfill.Kernel.exp.NCRUM) } Kernel.exp.NCRUM.list<-list() for(i in 1:group.num){Kernel.exp.NCRUM.list[[i]]<-eval(parse(text=paste("Kernel.exp.NCRUM_g",i,sep='')))} ############################################################# ##########060719update: Multiple Group End################### ############################################################# ############################################################# ##########060719update: Multiple Group####################### ############################################################# PImat.groupName<-paste("PImat_g",c(1:group.num),sep='') Reparm.multigroup<-array(0,dim = c(nr,nclass,group.num)) #Produce Stan code for PImat parameter for(loopi in 1:nr){ for( loopc in 1:nclass){ for (loopg in 1:group.num){ Reparm.multigroup[loopi,loopc,loopg]<-paste(' PImat_g',loopg,'[',loopi,',',loopc,']=inv_logit(',paste(Kernel.exp.NCRUM.list[[loopg]][loopi,loopc]),');\n',sep='') } } } ############################################################# ##########060319update: Multiple Group ###################### ############################################################# piTransf.multigroup<-array(0,dim = c(nr,1,group.num)) rTransf.multigroup<-array(0,dim = c(length(rTransf),1,group.num)) for(i in 1:group.num){ tempfill.piTransf<-piTransf tempfill.rTransf<-rTransf tempfill.piTransf<-str_replace_all(tempfill.piTransf,"piTransf",paste("piParm_g",i,sep='')) #tempfill.rTransf<-str_replace_all(tempfill.rTransf,"rTransf",paste("rTransf_g",i,sep='')) for (j in 1:length(noneOneallrParm)){ for( z in 1:length(tempfill.rTransf)){ tempfill.rTransf[z]<-str_replace_all(tempfill.rTransf[z],noneOneallrParm[j],paste(noneOneallrParm[j],"_g",i,sep='')) } } temp.piTransf<-piTransf fixedItem.vector_rTransf<-NULL; for(j in fixedItem.vector){ fixedItem.vector_rTransf<-c(fixedItem.vector_rTransf, grep(j,unlist(lapply(strsplit(unlist(lapply(strsplit(rTransf,"="),function(x){x[[1]]})),"_"),function(x){x[[1]]})))) } temp.rTransf<-rTransf[fixedItem.vector_rTransf] #060819update:a crtical step to handle piTransf lx_0 doesnt have group but some main effects have grouped difference Kernel.exp.fill<-Kernel.exp Kernel.exp.fill<-unique(Kernel.exp.fill) for (w in 1:nr){ for( j in 1:nclass){ for (z in 1:nr){ #Kernel.exp.fill[w,j]<-str_replace_all(Kernel.exp.fill[w,j],paste(intercept[z],"\\+",sep=''),"") Kernel.exp.fill[w,j]<-str_replace_all(Kernel.exp.fill[w,j],intercept[z],"") } } } Kernel.exp.fill<-strsplit(paste(unlist(c(Kernel.exp.fill[fixedItem.vector,])),collapse="+"),"\\+")[[1]] Kernel.exp.fill<-Kernel.exp.fill[Kernel.exp.fill!=''] Kernel.exp.fill<-unique(Kernel.exp.fill) Kernel.exp.fill<-c(Kernel.exp.fill,intercept[fixedItem.vector]) for(j in 1:length(temp.piTransf)){ temp.piTransf[j]<-str_replace_all(temp.piTransf[j],"piParm",paste("piParm_g",i,"",sep='')) } for(z in 1:length(Kernel.exp.fill)){ for(j in 1:length(temp.piTransf)){ temp.piTransf[j]<-str_replace_all(temp.piTransf[j],Kernel.exp.fill[z],paste(Kernel.exp.fill[z] ,"_g",i,sep='')) } } for(j in 1:length(temp.rTransf)){ # temp.piTransf<-piTransf temp.rTransf[j]<-str_replace_all(temp.rTransf[j],"\\)",paste("_g",i,")",sep='')) temp.rTransf[j]<-str_replace_all(temp.rTransf[j],"\\+",paste("_g",i,"+",sep='')) } for (j in 1:length(noneOneallrParm)){ for( z in 1:length(temp.rTransf)){ temp.rTransf[z]<-str_replace_all(temp.rTransf[z],noneOneallrParm[j],paste(noneOneallrParm[j],"_g",i,sep='')) } } tempfill.piTransf[]<-temp.piTransf tempfill.rTransf[fixedItem.vector_rTransf]<-temp.rTransf piTransf.multigroup[,,i]<-tempfill.piTransf rTransf.multigroup[,,i]<-tempfill.rTransf } ############################################################# ##########060319update: Multiple Group ###################### ############################################################# intercept.multigroup<-array(0,dim = c(nr,1,group.num)) mainEff.multigroup<-array(0,dim = c(numMainEffect,1,group.num)) interaction.multigroup<-array(0,dim = c(length(name.inter),1,group.num)) #Group Invariant Parameter Name fixedParmName<-NULL if(!is.na(fixeditem.vector)[1]){ for(i in fixeditem.vector){ fixedParmName<-c(fixedParmName,out[[3]][[i]]) } fixedParmName<-c(fixedParmName,out[[5]][fixeditem.vector]) } #Group Variant Parameter Name freeParmName<-c(out[[5]],itemParmName)[!c(out[[5]],itemParmName)%in%fixedParmName] for(i in 1:group.num){ tempfill.intercept<-out[[5]] tempfill.mainEff<-itemParmName[1:numMainEffect] tempfill.interaction<-name.inter temp.intercept<-tempfill.intercept[fixedItem.vector] temp.mainEff<-tempfill.mainEff[!tempfill.mainEff%in%fixedParmName] temp.interaction<-tempfill.interaction[!tempfill.interaction%in%fixedParmName] for(j in 1:length(temp.intercept)){ temp.intercept[j]<-paste(temp.intercept[j],"_g",i,sep='') } for(j in 1:length(temp.mainEff)){ temp.mainEff[j]<-paste(temp.mainEff[j],"_g",i,sep='') } for(j in 1:length(temp.interaction)){ temp.interaction[j]<-paste(temp.interaction[j],"_g",i,sep='') } tempfill.intercept[fixedItem.vector]<-temp.intercept tempfill.mainEff[!tempfill.mainEff%in%fixedParmName]<-temp.mainEff tempfill.interaction[!tempfill.interaction%in%fixedParmName]<-temp.interaction intercept.multigroup[,,i]<-tempfill.intercept mainEff.multigroup[,,i]<-tempfill.mainEff interaction.multigroup[,,i]<-tempfill.interaction } Constrain.List<-paste(' real',unique(mainEff.multigroup),';\n ') Constrain.List<-c(Constrain.List,paste(' real',unique(interaction.multigroup),';\n ')) Unconstrain.List<-paste(' real',unique(c(intercept.multigroup)),';\n ') ############################################################# ##########060319update: Multiple Group End################### ############################################################# Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) update.Parmprior<-Parmprior fix.Parmprior<-NULL update.Parmprior<-update.Parmprior[update.Parmprior!=''] ############################################################# #####060319update:Change from update.Parmprior&fix ###### ############################################################# update.Parmprior.multiGroup<-NULL for(i in 1:length(update.Parmprior)){ for (j in 1:length(freeParmName)){ for (z in 1:group.num){ if(sum(grepl(freeParmName[j],update.Parmprior[i]))>=1){ temp.update.Parmprior<-str_replace_all(update.Parmprior[i],freeParmName[j],paste(freeParmName[j],"_g",z,sep='')) update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup, temp.update.Parmprior) } } } } update.Parmprior.multiGroup<-unique(update.Parmprior.multiGroup) update.Parmprior.multiGroup<-c(" //Prior\n",paste(' ',fixedParmName,'~normal(0,5)',';\n',sep=''),update.Parmprior.multiGroup ) if(class.equal){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup,paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep='') ) }else{ for(i in 1:group.num){ update.Parmprior.multiGroup<-c(update.Parmprior.multiGroup, paste(' Vc_g',i,'~dirichlet(rep_vector(2.0, Nc));\n',sep='') ) } } ##therefore we can use: fix.Parmprior,update.Parmprior ############################################################# ############################################################# ############################################################# #####060419update:Likelihood Add PImat_g##################### ############################################################# PImat.likelihood1<-NULL PImat.likelihood0<-NULL Vc.likelihood<-NULL for(loopg in 1:group.num){ temp.PImat.likelihood1<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(1|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood1<-paste(PImat.likelihood1,temp.PImat.likelihood1,sep='\n') temp.PImat.likelihood0<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsI[iteri]=bernoulli_lpmf(0|PImat_g',loopg,'[iteri,iterc]);',sep=''), sep='\n') PImat.likelihood0<-paste(PImat.likelihood0,temp.PImat.likelihood0,sep='\n') } if(!class.equal){ for(loopg in 1:group.num){ temp.Vc.likelihood<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsC[iterc]=log(Vc_g',loopg,'[iterc])+sum(contributionsI);',sep=''), sep='\n') Vc.likelihood<-paste(Vc.likelihood,temp.Vc.likelihood,sep='\n') } } #Likelihood Stan code if(class.equal){ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '} contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ',sep='')}else{ Likelihood<-paste(' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, '}\n', Vc.likelihood ,' } target+=log_sum_exp(contributionsC); } ',sep='') } ############################################################# #####060419update:Likelihood Add PImat_g end################ ############################################################# #Data Specification data.spec<-' data{ int Np; int Ni; int Nc; matrix[Np, Ni] Y; vector[Np] GroupID; } ' ############################################################# #####060419update: Stan script############################## ############################################################# Constrain.List<-unique(Constrain.List);Unconstrain.List<-unique(Unconstrain.List) #Parameter Specification if(class.equal){parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='')}else{ parm.spec<-paste(c(' parameters{\n ', paste(paste(' simplex[Nc] Vc_g',1:group.num, ";",sep=''),"\n"), paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') } #To set rParm to real value rParm.real<- paste('real', str_replace_all(unique(unlist(lapply(strsplit(rTransf.multigroup,"="),function(x){x[[1]]})))," ",""),';\n',sep=' ') #Reparameter Specification tranrTransf.spec<-paste(c(' transformed parameters{ ', paste(' matrix[Ni, Nc] PImat_g',1:group.num,';\n',sep=''), paste(' vector[Ni] piParm_g',1:group.num,';\n',sep=''), rParm.real, unique(c(piTransf.multigroup)), #060419update unique(c(rTransf.multigroup)), #060419update paste0(c(Reparm.multigroup)),'}\n'),collapse='') #Model Specification update052619 model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,update.Parmprior.multiGroup,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] #Generated Quantities Specification IC.generatedquantities<-NULL if(!class.equal){ for(loopg in 1:group.num){ temp.IC.generatedquantities<-paste(paste(' if (GroupID[iterp]==',loopg,')',sep=''), paste(' contributionsIC[iteri,iterc]=log(Vc_g',loopg,'[iterc])+contributionsI[iteri];',sep=''), sep='\n') IC.generatedquantities<-paste(IC.generatedquantities,temp.IC.generatedquantities,sep='\n') } } if(class.equal){ generatedQuantities.spec<-paste(' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)' ,PImat.likelihood1,'\n', ' else' ,PImat.likelihood0, ' contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ',sep='')}else{ generatedQuantities.spec <- paste( ' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ if (Y[iterp,iteri] == 1)', PImat.likelihood1,'\n', ' else', PImat.likelihood0,'\n', " ",IC.generatedquantities,'\n', ' posteriorIC[iteri,iterc]=contributionsI[iteri];','\n }\n', ' log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ', sep = '' ) } if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file=filename, append=FALSE) cat( paste(c(' ', data.spec,parm.spec,tranrTransf.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\data.RData") StanORDM.run<-function(Qmatrix,response.matrix, script.path=NA,save.path=getwd(),save.name="ORDM_uninf", iter=1000,warmup = 0, chain.num=3,init.list='random',control.list=NA){ rstan.detect<-tryCatch(library("rstan"),error=function(e){"rstan is not loaded properly. See https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started for details."}) if(length(rstan.detect)==1){ stop() } Cdm.init<-F if(init.list=='cdm'){ Cdm.init<-T Install.package(c("CDM","stringr")) trueParmName<-Parm.name(Qmatrix=Qmatrix)$parm.name Classp.exp1<-Parm.name(Qmatrix=Qmatrix)$class.expression mod1<-gdina( data =respMatrix, q.matrix = Qmatrix , maxit=700,link = "logit",progress=F) CDMresult<-as.data.frame(coef(mod1)) library(stringr) CDM.parm.name<-paste(paste(paste('l',CDMresult[,3],sep=''),'_',sep=''),str_count(CDMresult$partype.attr,"Attr"),sep='') CDM.parm.name<-paste(CDM.parm.name, unlist(lapply(strsplit(unlist(lapply(strsplit(CDMresult$partype.attr, 'Attr', fixed=FALSE),function(x){paste(x,collapse="")})),'-'),function(x){paste(x,collapse="")})), sep='') CDM.parm.est<-CDMresult$est parm.ini<-round(CDM.parm.est[match(trueParmName,CDM.parm.name)],4) CDM.prop.est<-mod1$attribute.patt prop.ini<-CDM.prop.est[match(Classp.exp1,rownames(CDM.prop.est)),1] inilist1<-paste('list(',paste(noquote(paste(noquote(unlist(list( paste(paste('Vc=c(',paste((prop.ini),collapse=','),')',collapse=',')))) ))),collapse=',') ,')',collapse='') inilist1<-eval(parse(n =2000000 ,text=inilist1)) for( i in 2:chain.num){ temp.text<-paste('inilist',i,"<-inilist1",sep='') eval(parse(text=(temp.text))) } temp.text<-paste('init.list<-list(',paste(paste('inilist',1:chain.num,sep=''),collapse = ","),')',sep='') eval(parse(text=(temp.text))) } data.list<-Generate.datalist(Qmatrix,response.matrix) if(is.na(control.list)){control.list<-list(adapt_delta=0.82)} if(is.na(script.path)==T){ options(warn=-1) StanORDM.script(Qmatrix,save.path=save.path,save.name=save.name) script.path<-paste(paste(save.path,save.name,sep='/'),'.stan',sep='') options(warn=0) compiled_model<-stan_model(script.path) }else{ compiled_model<-stan_model(script.path) } if(Cdm.init==T){ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) }else{ estimated_model<-tryCatch(sampling(compiled_model, data = data.list, iter = iter, init = init.list, warmup = warmup, chains=chain.num, control=control.list), error=function(e){"The estimation process is terminated with errors"}) } estimated_model } #SEPERATION# #' @title Generate Stan code and Run the estimation for ORDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param save.path save the .stan file to somewhere; the default path is getwd() #' @param save.name name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\Data") StanORDM.script<-function(Qmatrix,scale.num,save.path=getwd(),save.name="ORDM_uninf"){ nc<-ncol(Qmatrix) nr<-nrow(Qmatrix) temp.table.col<-unique(apply(combn(rep(c(0,1),nc),nc),2,function(x){paste(x,collapse = "")})) temp.table.col<-temp.table.col[order(temp.table.col)] temp.table<-matrix(0,nr,length(temp.table.col)) colnames(temp.table)<-temp.table.col rownames(temp.table)<-paste('item',c(1:nr),sep='') temp.table<-as.data.frame(temp.table) for (i in 1:nr){ temp.table[i,]<-paste('l',i,'_0',sep='') } intercept<-temp.table[,1] #Generate attribute combinations comb.generator<-function(x.vector){ if(length(x.vector)>1){ temp.attr<-x.vector temp.attr.sav<-NULL for(i in 1:length(temp.attr)){ temp.1<-combn(temp.attr,i) temp.2<-apply(temp.1,2,function(x){paste(x,collapse = "")}) temp.attr.sav<-c(temp.attr.sav,temp.2) } } if(length(x.vector)==1){temp.attr.sav<-x.vector} temp.attr.sav } #vectors needed for combination.generator Item.load.id<-list() for ( i in 1:nr){ Item.load.id[[i]]<-grep('1',Qmatrix[i,])} Attr.load.id<-list() attr.load.id<-matrix(0,length(temp.table.col),nc) for ( i in 1:length(temp.table.col)){ attr.load.id[i,]<-unlist(strsplit(temp.table.col[i],split='')) Attr.load.id[[i]]<-grep('1',attr.load.id[i,]) } #Generate Combination for both Item.load and Attr.load Item.Comb<-list() for ( i in 1:nr){ Item.Comb[[i]]<-comb.generator(Item.load.id[[i]]) } Attr.Comb<-list() for ( i in 2:length(temp.table.col)){ Attr.Comb[[1]]<-0 Attr.Comb[[i]]<-comb.generator(Attr.load.id[[i]]) } constraints.list<-list() nway.inter.list<-list() for(i in 1:nr){ for(a in 2:length(temp.table.col)){ ifzero<-as.numeric(paste(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],collapse='')) if((!is.na(ifzero))){ temp.table[i,a]<-paste(c(temp.table[i,a], paste("S","l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='',collapse='') ),collapse='') if(a==length(temp.table.col)){ nway.inter.list[[i]]<-nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]) constraints.list[[i]]<-paste("l",i,"_",nchar(Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])]),Item.Comb[[i]][Item.Comb[[i]]%in%(Attr.Comb[[a]])],sep='') } } } } #Create Lambda Table Lamda.Table<-temp.table for(i in 1:nr){ for(a in 1:length(Lamda.Table)){ t.ref<-unique(as.character(Lamda.Table[i,])) pos<-c(1:length(t.ref))[Lamda.Table[i,a]==t.ref] temp.table[i,a]<-paste("t",i,"_",pos,sep='')}} #Generate LCDM specification out<-list() out[[1]]<-Lamda.Table out[[2]]<-temp.table out[[3]]<-constraints.list out[[4]]<-nway.inter.list out[[5]]<-intercept OUTPUT<-out nclass<-ncol(OUTPUT[[1]]);Nc<-nclass #Produce kernel expressions across items and attributes Kernel.exp<-OUTPUT[[1]] for (i in 1:nrow(OUTPUT[[1]])){ for ( j in 1:ncol(OUTPUT[[1]])){ if(sum(grep('S',OUTPUT[[1]][i,j]))!=0){Kernel.exp[i,j]<-gsub('S','+',OUTPUT[[1]][i,j])} } } #Monotonicity constraint in terms of the interaction terms of the item effects Constrain.List1<-NULL name.inter<-unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2] numway.inter<-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2] subname.inter<-substr((unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2]), (nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])-unlist(OUTPUT[[4]])[unlist(OUTPUT[[4]])>=2]+1), nchar(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])>=2])) if(length(name.inter)!=0){ for (inter in 1: length(name.inter)){ temp.nw<-numway.inter[inter] temp.nm<-name.inter[inter] temp.subnm<-strsplit(subname.inter[inter],split='')[[1]] temp.sel<-paste(unlist(strsplit(temp.nm,split = '_'))[1],"_",(1:(temp.nw-1)),sep='') first.sel<-unlist(OUTPUT[[3]])[grep(paste((temp.sel),collapse="|"),unlist(OUTPUT[[3]]))] second.sel<-sub(".*_.", "", first.sel) for (sel in 1:length(temp.subnm)){ SEL<-second.sel[sel] Constrain.List1<-rbind( paste(temp.nm,">-(0", paste("+",first.sel[grep(SEL,second.sel)], sep='',collapse=''),")",sep=''),Constrain.List1) } } Constrain.List1<-as.character(Constrain.List1) }else{ Constrain.List1<-NULL } itemParmName<-c(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1],unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==2],OUTPUT[[5]]) numMainEffect<-length(unlist(OUTPUT[[3]])[unlist(OUTPUT[[4]])==1]) Constrain.List<-paste(' real',itemParmName[1:numMainEffect],';\n ') Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Reparm<-array(rep(0,nr*nclass*(scale.num)),dim = c(nr,nclass,(scale.num))) #step parameters StepParmName<-paste('step',1:scale.num-1,sep='') StepParmName[1]<-0 CulStepParmName<-StepParmName ##########0529update Kernal.exp is updated to a list so that list[[?]] is a PImat for a cate Kernel.exp.list<-list(scale.num) for(i in 1:scale.num){ Kernel.exp.list[[i]]<-Kernel.exp } Kernel.exp.list.copy<-Kernel.exp.list for(i in 1:scale.num){ for (j in 1:nr){ for( z in 1:nclass){ Kernel.exp.list[[i]][j,z]<-str_replace_all(Kernel.exp.list[[i]][j,z],Kernel.exp.list.copy[[i]][j,1], paste(Kernel.exp.list.copy[[i]][j,1],"r",i-1,sep='')) } } } for(i in 1:scale.num){ for (j in 1:nr){ for( z in 2:nclass){ Kernel.exp.list[[i]][j,z]<-str_replace_all(Kernel.exp.list[[i]][j,z],"\\+", paste("+",i-1,"*",sep='')) } } } Kernel.exp.list[[1]]<-0 for(loopi in 1:nr){ for( loopc in 1:nclass){ for( loops in 1:1){Reparm[loopi,loopc,loops]<-paste(' PImat[',loopi,',',loopc,'][',loops,']=0;\n',sep='')} for( loops in 2:scale.num){ Reparm[loopi,loopc,loops]<-paste(' PImat[',loopi,',',loopc,'][',loops,']=', paste(Kernel.exp.list[[loops]][loopi,loopc]),';\n',sep='') } } } itemParmName<-c(unlist(OUTPUT[[3]])) for( i in 2:scale.num){ itemParmName<-c(itemParmName,c(Kernel.exp.list[[i]][,1])) } Unconstrain.List<-paste(' real',itemParmName[-(1:numMainEffect)],';\n ') Modelcontainer<-paste(' vector[Nc] contributionsC;\n',' vector[Ni] contributionsI;\n\n',sep='') Parmprior<-paste(c(paste(' //Prior\n'),paste(' ',itemParmName,'~normal(0,5)',';\n',sep=''),paste(' Vc~dirichlet(rep_vector(2.0, Nc));',sep=''))) #Likelihood Stan code Likelihood<-' \n //Likelihood for (iterp in 1:Np){ for (iterc in 1:Nc){ for (iteri in 1:Ni){ contributionsI[iteri]= categorical_lpmf(Y[iterp,iteri]| softmax(((PImat[iteri,iterc])))); } contributionsC[iterc]=log(Vc[iterc])+sum(contributionsI); } target+=log_sum_exp(contributionsC); } ' data.spec<-' data{ int Np; int Ni; int Nc; int Ns; matrix[Np, Ni] Y; } ' parm.spec<-paste(c(' parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List), '}\n'),collapse='') #Parameter Specification parm.spec<-paste(c('parameters{ simplex[Nc] Vc;\n ',paste0(Constrain.List),paste0(Unconstrain.List),'}\n'),collapse='') #Reparameter Specification transparm.spec<-paste(c('transformed parameters{ vector[Ns] PImat[Ni, Nc];\n', paste0(unlist(Reparm)),'}\n'),collapse='') #Model Specification model.spec<-paste(c('\nmodel {\n',paste(c(Modelcontainer,Parmprior,Likelihood),sep=''),'\n}',sep='')) model.spec<-model.spec[!startsWith(str_remove_all(model.spec," "),"~")] generatedQuantities.spec<-' \n generated quantities { vector[Ni] log_lik[Np]; vector[Ni] contributionsI; matrix[Ni,Nc] contributionsIC; matrix[Ni,Nc] posteriorIC; matrix[Np,Nc] posteriorPC; //Posterior for (iterp in 1:Np){ for (iteri in 1:Ni){ for (iterc in 1:Nc){ contributionsI[iteri]= categorical_lpmf(Y[iterp,iteri]| softmax(((PImat[iteri,iterc])))); contributionsIC[iteri,iterc]=log(Vc[iterc])+contributionsI[iteri]; posteriorIC[iteri,iterc]=contributionsI[iteri]; } log_lik[iterp,iteri]=log_sum_exp(contributionsIC[iteri,]); } for (iterc in 1:Nc){posteriorPC[iterp,iterc]=prod(exp(posteriorIC[,iterc]));} } } ' if (.Platform$OS.type == "unix") { filename = paste(paste(save.path,save.name,sep='/'),'.stan',sep='') }else{ filename = paste(paste(save.path,save.name,sep='\\'),'.stan',sep='') } sink(file= filename,append=FALSE) cat( paste(c(' ', data.spec,parm.spec,transparm.spec,model.spec,generatedQuantities.spec) )) sink(NULL) } #SEPERATION# StanPosterior.attribute<-function(stan.model,all.only=T){ nc=length(StanPosterior.class(stan.model)) class.name<-unique(apply(combn(rep(c(0,1),log2(nc)),log2(nc)),2,function(x){paste(x,collapse = "")})) posteriorAttr<-get_posterior_mean(stan.model,pars = c("posteriorPC")) rowNames<-row.names(posteriorAttr) newRowNames<-rowNames for(i in 1: length(newRowNames)){ id<-str_match_all(rowNames[i],"[0-9]+")[[1]][,1] newRowNames[i]<-paste(paste("p",id[1],sep=''),paste("c",class.name[as.numeric(as.character(id[2]))],sep=''),sep='') } np=nrow(posteriorAttr)/nc if(all.only){ posteriorAttr<-unlist(lapply(1:np,function(x){which.max(matrix(posteriorAttr[,ncol(posteriorAttr)],np,nc,byrow = T)[x,])})) posteriorAttr<-class.name[posteriorAttr] } if(!all.only){ row.names(posteriorAttr)<-newRowNames } posteriorAttr } #SEPERATION# StanPosterior.class<-function(stan.model,all.only=T){ res<-get_posterior_mean(stan.model,pars = c("Vc")) nc<-nrow(res) class.name<-unique(apply(combn(rep(c(0,1),log2(nc)),log2(nc)),2,function(x){paste(x,collapse = "")})) rownames(res)<-class.name if(all.only){ res<-res[,ncol(res)] } res } #SEPERATION# #' @title Generate Stan code and Run the estimation for LCDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param savepath save the .stan file to somewhere; the default path is getwd() #' @param savename name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\data.RData") ;Qmatrix<-cbind(Qmatrix,rep(1,9));Qmatrix[1,1]<-0 StanPrior.show<-function(script.path){ Install.package("plyr") Install.package('stringr') readStan.script<- readLines(script.path) readStan.script<-readStan.script[grep("~",readStan.script)] readStan.script<-str_remove_all(readStan.script," ") readStan.script<-str_remove_all(readStan.script,";") readStan.script } #SEPERATION# #' @title Generate Stan code and Run the estimation for LCDM #' #' @description #' The StanLCDM.script Function to automate Stan code geneartion for LCDMs with binary resposnes #' #' @param Qmatrix the Q-matrix specified for the LCDM #' @param savepath save the .stan file to somewhere; the default path is getwd() #' @param savename name the .stan #' @return a. stan file saved at the specified path #' #' @author {Zhehan Jiang, University of Alabama, \email{zjiang17@@ua.edu}} #' #' @export #loading needed packages #load("D:\\Dropbox\\Stan\\R\\data.RData") ;Qmatrix<-cbind(Qmatrix,rep(1,9));Qmatrix[1,1]<-0 StanPrior.update<-function(priorUpdate.matrix,script.path,save.path=getwd(),save.name=NA){ Install.package("plyr") Install.package('stringr') options(warn=-1) if(is.na(save.name)){ saveUpdate.name=paste(format(Sys.time(), "%a%b%Y"),'_priorUpdate.stan',sep='')}else{ saveUpdate.name=paste(paste(save.name,sep='/'),'.stan',sep='') } readStan.script<- readLines(script.path) for(i in 1:nrow(priorUpdate.matrix)){ beReplaced<-readStan.script[str_detect(str_remove_all(str_remove_all(readStan.script," "),";"),priorUpdate.matrix[i,1])][grep("~",readStan.script[str_detect(str_remove_all(str_remove_all(readStan.script," "),";"),priorUpdate.matrix[i,1])]) ] beReplaced<-str_remove_all(str_remove_all(beReplaced," "),";") readStan.script[which(str_remove_all(str_remove_all(readStan.script," "),";")==beReplaced)]<-paste(" ", paste(priorUpdate.matrix[i,1],priorUpdate.matrix[i,2],sep="~") , ';//UpdatedPrior ') } options(warn=0) filename = paste(paste(save.path,saveUpdate.name,sep='/'),sep='') writeLines(readStan.script,filename) }